我正在使用Postgres9.5,并且刚刚安装了PostGIS以获得一些扩展功能。我有一个带有(x,y)点的表,我想找到适合最大点数的矩形。约束条件是矩形边长是固定的。到目前为止,我正在计算框中没有旋转的点数。我的点以原点 (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)组成的辅助表,然后交叉连接到我的数据,并对每个角度的点进行计数,而按角度分组。然后,我可以选择哪个角度可以使矩形内的点最多。
是否有更有效,更优雅的方法 (例如使用Postgres Geometric数据类型或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
不幸的是,“按点旋转”功能不适用于Polygons。
最后,我生成了矩形顶点,旋转了这些顶点,然后将矩形(常数)的面积与通过包含测试点而形成的4个三角形的面积进行比较。
此技术基于简约的答案:
做成三角形。假设abcd是矩形,x是点,则如果area(abx)+area(bcx)+area(cdx)+area(dax) equals area(abcd)该点在其内部。
area(abx)+area(bcx)+area(cdx)+area(dax) equals area(abcd)
矩形定义为
甲 左下方(-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度。
x=10
y=20
这是代码。检查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度的旋转包括测试点(用表示*)
*