unit Ellipse; interface uses QuickDraw, Palettes, PrintTraps, 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.570796; 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: Real; xsave, y, ymin, ymax: Integer; aMinMax: TMinMax; TXList: array[0..maxY] of TMinMax; procedure GetMinMax (yValue: Integer; var xMinMax: TMinMax); var j1, j2, yr: Real; 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; var x1, y1, x2, y2, xy: extended; str1, str2, str3: str255; begin with moments do begin if BitCount = 0 then exit(GetMoments); n := bitcount; x1 := xsum / n; y1 := ysum / n; x2 := x2sum / n; y2 := y2sum / n; xy := xysum / n; xm := x1; ym := y1; u20 := x2 - sqr(x1); u02 := y2 - sqr(y1); u11 := xy - x1 * y1; exit(GetMoments); RealToString(u20, 8, 2, str1); RealToString(u02, 8, 2, str2); RealToString(u11, 8, 2, str3); ShowMessage(concat('u20=', str1, cr, 'u02=', str2, cr, 'u11=', str3)); wait(300); end; end; procedure GetEllipseParam (var Major, Minor, angle, xxcenter, yycenter: extended); {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. } var a11, a12, a22, m4, z, scale, tmp: extended; xoffset, yoffset, width, height: integer; begin with info^, info^.RoiRect do begin if PixelCount^[mCount] = 1 then begin major := 0.564; Minor := 0.564; angle := 0.0; xxcenter := left + 0.5; yycenter := top + 0.5; exit(GetEllipseParam); end; width := right - left; height := bottom - top; if width = 1 then begin major := 0.6 * height; Minor := 0.564; angle := 90.0; xxcenter := left + 0.5; yycenter := top + height / 2.0; exit(GetEllipseParam); end; if height = 1 then begin major := 0.6 * width; Minor := 0.564; angle := 0.0; xxcenter := left + width / 2.0; yycenter := top + 0.5; exit(GetEllipseParam); end; end; {with} GetMoments; with moments do begin m4 := 4.0 * abs(u02 * u20 - sqr(u11)); if m4 = 0.0 then m4 := 0.001; a11 := u02 / m4; a12 := u11 / m4; a22 := u20 / m4; xoffset := round(xm); yoffset := round(ym); end; if (a11 - a22) <> 0 then theta := 0.5 * arctan(2.0 * a12 / (a11 - a22)) else theta := 0.0; 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; if sin(theta) <> 0.0 then z := a12 * cos(theta) / sin(theta) else z := 0.0; major := sqrt(1.0 / abs(a22 + z)); minor := sqrt(1.0 / abs(a11 - z)); scale := sqrt(BitCount / (pi * major * minor)); { equalize areas } major := major * scale; minor := minor * scale; angle := 180.0 * theta / pi; 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); var x: integer; xe, ye: extended; begin for x := 0 to width - 1 do if MaskLine[x] = BlackIndex then begin xsum := xsum + x; ysum := ysum + y; xe := x; ye := y; x2sum := x2sum + xe * xe; y2sum := y2sum + ye * ye; xysum := xysum + xe * ye; bitCount := bitCount + 1; end; end; procedure ResetSums; begin xsum := 0; ysum := 0; x2sum := 0.0; y2sum := 0.0; xysum := 0.0; bitCount := 0; end; end.