unit Ellipse;

interface

	uses
		QuickDraw, OSIntf, PickerIntf, PrintTraps, ToolIntf, 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^.osroiRect, 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^.osroiRect, MeasurementsP^ do begin
				if PixelCount[nRegions] = 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 osRoiRect do begin
						xxCenter := left + xoffset;
						yyCenter := top + yoffset;
					end;
				SaveRect := osroiRect;
			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.