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<Point2D> POINT_X_CMP =
sascha@1794:         new Comparator<Point2D>() {
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<Point2D []> FIRST_POINT_X =
sascha@1794:         new Comparator<Point2D []>() {
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<Point2D> points;
sascha@1794: 
sascha@1794:     public Polygon2D() {
sascha@1795:         points = new ArrayList<Point2D>();
sascha@1795:     }
sascha@1795: 
sascha@1800:     public Polygon2D(List<Point2D> 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<Point2D> 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<Point2D> 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<Point2D []> splitByNaNs(
sascha@1794:         double [] xs,
sascha@1794:         double [] ys
sascha@1794:     ) {
sascha@1794:         List<Point2D []> parts = new ArrayList<Point2D []>();
sascha@1794: 
sascha@1794:         List<Point2D> tmp = new ArrayList<Point2D>(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<Point2D []> As,
sascha@1795:         List<Point2D []> 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<Polygon2D> positives,
sascha@1795:         List<Polygon2D> negatives
sascha@1795:     ) {
sascha@1798:         List<Point2D> apoints = new ArrayList<Point2D>();
sascha@1798:         List<Point2D> bpoints = new ArrayList<Point2D>();
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<Point2D>(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<Point2D>(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<Point2D>(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<Polygon2D> positives,
sascha@1794:         List<Polygon2D> negatives
sascha@1794:     ) {
sascha@1794:         if (xAs.length == 0 || xBs.length == 0) {
sascha@1794:             return;
sascha@1794:         }
sascha@1794: 
sascha@1794:         List<Point2D []> splAs = splitByNaNs(xAs, yAs);
sascha@1794:         List<Point2D []> 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 :