File: fit_ellipse.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (220 lines) | stat: -rw-r--r-- 9,088 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
;+
; 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 > 3
    arrayYSize = maxY-minY+1 > 3
    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