/* * 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 static java.lang.Math.*; import java.util.Arrays; import java.util.ArrayList; import java.util.List; import whitebox.geospatialfiles.WhiteboxRaster; import whitebox.geospatialfiles.ShapeFile; 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; /** * Can't find * @author Dr. John Lindsay email: jlindsay@uoguelph.ca */ public class LocateConjugatePrincipalPoint { private WhiteboxRaster referenceImage = null; private WhiteboxRaster transformedImage = null; private XYAndDirection[][] offsets; private static double epsilon = 1.0; private int maxNeighbourhoodSize = 400; private int[] numCellsInAnnulus; private double referenceNoData; private double transformedNoData; public static void main(String[] args) { LocateConjugatePrincipalPoint lcpp = new LocateConjugatePrincipalPoint(); lcpp.run(); } /** * Used to execute this plugin tool. */ private void run() { long startTime = System.nanoTime(); ShapeFile output = null; ShapeFile output2 = null; int progress, oldProgress; boolean conductFineSearch = false; int j; int refNeighbourhoodStart = 40; int refNeighbourhoodStep = 20; maxNeighbourhoodSize = 500; epsilon = 1.2; int polyOrder = 2; StringBuilder str; KdTree<Double> controlPointTree = new KdTree.SqrEuclid<>(2, new Integer(2000)); try { //String referenceFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-82_Blue_clipped_int.dep"; //String ppFile1 = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-82 principal point.shp"; String ppFile1 = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/test point3.shp"; // String transformedFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/Guelph_A19409-83_Blue_clipped_int.dep"; //String outputFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/pp mapped.shp"; //String outputFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/test point2 mapped.shp"; // String outputFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/tmp4.shp"; // String referenceTiePoints = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/82 tie points.shp"; // String transformedTiePoints = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/83 tie points.shp"; String referenceFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus_C6430-74072-L9_254_Blue_clipped.dep"; String transformedFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/GuelphCampus_C6430-74072-L9_253_Blue_clipped.dep"; String referenceTiePoints = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/campus 254 tie points.shp"; String transformedTiePoints = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/campus 253 tie points.shp"; String outputFile = "/Users/johnlindsay/Documents/Teaching/GEOG2420/airphotos/tmp6.shp"; DBFField[] fields = new DBFField[1]; fields[0] = new DBFField(); fields[0].setName("r1"); fields[0].setDataType(DBFField.DBFDataType.NUMERIC); fields[0].setDecimalCount(4); fields[0].setFieldLength(10); // fields[1] = new DBFField(); // fields[1].setName("r2"); // fields[1].setDataType(DBFField.DBFDataType.NUMERIC); // fields[1].setDecimalCount(4); // fields[1].setFieldLength(10); output = new ShapeFile(outputFile, ShapeType.POINT, fields); fields = new DBFField[1]; fields[0] = new DBFField(); fields[0].setName("r1"); fields[0].setDataType(DBFField.DBFDataType.NUMERIC); fields[0].setDecimalCount(4); fields[0].setFieldLength(10); output2 = new ShapeFile(outputFile.replace(".shp", "_2.shp"), ShapeType.POINT, fields); referenceImage = new WhiteboxRaster(referenceFile, "r"); referenceImage.setForceAllDataInMemory(true); int rows1 = referenceImage.getNumberRows(); int cols1 = referenceImage.getNumberColumns(); referenceNoData = referenceImage.getNoDataValue(); // WhiteboxRaster outputRaster = new WhiteboxRaster(outputFile.replace(".shp", ".dep"), "rw", referenceFile, WhiteboxRaster.DataType.FLOAT, 0); // epsilon = 5; // double[] data; // oldProgress = -1; // progress = 0; // for (int r = 0; r < rows1; r++) { // data = referenceImage.getRowValues(r); // double[][] data2 = new double[cols1][2]; // for (int c = 0; c < cols1; c++) { // data2[c][0] = c; // data2[c][1] = data[c]; // } // double[][] filteredData = douglasPeuckerFilter(data2, 0, cols1 - 1); // for (i = 1; i < filteredData.length; i++) { // j = (int)filteredData[i - 1][0]; // int k = (int)filteredData[i][0]; // double s = k - j; // double startZ = filteredData[i - 1][1]; // double endZ = filteredData[i][1]; // for (int m = j; m <= k; m++) { // double outZ = startZ + ((m - j) / s) * (endZ - startZ); // outputRaster.setValue(r, m, outZ); // } // //outputRaster.setValue(r, (int)filteredData[i][0], 1); // } // progress = (int)(100f * r / (rows1 - 1)); // if (progress > oldProgress) { // System.out.println(progress + "%"); // oldProgress = progress; // } // } // outputRaster.close(); transformedImage = new WhiteboxRaster(transformedFile, "r"); transformedImage.setForceAllDataInMemory(true); transformedNoData = transformedImage.getNoDataValue(); ShapeFile pp1 = new ShapeFile(ppFile1); if (pp1.getShapeType().getBaseType() != ShapeType.POINT && pp1.getShapeType().getBaseType() != ShapeType.MULTIPOINT) { //showFeedback("The input shapefile must be of a 'POINT' or 'MULTIPOINT' data type."); return; } ShapeFileRecord record = pp1.getRecord(0); double[][] point; // = record.getGeometry().getPoints(); ShapeFile refTiePoints = new ShapeFile(referenceTiePoints); if (refTiePoints.getShapeType().getBaseType() != ShapeType.POINT) { //showFeedback("The input shapefile must be of a 'POINT' or 'MULTIPOINT' data type."); return; } ShapeFile transTiePoints = new ShapeFile(transformedTiePoints); if (transTiePoints.getShapeType().getBaseType() != ShapeType.POINT) { //showFeedback("The input shapefile must be of a 'POINT' or 'MULTIPOINT' data type."); return; } int numTiePoints = refTiePoints.getNumberOfRecords(); if (transTiePoints.getNumberOfRecords() != numTiePoints) { return; } calculateOffsets(); conductFineSearch = true; ArrayList<XYPoint> tiePointsRef = new ArrayList<>(); ArrayList<XYPoint> tiePointsTransform = new ArrayList<>(); for (int r = 0; r < refTiePoints.getNumberOfRecords(); r++) { double[][] refPoint = refTiePoints.getRecord(r).getGeometry().getPoints(); int refCol = referenceImage.getColumnFromXCoordinate(refPoint[0][0]); int refRow = referenceImage.getRowFromYCoordinate(refPoint[0][1]); point = transTiePoints.getRecord(r).getGeometry().getPoints(); int transCol = transformedImage.getColumnFromXCoordinate(point[0][0]); int transRow = transformedImage.getRowFromYCoordinate(point[0][1]); RowPriorityGridCell gc = findPixelMatch(refCol, refRow, transCol, transRow, conductFineSearch, refNeighbourhoodStart, refNeighbourhoodStep, 30, 1.0); System.out.println("Control Point " + (r + 1) + ": " + gc.z); int matchedCol = gc.col; int matchedRow = gc.row; double matchedCorrelation = gc.z; if (matchedCorrelation >= 0.95) { double x2 = transformedImage.getXCoordinateFromColumn(matchedCol); double y2 = transformedImage.getYCoordinateFromRow(matchedRow); whitebox.geospatialfiles.shapefile.Point PP = new whitebox.geospatialfiles.shapefile.Point(x2, y2); Object[] rowData = new Object[1]; rowData[0] = new Double(matchedCorrelation); output.addRecord(PP, rowData); PP = new whitebox.geospatialfiles.shapefile.Point(refPoint[0][0], refPoint[0][1]); rowData = new Object[2]; rowData[0] = new Double(matchedCorrelation); rowData[1] = new Double(0.0); output2.addRecord(PP, rowData); tiePointsRef.add(new XYPoint(refPoint[0][0], refPoint[0][1])); tiePointsTransform.add(new XYPoint(x2, y2)); } else { System.out.println("No suitable match could be located."); } } conductFineSearch = false; List<KdTree.Entry<Double>> results; int newPolyOrder = polyOrder; if (newPolyOrder == 4 && tiePointsRef.size() < 15) { newPolyOrder--; } if (newPolyOrder == 3 && tiePointsRef.size() < 10) { newPolyOrder--; } if (newPolyOrder == 2 && tiePointsRef.size() < 6) { newPolyOrder--; } numTiePoints = 0; for (XYPoint tie : tiePointsRef) { double[] entry = {tie.x, tie.y}; controlPointTree.addPoint(entry, (double) numTiePoints); numTiePoints++; } PolynomialLeastSquares2DFitting pls = new PolynomialLeastSquares2DFitting( tiePointsRef, tiePointsTransform, newPolyOrder); double rmse = pls.getOverallRMSE(); System.out.println("\nRMSE: " + rmse); double north = transformedImage.getNorth(); double south = transformedImage.getSouth(); double east = transformedImage.getEast(); double west = transformedImage.getWest(); int totalPointsSearched = 0; int interval = 1000; double intervalSteps = 1.5; int loopNum = 1; do { System.out.println("Interval: " + interval); oldProgress = -1; for (int r = 0; r < rows1; r += interval) { for (int c = 0; c < cols1; c += interval) { if (referenceImage.getValue(r, c) != referenceNoData) { double refXCoord = referenceImage.getXCoordinateFromColumn(c); double refYCoord = referenceImage.getYCoordinateFromRow(r); double[] entry = {refXCoord, refYCoord}; int numNearestNeighbours = 15; if (numTiePoints < 15) { numNearestNeighbours = numTiePoints; } results = controlPointTree.nearestNeighbor(entry, numNearestNeighbours, true); j = results.size(); double[] X1 = new double[j]; double[] Y1 = new double[j]; double[] X2 = new double[j]; double[] Y2 = new double[j]; for (int k = 0; k < j; k++) { double val = results.get(k).value; X1[k] = tiePointsRef.get((int) val).x; Y1[k] = tiePointsRef.get((int) val).y; X2[k] = tiePointsTransform.get((int) val).x; Y2[k] = tiePointsTransform.get((int) val).y; } int count = 0; double scaleFactor = 0; for (int k = 0; k < j; k++) { double x1Ref = X1[k]; double y1Ref = Y1[k]; double x1Tr = X2[k]; double y1Tr = Y2[k]; for (int m = k + 1; m < j; m++) { double x2Ref = X1[m]; double y2Ref = Y1[m]; double x2Tr = X2[m]; double y2Tr = Y2[m]; double dist1 = sqrt((x2Ref - x1Ref) * (x2Ref - x1Ref) + (y2Ref - y1Ref) * (y2Ref - y1Ref)); double dist2 = sqrt((x2Tr - x1Tr) * (x2Tr - x1Tr) + (y2Tr - y1Tr) * (y2Tr - y1Tr)); if (dist1 > 0) { scaleFactor += dist2 / dist1; count++; } } } scaleFactor = scaleFactor / count; pls = new PolynomialLeastSquares2DFitting(X1, Y1, X2, Y2, 1); rmse = pls.getOverallRMSE(); XYPoint transCoords = pls.getForwardCoordinates( refXCoord, refYCoord); if (transCoords.x <= east && transCoords.x >= west && transCoords.y >= south && transCoords.y <= north) { totalPointsSearched++; int transCol = transformedImage.getColumnFromXCoordinate(transCoords.x); int transRow = transformedImage.getRowFromYCoordinate(transCoords.y); int searchWindowRadius = (int) rmse * 2; if (searchWindowRadius < 80) { searchWindowRadius = 80; } RowPriorityGridCell gc = findPixelMatch(c, r, transCol, transRow, conductFineSearch, refNeighbourhoodStart, refNeighbourhoodStep, searchWindowRadius, scaleFactor); int matchedCol = gc.col; int matchedRow = gc.row; double matchedCorrelation = gc.z; if (matchedCorrelation >= 0.95) { double x2 = transformedImage.getXCoordinateFromColumn(matchedCol); double y2 = transformedImage.getYCoordinateFromRow(matchedRow); whitebox.geospatialfiles.shapefile.Point PP = new whitebox.geospatialfiles.shapefile.Point(x2, y2); Object[] rowData = new Object[1]; rowData[0] = new Double(matchedCorrelation); output.addRecord(PP, rowData); PP = new whitebox.geospatialfiles.shapefile.Point(refXCoord, refYCoord); rowData = new Object[1]; rowData[0] = new Double(matchedCorrelation); output2.addRecord(PP, rowData); tiePointsRef.add(new XYPoint(refXCoord, refYCoord)); tiePointsTransform.add(new XYPoint(x2, y2)); entry = new double[]{refXCoord, refYCoord}; controlPointTree.addPoint(entry, (double) numTiePoints); numTiePoints++; newPolyOrder = polyOrder; if (newPolyOrder == 4 && tiePointsRef.size() < 15) { newPolyOrder--; } if (newPolyOrder == 3 && tiePointsRef.size() < 10) { newPolyOrder--; } if (newPolyOrder == 2 && tiePointsRef.size() < 6) { newPolyOrder--; } //System.out.println("Num. tie points: " + tiePointsRef.size() + " of " + totalPointsSearched + " (" + (100f * tiePointsRef.size() / totalPointsSearched) + "%)" + " , RMSE: " + rmse); } } } } progress = (int) ((100.0 * r) / rows1); if (progress > oldProgress) { System.out.println("Loop " + loopNum + " " + progress + "%" + ", Num. tie points: " + tiePointsRef.size() + " of " + totalPointsSearched + " (" + (100f * tiePointsRef.size() / totalPointsSearched) + "%)"); oldProgress = progress; } } loopNum++; interval = (int) (interval / intervalSteps); } while (interval >= 200); referenceImage.close(); transformedImage.close(); output.write(); output2.write(); System.out.println("\nOperation complete!"); long endTime = System.nanoTime(); double duration = (endTime - startTime); int secs = (int) (duration / 1000000000); // int days = secs / 86400; // secs = secs - 86400 * days; int hours = secs / 3600; secs = secs - 3600 * hours; int minutes = secs / 60; secs = secs - minutes * 60; int seconds = secs; str = new StringBuilder(); str.append("Duration: "); // if (days > 0) { // str.append(days).append(" days, "); // } if (hours > 0) { str.append(hours).append(" hours, "); } if (minutes > 0) { str.append(minutes).append(" minutes, "); } if (seconds > 0) { str.append(seconds).append(" seconds, "); } System.out.println(str.toString()); } catch (Exception e) { if (output != null && output2 != null) { try { output.write(); output2.write(); } catch (Exception e2) { } } e.printStackTrace(); } } private RowPriorityGridCell findPixelMatch(int refCol, int refRow, int transCol, int transRow, boolean conductFineSearch, int neighbourhoodStart, int neighbourhoodStep, int searchWindowRadius, double scaleFactor) { int i, a, row, col, n; double z, M, Q, w1, w2, w3; boolean[] coarsereferenceRings; boolean[] annulusVisited = new boolean[maxNeighbourhoodSize + 1]; int referenceRadiusN = 0; int referenceRadius = neighbourhoodStart; double[] referenceMeans = new double[maxNeighbourhoodSize + 1]; double[] referenceVariances = new double[maxNeighbourhoodSize + 1]; double[] referenceLumped = new double[maxNeighbourhoodSize + 1]; double[][] filterData1 = new double[maxNeighbourhoodSize + 1][2]; double[][] filterData2 = new double[maxNeighbourhoodSize + 1][2]; double[][] filterData3 = new double[maxNeighbourhoodSize + 1][2]; M = 0; Q = 0; boolean flag = true; do { for (i = 1; i <= referenceRadius; i++) { if (!annulusVisited[i]) { double total = 0; n = 0; M = 0; Q = 0; double previousZ = referenceNoData; double totalDiff = 0; for (a = 0; a < numCellsInAnnulus[i]; a++) { row = refRow + offsets[i][a].y; col = refCol + offsets[i][a].x; z = referenceImage.getValue(row, col); if (z != referenceNoData) { total += z; n++; if (a > 0) { M = M + (z - M) / (a + 1); Q = Q + (a * (z - M) * (z - M)) / (a + 1); } else { M = z; Q = 0; } if (previousZ != referenceNoData) { totalDiff += abs(z - previousZ); } } previousZ = z; } if (n > 1) { referenceMeans[i] = total / n; referenceVariances[i] = sqrt(Q / (n - 1)); referenceLumped[i] = totalDiff / (n - 1); } else { referenceMeans[i] = 0; referenceVariances[i] = 0; referenceLumped[i] = 0; } filterData1[i][0] = i; filterData1[i][1] = referenceMeans[i]; filterData2[i][0] = i; filterData2[i][1] = referenceVariances[i]; filterData3[i][0] = i; filterData3[i][1] = referenceLumped[i]; annulusVisited[i] = true; } } double[][] newData1 = douglasPeuckerFilter(filterData1, 1, referenceRadius); double[][] newData2 = douglasPeuckerFilter(filterData2, 1, referenceRadius); double oldEpsilon = epsilon; // epsilon = 0.5; // double[][] newData3 = douglasPeuckerFilter(filterData3, 1, referenceRadius); // epsilon = oldEpsilon; referenceRadiusN = 0; coarsereferenceRings = new boolean[referenceRadius + 1]; for (i = 0; i < newData1.length; i++) { coarsereferenceRings[(int) newData1[i][0]] = true; referenceRadiusN++; } for (i = 0; i < newData2.length; i++) { if (coarsereferenceRings[(int) newData2[i][0]] == false) { coarsereferenceRings[(int) newData2[i][0]] = true; referenceRadiusN++; } } w1 = (double) newData1.length / (newData1.length + newData2.length); w2 = (double) newData2.length / (newData1.length + newData2.length); // w1 = (double) newData1.length / (newData1.length + newData2.length + newData3.length); // w2 = (double) newData2.length / (newData1.length + newData2.length + newData3.length); // w3 = (double) newData3.length / (newData1.length + newData2.length + newData3.length); if (newData1.length > 8 && newData2.length > 8 && referenceRadiusN > 12) { // if (newData1.length > 8 && newData2.length > 8 && newData3.length > 8 && referenceRadiusN > 12) { // there must be enough information on both the means data // and the variance data. flag = false; } else { referenceRadius += neighbourhoodStep; if (referenceRadius > maxNeighbourhoodSize) { referenceRadius = maxNeighbourhoodSize; flag = false; } } } while (flag); double referenceMean = 0; double referenceVariance = 0; double referenceLump = 0; double referenceMeanDetailed = 0; double referenceVarianceDetailed = 0; double referenceLumpDetailed = 0; for (a = 1; a < referenceRadius + 1; a++) { referenceMeanDetailed += referenceMeans[a]; referenceVarianceDetailed += referenceVariances[a]; referenceLumpDetailed += referenceLumped[a]; if (coarsereferenceRings[a]) { referenceMean += referenceMeans[a]; referenceVariance += referenceVariances[a]; referenceLump += referenceLumped[a]; } } referenceMean = referenceMean / referenceRadiusN; referenceVariance = referenceVariance / referenceRadiusN; referenceLump = referenceLump / referenceRadiusN; referenceMeanDetailed = referenceMeanDetailed / referenceRadius; referenceVarianceDetailed = referenceVarianceDetailed / referenceRadius; referenceLumpDetailed = referenceLumpDetailed / referenceRadius; double[] referenceMeanDeviates = new double[referenceRadius + 1]; double[] referenceVarianceDeviates = new double[referenceRadius + 1]; double[] referenceLumpDeviates = new double[referenceRadius + 1]; double[] referenceMeanDeviatesDetailed = new double[referenceRadius + 1]; double[] referenceVarianceDeviatesDetailed = new double[referenceRadius + 1]; double[] referenceLumpDeviatesDetailed = new double[referenceRadius + 1]; double sqrDev1 = 0; double sqrDev2 = 0; double sqrDev3 = 0; double sqrDev1Detailed = 0; double sqrDev2Detailed = 0; double sqrDev3Detailed = 0; for (a = 1; a < referenceRadius + 1; a++) { referenceMeanDeviatesDetailed[a] = referenceMeans[a] - referenceMeanDetailed; referenceVarianceDeviatesDetailed[a] = referenceVariances[a] - referenceVarianceDetailed; referenceLumpDeviatesDetailed[a] = referenceLumped[a] - referenceLumpDetailed; sqrDev1Detailed += (referenceMeans[a] - referenceMeanDetailed) * (referenceMeans[a] - referenceMeanDetailed); sqrDev2Detailed += (referenceVariances[a] - referenceVarianceDetailed) * (referenceVariances[a] - referenceVarianceDetailed); sqrDev3Detailed += (referenceLumped[a] - referenceLumpDetailed) * (referenceLumped[a] - referenceLumpDetailed); referenceMeanDeviates[a] = referenceMeans[a] - referenceMean; referenceVarianceDeviates[a] = referenceVariances[a] - referenceVariance; referenceLumpDeviates[a] = referenceLumped[a] - referenceLump; if (coarsereferenceRings[a]) { sqrDev1 += (referenceMeans[a] - referenceMean) * (referenceMeans[a] - referenceMean); sqrDev2 += (referenceVariances[a] - referenceVariance) * (referenceVariances[a] - referenceVariance); sqrDev3 += (referenceLumped[a] - referenceLump) * (referenceLumped[a] - referenceLump); } } double maxCorrelations = 0; int maxCorrelationRow = -1, maxCorrelationCol = -1; for (int row2 = transRow - searchWindowRadius; row2 <= transRow + searchWindowRadius; row2++) { for (int col2 = transCol - searchWindowRadius; col2 <= transCol + searchWindowRadius; col2++) { // find the means and variances of each annulus double[] means = new double[referenceRadius + 1]; double[] variances = new double[referenceRadius + 1]; double[] lumps = new double[referenceRadius + 1]; for (i = 1; i <= referenceRadius; i++) { if (coarsereferenceRings[i]) { double total = 0; n = 0; M = 0; Q = 0; double previousZ = transformedNoData; double totalDiff = 0; int scaled_i = (int) round(i * scaleFactor); for (a = 0; a < numCellsInAnnulus[scaled_i]; a++) { row = row2 + offsets[scaled_i][a].y; col = col2 + offsets[scaled_i][a].x; z = transformedImage.getValue(row, col); if (z != transformedNoData) { total += z; n++; if (a > 0) { M = M + (z - M) / (a + 1); Q = Q + (a * (z - M) * (z - M)) / (a + 1); } else { M = z; Q = 0; } if (previousZ != transformedNoData) { totalDiff += abs(z - previousZ); } } previousZ = z; } if (n > 1) { means[i] = total / n; variances[i] = sqrt(Q / (n - 1)); lumps[i] = totalDiff / (n - 1); } else { means[i] = 0; variances[i] = 0; lumps[i] = 0; } } } double sampleMean = 0; double sampleVariance = 0; double sampleLump = 0; for (a = 1; a < referenceRadius + 1; a++) { if (coarsereferenceRings[a]) { sampleMean += means[a]; sampleVariance += variances[a]; sampleLump += lumps[a]; } } sampleMean = sampleMean / referenceRadiusN; sampleVariance = sampleVariance / referenceRadiusN; sampleLump = sampleLump / referenceRadiusN; double sampleSqrdDev1 = 0; double sampleSqrdDev2 = 0; double sampleSqrdDev3 = 0; // correlate the reference and sample means and variances double cov1 = 0; double cov2 = 0; double cov3 = 0; for (a = 1; a < referenceRadius + 1; a++) { if (coarsereferenceRings[a]) { cov1 += (means[a] - sampleMean) * referenceMeanDeviates[a]; cov2 += (variances[a] - sampleVariance) * referenceVarianceDeviates[a]; cov3 += (lumps[a] - sampleLump) * referenceLumpDeviates[a]; sampleSqrdDev1 += (means[a] - sampleMean) * (means[a] - sampleMean); sampleSqrdDev2 += (variances[a] - sampleVariance) * (variances[a] - sampleVariance); sampleSqrdDev3 += (lumps[a] - sampleLump) * (lumps[a] - sampleLump); } } double r1 = cov1 / (Math.sqrt(sqrDev1 * sampleSqrdDev1)); double r2 = cov2 / (Math.sqrt(sqrDev2 * sampleSqrdDev2)); // double r3 = cov3 / (Math.sqrt(sqrDev3 * sampleSqrdDev3)); if (!conductFineSearch) { if (r1 * w1 + r2 * w2 > maxCorrelations) { maxCorrelations = r1 * w1 + r2 * w2; maxCorrelationCol = col2; maxCorrelationRow = row2; } } else { if (r1 * w1 + r2 * w2 > 0.9) { // conduct a detailed correlation for (i = 1; i <= referenceRadius; i++) { if (!coarsereferenceRings[i]) { double total = 0; n = 0; double previousZ = transformedNoData; double totalDiff = 0; M = 0; Q = 0; int scaled_i = (int) round(i * scaleFactor); for (a = 0; a < numCellsInAnnulus[scaled_i]; a++) { row = row2 + offsets[scaled_i][a].y; col = col2 + offsets[scaled_i][a].x; z = transformedImage.getValue(row, col); if (z != transformedNoData) { total += z; n++; if (a > 0) { M = M + (z - M) / (a + 1); Q = Q + (a * (z - M) * (z - M)) / (a + 1); } else { M = z; Q = 0; } if (previousZ != transformedNoData) { totalDiff += abs(z - previousZ); } } previousZ = z; } if (n > 1) { means[i] = total / n; variances[i] = sqrt(Q / (n - 1)); lumps[i] = totalDiff / (n - 1); } else { means[i] = 0; variances[i] = 0; lumps[i] = 0; } } } sampleMean = 0; sampleVariance = 0; sampleLump = 0; for (a = 1; a < referenceRadius + 1; a++) { sampleMean += means[a]; sampleVariance += variances[a]; sampleLump += lumps[a]; } sampleMean = sampleMean / referenceRadius; sampleVariance = sampleVariance / referenceRadius; sampleLump = sampleLump / referenceRadius; // correlate the reference and sample means and variances cov1 = 0; cov2 = 0; cov3 = 0; sampleSqrdDev1 = 0; sampleSqrdDev2 = 0; sampleSqrdDev3 = 0; for (a = 1; a < referenceRadius + 1; a++) { cov1 += (means[a] - sampleMean) * referenceMeanDeviatesDetailed[a]; cov2 += (variances[a] - sampleVariance) * referenceVarianceDeviatesDetailed[a]; cov3 += (lumps[a] - sampleLump) * referenceLumpDeviatesDetailed[a]; sampleSqrdDev1 += (means[a] - sampleMean) * (means[a] - sampleMean); sampleSqrdDev2 += (variances[a] - sampleVariance) * (variances[a] - sampleVariance); sampleSqrdDev3 += (lumps[a] - sampleLump) * (lumps[a] - sampleLump); } r1 = cov1 / (Math.sqrt(sqrDev1Detailed * sampleSqrdDev1)); r2 = cov2 / (Math.sqrt(sqrDev2Detailed * sampleSqrdDev2)); // r3 = cov3 / (Math.sqrt(sqrDev3Detailed * sampleSqrdDev3)); if (r1 * w1 + r2 * w2 > maxCorrelations) { maxCorrelations = r1 * w1 + r2 * w2; maxCorrelationCol = col2; maxCorrelationRow = row2; } } } } } RowPriorityGridCell retCell = new RowPriorityGridCell(maxCorrelationRow, maxCorrelationCol, maxCorrelations); return retCell; } private void calculateOffsets() { // calculate the offsets int i, j, x, y, row, col, a, b; double dist; int maxOffset = maxNeighbourhoodSize * 2; //will accomodate a scalefactor of 2 numCellsInAnnulus = new int[maxOffset + 1]; x = maxOffset + 1; y = maxOffset + 1; for (row = 0; row <= maxOffset * 2 + 1; row++) { a = row - y; for (col = 0; col <= maxOffset * 2 + 1; col++) { b = col - x; dist = sqrt(a * a + b * b); i = (int) (round(dist)); if (i <= maxOffset) { numCellsInAnnulus[i]++; } } } offsets = new XYAndDirection[maxOffset + 1][]; for (i = 1; i <= maxOffset; i++) { offsets[i] = new XYAndDirection[numCellsInAnnulus[i]]; for (j = 0; j < numCellsInAnnulus[i]; j++) { offsets[i][j] = new XYAndDirection(); } } int[] currentNumInAnnulus = new int[maxOffset + 1]; for (row = 0; row <= maxOffset * 2 + 1; row++) { a = row - y; for (col = 0; col <= maxOffset * 2 + 1; col++) { b = col - x; dist = sqrt(a * a + b * b); i = (int) (round(dist)); if (i <= maxOffset && i > 0) { offsets[i][currentNumInAnnulus[i]].x = b; offsets[i][currentNumInAnnulus[i]].y = a; offsets[i][currentNumInAnnulus[i]].direction = atan2(-a, b); currentNumInAnnulus[i]++; } } } for (i = 1; i <= maxOffset; i++) { Arrays.sort(offsets[i]); } } private static double[][] douglasPeuckerFilter(double[][] points, int startIndex, int endIndex) { double dmax = 0; int index = 0; double a = endIndex - startIndex; double b = points[endIndex][1] - points[startIndex][1]; double c = -(b * startIndex - a * points[startIndex][1]); double norm = sqrt(a * a + b * b); for (int i = startIndex + 1; i < endIndex; i++) { double distance = abs(b * i - a * points[i][1] + c) / norm; if (distance > dmax) { index = i; dmax = distance; } } if (dmax >= epsilon) { double[][] recursiveResult1 = douglasPeuckerFilter(points, startIndex, index); double[][] recursiveResult2 = douglasPeuckerFilter(points, index, endIndex); double[][] result = new double[(recursiveResult1.length - 1) + recursiveResult2.length][2]; int result1Length = recursiveResult1.length; for (int i = 0; i < result1Length - 1; i++) { result[i][0] = recursiveResult1[i][0]; result[i][1] = recursiveResult1[i][1]; } for (int i = 0; i < recursiveResult2.length; i++) { result[result1Length + i - 1][0] = recursiveResult2[i][0]; result[result1Length + i - 1][1] = recursiveResult2[i][1]; } return result; } else { double[][] ret = { {points[startIndex][0], points[startIndex][1]}, {points[endIndex][0], points[endIndex][1]} }; return ret; } } class XYAndDirection implements Comparable<XYAndDirection> { private int x = Integer.MIN_VALUE; private int y = Integer.MIN_VALUE; private double direction = Double.NEGATIVE_INFINITY; @Override public int compareTo(XYAndDirection o) { final int BEFORE = -1; final int EQUAL = 0; final int AFTER = 1; if (this.direction > o.direction) { return BEFORE; } else if (this.direction < o.direction) { return AFTER; } if (this.x < o.x) { return BEFORE; } else if (this.x > o.x) { return AFTER; } if (this.y < o.y) { return BEFORE; } else if (this.y > o.y) { return AFTER; } return EQUAL; } } }