Celestin Apprentice 2/Apprentice-Release2.iso/Source Code/Pascal/Applications/NIH Image 1.55/Source/Ellipse.p

Text File | 1993-06-17 | 8.7 KB | 331 lines | [TEXT/PJMM]

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}
{  [phone]}


{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: 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;
    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;
            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, xoffset, yoffset: extended;
            width, height: integer;
            str1, str2, str3: str255;
    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.000001;
                a11 := u02 / m4;
                a12 := u11 / m4;
                a22 := u20 / m4;
                xoffset := xm;
                yoffset := ym;
            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 / (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.