In this example, NLREG is used to fit an ellipse to a roughly
elliptical pattern of data points (i.e., "elliptical regression").
The ellipse may be shifted from the origin, the semi-major and
semi-minor axis lengths must be determined, and the ellipse may
be tilted at an angle. So there are a total of 5 parameters that
must be calculated:
- The X coordinate of the center of the ellipse.
- The Y coordinate of the center of the ellipse.
- The distance from the center to the perimenter along the X axis.
- The distance from the center to the perimenter along the Y axis.
- The tilt angle of the ellipse (counter-clockwise from X axis).
Here is an example showing a set of data points and a general
ellipse fitted to the points.
The ellipse computed by this example minimizes the
sum of the squared distances from the the perimenter of the elipse
to the data points along a radial line extending from the center of
the ellipse to each data point.
The key formula used in this example is the polar equation for
an ellipse:
Radius = sqrt((xdim2 * ydim2) / ((xdim*sin(θ))2 + (ydim*cos(θ))2))
where Radius is the distance from the center (which in this formula is assumed to be (0,0)),
xdim is the distance from the center to the perimenter along the X axis,
ydim is the distance from the center to the perimenter along the Y axis,
and θ is the radial sweep angle counter-clockwise from the X axis.
Using the polar equation of the ellipse, we can compute the distance from the center to a point on
the perimeter for any angle. We also need to know the distance from the center to a given (x,y)
data point which we can calculate using Pythagorean distance formula:
Distance = sqrt((X-Xcenter)2 + (Y-Ycenter)2)
Once we know these two distances, the goal of the optimization is to minimize the sum of the
squared differences in distances (Distance - Radius). Here is the NLREG program that performs
this:
Title "Ellipse fitted to data points";
Variables X, Y;
/*
* Definite parameters to be calculated.
* Specifying reasonable starting values greatly increases the chances for convergence.
*/
Parameter Xcenter = 2; /* X coordinate of center of ellipse */
Parameter Ycenter = 2; /* Y coordinate of center of ellipse */
Parameter TiltAngle = 0; /* Rotation of ellipse in counter-clockwise direction (radians) */
Parameter Xdim; /* Radius of ellipse along X axis */
Parameter Ydim; /* Radius of ellipse along Y axis */
/*
* Work variables.
*/
double DataDistance, DataAngle, Deviation, Angle, r, EllipseX, EllipseY;
/*
* Compute the polar rotation angle and distance of this (x,y) data point
* relative to the estimated center of the ellipse.
*/
DataAngle = rtopa(x-Xcenter, y-Ycenter);
DataDistance = sqrt((x - Xcenter)^2 + (y - Ycenter)^2);
/*
* Compute the angle for the point on the ellipse with the tilt angle.
*/
Angle = DataAngle - TiltAngle;
/*
* Compute the radius of the ellipse (distance from center to perimeter) for
* this data angle. (Uses polar coordinate equation for an ellipse.)
*/
r = sqrt((Xdim^2 * Ydim^2) / ((Xdim*sin(Angle))^2 + (Ydim*cos(Angle))^2));
/*
* Compute the difference between the distance for the data point and the ellipse.
*/
Deviation = DataDistance - r;
/*
* Minimize the sum of squared deviations.
*/
Function Deviation;
/*
* X, Y data values.
*/
Data;
[ data goes here ]
In order to achieve convergence of this program, it is necessary to specify reasonable
starting values for the parameters. In particular, it is important to specify starting
values for the ellipse center (Xcenter,Ycenter) that are within the cluster of data
points (the mean X and Y values would work well).