我可以使用 Postgres 函数在固定大小的旋转矩形内查找点吗?
Can I use the Postgres functions to find points inside a rotating rectangle of fixed size?
我使用的是 Postgres 9.5,我刚刚安装了 PostGIS 以获得一些扩展功能。我有一个带有 (x,y) 点的 table,我想找到适合最大点数的矩形。约束是矩形的边长是固定的。到目前为止,我正在计算盒子里没有旋转的点数。我的观点以原点为中心,(0,0).
SELECT Sum(CASE
WHEN x > -5
AND x < 5
AND y > -10
AND y < 10 THEN 1
ELSE 0
END) AS inside_points,
Count(1) AS total_points
FROM track_t;
这个查询给出了原点 (0,0) 和长度 x = 10 和 y = 20.
从这里我将创建一个旋转矩形角点(角度、x1、y1、x2、y2)的助手 table,然后交叉连接到我的数据,并计算每个角度的点数,同时按角度分组。然后我可以 select 哪个角度给我矩形内的点数最多。
但这似乎有点过时,而且可能性能不佳。此外,旋转矩形内的计数点不是 trivial calculation.
有没有更高效更优雅的方法,也许用Postgres Geometric Datatypes or PostGIS Box2D, 旋转固定边长的矩形,然后计数里面的点数? 几何函数看起来不错,但它们似乎提供最小边界框而不是相反。
除了 Postgresql 之外,我还在使用 Python 框架,以防 SQL 无法完成这项工作。
更新:我试过的一件事是使用 Geometric Types,特别是 BOX
SELECT deg, Box(Point(-5, -10), Point(5, 10)) * Point(1, Radians(deg))
FROM Generate_series(0, 360, 90) AS deg
不幸的是,Rotate function by a Point 。
我最终生成了矩形顶点,旋转这些顶点,然后将矩形的面积(常量)与包含测试点的 4 个三角形的面积进行比较。
此技术基于 parsimonious answer:
Make triangle. Suppose, abcd is the rectangle and x is the point then if area(abx)+area(bcx)+area(cdx)+area(dax) equals area(abcd)
then the point is inside it.
矩形由
定义
A左下(-x/2,-y/2)
B左上(-x/2,+y/2)
C右上(+x/2,+y/2)
D右下(+x/2,-y/2)
此代码然后检查点 (qx,qy) 是否在宽度 x=10
和高度 y=20
的矩形内,该矩形围绕原点 (0,0) 旋转了一个角度范围为 0 到 180,以 10 度为单位。
这是代码。检查750k点需要9分钟,所以有一定的改进空间。此外,一旦我升级到 9.6
,它就可以并行化
with t as (select 10*0.5 as x, 20*0.5 as y, 17.0 as qx, -3.0 as qy)
select
z.angle
-- ABC area
--,abs(0.5*(z.ax*(z.by-z.cy)+z.bx*(z.cy-z.ay)+z.cx*(z.ay-z.by)))
-- CDA area
--,abs(0.5*(z.cx*(z.dy-z.ay)+z.dx*(z.ay-z.cy)+z.ax*(z.cy-z.dy)))
-- ABCD area
,abs(0.5*(z.ax*(z.by-z.cy)+z.bx*(z.cy-z.ay)+z.cx*(z.ay-z.by))) + abs(0.5*(z.cx*(z.dy-z.ay)+z.dx*(z.ay-z.cy)+z.ax*(z.cy-z.dy))) as abcd_area
-- ABQ area
--,abs(0.5*(z.ax*(z.by-z.qx)+z.bx*(z.qy-z.ay)+z.qx*(z.ay-z.by)))
-- BCQ area
--,abs(0.5*(z.bx*(z.cy-z.qx)+z.cx*(z.qy-z.by)+z.qx*(z.by-z.cy)))
-- CDQ area
--,abs(0.5*(z.cx*(z.dy-z.qx)+z.dx*(z.qy-z.cy)+z.qx*(z.cy-z.dy)))
-- DAQ area
--,abs(0.5*(z.dx*(z.ay-z.qx)+z.ax*(z.qy-z.dy)+z.qx*(z.dy-z.ay)))
-- total area of triangles with question point (ABQ + BCQ + CDQ + DAQ)
,abs(0.5*(z.ax*(z.by-z.qx)+z.bx*(z.qy-z.ay)+z.qx*(z.ay-z.by)))
+ abs(0.5*(z.bx*(z.cy-z.qx)+z.cx*(z.qy-z.by)+z.qx*(z.by-z.cy)))
+ abs(0.5*(z.cx*(z.dy-z.qx)+z.dx*(z.qy-z.cy)+z.qx*(z.cy-z.dy)))
+ abs(0.5*(z.dx*(z.ay-z.qx)+z.ax*(z.qy-z.dy)+z.qx*(z.dy-z.ay))) as point_area
from
(
SELECT
a.id as angle
-- bottom left (A)
,(-t.x) * cos(radians(a.id)) - (-t.y) * sin(radians(a.id)) as ax
,(-t.x) * sin(radians(a.id)) + (-t.y) * cos(radians(a.id)) as ay
--top left (B)
,(-t.x) * cos(radians(a.id)) - (t.y) * sin(radians(a.id)) as bx
,(-t.x) * sin(radians(a.id)) + (t.y) * cos(radians(a.id)) as by
--top right (C)
,(t.x) * cos(radians(a.id)) - (t.y) * sin(radians(a.id)) as cx
,(t.x) * sin(radians(a.id)) + (t.y) * cos(radians(a.id)) as cy
--bottom right (D)
,(t.x) * cos(radians(a.id)) - (-t.y) * sin(radians(a.id)) as dx
,(t.x) * sin(radians(a.id)) + (-t.y) * cos(radians(a.id)) as dy
-- point to check (Q)
,t.qx as qx
,t.qy as qy
FROM generate_series(0,180,10) AS a(id), t
) z
;
结果是
angle;abcd_area;point_area
0;200;340
10;200;360.6646055963
20;200;373.409049054212
30;200;377.846096908265
40;200;373.84093170467
50;200;361.515248361426
60;200;341.243556529821
70;200;313.641801308188
80;200;279.548648061772
90;200;240
*100;200;200*
*110;200;200*
*120;200;200*
*130;200;200*
*140;200;200*
150;200;237.846096908265
160;200;277.643408923024
170;200;312.04311584956
180;200;340
其中100度、110度、120度、130度和140度的旋转包括测试点(用*
表示)
我使用的是 Postgres 9.5,我刚刚安装了 PostGIS 以获得一些扩展功能。我有一个带有 (x,y) 点的 table,我想找到适合最大点数的矩形。约束是矩形的边长是固定的。到目前为止,我正在计算盒子里没有旋转的点数。我的观点以原点为中心,(0,0).
SELECT Sum(CASE
WHEN x > -5
AND x < 5
AND y > -10
AND y < 10 THEN 1
ELSE 0
END) AS inside_points,
Count(1) AS total_points
FROM track_t;
这个查询给出了原点 (0,0) 和长度 x = 10 和 y = 20.
从这里我将创建一个旋转矩形角点(角度、x1、y1、x2、y2)的助手 table,然后交叉连接到我的数据,并计算每个角度的点数,同时按角度分组。然后我可以 select 哪个角度给我矩形内的点数最多。
但这似乎有点过时,而且可能性能不佳。此外,旋转矩形内的计数点不是 trivial calculation.
有没有更高效更优雅的方法,也许用Postgres Geometric Datatypes or PostGIS Box2D, 旋转固定边长的矩形,然后计数里面的点数? 几何函数看起来不错,但它们似乎提供最小边界框而不是相反。
除了 Postgresql 之外,我还在使用 Python 框架,以防 SQL 无法完成这项工作。
更新:我试过的一件事是使用 Geometric Types,特别是 BOX
SELECT deg, Box(Point(-5, -10), Point(5, 10)) * Point(1, Radians(deg))
FROM Generate_series(0, 360, 90) AS deg
不幸的是,Rotate function by a Point
我最终生成了矩形顶点,旋转这些顶点,然后将矩形的面积(常量)与包含测试点的 4 个三角形的面积进行比较。
此技术基于 parsimonious answer:
Make triangle. Suppose, abcd is the rectangle and x is the point then if
area(abx)+area(bcx)+area(cdx)+area(dax) equals area(abcd)
then the point is inside it.
矩形由
定义A左下(-x/2,-y/2)
B左上(-x/2,+y/2)
C右上(+x/2,+y/2)
D右下(+x/2,-y/2)
此代码然后检查点 (qx,qy) 是否在宽度 x=10
和高度 y=20
的矩形内,该矩形围绕原点 (0,0) 旋转了一个角度范围为 0 到 180,以 10 度为单位。
这是代码。检查750k点需要9分钟,所以有一定的改进空间。此外,一旦我升级到 9.6
,它就可以并行化with t as (select 10*0.5 as x, 20*0.5 as y, 17.0 as qx, -3.0 as qy)
select
z.angle
-- ABC area
--,abs(0.5*(z.ax*(z.by-z.cy)+z.bx*(z.cy-z.ay)+z.cx*(z.ay-z.by)))
-- CDA area
--,abs(0.5*(z.cx*(z.dy-z.ay)+z.dx*(z.ay-z.cy)+z.ax*(z.cy-z.dy)))
-- ABCD area
,abs(0.5*(z.ax*(z.by-z.cy)+z.bx*(z.cy-z.ay)+z.cx*(z.ay-z.by))) + abs(0.5*(z.cx*(z.dy-z.ay)+z.dx*(z.ay-z.cy)+z.ax*(z.cy-z.dy))) as abcd_area
-- ABQ area
--,abs(0.5*(z.ax*(z.by-z.qx)+z.bx*(z.qy-z.ay)+z.qx*(z.ay-z.by)))
-- BCQ area
--,abs(0.5*(z.bx*(z.cy-z.qx)+z.cx*(z.qy-z.by)+z.qx*(z.by-z.cy)))
-- CDQ area
--,abs(0.5*(z.cx*(z.dy-z.qx)+z.dx*(z.qy-z.cy)+z.qx*(z.cy-z.dy)))
-- DAQ area
--,abs(0.5*(z.dx*(z.ay-z.qx)+z.ax*(z.qy-z.dy)+z.qx*(z.dy-z.ay)))
-- total area of triangles with question point (ABQ + BCQ + CDQ + DAQ)
,abs(0.5*(z.ax*(z.by-z.qx)+z.bx*(z.qy-z.ay)+z.qx*(z.ay-z.by)))
+ abs(0.5*(z.bx*(z.cy-z.qx)+z.cx*(z.qy-z.by)+z.qx*(z.by-z.cy)))
+ abs(0.5*(z.cx*(z.dy-z.qx)+z.dx*(z.qy-z.cy)+z.qx*(z.cy-z.dy)))
+ abs(0.5*(z.dx*(z.ay-z.qx)+z.ax*(z.qy-z.dy)+z.qx*(z.dy-z.ay))) as point_area
from
(
SELECT
a.id as angle
-- bottom left (A)
,(-t.x) * cos(radians(a.id)) - (-t.y) * sin(radians(a.id)) as ax
,(-t.x) * sin(radians(a.id)) + (-t.y) * cos(radians(a.id)) as ay
--top left (B)
,(-t.x) * cos(radians(a.id)) - (t.y) * sin(radians(a.id)) as bx
,(-t.x) * sin(radians(a.id)) + (t.y) * cos(radians(a.id)) as by
--top right (C)
,(t.x) * cos(radians(a.id)) - (t.y) * sin(radians(a.id)) as cx
,(t.x) * sin(radians(a.id)) + (t.y) * cos(radians(a.id)) as cy
--bottom right (D)
,(t.x) * cos(radians(a.id)) - (-t.y) * sin(radians(a.id)) as dx
,(t.x) * sin(radians(a.id)) + (-t.y) * cos(radians(a.id)) as dy
-- point to check (Q)
,t.qx as qx
,t.qy as qy
FROM generate_series(0,180,10) AS a(id), t
) z
;
结果是
angle;abcd_area;point_area
0;200;340
10;200;360.6646055963
20;200;373.409049054212
30;200;377.846096908265
40;200;373.84093170467
50;200;361.515248361426
60;200;341.243556529821
70;200;313.641801308188
80;200;279.548648061772
90;200;240
*100;200;200*
*110;200;200*
*120;200;200*
*130;200;200*
*140;200;200*
150;200;237.846096908265
160;200;277.643408923024
170;200;312.04311584956
180;200;340
其中100度、110度、120度、130度和140度的旋转包括测试点(用*
表示)