--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Fri Sep 28 12:14:00 2012 +0200
@@ -0,0 +1,416 @@
+/*
+ * Copyright (c) 2010 by Intevation GmbH
+ *
+ * This program is free software under the LGPL (>=v2.1)
+ * Read the file LGPL.txt coming with the software for details
+ * or visit http://www.gnu.org/licenses/ if it does not exist.
+ */
+
+package de.intevation.gnv.math;
+
+import com.vividsolutions.jts.algorithm.CGAlgorithms;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+import com.vividsolutions.jts.geom.Geometry;
+import com.vividsolutions.jts.geom.GeometryFactory;
+import com.vividsolutions.jts.geom.LinearRing;
+import com.vividsolutions.jts.geom.Point;
+import com.vividsolutions.jts.geom.Polygon;
+
+import com.vividsolutions.jts.index.ItemVisitor;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+
+import org.apache.log4j.Logger;
+
+/**
+ * Spans a rectangle of points to be used in spatial indexing.
+ *
+ * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
+ */
+public class GridCell
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(GridCell.class);
+
+    /**
+     * Callback for {@link com.vividsolutions.jts.index.SpatialIndex}
+     * to find a GridCell which contains a given point.
+     */
+    public static final class CellFinder
+    implements                ItemVisitor
+    {
+        /**
+         * Stores the found GridCell.
+         */
+        public    GridCell found;
+
+        /**
+         * The query point.
+         */
+        protected Point    point;
+
+        /**
+         * Default constructor.
+         */
+        public CellFinder() {
+        }
+
+        /**
+         * Prepares a spatial index lookup.
+         * @param center The query point.
+         */
+        public void prepare(Coordinate center) {
+            found = null;
+            point = GEOMETRY_FACTORY.createPoint(center);
+        }
+
+        /**
+         * Called by the spatial index as a 2nd order filter.
+         * @param item The GridCell to test.
+         */
+        public void visitItem(Object item) {
+            if (found == null) {
+                GridCell cell = (GridCell)item;
+                if (cell.contains(point)) {
+                    found = cell;
+                }
+            }
+        }
+    } // class CellFinder
+
+    /**
+     * The first point.
+     */
+    public Point2d p1;
+    /**
+     * The second point.
+     */
+    public Point2d p2;
+    /**
+     * The third point.
+     */
+    public Point2d p3;
+    /**
+     * The fourth point.
+     */
+    public Point2d p4;
+
+    /**
+     * Polygon created from the four points.
+     */
+    protected Polygon polygon;
+
+    /**
+     * Geometry factory to create the query points and the polygons.
+     */
+    public static final GeometryFactory GEOMETRY_FACTORY
+        = new GeometryFactory();
+
+    /**
+     * Default constructor.
+     */
+    public GridCell() {
+    }
+
+    /**
+     * Constructor to create a GridCell out of four given points..
+     * @param p1 The first point.
+     * @param p2 The second point.
+     * @param p3 The thrid point.
+     * @param p4 The fourth point.
+     */
+    public GridCell(Point2d p1, Point2d p2, Point2d p3, Point2d p4) {
+        this.p1 = p1;
+        this.p2 = p2;
+        this.p3 = p3;
+        this.p4 = p4;
+        createPolygon();
+    }
+
+    /**
+     * Creates the polygon from the four points.
+     */
+    protected void createPolygon() {
+        Coordinate [] coords = new Coordinate [] { p1, p2, p3, p4, p1 };
+        if (!CGAlgorithms.isCCW(coords)) {
+            for (int i = 0, j = coords.length-1; i < j; ++i, --j) {
+                Coordinate c = coords[i];
+                coords[i] = coords[j];
+                coords[j] = c;
+            }
+        }
+        LinearRing shell = GEOMETRY_FACTORY
+            .createLinearRing(coords);
+
+        if (!shell.isValid()) {
+            log.warn("linear ring is not valid");
+        }
+
+        polygon = GEOMETRY_FACTORY.createPolygon(shell, null);
+    }
+
+    /**
+     * Returns the envelope of the four point polygon.
+     * @return the envelope.
+     */
+    public Envelope getEnvelope() {
+        return polygon.getEnvelopeInternal();
+    }
+
+    /**
+     * Test if a given point is inside this grid cell.
+     * @param coord
+     * @return true, if this grid cell contains the given point - otherwise
+     * false.
+     */
+    public boolean contains(Geometry coord) {
+        return polygon.contains(coord);
+    }
+
+    /**
+     * Converts a list of points to a list of grid cells.
+     * @param points This list of points.
+     * @return This list of grid cells.
+     */
+    public static List<GridCell> pointsToGridCells(
+        List<? extends Point2d> points
+    ) {
+        return pointsToGridCells(points, null);
+    }
+
+    /**
+     * Converts a list of points to a list of grid cells.
+     * Points that are not in a relevant area are ignored.
+     * @param points The list of points.
+     * @param relevantArea The relevant area.
+     * @return The list of grid cells.
+     */
+    public static List<GridCell> pointsToGridCells(
+        List<? extends Point2d> points,
+        Envelope                relevantArea
+    ) {
+        return pointsToGridCells(points, relevantArea, 0);
+    }
+
+    private static final int NEIGHBORS [][] = {
+        { -1, -1 }, // 0
+        { -1,  0 }, // 1
+        { -1, +1 }, // 2
+        {  0, +1 }, // 3
+        { +1, +1 }, // 4
+        { +1,  0 }, // 5
+        { +1, -1 }, // 6
+        {  0, -1 }  // 7
+    };
+
+    /**
+     * Generate points by extrapolating border points.
+     * @param rows (i, j) indexed map of points.
+     * @param minI min known i.
+     * @param maxI max known i.
+     * @param minJ min known j.
+     * @param maxJ max known j.
+     * @param rounds Deternine how many extra rings should be generated.
+     * @param relevantArea The relevant area.
+     * @return number of newly generated points.
+     */
+    public static int extrapolate(
+        HashMap<Integer, HashMap<Integer, Point2d>> rows,
+        int minI, int maxI,
+        int minJ, int maxJ,
+        int rounds,
+        Envelope relevantArea
+    ) {
+        Point2d [] neighbors = new Point2d[NEIGHBORS.length];
+        Point2d [] positions = new Point2d[NEIGHBORS.length];
+
+        int total = 0;
+
+        ArrayList<ArrayList<IJKey>> prio =
+            new ArrayList<ArrayList<IJKey>>(NEIGHBORS.length);
+
+        for (int i = 0; i < NEIGHBORS.length; ++i) {
+            prio.add(new ArrayList<IJKey>());
+        }
+
+        while (rounds-- > 0) {
+            for (int i = minI; i <= maxI; ++i) {
+                for (int j = minJ; j <= maxJ; ++j) {
+                    Point2d p = get(rows, i, j);
+                    if (p != null) {
+                        continue;
+                    }
+
+                    int count = 0;
+
+                    for (int k = 0; k < neighbors.length; ++k) {
+                        neighbors[k] = positions[k] = null;
+                        int dij [] = NEIGHBORS[k];
+                        Point2d n1 = get(rows, i+  dij[0], j+  dij[1]);
+                        Point2d n2 = get(rows, i+2*dij[0], j+2*dij[1]);
+                        if (n1 != null && n2 != null) {
+                            ++count;
+                        }
+                    }
+
+                    if (count > 0) {
+                        prio.get(count-1).add(new IJKey(i, j));
+                    }
+                } // for all columns
+            } // for all rows
+
+            int N = 0;
+
+            for (int l = NEIGHBORS.length-1; l >= 0; --l) {
+                ArrayList<IJKey> list = prio.get(l);
+                for (IJKey ij: list) {
+                    int i = ij.i;
+                    int j = ij.j;
+                    for (int k = 0; k < neighbors.length; ++k) {
+                        neighbors[k] = positions[k] = null;
+                        int dij [] = NEIGHBORS[k];
+                        Point2d n1 = get(rows, i+  dij[0], j+  dij[1]);
+                        Point2d n2 = get(rows, i+2*dij[0], j+2*dij[1]);
+                        if (n1 != null && n2 != null) {
+                            neighbors[k] = n1;
+                            positions[k] = n1.extrapolate(-1d, n2);
+                        }
+                    }
+
+                    Point2d avg = Point2d.average(positions);
+
+                    if (avg != null && avg.near(positions)
+                    && (relevantArea == null
+                        || relevantArea.contains(avg.x, avg.y))) {
+                        avg.i = i;
+                        avg.j = j;
+                        avg.inverseDistanceWeighting(neighbors);
+                        put(rows, avg, i, j);
+                    }
+                }
+                N += list.size();
+                list.clear();
+            }
+
+            if (N == 0) {
+                break;
+            }
+            total += N;
+        } // for all rounds
+
+        return total;
+    }
+
+    /**
+     * Converts a list of points to a list of grid cells.
+     * @param points The list of points.
+     * @param relevantArea The relevant area. If a point is
+     * not inside this area it is ignored during the build process.
+     * @param extrapolationRounds Number of extra point rings.
+     * 0 = no extrpolation.
+     * @return The list of grid cells.
+     */
+    public static List<GridCell> pointsToGridCells(
+        List<? extends Point2d> points,
+        Envelope                relevantArea,
+        int                     extrapolationRounds
+    ) {
+        int minI = Integer.MAX_VALUE;
+        int maxI = Integer.MIN_VALUE;
+        int minJ = Integer.MAX_VALUE;
+        int maxJ = Integer.MIN_VALUE;
+
+        HashMap<Integer, HashMap<Integer, Point2d>> rows =
+            new HashMap<Integer, HashMap<Integer, Point2d>>();
+
+        int culled = 0;
+
+        for (Point2d p: points) {
+
+            if (relevantArea != null && !relevantArea.contains(p.x, p.y)) {
+                ++culled;
+                continue;
+            }
+
+            if (p.i < minI) minI = p.i;
+            if (p.i > maxI) maxI = p.i;
+            if (p.j < minJ) minJ = p.j;
+            if (p.j > maxJ) maxJ = p.j;
+
+            HashMap<Integer, Point2d> row = rows.get(p.i);
+
+            if (row == null) {
+                rows.put(p.i, row = new HashMap<Integer, Point2d>());
+            }
+
+            row.put(p.j, p);
+        }
+
+        ArrayList<GridCell> cells = new ArrayList<GridCell>(points.size());
+
+        int extrapolated = extrapolate(
+            rows,
+            minI, maxI,
+            minJ, maxJ,
+            extrapolationRounds,
+            relevantArea);
+
+        for (int i = minI + 1; i <= maxI; ++i) {
+            HashMap<Integer, Point2d> row1 = rows.get(i-1);
+            HashMap<Integer, Point2d> row2 = rows.get(i);
+
+            if (row1 == null || row2 == null) {
+                continue;
+            }
+
+            for (int j = minJ + 1; j < maxJ; ++j) {
+                Point2d p1 = row1.get(j-1);
+                Point2d p2 = row1.get(j);
+                Point2d p3 = row2.get(j);
+                Point2d p4 = row2.get(j-1);
+
+                if (p1 != null && p2 != null && p3 != null && p4 != null) {
+                    cells.add(new GridCell(p1, p2, p3, p4));
+                }
+            }
+        } // for all rows
+
+        if (log.isDebugEnabled()) {
+            log.debug("culled points: " + culled);
+            log.debug("extrapolated points: " + extrapolated);
+            log.debug("min/max i: " + minI + " / " + maxI);
+            log.debug("min/max j: " + minJ + " / " + maxJ);
+            log.debug("cells found: " + cells.size());
+        }
+
+        return cells;
+    }
+
+    private static Point2d get(
+        HashMap<Integer, HashMap<Integer, Point2d>> rows,
+        int i, int j
+    ) {
+        HashMap<Integer, Point2d> row = rows.get(i);
+        return row != null ? row.get(j) : null;
+    }
+
+    private static void put(
+        HashMap<Integer, HashMap<Integer, Point2d>> rows,
+        Point2d point,
+        int i, int j
+    ) {
+        Integer I = Integer.valueOf(i);
+        HashMap<Integer, Point2d> row = rows.get(I);
+        if (row == null) {
+            rows.put(I, row = new HashMap<Integer, Point2d>());
+        }
+        row.put(j, point);
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :