--- a/doodle/tk/geometry.d	Sun Aug 30 01:33:34 2009 +0930
+++ b/doodle/tk/geometry.d	Sun Aug 30 01:34:14 2009 +0930
@@ -6,6 +6,18 @@
     import doodle.tk.misc;
 }
 
+// In doodle x and y increase right/east and up/north respectively.
+
+// TODO explain the strategy for ensuring numerical stability
+// and the division of responsibility between users of these
+// types and the types themselves.
+//
+// Explain how numerical instability is handled. The current policy
+// is to correct bad user input (eg a gradient with miniscule length) and
+// print warnings rather than have assertions and cause crashes.
+
+// A location in 2D space
+
 struct Point {
     static immutable Point DEFAULT;
 
@@ -30,7 +42,7 @@
         return Vector(_x - p._x, _y - p._y);
     }
 
-    string toString() {
+    string toString() const {
         return std.string.format("(%f, %f)", _x, _y);
     }
 
@@ -50,6 +62,8 @@
     return Point(max(a.x, b.x), max(a.y, b.y));
 }
 
+// The displacement between two locations in 2D space
+
 struct Vector {
     static Vector DEFAULT;
 
@@ -86,7 +100,7 @@
         return sqrt(_x * _x + _y * _y);
     }
 
-    string toString() {
+    string toString() const {
         return std.string.format("[%f, %f]", _x, _y);
     }
 
@@ -98,6 +112,11 @@
     }
 }
 
+// A rectangle in 2D space.
+// Internally represented by:
+//   a point defining the bottom left corner
+//   a vector defining the displacement to the upper right corner
+
 struct Rectangle {
     static Rectangle DEFAULT;
 
@@ -118,31 +137,18 @@
         this(corner1.x, corner1.y, corner.x - corner1.x, corner.y - corner1.y);
     }
 
-    Point position() const {
-        return _position;
-    }
+    alias position min_corner;
+    Point position() const { return _position; }
 
-    Vector size() const {
-        return _size;
-    }
-
-    alias position min_corner;
+    Vector size() const { return _size; }
 
-    Point max_corner() const {
-        return _position + _size;
-    }
+    Point max_corner() const { return _position + _size; }
 
-    bool valid() const {
-        return _size.x > 0.0 & _size.y > 0.0;
-    }
+    bool valid() const { return _size.x > 0.0 & _size.y > 0.0; }
 
-    bool invalid() const {
-        return !valid();
-    }
+    bool invalid() const { return !valid(); }
 
-    double area() const {
-        return _size.x * _size.y;
-    }
+    double area() const { return _size.x * _size.y; }
 
     // Intersection
     Rectangle opAnd(in Rectangle r) const {
@@ -172,37 +178,45 @@
         }
         else {
             return Rectangle(min_extents(min_corner(), r.min_corner()),
-                              max_extents(max_corner(), r.max_corner()));
+                             max_extents(max_corner(), r.max_corner()));
         }
     }
 
-    // TODO make these free functions
+    // TODO review these functions.
+    // Want a clear and simple set.
+    // Consider making them free functions
 
     Rectangle moved(in Vector displacement) const {
         return Rectangle(_position + displacement, _size);
     }
 
+    Rectangle repositioned(in Point new_position) const {
+        return Rectangle(new_position, _size);
+    }
+
+    // Operations about the bottom left corner
+
     Rectangle expanded(in Vector expand_amount) const {
         return Rectangle(_position, _size + expand_amount);
     }
 
+    Rectangle shrunk(in Vector shrink_amount) const {
+        return Rectangle(_position, _size - shrink_amount);
+    }
+
+    // Operations about the centre
+
     Rectangle feathered(double amount) const {
         assert(amount >= 0.0);
         return Rectangle(Point(_position.x - amount, _position.y - amount),
                          Vector(_size.x + 2.0 * amount, _size.y + 2.0 * amount));
     }
 
-    Rectangle shrunk(in Vector shrink_amount) const {
-        return Rectangle(_position, _size - shrink_amount);
-    }
-
     Rectangle resized(in Vector new_size) const {
         return Rectangle(_position, new_size);
     }
 
-    Rectangle repositioned(in Point new_position) const {
-        return Rectangle(new_position, _size);
-    }
+    //
 
     void get_quantised(out int x, out int y, out int w, out int h) const {
         x = cast(int)floor(_position.x);
@@ -213,11 +227,9 @@
 
     //
 
-    Point centre() const {
-        return _position + _size / 2.0;
-    }
+    Point centre() const { return _position + _size / 2.0; }
 
-    string toString() {
+    string toString() const {
         return std.string.format("{%s, %s}", _position, _size);
     }
 
@@ -234,17 +246,129 @@
     }
 }
 
