Example usage for org.apache.commons.math3.linear RealMatrix getColumnVector

List of usage examples for org.apache.commons.math3.linear RealMatrix getColumnVector

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear RealMatrix getColumnVector.

Prototype

RealVector getColumnVector(int column) throws OutOfRangeException;

Source Link

Document

Get the entries at the given column index as a vector.

Usage

From source file:eagle.security.userprofile.model.kde.UserProfileKDEModeler.java

private void computeStats(RealMatrix m) {
    if (m.getColumnDimension() != this.cmdTypes.length) {
        LOG.error("Please fix the commands list in config file");
    }//  w  ww.ja v  a2 s  .  co  m

    statistics = new UserCommandStatistics[m.getColumnDimension()];

    for (int i = 0; i < m.getColumnDimension(); i++) {
        UserCommandStatistics stats = new UserCommandStatistics();
        stats.setCommandName(this.cmdTypes[i]);
        RealVector colData = m.getColumnVector(i);
        StandardDeviation deviation = new StandardDeviation();
        double stddev = deviation.evaluate(colData.toArray());

        if (LOG.isDebugEnabled())
            LOG.debug("Stddev is NAN ? " + (Double.isNaN(stddev) ? "yes" : "no"));
        if (stddev <= lowVarianceVal)
            stats.setLowVariant(true);
        else
            stats.setLowVariant(false);

        stats.setStddev(stddev);
        Mean mean = new Mean();
        double mu = mean.evaluate(colData.toArray());
        if (LOG.isDebugEnabled())
            LOG.debug("mu is NAN ? " + (Double.isNaN(mu) ? "yes" : "no"));

        stats.setMean(mu);
        statistics[i] = stats;
    }
}

From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java

@Override
public GaussianFactor division(GaussianFactor other) {
    if (!scope.contains(other.getVariables().getVariableIds())) {
        throw new FactorOperationException("Divisor scope " + other.getVariables() + " is not a subset of"
                + " this factor's scope " + scope);
    }//  w w  w . j a  va 2 s.  c o  m

    int[] otherMapping = scope.createContinuousVariableMapping(other.getVariables());

    RealMatrix newPrecisionMatrix = precisionMatrix.copy();

    RealMatrix otherPrecisionMatrix = other.getPrecisionMatrix();

    for (int i = 0; i < scope.size(); i++) {
        RealVector column = newPrecisionMatrix.getColumnVector(i);

        if (otherMapping[i] >= 0) {
            column = column.subtract(padVector(otherPrecisionMatrix.getColumnVector(otherMapping[i]),
                    scope.size(), otherMapping));
            newPrecisionMatrix.setColumnVector(i, column);
        }
    }

    RealVector newScaledMeanVector = scaledMeanVector.copy();
    RealVector otherScaledMeanVector = other.getScaledMeanVector();
    newScaledMeanVector = newScaledMeanVector
            .subtract(padVector(otherScaledMeanVector, scope.size(), otherMapping));

    double newNormalizationConstant = normalizationConstant - other.getNormalizationConstant();

    return new CanonicalGaussianFactor(scope, newPrecisionMatrix, newScaledMeanVector,
            newNormalizationConstant);
}

From source file:eagle.security.userprofile.model.eigen.UserProfileEigenModeler.java

private void computeStats(RealMatrix m) {

    if (m.getColumnDimension() != this.cmdTypes.length) {
        LOG.error("Please fix the commands list in config file");
        return;//from  w  ww .  j a  v  a  2s  . c  o m
    }
    statistics = new UserCommandStatistics[m.getColumnDimension()];
    for (int i = 0; i < m.getColumnDimension(); i++) {
        UserCommandStatistics stats = new UserCommandStatistics();
        stats.setCommandName(this.cmdTypes[i]);
        RealVector colData = m.getColumnVector(i);
        StandardDeviation deviation = new StandardDeviation();
        double stddev = deviation.evaluate(colData.toArray());
        //LOG.info("stddev is nan?" + (stddev == Double.NaN? "yes":"no"));
        if (stddev <= lowVarianceVal)
            stats.setLowVariant(true);
        else
            stats.setLowVariant(false);
        stats.setStddev(stddev);
        Mean mean = new Mean();
        double mu = mean.evaluate(colData.toArray());
        //LOG.info("mu is nan?" + (mu == Double.NaN? "yes":"no"));
        stats.setMean(mu);
        statistics[i] = stats;
    }
}

