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

Introduction

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

Prototype

public ArrayRealVector ebeMultiply(ArrayRealVector v) throws IllegalArgumentException 

Source Link

Document

Element-by-element multiplication.

Usage

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

/**
 * least-squares using algebraic cost function
 * <p>/*from  w  ww .j  a v a  2s.c  o m*/
 * 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;
}

---