From 9ee3d2265f2eaa911c732bd92b6a95bc70320b5e Mon Sep 17 00:00:00 2001 From: Diego Date: Mon, 29 May 2023 17:33:45 +0200 Subject: [PATCH 1/3] First draft of structural similarity class --- .../OmsStructuralSimilarity.java | 409 ++++++++++++++++++ 1 file changed, 409 insertions(+) create mode 100644 hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java diff --git a/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java new file mode 100644 index 000000000..a157d72c5 --- /dev/null +++ b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java @@ -0,0 +1,409 @@ +package org.hortonmachine.hmachine.modules.statistics.kerneldensity; + +import oms3.annotations.Execute; +import org.geotools.coverage.grid.GridCoverage2D; +import org.hortonmachine.gears.libs.exceptions.ModelsIllegalargumentException; +import org.hortonmachine.gears.libs.exceptions.ModelsRuntimeException; +import org.hortonmachine.gears.libs.modules.HMModel; +import org.hortonmachine.gears.libs.modules.HMRaster; +import org.jaitools.media.jai.kernel.KernelFactory; + +import javax.media.jai.KernelJAI; +import java.io.IOException; +import java.util.stream.IntStream; + +public class OmsStructuralSimilarity extends HMModel { + + + public GridCoverage2D inMap1 = null; + + public GridCoverage2D inMap2 = null; + + public double pK1 = 0.01; + public double pK2 = 0.03; + + public double pRelevanceMean = 1; + public double pRelevanceVariance = 1; + public double pRelevancePattern = 1; + + public int pKernel = 3; + + public int pRadius = 10; + + public boolean doConstant = false; + + public GridCoverage2D outStructuralSimilarity = null; + + public GridCoverage2D outMeanSimilarity = null; + + public GridCoverage2D outVarianceSimilarity = null; + + public GridCoverage2D outPatternSimilarity = null; + + + private volatile boolean errorOccurred = false; + private volatile String errorMessage; + + + @Execute + public void process() throws Exception { + + GridCoverage2D meanMap1 = windowMean(inMap1); + GridCoverage2D meanMap2 = windowMean(inMap2); + + GridCoverage2D varianceMap1 = windowVariance(inMap1,meanMap1); + GridCoverage2D varianceMap2 = windowVariance(inMap2,meanMap2); + + GridCoverage2D covarianceMap = windowCovariance(inMap1,inMap2,meanMap1,meanMap2); + + try (HMRaster meanRaster1 = HMRaster.fromGridCoverage(meanMap1); + HMRaster meanRaster2 = HMRaster.fromGridCoverage(meanMap2); + HMRaster varianceRaster1 = HMRaster.fromGridCoverage(varianceMap1); + HMRaster varianceRaster2 = HMRaster.fromGridCoverage(varianceMap2); + HMRaster covarianceRaster = HMRaster.fromGridCoverage(covarianceMap) + ) { + int cols = meanRaster1.getCols(); + int rows = meanRaster1.getRows(); + + // TODO: replace 0.0 by the value range among the two comparable rasters. + double c1 = Math.pow(pK1 * 0.0,2); + double c2 = Math.pow(pK2 * 0.0,2); + double c3 = 0.5*c2; + + HMRaster outputSSRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + HMRaster outputMSRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + HMRaster outputVSRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + HMRaster outputPSRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + + pm.beginTask("Estimating structural similarity...", cols ); + + IntStream.range(0, rows).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = 0; c < cols; c++ ) { + + double meanValue1 = meanRaster1.getValue(c, r); + double meanValue2 = meanRaster2.getValue(c, r); + double varianceValue1 = varianceRaster1.getValue(c, r); + double varianceValue2 = varianceRaster2.getValue(c, r); + double covarianceValue = covarianceRaster.getValue(c, r); + + if (meanRaster1.isNovalue(meanValue1) || meanRaster2.isNovalue(meanValue2) || varianceRaster1.isNovalue(varianceValue1) || varianceRaster2.isNovalue(varianceValue2) || covarianceRaster.isNovalue(covarianceValue)) { + continue; + } + + double meanSimilarity = (2*meanValue1*meanValue2 + c1)/(meanValue1*meanValue1 + meanValue2*meanValue2 + c1 ); + double varianceSimilarity = (2*Math.sqrt(varianceValue1)*Math.sqrt(varianceValue2) +c2)/(varianceValue1 + varianceValue2 + c2); + double patternSimilarity = (covarianceValue + c3)/(Math.sqrt(varianceValue1)*Math.sqrt(varianceValue2) + c3); + double structuralSimilarity = Math.pow(meanSimilarity,pRelevanceMean)*Math.pow(varianceSimilarity,pRelevanceVariance)*Math.pow(patternSimilarity,pRelevancePattern); + + try { + outputSSRaster.setValue(c, r, structuralSimilarity); + outputMSRaster.setValue(c, r, meanSimilarity); + outputVSRaster.setValue(c, r, varianceSimilarity); + outputPSRaster.setValue(c, r, patternSimilarity); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } + } + pm.worked(1); + }); + + pm.done(); + + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); + } + + outStructuralSimilarity = outputSSRaster.buildCoverage(); + outMeanSimilarity = outputVSRaster.buildCoverage(); + outVarianceSimilarity = outputVSRaster.buildCoverage(); + outPatternSimilarity = outputPSRaster.buildCoverage(); + } + + } + + + private GridCoverage2D windowCovariance(GridCoverage2D inMap1, GridCoverage2D inMap2, GridCoverage2D meanMap1, GridCoverage2D meanMap2) throws Exception { + + checkNull(inMap1); + checkNull(inMap2); + checkNull(meanMap1); + checkNull(meanMap2); + + try (HMRaster inRaster1 = HMRaster.fromGridCoverage(inMap1); + HMRaster meanRaster1 = HMRaster.fromGridCoverage(meanMap1); + HMRaster inRaster2 = HMRaster.fromGridCoverage(inMap2); + HMRaster meanRaster2 = HMRaster.fromGridCoverage(meanMap2) + ) { + + int cols = inRaster1.getCols(); + int rows = inRaster1.getRows(); + + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + float[] kernelData = kernel.getKernelData(); + + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { + + double inputValue1 = inRaster1.getValue(c, r); + if (inRaster1.isNovalue(inputValue1)) { + continue; + } + double inputValue2 = inRaster2.getValue(c, r); + if (inRaster2.isNovalue(inputValue2)) { + continue; + } + + double meanValue1 = meanRaster1.getValue(c,r); + double meanValue2 = meanRaster2.getValue(c,r); + if (meanRaster1.isNovalue(meanValue1) || meanRaster2.isNovalue(meanValue2)) { + continue; + } + + if (doConstant) { + inputValue1 = 1.0; + inputValue2 = 1.0; + } + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value1 = inRaster1.getValue(c + kc, r + kr); + if (inRaster1.isNovalue(value1)) { + value1 = 0; + } + double value2 = inRaster2.getValue(c + kc, r + kr); + if (inRaster2.isNovalue(value2)) { + value2 = 0; + } + try { + outputValue = outputValue + kernelData[k++] * (value1 - meanValue1) * (value2 - meanValue2); + } catch (Exception e) { + throw new RuntimeException(e); + } + } + } + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } + } + pm.worked(1); + }); + pm.done(); + + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); + } + + return outputRaster.buildCoverage(); + } + } + + private GridCoverage2D windowVariance(GridCoverage2D inMap, GridCoverage2D meanMap) throws Exception { + + checkNull(inMap); + checkNull(meanMap); + + try (HMRaster inRaster = HMRaster.fromGridCoverage(inMap); HMRaster meanRaster = HMRaster.fromGridCoverage(meanMap)) { + + int cols = inRaster.getCols(); + int rows = inRaster.getRows(); + + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap).build(); + float[] kernelData = kernel.getKernelData(); + + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { + double inputValue = inRaster.getValue(c, r); + if (inRaster.isNovalue(inputValue)) { + continue; + } + + double meanValue = meanRaster.getValue(c,r); + if (meanRaster.isNovalue(meanValue)) { + continue; + } + + if (doConstant) + inputValue = 1.0; + + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value = inRaster.getValue(c + kc, r + kr); + if (inRaster.isNovalue(value)) { + value = 0; + } + try { + outputValue = outputValue + kernelData[k++] * Math.pow((value - meanValue),2); + } catch (Exception e) { + throw new RuntimeException(e); + } + } + } + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } + } + pm.worked(1); + }); + pm.done(); + + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); + } + + return outputRaster.buildCoverage(); + } + } + + + private GridCoverage2D windowMean(GridCoverage2D inMap) throws Exception { + + checkNull(inMap); + + try (HMRaster inRaster = HMRaster.fromGridCoverage(inMap)) { + + int cols = inRaster.getCols(); + int rows = inRaster.getRows(); + + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap).build(); + float[] kernelData = kernel.getKernelData(); + + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { + double inputValue = inRaster.getValue(c, r); + if (inRaster.isNovalue(inputValue)) { + continue; + } + + if (doConstant) + inputValue = 1.0; + + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value = inRaster.getValue(c + kc, r + kr); + if (inRaster.isNovalue(value)) { + value = 0; + } + try { + outputValue = outputValue + kernelData[k++] * value; + } catch (Exception e) { + throw new RuntimeException(e); + } + } + } + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } + } + pm.worked(1); + }); + pm.done(); + + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); + } + + return outputRaster.buildCoverage(); + } + } + + private KernelFactory.ValueType getKernelType(){ + KernelFactory.ValueType type = KernelFactory.ValueType.EPANECHNIKOV; + switch( pKernel ) { + case 0: + type = KernelFactory.ValueType.BINARY; + break; + case 1: + type = KernelFactory.ValueType.COSINE; + break; + case 2: + type = KernelFactory.ValueType.DISTANCE; + break; + case 4: + type = KernelFactory.ValueType.GAUSSIAN; + break; + case 5: + type = KernelFactory.ValueType.INVERSE_DISTANCE; + break; + case 6: + type = KernelFactory.ValueType.QUARTIC; + break; + case 7: + type = KernelFactory.ValueType.TRIANGULAR; + break; + case 8: + type = KernelFactory.ValueType.TRIWEIGHT; + break; + } + return type; + } + + + + + + public static int getCodeForType( KernelFactory.ValueType type ) { + switch( type ) { + case BINARY: + return 0; + case COSINE: + return 1; + case DISTANCE: + return 2; + case EPANECHNIKOV: + return 3; + case GAUSSIAN: + return 4; + case INVERSE_DISTANCE: + return 5; + case QUARTIC: + return 6; + case TRIANGULAR: + return 7; + case TRIWEIGHT: + return 8; + default: + throw new ModelsIllegalargumentException("No kernel type: " + type, "OmsKernelDensity"); + } + } + + +} From ace0ced99f33c4cc0b6a0828310382ceb2ddc343 Mon Sep 17 00:00:00 2001 From: Diego Date: Tue, 30 May 2023 10:17:01 +0200 Subject: [PATCH 2/3] Moved null checks and GridCoverage conversion to HMRaster --- .../OmsStructuralSimilarity.java | 352 ++++++++---------- 1 file changed, 165 insertions(+), 187 deletions(-) diff --git a/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java index a157d72c5..90e1f7e47 100644 --- a/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java +++ b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java @@ -48,20 +48,21 @@ public class OmsStructuralSimilarity extends HMModel { @Execute public void process() throws Exception { - GridCoverage2D meanMap1 = windowMean(inMap1); - GridCoverage2D meanMap2 = windowMean(inMap2); + checkNull(inMap1); + checkNull(inMap2); + + try(HMRaster inRaster1 = HMRaster.fromGridCoverage(inMap1); + HMRaster inRaster2 = HMRaster.fromGridCoverage(inMap2) + ){ - GridCoverage2D varianceMap1 = windowVariance(inMap1,meanMap1); - GridCoverage2D varianceMap2 = windowVariance(inMap2,meanMap2); + HMRaster meanRaster1 = windowMean(inRaster1); + HMRaster meanRaster2 = windowMean(inRaster2); - GridCoverage2D covarianceMap = windowCovariance(inMap1,inMap2,meanMap1,meanMap2); + HMRaster varianceRaster1 = windowVariance(inRaster1,meanRaster1); + HMRaster varianceRaster2 = windowVariance(inRaster2,meanRaster2); + + HMRaster covarianceRaster = windowCovariance(inRaster1,inRaster2,meanRaster1,meanRaster2); - try (HMRaster meanRaster1 = HMRaster.fromGridCoverage(meanMap1); - HMRaster meanRaster2 = HMRaster.fromGridCoverage(meanMap2); - HMRaster varianceRaster1 = HMRaster.fromGridCoverage(varianceMap1); - HMRaster varianceRaster2 = HMRaster.fromGridCoverage(varianceMap2); - HMRaster covarianceRaster = HMRaster.fromGridCoverage(covarianceMap) - ) { int cols = meanRaster1.getCols(); int rows = meanRaster1.getRows(); @@ -126,223 +127,204 @@ public void process() throws Exception { } - private GridCoverage2D windowCovariance(GridCoverage2D inMap1, GridCoverage2D inMap2, GridCoverage2D meanMap1, GridCoverage2D meanMap2) throws Exception { + private HMRaster windowCovariance(HMRaster inRaster1, HMRaster inRaster2, HMRaster meanRaster1, HMRaster meanRaster2) throws Exception { - checkNull(inMap1); - checkNull(inMap2); - checkNull(meanMap1); - checkNull(meanMap2); - try (HMRaster inRaster1 = HMRaster.fromGridCoverage(inMap1); - HMRaster meanRaster1 = HMRaster.fromGridCoverage(meanMap1); - HMRaster inRaster2 = HMRaster.fromGridCoverage(inMap2); - HMRaster meanRaster2 = HMRaster.fromGridCoverage(meanMap2) - ) { + int cols = inRaster1.getCols(); + int rows = inRaster1.getRows(); - int cols = inRaster1.getCols(); - int rows = inRaster1.getRows(); + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); - KernelFactory.ValueType type = getKernelType(); - KernelJAI kernel = KernelFactory.createCircle(pRadius, type); - HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); - float[] kernelData = kernel.getKernelData(); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); - pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + float[] kernelData = kernel.getKernelData(); - IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { - if(errorOccurred) { - return; - } - for( int c = pRadius; c < cols - pRadius; c++ ) { + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); - double inputValue1 = inRaster1.getValue(c, r); - if (inRaster1.isNovalue(inputValue1)) { - continue; - } - double inputValue2 = inRaster2.getValue(c, r); - if (inRaster2.isNovalue(inputValue2)) { - continue; - } + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { - double meanValue1 = meanRaster1.getValue(c,r); - double meanValue2 = meanRaster2.getValue(c,r); - if (meanRaster1.isNovalue(meanValue1) || meanRaster2.isNovalue(meanValue2)) { - continue; - } + double inputValue1 = inRaster1.getValue(c, r); + if (inRaster1.isNovalue(inputValue1)) { + continue; + } + double inputValue2 = inRaster2.getValue(c, r); + if (inRaster2.isNovalue(inputValue2)) { + continue; + } - if (doConstant) { - inputValue1 = 1.0; - inputValue2 = 1.0; - } - int k = 0; - double outputValue = 0.0; - for( int kr = -pRadius; kr <= pRadius; kr++ ) { - for( int kc = -pRadius; kc <= pRadius; kc++ ) { - double value1 = inRaster1.getValue(c + kc, r + kr); - if (inRaster1.isNovalue(value1)) { - value1 = 0; - } - double value2 = inRaster2.getValue(c + kc, r + kr); - if (inRaster2.isNovalue(value2)) { - value2 = 0; - } - try { - outputValue = outputValue + kernelData[k++] * (value1 - meanValue1) * (value2 - meanValue2); - } catch (Exception e) { - throw new RuntimeException(e); - } + double meanValue1 = meanRaster1.getValue(c,r); + double meanValue2 = meanRaster2.getValue(c,r); + if (meanRaster1.isNovalue(meanValue1) || meanRaster2.isNovalue(meanValue2)) { + continue; + } + + if (doConstant) { + inputValue1 = 1.0; + inputValue2 = 1.0; + } + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value1 = inRaster1.getValue(c + kc, r + kr); + if (inRaster1.isNovalue(value1)) { + value1 = 0; + } + double value2 = inRaster2.getValue(c + kc, r + kr); + if (inRaster2.isNovalue(value2)) { + value2 = 0; + } + try { + outputValue = outputValue + kernelData[k++] * (value1 - meanValue1) * (value2 - meanValue2); + } catch (Exception e) { + throw new RuntimeException(e); } - } - try { - outputRaster.setValue(c, r, outputValue); - } catch (IOException e) { - errorOccurred = true; - errorMessage = e.getLocalizedMessage(); } } - pm.worked(1); - }); - pm.done(); - - if (errorOccurred) { - throw new ModelsRuntimeException(errorMessage, this); + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } } + pm.worked(1); + }); + pm.done(); - return outputRaster.buildCoverage(); + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); } - } - private GridCoverage2D windowVariance(GridCoverage2D inMap, GridCoverage2D meanMap) throws Exception { + return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + } - checkNull(inMap); - checkNull(meanMap); + private HMRaster windowVariance(HMRaster inRaster, HMRaster meanRaster) throws Exception { - try (HMRaster inRaster = HMRaster.fromGridCoverage(inMap); HMRaster meanRaster = HMRaster.fromGridCoverage(meanMap)) { + int cols = inRaster.getCols(); + int rows = inRaster.getRows(); - int cols = inRaster.getCols(); - int rows = inRaster.getRows(); + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + float[] kernelData = kernel.getKernelData(); - KernelFactory.ValueType type = getKernelType(); - KernelJAI kernel = KernelFactory.createCircle(pRadius, type); - HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap).build(); - float[] kernelData = kernel.getKernelData(); + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); - pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { + double inputValue = inRaster.getValue(c, r); + if (inRaster.isNovalue(inputValue)) { + continue; + } - IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { - if(errorOccurred) { - return; + double meanValue = meanRaster.getValue(c,r); + if (meanRaster.isNovalue(meanValue)) { + continue; } - for( int c = pRadius; c < cols - pRadius; c++ ) { - double inputValue = inRaster.getValue(c, r); - if (inRaster.isNovalue(inputValue)) { - continue; - } - double meanValue = meanRaster.getValue(c,r); - if (meanRaster.isNovalue(meanValue)) { - continue; - } + if (doConstant) + inputValue = 1.0; - if (doConstant) - inputValue = 1.0; - - int k = 0; - double outputValue = 0.0; - for( int kr = -pRadius; kr <= pRadius; kr++ ) { - for( int kc = -pRadius; kc <= pRadius; kc++ ) { - double value = inRaster.getValue(c + kc, r + kr); - if (inRaster.isNovalue(value)) { - value = 0; - } - try { - outputValue = outputValue + kernelData[k++] * Math.pow((value - meanValue),2); - } catch (Exception e) { - throw new RuntimeException(e); - } + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value = inRaster.getValue(c + kc, r + kr); + if (inRaster.isNovalue(value)) { + value = 0; + } + try { + outputValue = outputValue + kernelData[k++] * Math.pow((value - meanValue),2); + } catch (Exception e) { + throw new RuntimeException(e); } - } - try { - outputRaster.setValue(c, r, outputValue); - } catch (IOException e) { - errorOccurred = true; - errorMessage = e.getLocalizedMessage(); } } - pm.worked(1); - }); - pm.done(); - - if (errorOccurred) { - throw new ModelsRuntimeException(errorMessage, this); + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } } + pm.worked(1); + }); + pm.done(); - return outputRaster.buildCoverage(); + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); } - } - - private GridCoverage2D windowMean(GridCoverage2D inMap) throws Exception { + return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + } - checkNull(inMap); - try (HMRaster inRaster = HMRaster.fromGridCoverage(inMap)) { + private HMRaster windowMean(HMRaster inRaster) throws Exception { - int cols = inRaster.getCols(); - int rows = inRaster.getRows(); + int cols = inRaster.getCols(); + int rows = inRaster.getRows(); - KernelFactory.ValueType type = getKernelType(); - KernelJAI kernel = KernelFactory.createCircle(pRadius, type); - HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap).build(); - float[] kernelData = kernel.getKernelData(); + KernelFactory.ValueType type = getKernelType(); + KernelJAI kernel = KernelFactory.createCircle(pRadius, type); + HMRaster outputRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); + float[] kernelData = kernel.getKernelData(); - pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); + pm.beginTask("Estimating kernel density...", cols - 2 * pRadius); - IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { - if(errorOccurred) { - return; + IntStream.range(pRadius, rows - pRadius).parallel().forEach(r -> { + if(errorOccurred) { + return; + } + for( int c = pRadius; c < cols - pRadius; c++ ) { + double inputValue = inRaster.getValue(c, r); + if (inRaster.isNovalue(inputValue)) { + continue; } - for( int c = pRadius; c < cols - pRadius; c++ ) { - double inputValue = inRaster.getValue(c, r); - if (inRaster.isNovalue(inputValue)) { - continue; - } - if (doConstant) - inputValue = 1.0; - - int k = 0; - double outputValue = 0.0; - for( int kr = -pRadius; kr <= pRadius; kr++ ) { - for( int kc = -pRadius; kc <= pRadius; kc++ ) { - double value = inRaster.getValue(c + kc, r + kr); - if (inRaster.isNovalue(value)) { - value = 0; - } - try { - outputValue = outputValue + kernelData[k++] * value; - } catch (Exception e) { - throw new RuntimeException(e); - } + if (doConstant) + inputValue = 1.0; + + int k = 0; + double outputValue = 0.0; + for( int kr = -pRadius; kr <= pRadius; kr++ ) { + for( int kc = -pRadius; kc <= pRadius; kc++ ) { + double value = inRaster.getValue(c + kc, r + kr); + if (inRaster.isNovalue(value)) { + value = 0; + } + try { + outputValue = outputValue + kernelData[k++] * value; + } catch (Exception e) { + throw new RuntimeException(e); } - } - try { - outputRaster.setValue(c, r, outputValue); - } catch (IOException e) { - errorOccurred = true; - errorMessage = e.getLocalizedMessage(); } } - pm.worked(1); - }); - pm.done(); - - if (errorOccurred) { - throw new ModelsRuntimeException(errorMessage, this); + try { + outputRaster.setValue(c, r, outputValue); + } catch (IOException e) { + errorOccurred = true; + errorMessage = e.getLocalizedMessage(); + } } + pm.worked(1); + }); + pm.done(); - return outputRaster.buildCoverage(); + if (errorOccurred) { + throw new ModelsRuntimeException(errorMessage, this); } + + return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + } private KernelFactory.ValueType getKernelType(){ @@ -375,11 +357,7 @@ private KernelFactory.ValueType getKernelType(){ } return type; } - - - - - + public static int getCodeForType( KernelFactory.ValueType type ) { switch( type ) { case BINARY: From ed050b92824938a14cbbede5caa853e9db1c151c Mon Sep 17 00:00:00 2001 From: Diego Date: Tue, 30 May 2023 11:11:21 +0200 Subject: [PATCH 3/3] correct output of window functions, include raster ranges in calculation --- .../OmsStructuralSimilarity.java | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java index 90e1f7e47..af9a86396 100644 --- a/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java +++ b/hmachine/src/main/java/org/hortonmachine/hmachine/modules/statistics/kerneldensity/OmsStructuralSimilarity.java @@ -12,6 +12,8 @@ import java.io.IOException; import java.util.stream.IntStream; +import static org.hortonmachine.gears.modules.r.summary.OmsRasterSummary.getMinMax; + public class OmsStructuralSimilarity extends HMModel { @@ -63,12 +65,16 @@ public void process() throws Exception { HMRaster covarianceRaster = windowCovariance(inRaster1,inRaster2,meanRaster1,meanRaster2); + //TODO: check that maps to be compared are of the same dimensions. int cols = meanRaster1.getCols(); int rows = meanRaster1.getRows(); + + double[] range1 = getMinMax(inMap1); + double[] range2 = getMinMax(inMap2); + double range = Math.max(range1[1],range2[1]) - Math.min(range1[0],range2[0]); - // TODO: replace 0.0 by the value range among the two comparable rasters. - double c1 = Math.pow(pK1 * 0.0,2); - double c2 = Math.pow(pK2 * 0.0,2); + double c1 = Math.pow(pK1 * range,2); + double c2 = Math.pow(pK2 * range,2); double c3 = 0.5*c2; HMRaster outputSSRaster = new HMRaster.HMRasterWritableBuilder().setTemplate(inMap1).build(); @@ -201,7 +207,7 @@ private HMRaster windowCovariance(HMRaster inRaster1, HMRaster inRaster2, HMRast throw new ModelsRuntimeException(errorMessage, this); } - return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + return outputRaster; } private HMRaster windowVariance(HMRaster inRaster, HMRaster meanRaster) throws Exception { @@ -264,7 +270,7 @@ private HMRaster windowVariance(HMRaster inRaster, HMRaster meanRaster) throws E throw new ModelsRuntimeException(errorMessage, this); } - return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + return outputRaster; } @@ -323,7 +329,7 @@ private HMRaster windowMean(HMRaster inRaster) throws Exception { throw new ModelsRuntimeException(errorMessage, this); } - return HMRaster.fromGridCoverage(outputRaster.buildCoverage()); + return outputRaster; } @@ -357,7 +363,7 @@ private KernelFactory.ValueType getKernelType(){ } return type; } - + public static int getCodeForType( KernelFactory.ValueType type ) { switch( type ) { case BINARY: