;+
; 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