Fitting a plane ellipse in 3space

Directly to the input form

• The job of fitting a plane ellipse to 3D data occurs for example in computing the orbits of satellites or other celestial objects from measurements. It is a special application case of orthogonal least distance fitting. Due to the high nonlinearity of the problem the proper choice of model equations and the provision of good initial guesses for the parameters is important.
• The most obvious choice of the model considers the ellipse as the result of the intersection of a plane
  H : a₁x+a₂y+a₃z+a₄ = 0
  and a quadric
  Q : ∑ _{i,j: i+j=0,1,2,}c_i,jxⁱy^j = 0
  making the models locally unique by a normalization constraint on the coefficients. With N data points this results in an optimization problem with 3N+14 variables and 2N+2 nonlinear equality constraints (the constraints that the fitpoints lie on H and Q and the normalization of the coefficients a and c.) Unfortunately even worse, the Jacobian of these constraints will be illconditioned if the data points (and hence finally the fitpoints) will be dense in space.
• With this in mind, a completely different approach is taken here: we consider the ellipse as the image of the unit circle in the (x,y,z=0)-plane: using
  1. scaling of the x and y axis:
     x → λ₁ x and y → λ₂ y
     producing an ellipse.
  2. a rotation of this ellipse, defined by an orthonormal matrix V
     (x',y')^T = V (x,y)^T
     Since in this special case V will have a quite special form, namely
     [ α , β ; -β , α ]
     with α² + β² = 1
     this adds one more free parameter to the model.
  3. We shift the center of the ellipse from (0,0,0) to (xc,yc,0)
  4. next we rotate the 3D axis system such that the z- axis points into the direction of the intended hyperplanes normal (the vector (a₁, a₂, a₃) above with normalization to length one). This is achieved by a 3 × 3 orthonormal matrix W which has this normal as last column and the two first columns of the Householder reflector which transforms the normal to the vector (0,0,1)^T as column one and two. This adds two more degrees of freedom because of the normalization to length one.
     Warning: this reflector is not continuous as a function of the normal. Hence we may exchange two axis in the data in order to get the fitted plane flat relative to the new (x,y,z=0)-plane, i.e.
     |a₃| >= |a₁|,|a₂| !
  5. next we shift the rotated plane in the direction of the normal by the amount -a₄. This defines the plane in which the fitpoints will lie.
  6. N angles along the unit circle define the positions of the fitpoints.
  7. The objective of the minimization is the sum of squares of distances between the fitpoints and the datapoints.
  This gives a total of N+10 variables with two normalization constraints and hence is a minimal representation. There are 2 nonlinear equality constraints (the normalization of the plane normal and of the columns of V) and 2N bound constraints (lower bound 0 and upper bound 2 π for the N angles).
• The initial guesses for the parameters are obtained from the data points the other way around.
  1. first we fit in the sense of orthogonal distance regression a plane to the data points. This is achieved by a the computation of the SVD decomposition of the matrix with the rows [x_i,y_i,z_i,1] (see also this application here in this chapter). The parameters of this plane in Hesse normal form define the inputs for the normal and the plane shift.
  2. next we shift the plane by the amount a₄ in the direction of the normal. Now the point (0,0,0) is an element of the shifted plane and we rotate it into the standard coordinate system using the matrix W defined above in transposed form.
  3. The data are dealt with in the same manner. We project them now orthogonally to the (x,y,z=0) plane simply by neglecting the third component.
  4. Now we fit an ellipse to these plane data, again using the SVD, this time using the matrix with rows
     [x_i², x_iy _i, y_i²,x_i,y _i,1]
     obtaining the model of a second order curve (hopefully an ellipse) in the plane.
  5. next we rewrite this in the form
     (x-xc,y-yc) A (x-xc,y-yc)^T = 1
     checking of course the positive definiteness of the 2 × 2 matrix A
  6. Finally we decompose with orthonormal V and diagonal Λ
     A = V Λ V^T ,
     getting intial guesses for the parameter λ₁, λ₂, α , β .
  7. The code behind is a f77-translation of svd.f from linpack and the authors optimization code donlp2.

Input

• You may make an experiment, specifying the parameters for generating the data points and disturbing these subsequently, in order to get an impression about the sensitivity of the problem or
• You may fit data of your own.
• In the case of the experiment, you define the number of data points, the parameter settings (normalization of the parameters is done internally) and the amount of perturbation. This is done componentwise using a formula
  value_disturbed = value_true × (1+2(0.5-rand())*errlev)
  rand() is a quasi random number generator with values in [0,1].
• In the case of data of your own the number of points N, a table with these points (simply a sequence of 3N numbers)
• For both cases you finally must specify two angles for your view point for the 3D plot of curve and data.

Output

• A text file with information about the fits in computing the initial guesses. Maybe simply the message that your data could not be fitted by an ellipse. The fit error (in the euclidean norm) is simply the smallest singular value of the fit matrix (since the system is homogeneous). Since an absolute value has no information, we provide the quotient with the largest singular value. This is meant by ''norm(relative)''. This is followed by the results of the optimizer. The order of the variables in the initial and final value is:
  1. the N angles of the pictures of the data points and the fitpoints (measured on the unit circle!)
  2. the eigenvalues of the matrix A. The reciprocals of their square roots are the half axes of the projected ellipse in the (x,y,0)-plane
  3. the first column of the matrix V (rotation in the (x,y,0)-plane)
  4. the center of the projected ellipse in the (x,y,0)-plane
  5. the four parameters of the Hesse normal form of the plane
  fx is the sum of squares of the final distances of the data points to the ellipse (in 3D).
  Then there is the termination message of the optimizer, the values of the error measures (see the desciption of donlp2 in the section constrained optimization) and the number of function and gradient evaluations.
• There is another form of representation of the ellipse, namely the two half axes, three Euler angles for defining the total rotation, i.e. the factorization of the product
  W * [ V, (0,0)^T ; 0 , 0 , 1 ] = R₁R₂R₃
  into three plane rotations and the shift of (0,0,0) to the ellipse center. This makes the generation of the ellipse as curve for plotting easier and we provide you with the program text.
• A 3D plot showing your data points and the ellipse. Sometimes this plot may look strange: carefully look on the scales of the three axis: these are sometimes very different!

Questions ?!

• What can you observe about the sensitivity of the solution if you restrict the choice of the angles or require a large errlev?.
• Why do we need at least 8 data points?
• If you make the quotient λ₁/λ₂ very large or very small and errlev moderate large, what might occur and why?
• We require in the input form λ₁ and λ₂ not less 0.001 and the quotient of the two less or equal 1000. This has a practical reason. Which one?

to the input form: 2D3D ellipse fit