Example usage for org.apache.commons.math.linear ArrayRealVector mapMultiplyToSelf

Introduction

In this page you can find the example usage for org.apache.commons.math.linear ArrayRealVector mapMultiplyToSelf.

Prototype

@Override
public RealVector mapMultiplyToSelf(double d) 

Source Link

Usage

From source file:uk.ac.diamond.scisoft.analysis.fitting.AngleDerivativeFunction.java

/**
 * least-squares using algebraic cost function
 * <p>/*www  .  j  a v  a 2 s. com*/
 * This uses the method in "Numerically stable direct least squares fitting of ellipses"
 * by R. Halir and J. Flusser, Proceedings of the 6th International Conference in Central Europe
 * on Computer Graphics and Visualization. WSCG '98. CZ, Pilsen 1998, pp 125-132. 
 * <p>
 * @param x
 * @param y
 * @return geometric parameters
 */
private static double[] quickfit(AbstractDataset x, AbstractDataset y) {
    final AbstractDataset xx = Maths.square(x);
    final AbstractDataset yy = Maths.square(y);
    final AbstractDataset xxx = Maths.multiply(xx, x);
    final AbstractDataset yyy = Maths.multiply(yy, y);
    final AbstractDataset xy = Maths.multiply(x, y);

    Matrix S1 = new Matrix(3, 3);
    S1.set(0, 0, LinearAlgebra.dotProduct(xx, xx).getDouble());
    S1.set(0, 1, LinearAlgebra.dotProduct(xxx, y).getDouble());
    S1.set(0, 2, LinearAlgebra.dotProduct(xx, yy).getDouble());

    S1.set(1, 0, S1.get(0, 1));
    S1.set(1, 1, S1.get(0, 2));
    S1.set(1, 2, LinearAlgebra.dotProduct(x, yyy).getDouble());

    S1.set(2, 0, S1.get(0, 2));
    S1.set(2, 1, S1.get(1, 2));
    S1.set(2, 2, LinearAlgebra.dotProduct(yy, yy).getDouble());

    Matrix S2 = new Matrix(3, 3);
    S2.set(0, 0, ((Number) xxx.sum()).doubleValue());
    S2.set(0, 1, LinearAlgebra.dotProduct(xx, y).getDouble());
    S2.set(0, 2, ((Number) xx.sum()).doubleValue());

    S2.set(1, 0, S2.get(0, 1));
    S2.set(1, 1, LinearAlgebra.dotProduct(x, yy).getDouble());
    S2.set(1, 2, ((Number) xy.sum()).doubleValue());

    S2.set(2, 0, S2.get(1, 1));
    S2.set(2, 1, ((Number) yyy.sum()).doubleValue());
    S2.set(2, 2, ((Number) yy.sum()).doubleValue());

    Matrix S3 = new Matrix(3, 3);
    S3.set(0, 0, S2.get(0, 2));
    S3.set(0, 1, S2.get(1, 2));
    S3.set(0, 2, ((Number) x.sum()).doubleValue());

    S3.set(1, 0, S3.get(0, 1));
    S3.set(1, 1, S2.get(2, 2));
    S3.set(1, 2, ((Number) y.sum()).doubleValue());

    S3.set(2, 0, S3.get(0, 2));
    S3.set(2, 1, S3.get(1, 2));
    S3.set(2, 2, x.getSize());

    Matrix T = S3.solve(S2.transpose()).uminus();

    Matrix M = S1.plus(S2.times(T));

    Matrix Cinv = new Matrix(new double[] { 0, 0, 0.5, 0, -1.0, 0, 0.5, 0, 0 }, 3);
    Matrix Mp = Cinv.times(M);

    //      System.err.println("M " + Arrays.toString(Mp.getRowPackedCopy()));
    Matrix V = Mp.eig().getV();
    //      System.err.println("V " + Arrays.toString(V.getRowPackedCopy()));
    double[][] mv = V.getArray();
    ArrayRealVector v1 = new ArrayRealVector(mv[0]);
    ArrayRealVector v2 = new ArrayRealVector(mv[1]);
    ArrayRealVector v3 = new ArrayRealVector(mv[2]);

    v1.mapMultiplyToSelf(4);

    ArrayRealVector v = v1.ebeMultiply(v3).subtract(v2.ebeMultiply(v2));
    double[] varray = v.getData();
    int i = 0;
    for (; i < 3; i++) {
        if (varray[i] > 0)
            break;
    }
    if (i == 3) {
        throw new IllegalArgumentException("Could not find solution that satifies constraint");
    }

    v = new ArrayRealVector(new double[] { mv[0][i], mv[1][i], mv[2][i] });
    varray = v.getDataRef();
    final double ca = varray[0];
    final double cb = varray[1];
    final double cc = varray[2];
    Array2DRowRealMatrix mt = new Array2DRowRealMatrix(T.getArray(), false);
    varray = mt.operate(varray);
    final double cd = varray[0];
    final double ce = varray[1];
    final double cf = varray[2];

    //      System.err.println(String.format("Algebraic: %g, %g, %g, %g, %g, %g", ca, cb, cc, cd, ce, cf));
    final double disc = cb * cb - 4. * ca * cc;
    if (disc >= 0) {
        throw new IllegalArgumentException("Solution is not an ellipse");
    }

    if (cb == 0) {
        throw new IllegalArgumentException("Solution is a circle");
    }

    double[] qparameters = new double[PARAMETERS];
    qparameters[3] = (2. * cc * cd - cb * ce) / disc;
    qparameters[4] = (2. * ca * ce - cb * cd) / disc;

    final double sqrt = Math.sqrt((ca - cc) * (ca - cc) + cb * cb);
    qparameters[0] = -2. * (ca * ce * ce + cc * cd * cd + cf * cb * cb - cb * cd * ce - 4. * ca * cc * cf)
            / disc;
    qparameters[1] = qparameters[0] / (ca + cc + sqrt);
    qparameters[0] /= (ca + cc - sqrt);
    qparameters[0] = Math.sqrt(qparameters[0]);
    qparameters[1] = Math.sqrt(qparameters[1]);
    if (cb == 0) {
        qparameters[2] = 0.;
    } else {
        qparameters[2] = 0.5 * Math.atan2(cb, ca - cc);
    }
    if (qparameters[0] < qparameters[1]) {
        final double t = qparameters[0];
        qparameters[0] = qparameters[1];
        qparameters[1] = t;
    } else {
        qparameters[2] += Math.PI * 0.5;
    }

    return qparameters;
}

---