Fit an Ellipse to a Region of Interest

QUESTION: I have segmented my image into a set of distinct regions of interest. (Well, irregular blobs, really.) I would like to try to characterize the shape of these blobs in some way. I thought one way would be to fit an ellipse to each blob and use the ellipse parameters as substitutes for general shape parameters. But I can't figure out how to fit an ellipse to these irregular shapes. Do you have any ideas?

ANSWER: Not really. But when I have really tough questions I usually ask one of the regulars on the IDL newsgroup (comp.lang.idl-pvwave). They have all the answers. :-)

In this case, since this seems like a math question, I ran to Craig Markwardt (craigmnet@cow.physics.wisc.edu), the math guy. He provided me with some equations that I combined with a neat ellipse display program written by Wayne Landsman for the NASA Goddard IDL program library to come up with a "center of mass" ellipse fitting program, named Fit_Ellipse.

Here is a short (and probably inaccurate, because I don't really understand the math) explanation for how it works. (I'll be honest, I'm less concerned with how it works than I am with the fact that it does work.)

The method I am using here starts with the indices of the ROI and is a "mass density" method, in which each pixel has a "weight" that can be used to determine the center of mass and mass distribution of the ROI. The ellipse is fitted to this mass distribution.

The first step is create an array or image with just the pixels in the ROI having positive values. This array is used to calculate the center of mass of the ROI. In the code below, the variable indices contains a 1D vector of indices that describes the ROI. The variables xsize and ysize are the size of the window or array that the ROI came from. The code looks like this:

   array = BytArr(xsize, ysize)
   array[indices] = 255B
   totalMass = Total(array)
   xcm = Total( Total(array, 2) * Indgen(xsize) ) / totalMass
   ycm = Total( Total(array, 1) * Indgen(ysize) ) / totalMass
   center = [xcm, ycm]

The next step is to locate each ROI pixel in the array with respect to the center of mass (the variable center, above). In other words, we want to make the center of our coordinate system be the center of mass. The code looks like this:

   x = Findgen(xsize)
   y = Findgen(ysize)
   xx = (x # (y * 0 + 1)) - xcm
   yy = ((x * 0 + 1) # y) - ycm

Next, we calculate the mass distribution tensor for these pixels.

  npts = N_Elements(indices)
  i11 = Total(yy[indices]^2) / npts
  i22 = Total(xx[indices]^2) / npts
  i12 = -Total(xx[indices] * yy[indices]) / npts
  tensor = [ [i11, i12], [i12,i22] ]

We use the tensor to calculate the eigenvalues and eigenvectors. The EIGENQL command is used.

   evals = Eigenql(tensor, Eigenvectors=evecs)

By some trick of the sorcerer's art, the eigenvalues lead us directly to the semi-major and semi-minor axes of the ellipse we are trying to find.

   semimajor = Sqrt(evals[0]) * 2.0
   semiminor = Sqrt(evals[1]) * 2.0

But, to draw the ellipse we need the major and minor axes. Easy. (Even I knew how to do this part.)

   major = semimajor * 2.0
   minor = semiminor * 2.0
   semiAxes = [semimajor, semiminor]
   axes = [major, minor]

The orientation of the ellipse is found from the first eigenvector. The orientation is calculated in degrees counter-clockwise from the X axis.

   evec = evecs[*,0]
   orientation = ATAN(evec[1], evec[0]) * 180. / !Pi - 90.0

Finally, we calculate the points that make up the ellipse itself. In practice you might need 100-150 points for a smooth ellipse. Here I use 120.

   npoints = 120
   phi = 2 * !PI * (Findgen(npoints) / (npoints-1))

   ; Position angle in radians.

   ang = orientation / !RADEG

   ; Sin and cos of angle.

   cosang = Cos(ang)
   sinang = Sin(ang)

   ; Parameterized equation of ellipse.

   x =  semimajor * Cos(phi)
   y =  semiminor * Sin(phi)

    ; Rotate to desired position angle.

   xprime = xcm + (x * cosang) - (y * sinang)
   yprime = ycm + (x * sinang) + (y * cosang)

   ; Extract the ellipse points.

   pts = FltArr(2, N_Elements(xprime))
   pts[0,*] = xprime
   pts[1,*] = yprime

All that is left is to draw the ellipse.

   IDL> Plots, pts, /Device

You see the result in the figure below.

The Fit_Ellipse program has been incorporated into the larger Blob_Analyzer program to make ROI or blob analysis even easier.