From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java

@Override
public GaussianFactor product(GaussianFactor other) {
    Scope newScope = scope.union(other.getVariables());

    int newFactorSize = newScope.size();

    int[] thisMapping = newScope.createContinuousVariableMapping(scope);
    int[] otherMapping = newScope.createContinuousVariableMapping(other.getVariables());

    RealMatrix newPrecisionMatrix = new Array2DRowRealMatrix(newScope.size(), newScope.size());

    RealMatrix otherPrecisionMatrix = other.getPrecisionMatrix();

    for (int i = 0; i < newFactorSize; i++) {
        RealVector column = new ArrayRealVector(newFactorSize);
        if (thisMapping[i] >= 0) {
            column = column.add(/*from ww w  .j av a2  s  .  c om*/
                    padVector(precisionMatrix.getColumnVector(thisMapping[i]), newFactorSize, thisMapping));
        }

        if (otherMapping[i] >= 0) {
            column = column.add(padVector(otherPrecisionMatrix.getColumnVector(otherMapping[i]), newFactorSize,
                    otherMapping));
        }
        newPrecisionMatrix.setColumnVector(i, column);
    }

    RealVector newScaledMeanVector = padVector(scaledMeanVector, newScope.size(), thisMapping);
    RealVector otherScaledMeanVector = other.getScaledMeanVector();
    newScaledMeanVector = newScaledMeanVector
            .add(padVector(otherScaledMeanVector, newFactorSize, otherMapping));

    double newNormalizationConstant = normalizationConstant + other.getNormalizationConstant();

    return new CanonicalGaussianFactor(newScope, newPrecisionMatrix, newScaledMeanVector,
            newNormalizationConstant);
}

From source file:edu.cmu.tetrad.data.CovarianceMatrixOnTheFly.java

/**
 * Constructs a new covariance matrix from the given data set.
 *
 * @throws IllegalArgumentException if this is not a continuous data set.
 *//* www  . j  a  v  a  2 s .  c om*/
public CovarianceMatrixOnTheFly(DataSet dataSet) {
    if (!dataSet.isContinuous()) {
        throw new IllegalArgumentException("Not a continuous data set.");
    }

    this.variables = Collections.unmodifiableList(dataSet.getVariables());
    this.sampleSize = dataSet.getNumRows();

    if (dataSet instanceof BoxDataSet) {

        DataBox box = ((BoxDataSet) dataSet).getDataBox();

        if (box instanceof VerticalDoubleDataBox) {
            if (!dataSet.getVariables().equals(variables))
                throw new IllegalArgumentException();

            vectors = ((VerticalDoubleDataBox) box).getVariableVectors();

            //                final TetradMatrix doubleData = dataSet.getDoubleData();
            TetradVector means = DataUtils.means(vectors);
            DataUtils.demean(vectors, means);

            //                DataUtils.remean(doubleData, means);
        }

    }

    if (vectors == null) {
        final TetradMatrix doubleData = dataSet.getDoubleData();
        TetradVector means = DataUtils.means(doubleData);
        DataUtils.demean(doubleData, means);

        final RealMatrix realMatrix = doubleData.getRealMatrix();

        vectors = new double[variables.size()][];

        for (int i = 0; i < variables.size(); i++) {
            vectors[i] = realMatrix.getColumnVector(i).toArray();
        }

        DataUtils.remean(doubleData, means);
    }

    this.variances = new double[variables.size()];

    class VarianceTask extends RecursiveTask<Boolean> {
        private int chunk;
        private int from;
        private int to;

        public VarianceTask(int chunk, int from, int to) {
            this.chunk = chunk;
            this.from = from;
            this.to = to;
        }

        @Override
        protected Boolean compute() {
            if (to - from <= chunk) {
                for (int i = from; i < to; i++) {
                    double d = 0.0D;

                    double[] v1 = vectors[i];

                    for (int k = 0; k < sampleSize; ++k) {
                        d += v1[k] * v1[k];
                    }

                    double v = d;
                    v /= sampleSize;

                    variances[i] = v;

                    if (v == 0)
                        System.out.println("Zero variance! " + variables.get(i));
                }

                return true;
            } else {
                final int numIntervals = 4;

                int step = (to - from) / numIntervals + 1;

                List<VarianceTask> tasks = new ArrayList<VarianceTask>();

                for (int i = 0; i < numIntervals; i++) {
                    tasks.add(new VarianceTask(chunk, from + i * step, Math.min(from + (i + 1) * step, to)));
                }

                invokeAll(tasks);

                return true;
            }
        }
    }

    int _chunk = variables.size() / NTHREADS + 1;
    final int chunk = _chunk < minChunk ? minChunk : _chunk;

    ForkJoinPoolInstance.getInstance().getPool().invoke(new VarianceTask(chunk, 0, variables.size()));

    System.out.println("Done with variances.");

}

