# 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

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.

Copyright © 2003 David W. Fanning

Last Updated 3 December 2003