EllipseFitter.java

package ij.process;
import ij.*;
import ij.gui.*;
import java.awt.*;
import ij.plugin.filter.*;

/*
Best-fitting ellipse routines by:

  Bob Rodieck
  Department of Ophthalmology, RJ-10
  University of Washington, 
  Seattle, WA, 98195

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. 

public class Ellipse_Fitter implements PlugInFilter {
    public int setup(String arg, ImagePlus imp) {
        return DOES_ALL;
    }
    public void run(ImageProcessor ip) {
        EllipseFitter ef = new EllipseFitter();
        ef.fit(ip);
        IJ.write(IJ.d2s(ef.major)+" "+IJ.d2s(ef.minor)+" "+IJ.d2s(ef.angle)+" "+IJ.d2s(ef.xCenter)+" "+IJ.d2s(ef.yCenter));
        ef.drawEllipse(ip);
    }
}
*/


public class EllipseFitter {

    static final double HALFPI = 1.5707963267949;
    
    public double xCenter;

    public double  yCenter;
    
    public double major;
    
    public double minor;
    
    public double angle;
    
    public double theta;
    
    public int[] xCoordinates;
    public int[] yCoordinates;
    public int nCoordinates = 0;

    
    private int bitCount;
    private double  xsum, ysum, x2sum, y2sum, xysum;
    private byte[] mask;
    private int left, top, width, height;
    private double   n;
    private double   xm, ym;   //mean values
    private double   u20, u02, u11;  //central moments
    private ImageProcessor ip;
    private double pw, ph;
    private boolean record;

    public void fit(ImageProcessor ip, ImageStatistics stats) {
        this.ip = ip;
        mask = ip.getMaskArray();
        Rectangle r = ip.getRoi();
        this.pw = stats.pw;
        this.ph = stats.ph;
        left = r.x;
        top = r.y;
        width = r.width;
        height = r.height;
        getEllipseParam();
    }
    
    void getEllipseParam() {
        double    sqrtPi = 1.772453851;
        double    a11, a12, a22, m4, z, scale, tmp, xoffset, yoffset;
        double    RealAngle;

        if (mask==null) {
            major = (width*2) / sqrtPi;
            minor = (height*2) / sqrtPi; // * Info->PixelAspectRatio;
            angle = 0.0;
            theta = 0.0;
            if (major < minor) {
                tmp = major;
                major = minor;
                minor = tmp;
                angle = 90.0;
                theta = Math.PI/2.0;
            }
            xCenter = left + width / 2.0;
            yCenter = top + height / 2.0;
            return;
        }

        computeSums();
        getMoments();
        m4 = 4.0 * Math.abs(u02 * u20 - u11 * u11);
        if (m4 < 0.000001)
            m4 = 0.000001;
        a11 = u02 / m4;
        a12 = u11 / m4;
        a22 = u20 / m4;
        xoffset = xm;
        yoffset = ym;

        tmp = a11 - a22;
        if (tmp == 0.0)
            tmp = 0.000001;
        theta = 0.5 * Math.atan(2.0 * a12 / tmp);
        if (theta < 0.0)
            theta += HALFPI;
        if (a12 > 0.0)
            theta += HALFPI;
        else if (a12 == 0.0) {
            if (a22 > a11) {
                theta = 0.0;
                tmp = a22;
                a22 = a11;
                a11 = tmp;
            } else if (a11 != a22)
                theta = HALFPI;
        }
        tmp = Math.sin(theta);
        if (tmp == 0.0)
            tmp = 0.000001;
        z = a12 * Math.cos(theta) / tmp;
        major = Math.sqrt (1.0 / Math.abs(a22 + z));
        minor = Math.sqrt (1.0 / Math.abs(a11 - z));
        scale = Math.sqrt (bitCount / (Math.PI * major * minor)); //equalize areas
        major = major*scale*2.0;
        minor = minor*scale*2.0;
        angle = 180.0 * theta / Math.PI;
        if (angle == 180.0)
            angle = 0.0;
        if (major < minor) {
            tmp = major;
            major = minor;
            minor = tmp;
        }
        xCenter = left + xoffset + 0.5;
        yCenter = top + yoffset + 0.5;
    }