From source file:edu.cmu.tetrad.data.CovarianceMatrix.java

/**
 * Constructs a new covariance matrix from the given data set.
 *
 * @throws IllegalArgumentException if this is not a continuous data set.
 *//* w ww. ja va 2  s  .c o m*/
public CovarianceMatrix(DataSet dataSet) {
    if (!dataSet.isContinuous()) {
        throw new IllegalArgumentException("Not a continuous data set.");
    }

    this.matrix = new TetradMatrix(dataSet.getNumColumns(), dataSet.getNumColumns());

    this.variables = Collections.unmodifiableList(dataSet.getVariables());
    this.sampleSize = dataSet.getNumRows();

    if (dataSet instanceof BoxDataSet) {

        DataBox box = ((BoxDataSet) dataSet).getDataBox().copy();

        if (box instanceof VerticalDoubleDataBox) {
            if (!dataSet.getVariables().equals(variables))
                throw new IllegalArgumentException();

            vectors = ((VerticalDoubleDataBox) box).getVariableVectors();

            //                final TetradMatrix doubleData = dataSet.getDoubleData();

            //                DataUtils.remean(doubleData, means);
        }

    }

    if (vectors == null) {
        final TetradMatrix doubleData = dataSet.getDoubleData().copy();
        TetradVector means = DataUtils.means(doubleData);
        DataUtils.demean(doubleData, means);

        final RealMatrix realMatrix = doubleData.getRealMatrix();

        vectors = new double[variables.size()][];

        for (int i = 0; i < variables.size(); i++) {
            vectors[i] = realMatrix.getColumnVector(i).toArray();
        }
    }

    TetradVector means = DataUtils.means(vectors);

    int NTHREADS = Runtime.getRuntime().availableProcessors() * 10;
    int _chunk = variables.size() / NTHREADS + 1;
    int minChunk = 100;
    final int chunk = _chunk < minChunk ? minChunk : _chunk;

    class VarianceTask extends RecursiveTask<Boolean> {
        private int chunk;
        private int from;
        private int to;

        public VarianceTask(int chunk, int from, int to) {
            this.chunk = chunk;
            this.from = from;
            this.to = to;
        }

        @Override
        protected Boolean compute() {
            if (to - from <= chunk) {
                for (int i = from; i < to; i++) {
                    double d = 0.0D;

                    int count = 0;

                    double[] v1 = vectors[i];

                    for (int k = 0; k < sampleSize; ++k) {
                        if (Double.isNaN(v1[k])) {
                            continue;
                        }

                        d += v1[k] * v1[k];
                        count++;
                    }

                    double v = d;
                    v /= (count - 1);

                    matrix.set(i, i, v);

                    if (v == 0) {
                        System.out.println("Zero variance! " + variables.get(i));
                    }
                }

                return true;
            } else {
                int mid = (to + from) / 2;

                VarianceTask left = new VarianceTask(chunk, from, mid);
                VarianceTask right = new VarianceTask(chunk, mid, to);

                left.fork();
                right.compute();
                left.join();

                return true;
            }
        }
    }

    class RestOfThemTask extends RecursiveTask<Boolean> {
        private int chunk;
        private int from;
        private int to;

        public RestOfThemTask(int chunk, int from, int to) {
            this.chunk = chunk;
            this.from = from;
            this.to = to;
        }

        @Override
        protected Boolean compute() {
            if (to - from <= chunk) {
                for (int i = from; i < to; i++) {
                    for (int j = 0; j < i; j++) {

                        double d = 0.0D;

                        double[] v1 = vectors[i];
                        double[] v2 = vectors[j];
                        int count = 0;

                        for (int k = 0; k < sampleSize; k++) {
                            if (Double.isNaN(v1[k]))
                                continue;
                            if (Double.isNaN(v2[k]))
                                continue;

                            d += v1[k] * v2[k];
                            count++;
                        }

                        double v = d;
                        v /= (count - 1);

                        matrix.set(i, j, v);
                        matrix.set(j, i, v);
                    }
                }

                return true;
            } else {
                int mid = (to + from) / 2;

                RestOfThemTask left = new RestOfThemTask(chunk, from, mid);
                RestOfThemTask right = new RestOfThemTask(chunk, mid, to);

                left.fork();
                right.compute();
                left.join();

                return true;
            }
        }
    }

    VarianceTask task = new VarianceTask(chunk, 0, variables.size());
    ForkJoinPoolInstance.getInstance().getPool().invoke(task);

    RestOfThemTask task2 = new RestOfThemTask(chunk, 0, variables.size());
    ForkJoinPoolInstance.getInstance().getPool().invoke(task2);

    DataUtils.demean(vectors, means);

    this.variables = Collections.unmodifiableList(dataSet.getVariables());
    this.sampleSize = dataSet.getNumRows();
}

