Java tutorial
/* * Copyright (C) 2013 Dr. John Lindsay <jlindsay@uoguelph.ca> * * 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, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package whitebox.stats; import java.util.ArrayList; import org.apache.commons.math3.linear.*; import whitebox.structures.XYPoint; /** * * @author johnlindsay */ public class PolynomialLeastSquares2DFitting { private int polyOrder = 1; private double[] forwardRegressCoeffX; private double[] forwardRegressCoeffY; private double[] backRegressCoeffX; private double[] backRegressCoeffY; private int numCoefficients; private double[] xCoords1; private double[] yCoords1; private double[] xCoords2; private double[] yCoords2; private double[] residualsXY; private double[] residualsOrientation; //private boolean[] useGCP; private double xMin1; private double yMin1; private double xMin2; private double yMin2; private double overallRMSE = 0.0; public PolynomialLeastSquares2DFitting() { } public PolynomialLeastSquares2DFitting(double[] X1, double[] Y1, double[] X2, double[] Y2, int polyOrder) { this.polyOrder = polyOrder; addData(X1, Y1, X2, Y2); } public PolynomialLeastSquares2DFitting(ArrayList<XYPoint> data1, ArrayList<XYPoint> data2, int polyOrder) { this.polyOrder = polyOrder; double[] X1 = new double[data1.size()]; double[] Y1 = new double[data1.size()]; double[] X2 = new double[data2.size()]; double[] Y2 = new double[data2.size()]; int i = 0; for (XYPoint xy : data1) { X1[i] = xy.x; Y1[i] = xy.y; i++; } i = 0; for (XYPoint xy : data2) { X2[i] = xy.x; Y2[i] = xy.y; i++; } addData(X1, Y1, X2, Y2); } // properties public int getPolyOrder() { return polyOrder; } public void setPolyOrder(int polyOrder) { if (polyOrder < 1) { polyOrder = 1; } if (polyOrder > 10) { polyOrder = 10; } this.polyOrder = polyOrder; } public double[] getForwardRegressCoeffX() { return forwardRegressCoeffX; } public double[] getForwardRegressCoeffY() { return forwardRegressCoeffY; } public double[] getBackRegressCoeffX() { return backRegressCoeffX; } public double[] getBackRegressCoeffY() { return backRegressCoeffY; } public double[] getResidualsXY() { return residualsXY; } public double[] getResidualsOrientation() { return residualsOrientation; } public double getOverallRMSE() { return overallRMSE; } // methods public void addData(ArrayList<XYPoint> data1, ArrayList<XYPoint> data2) { double[] X1 = new double[data1.size()]; double[] Y1 = new double[data1.size()]; double[] X2 = new double[data2.size()]; double[] Y2 = new double[data2.size()]; int i = 0; for (XYPoint xy : data1) { X1[i] = xy.x; Y1[i] = xy.y; i++; } i = 0; for (XYPoint xy : data2) { X2[i] = xy.x; Y2[i] = xy.y; i++; } addData(X1, Y1, X2, Y2); } public final void addData(double[] X1, double[] Y1, double[] X2, double[] Y2) { int n = X1.length; if (Y1.length != n || X2.length != n || Y2.length != n) { return; } xCoords1 = new double[n]; yCoords1 = new double[n]; xCoords2 = new double[n]; yCoords2 = new double[n]; System.arraycopy(X1, 0, xCoords1, 0, n); System.arraycopy(Y1, 0, yCoords1, 0, n); System.arraycopy(X2, 0, xCoords2, 0, n); System.arraycopy(Y2, 0, yCoords2, 0, n); xMin1 = Double.POSITIVE_INFINITY; yMin1 = Double.POSITIVE_INFINITY; xMin2 = Double.POSITIVE_INFINITY; yMin2 = Double.POSITIVE_INFINITY; for (int i = 0; i < n; i++) { if (X1[i] < xMin1) { xMin1 = X1[i]; } if (Y1[i] < yMin1) { yMin1 = Y1[i]; } if (X2[i] < xMin2) { xMin2 = X2[i]; } if (Y2[i] < yMin2) { yMin2 = Y2[i]; } } calculateEquations(); } public void calculateEquations() { try { int m, i, j, k; int n = xCoords2.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } // for (i = 0; i < n; i++) { // xCoords1[i] -= xMin1; // yCoords1[i] -= yMin1; // xCoords2[i] -= xMin2; // yCoords2[i] -= yMin2; // } // Solve the forward transformation equations double[][] forwardCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { forwardCoefficientMatrix[i][m] = Math.pow(xCoords1[i], j) * Math.pow(yCoords1[i], k); m++; } } } RealMatrix coefficients = new Array2DRowRealMatrix(forwardCoefficientMatrix, false); DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); //DecompositionSolver solver = new QRDecomposition(coefficients).getSolver(); // do the x-coordinate first RealVector constants = new ArrayRealVector(xCoords2, false); RealVector solution = solver.solve(constants); forwardRegressCoeffX = new double[n]; for (int a = 0; a < numCoefficients; a++) { forwardRegressCoeffX[a] = solution.getEntry(a); } double[] residualsX = new double[n]; double SSresidX = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * forwardRegressCoeffX[j]; } residualsX[i] = xCoords2[i] - yHat; SSresidX += residualsX[i] * residualsX[i]; } double sumX = 0; double SSx = 0; for (i = 0; i < n; i++) { SSx += xCoords2[i] * xCoords2[i]; sumX += xCoords2[i]; } double varianceX = (SSx - (sumX * sumX) / n) / n; double SStotalX = (n - 1) * varianceX; double rsqX = 1 - SSresidX / SStotalX; // now the y-coordinate constants = new ArrayRealVector(yCoords2, false); solution = solver.solve(constants); forwardRegressCoeffY = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { forwardRegressCoeffY[a] = solution.getEntry(a); } double[] residualsY = new double[n]; residualsXY = new double[n]; residualsOrientation = new double[n]; double SSresidY = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * forwardRegressCoeffY[j]; } residualsY[i] = yCoords2[i] - yHat; SSresidY += residualsY[i] * residualsY[i]; residualsXY[i] = Math.sqrt(residualsX[i] * residualsX[i] + residualsY[i] * residualsY[i]); residualsOrientation[i] = Math.atan2(residualsY[i], residualsX[i]); } double sumY = 0; double sumR = 0; double SSy = 0; double SSr = 0; for (i = 0; i < n; i++) { SSy += yCoords2[i] * yCoords2[i]; SSr += residualsXY[i] * residualsXY[i]; sumY += yCoords2[i]; sumR += residualsXY[i]; } double varianceY = (SSy - (sumY * sumY) / n) / n; double varianceResiduals = (SSr - (sumR * sumR) / n) / n; double SStotalY = (n - 1) * varianceY; double rsqY = 1 - SSresidY / SStotalY; overallRMSE = Math.sqrt(varianceResiduals); //System.out.println("y-coordinate r-square: " + rsqY); // // Print the residuals. // System.out.println("\nResiduals:"); // for (i = 0; i < n; i++) { // System.out.println("Point " + (i + 1) + "\t" + residualsX[i] // + "\t" + residualsY[i] + "\t" + residualsXY[i]); // } // Solve the backward transformation equations double[][] backCoefficientMatrix = new double[n][numCoefficients]; for (i = 0; i < n; i++) { m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { backCoefficientMatrix[i][m] = Math.pow(xCoords2[i], j) * Math.pow(yCoords2[i], k); m++; } } } coefficients = new Array2DRowRealMatrix(backCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); solver = new QRDecomposition(coefficients).getSolver(); // do the x-coordinate first constants = new ArrayRealVector(xCoords1, false); solution = solver.solve(constants); backRegressCoeffX = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { backRegressCoeffX[a] = solution.getEntry(a); } // now the y-coordinate constants = new ArrayRealVector(yCoords1, false); solution = solver.solve(constants); backRegressCoeffY = new double[n]; for (int a = 0; a < numCoefficients; a++) { backRegressCoeffY[a] = solution.getEntry(a); } } catch (Exception e) { e.printStackTrace(); // showFeedback("Error in ImageRectificationDialog.calculateEquations: " // + e.getMessage()); } } public XYPoint getForwardCoordinates(XYPoint point) { return getForwardCoordinates(point.x, point.y); } public XYPoint getForwardCoordinates(double x, double y) { XYPoint ret; int j, k, m; double x_transformed = 0; //mapXMin; double y_transformed = 0; //mapYMin; double term; m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { term = Math.pow(x, j) * Math.pow(y, k); x_transformed += term * forwardRegressCoeffX[m]; y_transformed += term * forwardRegressCoeffY[m]; m++; } } ret = new XYPoint(x_transformed, y_transformed); return ret; } public XYPoint getBackwardCoordinates(XYPoint point) { return getBackwardCoordinates(point.x, point.y); } public XYPoint getBackwardCoordinates(double x, double y) { XYPoint ret; int j, k, m; double x_transformed = 0; //imageXMin; double y_transformed = 0; //imageYMin; double term; m = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { term = Math.pow(x, j) * Math.pow(y, k); x_transformed += term * backRegressCoeffX[m]; y_transformed += term * backRegressCoeffY[m]; m++; } } ret = new XYPoint(x_transformed, y_transformed); return ret; } }