Menu
▾
▴
[Geotools-commits] svn - r38814 - in trunk/modules/unsupported/process-raster/src: main/java/org/geotools/process/raster/surface test/java/org/geotools/process/raster/surface
Author: mdavis Date: 2012-06-18 11:34:42 -0700 (Mon, 18 Jun 2012) New Revision: 38814 Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapProcess.java trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapSurface.java trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/HeatmapProcessTest.java Log: Added Heatmap Rendering Transformation process Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapProcess.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapProcess.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapProcess.java 2012-06-18 18:34:42 UTC (rev 38814) @@ -0,0 +1,403 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * 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.raster.surface; + +import javax.measure.unit.NonSI; +import javax.measure.unit.SI; +import javax.measure.unit.Unit; + +import org.geotools.coverage.CoverageFactoryFinder; +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.coverage.grid.GridCoverageFactory; +import org.geotools.data.Query; +import org.geotools.data.simple.SimpleFeatureCollection; +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.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.process.gs.GSProcess; +import org.geotools.referencing.CRS; +import org.opengis.coverage.grid.GridCoverage; +import org.opengis.coverage.grid.GridGeometry; +import org.opengis.feature.simple.SimpleFeature; +import org.opengis.filter.Filter; +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.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Geometry; + +/** + * A Process that uses a {@link HeatmapSurface} to compute a heatmap surface over a set of + * irregular data points as a {@link GridCoverage}. + * Heatmaps are known more formally as <i>Multivariate Kernel Density Estimation</i>. + * <p> + * The appearance of the heatmap is controlled by the kernel radius, + * which determines the "radius of influence" of input points. + * The radius is specified by the radiusPixels parameter, + * which is in output pixels. + * Using pixels allows easy estimation of a value which will give a visually + * effective result, + * and ensures the heatmap appearance changes to match the zoom level. + * <p> + * By default each input point has weight 1. + * Optionally the weights of points may be supplied by an attribute specified by the <code>weightAttr</code> parameter. + * <p> + * All geometry types are allowed as input. For non-point geometries the centroid is used. + * <p> + * To improve performance, the surface grid can be computed at a lower resolution than the requested + * output image using the <code>pixelsPerCell</code> parameter. + * The grid is upsampled to match the required image size. Upsampling uses Bilinear + * Interpolation to maintain visual quality. This gives a large improvement in performance, with + * minimal impact on visual quality for small cell sizes (for instance, 10 pixels or less). + * <p> + * To ensure that the computed surface is stable (i.e. does not display obvious edge artifacts under + * zooming and panning), the data extent is expanded to be larger than the specified output + * extent. The expansion distance is equal to the size of <code>radiusPixels</code> in the input + * CRS. + * + * <h3>Parameters</h3> <i>M = mandatory, O = optional</i> + * <ul> + * <li><b>data</b> (M) - the FeatureCollection containing the point observations + * <li><b>radiusPixels</b> (M)- the density kernel radius, in pixels + * <li><b>weightAttr</b> (M)- the feature type attribute containing the observed surface value + * <li><b>pixelsPerCell</b> (O) - The pixels-per-cell value determines the resolution of the + * computed grid. Larger values improve performance, but degrade appearance. (Default = 1) + * <li><b>outputBBOX</b> (M) - The georeferenced bounding box of the output area + * <li><b>outputWidth</b> (M) - The width of the output raster + * <li><b>outputHeight</b> (M) - The height of the output raster + * </ul> + * The output of the process is a {@linkplain GridCoverage2D} with a single band, with cell values + * in the range [0, 1]. + * <p> + * Computation of the surface takes places in the CRS of the output. + * If the data CRS is different to the output CRS, the input points are transformed into the output CRS. + * + * <h3>Using the process as a Rendering Transformation</h3> + * + * This process can be used as a RenderingTransformation, since it implements the + * <tt>invertQuery(... Query, GridGeometry)</tt> method. In this case the <code>queryBuffer</code> + * parameter should be specified to expand the query extent appropriately. The output raster + * parameters may be provided from the request extents, using the following SLD environment + * variables: + * <ul> + * <li><b>outputBBOX</b> - env var = <tt>wms_bbox</tt> + * <li><b>outputWidth</b> - env var = <tt>wms_width</tt> + * <li><b>outputHeight</b> - env var = <tt>wms_height</tt> + * </ul> + * When used as an Rendering Transformation the data query is rewritten to expand the query BBOX, to + * ensure that enough data points are queried to make the computed surface stable under panning and + * zooming. + * + * <p> + * + * @author Martin Davis - OpenGeo + * + */ +@DescribeProcess(title = "Heatmap", description = "Computes a heatmap surface over a set of irregular data points as a GridCoverage.") +public class HeatmapProcess implements GSProcess { + + @DescribeResult(name = "result", description = "The heat map surface as a raster") + public GridCoverage2D execute( + + // process data + @DescribeParameter(name = "data", description = "Features containing the data points") SimpleFeatureCollection obsFeatures, + + // process parameters + @DescribeParameter(name = "radiusPixels", description = "Radius to use for the kernel, in pixels") Integer argRadiusPixels, + @DescribeParameter(name = "weightAttr", description = "Featuretype attribute containing the point weight value", min = 0, max = 1) String valueAttr, + @DescribeParameter(name = "pixelsPerCell", description = "Number of pixels per grid cell (default = 1)", min = 0, max = 1) Integer argPixelsPerCell, + + // output image parameters + @DescribeParameter(name = "outputBBOX", description = "Georeferenced bounding box of the output") ReferencedEnvelope argOutputEnv, + @DescribeParameter(name = "outputWidth", description = "Width of the output raster") Integer argOutputWidth, + @DescribeParameter(name = "outputHeight", description = "Height of the output raster") Integer argOutputHeight, + + ProgressListener monitor) throws ProcessException { + + /** + * -------- Extract required information from process arguments ------------- + */ + int pixelsPerCell = 1; + if (argPixelsPerCell != null && argPixelsPerCell > 1) { + pixelsPerCell = argPixelsPerCell; + } + int outputWidth = argOutputWidth; + int outputHeight = argOutputHeight; + int gridWidth = outputWidth; + int gridHeight = outputHeight; + if (pixelsPerCell > 1) { + gridWidth = outputWidth / pixelsPerCell; + gridHeight = outputHeight / pixelsPerCell; + } + + /** + * Compute transform to convert input coords into output CRS + */ + CoordinateReferenceSystem srcCRS = obsFeatures.getSchema().getCoordinateReferenceSystem(); + CoordinateReferenceSystem dstCRS = argOutputEnv.getCoordinateReferenceSystem(); + MathTransform trans = null; + try { + trans = CRS.findMathTransform(srcCRS, dstCRS); + } catch (FactoryException e) { + throw new ProcessException(e); + } + + //------------ Kernel Radius + /* + * // not used for now - only pixel radius values are supported double + * distanceConversionFactor = distanceConversionFactor(srcCRS, dstCRS); double dstRadius = + * argRadius * distanceConversionFactor; + */ + int radiusCells = 100; + if (argRadiusPixels != null) + radiusCells = argRadiusPixels; + if (pixelsPerCell > 1) { + radiusCells /= pixelsPerCell; + } + + + /** + * -------------- Extract the input observation points ----------- + */ + HeatmapSurface heatMap = new HeatmapSurface(radiusCells, argOutputEnv, gridWidth, + gridHeight); + try { + extractPoints(obsFeatures, valueAttr, trans, heatMap); + } catch (CQLException e) { + throw new ProcessException(e); + } + + /** + * --------------- Do the processing ------------------------------ + */ + // Stopwatch sw = new Stopwatch(); + // compute the heatmap at the specified resolution + float[][] heatMapGrid = heatMap.computeSurface(); + + // flip now, since grid size may be smaller + heatMapGrid = flipXY(heatMapGrid); + + // upsample to output resolution if necessary + float[][] outGrid = heatMapGrid; + if (pixelsPerCell > 1) + outGrid = upsample(heatMapGrid, -999, outputWidth, outputHeight); + + // convert to the GridCoverage2D required for output + GridCoverageFactory gcf = CoverageFactoryFinder.getGridCoverageFactory(null); + GridCoverage2D gridCov = gcf.create("Process Results", outGrid, argOutputEnv); + + // System.out.println("************** Heatmap computed in " + sw.getTimeString()); + + return gridCov; + } + + /* + * An approximate value for the length of a degree at the equator in meters. This doesn't have + * to be precise, since it is only used to convert values which are themselves rough + * approximations. + */ + private static final double METRES_PER_DEGREE = 111320; + + private static double distanceConversionFactor(CoordinateReferenceSystem srcCRS, + CoordinateReferenceSystem dstCRS) { + Unit<?> srcUnit = srcCRS.getCoordinateSystem().getAxis(0).getUnit(); + Unit<?> dstUnit = dstCRS.getCoordinateSystem().getAxis(0).getUnit(); + if (srcUnit == dstUnit) { + return 1; + } else if (srcUnit == NonSI.DEGREE_ANGLE && dstUnit == SI.METER) { + return METRES_PER_DEGREE; + } else if (srcUnit == SI.METER && dstUnit == NonSI.DEGREE_ANGLE) { + return 1.0 / METRES_PER_DEGREE; + } + throw new IllegalStateException("Unable to convert distances from " + srcUnit + " to " + + dstUnit); + } + + /** + * Flips an XY matrix along the X=Y axis, and inverts the Y axis. Used to convert from + * "map orientation" into the "image orientation" used by GridCoverageFactory. The surface + * interpolation is done on an XY grid, with Y=0 being the bottom of the space. GridCoverages + * are stored in an image format, in a YX grid with Y=0 being the top. + * + * @param grid the grid to flip + * @return the flipped grid + */ + private float[][] flipXY(float[][] grid) { + int xsize = grid.length; + int ysize = grid[0].length; + + float[][] grid2 = new float[ysize][xsize]; + for (int ix = 0; ix < xsize; ix++) { + for (int iy = 0; iy < ysize; iy++) { + int iy2 = ysize - iy - 1; + grid2[iy2][ix] = grid[ix][iy]; + } + } + return grid2; + } + + private float[][] upsample(float[][] grid, float noDataValue, int width, int height) { + BilinearInterpolator bi = new BilinearInterpolator(grid, noDataValue); + float[][] outGrid = bi.interpolate(width, height, false); + return outGrid; + } + + /** + * Given a target query and a target grid geometry returns the query to be used to read the + * input data of the process involved in rendering. In this process this method is used to: + * <ul> + * <li>determine the extent & CRS of the output grid + * <li>expand the query envelope to ensure stable surface generation + * <li>modify the query hints to ensure point features are returned + * </ul> + * Note that in order to pass validation, all parameters named here must also appear in the + * parameter list of the <tt>execute</tt> method, even if they are not used there. + * + * @param argRadiusPixels the feature type attribute that contains the observed surface value + * @param targetQuery the query used against the data source + * @param targetGridGeometry the grid geometry of the destination image + * @return The transformed query + */ + public Query invertQuery( + @DescribeParameter(name = "radiusPixels", description = "Radius to use for the kernel", min = 0, max = 1) Integer argRadiusPixels, + // output image parameters + @DescribeParameter(name = "outputBBOX", description = "Georeferenced bounding box of the output") ReferencedEnvelope argOutputEnv, + @DescribeParameter(name = "outputWidth", description = "Width of the output raster") Integer argOutputWidth, + @DescribeParameter(name = "outputHeight", description = "Height of the output raster") Integer argOutputHeight, + + Query targetQuery, GridGeometry targetGridGeometry) throws ProcessException { + + // TODO: handle different CRSes in input and output + + int radiusPixels = argRadiusPixels > 0 ? argRadiusPixels : 0; + // input parameters are required, so should be non-null + double queryBuffer = radiusPixels / pixelSize(argOutputEnv, argOutputWidth, argOutputHeight); + /* + * if (argQueryBuffer != null) { queryBuffer = argQueryBuffer; } + */ + targetQuery.setFilter(expandBBox(targetQuery.getFilter(), queryBuffer)); + + // clear properties to force all attributes to be read + // (required because the SLD processor cannot see the value attribute specified in the + // transformation) + // TODO: set the properties to read only the specified value attribute + targetQuery.setProperties(null); + + // set the decimation hint to ensure points are read + Hints hints = targetQuery.getHints(); + hints.put(Hints.GEOMETRY_DISTANCE, 0.0); + + return targetQuery; + } + + private double pixelSize(ReferencedEnvelope outputEnv, int outputWidth, int outputHeight) + { + // error-proofing + if (outputEnv.getWidth() <= 0) return 0; + // assume view is isotropic + return outputWidth / outputEnv.getWidth(); + } + + private Filter expandBBox(Filter filter, double distance) { + return (Filter) filter.accept(new BBOXExpandingFilterVisitor(distance, distance, distance, + distance), null); + } + + public static void extractPoints(SimpleFeatureCollection obsPoints, String attrName, + MathTransform trans, HeatmapSurface heatMap) 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; + try { + while (obsIt.hasNext()) { + SimpleFeature feature = obsIt.next(); + + try { + // get the weight value, if any + double val = 1; + 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); + } catch (Exception e) { + // just carry on for now (debugging) + // throw new ProcessException("Expression " + attrExpr + + // " failed to evaluate to a numeric value", e); + } + } + } finally { + obsIt.close(); + } + } + + /** + * Gets a point to represent the Geometry. If the Geometry is a point, this is returned. + * Otherwise, the centroid is used. + * + * @param g the geometry to find a point for + * @return a point representing the Geometry + */ + private static Coordinate getPoint(Geometry g) { + if (g.getNumPoints() == 1) + return g.getCoordinate(); + return g.getCentroid().getCoordinate(); + } + + /** + * Gets the value for a point from the supplied attribute. + * The value is checked for validity, + * and a default of 1 is used if necessary. + * + * @param feature the feature to extract the value from + * @param attrExpr the expression specifying the attribute to read + * @return the value for the point + */ + private static double getPointValue(SimpleFeature feature, Expression attrExpr) { + Double valObj = attrExpr.evaluate(feature, Double.class); + if (valObj != null) { + return valObj; + } + return 1; + } +} Added: trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapSurface.java =================================================================== --- trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapSurface.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/main/java/org/geotools/process/raster/surface/HeatmapSurface.java 2012-06-18 18:34:42 UTC (rev 38814) @@ -0,0 +1,266 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2011, Open Source Geospatial Foundation (OSGeo) + * (C) 2008-2011 TOPP - www.openplans.org. + * + * 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.raster.surface; + +import com.vividsolutions.jts.geom.Envelope; + +/** + * Computes a Heat Map surface from a set of irregular data points, each containing a positive + * height value. The nature of the surface is determined by a kernelRadius value, which indicates + * how far out each data points "spreads". + * <p> + * The Heatmap surface is computed as a grid (raster) of values representing the surface. For + * stability, the compute grid is expanded by the kernel radius on all four sides. This avoids + * "edge effects" from distorting the surface within the requested envelope. + * <p> + * The values in the output surface are normalized to lie in the range [0, 1]. + * + * @author Martin Davis, OpenGeo + * + */ +public class HeatmapSurface { + /** + * Number of iterations of box blur to approximate a Gaussian blur + */ + private static final int GAUSSIAN_APPROX_ITER = 4; + + private Envelope srcEnv; + + private int xSize; + + private int ySize; + + private GridTransform gridTrans; + + private float[][] grid; + + private int kernelRadiusGrid; + + /** + * Creates a new heatmap surface. + * + * @param kernelRadius the kernel radius, in grid units + * @param srcEnv the envelope defining the data space + * @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) { + // radius must be non-negative + this.kernelRadiusGrid = Math.max(kernelRadius, 0); + + this.srcEnv = srcEnv; + this.xSize = xSize; + this.ySize = ySize; + + init(); + } + + private void init() { + gridTrans = new GridTransform(srcEnv, xSize, ySize); + + int xSizeExp = xSize + 2 * kernelRadiusGrid; + int ySizeExp = ySize + 2 * kernelRadiusGrid; + + grid = new float[xSizeExp][ySizeExp]; + } + + /** + * Adds a new data point to the surface. Data points can be coincident. + * + * @param x the X ordinate of the point + * @param y the Y ordinate of the point + * @param value the data value of the point + */ + public void addPoint(double x, double y, double value) + { + /** + * Input points are converted to grid space, and offset by the grid expansion offset + */ + int gi = gridTrans.i(x) + kernelRadiusGrid; + int gj = gridTrans.j(y) + kernelRadiusGrid; + + // check if point falls outside grid - skip it if so + if (gi < 0 || gi > grid.length || gj < 0 || gj > grid[0].length) + return; + + grid[gi][gj] += value; + // System.out.println("data[" + gi + ", " + gj + "] <- " + value); + } + + /** + * Computes a grid representing the heatmap surface. The grid is structured as an XY matrix, + * with (0,0) being the bottom left corner of the data space + * + * @return a grid representing the surface + */ + public float[][] computeSurface() { + + computeHeatmap(grid, kernelRadiusGrid); + + float[][] gridOut = extractGrid(grid, kernelRadiusGrid, kernelRadiusGrid, xSize, ySize); + + return gridOut; + } + + private float[][] extractGrid(float[][] grid, int xBase, int yBase, int xSize, int ySize) { + float[][] gridExtract = new float[xSize][ySize]; + for (int i = 0; i < xSize; i++) { + for (int j = 0; j < ySize; j++) { + gridExtract[i][j] = grid[xBase + i][yBase + j]; + } + } + return gridExtract; + } + + private float[][] computeHeatmap(float[][] grid, int kernelRadius) { + int xSize = grid.length; + int ySize = grid[0].length; + + 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 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; + } + + /** + * DON'T USE This method is too simplistic to determine normalization factor. Would need to use + * a full 2D grid and smooth it to get correct value + * + * @param baseBoxKernelRadius + * @param radiusIncBreak + */ + private void testNormalizeFactor(int baseBoxKernelRadius, int radiusIncBreak) { + double val = 1.0; + 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++; + + int dia = 2 * boxKernelRadius + 1; + float kernelVal = kernelVal(boxKernelRadius); + System.out.println(boxKernelRadius + " kernel val = " + kernelVal); + + if (count == 0) { + val = val * 1 * kernelVal; + } else { + val = val * dia * kernelVal; + } + System.out.println("norm val = " + val); + if (count == 0) { + val = val * 1 * kernelVal; + } else { + val = val * dia * kernelVal; + } + } + System.out.println("norm factor = " + val); + } + + /** + * Normalizes grid values to range [0,1] + * + * @param grid + */ + 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++) { + if (grid[i][j] > max) + max = grid[i][j]; + } + } + + float normFactor = 1.0f / max; + + for (int i = 0; i < grid.length; i++) { + for (int j = 0; j < grid[0].length; j++) { + grid[i][j] *= normFactor; + } + } + } + + private float kernelVal(int kernelRadius) { + // This kernel function has been confirmed to integrate to 1 over the full radius + float val = (float) (1.0f / (2 * kernelRadius + 1)); + return val; + } + + private void boxBlur(int kernelRadius, float[][] input, float[][] output) { + int width = input.length; + int height = input[0].length; + + // init moving average total + float kernelVal = kernelVal(kernelRadius); + // System.out.println("boxblur: radius = " + kernelRadius + " kernel val = " + kernelVal); + + for (int j = 0; j < height; j++) { + + double tot = 0.0; + + for (int i = -kernelRadius; i <= kernelRadius; i++) { + if (i < 0 || i >= width) + continue; + tot += kernelVal * input[i][j]; + } + + // System.out.println(tot); + + output[j][0] = (float) tot; + + for (int i = 1; i < width; i++) { + + // update box running total + int iprev = i - 1 - kernelRadius; + if (iprev >= 0) + tot -= kernelVal * input[iprev][j]; + + int inext = i + kernelRadius; + if (inext < width) + tot += kernelVal * input[inext][j]; + + output[j][i] = (float) tot; + // if (i==49 && j==147) System.out.println("val[ " + i + ", " + j + "] = " + tot); + + } + } + } +} Added: trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/HeatmapProcessTest.java =================================================================== --- trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/HeatmapProcessTest.java (rev 0) +++ trunk/modules/unsupported/process-raster/src/test/java/org/geotools/process/raster/surface/HeatmapProcessTest.java 2012-06-18 18:34:42 UTC (rev 38814) @@ -0,0 +1,127 @@ +/* + * GeoTools - The Open Source Java GIS Toolkit + * http://geotools.org + * + * (C) 2002-2011, 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.raster.surface; + +import static org.junit.Assert.assertTrue; + +import java.awt.geom.Point2D; + +import org.geotools.coverage.grid.GridCoverage2D; +import org.geotools.data.simple.SimpleFeatureCollection; +import org.geotools.feature.FeatureCollections; +import org.geotools.feature.simple.SimpleFeatureBuilder; +import org.geotools.feature.simple.SimpleFeatureTypeBuilder; +import org.geotools.geometry.jts.ReferencedEnvelope; +import org.geotools.referencing.crs.DefaultGeographicCRS; +import org.junit.Test; +import org.opengis.feature.simple.SimpleFeatureType; +import org.opengis.util.ProgressListener; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Geometry; +import com.vividsolutions.jts.geom.GeometryFactory; +import com.vividsolutions.jts.geom.MultiPoint; +import com.vividsolutions.jts.geom.impl.PackedCoordinateSequenceFactory; + +/** + * @author Martin Davis - OpenGeo + * + */ +public class HeatmapProcessTest { + + /** + * A test of a simple surface, validating that the process + * can be invoked and return a reasonable result in a simple situation. + * + * @throws Exception + */ + @Test + public void testSimpleSurface() { + + ReferencedEnvelope bounds = new ReferencedEnvelope(0, 10, 0, 10, DefaultGeographicCRS.WGS84); + Coordinate[] data = new Coordinate[] { + new Coordinate(4, 4), + new Coordinate(4, 6) + }; + SimpleFeatureCollection fc = createPoints(data, bounds); + + ProgressListener monitor = null; + + HeatmapProcess process = new HeatmapProcess(); + GridCoverage2D cov = process.execute(fc, // data + 20, //radius + null, // weightAttr + 1, // pixelsPerCell + 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); + + } + + private float coverageValue(GridCoverage2D cov, double x, double y) + { + float[] covVal = new float[1]; + Point2D worldPos = new Point2D.Double(x, y); + cov.evaluate(worldPos, covVal); + return covVal[0]; + } + + private SimpleFeatureCollection createPoints(Coordinate[] pts, ReferencedEnvelope bounds) + { + + SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder(); + tb.setName("data"); + tb.setCRS(bounds.getCoordinateReferenceSystem()); + tb.add("shape", MultiPoint.class); + tb.add("value", Double.class); + + SimpleFeatureType type = tb.buildFeatureType(); + SimpleFeatureBuilder fb = new SimpleFeatureBuilder(type); + SimpleFeatureCollection fc = FeatureCollections.newCollection(); + + GeometryFactory factory = new GeometryFactory(new PackedCoordinateSequenceFactory()); + + for (Coordinate p : pts) { + Geometry point = factory.createPoint(p); + fb.add(point); + fb.add(p.z); + fc.add(fb.buildFeature(null)); + } + + return fc; + } + +} |
