如何计算任意三角形与正方形的交点面积？
blu*_*nut 2 geometry polygons computational-geometry
所以,今天我一直在努力解决一个坦率的现在真气的问题.
给定平面上三角形的一组顶点(仅3个点,6个自由参数),我需要计算该三角形与{0,0}和{1,1}定义的单位平方的交点区域.(我之所以选择这个,是因为2D中的任何方形都可以转换为此方,同一个转换可以移动3个顶点).
所以,现在问题简化为只有6个参数,3个点......我认为这个问题很短,我愿意编写完整的解决方案/找到完整的解决方案.
(如果可能的话,我想在GPU上运行超过200万个三角形,每个<0.5秒.至于需要简化/没有数据结构/库)
就我对解决方案的尝试而言,我已经...列出了我想出的方法,其中没有一个看起来很快或者......特定于好的情况(太笼统).
选项1:找到封闭的多边形,它可以是从三角形到6格的任何东西.通过在我发现的O(n)时间算法中使用一些凸多边形的交叉来做到这一点.然后我会按CW或CCw顺序对这些交点(新顶点,最多7个O(n log n))进行排序,这样我就可以在点上运行一个简单的区域算法(基于格林函数)( O(n)再次).对于与另一个m-gon相交的任意凸n-gon,这是我能提供的最快速度.但是......我的问题绝对不是那么复杂,它是一个特例,所以它应该有一个更好的解决方案......
选项2:因为我知道它是一个三角形和单位正方形,所以我可以用更强力的方式找到交叉点列表(而不是使用一些算法......坦率地说实施起来有点令人沮丧,如上所列)
只有19分要检查.4个点是三角形内的正方形角.正方形内有三角三角形.然后对于三角形的每条线,每条线将与正方形相交4条线(例如,y = 0,y = 1,x = 0,x = 1条线).那是另外12点.所以,12 + 3 + 4 = 19点检查.一旦我有了这个交叉点,最多6个,最少3个点,我可以跟进我能想到的两种方法中的一种.
2a:通过增加x值对它们进行排序,然后简单地将形状分解为其子三角形/ 4-gon形状,每个形状都有一个基于限制顶部和底部线的简单公式.总结一下这些地区.
或2b:再以某种循环方式对交点进行排序,然后根据绿色函数计算区域.
不幸的是,就我所知,这仍然是最复杂的.为了找到交点,我可以开始分解所有的情况,因为我知道它的正方形只有0和1,这使得数学删除了一些条款..但它不一定简单.
选项3:根据各种条件开始分离问题.例如.正方形内的0,1,2或3个三角形点.然后针对每种情况,遍历所有可能数量的交叉点,然后针对每个多边形形状的情况,唯一地记下区域解决方案.
选项4:一些具有重载步骤功能的配方.这是我最想要的那个,我怀疑它会有点......很大,但也许我很乐观这是可能的,并且一旦我有了公式,这将是计算最快的运行时间.
---总的来说,我知道可以使用一些高级库(例如限幅器)来解决它.我还意识到,在使用各种数据结构(链表,然后对其进行排序)时,编写通用解决方案并不是那么难.如果我只需要这样做几次,所有这些情况都可以.但是,因为我需要将它作为一个图像处理步骤运行,每张图像大约> 9*1024*1024次,我正在拍摄图像...让我说1 fps(技术上我会想要推动这个速度尽可能快地上升,但是下限是1秒来计算这些三角形交叉区域问题中的900万个).这在CPU上可能是不可能的,这很好,我可能最终会在Cuda中实现它,但我确实想要在这个问题上推动速度限制.
编辑:所以,我最终选择了2b.由于可能只有19个交叉点,其中最多6个将定义形状,我首先找到3到6个顶点.然后我按循环(CCW)顺序对它们进行排序.然后我通过计算该多边形的面积来找到该区域.
这是我写的测试代码(这是为了Igor,但应该是可读的伪代码)不幸的是它有点长啰嗦,但是......我认为除了我糟糕的排序算法(不应该超过20掉期) ,所以编写更好的排序没有那么多开销)...除了排序之外,我认为我不能让它更快.尽管如此,我对选择此选项时可能有的任何建议或疏忽持开放态度.
function calculateAreaUnitSquare(xPos, yPos)
wave xPos
wave yPos
// First, make array of destination. Only 7 possible results at most for this geometry.
Make/o/N=(7) outputVertexX = NaN
Make/o/N=(7) outputVertexY = NaN
variable pointsfound = 0
// Check 4 corners of square
// Do this by checking each corner against the parameterized plane described by basis vectors p2-p0 and p1-p0.
// (eg. project onto point - p0 onto p2-p0 and onto p1-p0. Using appropriate parameterization scaling (not unit).
// Once we have the parameterizations, then it's possible to check if it is inside the triangle, by checking that u and v are bounded by u>0, v>0 1-u-v > 0
variable denom = yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]*xPos[2]+yPos[1]*xPos[2]+xPos[0]*yPos[2]-xPos[1]*yPos[2]
//variable u00 = yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]*Xx+yPos[1]*Xx+xPos[0]*Yx-xPos[1]*Yx
//variable v00 = -yPos[2]*Xx+yPos[0]*(Xx-xPos[2])+xPos[0]*(yPos[2]-Yx)+yPos[2]*Yx
variable u00 = (yPos[0]*xPos[1]-xPos[0]*yPos[1])/denom
variable v00 = (yPos[0]*(-xPos[2])+xPos[0]*(yPos[2]))/denom
variable u01 =(yPos[0]*xPos[1]-xPos[0]*yPos[1]+xPos[0]-xPos[1])/denom
variable v01 =(yPos[0]*(-xPos[2])+xPos[0]*(yPos[2]-1)+xPos[2])/denom
variable u11 = (yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]+yPos[1]+xPos[0]-xPos[1])/denom
variable v11 = (-yPos[2]+yPos[0]*(1-xPos[2])+xPos[0]*(yPos[2]-1)+xPos[2])/denom
variable u10 = (yPos[0]*xPos[1]-xPos[0]*yPos[1]-yPos[0]+yPos[1])/denom
variable v10 = (-yPos[2]+yPos[0]*(1-xPos[2])+xPos[0]*(yPos[2]))/denom
if(u00 >= 0 && v00 >=0 && (1-u00-v00) >=0)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
if(u01 >= 0 && v01 >=0 && (1-u01-v01) >=0)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
if(u10 >= 0 && v10 >=0 && (1-u10-v10) >=0)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
if(u11 >= 0 && v11 >=0 && (1-u11-v11) >=0)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
// Check 3 points for triangle. This is easy, just see if its bounded in the unit square. if it is, add it.
variable i = 0
for(i=0; i<3; i+=1)
if(xPos[i] >= 0 && xPos[i] <= 1 )
if(yPos[i] >=0 && yPos[i] <=1)
if(!((xPos[i] == 0 || xPos[i] == 1) && (yPos[i] == 0 || yPos[i] == 1) ))
outputVertexX[pointsfound] = xPos[i]
outputVertexY[pointsfound] = yPos[i]
pointsfound+=1
endif
endif
endif
endfor
// Check intersections.
// Procedure is: loop over 3 lines of triangle.
// For each line
// Check if vertical
// If not vertical, find y intercept with x=0 and x=1 lines.
// if y intercept is between 0 and 1, then add the point
// Check if horizontal
// if not horizontal, find x intercept with y=0 and y=1 lines
// if x intercept is between 0 and 1, then add the point
for(i=0; i<3; i+=1)
variable iN = mod(i+1,3)
if(xPos[i] != xPos[iN])
variable tx0 = xPos[i]/(xPos[i] - xPos[iN])
variable tx1 = (xPos[i]-1)/(xPos[i] - xPos[iN])
if(tx0 >0 && tx0 < 1)
variable yInt = (yPos[iN]-yPos[i])*tx0+yPos[i]
if(yInt > 0 && yInt <1)
outputVertexX[pointsfound] = 0
outputVertexY[pointsfound] = yInt
pointsfound+=1
endif
endif
if(tx1 >0 && tx1 < 1)
yInt = (yPos[iN]-yPos[i])*tx1+yPos[i]
if(yInt > 0 && yInt <1)
outputVertexX[pointsfound] = 1
outputVertexY[pointsfound] = yInt
pointsfound+=1
endif
endif
endif
if(yPos[i] != yPos[iN])
variable ty0 = yPos[i]/(yPos[i] - yPos[iN])
variable ty1 = (yPos[i]-1)/(yPos[i] - yPos[iN])
if(ty0 >0 && ty0 < 1)
variable xInt = (xPos[iN]-xPos[i])*ty0+xPos[i]
if(xInt > 0 && xInt <1)
outputVertexX[pointsfound] = xInt
outputVertexY[pointsfound] = 0
pointsfound+=1
endif
endif
if(ty1 >0 && ty1 < 1)
xInt = (xPos[iN]-xPos[i])*ty1+xPos[i]
if(xInt > 0 && xInt <1)
outputVertexX[pointsfound] = xInt
outputVertexY[pointsfound] = 1
pointsfound+=1
endif
endif
endif
endfor
// Now we have all 6 verticies that we need. Next step: find the lowest y point of the verticies
// if there are multiple with same low y point, find lowest X of these.
// swap this vertex to be first vertex.
variable lowY = 1
variable lowX = 1
variable m = 0;
for (i=0; i<pointsfound ; i+=1)
if (outputVertexY[i] < lowY)
m=i
lowY = outputVertexY[i]
lowX = outputVertexX[i]
elseif(outputVertexY[i] == lowY)
if(outputVertexX[i] < lowX)
m=i
lowY = outputVertexY[i]
lowX = outputVertexX[i]
endif
endif
endfor
outputVertexX[m] = outputVertexX[0]
outputVertexY[m] = outputVertexY[0]
outputVertexX[0] = lowX
outputVertexY[0] = lowY
// now we have the bottom left corner point, (bottom prefered).
// calculate the cos(theta) of unit x hat vector to the other verticies
make/o/N=(pointsfound) angles = (p!=0)?( (outputVertexX[p]-lowX) / sqrt( (outputVertexX[p]-lowX)^2+(outputVertexY[p]-lowY)^2) ) : 0
// Now sort the remaining verticies based on this angle offset. This will orient the points for a convex polygon in its maximal size / ccw orientation
// (This sort is crappy, but there will be in theory, at most 25 swaps. Which in the grand sceme of operations, isn't so bad.
variable j
for(i=1; i<pointsfound; i+=1)
for(j=i+1; j<pointsfound; j+=1)
if( angles[j] > angles[i] )
variable tempX = outputVertexX[j]
variable tempY = outputVertexY[j]
outputVertexX[j] = outputVertexX[i]
outputVertexY[j] =outputVertexY[i]
outputVertexX[i] = tempX
outputVertexY[i] = tempY
variable tempA = angles[j]
angles[j] = angles[i]
angles[i] = tempA
endif
endfor
endfor
// Now the list is ordered!
// now calculate the area given a list of CCW oriented points on a convex polygon.
// has a simple and easy math formula : http://www.mathwords.com/a/area_convex_polygon.htm
variable totA = 0
for(i = 0; i<pointsfound; i+=1)
totA += outputVertexX[i]*outputVertexY[mod(i+1,pointsfound)] - outputVertexY[i]*outputVertexX[mod(i+1,pointsfound)]
endfor
totA /= 2
return totA
结束
我认为Cohen-Sutherland线裁剪算法在这里是你的朋友.
首先检查三角形的边界框与正方形,以捕捉琐碎的情况(正方形内的三角形,正方形外的三角形).
接下来检查方块完全位于三角形内的情况.
接下来考虑您的三角形顶点A,B并按C顺时针顺序.剪切线段AB,BC并CA对着方块.它们要么被改变,要么它们位于广场内,要么发现它们位于外面,在这种情况下它们可以被忽略.
您现在有一个最多三个线段的有序列表,用于定义一些边交叉多边形.很容易弄清楚如何从一个边到另一个边遍历以找到交叉多边形的其他边.考虑一个线段(的端点e)对未来的开始(s)
- 如果
e与s三角形顶点位于正方形内的情况一致,则不需要遍历. - 如果
e和s不同,那么我们需要绕着正方形的边界顺时针移动.
请注意,此遍历将按顺时针顺序排列,因此无需计算交点形状的顶点,将它们按顺序排序,然后计算区域.无需存储顶点即可随意计算区域.
请考虑以下示例:
在第一种情况下:
- 我们剪切线条
AB,BC并CA对着正方形,产生线段ab>ba和ca>ac ab>ba形成交叉多边形的第一条边- 从
ba,遍历到ca:ba依赖于y=1,而ca不是,所以下一个边缘是ca>(1,1) (1,1)并且ca两者都在x=1,所以下一个边缘是(1,1)>ca- 下一个边缘是我们已经拥有的线段,
ca>ac ac并且ab是重合的,因此不需要遍历(您可能只是计算退化边缘的区域并避免在这些情况下的分支)
在第二种情况下,将三角形边缘剪切到正方形给我们ab>ba,bc>cb和ca>ac.这些段之间的遍历是微不足道的,因为起点和终点位于相同的方形边缘上.
在第三种情况下,从穿越ba到ca经过两个正方形顶点,但它仍然是在它们所在的广场边缘比较简单的事情:
ba谎言y=1,ca不,所以下一个顶点是(1,1)(1,1)谎言x=1,ca不,所以下一个顶点是(1,0)(1,0)位于y=0一样,ca,那么接下来的顶点ca.