所以，我今天整天都在与一个坦率地说现在令人愤怒的问题作斗争。
给定平面上三角形的一组顶点（只有 3 个点，6 个自由参数），我需要计算这个三角形与 {0,0} 和 {1,1} 定义的单位正方形的交点面积。（我之所以选择这个，是因为 2D 中的任何正方形都可以转换为这个，并且相同的转换可以移动 3 个顶点）。
所以，现在问题被简化为只有 6 个参数，3 个点......我认为这足够短，我愿意编写完整的解决方案/找到完整的解决方案。
（如果可能的话，我希望它在 GPU 上每 <0.5 秒运行超过 200 万个三角形。至于需要简化/没有数据结构/库）
就我对解决方案的尝试而言，我已经......列出了我想出的方法，其中没有一个看起来很快或......特定于好的案例（太笼统）。
选项 1：找到封闭的多边形，它可以是从三角形到 6 边形的任何东西。通过使用我发现的 O(n) 时间算法中的一些凸多边形的交集来做到这一点。然后我会按照 CW 或 CCw 顺序对这些交点（新顶点，最多 7 个 O(n log n) ）进行排序，这样我就可以在这些点上运行一个简单的区域算法（基于 Greens 函数）（ O(n) 再次）。对于与另一个 m-gon 相交的任意凸 n-gon，这是我所能得到的最快速度。但是...我的问题绝对不是那么复杂，它是一个特例，所以它应该有更好的解决方案...
选项2：因为我知道它是一个三角形和单位正方形，所以我可以简单地以更强力的方式找到交点列表（而不是使用某种算法......坦率地说，实现起来有点令人沮丧，如上所列）
只有 19 点需要检查。4个点是三角形内正方形的角。3个点是正方形内的三角形。然后对于三角形的每条线，每条线将与正方形的 4 条线相交（例如 y=0、y=1、x=0、x=1 线）。那是另外12点。所以，12+3+4 = 19 点检查。一旦我拥有最多 6 个，至少 3 个进行此交集的点，我就可以使用我能想到的两种方法之一进行跟进。
2a：通过增加 x 值对它们进行排序，并将形状简单地分解为其子三角形 / 4 边形，每个形状都有一个基于限制顶线和底线的简单公式。总结领域。
或2b：再次以某种循环方式对交叉点进行排序，然后根据greens函数计算面积。
不幸的是，据我所知，这仍然很复杂。为了找到交点，我可以开始进一步分解所有的情况，因为我知道正方形只有 0 和 1，这使得数学放弃了一些术语.. 但这并不一定很简单。
选项3：开始根据各种条件分离问题。例如。正方形内有 0、1、2 或 3 个三角形点。然后对于每种情况，遍历所有可能数量的交点，然后对于每种多边形形状的情况，唯一地写下面积解决方案。
选项 4：一些带有重质阶跃函数的公式。这可能是我最想要的，我怀疑它会有点……大，但也许我乐观地认为这是可能的，一旦我有了公式，这将是最快的计算运行时间。
--- 总的来说，我知道可以使用一些高级库（例如 Clipper）来解决它。我还意识到，在使用各种数据结构（链表，然后对其进行排序）时，编写通用解决方案并不难。如果我只需要这样做几次，所有这些情况都可以。但是，由于我需要将其作为图像处理步骤运行，因此每张图像的顺序为 >9 * 1024*1024 次，并且我正在以 .. 的速度拍摄图像，比如说 1 fps（从技术上讲，我想提高这个速度尽可能快，但计算 900 万个三角形相交面积问题的下限是 1 秒）。这在 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
结尾