From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java

@Test
public void testSetColumnVector() {
    RealMatrix m = new BigSparseRealMatrix(subTestData);
    RealVector mColumn3 = columnToVector(subColumn3);
    Assert.assertNotSame(mColumn3, m.getColumnVector(1));
    m.setColumnVector(1, mColumn3);/* www. j av  a2s.co m*/
    Assert.assertEquals(mColumn3, m.getColumnVector(1));
    try {
        m.setColumnVector(-1, mColumn3);
        Assert.fail("Expecting OutOfRangeException");
    } catch (OutOfRangeException ex) {
        // expected
    }
    try {
        m.setColumnVector(0, new ArrayRealVector(5));
        Assert.fail("Expecting MatrixDimensionMismatchException");
    } catch (MatrixDimensionMismatchException ex) {
        // expected
    }
}

From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java

@Test
public void testGetColumnVector() {
    RealMatrix m = new BigSparseRealMatrix(subTestData);
    RealVector mColumn1 = columnToVector(subColumn1);
    RealVector mColumn3 = columnToVector(subColumn3);
    Assert.assertEquals("Column1", mColumn1, m.getColumnVector(1));
    Assert.assertEquals("Column3", mColumn3, m.getColumnVector(3));
    try {// w w w . ja v  a2 s .  c o  m
        m.getColumnVector(-1);
        Assert.fail("Expecting OutOfRangeException");
    } catch (OutOfRangeException ex) {
        // expected
    }
    try {
        m.getColumnVector(4);
        Assert.fail("Expecting OutOfRangeException");
    } catch (OutOfRangeException ex) {
        // expected
    }
}

From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java

