List of usage examples for org.apache.commons.math3.linear RealVector getEntry
public abstract double getEntry(int index) throws OutOfRangeException;
From source file:org.rhwlab.BHCnotused.StdGaussianWishartGaussian.java
private void printVector(PrintStream stream, RealVector v) { boolean first = true; stream.print("("); for (int i = 0; i < v.getDimension(); ++i) { if (!first) { stream.printf(",%d", (int) v.getEntry(i)); } else {// www .j a v a 2 s. c om stream.printf("%d", (int) v.getEntry(i)); } first = false; } stream.print(")"); }
From source file:org.rhwlab.BHCnotused.ThreadedAlgorithm.java
public void init(int seg) throws Exception { clusters = new ArrayList<>(); // build the initial clusters with one data point in each cluster for (int n = 0; n < source.getK(); ++n) { RealVector v = source.getCenter(n); Dfp[] z = new Dfp[v.getDimension()]; for (int i = 0; i < z.length; ++i) { z[i] = field.newDfp(v.getEntry(i)); }//from w ww . ja v a2 s . c o m LabeledFieldVector fv = new LabeledFieldVector(z, n); Cluster cluster = new Cluster(new GaussianGIWPrior(fv)); clusters.add(cluster); } // make all possible pairings of initial clusters pairs = new HashMap<>(); for (int i = 0; i < clusters.size() - 1; ++i) { HashMap<Cluster, Cluster> map = new HashMap<>(); MergeAction merge = new MergeAction(clusters, i, i + 1, map); ForkJoinPool pool = new ForkJoinPool(); pool.invoke(merge); Cluster clusterI = clusters.get(i); pairs.put(clusterI, map); System.out.printf("Cluster %d paired\n", i); } }
From source file:org.rhwlab.dispim.datasource.Segmentation.java
public Segmentation(VoxelDataSource source, double thresh, BoundingBox box) { this.source = source; this.thresh = thresh; segmentIndex = new ArrayList<>(); mins = new double[source.getD()]; maxs = new double[source.getD()]; for (int d = 0; d < source.getD(); ++d) { mins[d] = Double.MAX_VALUE; maxs[d] = 0.0;/*from w ww. j ava2s . c om*/ } for (long i = 0; i < source.getN(); ++i) { // process each voxel in the segmented tiff Voxel segVox = source.get(i); if (segVox.getIntensity() >= thresh && box.isWithin(segVox)) { // do not include the voxel in the segmentation if less than the threshold or outside the bounding box for (int d = 0; d < mins.length; ++d) { RealVector v = segVox.coords; double e = v.getEntry(d); if (e < mins[d]) { mins[d] = e; } if (e > maxs[d]) { maxs[d] = e; } } segmentIndex.add(i); } } }
From source file:org.rhwlab.dispim.nucleus.BHCNucleusData.java
public String getRadiusLabel(int i) { // find the adjusted eigenvector closest to the original unadjusted eigenvector for dimension i // this keeps the order of the adjusted eigenvectors the same as the original unadjusted eigenvectors // the eigendecomposition returns the eigenvectors sorted by eigenvalue // this procedure puts them back in their original order int adjustedI = 0; double maxD = 0.0; RealVector aV = eigenA.getEigenvector(i); for (int j = 0; j < A.getColumnDimension(); ++j) { RealVector v = adjustedEigenA.getEigenvector(j); double d = Math.abs(v.dotProduct(aV)); if (d > maxD) { maxD = d;//w w w . j a v a 2s . com adjustedI = j; } } RealVector v = adjustedEigenA.getEigenvector(adjustedI); // double eigenVal = adjustedEigenA.getRealEigenvalue(adjustedI); // double r = 1.0/Math.sqrt(Ace3D_Frame.R*eigenVal); double r = this.getRadius(adjustedI); return String.format("%4.1f(%.2f,%.2f,%.2f)", r, v.getEntry(0), v.getEntry(1), v.getEntry(2)); }
From source file:org.rhwlab.dispim.nucleus.NucleusData.java
public Ellipse2d ellipse(int xi, int yi, int zi, Coeff coef) { Array2DRowRealMatrix Q = new Array2DRowRealMatrix(3, 3); Q.setEntry(0, 0, coef.A);//from w w w . j av a 2 s .com Q.setEntry(1, 0, coef.B / 2); Q.setEntry(0, 1, coef.B / 2); Q.setEntry(1, 1, coef.C); Q.setEntry(2, 0, coef.D / 2); Q.setEntry(0, 2, coef.D / 2); Q.setEntry(2, 1, coef.E / 2); Q.setEntry(1, 2, coef.E / 2); Q.setEntry(2, 2, coef.F); EigenDecomposition ed = new EigenDecomposition(Q); double detQ = ed.getDeterminant(); RealMatrix rm = new Array2DRowRealMatrix(2, 2); rm.setEntry(0, 0, coef.A); rm.setEntry(1, 1, coef.C); rm.setEntry(0, 1, coef.B / 2.0); rm.setEntry(1, 0, coef.B / 2.0); EigenDecomposition eigenDecomp = new EigenDecomposition(rm); double detA33 = eigenDecomp.getDeterminant(); double[] eigenValues = eigenDecomp.getRealEigenvalues(); //System.out.printf("Eigenvalues: %f,%f\n",eigenValues[0],eigenValues[1]) ; RealVector eigenvector0 = eigenDecomp.getEigenvector(0); Ellipse2d e = new Ellipse2d(); double cot2theta = (coef.A - coef.C) / coef.B; if (Double.isFinite(cot2theta)) { double d = Math.sqrt(1.0 + cot2theta * cot2theta); double cos2theta = cot2theta / d; e.cosine = Math.sqrt((1.0 + cos2theta) / 2.0); e.sine = Math.sqrt((1.0 - cos2theta) / 2.0); } else { e.cosine = 1.0; e.sine = 0.0; } double dd = coef.B * coef.B - 4.0 * coef.A * coef.C; double xc = (2.0 * coef.C * coef.D - coef.B * coef.E) / dd; double yc = (2.0 * coef.A * coef.E - coef.B * coef.D) / dd; // System.out.printf("dd=%f,xc=%f,xcn=%f,yc=%f,ycn%f\n",dd,xc,xcn,yc,ycn); double[] ce = this.getCenter(); e.x = ce[xi] + xc; e.y = ce[yi] + yc; // double f = -detQ/detA33 * 2; double f = -detQ / detA33; double a = eigenValues[0] / f; double b = eigenValues[1] / f; if (a <= 0.0 || b <= 0.0) { return null; } e.a = 1.0 / Math.sqrt(a); e.b = 1.0 / Math.sqrt(b); //System.out.printf("detQ=%e,detA33=%e,f=%f,a=%e,b=%e\n",detQ,detA33,f,e.a,e.b); //System.out.printf("eigenValues (%f,%f)\n",eigenValues[0],eigenValues[1]); e.cosine = eigenvector0.getEntry(0); e.sine = eigenvector0.getEntry(1); e.low[xi] = (long) (e.x - e.a); e.low[yi] = (long) (e.y - e.b); e.low[zi] = 0; e.high[xi] = (long) (e.x + e.a); e.high[yi] = (long) (e.y + e.b); e.high[zi] = 0; //System.out.printf("Ellipse: a=%f,b=%f,x=%f,y=%f\n", e.a,e.b,e.x,e.y); return e; }
From source file:playground.sergioo.facilitiesGenerator2012.WorkFacilitiesGeneration.java
private static Set<PointPerson> getPCATransformation(Collection<PointPerson> points) { RealMatrix pointsM = new Array2DRowRealMatrix(points.iterator().next().getDimension(), points.size()); int k = 0;// ww w.j a v a 2s . c o m for (PointND<Double> point : points) { for (int f = 0; f < point.getDimension(); f++) pointsM.setEntry(f, k, point.getElement(f)); k++; } RealMatrix means = new Array2DRowRealMatrix(pointsM.getRowDimension(), 1); for (int r = 0; r < means.getRowDimension(); r++) { double mean = 0; for (int c = 0; c < pointsM.getColumnDimension(); c++) mean += pointsM.getEntry(r, c) / pointsM.getColumnDimension(); means.setEntry(r, 0, mean); } RealMatrix deviations = new Array2DRowRealMatrix(pointsM.getRowDimension(), pointsM.getColumnDimension()); for (int r = 0; r < deviations.getRowDimension(); r++) for (int c = 0; c < deviations.getColumnDimension(); c++) deviations.setEntry(r, c, pointsM.getEntry(r, c) - means.getEntry(r, 0)); RealMatrix covariance = deviations.multiply(deviations.transpose()) .scalarMultiply(1 / (double) pointsM.getColumnDimension()); EigenDecomposition eigenDecomposition = new EigenDecomposition(covariance, 0); RealMatrix eigenVectorsT = eigenDecomposition.getVT(); RealVector eigenValues = new ArrayRealVector(eigenDecomposition.getD().getRowDimension()); for (int r = 0; r < eigenDecomposition.getD().getRowDimension(); r++) eigenValues.setEntry(r, eigenDecomposition.getD().getEntry(r, r)); for (int i = 0; i < eigenValues.getDimension(); i++) { for (int j = i + 1; j < eigenValues.getDimension(); j++) if (eigenValues.getEntry(i) < eigenValues.getEntry(j)) { double tempValue = eigenValues.getEntry(i); eigenValues.setEntry(i, eigenValues.getEntry(j)); eigenValues.setEntry(j, tempValue); RealVector tempVector = eigenVectorsT.getRowVector(i); eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(j)); eigenVectorsT.setRowVector(j, tempVector); } eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(i).mapMultiply(Math.sqrt(1 / eigenValues.getEntry(i)))); } RealVector standardDeviations = new ArrayRealVector(pointsM.getRowDimension()); for (int r = 0; r < covariance.getRowDimension(); r++) standardDeviations.setEntry(r, Math.sqrt(covariance.getEntry(r, r))); double zValue = standardDeviations.dotProduct(new ArrayRealVector(pointsM.getRowDimension(), 1)); RealMatrix zScore = deviations.scalarMultiply(1 / zValue); pointsM = eigenVectorsT.multiply(zScore); Set<PointPerson> pointsC = new HashSet<PointPerson>(); k = 0; for (PointPerson point : points) { PointPerson pointC = new PointPerson(point.getId(), point.getOccupation(), new Double[] { pointsM.getEntry(0, k), pointsM.getEntry(1, k) }, point.getPlaceType()); pointC.setWeight(point.getWeight()); pointsC.add(pointC); k++; } return pointsC; }
From source file:plugins.ImageRectificationPanel.java
public void calculateEquations(double[] imageX, double[] imageY, double[] mapX, double[] mapY) { try {// w w w . j a v a2s. com int m, i, j, k; int n = mapX.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++) { imageX[i] -= imageXMin; imageY[i] -= imageYMin; mapX[i] -= mapXMin; mapY[i] -= mapYMin; } // 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(imageX[i], j) * Math.pow(imageY[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(mapX, 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] = mapX[i] - yHat; SSresidX += residualsX[i] * residualsX[i]; } double sumX = 0; double SSx = 0; for (i = 0; i < n; i++) { SSx += mapX[i] * mapX[i]; sumX += mapX[i]; } double varianceX = (SSx - (sumX * sumX) / n) / n; double SStotalX = (n - 1) * varianceX; double rsqX = 1 - SSresidX / SStotalX; //System.out.println("x-coordinate r-square: " + rsqX); // now the y-coordinate constants = new ArrayRealVector(mapY, 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]; 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] = mapY[i] - yHat; SSresidY += residualsY[i] * residualsY[i]; residualsXY[i] = Math.sqrt(residualsX[i] * residualsX[i] + residualsY[i] * residualsY[i]); } double sumY = 0; double sumR = 0; double SSy = 0; double SSr = 0; for (i = 0; i < n; i++) { SSy += mapY[i] * mapY[i]; SSr += residualsXY[i] * residualsXY[i]; sumY += mapY[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(mapX[i], j) * Math.pow(mapY[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(imageX, 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(imageY, false); solution = solver.solve(constants); backRegressCoeffY = new double[n]; for (int a = 0; a < numCoefficients; a++) { backRegressCoeffY[a] = solution.getEntry(a); } } catch (OutOfMemoryError oe) { myHost.showFeedback("An out-of-memory error has occurred during operation."); } catch (Exception e) { myHost.showFeedback("An error has occurred during operation. See log file for details."); myHost.logException("Error in ImageRectification", e); } }
From source file:plugins.SURFPixelMatching.java
private void run() { try {/*from w w w . j a v a 2 s .c o m*/ // variables int a, b, c, i, j, k, n, r; int row, col; int progress, oldProgress; double x, y, z, newX, newY; double x2, y2; double north, south, east, west; double newNorth, newSouth, newEast, newWest; double rightNodata; double leftNodata; Object[] rowData; whitebox.geospatialfiles.shapefile.Point outputPoint; ShapeFile rightTiePoints; ShapeFile leftTiePoints; ShapeFile rightFiducials; ShapeFile leftFiducials; XYPoint xyPoint; ArrayList<XYPoint> leftTiePointsList = new ArrayList<>(); ArrayList<XYPoint> rightTiePointsList = new ArrayList<>(); float balanceValue = 0.81f; float threshold = 0.004f; int octaves = 4; double maxAllowableRMSE = 1.0; double matchingThreshold = 0.6; // left image //String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus_C6430-74072-L9_253_Blue_clipped.dep"; //String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-82_Blue.dep"; //String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus 253.dep"; //String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus 253 epipolar.dep"; String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-82_Blue low res.dep"; //String leftImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/tmp1.dep"; // right image //String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus_C6430-74072-L9_254_Blue_clipped.dep"; //String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-83_Blue.dep"; //String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus 254.dep"; //String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus 254 epipolar.dep"; String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-83_Blue low res.dep"; //String rightImageName = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus 253.dep"; String leftOutputHeader = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/tmp4.shp"; String rightOutputHeader = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/tmp4_2.shp"; DBFField[] fields = new DBFField[5]; fields[0] = new DBFField(); fields[0].setName("FID"); fields[0].setDataType(DBFField.DBFDataType.NUMERIC); fields[0].setDecimalCount(4); fields[0].setFieldLength(10); fields[1] = new DBFField(); fields[1].setName("ORIENT"); fields[1].setDataType(DBFField.DBFDataType.NUMERIC); fields[1].setDecimalCount(4); fields[1].setFieldLength(10); fields[2] = new DBFField(); fields[2].setName("SCALE"); fields[2].setDataType(DBFField.DBFDataType.NUMERIC); fields[2].setDecimalCount(4); fields[2].setFieldLength(10); fields[3] = new DBFField(); fields[3].setName("LAPLAC"); fields[3].setDataType(DBFField.DBFDataType.NUMERIC); fields[3].setDecimalCount(4); fields[3].setFieldLength(10); fields[4] = new DBFField(); fields[4].setName("RESID"); fields[4].setDataType(DBFField.DBFDataType.NUMERIC); fields[4].setDecimalCount(4); fields[4].setFieldLength(10); // read the input data WhiteboxRaster leftImage = new WhiteboxRaster(leftImageName, "r"); // leftImage.setForceAllDataInMemory(true); int nRowsLeft = leftImage.getNumberRows(); int nColsLeft = leftImage.getNumberColumns(); leftNodata = leftImage.getNoDataValue(); WhiteboxRaster rightImage = new WhiteboxRaster(rightImageName, "r"); //rightImage.setForceAllDataInMemory(true); // int nRowsRight = rightImage.getNumberRows(); // int nColsRight = rightImage.getNumberColumns(); rightNodata = rightImage.getNoDataValue(); // ArrayList<InterestPoint> interest_points; // double threshold = 600; // double balanceValue = 0.9; // int octaves = 4; // ISURFfactory mySURF = SURF.createInstance(leftImage, balanceValue, threshold, octaves); // IDetector detector = mySURF.createDetector(); // interest_points = detector.generateInterestPoints(); // System.out.println("Interest points generated"); // IDescriptor descriptor = mySURF.createDescriptor(interest_points); // descriptor.generateAllDescriptors(); System.out.println("Performing SURF analysis on left image..."); Surf leftSurf = new Surf(leftImage, balanceValue, threshold, octaves); // if (leftSurf.getNumberOfPoints() > 500000) { // System.out.println("Number of points exceeds limit, reset threshold: " + leftSurf.getNumberOfPoints()); // return; // } if (leftSurf.getNumberOfPoints() == 0) { System.out .println("Number of points equals zero, reset threshold: " + leftSurf.getNumberOfPoints()); return; } System.out.println("Performing SURF analysis on right image..."); Surf rightSurf = new Surf(rightImage, balanceValue, threshold, octaves); if (rightSurf.getNumberOfPoints() == 0) { System.out .println("Number of points equals zero, reset threshold: " + leftSurf.getNumberOfPoints()); return; } System.out.println("Matching points of interest..."); Map<SURFInterestPoint, SURFInterestPoint> matchingPoints = leftSurf.getMatchingPoints(rightSurf, matchingThreshold, false); int numTiePoints = matchingPoints.size(); if (numTiePoints < 3) { System.err.println( "The number of potential tie points is less than 3. Adjust your threshold parameters and retry."); return; } System.out.println(numTiePoints + " potential tie points located"); System.out.println("Trimming outlier tie points..."); boolean flag; do { flag = false; leftTiePointsList.clear(); rightTiePointsList.clear(); i = 0; for (SURFInterestPoint point : matchingPoints.keySet()) { x = point.getX(); y = point.getY(); SURFInterestPoint target = matchingPoints.get(point); x2 = target.getX(); y2 = target.getY(); leftTiePointsList.add(new XYPoint(x, y)); rightTiePointsList.add(new XYPoint(x2, y2)); i++; } PolynomialLeastSquares2DFitting overallFit = new PolynomialLeastSquares2DFitting(leftTiePointsList, rightTiePointsList, 1); double maxDist = 0; SURFInterestPoint mostInfluentialPoint = null; i = 0; for (SURFInterestPoint point : matchingPoints.keySet()) { leftTiePointsList.clear(); rightTiePointsList.clear(); for (SURFInterestPoint point2 : matchingPoints.keySet()) { if (point2 != point) { x = point2.getX(); y = point2.getY(); SURFInterestPoint target = matchingPoints.get(point2); x2 = target.getX(); y2 = target.getY(); leftTiePointsList.add(new XYPoint(x, y)); rightTiePointsList.add(new XYPoint(x2, y2)); } } PolynomialLeastSquares2DFitting newFit = new PolynomialLeastSquares2DFitting(leftTiePointsList, rightTiePointsList, 1); x = point.getX(); y = point.getY(); XYPoint pt1 = overallFit.getForwardCoordinates(x, y); XYPoint pt2 = newFit.getForwardCoordinates(x, y); double dist = pt1.getSquareDistance(pt2); if (dist > maxDist) { maxDist = dist; mostInfluentialPoint = point; } } if (maxDist > 10 && mostInfluentialPoint != null) { matchingPoints.remove(mostInfluentialPoint); flag = true; } System.out.println(maxDist); } while (flag); int numPoints = matchingPoints.size(); // create homogeneous points matrices double[][] leftPoints = new double[3][numPoints]; double[][] rightPoints = new double[3][numPoints]; i = 0; for (SURFInterestPoint point : matchingPoints.keySet()) { leftPoints[0][i] = point.getX(); leftPoints[1][i] = point.getY(); leftPoints[2][i] = 1; SURFInterestPoint target = matchingPoints.get(point); rightPoints[0][i] = target.getX(); rightPoints[1][i] = target.getY(); rightPoints[2][i] = 1; i++; } double[][] normalizedLeftPoints = Normalize2DHomogeneousPoints.normalize(leftPoints); RealMatrix Tl = MatrixUtils.createRealMatrix(Normalize2DHomogeneousPoints.T); double[][] normalizedRightPoints = Normalize2DHomogeneousPoints.normalize(rightPoints); RealMatrix Tr = MatrixUtils.createRealMatrix(Normalize2DHomogeneousPoints.T); RealMatrix pnts1 = MatrixUtils.createRealMatrix(normalizedLeftPoints); RealMatrix pnts2 = MatrixUtils.createRealMatrix(normalizedRightPoints); RealMatrix A = MatrixUtils.createRealMatrix(buildA(normalizedLeftPoints, normalizedRightPoints)); //RealMatrix ata = A.transpose().multiply(A); SingularValueDecomposition svd = new SingularValueDecomposition(A); RealMatrix V = svd.getV(); RealVector V_smallestSingularValue = V.getColumnVector(8); RealMatrix F = MatrixUtils.createRealMatrix(3, 3); for (i = 0; i < 9; i++) { F.setEntry(i / 3, i % 3, V_smallestSingularValue.getEntry(i)); } for (i = 0; i < V.getRowDimension(); i++) { System.out.println(V.getRowVector(i).toString()); } SingularValueDecomposition svd2 = new SingularValueDecomposition(F); RealMatrix U = svd2.getU(); RealMatrix S = svd2.getS(); V = svd2.getV(); RealMatrix m = MatrixUtils.createRealMatrix( new double[][] { { S.getEntry(1, 1), 0, 0 }, { 0, S.getEntry(2, 2), 0 }, { 0, 0, 0 } }); F = U.multiply(m).multiply(V).transpose(); // Denormalise F = Tr.transpose().multiply(F).multiply(Tl); for (i = 0; i < F.getRowDimension(); i++) { System.out.println(F.getRowVector(i).toString()); } svd2 = new SingularValueDecomposition(F); //[U,D,V] = svd(F,0); RealMatrix e1 = svd2.getV().getColumnMatrix(2); //hnormalise(svd2.getV(:,3)); RealMatrix e2 = svd2.getU().getColumnMatrix(2); //e2 = hnormalise(U(:,3)); e1.setEntry(0, 0, (e1.getEntry(0, 0) / e1.getEntry(2, 0))); e1.setEntry(1, 0, (e1.getEntry(1, 0) / e1.getEntry(2, 0))); e1.setEntry(2, 0, 1); e2.setEntry(0, 0, (e2.getEntry(0, 0) / e2.getEntry(2, 0))); e2.setEntry(1, 0, (e2.getEntry(1, 0) / e2.getEntry(2, 0))); e2.setEntry(2, 0, 1); System.out.println(""); // boolean[] removeTiePoint = new boolean[numTiePoints]; // double[] residuals = null; // double[] residualsOrientation = null; // boolean flag; // do { // // perform the initial tie point transformation // leftTiePointsList.clear(); // rightTiePointsList.clear(); // int numPointsIncluded = 0; // i = 0; // for (SURFInterestPoint point : matchingPoints.keySet()) { // if (i < numTiePoints && !removeTiePoint[i]) { // x = point.getX(); // y = point.getY(); // // SURFInterestPoint target = matchingPoints.get(point); // x2 = target.getX(); // y2 = target.getY(); // // leftTiePointsList.add(new XYPoint(x, y)); // rightTiePointsList.add(new XYPoint(x2, y2)); // // numPointsIncluded++; // } // i++; // } // // PolynomialLeastSquares2DFitting plsFit = new PolynomialLeastSquares2DFitting( // leftTiePointsList, rightTiePointsList, 1); // // double rmse = plsFit.getOverallRMSE(); // System.out.println("RMSE: " + rmse + " with " + numPointsIncluded + " points included."); // // flag = false; // residuals = plsFit.getResidualsXY(); // residualsOrientation = plsFit.getResidualsOrientation(); // if (rmse > maxAllowableRMSE) { // i = 0; // for (k = 0; k < numTiePoints; k++) { // if (!removeTiePoint[k]) { // if (residuals[i] > 3 * rmse) { // removeTiePoint[k] = true; // flag = true; // } // i++; // } // } // } // } while (flag); // // i = 0; // for (k = 0; k < numTiePoints; k++) { // if (!removeTiePoint[k]) { // i++; // } // } System.out.println(numPoints + " tie points remain."); System.out.println("Outputing tie point files..."); ShapeFile leftOutput = new ShapeFile(leftOutputHeader, ShapeType.POINT, fields); ShapeFile rightOutput = new ShapeFile(rightOutputHeader, ShapeType.POINT, fields); i = 0; k = 0; for (SURFInterestPoint point : matchingPoints.keySet()) { // if (i < numTiePoints && !removeTiePoint[i]) { x = leftImage.getXCoordinateFromColumn((int) point.getX()); y = leftImage.getYCoordinateFromRow((int) point.getY()); SURFInterestPoint target = matchingPoints.get(point); x2 = rightImage.getXCoordinateFromColumn((int) target.getX()); y2 = rightImage.getYCoordinateFromRow((int) target.getY()); outputPoint = new whitebox.geospatialfiles.shapefile.Point(x, y); rowData = new Object[5]; rowData[0] = new Double(k + 1); rowData[1] = new Double(point.getOrientation()); rowData[2] = new Double(point.getScale()); rowData[3] = new Double(point.getLaplacian()); rowData[4] = new Double(0); //residuals[k]); leftOutput.addRecord(outputPoint, rowData); outputPoint = new whitebox.geospatialfiles.shapefile.Point(x2, y2); rowData = new Object[5]; rowData[0] = new Double(k + 1); rowData[1] = new Double(target.getOrientation()); rowData[2] = new Double(target.getScale()); rowData[3] = new Double(target.getLaplacian()); rowData[4] = new Double(0); //residuals[k]); rightOutput.addRecord(outputPoint, rowData); k++; // } i++; } leftOutput.write(); rightOutput.write(); System.out.println("Done!"); } catch (Exception e) { e.printStackTrace(); } }
From source file:plugins.TrendSurface.java
public double calculateEquation(double[] X, double[] Y, double[] Z) { try {/* w ww.j a v a 2s . c o m*/ int m, i, j, k; int n = Z.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } // 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(X[i], j) * Math.pow(Y[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(Z, false); RealVector solution = solver.solve(constants); regressCoefficents = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { regressCoefficents[a] = solution.getEntry(a); } double[] residuals = new double[n]; double SSresid = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * regressCoefficents[j]; } residuals[i] = Z[i] - yHat; SSresid += residuals[i] * residuals[i]; } double sum = 0; double SS = 0; for (i = 0; i < n; i++) { SS += Z[i] * Z[i]; sum += Z[i]; } double variance = (SS - (sum * sum) / n) / n; double SStotal = (n - 1) * variance; double rsq = 1 - SSresid / SStotal; return rsq; } catch (DimensionMismatchException | NoDataException | NullArgumentException | OutOfRangeException e) { showFeedback("Error in TrendSurface.calculateEquation: " + e.toString()); return -1; } }
From source file:plugins.TrendSurfaceVectorPoints.java
public double calculateEquation(double[] X, double[] Y, double[] Z) { try {/*from www .j ava2 s . c o m*/ int m, i, j, k; int n = Z.length; // How many coefficients are there? numCoefficients = 0; for (j = 0; j <= polyOrder; j++) { for (k = 0; k <= (polyOrder - j); k++) { numCoefficients++; } } // 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(X[i], j) * Math.pow(Y[i], k); m++; } } } RealMatrix coefficients = new Array2DRowRealMatrix(forwardCoefficientMatrix, false); //DecompositionSolver solver = new SingularValueDecomposition(coefficients).getSolver(); DecompositionSolver solver = new QRDecomposition(coefficients).getSolver(); // do the z-coordinate RealVector constants = new ArrayRealVector(Z, false); RealVector solution = solver.solve(constants); regressCoefficents = new double[numCoefficients]; for (int a = 0; a < numCoefficients; a++) { regressCoefficents[a] = solution.getEntry(a); } double[] residuals = new double[n]; double SSresid = 0; for (i = 0; i < n; i++) { double yHat = 0.0; for (j = 0; j < numCoefficients; j++) { yHat += forwardCoefficientMatrix[i][j] * regressCoefficents[j]; } residuals[i] = Z[i] - yHat; SSresid += residuals[i] * residuals[i]; } double sum = 0; double SS = 0; for (i = 0; i < n; i++) { SS += Z[i] * Z[i]; sum += Z[i]; } double variance = (SS - (sum * sum) / n) / n; double SStotal = (n - 1) * variance; double rsq = 1 - SSresid / SStotal; return rsq; } catch (DimensionMismatchException | NoDataException | NullArgumentException | OutOfRangeException e) { showFeedback("Error in TrendSurface.calculateEquation: " + e.toString()); return -1; } }