Added Robust Line intersector. Need more tests still

algorithm/cga/big-double-cga.go → algorithm/big_cga/big_cga.go

@@ -1,4 +1,4 @@
-package cga
+package big_cga
 
 import (
 	"github.com/twpayne/go-geom"
@@ -11,38 +11,54 @@ import (
  */
 var dp_safe_epsilon = 1e-15
 
+type Orientation int
+
+const (
+	CLOCKWISE Orientation = iota - 1
+	COLLINEAR
+	COUNTER_CLOCKWISE
+)
+
+var orientationLabels = [3]string{"CLOCKWISE", "COLLINEAR", "COUNTER_CLOCKWISE"}
+
+func (o Orientation) String() string {
+	return orientationLabels[int(o+1)]
+}
+
 /**
- * Returns the index of the direction of the point <code>q</code> relative to
- * a vector specified by <code>p1-p2</code>.
+ * Returns the index of the direction of the point <code>point</code> relative to
+ * a vector specified by <code>vectorOrigin-vectorEnd</code>.
  *
- * @param p1 the origin point of the vector
- * @param p2 the final point of the vector
- * @param q the point to compute the direction to
+ * @param vectorOrigin the origin point of the vector
+ * @param vectorEnd the final point of the vector
+ * @param point the point to compute the direction to
  *
- * @return 1 if q is counter-clockwise (left) from p1-p2
- * @return -1 if q is clockwise (right) from p1-p2
- * @return 0 if q is collinear with p1-p2
+ * @return COUNTER_CLOCKWISE if point is counter-clockwise (left) from vectorOrigin-vectorEnd
+ * @return CLOCKWISE if point is clockwise (right) from vectorOrigin-vectorEnd
+ * @return COLLINEAR if point is collinear with vectorOrigin-vectorEnd
  */
-func OrientationIndex(p1, p2, q geom.Coord) int {
+func OrientationIndex(vectorOrigin, vectorEnd, point geom.Coord) Orientation {
 	// fast filter for orientation index
 	// avoids use of slow extended-precision arithmetic in many cases
-	index := OrientationIndexFilter(p1, p2, q)
+	index := orientationIndexFilter(vectorOrigin, vectorEnd, point)
 	if index <= 1 {
 		return index
 	}
 
 	var dx1, dy1, dx2, dy2 big.Float
 
 	// normalize coordinates
-	dx1.SetFloat64(p2.X()).Add(&dx1, big.NewFloat(-p1.X()))
-	dy1.SetFloat64(p2.Y()).Add(&dy1, big.NewFloat(-p1.Y()))
-	dx2.SetFloat64(q.X()).Add(&dx2, big.NewFloat(-p2.X()))
-	dy2.SetFloat64(q.Y()).Add(&dy2, big.NewFloat(-p2.Y()))
+	dx1.SetFloat64(vectorEnd.X()).Add(&dx1, big.NewFloat(-vectorOrigin.X()))
+	dy1.SetFloat64(vectorEnd.Y()).Add(&dy1, big.NewFloat(-vectorOrigin.Y()))
+	dx2.SetFloat64(point.X()).Add(&dx2, big.NewFloat(-vectorEnd.X()))
+	dy2.SetFloat64(point.Y()).Add(&dy2, big.NewFloat(-vectorEnd.Y()))
 
-	// calculate determinant.  Calculation takes place in dx1
-	dx1.Mul(&dx1, &dy2).Sub(&dx1, dy1.Mul(&dx1, &dx2))
+	// calculate determinant.  Calculation takes place in dx1 for performance
+	dx1.Mul(&dx1, &dy2)
+	dy1.Mul(&dy1, &dx2)
+	dx1.Sub(&dx1, &dy1)
 
-	return bigSignum(dx1)
+	return Orientation(bigSignum(dx1))
 }
 
 /**
@@ -65,7 +81,7 @@ func OrientationIndex(p1, p2, q geom.Coord) int {
  * @return the orientation index if it can be computed safely
  * @return i > 1 if the orientation index cannot be computed safely
  */
-func OrientationIndexFilter(pa, pb, pc geom.Coord) int {
+func orientationIndexFilter(pa, pb, pc geom.Coord) Orientation {
 	var detsum float64
 
 	detleft := (pa.X() - pc.X()) * (pb.Y() - pc.Y())
@@ -96,18 +112,25 @@ func OrientationIndexFilter(pa, pb, pc geom.Coord) int {
 	return 2
 }
 
-func signum(x float64) int {
+func signum(x float64) Orientation {
 	if x > 0 {
-		return 1
+		return COUNTER_CLOCKWISE
 	}
 	if x < 0 {
-		return -1
+		return CLOCKWISE
 	}
-	return 0
+	return COLLINEAR
 }
-func bigSignum(x big.Float) int {
+func bigSignum(x big.Float) Orientation {
 	if x.IsInf() {
-		return 0
+		return COLLINEAR
+	}
+	switch x.Sign() {
+	case -1:
+		return CLOCKWISE
+	case 0:
+		return COLLINEAR
+	default:
+		return COUNTER_CLOCKWISE
 	}
-	return x.Sign()
 }

algorithm/big_cga/big_cga_test.go

@@ -0,0 +1,57 @@
+package big_cga_test
+
+import (
+	"github.com/twpayne/go-geom"
+	"github.com/twpayne/go-geom/algorithm/big_cga"
+	"testing"
+)
+
+func TestOrientationIndex(t *testing.T) {
+	for i, testData := range []struct {
+		vectorOrigin, vectorEnd, point geom.Coord
+		result                         big_cga.Orientation
+	}{
+
+		{
+			vectorOrigin: geom.Coord{-1.0, -1.0},
+			vectorEnd:    geom.Coord{1.0, 1.0},
+			point:        geom.Coord{0, 0},
+			result:       big_cga.COLLINEAR,
+		},
+		{
+			vectorOrigin: geom.Coord{1.0, 1.0},
+			vectorEnd:    geom.Coord{-1.0, -1.0},
+			point:        geom.Coord{0, 0},
+			result:       big_cga.COLLINEAR,
+		},
+		{
+			vectorOrigin: geom.Coord{10.0, 10.0},
+			vectorEnd:    geom.Coord{20.0, 20.0},
+			point:        geom.Coord{10.0, 20.0},
+			result:       big_cga.COUNTER_CLOCKWISE,
+		},
+		{
+			vectorOrigin: geom.Coord{10.0, 10.0},
+			vectorEnd:    geom.Coord{20.0, 20.0},
+			point:        geom.Coord{20.0, 10.0},
+			result:       big_cga.CLOCKWISE,
+		},
+		{
+			vectorOrigin: geom.Coord{10.0, 20.0},
+			vectorEnd:    geom.Coord{20.0, 10.0},
+			point:        geom.Coord{10.0, 10.0},
+			result:       big_cga.CLOCKWISE,
+		},
+		{
+			vectorOrigin: geom.Coord{10.0, 20.0},
+			vectorEnd:    geom.Coord{20.0, 10.0},
+			point:        geom.Coord{20.0, 20.00},
+			result:       big_cga.COUNTER_CLOCKWISE,
+		},
+	} {
+		orientationIndex := big_cga.OrientationIndex(testData.vectorOrigin, testData.vectorEnd, testData.point)
+		if orientationIndex != testData.result {
+			t.Errorf("Test %v Failed. Expected: %v (%v) but was %v (%v) : TestData: %v", i+1, testData.result, int(testData.result), orientationIndex, int(orientationIndex), testData)
+		}
+	}
+}

algorithm/central_endpoint/intersector.go

@@ -0,0 +1,77 @@
+package central_endpoint
+
+import (
+	"github.com/twpayne/go-geom"
+	"math"
+)
+
+/*
+ * Computes an approximate intersection of two line segments by taking the most central of the endpoints of the segments.
+
+ * This is effective in cases where the segments are nearly parallel and should intersect at an endpoint.
+ * It is also a reasonable strategy for cases where the endpoint of one segment lies on or almost on the interior of another one.
+ * Taking the most central endpoint ensures that the computed intersection point lies in the envelope of the segments.
+
+ * Also, by always returning one of the input points, this should result  in reducing segment fragmentation.
+ * Intended to be used as a last resort for  computing ill-conditioned intersection situations which cause other methods to fail.
+ *
+ */
+func GetIntersection(line1End1, line1End2, line2End1, line2End2 geom.Coord) geom.Coord {
+	intersector := centralEndpointIntersector{
+		line1End1: line1End1,
+		line1End2: line1End2,
+		line2End1: line2End1,
+		line2End2: line2End2}
+	intersector.compute()
+	return intersector.intPt
+}
+
+type centralEndpointIntersector struct {
+	line1End1, line1End2, line2End1, line2End2, intPt geom.Coord
+}
+
+func (intersector *centralEndpointIntersector) compute() {
+	pts := [4]geom.Coord{intersector.line1End1, intersector.line1End2, intersector.line2End1, intersector.line2End2}
+	centroid := average(pts)
+	intersector.intPt = findNearestPoint(centroid, pts)
+}
+
+func average(pts [4]geom.Coord) geom.Coord {
+	n := float64(len(pts))
+	avg := geom.Coord{0, 0}
+
+	for i := 0; i < len(pts); i++ {
+		avg[0] += pts[i][0]
+		avg[1] += pts[i][1]
+	}
+	if n > 0 {
+		avg[0] = avg[0] / n
+		avg[1] = avg[1] / n
+	}
+	return avg
+}
+
+/**
+ * Determines a point closest to the given point.
+ *
+ * @param p the point to compare against
+ * @param p1 a potential result point
+ * @param p2 a potential result point
+ * @param q1 a potential result point
+ * @param q2 a potential result point
+ * @return the point closest to the input point p
+ */
+func findNearestPoint(p geom.Coord, pts [4]geom.Coord) geom.Coord {
+	minDist := math.MaxFloat64
+	result := geom.Coord{}
+	for i := 0; i < len(pts); i++ {
+		dist := p.Distance2D(pts[i])
+
+		// always initialize the result
+		if i == 0 || dist < minDist {
+			minDist = dist
+			result = pts[i]
+		}
+	}
+	return result
+}

algorithm/central_endpoint/intersector_test.go

@@ -0,0 +1,43 @@
+package central_endpoint_test
+
+import (
+	"github.com/twpayne/go-geom"
+	"github.com/twpayne/go-geom/algorithm/central_endpoint"
+	"testing"
+)
+
+func TestGetIntersection(t *testing.T) {
+	for i, tc := range []struct {
+		p1, p2, p3, p4 geom.Coord
+		layout         geom.Layout
+		result         geom.Coord
+	}{
+		{
+			p1: geom.Coord{-1, 0}, p2: geom.Coord{1, 0}, p3: geom.Coord{0, -1}, p4: geom.Coord{0, 1},
+			layout: geom.XY,
+			result: geom.Coord{-1.0, 0.0},
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{10, 20}, p4: geom.Coord{20, 10},
+			layout: geom.XY,
+			result: geom.Coord{10.0, 10.0},
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{20, 20}, p4: geom.Coord{10, 10},
+			layout: geom.XY,
+			result: geom.Coord{10.0, 10.0},
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{30, 20}, p4: geom.Coord{20, 10},
+			layout: geom.XY,
+			result: geom.Coord{20.0, 20.0},
+		},
+	} {
+		calculatedResult := central_endpoint.GetIntersection(tc.p1, tc.p2, tc.p3, tc.p4)
+
+		if !calculatedResult.Equal(tc.layout, tc.result) {
+			t.Errorf("Test '%v' failed: expected \n%v but was \n%v", i+1, tc.result, calculatedResult)
+		}
+
+	}
+}

algorithm/hcoords/hcoords.go

@@ -0,0 +1,43 @@
+package hcoords
+
+import (
+	"fmt"
+	"github.com/twpayne/go-geom"
+	"math"
+)
+
+/**
+ * Computes the (approximate) intersection point between two line segments
+ * using homogeneous coordinates.
+ * <p>
+ * Note that this algorithm is not numerically stable; i.e. it can produce intersection points which
+ * lie outside the envelope of the line segments themselves.  In order to increase the precision of the calculation
+ * input points should be normalized before passing them to this routine.
+ */
+func GetIntersection(p1, p2, q1, q2 geom.Coord) (geom.Coord, error) {
+	// unrolled computation
+	px := p1[1] - p2[1]
+	py := p2[0] - p1[0]
+	pw := p1[0]*p2[1] - p2[0]*p1[1]
+
+	qx := q1[1] - q2[1]
+	qy := q2[0] - q1[0]
+	qw := q1[0]*q2[1] - q2[0]*q1[1]
+
+	x := py*qw - qy*pw
+	y := qx*pw - px*qw
+	w := px*qy - qx*py
+
+	xInt := x / w
+	yInt := y / w
+
+	if math.IsNaN(xInt) || math.IsNaN(xInt) || math.IsNaN(yInt) || math.IsNaN(yInt) {
+		return nil, fmt.Errorf("intersection cannot be calculated using the h-coords implementation")
+	}
+
+	if math.IsInf(xInt, 0) || math.IsInf(xInt, 0) || math.IsInf(yInt, 0) || math.IsInf(yInt, 0) {
+		return nil, fmt.Errorf("intersection cannot be calculated using the h-coords implementation")
+	}
+
+	return geom.Coord{xInt, yInt}, nil
+}

algorithm/hcoords/hcoords_test.go

@@ -0,0 +1,51 @@
+package hcoords_test
+
+import (
+	"github.com/twpayne/go-geom"
+	"github.com/twpayne/go-geom/algorithm/hcoords"
+	"testing"
+)
+
+func TestGetIntersection(t *testing.T) {
+	for i, tc := range []struct {
+		p1, p2, p3, p4 geom.Coord
+		layout         geom.Layout
+		result         geom.Coord
+		err            bool
+	}{
+		{
+			p1: geom.Coord{-1, 0}, p2: geom.Coord{1, 0}, p3: geom.Coord{0, -1}, p4: geom.Coord{0, 1},
+			layout: geom.XY,
+			result: geom.Coord{0.0, 0.0},
+			err:    false,
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{10, 20}, p4: geom.Coord{20, 10},
+			layout: geom.XY,
+			result: geom.Coord{15.0, 15.0},
+			err:    false,
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{20, 20}, p4: geom.Coord{10, 10},
+			layout: geom.XY,
+			result: nil,
+			err:    true,
+		},
+		{
+			p1: geom.Coord{10, 10}, p2: geom.Coord{20, 20}, p3: geom.Coord{30, 20}, p4: geom.Coord{20, 10},
+			layout: geom.XY,
+			result: nil,
+			err:    true,
+		},
+	} {
+		calculatedResult, err := hcoords.GetIntersection(tc.p1, tc.p2, tc.p3, tc.p4)
+
+		if !calculatedResult.Equal(tc.layout, tc.result) {
+			t.Errorf("Test '%v' failed: expected \n%v but was \n%v", i+1, tc.result, calculatedResult)
+		}
+
+		if tc.err != (err != nil) {
+			t.Errorf("Test '%v' failed: expected \n%v but was \n%v", i+1, tc.result, calculatedResult)
+		}
+	}
+}