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 Sascha L. Teichmann 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 pointsToGridCells( sascha@514: List 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 pointsToGridCells( sascha@517: List 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> 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> prio = sascha@593: new ArrayList>(NEIGHBORS.length); sascha@593: sascha@593: for (int i = 0; i < NEIGHBORS.length; ++i) { sascha@593: prio.add(new ArrayList()); 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 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 pointsToGridCells( sascha@593: List 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> rows = sascha@514: new HashMap>(); 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 row = rows.get(p.i); sascha@514: sascha@514: if (row == null) { sascha@514: rows.put(p.i, row = new HashMap()); sascha@514: } sascha@514: sascha@514: row.put(p.j, p); sascha@514: } sascha@514: sascha@514: ArrayList cells = new ArrayList(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 row1 = rows.get(i-1); sascha@514: HashMap 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> rows, sascha@593: int i, int j sascha@593: ) { sascha@593: HashMap 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> rows, sascha@593: Point2d point, sascha@593: int i, int j sascha@593: ) { sascha@593: Integer I = Integer.valueOf(i); sascha@593: HashMap row = rows.get(I); sascha@593: if (row == null) { sascha@593: rows.put(I, row = new HashMap()); 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 :