/**
* Creates a correction from a set of objects whose positions should be the same in each channel.
* 
* @param imageObjects                  A Vector containing all the ImageObjects to be used for the correction
*                                      or in the order it appears in a multiwavelength image file.
* @return                              A Correction object that can be used to correct the positions of other objects based upon the standards provided.
*//*from w ww  . ja v  a 2  s . co  m*/
public Correction getCorrection(java.util.List<ImageObject> imageObjects) {

    int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM);

    int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM);

    if (!this.parameters.hasKeyAndTrue(DET_CORR_PARAM)) {
        try {
            return Correction.readFromDisk(FileUtils.getCorrectionFilename(this.parameters));
        } catch (java.io.IOException e) {

            java.util.logging.Logger
                    .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME)
                    .severe("Exception encountered while reading correction from disk: ");
            e.printStackTrace();

        } catch (ClassNotFoundException e) {

            java.util.logging.Logger
                    .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME)
                    .severe("Exception encountered while reading correction from disk: ");
            e.printStackTrace();

        }

        return null;
    }

    int numberOfPointsToFit = this.parameters.getIntValueForKey(NUM_POINT_PARAM);

    RealMatrix correctionX = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters);
    RealMatrix correctionY = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters);
    RealMatrix correctionZ = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters);

    RealVector distanceCutoffs = new ArrayRealVector(imageObjects.size(), 0.0);

    RealVector ones = new ArrayRealVector(numberOfPointsToFit, 1.0);

    RealVector distancesToObjects = new ArrayRealVector(imageObjects.size(), 0.0);

    RealMatrix allCorrectionParametersMatrix = new Array2DRowRealMatrix(numberOfPointsToFit,
            numberOfCorrectionParameters);

    for (int i = 0; i < imageObjects.size(); i++) {

        RealVector ithPos = imageObjects.get(i).getPositionForChannel(referenceChannel);

        for (int j = 0; j < imageObjects.size(); j++) {

            double d = imageObjects.get(j).getPositionForChannel(referenceChannel).subtract(ithPos).getNorm();

            distancesToObjects.setEntry(j, d);

        }

        //the sorting becomes a bottleneck once the number of points gets large

        //reverse comparator so we can use the priority queue and get the max element at the head

        Comparator<Double> cdReverse = new Comparator<Double>() {

            public int compare(Double o1, Double o2) {

                if (o1.equals(o2))
                    return 0;
                if (o1 > o2)
                    return -1;
                return 1;
            }

        };

        PriorityQueue<Double> pq = new PriorityQueue<Double>(numberOfPointsToFit + 2, cdReverse);

        double maxElement = Double.MAX_VALUE;

        for (int p = 0; p < numberOfPointsToFit + 1; p++) {

            pq.add(distancesToObjects.getEntry(p));

        }

        maxElement = pq.peek();

        for (int p = numberOfPointsToFit + 1; p < distancesToObjects.getDimension(); p++) {

            double value = distancesToObjects.getEntry(p);

            if (value < maxElement) {

                pq.poll();

                pq.add(value);

                maxElement = pq.peek();

            }

        }

        double firstExclude = pq.poll();
        double lastDist = pq.poll();

        double distanceCutoff = (lastDist + firstExclude) / 2.0;

        distanceCutoffs.setEntry(i, distanceCutoff);

        RealVector xPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0);
        RealVector yPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0);
        RealVector zPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0);

        RealMatrix differencesToFit = new Array2DRowRealMatrix(numberOfPointsToFit,
                imageObjects.get(0).getPositionForChannel(referenceChannel).getDimension());

        int toFitCounter = 0;

        for (int j = 0; j < imageObjects.size(); j++) {
            if (distancesToObjects.getEntry(j) < distanceCutoff) {
                xPositionsToFit.setEntry(toFitCounter,
                        imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(0));
                yPositionsToFit.setEntry(toFitCounter,
                        imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(1));
                zPositionsToFit.setEntry(toFitCounter,
                        imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(2));

                differencesToFit.setRowVector(toFitCounter, imageObjects.get(j)
                        .getVectorDifferenceBetweenChannels(referenceChannel, channelToCorrect));

                toFitCounter++;
            }
        }

        RealVector x = xPositionsToFit.mapSubtractToSelf(ithPos.getEntry(0));
        RealVector y = yPositionsToFit.mapSubtractToSelf(ithPos.getEntry(1));

        allCorrectionParametersMatrix.setColumnVector(0, ones);
        allCorrectionParametersMatrix.setColumnVector(1, x);
        allCorrectionParametersMatrix.setColumnVector(2, y);
        allCorrectionParametersMatrix.setColumnVector(3, x.map(new Power(2)));
        allCorrectionParametersMatrix.setColumnVector(4, y.map(new Power(2)));
        allCorrectionParametersMatrix.setColumnVector(5, x.ebeMultiply(y));

        DecompositionSolver solver = (new QRDecomposition(allCorrectionParametersMatrix)).getSolver();

        RealVector cX = solver.solve(differencesToFit.getColumnVector(0));
        RealVector cY = solver.solve(differencesToFit.getColumnVector(1));
        RealVector cZ = solver.solve(differencesToFit.getColumnVector(2));

        correctionX.setRowVector(i, cX);
        correctionY.setRowVector(i, cY);
        correctionZ.setRowVector(i, cZ);

    }

    Correction c = new Correction(correctionX, correctionY, correctionZ, distanceCutoffs, imageObjects,
            referenceChannel, channelToCorrect);

    return c;

}

