--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/geom/VectorUtils.java	Mon Oct 31 16:04:59 2011 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/geom/VectorUtils.java	Mon Oct 31 17:05:14 2011 +0000
@@ -78,5 +78,72 @@
             && Math.abs(Y(a)-Y(b)) < EPSILON;
     }
 
+    public static final Point2D intersection(
+        Point2D p1, Point2D p2,
+        Point2D p3, Point2D p4
+    ) {
+        double x1 = X(p1);
+        double y1 = Y(p1);
+        double x2 = X(p2);
+        double y2 = Y(p2);
+        double x3 = X(p3);
+        double y3 = Y(p3);
+        double x4 = X(p4);
+        double y4 = Y(p4);
+
+        // Compute a1, b1, c1, where line joining points 1 and 2
+        // is "a1 x + b1 y + c1 = 0".
+        double a1 = y2 - y1;
+        double b1 = x1 - x2;
+        double c1 = x2*y1 - x1*y2;
+
+        // Compute r3 and r4.
+        double r3 = a1*x3 + b1*y3 + c1;
+        double r4 = a1*x4 + b1*y4 + c1;
+
+        if (r3 != 0d && r4 != 0d && r3*r4 >= 0) {
+            return null;
+        }
+
+        // Compute a2, b2, c2
+        double a2 = y4 - y3;
+        double b2 = x3 - x4;
+        double c2 = (x4 * y3) - (x3 * y4);
+
+        // Compute r1 and r2
+        double r1 = a2*x1 + b2*y1 + c2;
+        double r2 = a2*x2 + b2*y2 + c2;
+
+        if (r1 != 0d && r2 != 0d && r1*r2 >= 0) {
+            return null;
+        }
+
+        // Line segments intersect: compute intersection point.
+        double denom = a1*b2 - a2*b1;
+
+        if (denom == 0d) { // collinear
+            return null;
+        }
+
+        double offset = Math.abs(denom)/2d;
+
+        // The denom/2 is to get rounding instead of truncating. It
+        // is added or subtracted to the numerator, depending upon the
+        // sign of the numerator.
+        double num = b1*c2 - b2*c1;
+
+        double x = num < 0d
+            ? (num - offset)/denom
+            : (num + offset)/denom;
+
+        num = a2*c1 - a1*c2;
+
+        double y = num < 0d
+            ? (num - offset)/denom
+            : (num + offset)/denom;
+
+        return new Point2D.Double(x, y);
+    }
+
 }
 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :