//**************************************************************************
 // 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 yi.iu.processor.ellipse;

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

/**
 * Based on the code provided by Maurizlio Pilu
 *
 * http://www.hpl.hp.com/personal/Maurizio_Pilu/research/ellipse.htm
 */
public class ElliFit
    extends Applet {

  public static final int BOOKSTEIN = 0;
  public static final int TAUBIN = 1;
  public static final int FPF = 2;

  /**
   * Same as fit(int mode, Vector points).  Convert the java.awt.Point array
   * to ControlPoint first.
   * @param mode
   * @param points
   * @return
   */
  public static double[] fit(int mode, Point[] points) {
    Vector cp_vect = new Vector(points.length);
    for (int i = 0; i < points.length; i++) {
      cp_vect.add(new ControlPoint(points[i].x, points[i].y));
    }
    return (new CurvePanel()).fit(mode, cp_vect);
  }

  public static double[][] getEllipsePoints(double pvec[], int npts) {
    double points[][] = new double[3][npts + 1];
    (new CurvePanel()).draw_conic(pvec, npts, points);
    return points;
  }

  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 {
  private int mode = ElliFit.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 ElliFit.BOOKSTEIN:
      case ElliFit.TAUBIN:
      case ElliFit.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) {
            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];

    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 (ElliFit.TAUBIN):
        g.drawString(warn_taub_str, size().width / 18, size().height / 18);
        break;
    }

    double pvec[] = fit(mode, points);
    if (pvec == null)
      return;

    // ...and plot it
    draw_conic(pvec, npts, XY);
    for (int i = 1; i <= npts; i++) {
      if (XY[1][i] == -1 || (i < npts && 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]);
    }
  }

  /**
   * Fit ellipse to a set of point
   * @param mode
   * @param points
   * @return the parameter of ax^2+bxy+cy^2+dx+ey+f=0;
   */
  public double[] fit(int mode, Vector points) {
    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];

    switch (mode) {
      case (ElliFit.FPF):

        //System.out.println("FPF mode");
        Const[1][3] = -2;
        Const[2][2] = 1;
        Const[3][1] = -2;
        break;
      case (ElliFit.BOOKSTEIN):

        //System.out.println("BOOK mode");
        Const[1][1] = 2;
        Const[2][2] = 1;
        Const[3][3] = 2;

    }

    if (np < 6)
      return null;

    // 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 (ElliFit.BOOKSTEIN): // smallest eigenvalue
        for (i = 1; i <= 6; i++)
          if (d[i] < minev && Math.abs(d[i]) > zero)
            solind = i;
        break;
      case (ElliFit.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");

    return pvec;
  }
}

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(ElliFit.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(ElliFit.TAUBIN);
          }
          else if (cbox.equals(st_fpf_label)) {
            target.setCurveType(ElliFit.FPF);
          }
        }
        target.repaint();
        return (true);
      default:
        return (false);
    }
  }
}