calculate the distance from a point to a minor arc with a iterative a…

…lgorithm and add tests
limgrocompany · Jul 25, 2022 · 49b9e55 · 49b9e55
1 parent 2513144
commit 49b9e55
Show file tree

Hide file tree

Showing 3 changed files with 162 additions and 24 deletions.
diff --git a/src/Utility/PointToLineDistance.php b/src/Utility/PointToLineDistance.php
@@ -4,12 +4,17 @@
 
 namespace Location\Utility;
 
+use Location\Bearing\BearingSpherical;
 use Location\Coordinate;
 use Location\Distance\DistanceInterface;
+use Location\Distance\Vincenty;
+use Location\Exception\NotConvergingException;
 use Location\Line;
 
+use function PHPUnit\Framework\throwException;
+
 /**
- * Calculate the distance between a Line and a Point.
+ * Calculate the distance between a Line (minor arc of a Great Circle) and a Point.
  *
  * @author Marcus Jaschen <[email protected]>
  */
@@ -20,39 +25,55 @@ class PointToLineDistance
      */
     private $distanceCalculator;
 
-    public function __construct(DistanceInterface $distanceCalculator)
+    /**
+     * @var float
+     */
+    private $epsilon;
+
+    public function __construct(DistanceInterface $distanceCalculator, float $epsilon = 0.001)
     {
         $this->distanceCalculator = $distanceCalculator;
+        $this->epsilon = $epsilon;
     }
 
     public function getDistance(Coordinate $point, Line $line): float
     {
-        if ($line->getPoint1()->hasSameLocation($line->getPoint2())) {
+        if ($line->getPoint1()->hasSameLocation($line->getPoint2(), $this->epsilon)) {
             return $this->distanceCalculator->getDistance($point, $line->getPoint1());
         }
+        if ($point->hasSameLocation($line->getPoint1(), $this->epsilon)) {
+            return 0.0;
+        }
+        if ($point->hasSameLocation($line->getPoint2(), $this->epsilon)) {
+            return 0.0;
+        }
 
-        $pLat = deg2rad($point->getLat());
-        $pLng = deg2rad($point->getLng());
-
-        $l1Lat = deg2rad($line->getPoint1()->getLat());
-        $l1Lng = deg2rad($line->getPoint1()->getLng());
-        $l2Lat = deg2rad($line->getPoint2()->getLat());
-        $l2Lng = deg2rad($line->getPoint2()->getLng());
+        $iterationsCounter = 0;
+        $iterationLine = clone $line;
 
-        $deltal2l1Lat = $l2Lat - $l1Lat;
-        $deltal2l1Lng = $l2Lng - $l1Lng;
+        do {
+            $linePoint1 = $iterationLine->getPoint1();
+            $linePoint2 = $iterationLine->getPoint2();
+            $lineMidPoint = $iterationLine->getMidpoint();
 
-        $u = (($pLat - $l1Lat) * $deltal2l1Lat + ($pLng - $l1Lng) * $deltal2l1Lng) /
-            ($deltal2l1Lat ** 2 + $deltal2l1Lng ** 2);
+            $distancePointToLinePoint1 = $point->getDistance($linePoint1, $this->distanceCalculator);
+            $distancePointToLinePoint2 = $point->getDistance($linePoint2, $this->distanceCalculator);
+            $distancePointToLineMidPoint = $point->getDistance($lineMidPoint, $this->distanceCalculator);
 
-        if ($u <= 0) {
-            return $this->distanceCalculator->getDistance($point, $line->getPoint1());
-        }
+            if (($distancePointToLinePoint1 + $distancePointToLineMidPoint) <= ($distancePointToLineMidPoint + $distancePointToLinePoint2)) {
+                $iterationLine = new Line($linePoint1, $lineMidPoint);
+            } else {
+                $iterationLine = new Line($lineMidPoint, $linePoint2);
+            }
 
-        if ($u >= 1) {
-            return $this->distanceCalculator->getDistance($point, $line->getPoint2());
-        }
+            if (abs($distancePointToLinePoint1 - $distancePointToLinePoint2) < $this->epsilon) {
+                return $point->getDistance($iterationLine->getMidpoint(), $this->distanceCalculator);
+            }
 
-        return (new PerpendicularDistance())->getPerpendicularDistance($point, $line);
+            $iterationsCounter++;
+            if ($iterationsCounter > 100) {
+                throw new NotConvergingException('Calculation of Point to Minor Arc did not converge after 100 iterations.', 6391878367);
+            }
+        } while (true);
     }
 }
diff --git a/tests/Location/Utility/PointToLineDistanceTest.php b/tests/Location/Utility/PointToLineDistanceTest.php
@@ -52,6 +52,17 @@ public function testLinePoint1IsNearer(): void
         );
     }
 
