geogebra.common.kernel.algos.AlgoConicFivePoints.java Source code

Java tutorial

Introduction

Here is the source code for geogebra.common.kernel.algos.AlgoConicFivePoints.java

Source

/* 
GeoGebra - Dynamic Mathematics for Everyone
http://www.geogebra.org
    
This file is part of GeoGebra.
    
This program is free software; you can redistribute it and/or modify it 
under the terms of the GNU General Public License as published by 
the Free Software Foundation.
    
 */

/*
 * AlgoConicFivePoints.java
 *
 * Created on 15. November 2001, 21:37
 */

package geogebra.common.kernel.algos;

import geogebra.common.euclidian.EuclidianConstants;
import geogebra.common.kernel.Construction;
import geogebra.common.kernel.Kernel;
import geogebra.common.kernel.LocusEquation;
import geogebra.common.kernel.StringTemplate;
import geogebra.common.kernel.commands.Commands;
import geogebra.common.kernel.geos.GeoConic;
import geogebra.common.kernel.geos.GeoElement;
import geogebra.common.kernel.geos.GeoLine;
import geogebra.common.kernel.geos.GeoPoint;
import geogebra.common.kernel.geos.GeoVec3D;
import geogebra.common.kernel.kernelND.GeoConicND;
import geogebra.common.kernel.kernelND.GeoPointND;

import java.util.ArrayList;

/**
 *
 * @author Markus
 * @version
 */
public class AlgoConicFivePoints extends AlgoElement {

    // #4156 these are rather arbitrary tradeoffs between compatibility and
    // numeric stability
    private static final double MULTIPLIER_MIN = 0.001;
    private static final double MULTIPLIER_MAX = 1000;
    protected GeoPoint[] P; // input five points
    protected GeoConicND conic; // output
    private boolean criticalCase; // true when 5 points is on a parabola

    private double[][] A, B, C;
    private double e1, e2;
    private GeoVec3D[] line;
    private int i, j;

    public AlgoConicFivePoints(Construction cons, String label, GeoPointND[] inputP) {
        this(cons, inputP);
        conic.setLabel(label);
    }

    protected void setInputPoints() {
        input = P;
    }

    protected GeoPoint[] createPoints2D(GeoPointND[] inputP) {
        return (GeoPoint[]) inputP;
    }

    public AlgoConicFivePoints(Construction cons, GeoPointND[] inputP) {
        super(cons);
        this.P = createPoints2D(inputP);
        conic = newGeoConic(cons);

        setInputOutput(); // for AlgoElement

        line = new GeoVec3D[4];
        for (i = 0; i < 4; i++) {
            line[i] = new GeoLine(cons);
        }
        A = new double[3][3];
        B = new double[3][3];
        C = new double[3][3];
        checkCriticalCase();
        initCoords();
        compute();
        addIncidence();

        /*
         * moved into addIncidence() for (int i=0; i < P.length; i++) {
         * conic.addPointOnConic(P[i]); }
         */

    }

    /**
     * init Coords values
     */
    protected void initCoords() {
        // none here
    }

    /**
     * 
     * @param cons
     *            construction
     * @return output conic
     */
    protected GeoConicND newGeoConic(Construction cons) {
        return new GeoConic(cons);
    }

    private void checkCriticalCase() {
        criticalCase = false;

        for (int i = 0; i < P.length; i++) {
            if (P[i].getIncidenceList() == null)
                return;
        }

        ArrayList<GeoElement> firstList = P[0].getIncidenceList();

        for (int j = 0; j < firstList.size(); j++) {
            if (firstList.get(j).isGeoConic()) {
                GeoConic p = (GeoConic) firstList.get(j);
                if (p.getType() == GeoConic.CONIC_PARABOLA) {
                    criticalCase = true;
                    for (int i = 1; i < 5; i++) {
                        if (!P[i].getIncidenceList().contains(p)) {
                            criticalCase = false;
                            break;
                        }
                    }
                }
            }

            if (criticalCase) {
                break;
            }
        }

    }