    void computeSums () {
        xsum = 0.0;
        ysum = 0.0;
        x2sum = 0.0;
        y2sum = 0.0;
        xysum = 0.0;
        int bitcountOfLine;
        double   xe, ye;
        int xSumOfLine;
        for (int y=0; y<height; y++) {
            bitcountOfLine = 0;
            xSumOfLine = 0;
            int offset = y*width;
            for (int x=0; x<width; x++) {
                if (mask[offset+x] != 0) {
                    bitcountOfLine++;
                    xSumOfLine += x;
                    x2sum += x * x;
                }
            } 
            xsum += xSumOfLine;
            ysum += bitcountOfLine * y;
            ye = y;
            xe = xSumOfLine;
            xysum += xe*ye;
            y2sum += ye*ye*bitcountOfLine;
            bitCount += bitcountOfLine;
        }
    }

    void getMoments () {
        double   x1, y1, x2, y2, xy;

        if (bitCount == 0)
            return;

        x2sum += 0.08333333 * bitCount;
        y2sum += 0.08333333 * bitCount;
        n = bitCount;
        x1 = xsum/n;
        y1 = ysum / n;
        x2 = x2sum / n;
        y2 = y2sum / n;
        xy = xysum / n;
        xm = x1;
        ym = y1;
        u20 = x2 - (x1 * x1);
        u02 = y2 - (y1 * y1);
        u11 = xy - x1 * y1;
    }
    
    /* 
    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    
    */
    
    public void drawEllipse(ImageProcessor ip) {
        if (major==0.0 && minor==0.0)
            return;
        int xc = (int)Math.round(xCenter);
        int yc = (int)Math.round(yCenter);
        int maxY = ip.getHeight();
        int xmin, xmax;
        double sint, cost, rmajor2, rminor2, g11, g12, g22, k1, k2, k3;
        int x, xsave, ymin, ymax;
        int[] txmin = new int[maxY];
        int[] txmax = new int[maxY];
        double j1, j2, yr;

        sint = Math.sin(theta);
        cost = Math.cos(theta);
        rmajor2 = 1.0 / sqr(major/2);
        rminor2 = 1.0 / sqr(minor/2);
        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 = (int)Math.floor(Math.sqrt(Math.abs(k3 / k2)));
        if (ymax>maxY)
            ymax = maxY;
        if (ymax<1)
            ymax = 1;
        ymin = -ymax;
        // Precalculation and use of symmetry speed things up
        for (int y=0; y<=ymax; y++) {
            //GetMinMax(y, aMinMax);
            j2 = Math.sqrt(k2 * sqr(y) + k3);
            j1 = k1 * y;
            txmin[y] = (int)Math.round(j1 - j2);
            txmax[y] = (int)Math.round(j1 + j2);
        }
        if (record) {
            xCoordinates[nCoordinates] = xc + txmin[ymax - 1];
            yCoordinates[nCoordinates] = yc + ymin;
            nCoordinates++;
        } else
            ip.moveTo(xc + txmin[ymax - 1], yc + ymin);
        for (int y=ymin; y<ymax; y++) {
            x = y<0?txmax[-y]:-txmin[y];
            if (record) {
                xCoordinates[nCoordinates] = xc + x;
                yCoordinates[nCoordinates] = yc + y;
                nCoordinates++;
            } else
                ip.lineTo(xc + x, yc + y);
        }
        for (int y=ymax; y>ymin; y--) {
            x = y<0?txmin[-y]:-txmax[y];
            if (record) {
                xCoordinates[nCoordinates] = xc + x;
                yCoordinates[nCoordinates] = yc + y;
                nCoordinates++;
            } else
                ip.lineTo(xc + x, yc + y);
        }
    }
    
    public void makeRoi(ImageProcessor ip) {
        record = true;
        int size = ip.getHeight()*3;
        xCoordinates = new int[size];
        yCoordinates = new int[size];
        nCoordinates = 0;
        drawEllipse(ip);
        record = false;
    }
    
    private double sqr(double x) {
        return x*x;
    }

}