MySQL实现光线投射算法？

Question

我们需要找到一种快速而相当准确的方法，用于在 Google Maps 上的经纬度值和多边形中确定点是否在多边形内。经过一些研究 - 我们发现了一些关于 MySQL 几何扩展的帖子，并实现了这个功能 -

SELECT id, Contains( PolyFromText( 'POLYGON(".$polygonpath.")' ) , PointFromText( concat( \"POINT(\", latitude, \" \", longitude, \")\" ) ) ) AS
            CONTAINS
FROM tbl_points

然而，这对由大量点组成的多边形无效 :(

经过更多的研究，我发现了一种称为射线投射算法的标准算法，但在尝试在MySQL中开发查询之前，我想冒险看看是否已经有人通过这个算法或者是否有一个有用的链接来展示如何在MySQL/SQL-server中实现该算法。

所以，简而言之 - 问题是：

请问是否能提供Ray casting算法在MySQL/SQL-server中的实现？

其他细节:

• 多边形可能是凸的、凹的或复杂的。
• 目标是快速执行而非100%准确性。

- zarun

---

大约一年前，当我查看MySQL中的地理空间扩展时，它们的质量最多只能达到测试版水平。因此，我最终选择了PostgreSQL作为我的地理数据库。不过，我从未听说过有人尝试在数据库中实现光线投射。 - Eric J.

@EricJ. - 感谢您的回复。我希望我可以在这个问题上使用postGIS...但是由于它只是一个庞大系统中的一小部分，所以无法这样做。 :( - zarun

MySQL有cos/sin/tan函数，这会对你有所帮助吗？链接：http://dev.mysql.com/doc/refman/5.0/en/mathematical-functions.html - Johan

Answer 1

以下函数（Raycasting算法的MYSQL版本）让我大开眼界：

CREATE FUNCTION myWithin(p POINT, poly POLYGON) RETURNS INT(1) DETERMINISTIC 
BEGIN 
DECLARE n INT DEFAULT 0; 
DECLARE pX DECIMAL(9,6); 
DECLARE pY DECIMAL(9,6); 
DECLARE ls LINESTRING; 
DECLARE poly1 POINT; 
DECLARE poly1X DECIMAL(9,6); 
DECLARE poly1Y DECIMAL(9,6); 
DECLARE poly2 POINT; 
DECLARE poly2X DECIMAL(9,6); 
DECLARE poly2Y DECIMAL(9,6); 
DECLARE i INT DEFAULT 0; 
DECLARE result INT(1) DEFAULT 0; 
SET pX = X(p); 
SET pY = Y(p); 
SET ls = ExteriorRing(poly); 
SET poly2 = EndPoint(ls); 
SET poly2X = X(poly2); 
SET poly2Y = Y(poly2); 
SET n = NumPoints(ls); 
WHILE i<n DO 
SET poly1 = PointN(ls, (i+1)); 
SET poly1X = X(poly1); 
SET poly1Y = Y(poly1); 
IF ( ( ( ( poly1X <= pX ) && ( pX < poly2X ) ) || ( ( poly2X <= pX ) && ( pX < poly1X ) ) ) && ( pY > ( poly2Y - poly1Y ) * ( pX - poly1X ) / ( poly2X - poly1X ) + poly1Y ) ) THEN 
SET result = !result; 
END IF; 
SET poly2X = poly1X; 
SET poly2Y = poly1Y; 
SET i = i + 1; 
END WHILE; 
RETURN result; 
End; 

添加

  DELIMITER ;; 

在函数之前添加所需内容。 该函数的使用方法为：

 SELECT myWithin(point, polygon) as result;

在哪里

 point  = Point(lat,lng) 
 polygon = Polygon(lat1 lng1, lat2 lng2, lat3 lng3, .... latn lngn, lat1 lng1)

请注意多边形应该是封闭的（通常如果您检索标准kml或google地图数据，则它已经是封闭的，但请确保 - 注意lat1 lng1设置在末尾重复）

我在数据库中没有点和多边形作为几何字段，所以我不得不做一些类似的事情：

 Select myWithin(PointFromText( concat( "POINT(", latitude, " ", longitude, ")" ) ),PolyFromText( 'POLYGON((lat1 lng1, ..... latn lngn, lat1 lng1))' ) ) as result

我希望这可以帮助到某些人。

- zarun

Answer 2

我将编写自定义UDF，使用C或Delphi或其他高级工具实现射线投射算法： 编写UDF的链接
这是一个MySQL gis实现的源代码，它查找球面上的点（将其用作模板，以了解如何与MySQL中的空间数据类型交互）。
http://www.lenzg.net/archives/220-New-UDF-for-MySQL-5.1-provides-GIS-functions-distance_sphere-and-distance_spheroid.html 来自MySQL手册：
http://dev.mysql.com/doc/refman/5.0/en/adding-functions.html

MS Visual C++的UDF教程：
http://rpbouman.blogspot.com/2007/09/creating-mysql-udfs-with-microsoft.html

Delphi的UDF教程：
在Delphi中创建MySQL的UDF

射线投影算法的源代码
伪代码：http://rosettacode.org/wiki/Ray-casting_algorithm
drDobbs上的文章（注意文章顶部的代码链接）：http://drdobbs.com/cpp/184409586
Delphi（实际上是FreePascal）：http://www.cabiatl.com/mricro/raycast/

- Johan

---

请注意，升级MySQL时自定义UDF可能会出现故障。 - cEz

谢谢！那里有很多有用的链接。我正在逐个查看它们。 此外，当涉及到通知系统时，我遇到了一个问题，该系统应在发布重要内容时通知用户边界内的情况。因此 - 用户（点），事件（点）和边界（多边形）..我想将数百个这样的边界信息传递给函数然后检查是否发生的点在边界内或不在边界内会非常缓慢-然后再次检查所有用户是否在有效边界内并通知他们 ;( - zarun

@zarun，诀窍在于确保您可以使用索引。索引允许您消除99.99％的搜索空间。这就是它们如何加快速度的。 - Johan

@Johan - 谢谢伙计..我承认你的建议。但是我的回复中所说的是关于“使用C、Delphi或其他高级工具实现射线投射算法的UDF”，也就是说，将值传递到php函数中会很慢（我正在为上述问题使用mysql/php）。但没关系，谢谢 :) - zarun

你不能用php编写UDF，它必须编译成DLL。那个DLL你要附加到MySQL上。这只留下生成本机x86代码的工具。 - Johan

谢谢，我明白了。其实我只是在回复原始评论而已。另外，我现在有一个解决方案了 - 在MYSQL中实现了Raycasting！感谢提供的资源！ - zarun

Answer 3

为了保险起见，有一个接受MULTIPOLYGON作为输入的MySQL函数：http://forums.mysql.com/read.php?23,286574,286574

DELIMITER $$

CREATE DEFINER=`root`@`localhost` FUNCTION `GISWithin`(pt POINT, mp MULTIPOLYGON) RETURNS int(1)
    DETERMINISTIC
BEGIN 

DECLARE str, xy TEXT; 
DECLARE x, y, p1x, p1y, p2x, p2y, m, xinters DECIMAL(16, 13) DEFAULT 0; 
DECLARE counter INT DEFAULT 0; 
DECLARE p, pb, pe INT DEFAULT 0; 

SELECT MBRWithin(pt, mp) INTO p; 

IF p != 1 OR ISNULL(p) THEN 
    RETURN p; 
END IF; 

SELECT X(pt), Y(pt), ASTEXT(mp) INTO x, y, str; 
SET str = REPLACE(str, 'POLYGON((',''); 
SET str = REPLACE(str, '))', ''); 
SET str = CONCAT(str, ','); 

SET pb = 1; 
SET pe = LOCATE(',', str); 
SET xy = SUBSTRING(str, pb, pe - pb); 
SET p = INSTR(xy, ' '); 
SET p1x = SUBSTRING(xy, 1, p - 1); 
SET p1y = SUBSTRING(xy, p + 1); 
SET str = CONCAT(str, xy, ','); 

WHILE pe > 0 DO 
    SET xy = SUBSTRING(str, pb, pe - pb); 
    SET p = INSTR(xy, ' '); 
    SET p2x = SUBSTRING(xy, 1, p - 1); 
    SET p2y = SUBSTRING(xy, p + 1); 

    IF p1y < p2y THEN SET m = p1y; ELSE SET m = p2y; END IF; 

    IF y > m THEN 
        IF p1y > p2y THEN SET m = p1y; ELSE SET m = p2y; END IF; 
        IF y <= m THEN 
            IF p1x > p2x THEN SET m = p1x; ELSE SET m = p2x; END IF; 
            IF x <= m THEN 
                IF p1y != p2y THEN 
                    SET xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x; 
                END IF; 
                IF p1x = p2x OR x <= xinters THEN 
                    SET counter = counter + 1; 
                END IF; 
            END IF; 
        END IF; 
    END IF; 

    SET p1x = p2x; 
    SET p1y = p2y; 
    SET pb = pe + 1; 
    SET pe = LOCATE(',', str, pb); 
END WHILE; 

RETURN counter % 2; 

END

- Igor

Answer 4

回复zarun有关在多边形内查找经纬度的功能。

我有一个包含经纬度信息的属性表。因此，我需要获取其经纬度位于多边形lats/longs（我从Google API中获取）内的记录。起初我不知道如何使用Zarun函数。以下是解决方案查询：

• 表格：属性
• 字段：id，纬度，经度，床位等...
• 查询：

SELECT id
FROM properties
WHERE myWithin(
    PointFromText(concat( "POINT(", latitude, " ", longitude, ")")),
    PolyFromText('POLYGON((37.628134 -77.458334,37.629867 -77.449021,37.62324 -77.445416,37.622424 -77.457819,37.628134 -77.458334))' )
) = 1 limit 0,50;

希望这可以为像我这样的新手节省时间 ;)

- Nasir Hussain

---

谢谢！我没想到要发布如何使用它的说明！:-D - zarun

Answer 5

我想在一个多边形表中使用上述存储过程，以下是执行此操作的命令。

使用ogr2ogr将包含多边形的shapefile导入到mysql后，可以按如下方式执行：

ogr2ogr -f "mysql" MYSQL:"mydbname,host=localhost,user=root,password=mypassword,port=3306" -nln "mytablename" -a_srs "EPSG:4326" /path/to/shapefile.shp

接下来，您可以使用MBRwithin对表进行预过滤，然后使用mywithin完成如下操作

DROP TEMPORARY TABLE IF EXISTS POSSIBLE_POLYS;
CREATE TEMPORARY TABLE POSSIBLE_POLYS(OGR_FID INT,SHAPE POLYGON);
INSERT INTO POSSIBLE_POLYS (OGR_FID, SHAPE) SELECT mytablename.OGR_FID,mytablename.SHAPE FROM mytablename WHERE MBRwithin(@testpoint,mytablename.SHAPE);

DROP TEMPORARY TABLE IF EXISTS DEFINITE_POLY;
CREATE TEMPORARY TABLE DEFINITE_POLY(OGR_FID INT,SHAPE POLYGON);
INSERT INTO DEFINITE_POLY (OGR_FID, SHAPE) SELECT POSSIBLE_POLYS.OGR_FID,POSSIBLE_POLYS.SHAPE FROM POSSIBLE_POLYS WHERE mywithin(@testpoint,POSSIBLE_POLYS.SHAPE);

其中@testpoint是从以下位置创建的，例如：

SET @longitude=120;
SET @latitude=-30;
SET @testpoint =(PointFromText( concat( "POINT(", @longitude, " ", @latitude, ")" ) ));

- Mr Purple

Answer 6

这是一个适用于MULTIPOLYGONS的版本（Zarun的版本只适用于POLYGONS）。

CREATE FUNCTION GISWithin(p POINT, multipoly MULTIPOLYGON) RETURNS INT(1) DETERMINISTIC 
BEGIN 
DECLARE n,i,m,x INT DEFAULT 0; 
DECLARE pX,pY,poly1X,poly1Y,poly2X,poly2Y DECIMAL(13,10); 
DECLARE ls LINESTRING; 
DECLARE poly MULTIPOLYGON;
DECLARE poly1,poly2 POINT; 
DECLARE result INT(1) DEFAULT 0; 
SET pX = X(p); 
SET pY = Y(p); 
SET m = NumGeometries(multipoly);
WHILE x<m DO 
SET poly = GeometryN(multipoly,x);
SET ls = ExteriorRing(poly); 
SET poly2 = EndPoint(ls); 
SET poly2X = X(poly2); 
SET poly2Y = Y(poly2); 
SET n = NumPoints(ls); 
WHILE i<n DO 
SET poly1 = PointN(ls, (i+1)); 
SET poly1X = X(poly1); 
SET poly1Y = Y(poly1); 
IF ( ( ( ( poly1X <= pX ) && ( pX < poly2X ) ) || ( ( poly2X <= pX ) && ( pX < poly1X ) ) ) && ( pY > ( poly2Y - poly1Y ) * ( pX - poly1X ) / ( poly2X - poly1X ) + poly1Y ) ) THEN 
SET result = !result; 
END IF; 
SET poly2X = poly1X; 
SET poly2Y = poly1Y; 
SET i = i + 1; 
END WHILE; 
SET x = x + 1;
END WHILE;
RETURN result; 
End; 

- daamsie

Answer 7

从MySQL5.6.1及以上版本开始，它现在是一个空间扩展。请参见文档上的function_st-contains。

- Nasser Al-Wohaibi

---

在StackOverflow上，仅仅一个链接并不是一个答案。请将链接中的信息（例如摘要和代码示例）融入到帖子中。 - M. Mimpen