ellifit.java

It help you to find the center of an ellipse

//**************************************************************************
// This java code is  an interactive demo of the first ellipse-specific
// direct fitting method presented in the papers:

//    M. Pilu, A. Fitzgibbon, R.Fisher ``Ellipse-specific Direct
//    least-square Fitting '' , IEEE International Conference on Image
//    Processing, Lausanne, September 1996. (poscript) (HTML) 

//    A. Fitzgibbon, M. Pilu , R.Fisher ``Direct least-square fitting of
//    Ellipses '' , International Conference on Pattern Recognition, Vienna,
//    August 1996. (poscript) - Extended version available as DAI Research
//    Paper #794

// The demo can be tried out at   
//   http://www.dai.ed.ac.uk/students/maurizp/ElliFitDemo/demo.html

// The code was written by Maurizio Pilu , University of Edinburgh.  The
// applet's graphic interface was much inspired by the Curve Applet
// written by Michael Heinrichs at SFU, Vancouver:
// http://fas.sfu.ca:80/1/cs/people/GradStudents/heinrica/personal/curve.html 

// Some math routines are from the "Numerical Recipes in C" by
// Press/Teukolsky/Vettering/Flannery, Cambridge Uiniversity Press,
// Second Edition (1988). PLEASE READ COPYRIGHT ISSUES ON THE NUMERICAL
// RECIPES BOOK.

// NOTE: Some parts of the program are rather scruffy. The author
//       will tidy it up whan he has some spare time.

// DISCLAIMER: The authors and the department assume no responsabilities
//             whatsoever for any wrong use of this code. 

// COPYRIGHT: Any commercial use of the code and the method is forbidden
//            without written authorization from the authors.
//**************************************************************************


import java.awt.*;
import java.applet.*;
import java.lang.Double;
import java.lang.Math;
import java.util.Vector;



public class ElliFit extends Applet {
    public void init() {
	setLayout(new BorderLayout());
	CurvePanel dp = new CurvePanel();
	add("Center", dp);
	add("South",new CurveControls(dp));
	add("North",new CurveControls2(dp));
    }

    public boolean handleEvent(Event e) {
	switch (e.id) {
	  case Event.WINDOW_DESTROY:
	    System.exit(0);
	    return true;
	  default:
	    return false;
	}
    }

    public static void main(String args[]) {
	Frame f = new Frame("ElliFIt");
	ElliFit curve = new ElliFit();
	curve.init();
	curve.start();

	f.add("Center", curve);
	f.show();
    }
}




class ControlPoint extends Object {
    public int x;
    public int y;
    public static final int PT_SIZE = 2;

    public ControlPoint(int a, int b) {
	x = a;
	y = b;
    }

    public boolean within(int a, int b) {
	if (a >= x - PT_SIZE && 
	    b >= y - PT_SIZE &&
	    a <= x + PT_SIZE && 
	    b <= y + PT_SIZE)
	    return true;
	else
	    return false;
    }
}



class CurvePanel extends Panel {
    public static final int BOOKSTEIN = 0;
    public static final int TAUBIN = 1;
    public static final int FPF = 2;
    private int	   mode = BOOKSTEIN;

    public static final int ADD = 0;
    public static final int MOVE = 1;
    public static final int DELETE = 2;
    private int	  action = ADD;
    private int   i;

    public String warn_taub_str = "Sorry! Taubin method not yet implemented.";

    private Vector points = new Vector(16,4);

    // If a control point is being moved, this is the index into the list
    // of the moving point.  Otherwise it contains -1
    private int moving_point;
    private int precision;
    private static double BooksteinConstraint[][] = new double[7][7];
    private static double FPFConstraint[][] = new double[7][7];
    private static double TaubinConstraint[][] = new double[7][7];

    public CurvePanel() {
	setBackground(Color.red);
    }

    private void initConstraint() {
      // INitialize Constraint matrices
      // Here Initialization $$$

      }
   
    public void setAction(int action) {
    // Change the action type
	switch (action) {
	  case ADD:
	  case MOVE:
	  case DELETE:
	    this.action = action;
	    break;
	  default:
	    throw new IllegalArgumentException();
	}
    }

    public void setCurveType(int mode) {
    // Change the curve display type
	switch (mode) {
	  case BOOKSTEIN:
	  case TAUBIN:
	  case FPF:
	    this.mode = mode;
	    break;
	  default:
	    throw new IllegalArgumentException();
	}
    }

