unit Ellipse;

interface

	uses
		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, <x, centerx-0.5, centerx+0.5>]-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.