unit Ellipse; interface uses Types, QuickDraw, Palettes, globals, Utilities, graphics; procedure DrawEllipse; procedure GetEllipseParam (var Major, Minor, angle, xxcenter, yycenter: extended); procedure ComputeSums (y, width: integer; var MaskLine: LineType); procedure ResetSums; {Best-fitting ellipse routines by:} { Bob Rodieck} { Department of Ophthalmology, RJ-10} { University of Washington, } { Seattle, WA, 98195} { (206) 543-2449} {Notes on best-fitting ellipse:} { Consider some arbitrarily shaped closed profile, which we wish to} { characterize in a quantitative manner by a series of terms, each } { term providing a better approximation to the shape of the profile. } { Assume also that we wish to include the orientation of the profile } { (i.e. which way is up) in our characterization. } { One approach is to view the profile as formed by a series harmonic } { components, much in the same way that one can decompose a waveform} { over a fixed interval into a series of Fourier harmonics over that } { interval. From this perspective the first term is the mean radius,} { or some related value (i.e. the area). The second term is the } { magnitude and phase of the first harmonic, which is equivalent to the} { best-fitting ellipse. } { What constitutes the best-fitting ellipse? First, it should have the} { same area. In statistics, the measure that attempts to characterize some} { two-dimensional distribution of data points is the 'ellipse of } { concentration' (see Cramer, Mathematical Methods of Statistics, } { Princeton Univ. Press, 945, page 283). This measure equates the second} { order central moments of the ellipse to those of the distribution, } { and thereby effectively defines both the shape and size of the ellipse. } { This technique can be applied to a profile by assuming that it constitutes} { a uniform distribution of points bounded by the perimeter of the profile.} { For most 'blob-like' shapes the area of the ellipse is close to that} { of the profile, differing by no more than about 4%. We can then make} { a small adjustment to the size of the ellipse, so as to give it the } { same area as that of the profile. This is what is done here, and } { therefore this is what we mean by 'best-fitting'. } { For a real pathologic case, consider a dumbell shape formed by two small} { circles separated by a thin line. Changing the distance between the} { circles alters the second order moments, and thus the size of the ellipse } { of concentration, without altering the area of the profile. } implementation const HalfPi = 1.5707963267949; type TMoments= record n: extended; xm, ym, { mean values } u20, u02, u11: extended; { central moments } end; var BitCount, xsum, ysum: LongInt; x2sum, y2sum, xysum: extended; Moments: TMoments; gMajor, gMinor, Theta: extended; gxCenter, gyCenter: integer; SaveRect: rect; procedure DrawEllipse; { basic equations: } { a: major axis} { b: minor axis} { t: theta, angle of major axis, clockwise with respect to x axis. } { g11*x^2 + 2*g12*x*y + g22*y^2 = 1 -- equation of ellipse} { g11:= ([cos(t)]/a)^2 + ([sin(t)]/b)^2} { g12:= (1/a^2 - 1/b^2) * sin(t) * cos(t)} { g22:= ([sin(t)]/a)^2 + ([cos(t)]/b)^2} { solving for x: x:= k1*y ± sqrt( k2*y^2 + k3 )} { where: k1:= -g12/g11} { k2:= (g12^2 - g11*g22)/g11^2} { k3:= 1/g11} { ymax or ymin occur when there is a single value for x, that is when: } { k2*y^2 + k3 = 0 } const maxY = 1000; type TMinMax = record xmin, xmax: Integer; end; var sint, cost, rmajor2, rminor2, g11, g12, g22, k1, k2, k3: extended; xsave, y, ymin, ymax: Integer; aMinMax: TMinMax; TXList: array[0..maxY] of TMinMax; procedure GetMinMax (yValue: Integer; var xMinMax: TMinMax); var j1, j2, yr: extended; begin yr := yvalue; j2 := sqrt(k2 * sqr(yr) + k3); j1 := k1 * yr; with xMinMax do begin xmin := round(j1 - j2); xmax := round(j1 + j2); end; end; procedure Plot (x: Integer); begin MoveTo(gxCenter + xsave, gyCenter + y); LineTo(gxCenter + x, gyCenter + y); xsave := x; end; begin if not EqualRect(info^.RoiRect, SaveRect) then exit(DrawEllipse); sint := sin(Theta); cost := cos(Theta); rmajor2 := 1.0 / sqr(gMajor); rminor2 := 1.0 / sqr(gMinor); g11 := rmajor2 * sqr(cost) + rminor2 * sqr(sint); g12 := (rmajor2 - rminor2) * sint * cost; g22 := rmajor2 * sqr(sint) + rminor2 * sqr(cost); k1 := -g12 / g11; k2 := (sqr(g12) - g11 * g22) / sqr(g11); k3 := 1.0 / g11; ymax := Trunc(sqrt(abs(k3 / k2))); if ymax > maxy then ymax := maxy; ymin := -ymax; { Precalculation and use of symmetry speed things up } for y := 0 to ymax do begin GetMinMax(y, aMinMax); TXList[y] := aMinMax; end; xsave := TXList[ymax - 1].xmin; { i.e. abs(ymin+1) } for y := ymin to ymax - 1 do with TXList[abs(y)] do if y < 0 then Plot(xmax) else Plot(-xmin); for y := ymax downto ymin + 1 do with TXList[abs(y)] do if y < 0 then Plot(xmin) else Plot(-xmax); end; { TraceOval } procedure GetMoments;{changed n_} var x1, y1, x2, y2, xy: extended; begin with moments, Info^ do begin if BitCount = 0 then exit(GetMoments); x2sum := x2sum + 0.08333333* BitCount; {NIntegrate[x^2, ]-center^2 = 0.08333333} y2sum := y2sum + 0.08333333* BitCount; {=correction when the mass of a pixel is seen as an area instead of a point} n := bitcount; x1 := xsum / n; y1 := ysum * PixelAspectRatio/ n; x2 := x2sum / n; y2 := y2sum * sqr(PixelAspectRatio)/ n; xy := xysum * PixelAspectRatio / n; xm := x1; ym := y1; u20 := x2 - sqr(x1); u02 := y2 - sqr(y1); u11 := xy - x1 * y1; end; end; procedure GetEllipseParam (var Major, Minor, angle, xxcenter, yycenter: extended);{changed n_} {Return the parameters of an ellipse that has the same second-order } {moments as those specified by 'm'. } {See Cramer, Mathematical Methods of Statistics, } {Princeton Univ. Press 1945, page 283.} {The elliptical parameters returned specify an ellipse that} {has the the same second order moments as that of the profile} {that generated the moments. This ellipse need not have the same} {area as that of the profile, although its area will be close to} {that of the profile. In order to refine our measure, we scale} {the major and minor axes so as to make the area equal to that} {of the profile. } const sqrtPi = 1.772453851; var a11, a12, a22, m4, z, scale, tmp, xoffset, yoffset: extended; width, height: integer; str1, str2, str3: str255; RealAngle: real; begin with info^, info^.RoiRect do begin width := right - left; height := bottom - top; if RoiType = RectRoi then begin major := width / sqrtPi; minor := height / sqrtPi * PixelAspectRatio; angle := 0.0; if major < minor then begin tmp := major; major := minor; minor := tmp; angle := 90.0; end; xxcenter := left - 0.5 + width / 2.0; yycenter := top - 0.5 + height / 2.0; exit(GetEllipseParam); end; end; GetMoments; with moments do begin m4 := 4.0 * abs(u02 * u20 - sqr(u11)); if m4 <0.000001 then m4 := 0.000001; a11 := u02 / m4; a12 := u11 / m4; a22 := u20 / m4; xoffset := xm; yoffset := ym/Info^.PixelAspectRatio; end; tmp := a11 - a22; if tmp = 0.0 then tmp := 0.000001; theta := 0.5 * arctan(2.0 * a12 / tmp); if theta < 0.0 then theta := theta + halfpi; if a12 > 0.0 then theta := theta + halfpi else if a12 = 0 then begin if a22 > a11 then begin theta := 0; tmp := a22; a22 := a11; a11 := tmp; end else if a11 <> a22 then theta := halfpi; end; tmp := sin(theta); if tmp = 0.0 then tmp := 0.000001; z := a12 * cos(theta) / tmp; major := sqrt(1.0 / abs(a22 + z)); minor := sqrt(1.0 / abs(a11 - z)); scale := sqrt(BitCount * Info^.PixelAspectRatio/ (pi * major * minor )); {equalize areas } major := major * scale; minor := minor * scale; RealAngle := 180.0 * theta / pi;{force rounding by using real} angle := RealAngle; if angle = 180.0 then angle := 0.0; if major < minor then begin tmp := major; major := minor; minor := tmp; end; with info^ do begin with RoiRect do begin xxCenter := left + xoffset; yyCenter := top + yoffset; end; SaveRect := RoiRect; end; gxCenter := round(xxCenter); gyCenter := round(yyCenter); gMajor := major; gMinor := minor; end; procedure ComputeSums (y, width: integer; var MaskLine: LineType);{changed n_} var x: longint; BitcountOfLine: longint; xe, ye: extended; xSumOfLine: longint; begin BitcountOfLine := 0; xSumOfLine:= 0; for x := 0 to width - 1 do if MaskLine[x] = BlackIndex then begin BitcountOfLine := BitcountOfLine+1; xSumOfLine := xSumOfLine + x; x2sum := x2sum + sqr(x); end; xsum := xsum + xSumOfLine; ysum := ysum + BitcountOfLine * y; ye := y; xe := xSumOfLine; xysum := xysum + xe * ye; y2sum := y2sum + sqr(ye) * BitcountOfLine; bitCount := bitCount + BitcountOfLine; end; procedure ResetSums; begin xsum := 0; ysum := 0; x2sum := 0.0; y2sum := 0.0; xysum := 0.0; bitCount := 0; end; end.