diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/AbstractGeometryRasterizer.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/AbstractGeometryRasterizer.java new file mode 100644 index 0000000..adcc315 --- /dev/null +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/AbstractGeometryRasterizer.java @@ -0,0 +1,18 @@ +package org.geotools.process.vector; + +public abstract class AbstractGeometryRasterizer implements GeometryRasterizer { + + GridTransform trans; + RasterHandler handler; + + @Override + public void setTransformation(GridTransform trans) { + this.trans = trans; + } + + @Override + public void setHandler(RasterHandler handler) { + this.handler = handler; + } + +} diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/BoundingBoxGeometryRasterizer.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/BoundingBoxGeometryRasterizer.java new file mode 100644 index 0000000..9401659 --- /dev/null +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/BoundingBoxGeometryRasterizer.java @@ -0,0 +1,49 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2012, Open Source Geospatial Foundation (OSGeo) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library 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 + * Lesser General Public License for more details. + */ + +package org.geotools.process.vector; + +import com.vividsolutions.jts.geom.Envelope; +import com.vividsolutions.jts.geom.Geometry; + +/** + * Simple rasterizer which covers the entire bounding box of the geometry. + * + * @author Kevin Smith, OpenGeo + * + */ +public class BoundingBoxGeometryRasterizer extends AbstractGeometryRasterizer { + + @Override + public void rasterize(Geometry g, Object userData) { + final Envelope geomBounds = g.getEnvelopeInternal(); + + + final int minI = trans.safeI(geomBounds.getMinX()); + final int minJ = trans.safeJ(geomBounds.getMinY()); + final int maxI = trans.safeI(geomBounds.getMaxX()); + final int maxJ = trans.safeJ(geomBounds.getMaxY()); + + + for (int i = minI; i <= maxI; i++) { + for (int j = minJ; j <= maxJ; j++) { + handler.point(i, j, userData, null); + } + } + } + +} diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/CentroidGeometryRasterizer.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/CentroidGeometryRasterizer.java new file mode 100644 index 0000000..c3fefe3 --- /dev/null +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/CentroidGeometryRasterizer.java @@ -0,0 +1,43 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2012, Open Source Geospatial Foundation (OSGeo) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library 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 + * Lesser General Public License for more details. + */ + +package org.geotools.process.vector; + +import com.vividsolutions.jts.geom.Geometry; +import com.vividsolutions.jts.geom.Point; + +/** + * Simple geometry rasterizer which sets a single pixel at the centroid of the geometry. + * @author Kevin Smith, OpenGeo + * + */ +public class CentroidGeometryRasterizer extends AbstractGeometryRasterizer { + + @Override + public void rasterize(Geometry g, Object userData) { + Point p = g.getCentroid(); + + int i = trans.i(p.getX()); + int j = trans.j(p.getY()); + + if(i<0 || i>= trans.getXSize()) return; + if(j<0 || j>= trans.getYSize()) return; + + handler.point(i, j, userData, null); + } + +} diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GeometryRasterizer.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GeometryRasterizer.java new file mode 100644 index 0000000..5af76d8 --- /dev/null +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GeometryRasterizer.java @@ -0,0 +1,66 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2012, Open Source Geospatial Foundation (OSGeo) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library 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 + * Lesser General Public License for more details. + */ + +package org.geotools.process.vector; + +import com.vividsolutions.jts.geom.Geometry; + + +/** + * Interface for an implementation of a rasterization algorithm + * + * @author Kevin Smith, OpenGeo + * + */ +public interface GeometryRasterizer { + + /** + * Functor to abstract the raster implementation being written to + * + */ + public interface RasterHandler { + /** + * + * @param i the horizontal grid ordinate + * @param j the vertical grid ordinate + * @param userData the user data passed in with the geometry + * @param algorithmData extra information generated by the rasterization algorithm + */ + public void point(int i, int j, Object userData, Object algorithmData); + + } + + + /** + * Set the GridTransformation used to transform from spatial to pixel coordinates + * @param trans + */ + public void setTransformation(GridTransform trans); + + /** + * Set the handler used to write to the raster + * @param handler + */ + public void setHandler(RasterHandler handler); + + /** + * Rasterize a geometry + * @param g the geometry + * @param o user data to pass to the handler + */ + public void rasterize(Geometry g, Object o); +} diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GridTransform.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GridTransform.java index f2bd2f4..d7f55f9 100644 --- a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GridTransform.java +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/GridTransform.java @@ -2,7 +2,7 @@ * GeoTools - The Open Source Java GIS Toolkit * http://geotools.org * - * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2011-2012, Open Source Geospatial Foundation (OSGeo) * (C) 2008-2011 TOPP - www.openplans.org. * * This library is free software; you can redistribute it and/or @@ -31,19 +31,20 @@ import com.vividsolutions.jts.geom.Envelope; * must check that values are in an acceptable range. * * @author Martin Davis - OpenGeo + * @author Kevin Smith - OpenGeo * */ class GridTransform { - private Envelope env; + final private Envelope env; - private int xSize; + final private int xSize; - private int ySize; + final private int ySize; - private double dx; + final private double dx; - private double dy; + final private double dy; private boolean isClamped = true; @@ -58,41 +59,117 @@ class GridTransform { this.env = env; this.xSize = xSize; this.ySize = ySize; - dx = env.getWidth() / (xSize - 1); - dy = env.getHeight() / (ySize - 1); + dx = env.getWidth() / (xSize); + dy = env.getHeight() / (ySize); } + /** + * Returns a GridTransform with the same scale factors, but with an expanded envelope. Negative + * margins may be used to contract the envelope. + * + * For points within both envelopes or when clamping is off, the pixel coordinates in the new + * grid will be increased by the minI/J parameters given here, relative to the old grid. + * + * @param minI margin in pixels to expand the lower I/X bound by. + * @param minJ margin in pixels to expand the lower J/Y bound by. + * @param maxI margin in pixels to expand the upper I/X bound by. + * @param maxJ margin in pixels to expand the upper J/Y bound by. + * @return + */ + public GridTransform expand(int minI, int minJ, int maxI, int maxJ){ + int newXSize = this.xSize+minI+maxI; + int newYSize = this.ySize+minJ+maxJ; + + Envelope newEnv = new Envelope( + env.getMinX()-minI*dx, + env.getMaxX()+maxI*dx, + env.getMinY()-minJ*dy, + env.getMaxY()+maxJ*dy); + + GridTransform newGrid = new GridTransform(newEnv, newXSize, newYSize); + newGrid.setClamp(isClamped); // New Transform should have same clamped state as old one. + + return newGrid; + } + + /** + * Returns a GridTransform equivalent to this one, but with an expanded envelope. A negative + * margin may be used to contract the envelope. + * + * For points within both envelopes or when clamping is off, the pixel coordinates in the new + * grid will be increased by the margin parameter given here, relative to the old grid. + * + * @param margin margin in pixels to expand the bounds by. + * @return + */ + public GridTransform expand(int margin){ + + return this.expand(margin, margin, margin, margin); + } + + /** * Sets whether to clamp outputs from transform to input envelope. * Default is to clamp the outputs. * + * Clamped values will be set to -1 or x/ySize. + * * @param isClamped true if input is to be clamped */ public void setClamp(boolean isClamped) { this.isClamped = isClamped; } - + /** - * Computes the X ordinate of the i'th grid column. + * The envelope of the GridTransformation + * @return + */ + public Envelope getEnv() { + return env; + } + + + /** + * The width of the grid + * @return + */ + public int getXSize() { + return xSize; + } + + + /** + * The height of the grid + * @return + */ + public int getYSize() { + return ySize; + } + + + /** + * Computes the X ordinate of the centre of the Ith grid column. * @param i the index of a grid column * @return the X ordinate of the column */ public double x(int i) { if (i >= xSize - 1) return env.getMaxX(); - return env.getMinX() + i * dx; + // Get the centre of the cell + return env.getMinX() + (2*i+1) * dx / 2; } /** - * Computes the Y ordinate of the i'th grid row. + * Computes the Y ordinate of the centre of the Jth grid row. * @param j the index of a grid row * @return the Y ordinate of the row */ public double y(int j) { if (j >= ySize - 1) return env.getMaxY(); - return env.getMinY() + j * dy; + // Get the centre of the cell + return env.getMinY() + (2*j+1) * dy /2; } /** @@ -129,4 +206,31 @@ class GridTransform { return j; } + /** + * Computes the column index of an X ordinate. Out of bounds ordinates are set to the nearest edge. + * @param x + * @return + */ + public int safeI(double x){ + int i = i(x); + if(isClamped){ + if(i<0) return 0; + if(i>=getXSize()) return i-1; + } + return i; + } + + /** + * Computes the row index of a Y ordinate. Out of bounds ordinates are set to the nearest edge. + * @param y + * @return + */ + public int safeJ(double y){ + int j = j(y); + if(isClamped){ + if(j<0) return 0; + if(j>=getYSize()) return j-1; + } + return j; + } } \ No newline at end of file diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapProcess.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapProcess.java index 5d97611..10619a3 100644 --- a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapProcess.java +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapProcess.java @@ -2,7 +2,7 @@ * GeoTools - The Open Source Java GIS Toolkit * http://geotools.org * - * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2011-2012, Open Source Geospatial Foundation (OSGeo) * (C) 2008-2011 TOPP - www.openplans.org. * * This library is free software; you can redistribute it and/or @@ -30,12 +30,14 @@ import org.geotools.data.simple.SimpleFeatureIterator; import org.geotools.factory.Hints; import org.geotools.filter.text.cql2.CQLException; import org.geotools.filter.text.ecql.ECQL; +import org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer; import org.geotools.geometry.jts.ReferencedEnvelope; import org.geotools.process.ProcessException; import org.geotools.process.factory.DescribeParameter; import org.geotools.process.factory.DescribeProcess; import org.geotools.process.factory.DescribeResult; import org.geotools.referencing.CRS; +import org.geotools.referencing.operation.projection.ProjectionException; import org.opengis.coverage.grid.GridCoverage; import org.opengis.coverage.grid.GridGeometry; import org.opengis.feature.simple.SimpleFeature; @@ -44,9 +46,11 @@ import org.opengis.filter.expression.Expression; import org.opengis.referencing.FactoryException; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.TransformException; import org.opengis.util.ProgressListener; import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.geom.Geometry; /** @@ -129,12 +133,14 @@ public class HeatmapProcess implements VectorProcess { @DescribeParameter(name = "radiusPixels", description = "Radius of the density kernel in pixels") Integer argRadiusPixels, @DescribeParameter(name = "weightAttr", description = "Name of the attribute to use for data point weight", min = 0, max = 1) String valueAttr, @DescribeParameter(name = "pixelsPerCell", description = "Resolution at which to compute the heatmap (in pixels). Default = 1", min = 0, max = 1) Integer argPixelsPerCell, + @DescribeParameter(name = "rasterizeMode", description = "Method used to rasterize complex geometries (\"centroid\", \"envelope\"). Default = \"centroid\"", min = 0, max = 1) String argRasterizeMode, + @DescribeParameter(name = "normalizeMinimum", description = "Normalize so that the minimum value is 0. Default = false", min = 0, max = 1) Boolean argNormalizeMinimum, // output image parameters @DescribeParameter(name = "outputBBOX", description = "Bounding box of the output") ReferencedEnvelope argOutputEnv, @DescribeParameter(name = "outputWidth", description = "Width of output raster in pixels") Integer argOutputWidth, @DescribeParameter(name = "outputHeight", description = "Height of output raster in pixels") Integer argOutputHeight, - + ProgressListener monitor) throws ProcessException { /** @@ -152,7 +158,22 @@ public class HeatmapProcess implements VectorProcess { gridWidth = outputWidth / pixelsPerCell; gridHeight = outputHeight / pixelsPerCell; } - + + GeometryRasterizer rasterizer = new CentroidGeometryRasterizer(); + if(argRasterizeMode!=null) { + if(argRasterizeMode.equalsIgnoreCase("centroid")) { + rasterizer = new CentroidGeometryRasterizer(); + } else if(argRasterizeMode.equalsIgnoreCase("envelope")) { + rasterizer = new BoundingBoxGeometryRasterizer(); + } else { + throw new ProcessException(String.format("Unknown rasterizeMode for Heatmap: \"%s\"", argRasterizeMode)); + } + } + + boolean normalizeMinimum = false; + if(argNormalizeMinimum!=null) { + normalizeMinimum = argNormalizeMinimum; + } /** * Compute transform to convert input coords into output CRS */ @@ -183,9 +204,9 @@ public class HeatmapProcess implements VectorProcess { * -------------- Extract the input observation points ----------- */ HeatmapSurface heatMap = new HeatmapSurface(radiusCells, argOutputEnv, gridWidth, - gridHeight); + gridHeight, normalizeMinimum); try { - extractPoints(obsFeatures, valueAttr, trans, heatMap); + extractPoints(obsFeatures, valueAttr, trans, heatMap, rasterizer); } catch (CQLException e) { throw new ProcessException(e); } @@ -327,18 +348,26 @@ public class HeatmapProcess implements VectorProcess { } public static void extractPoints(SimpleFeatureCollection obsPoints, String attrName, - MathTransform trans, HeatmapSurface heatMap) throws CQLException { + MathTransform trans, final HeatmapSurface heatMap, GeometryRasterizer rasterizer) throws CQLException { Expression attrExpr = null; if (attrName != null) { attrExpr = ECQL.toExpression(attrName); } SimpleFeatureIterator obsIt = obsPoints.features(); - - double[] srcPt = new double[2]; - double[] dstPt = new double[2]; - - int i = 0; + + GeometryCoordinateSequenceTransformer transformer = new GeometryCoordinateSequenceTransformer(); + transformer.setMathTransform(trans); + + rasterizer.setTransformation(heatMap.gridTrans); + rasterizer.setHandler(new GeometryRasterizer.RasterHandler() { + + @Override + public void point(int i, int j, Object userData, Object algorithmData) { + heatMap.addPoint(i, j, (Float) userData); + } + }); + try { while (obsIt.hasNext()) { SimpleFeature feature = obsIt.next(); @@ -349,20 +378,16 @@ public class HeatmapProcess implements VectorProcess { if (attrExpr != null) { val = getPointValue(feature, attrExpr); } - - // get the point location from the geometry + Geometry geom = (Geometry) feature.getDefaultGeometry(); - Coordinate p = getPoint(geom); - srcPt[0] = p.x; - srcPt[1] = p.y; - trans.transform(srcPt, 0, dstPt, 0, 1); - Coordinate pobs = new Coordinate(dstPt[0], dstPt[1], val); - - heatMap.addPoint(pobs.x, pobs.y, val); + + geom=transformer.transform((Geometry) feature.getDefaultGeometry()); + rasterizer.rasterize(geom, new Float(val)); + } catch (ProjectionException e){ + // Keep trying with other geometries } catch (Exception e) { - // just carry on for now (debugging) - // throw new ProcessException("Expression " + attrExpr + - // " failed to evaluate to a numeric value", e); + throw new ProcessException("Expression " + attrExpr + + " failed to evaluate to a numeric value", e); } } } finally { diff --git a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapSurface.java b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapSurface.java index 96c85f1..e8dde44 100644 --- a/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapSurface.java +++ b/modules/unsupported/process-feature/src/main/java/org/geotools/process/vector/HeatmapSurface.java @@ -2,7 +2,7 @@ * GeoTools - The Open Source Java GIS Toolkit * http://geotools.org * - * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2011-2012, Open Source Geospatial Foundation (OSGeo) * (C) 2008-2011 TOPP - www.openplans.org. * * This library is free software; you can redistribute it and/or @@ -15,6 +15,7 @@ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. */ + package org.geotools.process.vector; import com.vividsolutions.jts.geom.Envelope; @@ -39,17 +40,19 @@ public class HeatmapSurface { */ private static final int GAUSSIAN_APPROX_ITER = 4; - private Envelope srcEnv; + final Envelope srcEnv; - private int xSize; + final int xSize; - private int ySize; + final int ySize; - private GridTransform gridTrans; + GridTransform gridTrans; private float[][] grid; - private int kernelRadiusGrid; + final int kernelRadiusGrid; + + final boolean normalizeMinimum; /** * Creates a new heatmap surface. @@ -59,35 +62,27 @@ public class HeatmapSurface { * @param xSize the width of the output grid * @param ySize the height of the output grid */ - public HeatmapSurface(int kernelRadius, Envelope srcEnv, int xSize, int ySize) { + public HeatmapSurface(int kernelRadius, Envelope srcEnv, int xSize, int ySize, boolean normalizeMinimum) { // radius must be non-negative this.kernelRadiusGrid = Math.max(kernelRadius, 0); this.srcEnv = srcEnv; this.xSize = xSize; this.ySize = ySize; + this.normalizeMinimum = normalizeMinimum; init(); } private void init() { gridTrans = new GridTransform(srcEnv, xSize, ySize); - /** - * Do NOT clamp transform output, since - * the actual target grid is larger than the source Envelope, - * due to the required buffering. - * This means that transform outputs MUST be checked for validity. - */ - gridTrans.setClamp(false); - - int xSizeExp = xSize + 2 * kernelRadiusGrid; - int ySizeExp = ySize + 2 * kernelRadiusGrid; + gridTrans=gridTrans.expand(kernelRadiusGrid); - grid = new float[xSizeExp][ySizeExp]; + grid = new float[gridTrans.getXSize()][gridTrans.getYSize()]; } /** - * Adds a new data point to the surface. Data points can be coincident. + * Adds a new data point to the surface in map space. Data points can be coincident. * * @param x the X ordinate of the point * @param y the Y ordinate of the point @@ -96,16 +91,28 @@ public class HeatmapSurface { public void addPoint(double x, double y, double value) { /** - * Input points are converted to grid space, and offset by the grid expansion offset + * Input points are converted to grid space */ - int gi = gridTrans.i(x) + kernelRadiusGrid; - int gj = gridTrans.j(y) + kernelRadiusGrid; - + int i = gridTrans.i(x); + int j = gridTrans.j(y); + + addPoint(i, j, value); + } + + /** + * Adds a new data point to the surface in pixel space. Data points can be coincident. + * + * @param i the I (X) ordinate of the point + * @param j the J (Y) ordinate of the point + * @param value the data value of the point + */ + public void addPoint(int i, int j, double value) + { // check if point falls outside grid - skip it if so - if (gi < 0 || gi > grid.length || gj < 0 || gj > grid[0].length) + if (i < 0 || i >= grid.length || j < 0 || j >= grid[0].length) return; - grid[gi][gj] += value; + grid[i][j] += value; // System.out.println("data[" + gi + ", " + gj + "] <- " + value); } @@ -141,26 +148,28 @@ public class HeatmapSurface { int baseBoxKernelRadius = kernelRadius / GAUSSIAN_APPROX_ITER; int radiusIncBreak = kernelRadius - baseBoxKernelRadius * GAUSSIAN_APPROX_ITER; - /** - * Since Box Blur is linearly separable, can implement it by doing 2 1-D box blurs in - * different directions. Using a flipped buffer grid allows the same code to compute each - * direction, as well as preserving input grid values. - */ - // holds flipped copy of first box blur pass - float[][] grid2 = new float[ySize][xSize]; - for (int count = 0; count < GAUSSIAN_APPROX_ITER; count++) { - int boxKernelRadius = baseBoxKernelRadius; + // if the kernel radius is 0, just skip the blur entirely. + if(kernelRadius>0){ /** - * If required, increment radius to ensure sum of radii equals total kernel radius + * Since Box Blur is linearly separable, can implement it by doing 2 1-D box blurs in + * different directions. Using a flipped buffer grid allows the same code to compute each + * direction, as well as preserving input grid values. */ - if (count < radiusIncBreak) - boxKernelRadius++; - // System.out.println(boxKernelRadius); - - boxBlur(boxKernelRadius, grid, grid2); - boxBlur(boxKernelRadius, grid2, grid); + // holds flipped copy of first box blur pass + float[][] grid2 = new float[ySize][xSize]; + for (int count = 0; count < GAUSSIAN_APPROX_ITER; count++) { + int boxKernelRadius = baseBoxKernelRadius; + /** + * If required, increment radius to ensure sum of radii equals total kernel radius + */ + if (count < radiusIncBreak) + boxKernelRadius++; + // System.out.println(boxKernelRadius); + + boxBlur(boxKernelRadius, grid, grid2); + boxBlur(boxKernelRadius, grid2, grid); + } } - // testNormalizeFactor(baseBoxKernelRadius, radiusIncBreak); normalize(grid); return grid; @@ -202,6 +211,8 @@ public class HeatmapSurface { System.out.println("norm factor = " + val); } + + /** * Normalizes grid values to range [0,1] * @@ -209,18 +220,28 @@ public class HeatmapSurface { */ private void normalize(float[][] grid) { float max = Float.NEGATIVE_INFINITY; - for (int i = 0; i < grid.length; i++) { - for (int j = 0; j < grid[0].length; j++) { + float min = Float.POSITIVE_INFINITY; + + // Want to normalize over just the visible portion + + for (int i = kernelRadiusGrid; i < kernelRadiusGrid+xSize; i++) { + for (int j = kernelRadiusGrid; j < kernelRadiusGrid+ySize; j++) { if (grid[i][j] > max) max = grid[i][j]; + if (grid[i][j] < min) + min = grid[i][j]; } } - float normFactor = 1.0f / max; + float normOffset=0; + if(normalizeMinimum) { + normOffset=min; + } + float normFactor = 1.0f / (max-normOffset); for (int i = 0; i < grid.length; i++) { for (int j = 0; j < grid[0].length; j++) { - grid[i][j] *= normFactor; + grid[i][j] = (grid[i][j]-normOffset)*normFactor; } } } diff --git a/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/GridTransformTest.java b/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/GridTransformTest.java new file mode 100644 index 0000000..6d1f3fe --- /dev/null +++ b/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/GridTransformTest.java @@ -0,0 +1,81 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2012, Open Source Geospatial Foundation (OSGeo) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; + * version 2.1 of the License. + * + * This library 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 + * Lesser General Public License for more details. + */ + +package org.geotools.process.vector; + +import static org.junit.Assert.*; + +import org.junit.Test; + +import com.vividsolutions.jts.geom.Envelope; + +public class GridTransformTest { + + @Test + public void testSimple() { + Envelope env = new Envelope(500, 1600, 400, 1400); + GridTransform trans = new GridTransform(env, 11, 10); + + assertEquals(2, trans.i(750)); + assertEquals(4, trans.j(850)); + + assertEquals(750.0, trans.x(2), 0.001); + assertEquals(850.0, trans.y(4), 0.001); + + } + + @Test + public void testClamp() { + Envelope env = new Envelope(500, 1600, 400, 1400); + GridTransform trans = new GridTransform(env, 11, 10); + + assertEquals(-1, trans.i(450)); + assertEquals(-1, trans.j(350)); + assertEquals(trans.getXSize(), trans.i(1650)); + assertEquals(trans.getYSize(), trans.j(1450)); + } + + @Test + public void testExpand() { + Envelope env = new Envelope(500, 1600, 400, 1400); + GridTransform trans = new GridTransform(env, 11, 10); + GridTransform trans2 = trans.expand(2, 1, 4, 3); + + assertEquals(new Envelope(300, 2000, 300, 1700), trans2.getEnv()); + assertEquals(17, trans2.getXSize()); + assertEquals(14, trans2.getYSize()); + + + // point inside + assertEquals(4, trans2.i(750)); + assertEquals(5, trans2.j(850)); + + // points inside the new margins + assertEquals(1, trans2.i(450)); + assertEquals(0, trans2.j(350)); + assertEquals(13, trans2.i(1650)); + assertEquals(11, trans2.j(1450)); + + // points outside the new margins + assertEquals(-1, trans.i(250)); + assertEquals(-1, trans.j(150)); + assertEquals(trans2.getXSize(), trans2.i(2050)); + assertEquals(trans2.getYSize(), trans2.j(1750)); + + } + +} diff --git a/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/HeatmapProcessTest.java b/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/HeatmapProcessTest.java index cc71c6a..1f09eab 100644 --- a/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/HeatmapProcessTest.java +++ b/modules/unsupported/process-feature/src/test/java/org/geotools/process/vector/HeatmapProcessTest.java @@ -17,6 +17,7 @@ package org.geotools.process.vector; import static org.junit.Assert.assertTrue; +import static org.junit.Assert.assertEquals; import java.awt.geom.Point2D; @@ -25,6 +26,7 @@ import org.geotools.data.simple.SimpleFeatureCollection; import org.geotools.feature.DefaultFeatureCollection; import org.geotools.feature.simple.SimpleFeatureBuilder; import org.geotools.feature.simple.SimpleFeatureTypeBuilder; +import org.geotools.geometry.jts.JTS; import org.geotools.geometry.jts.ReferencedEnvelope; import org.geotools.referencing.crs.DefaultGeographicCRS; import org.junit.Test; @@ -72,6 +74,8 @@ public class HeatmapProcessTest { 20, //radius null, // weightAttr 1, // pixelsPerCell + (String) null, // rasterizeMode + (Boolean) null, // normalizeMinimum bounds, // outputEnv 100, // outputWidth 100, // outputHeight @@ -96,6 +100,117 @@ public class HeatmapProcessTest { assertTrue(far < center1 / 1000); } + + /** + * Repeat testSimpleSurface with rasterizeMode=envelope. + * Should give the same result for point features. + * + * @throws Exception + */ + @Test + public void testSimpleSurfaceEnvelope() { + + ReferencedEnvelope bounds = new ReferencedEnvelope(0, 10, 0, 10, DefaultGeographicCRS.WGS84); + Coordinate[] data = new Coordinate[] { + new Coordinate(4, 4), + new Coordinate(4, 6), + // include a coordinate outside the heatmap buffer bounds, to ensure it is filtered correctly + new Coordinate(100, 100) + }; + SimpleFeatureCollection fc = createPoints(data, bounds); + + ProgressListener monitor = null; + + HeatmapProcess process = new HeatmapProcess(); + GridCoverage2D cov = process.execute(fc, // data + 20, //radius + null, // weightAttr + 1, // pixelsPerCell + "envelope", // rasterizeMode + (Boolean) null, // normalizeMinimum + bounds, // outputEnv + 100, // outputWidth + 100, // outputHeight + monitor // monitor) + ); + + // following tests are checking for an appropriate shape for the surface + + float center1 = coverageValue(cov, 4, 4); + float center2 = coverageValue(cov, 4, 6); + float midway = coverageValue(cov, 4, 5); + float far = coverageValue(cov, 9, 9); + + // peaks are roughly equal + float peakDiff = Math.abs(center1 - center2); + assert(peakDiff < center1 / 10); + + // dip between peaks + assertTrue(midway > center1 / 2); + + // surface is flat far away + assertTrue(far < center1 / 1000); + + } + + /** + * Test minimum normalization. + * + * @throws Exception + */ + @Test + public void testNormalize() { + + ReferencedEnvelope bounds = new ReferencedEnvelope(0, 10, 0, 10, DefaultGeographicCRS.WGS84); + + + SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder(); + tb.setName("data"); + tb.setCRS(bounds.getCoordinateReferenceSystem()); + tb.add("shape", Geometry.class); + tb.add("value", Double.class); + + SimpleFeatureType type = tb.buildFeatureType(); + SimpleFeatureBuilder fb = new SimpleFeatureBuilder(type); + DefaultFeatureCollection fc = new DefaultFeatureCollection(); + + GeometryFactory factory = new GeometryFactory(new PackedCoordinateSequenceFactory()); + + Geometry point = factory.createPoint(new Coordinate(4, 4)); + fb.add(point); + fb.add(1.0); + fc.add(fb.buildFeature(null)); + + fb.add(JTS.toGeometry(bounds)); + fb.add(1.0); + fc.add(fb.buildFeature(null)); + + ProgressListener monitor = null; + + HeatmapProcess process = new HeatmapProcess(); + GridCoverage2D cov = process.execute(fc, // data + 0, //radius + null, // weightAttr + 1, // pixelsPerCell + "envelope", // rasterizeMode + true, // normalizeMinimum + bounds, // outputEnv + 100, // outputWidth + 100, // outputHeight + monitor // monitor) + ); + + // following tests are checking for an appropriate shape for the surface + + float center1 = coverageValue(cov, 4, 4); + float far = coverageValue(cov, 9, 9); + + assertEquals(1.0f, center1, 0.0001f); + + // surface is flat far away + assertEquals(0.0f, far, 0.0001f); + + } private float coverageValue(GridCoverage2D cov, double x, double y) {