+    public function testPointIsSameLocationAsLinePoint1(): void
+    {
+        $point = new Coordinate(52.45, 13.05);
+
+        $linePoint1 = new Coordinate(52.45, 13.05);
+        $linePoint2 = new Coordinate(52.6, 13.12);
+        $line = new Line($linePoint1, $linePoint2);
+
+        $this->assertEquals(0, $this->pointToLineDistance->getDistance($point, $line));
+    }
+
     public function testLinePoint2IsNearer(): void
     {
         $point = new Coordinate(52.6001, 13.1201);
@@ -60,12 +71,36 @@ public function testLinePoint2IsNearer(): void
         $linePoint2 = new Coordinate(52.6, 13.12);
         $line = new Line($linePoint1, $linePoint2);
 
-        $this->assertEquals(
+        $this->assertEqualsWithDelta(
             $point->getDistance($linePoint2, $this->vincenty),
-            $this->pointToLineDistance->getDistance($point, $line)
+            $this->pointToLineDistance->getDistance($point, $line),
+            0.001
         );
     }
 
+    public function testPointIsSameLocationAsLinePoint2(): void
+    {
+        $point = new Coordinate(52.45, 13.05);
+
+        $linePoint1 = new Coordinate(52.5, 13.1);
+        $linePoint2 = new Coordinate(52.45, 13.05);
+        $line = new Line($linePoint1, $linePoint2);
+
+        $this->assertEquals(0, $this->pointToLineDistance->getDistance($point, $line));
+    }
+
+    public function testPointIsSameLocationAsLineMidPoint(): void
+    {
+
+        $linePoint1 = new Coordinate(52.5, 13.1);
+        $linePoint2 = new Coordinate(52.6, 13.12);
+        $line = new Line($linePoint1, $linePoint2);
+
+        $point = $line->getMidpoint();
+
+        $this->assertEquals(0, $this->pointToLineDistance->getDistance($point, $line));
+    }
+
     public function testDistanceIsCalculatedToSomewhereOnLine(): void
     {
         $point = new Coordinate(52.04, 13.01);
@@ -95,6 +130,18 @@ public function testDistanceMatchesPerpendicularDistance(): void
         $perpendicularDistance = $pdCalculator->getPerpendicularDistance($point, $line);
         $pointToLineDistance = $this->pointToLineDistance->getDistance($point, $line);
 
-        $this->assertEquals($pointToLineDistance, $perpendicularDistance);
+        // allowed delta is relatively large because the perpdendicular distance
+        // is calculated with a simple spherical model
+        $this->assertEqualsWithDelta($pointToLineDistance, $perpendicularDistance, 0.3);
+    }
+
+    public function testPointToLineDistanceCaseA()
+    {
+        $line = new Line(new Coordinate(55.9857757, 13.5475307), new Coordinate(55.9869533, 13.5509295));
+        $point0a = new Coordinate(55.9839105, 13.5465958);
+
+        $this->assertEqualsWithDelta(215.636, $point0a->getDistance($line->getPoint1(), new Vincenty()), 0.1);
+
+        $this->assertEqualsWithDelta(215.636, $this->pointToLineDistance->getDistance($point0a, $line), 0.1);
     }
 }
diff --git a/tests/Regression/Github/92/Issue92Test.php b/tests/Regression/Github/92/Issue92Test.php
@@ -0,0 +1,70 @@
+<?php
+
+declare(strict_types=1);
+
+namespace Location;
+
+use Location\Distance\Haversine;
+use Location\Distance\Vincenty;
+use Location\Factory\CoordinateFactory;
+use Location\Utility\PointToLineDistanceOld;
+use Location\Utility\PointToLineDistance;
+use PHPUnit\Framework\TestCase;
+
+class Issue92Test extends TestCase
+{
+    public function testForIssue92(): void
+    {
+        $ll1 = array (
+            'lat' => 55.98467000751469413444283418357372283935546875,
+            'lon' => 13.544430000324179985682349069975316524505615234375,
+        );
+        $ll2 = array (
+            'lat' => 55.98461833972562118333371472544968128204345703125,
+            'lon' => 13.54429500218709137016048771329224109649658203125,
+        );
+        $p = array (
+            'lat' => 55.97609299999999876717993174679577350616455078125,
+            'lon' => 13.5475989999999999469082467840053141117095947265625,
+        );
+
+        $pointToLineDistanceCalculator = new PointToLineDistance(new Vincenty());
+
+        $dist = $pointToLineDistanceCalculator->getDistance(
+            new Coordinate($p['lat'], $p['lon']),
+            new Line(
+                new Coordinate($ll1['lat'], $ll1['lon']),
+                new Coordinate($ll2['lat'], $ll2['lon'])
+            )
+        );
+
+        $this->assertEqualsWithDelta(971.37, $dist, 0.01);
+    }
+    public function testForIssue92RoundedValues(): void
+    {
+        $ll1 = array (
+            'lat' => 55.98467,
+            'lon' => 13.54443,
+        );
+        $ll2 = array (
+            'lat' => 55.98461,
+            'lon' => 13.54429,
+        );
+        $p = array (
+            'lat' => 55.97609,
+            'lon' => 13.54759,
+        );
+
+        $pointToLineDistanceCalculator = new PointToLineDistance(new Vincenty());
+
+        $dist = $pointToLineDistanceCalculator->getDistance(
+            new Coordinate($p['lat'], $p['lon']),
+            new Line(
+                new Coordinate($ll1['lat'], $ll1['lon']),
+                new Coordinate($ll2['lat'], $ll2['lon'])
+            )
+        );
+
+        $this->assertEqualsWithDelta(970.74, $dist, 0.01);
+    }
+}