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