ellifit.java

It help you to find the center of an ellipse

			   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;  

      // Generate normals linspace
      for (i=1, theta=0.0; i<=npts; i++, theta+=(pi/npts)) {
	u[1][i] = Math.cos(theta);
        u[2][i] = Math.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] = Math.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] /=  Math.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]<minev && Math.abs(d[i])>zero)	
		   solind = i;
	      break;
           case(FPF):
	      for (i=1; i<=6; i++)
		 if (d[i]<0 && Math.abs(d[i])>zero)	
		   solind = i;
          }
	// Now fetch the right solution
	for (int j=1;j<=6;j++)
	  pvec[j] = sol[j][solind];
	//pv(pvec,"the solution");

	// ...and plot it	
	draw_conic(pvec,npts,XY);	
	for (int i=1; i<npts; i++) {
	  if (XY[1][i]==-1 || XY[1][i+1]==-1  )
	    continue;
	  else
	    if (i<npts)
	      g.drawLine((int)XY[1][i],(int)XY[2][i],(int)XY[1][i+1],(int)XY[2][i+1] );
	    else
	      g.drawLine((int)XY[1][i],(int)XY[2][i],(int)XY[1][1],(int)XY[2][1] );
	}	
   }
}

class CurveControls extends Panel {
    CurvePanel target;
    Checkbox cb_add;
    Checkbox cb_move;
    Checkbox cb_delete;
    String st_add_label = "Add Points";
    String st_move_label = "Move Points";
    String st_delete_label = "Delete Points";

    public CurveControls(CurvePanel target) {
	this.target = target;
	setLayout(new FlowLayout(FlowLayout.CENTER));
	setBackground(Color.lightGray);
	Button clear = new Button("Clear");
	add("West",clear);

	CheckboxGroup action_group = new CheckboxGroup();
	add(cb_add = new Checkbox(st_add_label, action_group, true));
	add(cb_move = new Checkbox(st_move_label, action_group, false));
	add(cb_delete = new Checkbox(st_delete_label, action_group, false));
    }

    public void paint(Graphics g) {
	Rectangle r = bounds();

	g.setColor(Color.lightGray);
	g.draw3DRect(0, 0, r.width, r.height, false);
    }

    public boolean action(Event e, Object arg) {
	if (e.target instanceof Checkbox) {
	    String cbox = ((Checkbox)(e.target)).getLabel();
	    if (cbox.equals(st_add_label)) {
	    	target.setAction(CurvePanel.ADD);
	    } else if (cbox.equals(st_move_label)) {
	    	target.setAction(CurvePanel.MOVE);
	    } else if (cbox.equals(st_delete_label)) {
	    	target.setAction(CurvePanel.DELETE);
	    }
	} else if (e.target instanceof Button) {
	    String button = ((Button)(e.target)).getLabel();
	    if (button.equals("Clear")) {
	    	target.clearPoints();

		// After clearing the control points, put the user back into
		// ADD mode, since none of the other modes make any sense.
		cb_add.setState(true);
		cb_delete.setState(false);
		cb_move.setState(false);
	    	target.setAction(CurvePanel.ADD);
		target.repaint();
	    }
	}
	return true;
    }
}

class CurveControls2 extends Panel {
    CurvePanel target;
    Checkbox cb_hermite;
    Checkbox cb_bezier;
    Checkbox cb_bspline;
    String st_book_label = "Bookstein";
    String st_taubin_label = "Taubin";
    String st_fpf_label = "Ellipse-specific";
    Scrollbar scroll;
    Label precision_display;

    public CurveControls2(CurvePanel target) {
	this.target = target;
	setLayout(new FlowLayout(1));

	CheckboxGroup type_group = new CheckboxGroup();
	add(cb_hermite = new Checkbox(st_book_label, type_group, true));
	add(cb_bezier = new Checkbox(st_taubin_label, type_group, false));
	//	add(cb_bezier = new Checkbox("foo", type_group, false));
	add(cb_bspline = new Checkbox(st_fpf_label, type_group, false));

	// Set up the scrollbar: Value 15, Page 5, Min 2, Max 40
	//add(scroll = new Scrollbar(Scrollbar.HORIZONTAL, 15, 5, 2, 40));
	//target.setPrecision(scroll.getValue());
	//add(precision_display = new Label("Precision: "+
	//	String.valueOf(scroll.getValue()),Label.LEFT));
    }

    public void paint(Graphics g) {
	Rectangle r = bounds();

	g.setColor(Color.lightGray);
	g.draw3DRect(0, 0, r.width, r.height, false);
    }

    public boolean handleEvent(Event e) {
	switch (e.id) {
	  case Event.SCROLL_LINE_DOWN:
	  case Event.SCROLL_LINE_UP:
	  case Event.SCROLL_PAGE_DOWN:
	  case Event.SCROLL_PAGE_UP:  
		target.repaint();	
		return true;
	  case Event.ACTION_EVENT:
	  // Handle other action events
		if (e.target instanceof Checkbox) {
		   String cbox = ((Checkbox)(e.target)).getLabel();
		   if (cbox.equals(st_book_label)) {
		      target.setCurveType(CurvePanel.BOOKSTEIN);
		   } else if (cbox.equals(st_taubin_label)) {
		     // system.out.println("------------------------------------------");
		     //          g.drawString("Taubin not implemented",e.x+5,e.y-5);
		      // System.out.println("Taubin algorithm not yet implemented!");
		      target.setCurveType(CurvePanel.TAUBIN);
		   } else if (cbox.equals(st_fpf_label)) {
		      target.setCurveType(CurvePanel.FPF);
		   }
		}
		target.repaint();
		return(true);
	  default:
		return(false);
	}
    }
}