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