如何计算任意三角形与正方形的交点面积？

所以，我今天一整天都在努力解决一个坦率的现在真气的问题 .

给定平面上三角形的一组顶点（仅3个点，6个自由参数），我需要计算该三角形与{0,0}和{1,1}定义的单位平方的交点区域 . （我之所以选择这个是因为2D中的任何方形都可以转换为此方，并且相同的转换可以移动3个顶点） .

所以，现在问题简化为只有6个参数，3个点......我认为这个问题很短，我愿意编写完整的解决方案/找到完整的解决方案 .

（如果可能的话，我希望这能在每个<0.5秒的GPU上运行超过200万个三角形 . 至于需要简化/没有数据结构/库）

就我对解决方案的尝试而言，我......已经列出了我想出的方法，其中没有一个看起来很快或者......特定于好的情况（太笼统） .

选项1：找到封闭的多边形，它可以是从三角形到6-gon的任何东西 . 通过使用我发现的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：具有重载步骤功能的一些配方 . 这是我最想要的那个，我怀疑它会有点......很大，但也许我很乐观它有可能，并且一旦我有了公式，它将是计算最快的运行时间 .

---总的来说，我知道它可以使用一些高级库（例如clipper）来解决 . 我也意识到，在使用各种数据结构（链表，然后对其进行排序）时，编写通用解决方案并不是那么困难 . 所有这些情况都没问题，如果我只需要这样做几次 . 但是，因为我需要将它作为一个图像处理步骤运行，每张图像大约> 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

结束