From source file:org.akvo.caddisfly.sensor.colorimetry.strip.calibration.CalibrationCard.java

@NonNull
private static Mat doIlluminationCorrection(@NonNull Mat imgLab, @NonNull CalibrationData calData) {
    // create HLS image for homogeneous illumination calibration
    int pHeight = imgLab.rows();
    int pWidth = imgLab.cols();

    RealMatrix points = createWhitePointMatrix(imgLab, calData);

    // create coefficient matrix for all three variables L,A,B
    // the model for all three is y = ax + bx^2 + cy + dy^2 + exy + f
    // 6th row is the constant 1
    RealMatrix coefficient = new Array2DRowRealMatrix(points.getRowDimension(), 6);
    coefficient.setColumnMatrix(0, points.getColumnMatrix(0));
    coefficient.setColumnMatrix(2, points.getColumnMatrix(1));

    //create constant, x^2, y^2 and xy terms
    for (int i = 0; i < points.getRowDimension(); i++) {
        coefficient.setEntry(i, 1, Math.pow(coefficient.getEntry(i, 0), 2)); // x^2
        coefficient.setEntry(i, 3, Math.pow(coefficient.getEntry(i, 2), 2)); // y^2
        coefficient.setEntry(i, 4, coefficient.getEntry(i, 0) * coefficient.getEntry(i, 2)); // xy
        coefficient.setEntry(i, 5, 1d); // constant = 1
    }/*from   w  ww. ja v  a  2  s .  com*/

    // create vectors
    RealVector L = points.getColumnVector(2);
    RealVector A = points.getColumnVector(3);
    RealVector B = points.getColumnVector(4);

    // solve the least squares problem for all three variables
    DecompositionSolver solver = new SingularValueDecomposition(coefficient).getSolver();
    RealVector solutionL = solver.solve(L);
    RealVector solutionA = solver.solve(A);
    RealVector solutionB = solver.solve(B);

    // get individual coefficients
    float La = (float) solutionL.getEntry(0);
    float Lb = (float) solutionL.getEntry(1);
    float Lc = (float) solutionL.getEntry(2);
    float Ld = (float) solutionL.getEntry(3);
    float Le = (float) solutionL.getEntry(4);
    float Lf = (float) solutionL.getEntry(5);

    float Aa = (float) solutionA.getEntry(0);
    float Ab = (float) solutionA.getEntry(1);
    float Ac = (float) solutionA.getEntry(2);
    float Ad = (float) solutionA.getEntry(3);
    float Ae = (float) solutionA.getEntry(4);
    float Af = (float) solutionA.getEntry(5);

    float Ba = (float) solutionB.getEntry(0);
    float Bb = (float) solutionB.getEntry(1);
    float Bc = (float) solutionB.getEntry(2);
    float Bd = (float) solutionB.getEntry(3);
    float Be = (float) solutionB.getEntry(4);
    float Bf = (float) solutionB.getEntry(5);

    // compute mean (the luminosity value of the plane in the middle of the image)
    float L_mean = (float) (0.5 * La * pWidth + 0.5 * Lc * pHeight + Lb * pWidth * pWidth / 3.0
            + Ld * pHeight * pHeight / 3.0 + Le * 0.25 * pHeight * pWidth + Lf);
    float A_mean = (float) (0.5 * Aa * pWidth + 0.5 * Ac * pHeight + Ab * pWidth * pWidth / 3.0
            + Ad * pHeight * pHeight / 3.0 + Ae * 0.25 * pHeight * pWidth + Af);
    float B_mean = (float) (0.5 * Ba * pWidth + 0.5 * Bc * pHeight + Bb * pWidth * pWidth / 3.0
            + Bd * pHeight * pHeight / 3.0 + Be * 0.25 * pHeight * pWidth + Bf);

    // Correct image
    // we do this per row. We tried to do it in one block, but there is no speed difference.
    byte[] temp = new byte[imgLab.cols() * imgLab.channels()];
    int valL, valA, valB;
    int ii, ii3;
    float iiSq, iSq;
    int imgCols = imgLab.cols();
    int imgRows = imgLab.rows();

    // use lookup tables to speed up computation
    // create lookup tables
    float[] L_aii = new float[imgCols];
    float[] L_biiSq = new float[imgCols];
    float[] A_aii = new float[imgCols];
    float[] A_biiSq = new float[imgCols];
    float[] B_aii = new float[imgCols];
    float[] B_biiSq = new float[imgCols];

    float[] Lci = new float[imgRows];
    float[] LdiSq = new float[imgRows];
    float[] Aci = new float[imgRows];
    float[] AdiSq = new float[imgRows];
    float[] Bci = new float[imgRows];
    float[] BdiSq = new float[imgRows];

    for (ii = 0; ii < imgCols; ii++) {
        iiSq = ii * ii;
        L_aii[ii] = La * ii;
        L_biiSq[ii] = Lb * iiSq;
        A_aii[ii] = Aa * ii;
        A_biiSq[ii] = Ab * iiSq;
        B_aii[ii] = Ba * ii;
        B_biiSq[ii] = Bb * iiSq;
    }

    for (int i = 0; i < imgRows; i++) {
        iSq = i * i;
        Lci[i] = Lc * i;
        LdiSq[i] = Ld * iSq;
        Aci[i] = Ac * i;
        AdiSq[i] = Ad * iSq;
        Bci[i] = Bc * i;
        BdiSq[i] = Bd * iSq;
    }

    // We can also improve the performance of the i,ii term, if we want, but it won't make much difference.
    for (int i = 0; i < imgRows; i++) { // y
        imgLab.get(i, 0, temp);
        ii3 = 0;
        for (ii = 0; ii < imgCols; ii++) { //x
            valL = capValue(
                    Math.round((temp[ii3] & 0xFF)
                            - (L_aii[ii] + L_biiSq[ii] + Lci[i] + LdiSq[i] + Le * i * ii + Lf) + L_mean),
                    0, 255);
            valA = capValue(
                    Math.round((temp[ii3 + 1] & 0xFF)
                            - (A_aii[ii] + A_biiSq[ii] + Aci[i] + AdiSq[i] + Ae * i * ii + Af) + A_mean),
                    0, 255);
            valB = capValue(
                    Math.round((temp[ii3 + 2] & 0xFF)
                            - (B_aii[ii] + B_biiSq[ii] + Bci[i] + BdiSq[i] + Be * i * ii + Bf) + B_mean),
                    0, 255);

            temp[ii3] = (byte) valL;
            temp[ii3 + 1] = (byte) valA;
            temp[ii3 + 2] = (byte) valB;
            ii3 += 3;
        }
        imgLab.put(i, 0, temp);
    }

    return imgLab;
}