    /**
     * @author Tam
     * 
     *         for special cases of e.g. AlgoIntersectLineConic
     */
    private void addIncidence() {
        for (int i = 0; i < P.length; ++i) {
            P[i].addIncidence(conic, false);
        }
    }

    @Override
    public Commands getClassName() {
        return Commands.Conic;
    }

    @Override
    public int getRelatedModeID() {
        return EuclidianConstants.MODE_CONIC_FIVE_POINTS;
    }

    // for AlgoElement
    @Override
    protected void setInputOutput() {

        setInputPoints();

        super.setOutputLength(1);
        super.setOutput(0, conic);
        setDependencies(); // done by AlgoElement
    }

    public GeoConicND getConic() {
        return conic;
    }

    GeoPoint[] getPoints() {
        return P;
    }

    /**
     * Method created for LocusEqu project.
     * 
     * @return a copy of inner array so it cannot be manipulated from outside.
     */
    public GeoPoint[] getAllPoints() {
        GeoPoint[] copy = new GeoPoint[this.getPoints().length];
        System.arraycopy(this.getPoints(), 0, copy, 0, copy.length);
        return copy;
    }

    // compute conic through five points P[0] ... P[4]
    // with Plcker  method
    @Override
    public void compute() {
        // compute lines P0 P1, P2 P3,
        // P0 P2, P1 P3
        GeoVec3D.cross(P[0], P[1], line[0]);
        GeoVec3D.cross(P[2], P[3], line[1]);
        GeoVec3D.cross(P[0], P[2], line[2]);
        GeoVec3D.cross(P[1], P[3], line[3]);

        // compute degenerate conics A = line[0] u line[1],
        // B = line[2] u line[3]
        degCone(line[0], line[1], A);
        degCone(line[2], line[3], B);
        e1 = evalMatrix(B, P[4]);
        e2 = -evalMatrix(A, P[4]);

        // try to avoid tiny/huge value for matrix
        if (shouldInvert(e1) && shouldInvert(e2)) {

            double tmp = e1;

            e1 = 1 / e2;
            e2 = 1 / tmp;
        }

        linComb(A, B, e1, e2, C);

        /***
         * Perturbation method to estimate the error of "detS" of the conic
         * 
         * The following is a random perturbation method to estimate detS it is
         * commented out because (1) it is not deterministic (2) it still can't
         * solve the problem of #1294
         * 
         * Tam 6/9/2012
         */

        /*
         * // compute a perturbed Cpert kernel.setSilentMode(true); Ppert = new
         * GeoPoint2[5]; for (int i=0; i<5; i++) { Ppert[i] = new
         * GeoPoint2(P[i]);
         * 
         * }
         */

        /*
         * double maxDistSqr = 0; double temp; int maxDistP1 = 0; int maxDistP2
         * = 0;
         * 
         * 
         * for (int i=0; i<5; i++) { temp = P[0].distanceSqr(P[i]); if ( temp >
         * maxDistSqr) { maxDistP2 = i; maxDistSqr = temp; } } for (int i=0;
         * i<5; i++) { temp = P[i].distanceSqr(P[maxDistP2]); if (temp >
         * maxDistSqr) { maxDistP1 = i; maxDistSqr = temp; } }
         * 
         * if (maxDistSqr!=0) { double t = Kernel.EPSILON/Math.sqrt(maxDistSqr);
         * Ppert[maxDistP1].setCoords( P[maxDistP1].inhomX * (1+t) -
         * P[maxDistP2].inhomX * t, P[maxDistP1].inhomY * (1+t) -
         * P[maxDistP2].inhomY * t, 1 ); Ppert[maxDistP2].setCoords(
         * P[maxDistP2].inhomX * (1+t) - P[maxDistP1].inhomX * t,
         * P[maxDistP2].inhomY * (1+t) - P[maxDistP1].inhomY * t, 1 );
         */

        /*
         * perturbation for finding out deltaS
         */
        /*
         * delta = Kernel.MIN_PRECISION; int repetition = 5; for (int m=0; m<3;
         * m++) for (int n=0; n<3; n++) Cmin[m][n]=Double.POSITIVE_INFINITY;
         * 
         * for (int k=0; k<repetition; k++) { for (int i=0; i<5; i++) {
         * Ppert[i].randomizeForErrorEstimation(); }
         * 
         * GeoVec3D.cross(Ppert[0], Ppert[1], line[0]); GeoVec3D.cross(Ppert[2],
         * Ppert[3], line[1]); GeoVec3D.cross(Ppert[0], Ppert[2], line[2]);
         * GeoVec3D.cross(Ppert[1], Ppert[3], line[3]); degCone(line[0],
         * line[1], A); degCone(line[2], line[3], B); l = evalMatrix(B,
         * Ppert[4]); m = -evalMatrix(A, Ppert[4]); linComb(A, B, l, m, Cpert);
         * 
         * //calculate the estimation of error of detS delta = Math.min(delta,
         * Math.abs( Cpert[0][0]*Cpert[1][1] -
         * (Cpert[0][1]+Cpert[1][0])*(Cpert[0][1]+Cpert[1][0])/4)); for (int
         * m=0; m<3; m++) for (int n=m; n<3; n++)
         * Cmin[m][n]=Math.min(Cmin[m][n],
         * Math.abs((Cpert[m][n]+Cpert[n][m])/2));
         * 
         * }
         * 
         * conic.errDetS = delta; //TODO: this is not reasonable enough. if
         * (Math.abs(C[0][0])< Math.abs(Cmin[0][0])/10) C[0][0]=0; if
         * (Math.abs(C[1][1])< Math.abs(Cmin[1][1])/10) C[1][1]=0; if
         * (Math.abs(C[0][1] + C[1][0])< Math.abs( Cmin[0][1]) /5 ) { C[0][1]=0;
         * C[1][0]=0; }
         */
        /**
         * perturbation method ends
         */

        /***
         * Testing: use analytic method:
         * 
         * solving system of linear equations
         * 
         * uses:
         * 
         * import org.apache.commons.math.linear.AbstractFieldMatrix; import
         * org.apache.commons.math.linear.Array2DRowFieldMatrix; import
         * org.apache.commons.math.linear.Array2DRowRealMatrix; import
         * org.apache.commons.math.linear.FieldMatrix; import
         * org.apache.commons.math.linear.DecompositionSolver; import
         * org.apache.commons.math.linear.QRDecompositionImpl; import
         * org.apache.commons.math.linear.RealMatrix; import
         * org.apache.commons.math.linear.SingularValueDecompositionImpl; import
         * org.apache.commons.math.util.BigReal;
         */
        /*
         * RealMatrix coeffM = new Array2DRowRealMatrix(5, 6);
         * 
         * for (int i = 0; i<5; i++) { coeffM.setRow(i, new double[]
         * {P[i].inhomX * P[i].inhomX,P[i].inhomX * P[i].inhomY,P[i].inhomY *
         * P[i].inhomY, P[i].inhomX , P[i].inhomY, 1.0}); }
         * 
         * SingularValueDecompositionImpl solver = new
         * SingularValueDecompositionImpl(coeffM);
         * 
         * solver.getSolver().solve(new double[] {0,0,0,0,0});
         * 
         * //RealMatrix test1 = solver.getV().multiply(solver.getVT());
         * //RealMatrix test2 = solver.getVT().multiply(solver.getV());
         * 
         * int key = -1; double keysum = 1; for (i=0; i<6; i++) { double sum =
         * 0; for (j=0; j<5; j++) { sum += solver.getV().getEntry(i,j) *
         * solver.getV().getEntry(i,j); } if (sum < keysum) { key = i; keysum =
         * sum; } }
         * 
         * double[] xx = new double[6]; double[] v6 = new double[6]; if
         * (!Kernel.isZero(1-keysum)) { xx[5] = 1/(1-keysum); for (int j=0; j<5;
         * j++) { xx[j] = -xx[5]*solver.getV().getEntry(key,j); } }
         * 
         * for (int i=0; i<5; i++) { v6[i] = 0; }
         * 
         * v6[key] = xx[5];
         * 
         * for (int i=0; i<6; i++) { for (int j=0; j<5; j++) { v6[i] += xx[j] *
         * solver.getV().getEntry(i,j); } }
         * 
         * RealMatrix checkSol = new Array2DRowRealMatrix(6,1);
         * checkSol.setColumn(0, v6); RealMatrix check =
         * coeffM.multiply(checkSol); //double[] solution =
         * solver.getg.getV().getColumn(5); for (int i=0;
         * i<check.getRowDimension(); i++) {
         * System.out.println(check.getRow(1).toString()); }
         * 
         * conic.setMatrix(new double[]
         * {v6[0],v6[2],v6[5],v6[1]/2,v6[3]/2,v6[4]/2});
         */
        /***
         * solving system of linear equations test ends
         ***/

        /***
         * critical case: five points lie on an unstable conic now only for
         * parabola. Need more tests for: one line; two lines; one point; two
         * points
         */

        if (criticalCase) {
            conic.errDetS = Double.POSITIVE_INFINITY;
        } else {
            conic.errDetS = Kernel.MIN_PRECISION;
        }

        conic.setMatrix(C);
        // System.out.println(conic.getTypeString());

    }

