Fix @turf/nearest-point-on-line endpoint selection and degenerate input cases

packages/turf-nearest-point-on-line/README.md

@@ -6,6 +6,10 @@
 
 Returns the nearest point on a line to a given point.
 
+If any of the segments in the input line string are antipodal and therefore
+have an undefined arc, this function will instead return that the point lies
+on the line.
+
 ### Parameters
 
 *   `lines` **([Geometry][1] | [Feature][2]<([LineString][3] | [MultiLineString][4])>)** lines to snap to

packages/turf-nearest-point-on-line/index.ts

@@ -13,6 +13,10 @@ import { getCoord, getCoords } from "@turf/invariant";
 /**
  * Returns the nearest point on a line to a given point.
  *
+ * If any of the segments in the input line string are antipodal and therefore
+ * have an undefined arc, this function will instead return that the point lies
+ * on the line.
+ *
  * @function
  * @param {Geometry|Feature<LineString|MultiLineString>} lines lines to snap to
  * @param {Geometry|Feature<Point>|number[]} pt point to snap from
@@ -75,12 +79,10 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
       for (let i = 0; i < coords.length - 1; i++) {
         //start - start of current line section
         const start: Feature<Point, { dist: number }> = point(coords[i]);
-        start.properties.dist = distance(pt, start, options);
         const startPos = getCoord(start);
 
         //stop - end of current line section
         const stop: Feature<Point, { dist: number }> = point(coords[i + 1]);
-        stop.properties.dist = distance(pt, stop, options);
         const stopPos = getCoord(stop);
 
         // sectionLength
@@ -89,40 +91,28 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
         let wasEnd: boolean;
 
         // Short circuit if snap point is start or end position of the line
-        // segment or if start is equal to stop position.
-        if (startPos[0] === ptPos[0] && startPos[1] === ptPos[1]) {
-          [intersectPos, , wasEnd] = [startPos, undefined, false];
-        } else if (stopPos[0] === ptPos[0] && stopPos[1] === ptPos[1]) {
-          [intersectPos, , wasEnd] = [stopPos, undefined, true];
-        } else if (startPos[0] === stopPos[0] && startPos[1] === stopPos[1]) {
-          [intersectPos, , wasEnd] = [stopPos, undefined, true];
+        // Test the end position first for consistency in case they are
+        // coincident
+        if (stopPos[0] === ptPos[0] && stopPos[1] === ptPos[1]) {
+          [intersectPos, wasEnd] = [stopPos, true];
+        } else if (startPos[0] === ptPos[0] && startPos[1] === ptPos[1]) {
+          [intersectPos, wasEnd] = [startPos, false];
         } else {
           // Otherwise, find the nearest point the hard way.
-          [intersectPos, , wasEnd] = nearestPointOnSegment(
-            start.geometry.coordinates,
-            stop.geometry.coordinates,
-            getCoord(pt)
+          [intersectPos, wasEnd] = nearestPointOnSegment(
+            startPos,
+            stopPos,
+            ptPos
           );
         }
-        let intersectPt:
-          | Feature<
-              Point,
-              { dist: number; multiFeatureIndex: number; location: number }
-            >
-          | undefined;
-
-        if (intersectPos) {
-          intersectPt = point(intersectPos, {
-            dist: distance(pt, intersectPos, options),
-            multiFeatureIndex: multiFeatureIndex,
-            location: length + distance(start, intersectPos, options),
-          });
-        }
 
-        if (
-          intersectPt &&
-          intersectPt.properties.dist < closestPt.properties.dist
-        ) {
+        const intersectPt = point(intersectPos, {
+          dist: distance(pt, intersectPos, options),
+          multiFeatureIndex: multiFeatureIndex,
+          location: length + distance(start, intersectPos, options),
+        });
+
+        if (intersectPt.properties.dist < closestPt.properties.dist) {
           closestPt = {
             ...intersectPt,
             properties: {
@@ -143,12 +133,7 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
   return closestPt;
 }
 
-/*
- * Plan is to externalise these vector functions to a simple third party
- * library.
- * Possible candidate is @amandaghassaei/vector-math though having some import
- * issues.
- */
+// A simple Vector3 type for cartesian operations.
 type Vector = [number, number, number];
 
 function dot(v1: Vector, v2: Vector): number {
@@ -164,13 +149,13 @@ function cross(v1: Vector, v2: Vector): Vector {
   return [v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x];
 }
 
-function magnitude(v: Vector) {
+function magnitude(v: Vector): number {
   return Math.sqrt(Math.pow(v[0], 2) + Math.pow(v[1], 2) + Math.pow(v[2], 2));
 }
 
-function angle(v1: Vector, v2: Vector): number {
-  const theta = dot(v1, v2) / (magnitude(v1) * magnitude(v2));
-  return Math.acos(Math.min(Math.max(theta, -1), 1));
+function normalize(v: Vector): Vector {
+  const mag = magnitude(v);
+  return [v[0] / mag, v[1] / mag, v[2] / mag];
 }
 
 function lngLatToVector(a: Position): Vector {
@@ -185,7 +170,10 @@ function lngLatToVector(a: Position): Vector {
 
 function vectorToLngLat(v: Vector): Position {
   const [x, y, z] = v;
-  const lat = radiansToDegrees(Math.asin(z));
+  // Clamp the z-value to ensure that is inside the [-1, 1] domain as required
+  // by asin. Note that
+  const zClamp = Math.min(Math.max(z, -1), 1);
+  const lat = radiansToDegrees(Math.asin(zClamp));
   const lng = radiansToDegrees(Math.atan2(y, x));
 
   return [lng, lat];
@@ -195,7 +183,7 @@ function nearestPointOnSegment(
   posA: Position, // start point of segment to measure to
   posB: Position, // end point of segment to measure to
   posC: Position // point to measure from
-): [Position, boolean, boolean] {
+): [Position, boolean] {
   // Based heavily on this article on finding cross track distance to an arc:
   // https://gis.stackexchange.com/questions/209540/projecting-cross-track-distance-on-great-circle
 
@@ -206,62 +194,78 @@ function nearestPointOnSegment(
   const B = lngLatToVector(posB); // ... to posB
   const C = lngLatToVector(posC); // ... to posC
 
-  // Components of target point.
-  const [Cx, Cy, Cz] = C;
-
-  // Calculate coefficients.
-  const [D, E, F] = cross(A, B);
-  const a = E * Cz - F * Cy;
-  const b = F * Cx - D * Cz;
-  const c = D * Cy - E * Cx;
-
-  const f = c * E - b * F;
-  const g = a * F - c * D;
-  const h = b * D - a * E;
+  // The axis (normal vector) of the great circle plane containing the line segment
+  const segmentAxis = cross(A, B);
+
+  // Two degenerate cases exist for the segment axis cross product. The first is
+  // when vectors are aligned (within the bounds of floating point tolerance).
+  // The second is where vectors are antipodal (again within the bounds of
+  // tolerance. Both cases produce a [0, 0, 0] cross product which invalidates
+  // the rest of the algorithm, but each case must be handled separately:
+  // - The aligned case indicates coincidence of A and B. therefore this can be
+  //   an early return assuming the closest point is the end (for consistency).
+  // - The antipodal case is truly degenerate - an infinte number of great
+  //   circles are possible and one will always pass through C. However, given
+  //   that this case is both highly unlikely to occur in practice and that is
+  //   will usually be logically sound to return that the point is on the
+  //   segment, we choose to return the provided point.
+  if (segmentAxis[0] === 0 && segmentAxis[1] === 0 && segmentAxis[2] === 0) {
+    if (dot(A, B) > 0) {
+      return [[...posB], true];
+    } else {
+      return [[...posC], false];
+    }
+  }
 
-  const t = 1 / Math.sqrt(Math.pow(f, 2) + Math.pow(g, 2) + Math.pow(h, 2));
+  // The axis of the great circle passing through the segment's axis and the
+  // target point
+  const targetAxis = cross(segmentAxis, C);
 
-  // Vectors to the two points these great circles intersect.
-  const I1: Vector = [f * t, g * t, h * t];
-  const I2: Vector = [-1 * f * t, -1 * g * t, -1 * h * t];
+  // This cross product also has a degenerate case where the segment axis is
+  // coincident with or antipodal to the target point. In this case the point
+  // is equidistant to the entire segment. For consistency, we early return the
+  // endpoint as the matching point.
+  if (targetAxis[0] === 0 && targetAxis[1] === 0 && targetAxis[2] === 0) {
+    return [[...posB], true];
+  }
 
-  // Figure out which is the closest intersection to this segment of the great
-  // circle.
-  const angleAB = angle(A, B);
-  const angleAI1 = angle(A, I1);
-  const angleBI1 = angle(B, I1);
-  const angleAI2 = angle(A, I2);
-  const angleBI2 = angle(B, I2);
+  // The line of intersection between the two great circle planes
+  const intersectionAxis = cross(targetAxis, segmentAxis);
 
-  let I: Vector;
+  // Vectors to the two points these great circles intersect are the normalized
+  // intersection and its antipodes
+  const I1 = normalize(intersectionAxis);
+  const I2: Vector = [-I1[0], -I1[1], -I1[2]];
 
-  if (
-    (angleAI1 < angleAI2 && angleAI1 < angleBI2) ||
-    (angleBI1 < angleAI2 && angleBI1 < angleBI2)
-  ) {
-    I = I1;
-  } else {
-    I = I2;
-  }
+  // Figure out which is the closest intersection to this segment of the great circle
+  // Note that for points on a unit sphere, the dot product represents the
+  // cosine of the angle between the two vectors which monotonically increases
+  // the closer the two points are together and therefore determines proximity
+  const I = dot(C, I1) > dot(C, I2) ? I1 : I2;
 
   // I is the closest intersection to the segment, though might not actually be
-  // ON the segment.
-
-  // If angle AI or BI is greater than angleAB, I lies on the circle *beyond* A
-  // and B so use the closest of A or B as the intersection
-  if (angle(A, I) > angleAB || angle(B, I) > angleAB) {
-    if (
-      distance(vectorToLngLat(I), vectorToLngLat(A)) <=
-      distance(vectorToLngLat(I), vectorToLngLat(B))
-    ) {
-      return [vectorToLngLat(A), true, false];
-    } else {
-      return [vectorToLngLat(B), false, true];
-    }
+  // ON the segment. To test whether the closest intersection lies on the arc or
+  // not, we do a cross product comparison to check rotation around the unit
+  // circle defined by the great circle plane.
+  const segmentAxisNorm = normalize(segmentAxis);
+  const cmpAI = dot(cross(A, I), segmentAxisNorm);
+  const cmpIB = dot(cross(I, B), segmentAxisNorm);
+
+  // When both comparisons are positive, the rotation from A to I to B is in the
+  // same direction, implying that I lies between A and B
+  if (cmpAI >= 0 && cmpIB >= 0) {
+    return [vectorToLngLat(I), false];
   }
 
-  // As angleAI nor angleBI don't exceed angleAB, I is on the segment
-  return [vectorToLngLat(I), false, false];
+  // Finally process the case where the intersection is not on the segment,
+  // using the dot product with the original point to find the closest endpoint
+  if (dot(A, C) > dot(B, C)) {
+    // Clone position when returning as it is reasonable to not expect structural
+    // sharing on the returned Position in all return cases
+    return [[...posA], false];
+  } else {
+    return [[...posB], true];
+  }
 }
 
 export { nearestPointOnLine };

packages/turf-nearest-point-on-line/test.ts

@@ -14,6 +14,7 @@ import {
   point,
   featureCollection,
   round,
+  degreesToRadians,
 } from "@turf/helpers";
 import { nearestPointOnLine } from "./index.js";
 
@@ -517,3 +518,89 @@ test("turf-nearest-point-on-line -- issue 2808 redundant point support", (t) =>
 
   t.end();
 });
+
+test("turf-nearest-point-on-line -- issue 2934 correct endpoint chosen when in opposite hemisphere", (t) => {
+  // create a long line where the southern end point should be chosen, but the
+  // northern endpoint is closer to the projected great circles
+  const line = lineString([
+    [25, 88],
+    [25, -70],
+  ]);
+  const pt = point([-45, -88]);
+  const nearest = nearestPointOnLine(line, pt);
+
+  t.deepEqual(
+    truncate(nearest, { precision: 8 }).geometry.coordinates,
+    [25, -70],
+    "nearest point is in the Southern Hemisphere"
+  );
+
+  t.end();
+});
+
+test("turf-nearest-point-on-line -- correctly reports external intersection for large arcs", (t) => {
+  // create a test case where the arc is > 2Pi and the closest intersection
+  // point is far from the arc - this case was failing due to an angular
+  // comparison only working for small arcs; this test case checks this
+  // regression
+  const line = lineString([
+    [25, 80],
+    [25, -70],
+  ]);
+  const pt = point([-88, 13]);
+  const nearest = nearestPointOnLine(line, pt);
+
+  t.deepEqual(
+    nearest.geometry.coordinates,
+    [25, 80],
+    "nearest point is at the northern end"
+  );
+
+  t.end();
+});
+
+test("turf-nearest-point-on-line -- issue 2939 nearestPointOnSegment handles tiny segments", (t) => {
+  // create a test case where the line segment passed is small enough such that
+  // a cross product used to generate its normal in nearestPointOnSegment ends
+  // up as [0, 0, 0] but large enough so that the two points are not ===
+  //
+  // In v7.2.0 this test case was failing in browser (Chrome 141.0) but
+  // succeeding on Node (v22.20.0) due to a bit difference in the cosine
+  // implementation. To reproduce the same effect in Node, we monkey patch
+  // Math.cos to temporarily return the browser value as it was too difficult to
+  // find a repro case in node. This monkey patch version fails on v7.2.0 but
+  // passes with degenerate point fixes.
+
+  const lngA = 35.000102519989014;
+  const lngB = 35.00010251998902;
+  const line = lineString([
+    [lngA, 32.00010141921075],
+    [lngB, 32.00010141921075],
+  ]);
+
+  // WARN: Override cos function to give specific results for this test
+  const originalCos = Math.cos;
+  Math.cos = (x) => {
+    // The cosine of lngA gives a different result in the tested browsers than
+    // it does in node, therefore capturing here for this test run only.
+    if (x === degreesToRadians(lngA)) {
+      return originalCos(degreesToRadians(lngB));
+    } else {
+      return originalCos(x);
+    }
+  };
+
+  const pt = point([35.0005, 32.0005]);
+  const nearest = nearestPointOnLine(line, pt);
+
+  // WARN: Reset cos function
+  Math.cos = originalCos;
+
+  t.deepEqual(
+    nearest.geometry.coordinates,
+    [35.00010251998902, 32.00010141921075],
+    "nearest point should be the end point of the line string"
+  );
+
+  t.end();
+});