;+ ; NAME: ; Fit_Ellipse ; ; PURPOSE: ; ; This program fits an ellipse to an ROI given by a vector of ROI indices. ; ; AUTHOR: ; ; FANNING SOFTWARE CONSULTING ; David Fanning, Ph.D. ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: david@idlcoyote.com ; Coyote's Guide to IDL Programming: http://www.idlcoyote.com ; ; CATEGORY: ; ; Graphics, math. ; ; CALLING SEQUENCE: ; ; ellipsePts = Fit_Ellipse(indices) ; ; OPTIONAL INPUTS: ; ; indices - A 1D vector of pixel indices that describe the ROI. For example, ; the indices may be returned as a result of the WHERE function. ; ; OUTPUTS: ; ; ellipsePts - A 2-by-npoints array of the X and Y points that describe the ; fitted ellipse. The points are in the device coodinate system. ; ; INPUT KEYWORDS: ; ; NPOINTS - The number of points in the fitted ellipse. Set to 120 by default. ; ; SCALE - A two-element array that gives the scaling parameters for each X and Y pixel, respectively. ; Set to [1.0,1.0] by default. ; ; XSIZE - The X size of the window or array from which the ROI indices are taken. ; Set to !D.X_Size by default. ; ; YSIZE - The Y size of the window or array from which the ROI indices are taken. ; Set to !D.Y_Size by default. ; ; OUTPUT KEYWORDS: ; ; CENTER -- Set to a named variable that contains the X and Y location of the center ; of the fitted ellipse in device coordinates. ; ; ORIENTATION - Set to a named variable that contains the orientation of the major ; axis of the fitted ellipse. The direction is calculated in degrees ; counter-clockwise from the X axis. ; ; AXES - A two element array that contains the length of the major and minor ; axes of the fitted ellipse, respectively. ; ; SEMIAXES - A two element array that contains the length of the semi-major and semi-minor ; axes of the fitted ellipse, respectively. (This is simple AXES/2.) ; ; EXAMPLE: ; ; LoadCT, 0, /Silent ; image = BytArr(400, 300)+125 ; image[180:245, 125:175] = 255B ; indices = Where(image EQ 255) ; Window, XSize=400, YSize=300 ; TV, image ; PLOTS, Fit_Ellipse(indices, XSize=400, YSize=300), /Device, Color=cgColor('red') ; ; MODIFICATION HISTORY: ; ; Written by David W. Fanning, April 2002. Based on algorithms provided by Craig Markwardt ; and Wayne Landsman in his TVEllipse program. ; Added SCALE keyword and modified the algorithm to use memory more efficiently. ; I no longer have to make huge arrays. The arrays are only as big as the blob ; being fitted. 17 AUG 2008. DWF. ; Fixed small typo that caused blobs of indices with a longer X axis than Y axis ; to misrepresent the center of the ellipse. 23 February 2009. ;- ;******************************************************************************************; ; Copyright (c) 2008-2009, by Fanning Software Consulting, Inc. ; ; All rights reserved. ; ; ; ; Redistribution and use in source and binary forms, with or without ; ; modification, are permitted provided that the following conditions are met: ; ; ; ; * Redistributions of source code must retain the above copyright ; ; notice, this list of conditions and the following disclaimer. ; ; * Redistributions in binary form must reproduce the above copyright ; ; notice, this list of conditions and the following disclaimer in the ; ; documentation and/or other materials provided with the distribution. ; ; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ; ; contributors may be used to endorse or promote products derived from this ; ; software without specific prior written permission. ; ; ; ; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ; ; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ; ; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ; ; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ; ; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ; ; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ; ; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ; ; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ; ; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ; ; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ; ;******************************************************************************************; FUNCTION Fit_Ellipse, indices, $ AXES=axes, $ CENTER=center, $ NPOINTS=npoints, $ ORIENTATION=orientation, $ SEMIAXES=semiAxes, $ SCALE=scale, $ XSIZE=xsize, $ YSIZE=ysize ; The following method determines the "mass density" of the ROI and fits ; an ellipse to it. This is used to calculate the major and minor axes of ; the ellipse, as well as its orientation. The orientation is calculated in ; degrees counter-clockwise from the X axis. IF N_Elements(xsize) EQ 0 THEN xsize = !D.X_Size IF N_Elements(ysize) EQ 0 THEN ysize = !D.Y_Size IF N_Elements(npoints) EQ 0 THEN npoints = 120 IF N_Elements(scale) EQ 0 THEN scale = [1.0, 1.0] IF N_Elements(scale) EQ 1 THEN scale = [scale, scale] ; Fake indices for testing purposes. IF N_Elements(indices) EQ 0 THEN BEGIN xs = xsize / 4 xf = xsize / 4 * 2 ys = ysize / 4 yf = ysize / 4 * 2 array = BytArr(xsize, ysize) array[xs:xf, ys:yf] = 255B indices = Where(array EQ 255) ENDIF ; Convert the indices to COL/ROW coordinates. Find min and max values xyindices = Array_Indices([xsize,ysize], indices, /DIMENSIONS) minX = Min(xyindices[0,*], MAX=maxX) minY = Min(xyindices[1,*], MAX=maxY) cols = Reform(xyindices[0,*]) - minX rows = Reform(xyindices[1,*]) - minY ; Make an array large enough to hold the blob. arrayXSize = maxX-minX+1 arrayYSize = maxY-minY+1 array = BytArr(arrayXSize, arrayYSize) array[xyindices[0,*] - minX, xyindices[1,*] - minY] = 255B totalMass = Total(array) xcm = Total( Total(array, 2) * Indgen(arrayXSize) * scale[0] ) / totalMass ycm = Total( Total(array, 1) * Indgen(arrayYSize) * scale[1] ) / totalMass center = [xcm, ycm] ; Obtain the position of every pixel in the image, with the origin ; at the center of mass of the ROI. x = Findgen(arrayXSize) * scale[0] y = Findgen(arrayYSize) * scale[1] xx = (x # (y * 0 + 1)) - (xcm) yy = ((x * 0 + 1) # y) - (ycm) npts = N_Elements(indices) ; Calculate the mass distribution tensor. i11 = Total(yy[cols, rows]^2) / npts i22 = Total(xx[cols, rows]^2) / npts i12 = -Total(xx[cols, rows] * yy[cols, rows]) / npts tensor = [[ i11, i12],[i12,i22]] ; Find the eigenvalues and eigenvectors of the tensor. evals = Eigenql(tensor, Eigenvectors=evecs) ; The semi-major and semi-minor axes of the ellipse are obtained from the eigenvalues. semimajor = Sqrt(evals[0]) * 2.0 semiminor = Sqrt(evals[1]) * 2.0 ; We want the actual axes lengths. major = semimajor * 2.0 minor = semiminor * 2.0 semiAxes = [semimajor, semiminor] axes = [major, minor] ; The orientation of the ellipse is obtained from the first eigenvector. evec = evecs[*,0] ; Degrees counter-clockwise from the X axis. orientation = ATAN(evec[1], evec[0]) * 180. / !Pi - 90.0 ; Divide a circle into Npoints. 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 points to return. pts = FltArr(2, N_Elements(xprime)) pts[0,*] = xprime + minX pts[1,*] = yprime + minY center = center + [minX, minY] RETURN, pts END