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 plugins; import jopensurf.*; import static java.lang.Math.*; import java.util.Arrays; import java.util.ArrayList; import java.util.List; import java.util.Iterator; import java.util.Map; import org.apache.commons.math3.linear.*; import org.apache.commons.math3.linear.SingularValueDecomposition; import whitebox.geospatialfiles.WhiteboxRaster; import whitebox.geospatialfiles.ShapeFile; import whitebox.geospatialfiles.shapefile.PointsList; import whitebox.geospatialfiles.shapefile.PolyLine; import whitebox.geospatialfiles.shapefile.ShapeFileRecord; import whitebox.geospatialfiles.shapefile.ShapeType; import whitebox.geospatialfiles.shapefile.attributes.DBFField; import whitebox.structures.XYPoint; import whitebox.stats.PolynomialLeastSquares2DFitting; import whitebox.structures.KdTree; import whitebox.structures.RowPriorityGridCell; import photogrammetry.Normalize2DHomogeneousPoints; /** * * @author johnlindsay */ public class SURFPixelMatching { public static void main(String[] args) { SURFPixelMatching surf = new SURFPixelMatching(); surf.run(); } private void run() { try { // 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(); } } private double[][] buildA(double[][] points1, double[][] points2) { int numPoints = points1[0].length; double[][] result = new double[numPoints][9]; for (int i = 0; i < numPoints; i++) { double x1i = points1[0][i]; double x2i = points2[0][i]; double y1i = points1[1][i]; double y2i = points2[1][i]; double[] row = new double[] { x1i * x2i, y1i * x2i, x2i, x1i * y2i, y1i * y2i, y2i, x1i, y1i, 1 }; System.arraycopy(row, 0, result[i], 0, 9); } return result; } }