--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Fri Jan 29 13:28:28 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Sun Jan 31 21:50:15 2010 +0000
@@ -1,5 +1,7 @@
 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;
@@ -73,8 +75,21 @@
     }
 
     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(new Coordinate [] { p1, p2, p3, p4, p1 });
+            .createLinearRing(coords);
+
+        if (!shell.isValid()) {
+            log.warn("linear ring is not valid");
+        }
+
         polygon = GEOMETRY_FACTORY.createPolygon(shell, null);
     }
 
@@ -96,6 +111,111 @@
         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
+    };
+
+    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;
+    }
+
+    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;
@@ -129,25 +249,36 @@
 
         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) {
-                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 columns in row
+            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());
@@ -155,5 +286,26 @@
 
         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 :