Minor cleanup for radial sorting

- The radial sort 'less' function is extracted out into its own
  standalone function, making the parent function (`fixVertex`) not need
  to know about the details of the sort.

- A new `XY` method `lengthSq` is added, that computes the square of the
  length of the distance between the origin and the XY value.

- The `Cross` and `lengthSq` methods are used for radial sorting rather
  than reimplementing the logic for those computations.

geom/alg_exact_equals.go

@@ -32,7 +32,7 @@ func (c exactEqualsComparator) eq(a, b Coordinates) bool {
 		return false
 	}
 	asb := a.XY.Sub(b.XY)
-	if asb.Dot(asb) > c.toleranceSq {
+	if asb.lengthSq() > c.toleranceSq {
 		return false
 	}
 	if a.Type.Is3D() && a.Z != b.Z {

geom/alg_point_on_surface.go

@@ -29,7 +29,7 @@ func (n *nearestPointAccumulator) consider(candidate Point) {
 	}
 
 	delta := targetXY.Sub(candidateXY)
-	candidateDist := delta.Dot(delta)
+	candidateDist := delta.lengthSq()
 	if n.point.IsEmpty() || candidateDist < n.dist {
 		n.dist = candidateDist
 		n.point = candidate

geom/dcel_fixup.go

@@ -23,46 +23,12 @@ func (d *doublyConnectedEdgeList) fixVertex(v *vertexRecord) {
 
 	// Perform the sort.
 	if !alreadySorted {
-		// This solution is a reworking of
-		// https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
-		// to avoid using trigonometry.
 		sort.Slice(incidents, func(i, j int) bool {
-			// Sort edges in ascending order of their angle relative to the
-			// x-axis. This is a stricter sort than necessary but is easy to
-			// implement. We only really care that the edges are sorted
-			// relative to each other (we don't care about the starting point).
 			ei := incidents[i]
 			ej := incidents[j]
 			di := ei.seq.GetXY(1).Sub(ei.seq.GetXY(0))
 			dj := ej.seq.GetXY(1).Sub(ej.seq.GetXY(0))
-
-			if di.X >= 0 && dj.X < 0 {
-				return true
-			}
-			if di.X < 0 && dj.X >= 0 {
-				return false
-			}
-			if di.X == 0 && dj.X == 0 {
-				if di.Y >= 0 || dj.Y >= 0 {
-					return di.Y < dj.Y
-				}
-				return dj.Y < di.Y
-			}
-
-			// compute the cross product of vectors from origin
-			det := di.X*dj.Y - dj.X*di.Y
-			if det < 0 {
-				return false
-			}
-			if det > 0 {
-				return true
-			}
-
-			// points are on the same line from the center
-			// check which point is further from the center
-			d1 := di.X*di.X + di.Y*di.Y
-			d2 := dj.X*dj.X + dj.Y*dj.Y
-			return d1 < d2
+			return radialLess(di, dj)
 		})
 	}
 
@@ -75,6 +41,38 @@ func (d *doublyConnectedEdgeList) fixVertex(v *vertexRecord) {
 	}
 }
 
+// radialLess provides an ordering for sorting vectors radially around the origin.
+// This solution is a reworking of
+// https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
+// to avoid using trigonometry.
+func radialLess(di, dj XY) bool {
+	if di.X >= 0 && dj.X < 0 {
+		return true
+	}
+	if di.X < 0 && dj.X >= 0 {
+		return false
+	}
+	if di.X == 0 && dj.X == 0 {
+		if di.Y >= 0 || dj.Y >= 0 {
+			return di.Y < dj.Y
+		}
+		return dj.Y < di.Y
+	}
+
+	// Due to the previous checks, di and dj must be in different sides (LHS vs
+	// RHS) of the XY plane. Therefore the sign of the cross product can
+	// provide an ordering within each half.
+	if det := di.Cross(dj); det != 0 {
+		return det > 0
+	}
+
+	// Points are on the same line from the center.
+	// Check which point is further from the center.
+	li := di.lengthSq()
+	lj := dj.lengthSq()
+	return li < lj
+}
+
 // assignFaces populates the face list based on half edge loops.
 func (d *doublyConnectedEdgeList) assignFaces() {
 	// Find all cycles.

geom/type_line_string.go

@@ -3,7 +3,6 @@ package geom
 import (
 	"database/sql/driver"
 	"fmt"
-	"math"
 	"unsafe"
 
 	"github.com/peterstace/simplefeatures/rtree"
@@ -334,7 +333,7 @@ func (s LineString) Length() float64 {
 		xyA := s.seq.GetXY(i)
 		xyB := s.seq.GetXY(i + 1)
 		delta := xyA.Sub(xyB)
-		sum += math.Sqrt(delta.Dot(delta))
+		sum += delta.Length()
 	}
 	return sum
 }

geom/xy.go

@@ -98,7 +98,12 @@ func (w XY) Unit() XY {
 
 // Length treats XY as a vector, and returns its length.
 func (w XY) Length() float64 {
-	return math.Sqrt(w.Dot(w))
+	return math.Sqrt(w.lengthSq())
+}
+
+// lengthSq treads XY as a vector, and returns its squared length.
+func (w XY) lengthSq() float64 {
+	return w.Dot(w)
 }
 
 // Less gives an ordering on XYs. If two XYs have different X values, then the