From c81db372d79bf96acc8a4a2cdb2fe7f0430b57e4 Mon Sep 17 00:00:00 2001
From: Madhura Jayaratne <madhura.cj@gmail.com>
Date: Sat, 9 Jul 2011 19:18:03 +0530
Subject: [PATCH] Manually implemented PointOnSurface() function as it is not
 yet implemented in MySQL

---
 libraries/gis/pma_gis_multipolygon.php | 44 +++-----------------
 libraries/gis/pma_gis_polygon.php      | 57 ++++++++++++++++++++++++++
 2 files changed, 63 insertions(+), 38 deletions(-)

diff --git a/libraries/gis/pma_gis_multipolygon.php b/libraries/gis/pma_gis_multipolygon.php
index 31b9ce5729..4f57c1db59 100644
--- a/libraries/gis/pma_gis_multipolygon.php
+++ b/libraries/gis/pma_gis_multipolygon.php
@@ -376,7 +376,12 @@ class PMA_GIS_Multipolygon extends PMA_GIS_Geometry
             $row_data['parts'][$i]['isOuter'] = PMA_GIS_Polygon::isOuterRing($ring['points']);
         }
 
-        $this->getPointsOnSurface($row_data['parts']);
+        // Find points on surface for inner rings
+        foreach ($row_data['parts'] as $i => $ring) {
+            if (! $ring['isOuter']) {
+                $row_data['parts'][$i]['pointOnSurface'] = PMA_GIS_Polygon::getPointOnSurface($ring['points']);
+            }
+        }
 
         // Classify inner rings to their respective outer rings.
         foreach ($row_data['parts'] as $j => $ring1) {
@@ -429,43 +434,6 @@ class PMA_GIS_Multipolygon extends PMA_GIS_Geometry
         return $wkt;
     }
 
-    /**
-     * Attach to each ring a point on its surface.
-     *
-     * @param array $rings
-     */
-    private function getPointsOnSurface(&$rings)
-    {
-        $sql = 'SELECT ';
-        $inner_rings =  false;
-
-        foreach ($rings as $i => $ring) {
-            if (! $ring['isOuter']) {
-                $inner_rings = true;
-                // gis_data[0]['MULTIPOLYGON'][0][0] will have $ring
-                $gis_data = array(array('MULTIPOLYGON' => array(array($ring['points']))));
-                $gis_data[0]['MULTIPOLYGON'][0][0]['no_of_points'] = count($ring['points']);
-                $geom = $this->generateWkt($gis_data, 0);
-                $sql .= 'AsText(PointOnSurface(GeomFromText("' . $geom . '"))) as `P' . $i . '`, ';
-            }
-        }
-        if (! $inner_rings) {
-            return;
-        }
-
-        $sql = substr($sql, 0, strlen($sql) - 2);
-        $result = PMA_DBI_fetch_result($sql);
-
-        require_once './libraries/gis/pma_gis_point.php';
-        $point = PMA_GIS_Point::singleton();
-        foreach ($rings as $i => $ring) {
-            if (! $ring['isOuter']) {
-                $param = $point->generateParams($result[0]['P' . $i], 0);
-                $rings[$i]['pointOnSurface'] = $param[0]['POINT'];
-            }
-        }
-    }
-
     /**
      * Generate parameters for the GIS data editor from the value of the GIS column.
      *
diff --git a/libraries/gis/pma_gis_polygon.php b/libraries/gis/pma_gis_polygon.php
index f1c051ec40..19e2f53ad0 100644
--- a/libraries/gis/pma_gis_polygon.php
+++ b/libraries/gis/pma_gis_polygon.php
@@ -418,6 +418,63 @@ class PMA_GIS_Polygon extends PMA_GIS_Geometry
         }
     }
 
+    /**
+     * Returns a point that is guaranteed to be on the surface of the ring.
+     * (for simple closed rings)
+     *
+     * @param array $ring array of points forming the ring
+     */
+    public static function getPointOnSurface($ring)
+    {
+        // Find two consecutive distinct points.
+        for ($i = 0; $i < count($ring) - 1; $i++) {
+            if (($ring[$i]['x'] != $ring[$i + 1]['x'])
+                || ($ring[$i]['y'] != $ring[$i + 1]['y'])
+            ) {
+                $x0 = $ring[$i]['x'];
+                $x1 = $ring[$i + 1]['x'];
+                $y0 = $ring[$i]['y'];
+                $y1 = $ring[$i + 1]['y'];
+            }
+        }
+
+        if (! isset($x0)) {
+            return false;
+        }
+
+        // Find the mid point
+        $x2 = ($x0 + $x1) / 2;
+        $y2 = ($y0 + $y1) / 2;
+
+        // Always keep $epsilon < 1 to go with the reduction logic found later in this method
+        $epsilon = 0.1;
+        $denominator = sqrt(pow(($y1 - $y0), 2) + pow(($x0 - $x1), 2));
+        $pointA = array(); $pointB = array();
+
+        while (true) {
+            // Get the points on either sides of the line with a distance of epsilon to the mid point
+            $pointA['x'] = $x2 + ($epsilon * ($y1 - $y0)) / $denominator;
+            $pointA['y'] = $y2 + ($pointA['x'] - $x2) * ($x0 - $x1) / ($y1 - $y0);
+
+            $pointB['x'] = $x2 + ($epsilon * ($y1 - $y0)) / (0 - $denominator);
+            $pointB['y'] = $y2 + ($pointB['x'] - $x2) * ($x0 - $x1) / ($y1 - $y0);
+
+            // One of the points should be inside the polygon, unless epcilon chosen is too large
+            if (PMA_GIS_Polygon::isPointInsidePolygon($pointA, $ring)) {
+                return $pointA;
+            } else if (PMA_GIS_Polygon::isPointInsidePolygon($pointB, $ring)) {
+                return $pointB;
+            }
+
+            // If both are outside the polygon reduce the epsilon and recalculate the points
+            // (reduce exponentially for faster convergance)
+            else {
+                $epsilon = pow($epsilon, 2);
+            }
+
+        }
+    }
+
     /** Generate parameters for the GIS data editor from the value of the GIS column.
      *
      * @param string $value of the GIS column