view gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java @ 1085:ec512e7992c6

Improved data initialization and axis range calculation in vertical profiles (scalar and vector). gnv-artifacts/trunk@1187 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Ingo Weinzierl <ingo.weinzierl@intevation.de>
date Wed, 09 Jun 2010 17:11:54 +0000
parents 22c18083225e
children f953c9a559d8
line wrap: on
line source
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 :

http://dive4elements.wald.intevation.org