sascha@1794: package de.intevation.flys.geom; sascha@1794: sascha@1794: import java.io.Serializable; sascha@1794: sascha@1794: import java.awt.Shape; sascha@1794: sascha@1794: import java.awt.geom.Path2D; sascha@1794: import java.awt.geom.Point2D; sascha@1794: sascha@1794: import java.util.ArrayList; sascha@1794: import java.util.List; sascha@1794: import java.util.Arrays; sascha@1794: import java.util.Comparator; sascha@1794: import java.util.Collections; sascha@1794: sascha@1798: import de.intevation.flys.artifacts.math.Linear; sascha@1798: sascha@1799: import static de.intevation.flys.geom.VectorUtils.X; sascha@1799: import static de.intevation.flys.geom.VectorUtils.Y; sascha@1799: import static de.intevation.flys.geom.VectorUtils.EPSILON; sascha@1799: sascha@1794: public class Polygon2D sascha@1794: implements Serializable sascha@1794: { sascha@1794: public static final Comparator POINT_X_CMP = sascha@1794: new Comparator() { sascha@1794: public int compare(Point2D a, Point2D b) { sascha@1795: double d = X(a) - X(b); sascha@1794: if (d < 0d) return -1; sascha@1794: return d > 0d ? + 1 : 0; sascha@1794: } sascha@1794: }; sascha@1794: sascha@1794: public static final Comparator FIRST_POINT_X = sascha@1794: new Comparator() { sascha@1794: public int compare(Point2D [] a, Point2D [] b) { sascha@1794: if (a.length == 0) return -1; // should not happen. sascha@1794: if (b.length == 0) return +1; // should not happen. sascha@1795: double d = X(a[0]) - Y(b[0]); sascha@1794: if (d < 0d) return -1; sascha@1794: return d > 0d ? + 1 : 0; sascha@1794: } sascha@1794: }; sascha@1794: sascha@1795: protected List points; sascha@1794: sascha@1794: public Polygon2D() { sascha@1795: points = new ArrayList(); sascha@1795: } sascha@1795: sascha@1800: public Polygon2D(List points) { sascha@1799: this.points = points; sascha@1799: } sascha@1799: sascha@1794: public void add(double x, double y) { sascha@1795: points.add(new Point2D.Double(x, y)); sascha@1795: } sascha@1795: sascha@1798: protected static boolean addCheck(Point2D p, List points) { sascha@1795: switch (points.size()) { sascha@1795: case 0: sascha@1795: points.add(p); sascha@1795: return true; sascha@1795: case 1: sascha@1797: if (VectorUtils.epsilonEquals(points.get(0), p)) { sascha@1795: return false; sascha@1795: } sascha@1795: points.add(p); sascha@1795: return true; sascha@1795: default: sascha@1795: int L = points.size()-1; sascha@1795: Point2D last = points.get(L); sascha@1797: if (VectorUtils.epsilonEquals(last, p)) { sascha@1795: return false; sascha@1795: } sascha@1795: Point2D before = points.get(L-1); sascha@1797: if (VectorUtils.collinear(before, last, p)) { sascha@1795: points.set(L, p); sascha@1795: } sascha@1795: else { sascha@1795: points.add(p); sascha@1795: } sascha@1795: return true; sascha@1795: } sascha@1795: } sascha@1795: sascha@1798: public boolean addCheck(Point2D p) { sascha@1798: return addCheck(p, points); sascha@1798: } sascha@1798: sascha@1795: public void addReversed(List other) { sascha@1795: for (int i = other.size()-1; i >= 0; --i) { sascha@1795: addCheck(other.get(i)); sascha@1795: } sascha@1795: } sascha@1795: sascha@1794: public double area() { sascha@1794: double area = 0d; sascha@1794: sascha@1795: for (int i = 0, N = points.size(); i < N; ++i) { sascha@1794: int j = (i + 1) % N; sascha@1795: Point2D pi = points.get(i); sascha@1795: Point2D pj = points.get(j); sascha@1795: area += X(pi)*Y(pj); sascha@1795: area -= X(pj)*Y(pi); sascha@1794: } sascha@1794: sascha@1794: return 0.5d*area; sascha@1794: } sascha@1794: sascha@1794: public Shape toShape() { sascha@1794: Path2D.Double path = new Path2D.Double(); sascha@1794: sascha@1795: int N = points.size(); sascha@1794: sascha@1794: if (N > 0) { sascha@1795: Point2D p = points.get(0); sascha@1795: path.moveTo(X(p), Y(p)); sascha@1794: for (int i = 1; i < N; ++i) { sascha@1795: p = points.get(i); sascha@1795: path.lineTo(X(p), Y(p)); sascha@1794: } sascha@1794: path.closePath(); sascha@1794: } sascha@1794: sascha@1794: return path; sascha@1794: } sascha@1794: sascha@1794: protected static List splitByNaNs( sascha@1794: double [] xs, sascha@1794: double [] ys sascha@1794: ) { sascha@1794: List parts = new ArrayList(); sascha@1794: sascha@1794: List tmp = new ArrayList(xs.length); sascha@1794: sascha@1794: for (int i = 0; i < xs.length; ++i) { sascha@1794: double x = xs[i]; sascha@1794: double y = ys[i]; sascha@1794: sascha@1794: if (Double.isNaN(x) || Double.isNaN(y)) { sascha@1794: if (!tmp.isEmpty()) { sascha@1794: Point2D [] part = new Point2D[tmp.size()]; sascha@1794: parts.add(tmp.toArray(part)); sascha@1794: tmp.clear(); sascha@1794: } sascha@1794: } sascha@1794: else { sascha@1794: tmp.add(new Point2D.Double(x, y)); sascha@1794: } sascha@1794: } sascha@1794: sascha@1794: if (!tmp.isEmpty()) { sascha@1794: Point2D [] part = new Point2D[tmp.size()]; sascha@1794: parts.add(tmp.toArray(part)); sascha@1794: } sascha@1794: sascha@1794: return parts; sascha@1794: } sascha@1794: sascha@1795: protected static boolean removeNoneIntersecting( sascha@1795: List As, sascha@1795: List Bs sascha@1795: ) { sascha@1798: int B = Bs.size()-1; sascha@1795: OUTER: for (int i = 0; i < As.size();) { sascha@1795: Point2D [] a = As.get(i); sascha@1798: int lo = 0, hi = B; sascha@1795: while (lo <= hi) { sascha@1795: int mid = (lo + hi) >> 1; sascha@1795: Point2D [] b = Bs.get(mid); sascha@1795: if (X(a[0]) > X(b[b.length-1])) lo = mid+1; sascha@1795: else if (X(a[a.length-1]) < X(b[0])) hi = mid-1; sascha@1795: else { sascha@1795: // found: keep sascha@1795: ++i; sascha@1795: continue OUTER; sascha@1795: } sascha@1795: } sascha@1795: // not found: remove sascha@1795: As.remove(i); sascha@1795: } sascha@1795: sascha@1795: return As.isEmpty(); sascha@1795: } sascha@1795: sascha@1795: protected static void buildPolygons( sascha@1795: Point2D [] as, sascha@1795: Point2D [] bs, sascha@1795: Point2D [][] rest, sascha@1795: List positives, sascha@1795: List negatives sascha@1795: ) { sascha@1798: List apoints = new ArrayList(); sascha@1798: List bpoints = new ArrayList(); sascha@1798: sascha@1798: double ax = X(as[0]); sascha@1798: double bx = X(bs[0]); sascha@1798: sascha@1799: int ai = 1; sascha@1799: int bi = 1; sascha@1799: sascha@1799: boolean intersect = false; sascha@1798: sascha@1798: if (ax == bx) { sascha@1798: apoints.add(as[0]); sascha@1798: bpoints.add(bs[0]); sascha@1798: } sascha@1798: else if (ax > bx) { sascha@1798: apoints.add(as[0]); sascha@1798: double bx1; sascha@1798: while ((bx1 = X(bs[bi])) < ax) ++bi; sascha@1798: if (bx1 == ax) { sascha@1798: bpoints.add(bs[bi]); sascha@1798: } sascha@1798: else { // need to calculate start b point. sascha@1799: intersect = true; sascha@1798: double by1 = Linear.linear( sascha@1798: ax, sascha@1798: X(bs[bi-1]), bx1, sascha@1798: Y(bs[bi-1]), Y(bs[bi])); sascha@1798: sascha@1798: bpoints.add(new Point2D.Double(ax, by1)); sascha@1798: } sascha@1798: } sascha@1798: else { // bx > ax: Symmetric case sascha@1798: bpoints.add(bs[0]); sascha@1798: double ax1; sascha@1798: while ((ax1 = X(as[ai])) < bx) ++ai; sascha@1798: if (ax1 == bx) { sascha@1798: apoints.add(as[ai]); sascha@1798: } sascha@1798: else { // need to calculate start b point. sascha@1799: intersect = true; sascha@1798: double ay1 = Linear.linear( sascha@1798: bx, sascha@1798: X(as[ai-1]), ax1, sascha@1798: Y(as[ai-1]), Y(as[ai])); sascha@1798: sascha@1798: apoints.add(new Point2D.Double(bx, ay1)); sascha@1798: } sascha@1798: } sascha@1798: sascha@1798: // now we have a point in each list, decide if neg/pos. sascha@1798: boolean neg = Y(bpoints.get(0)) > Y(apoints.get(0)); sascha@1798: sascha@1799: // Continue with inner points sascha@1799: sascha@1799: Line line = new Line(); sascha@1799: sascha@1799: while (ai < as.length && bi < bs.length) { sascha@1799: double xan = X(as[ai]); sascha@1799: double xbn = X(bs[bi]); sascha@1799: if (xan == xbn) { sascha@1799: double yan = Y(as[ai]); sascha@1799: double ybn = Y(bs[ai]); sascha@1799: sascha@1799: if (neg) { sascha@1799: if (yan > ybn) { // intersection sascha@1801: Point2D ip = VectorUtils.intersection( sascha@1801: apoints.get(apoints.size()-1), as[ai], sascha@1801: bpoints.get(bpoints.size()-1), bs[bi]); sascha@1801: addCheck(ip, apoints); sascha@1801: addCheck(ip, bpoints); sascha@1801: Polygon2D p = new Polygon2D( sascha@1801: new ArrayList(apoints)); sascha@1801: p.addReversed(bpoints); sascha@1801: negatives.add(p); sascha@1801: apoints.clear(); sascha@1801: bpoints.clear(); sascha@1801: apoints.add(ip); sascha@1801: bpoints.add(ip); sascha@1801: neg = !neg; sascha@1799: } sascha@1799: else { // no intersection sascha@1799: addCheck(as[ai], apoints); sascha@1799: addCheck(bs[bi], bpoints); sascha@1799: } sascha@1799: } sascha@1799: else { // not neg sascha@1799: if (ybn > yan) { // intersection sascha@1801: Point2D ip = VectorUtils.intersection( sascha@1801: apoints.get(apoints.size()-1), as[ai], sascha@1801: bpoints.get(bpoints.size()-1), bs[bi]); sascha@1801: addCheck(ip, apoints); sascha@1801: addCheck(ip, bpoints); sascha@1801: Polygon2D p = new Polygon2D( sascha@1801: new ArrayList(apoints)); sascha@1801: p.addReversed(bpoints); sascha@1801: positives.add(p); sascha@1801: apoints.clear(); sascha@1801: bpoints.clear(); sascha@1801: apoints.add(ip); sascha@1801: bpoints.add(ip); sascha@1801: neg = !neg; sascha@1799: } sascha@1799: else { // no intersection sascha@1799: addCheck(as[ai], apoints); sascha@1799: addCheck(bs[bi], bpoints); sascha@1799: } sascha@1799: } sascha@1801: ++ai; sascha@1801: ++bi; sascha@1799: } sascha@1799: else if (xan > xbn) { sascha@1799: line.set(apoints.get(apoints.size()-1), as[ai]); sascha@1799: double dir = neg ? -1d: 1d; // XXX: correct sign? sascha@1799: while (bi < bs.length sascha@1799: && X(bs[bi]) < xan sascha@1799: && line.eval(bs[bi])*dir > EPSILON) sascha@1799: ++bi; sascha@1799: if (bi == bs.length) { sascha@1799: // b run out of points sascha@1799: // calculate Y(last_a, as[ai]) for X(bs[bi-1]) sascha@1799: double ay1 = Linear.linear( sascha@1799: X(bs[bi-1]), sascha@1799: X(apoints.get(apoints.size()-1)), X(as[ai]), sascha@1799: Y(apoints.get(apoints.size()-1)), Y(as[ai])); sascha@1799: addCheck(new Point2D.Double(X(bs[bi-1]), ay1), apoints); sascha@1799: addCheck(bs[bi-1], bpoints); sascha@1801: Polygon2D p = new Polygon2D( sascha@1801: new ArrayList(apoints)); sascha@1801: p.addReversed(bpoints); sascha@1801: apoints.clear(); sascha@1801: bpoints.clear(); sascha@1801: (neg ? negatives : positives).add(p); sascha@1799: break; sascha@1799: } sascha@1799: else { sascha@1799: // TODO: intersect line and/or X(bs[bi]) >= xan? sascha@1799: } sascha@1799: } sascha@1799: else { // xbn > xan sascha@1799: line.set(bpoints.get(bpoints.size()-1), bs[bi]); sascha@1799: // TODO: continue symmetric sascha@1799: } sascha@1799: } sascha@1799: sascha@1799: // TODO: Continue with closing segment sascha@1799: } sascha@1799: sascha@1799: public static final class Line { sascha@1799: sascha@1799: private double a; sascha@1799: private double b; sascha@1799: private double c; sascha@1799: sascha@1799: public Line() { sascha@1799: } sascha@1799: sascha@1799: public Line(Point2D p1, Point2D p2) { sascha@1799: set(p1, p2); sascha@1799: } sascha@1799: sascha@1799: public void set(Point2D p1, Point2D p2) { sascha@1799: Point2D p3 = sascha@1799: VectorUtils.normalize( sascha@1799: VectorUtils.sub(p1, p2)); sascha@1799: sascha@1799: Point2D n = VectorUtils.ortho(p3); sascha@1799: sascha@1799: a = X(n); sascha@1799: b = Y(n); sascha@1799: sascha@1799: // a*x + b*y + c = 0 sascha@1799: // c = -a*x -b*y sascha@1799: sascha@1799: c = -a*X(p1) - b*Y(p1); sascha@1799: } sascha@1799: sascha@1799: public double eval(Point2D p) { sascha@1799: return a*X(p) + b*Y(p) + c; sascha@1799: } sascha@1795: } sascha@1795: sascha@1794: public static void createPolygons( sascha@1794: double [] xAs, double [] yAs, sascha@1794: double [] xBs, double [] yBs, sascha@1794: List positives, sascha@1794: List negatives sascha@1794: ) { sascha@1794: if (xAs.length == 0 || xBs.length == 0) { sascha@1794: return; sascha@1794: } sascha@1794: sascha@1794: List splAs = splitByNaNs(xAs, yAs); sascha@1794: List splBs = splitByNaNs(xBs, yBs); sascha@1794: sascha@1794: // They feeded us with NaNs only. sascha@1794: if (splAs.isEmpty() || splBs.isEmpty()) { sascha@1794: return; sascha@1794: } sascha@1794: sascha@1794: // Sort each part by x to ensure that the first sascha@1794: // is the smallest. sascha@1794: for (Point2D [] splA: splAs) { sascha@1794: Arrays.sort(splA, POINT_X_CMP); sascha@1794: } sascha@1794: sascha@1794: for (Point2D [] splB: splBs) { sascha@1794: Arrays.sort(splB, POINT_X_CMP); sascha@1794: } sascha@1794: sascha@1794: // Now sort all parts by there first elements. sascha@1794: // Should be good enough to find overlapping regions. sascha@1794: Collections.sort(splAs, FIRST_POINT_X); sascha@1794: Collections.sort(splBs, FIRST_POINT_X); sascha@1794: sascha@1794: // Check if the two series intersect at all. sascha@1794: // If no then there will be no area between them. sascha@1794: sascha@1794: Point2D [] p1 = splAs.get(0); sascha@1794: Point2D [] p2 = splBs.get(splBs.size()-1); sascha@1794: sascha@1795: // Sort out the ranges that are not intersecting sascha@1795: // the ranges in the other series. sascha@1795: // We are going to merge them anyway sascha@1795: // so this is not strictly required. sascha@1795: // Keep it to recude cases. sascha@1795: if (removeNoneIntersecting(splAs, splBs) sascha@1795: || removeNoneIntersecting(splBs, splAs) sascha@1795: ) { sascha@1795: // They do not intersect at all. sascha@1794: return; sascha@1794: } sascha@1794: sascha@1794: // TODO: Intersect/split the two series parts. sascha@1794: } sascha@1794: } sascha@1794: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :