//**************************************************************************
// 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.
//**************************************************************************
package ellifit;
import java.awt.*;
import java.applet.*;
import java.lang.Double;
import java.util.Vector;
import javax.swing.JFrame;
import math.geom2d.Point2D;
import math.geom2d.conic.Ellipse2D;
import org.apache.commons.math3.util.FastMath;
public class ElliFit extends Applet {
/**
*
*/
private static final long serialVersionUID = 1L;
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[]) {
JFrame f = new JFrame("ElliFIt");
ElliFit curve = new ElliFit();
curve.init();
curve.start();
f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
f.setSize(400,400);
f.add("Center", curve);
f.setVisible(true);
}
}
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 = FPF; // 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.white);
// initialize
points.addElement(new ControlPoint(10,10));
points.addElement(new ControlPoint(10,15));
points.addElement(new ControlPoint(15,15));
points.addElement(new ControlPoint(20,55));
points.addElement(new ControlPoint(45,45));
points.addElement(new ControlPoint(40,45));
}
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 += FastMath.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*FastMath.abs(a[ip][iq]);
if (i > 4 && FastMath.abs(d[ip])+g == FastMath.abs(d[ip])
&& FastMath.abs(d[iq])+g == FastMath.abs(d[iq]))
a[ip][iq]=0.0;
else if (FastMath.abs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if (FastMath.abs(h)+g == FastMath.abs(h))
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
t=1.0/(FastMath.abs(theta)+FastMath.sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/FastMath.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]=FastMath.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=FastMath.abs((double)A[k][k]);
npivot=k;
for (i=k;i<=N;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)
{
mult=A[i][k];
for (j=2*N+1;j>=k;j--)
A[i][j]=A[i][j]-mult*A[k][j] ;
}
}
}
else
{ // printf("\n The matrix may be singular !!") ;
return(-1);
};
}
/** Copia il risultato nella matrice InvB ***/
for (k=1,p=1;k<=N;k++,p++)
for (j=N+2,q=1;j<=2*N+1;j++,q++)
InvB[p][q]=A[k][j];
return(0);
} /* End of INVERSE */
private void AperB(double _A[][], double _B[][], double _res[][],
int _righA, int _colA, int _righB, int _colB) {
int p,q,l;
for (p=1;p<=_righA;p++)
for (q=1;q<=_colB;q++)
{ _res[p][q]=0.0;
for (l=1;l<=_colA;l++)
_res[p][q]=_res[p][q]+_A[p][l]*_B[l][q];
}
}
private void A_TperB(double _A[][], double _B[][], double _res[][],
int _righA, int _colA, int _righB, int _colB) {
int p,q,l;
for (p=1;p<=_colA;p++)
for (q=1;q<=_colB;q++)
{ _res[p][q]=0.0;
for (l=1;l<=_righA;l++)
_res[p][q]=_res[p][q]+_A[l][p]*_B[l][q];
}
}
private void AperB_T(double _A[][], double _B[][], double _res[][],
int _righA, int _colA, int _righB, int _colB) {
int p,q,l;
for (p=1;p<=_colA;p++)
for (q=1;q<=_colB;q++)
{ _res[p][q]=0.0;
for (l=1;l<=_righA;l++)
_res[p][q]=_res[p][q]+_A[p][l]*_B[q][l];
}
}
public void pv(double v[], java.lang.String str) {
System.out.println("------------" + str + "--------------");
System.out.println(" " + v[1] + " " + v[2] + " " + v[3]
+ " " + v[4] + " " + v[5] + " " + v[6] );
System.out.println("------------------------------------------");
}
public void pm(double S[][], java.lang.String str) {
System.out.println("------------" + str + "--------------");
System.out.println(" " + S[1][1] + " " + S[1][2] + " " + S[1][3]
+ " " + S[1][4] + " " + S[1][5] + " " + S[1][6] );
System.out.println(" " + S[2][1] + " " + S[2][2] + " " + S[2][3]
+ " " + S[2][4] + " " + S[2][5] + " " + S[2][6] );
System.out.println(" " + S[3][1] + " " + S[3][2] + " " + S[3][3]
+ " " + S[3][4] + " " + S[3][5] + " " + S[3][6] );
System.out.println(" " + S[4][1] + " " + S[4][2] + " " + S[4][3]
+ " " + S[4][4] + " " + S[4][5] + " " + S[4][6] );
System.out.println(" " + S[5][1] + " " + S[5][2] + " " + S[5][3]
+ " " + S[5][4] + " " + S[5][5] + " " + S[5][6] );
System.out.println(" " + S[6][1] + " " + S[6][2] + " " + S[6][3]
+ " " + S[6][4] + " " + S[6][5] + " " + S[6][6] );
System.out.println("------------------------------------------");
}
public void draw_conic(double pvec[], int nptsk, double points[][]) {
int npts=nptsk/2;
double u[][] = new double[3][npts+1];
double Aiu[][] = new double[3][npts+1];
double L[][] = new double[3][npts+1];
double B[][] = new double[3][npts+1];
double Xpos[][] = new double[3][npts+1];
double Xneg[][] = new double[3][npts+1];
double ss1[][] = new double[3][npts+1];
double ss2[][] = new double[3][npts+1];
double lambda[] = new double[npts+1];
double uAiu[][] = new double[3][npts+1];
double A[][] = new double[3][3];
double Ai[][] = new double[3][3];
double Aib[][] = new double[3][2];
double b[][] = new double[3][2];
double r1[][] = new double[2][2];
double Ao, Ax, Ay, Axx, Ayy, Axy;
double pi = 3.14781;
double theta;
int i;
int j;
double kk;
Ao = pvec[6];
Ax = pvec[4];
Ay = pvec[5];
Axx = pvec[1];
Ayy = pvec[3];
Axy = pvec[2];
A[1][1] = Axx; A[1][2] = Axy/2;
A[2][1] = Axy/2; A[2][2] = Ayy;
b[1][1] = Ax; b[2][1] = Ay;
/*
System.out.println("\nDebug=");
System.out.println(Axx);
System.out.println(Ayy);
System.out.println(Ax);
System.out.println(Ay);
*/
// Generate normals linspace
for (i=1, theta=0.0; i<=npts; i++, theta+=(pi/npts)) {
u[1][i] = FastMath.cos(theta);
u[2][i] = FastMath.sin(theta); }
inverse(A,Ai,2);
AperB(Ai,b,Aib,2,2,2,1);
A_TperB(b,Aib,r1,2,1,2,1);
r1[1][1] = r1[1][1] - 4*Ao;
AperB(Ai,u,Aiu,2,2,2,npts);
for (i=1; i<=2; i++)
for (j=1; j<=npts; j++)
uAiu[i][j] = u[i][j] * Aiu[i][j];
for (j=1; j<=npts; j++) {
if ( (kk=(r1[1][1] / (uAiu[1][j]+uAiu[2][j]))) >= 0.0)
lambda[j] = FastMath.sqrt(kk);
else
lambda[j] = -1.0; }
// Builds up B and L
for (j=1; j<=npts; j++)
L[1][j] = L[2][j] = lambda[j];
for (j=1; j<=npts; j++) {
B[1][j] = b[1][1];
B[2][j] = b[2][1]; }
for (j=1; j<=npts; j++) {
ss1[1][j] = 0.5 * ( L[1][j] * u[1][j] - B[1][j]);
ss1[2][j] = 0.5 * ( L[2][j] * u[2][j] - B[2][j]);
ss2[1][j] = 0.5 * ( -L[1][j] * u[1][j] - B[1][j]);
ss2[2][j] = 0.5 * ( -L[2][j] * u[2][j] - B[2][j]); }
AperB(Ai,ss1,Xpos,2,2,2,npts);
AperB(Ai,ss2,Xneg,2,2,2,npts);
for (j=1; j<=npts; j++) {
if (lambda[j]==-1.0) {
points[1][j] = -1.0;
points[2][j] = -1.0;
points[1][j+npts] = -1.0;
points[2][j+npts] = -1.0;
}
else {
points[1][j] = Xpos[1][j];
points[2][j] = Xpos[2][j];
points[1][j+npts] = Xneg[1][j];
points[2][j+npts] = Xneg[2][j];
}
}
}
public void paint(Graphics g) {
int np = points.size(); // number of points
double D[][] = new double[np+1][7];
double S[][] = new double[7][7];
double Const[][] = new double[7][7];
double temp[][] = new double [7][7];
double L[][] = new double [7][7];
double C[][] = new double [7][7];
double invL[][] = new double [7][7];
double d[] = new double [7];
double V[][] = new double [7][7];
double sol[][] = new double [7][7];
double tx,ty;
int nrot=0;
int npts=50;
double XY[][] = new double[3][npts+1];
double pvec[] = new double[7];
g.setColor(getForeground());
g.setPaintMode();
// draw a border around the canvas
g.drawRect(0,0, size().width-1, size().height-1);
// draw the control points
for (int i=0; i < np; i++) {
ControlPoint p = (ControlPoint)points.elementAt(i);
g.drawOval(p.x-p.PT_SIZE, p.y-p.PT_SIZE, p.PT_SIZE*2, p.PT_SIZE*2);
//g.drawString(String.valueOf(i),p.x+p.PT_SIZE,p.y-p.PT_SIZE);
}
switch (mode) {
case(FPF):
//System.out.println("FPF mode");
Const[1][3]=-2;
Const[2][2]=1;
Const[3][1]=-2;
break;
case(TAUBIN):
g.drawString(warn_taub_str,size().width/18 , size().height/18 );
break;
case(BOOKSTEIN):
//System.out.println("BOOK mode");
Const[1][1]=2;
Const[2][2]=1;
Const[3][3]=2;
}
if (np<6)
return;
// Now first fill design matrix
for (int i=1; i <= np; i++)
{
tx = ((ControlPoint)points.elementAt(i-1)).x;
ty = ((ControlPoint)points.elementAt(i-1)).y;
D[i][1] = tx*tx;
D[i][2] = tx*ty;
D[i][3] = ty*ty;
D[i][4] = tx;
D[i][5] = ty;
D[i][6] = 1.0;
}
//pm(Const,"Constraint");
// Now compute scatter matrix S
A_TperB(D,D,S,np,6,np,6);
//pm(S,"Scatter");
choldc(S,6,L);
//pm(L,"Cholesky");
inverse(L,invL,6);
//pm(invL,"inverse");
AperB_T(Const,invL,temp,6,6,6,6);
AperB(invL,temp,C,6,6,6,6);
//pm(C,"The C matrix");
jacobi(C,6,d,V,nrot);
//pm(V,"The Eigenvectors"); /* OK */
//pv(d,"The eigevalues");
A_TperB(invL,V,sol,6,6,6,6);
//pm(sol,"The GEV solution unnormalized"); /* SOl */
// Now normalize them
for (int j=1;j<=6;j++) /* Scan columns */
{
double mod = 0.0;
for (int i=1;i<=6;i++)
mod += sol[i][j]*sol[i][j];
for (int i=1;i<=6;i++)
sol[i][j] /= FastMath.sqrt(mod);
}
//pm(sol,"The GEV solution"); /* SOl */
double zero=10e-20;
double minev=10e+20;
int solind=0;
switch (mode) {
case(BOOKSTEIN): // smallest eigenvalue
for (i=1; i<=6; i++)
if (d[i]zero)
solind = i;
break;
case(FPF):
for (i=1; i<=6; i++)
if (d[i]<0 && FastMath.abs(d[i])>zero)
solind = i;
}
// Now fetch the right solution
for (int j=1;j<=6;j++)
pvec[j] = sol[j][solind];
// he ellipse is defined through the 6 coefficients of its algebraic equation:
// Axx*x**2 + Axy*x*y + Ayy*y**2 + Ax*x + Ay*y + A0 = 0
pv(pvec,"the solution.");
double A0 = pvec[6];
double Ax = pvec[4];
double Ay = pvec[5];
double Axx = pvec[1];
double Ayy = pvec[3];
double Axy = pvec[2];
double aC= Axx; double bC = Axy/2;
double cC= Ayy;
double dC = Ax/2; double fC = Ay/2;
double gC=A0;
//ax^2+2bxy+cy^2+2dx+2fy+g=0
//b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
double num = bC*bC-aC*cC;
double x0=(cC*dC-bC*fC)/num;
double y0=(aC*fC-bC*dC)/num;
System.out.println("Center x="+Double.toString(x0)+" y="+Double.toString(y0));
for (int i=0; i < np; i++) {
ControlPoint p = (ControlPoint)points.elementAt(i);
// p.x, p.y - points X, Y
//wOval(p.x-p.PT_SIZE, p.y-p.PT_SIZE, p.PT_SIZE*2, p.PT_SIZE*2);
//g.drawString(String.valueOf(i),p.x+p.PT_SIZE,p.y-p.PT_SIZE);
}
//ellipse_angle_of_rotation
//b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
//return 0.5*np.arctan(2*b/(a-c))
double angle=0.5*FastMath.atan( (2*bC/(aC-cC)) );
System.out.println("Angle of rotation="+Double.toString(angle));
double up = 2*(aC*fC*fC+cC*dC*dC+gC*bC*bC-2*bC*dC*fC-aC*cC*gC);
double xx=(bC*bC-aC*cC);
double down1=xx*( (cC-aC)*FastMath.sqrt(1+4*bC*bC/((aC-cC)*(aC-cC)))-(cC+aC));
double down2=xx*( (aC-cC)*FastMath.sqrt(1+4*bC*bC/((aC-cC)*(aC-cC)))-(cC+aC));
double res1=FastMath.sqrt(up/down1);
double res2=FastMath.sqrt(up/down2);
double major=res1;
double minor=res2;
if (major
Ads help maintain this website.