    private boolean shouldInvert(double d) {
        return (!Kernel.isZero(d) && Math.abs(d) < MULTIPLIER_MIN) || Math.abs(d) > MULTIPLIER_MAX;
    }

    // compute degenerate conic from lines a, b
    // the result is written into A as a NON-SYMMETRIC Matrix
    final private static void degCone(GeoVec3D a, GeoVec3D b, double[][] A) {
        // A = a . b^t
        A[0][0] = a.x * b.x;
        A[0][1] = a.x * b.y;
        A[0][2] = a.x * b.z;
        A[1][0] = a.y * b.x;
        A[1][1] = a.y * b.y;
        A[1][2] = a.y * b.z;
        A[2][0] = a.z * b.x;
        A[2][1] = a.z * b.y;
        A[2][2] = a.z * b.z;
    }

    // computes P.A.P, where A is a (possibly not symmetric) 3x3 matrix
    final private static double evalMatrix(double[][] A, GeoPoint P) {
        return A[0][0] * P.x * P.x + A[1][1] * P.y * P.y + A[2][2] * P.z * P.z + (A[0][1] + A[1][0]) * P.x * P.y
                + (A[0][2] + A[2][0]) * P.x * P.z + (A[1][2] + A[2][1]) * P.y * P.z;
    }

    // computes the linear combination C = l * A + m * B
    final private void linComb(double[][] A, double[][] B, double l, double m, double[][] C) {
        for (i = 0; i < 3; i++) {
            for (j = 0; j < 3; j++) {
                C[i][j] = l * A[i][j] + m * B[i][j];
            }
        }
    }

    @Override
    final public String toString(StringTemplate tpl) {
        // Michael Borcherds 2008-03-30
        // simplified to allow better Chinese translation
        return getLoc().getPlain("ConicThroughABCDE", P[0].getLabel(tpl), P[1].getLabel(tpl), P[2].getLabel(tpl),
                P[3].getLabel(tpl), P[4].getLabel(tpl));
    }

    @Override
    public boolean isLocusEquable() {
        return true;
    }

    @Override
    public EquationElementInterface buildEquationElementForGeo(GeoElement geo, EquationScopeInterface scope) {
        return LocusEquation.eqnConicFivePoints(geo, this, scope);
    }
}