    public void clearPoints() {
	points.removeAllElements();
    }

    private int findPoint(int a, int b) {
    // Scan the list of control points to find out which (if any) point
    // contains the coordinates: a,b.
    // If a point is found, return the point's index, otherwise return -1
        int max = points.size();

	for(int i = 0; i < max; i++) {
	    ControlPoint pnt = (ControlPoint)points.elementAt(i);
	    if (pnt.within(a,b)) {
		return i;
	    }
	}
	return -1;
    }

    public boolean handleEvent(Event e) {
	switch (e.id) {
	  case Event.MOUSE_DOWN:
	    // How we handle a MOUSE_DOWN depends on the action mode
	    switch (action) {
	      case ADD:
		// Add a new control point at the specified location
		ControlPoint pnt;
		points.addElement(pnt = new ControlPoint(e.x, e.y));
		repaint();
		break;
	      case MOVE:
		// Attempt to select the point at the location specified.
		// If there is no point at the location, findPoint returns
		// -1 (i.e. there is no point to be moved)
		moving_point = findPoint(e.x, e.y);
		break;
	      case DELETE:
		// Delete a point if one has been clicked
		int delete_pt = findPoint(e.x, e.y);
		if(delete_pt >= 0) {
		   points.removeElementAt(delete_pt);
		   repaint();
		}
		break;
	      default:
	        throw new IllegalArgumentException();
	    }
	    return true;
	  case Event.MOUSE_UP:
	    // We only care about MOUSE_UP's if we've been moving a control
	    // point.  If so, drop the control point.
	    if (moving_point >=0) {
		moving_point = -1;
	        repaint();
	    }
	    return true;
	  case Event.MOUSE_DRAG:
	    // We only care about MOUSE_DRAG's while we are moving a control
	    // point.  Otherwise, do nothing.
	    if (moving_point >=0) {
		ControlPoint pnt = (ControlPoint) points.elementAt(moving_point);
		pnt.x = e.x;
		pnt.y = e.y;
		repaint();
	    }
	    return true;
	  case Event.WINDOW_DESTROY:
	    System.exit(0);
	    return true;
	  default:
	    return false;
	}
    }

    private void multMatrix(double m[][], double g[][], double mg[][]) {
    // This function performs the meat of the calculations for the
    // curve plotting.  Note that it is not a matrix multiplier in the
    // pure sense.  The first matrix is the curve matrix (each curve type
    // has its own matrix), and the second matrix is the geometry matrix
    // (defined by the control points).  The result is returned in the
    // third matrix.

	// First clear the return array
	for(int i=0; i<4; i++) 
	    for(int j=0; j<2; j++) 
		mg[i][j]=0;

	// Perform the matrix math
	for(int i=0; i<4; i++) 
	    for(int j=0; j<2; j++) 
		for(int k=0; k<4; k++) 
		    mg[i][j]=mg[i][j] + (m[i][k] * g[k][j]);
    }