-/*
+private {
+    // This is a commmon building block for intersection of lines and segments
+    // Two lines "a" and "b" are defined by two points each 
+    // "ua" and "ub" define the intersection as a fraction alone
+    // the two respetive line segments (FIXME this comment sucks)
+    // Influenced by http://ozviz.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
+    bool compute_intersection(in Point pa1, in Point pa2, out double ua,
+                              in Point pb1, in Point pb2, out double ub) {
+        //writefln("pa1: %s, pa2: %s, pb1: %s, pb2: %s", pa1, pa2, pb1, pb2);
+        double den = (pb2.y - pb1.y) * (pa2.x - pa1.x) - (pb2.x - pb1.x) * (pa2.y - pa1.y);
+
+        if (abs(den) < 1e-9) {          // TODO consolidate constants used for numerical stability
+            // Lines are parallel or nearly so
+            writefln("Warning, parallel lines!");
+            return false;
+        }
+        else {
+            // It will be safe to divide by den
+            double num_a = (pb2.x - pb1.x) * (pa1.y - pb1.y) - (pb2.y - pb1.y) * (pa1.x - pb1.x);
+            double num_b = (pa2.x - pa1.x) * (pa1.y - pb1.y) - (pa2.y - pa1.y) * (pa1.x - pb1.x);
+
+            ua = num_a / den;
+            ub = num_b / den;
+
+            {
+                // FIXME remove this debugging stuff
+                //writefln("ua: %s, ub: %s", ua, ub);
+                Point pa = pa1 + ua * (pa2 - pa1);
+                Point pb = pb1 + ub * (pb2 - pb1);
+                //writefln("pa: %s, pb: %s", pa, pb);
+                assert(pa1 + ua * (pa2 - pa1) == pb1 + ub * (pb2 - pb1));
+            }
+
+            return true;
+        }
+    }
+}
+
+// A line (notionally infinitely extending in both directions) in 2D space.
+// Internally represented by:
+//   a point at an arbitrary location along the line
+//   a vector defining the gradient of the line
+
+struct Line {
+    this(in Point p, in Vector g) {
+        _point = p;
+        _gradient = g;
+        assert(_gradient.length > 1e-6);        // FIXME how to best deal with this
+    }
+
+    this(in Point a, in Point b) {
+        _point = a;
+        _gradient = b - a;
+        assert(_gradient.length > 1e-6);        // FIXME as above
+    }
+
+    Point point() const { return _point; }
+    Vector gradient() const { return _gradient; }
+
+    private {
+        Point _point;           // Arbitrary point along line
+        Vector _gradient;
+    }
+}
+
+// Calculate the point "p" where lines "a" and "b" intersect.
+// Returns false if lines are parallel or too close for numerical stability
+bool intersection(in Line a, in Line b, out Point p) {
+    Point  pa = a.point;
+    Vector va = a.gradient;
+    double ua;
+    Point  pb = b.point;
+    Vector vb = b.gradient;
+    double ub;
+
+    if (compute_intersection(pa, pa + va, ua, pb, pb + vb, ub)) {
+        // We could just have easily evaluated for line b...
+        p = pa + ua * va;
+        // p = pb + ub * vb;
+        return true;
+    }
+    else {
+        return false;
+    }
+}
+
+// A line segment (has a beginning and an end) in 2D space.
+
 struct Segment {
+    this(in Point a, in Point b) {
+        _begin = a;
+        _end = b;
+    }
+
+    Point begin() const { return _begin; }
+    Point end() const { return _end; }
+
     private {
         Point _begin, _end;
     }
 }
 
-struct Line {
-    private {
-        Point _point;
-        Vector _gradien;
+bool intersection(in Segment a, in Segment b, out Point p) {
+    Point pa1 = a.begin;
+    Point pa2 = a.end;
+    double ua;
+    Point pb1 = b.begin;
+    Point pb2 = b.end;
+    double ub;
+
+    if (compute_intersection(pa1, pa2, ua, pb1, pb2, ub)) {
+        if (ua >= 0.0 && ua <= 1.0 &&     // inside of segment a
+            ub >= 0.0 && ub <= 1.0) {     // inside of segment b
+            // We could just have easily evaluated for line b...
+            p = pa1 + ua * (pa2 - pa1);
+            // p = pa2 + ub * (pb2 - pb1);
+            return true;
+        }
+        else {
+            return false;
+        }
+    }
+    else {
+        return false;
     }
 }
-*/