    private void ROTATE(double a[][], int i, int j, int k, int l, double tau, double s) 
      {
	double g,h;
	g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);
	a[k][l]=h+s*(g-h*tau);
      }
    
    private void jacobi(double a[][], int n, double d[] , double v[][], int nrot)      
      {
	int j,iq,ip,i;
	double tresh,theta,tau,t,sm,s,h,g,c;

	double b[] = new double[n+1];
	double z[] = new double[n+1];
	
	for (ip=1;ip<=n;ip++) {
	  for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;
	  v[ip][ip]=1.0;
	}
	for (ip=1;ip<=n;ip++) {
	  b[ip]=d[ip]=a[ip][ip];
	  z[ip]=0.0;
	}
	nrot=0;
	for (i=1;i<=50;i++) {
	  sm=0.0;
	  for (ip=1;ip<=n-1;ip++) {
	    for (iq=ip+1;iq<=n;iq++)
	      sm += Math.abs(a[ip][iq]);
	  }
	  if (sm == 0.0) {
	    /*    free_vector(z,1,n);
		  free_vector(b,1,n);  */
	    return;
	  }
	  if (i < 4)
	    tresh=0.2*sm/(n*n);
	  else
	    tresh=0.0;
	  for (ip=1;ip<=n-1;ip++) {
	    for (iq=ip+1;iq<=n;iq++) {
	      g=100.0*Math.abs(a[ip][iq]);
	      if (i > 4 && Math.abs(d[ip])+g == Math.abs(d[ip])
		  && Math.abs(d[iq])+g == Math.abs(d[iq]))
		a[ip][iq]=0.0;
	      else if (Math.abs(a[ip][iq]) > tresh) {
		h=d[iq]-d[ip];
		if (Math.abs(h)+g == Math.abs(h))
		  t=(a[ip][iq])/h;
		else {
		  theta=0.5*h/(a[ip][iq]);
		  t=1.0/(Math.abs(theta)+Math.sqrt(1.0+theta*theta));
		  if (theta < 0.0) t = -t;
		}
		c=1.0/Math.sqrt(1+t*t);
		s=t*c;
		tau=s/(1.0+c);
		h=t*a[ip][iq];
		z[ip] -= h;
		z[iq] += h;
		d[ip] -= h;
		d[iq] += h;
		a[ip][iq]=0.0;
		for (j=1;j<=ip-1;j++) {
		  ROTATE(a,j,ip,j,iq,tau,s);
		  }
		for (j=ip+1;j<=iq-1;j++) {
		  ROTATE(a,ip,j,j,iq,tau,s);
		  }
		for (j=iq+1;j<=n;j++) {
		  ROTATE(a,ip,j,iq,j,tau,s);
		  }
		for (j=1;j<=n;j++) {
		  ROTATE(v,j,ip,j,iq,tau,s);
		  }
		++nrot;
	      }
	    }
	  }
	  for (ip=1;ip<=n;ip++) {
	    b[ip] += z[ip];
	    d[ip]=b[ip];
	    z[ip]=0.0;
	  }
	}
	//printf("Too many iterations in routine JACOBI");
      }
    

    //  Perform the Cholesky decomposition    
    // Return the lower triangular L  such that L*L'=A  
    private void choldc(double a[][], int n, double l[][])
      {
	int i,j,k;
	double sum;
	double p[] = new double[n+1];
	
	for (i=1; i<=n; i++)  {
	  for (j=i; j<=n; j++)  {
	    for (sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];
	    if (i == j) {
	      if (sum<=0.0)  
		// printf("\nA is not poitive definite!");
		{}
	      else 
		p[i]=Math.sqrt(sum); }
	    else 
	      {
		a[j][i]=sum/p[i];
	      }
	  }
	}       
	for (i=1; i<=n; i++)  
	  for (j=i; j<=n; j++)  
	    if (i==j)
	      l[i][i] = p[i];
	    else
	      {
		l[j][i]=a[j][i];  
		l[i][j]=0.0;
	      }    
      }


    /********************************************************************/
    /**    Calcola la inversa della matrice  B mettendo il risultato   **/
    /**    in InvB . Il metodo usato per l'inversione e' quello di     **/
    /**    Gauss-Jordan.   N e' l'ordine della matrice .               **/
    /**    ritorna 0 se l'inversione  corretta altrimenti ritorna     **/
    /**    SINGULAR .                                                  **/
    /********************************************************************/
    int inverse(double TB[][], double InvB[][], int N) {  
      int k,i,j,p,q;
      double mult;
      double D,temp;
      double maxpivot;
      int npivot;
      double B[][] = new double [N+1][N+2];
      double A[][] = new double [N+1][2*N+2];
      double C[][] = new double [N+1][N+1];
      double eps = 10e-20;
      
      
      for(k=1;k<=N;k++)
	for(j=1;j<=N;j++)
	  B[k][j]=TB[k][j];
      
      for (k=1;k<=N;k++)
	{
	  for (j=1;j<=N+1;j++)
	    A[k][j]=B[k][j];
	  for (j=N+2;j<=2*N+1;j++)
	    A[k][j]=(float)0;
	  A[k][k-1+N+2]=(float)1;
	}
      for (k=1;k<=N;k++)
	{
	  maxpivot=Math.abs((double)A[k][k]);
	  npivot=k;
	  for (i=k;i<=N;i++)
	    if (maxpivot<Math.abs((double)A[i][k]))
	      {
		maxpivot=Math.abs((double)A[i][k]);
		npivot=i;
	      }
	  if (maxpivot>=eps)
	    {      if (npivot!=k)
		     for (j=k;j<=2*N+1;j++)
		       {
			 temp=A[npivot][j];
			 A[npivot][j]=A[k][j];
			 A[k][j]=temp;
		       } ;
		   D=A[k][k];
		   for (j=2*N+1;j>=k;j--)
		     A[k][j]=A[k][j]/D;
		   for (i=1;i<=N;i++)
		     {
		       if (i!=k)
			 {