How to calculate area of intersection of an arbitrary triangle with a square?

Question

So, I've been struggling with a frankly now infuriating problem all day today.

Given a set of verticies of a triangle on a plane (just 3 points, 6 free parameters), I need to calculate the area of intersection of this triangle with the unit square defined by {0,0} and {1,1}. (I choose this because any square in 2D can be transformed to this, and the same transformation can move the 3 vertices).

So, now the problem is simplified down to only 6 parameters, 3 points... which I think is short enough that I'd be willing to code up the full solution / find the full solution.

( I would like this to run on a GPU for literally more than 2 million triangles every <0.5 seconds, if possible. as for the need for simplification / no data structures / libraries)

In terms of my attempt at the solution, I've... got a list of ways I've come up with, none of which seem fast or ... specific to the nice case (too general).

Option 1: Find the enclosed polygon, it can be anything from a triangle up to a 6-gon. Do this by use of some intersection of convex polygon in O(n) time algorithms that I found. Then I would sort these intersection points (new vertices, up to 7 of them O(n log n) ), in either CW or CCw order, so that I can run a simple area algorithm on the points (based on Greens function) (O(n) again). This is the fastest i can come with for an arbitrary convex n-gon intersecting with another m-gon. However... my problem is definitely not that complex, its a special case, so it should have a better solution...

Option 2: Since I know its a triangle and unit square, i can simply find the list of intersection points in a more brute force way (rather than using some algorithm that is ... frankly a little frustrating to implement, as listed above)

There are only 19 points to check. 4 points are corners of square inside of triangle. 3 points are triangle inside square. And then for each line of the triangle, each will intersect 4 lines from the square (eg. y=0, y=1, x=0, x=1 lines). that is another 12 points. so, 12+3+4 = 19 points to check. Once I have the, at most 6, at fewest 3, points that do this intersection, i can then follow up with one of two methods that I can think of.

2a: Sort them by increasing x value, and simply decompose the shape into its sub triangle / 4-gon shapes, each with an easy formula based on the limiting top and bottom lines. sum up the areas.

or 2b: Again sort the intersection points in some cyclic way, and then calculate the area based on greens function.

Unfortunately, this still ends up being just as complex as far as I can tell. I can start breaking up all the cases a little more, for finding the intersection points, since i know its just 0s and 1s for the square, which makes the math drop out some terms.. but it's not necessarily simple.

Option 3: Start separating the problem based on various conditions. Eg. 0, 1, 2, or 3 points of triangle inside square. And then for each case, run through all possible number of intersections, and then for each of those cases of polygon shapes, write down the area solution uniquely.

Option 4: some formula with heaviside step functions. This is the one I want the most probably, I suspect it'll be a little... big, but maybe I'm optimistic that it is possible, and that it would be the fastest computationally run time once I have the formula.

--- Overall, I know that it can be solved using some high level library (clipper for instance). I also realize that writing general solutions isn't so hard when using data structures of various kinds (linked list, followed by sorting it). And all those cases would be okay, if I just needed to do this a few times. But, since I need to run it as an image processing step, on the order of >9 * 1024*1024 times per image, and I'm taking images at .. lets say 1 fps (technically I will want to push this speed up as fast as possible, but lower bound is 1 second to calculate 9 million of these triangle intersection area problems). This might not be possible on a CPU, which is fine, I'll probably end up implementing it in Cuda anyways, but I do want to push the limit of speed on this problem.

Edit: So, I ended up going with Option 2b. Since there are only 19 intersections possible, of which at most 6 will define the shape, I first find those 3 to 6 verticies. Then i sort them in a cyclic (CCW) order. And then I find the area by calculating the area of that polygon.

Here is my test code I wrote to do that (it's for Igor, but should be readable as pseudocode) Unfortunately it's a little long winded, but.. I think other than my crappy sorting algorithm (shouldn't be more than 20 swaps though, so not so much overhead for writing better sorting)... other than that sorting, I don't think I can make it any faster. Though, I am open to any suggestions or oversights I might have had in chosing this option.

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

end

Answer 1

I think the Cohen-Sutherland line-clipping algorithm is your friend here.

First off check the bounding box of the triangle against the square to catch the trivial cases (triangle inside square, triangle outside square).

Next check for the case where the square lies completely within the triangle.

Next consider your triangle vertices A, B and C in clockwise order. Clip the line segments AB, BC and CA against the square. They will either be altered such that they lie within the square or are found to lie outside, in which case they can be ignored.

You now have an ordered list of up to three line segments that define the some of the edges intersection polygon. It is easy to work out how to traverse from one edge to the next to find the other edges of the intersection polygon. Consider the endpoint of one line segment (e) against the start of the next (s)

• If e is coincident with s, as would be the case when a triangle vertex lies within the square, then no traversal is required.
• If e and s differ, then we need to traverse clockwise around the boundary of the square.

Note that this traversal will be in clockwise order, so there is no need to compute the vertices of the intersection shape, sort them into order and then compute the area. The area can be computed as you go without having to store the vertices.

Consider the following examples:

In the first case:

1. We clip the lines AB, BC and CA against the square, producing the line segments ab>ba and ca>ac
2. ab>ba forms the first edge of the intersection polygon
3. To traverse from ba to ca: ba lies on y=1, while ca does not, so the next edge is ca>(1,1)
4. (1,1) and ca both lie on x=1, so the next edge is (1,1)>ca
5. The next edge is a line segment we already have, ca>ac
6. ac and ab are coincident, so no traversal is needed (you might be as well just computing the area for a degenerate edge and avoiding the branch in these cases)

In the second case, clipping the triangle edges against the square gives us ab>ba, bc>cb and ca>ac. Traversal between these segments is trivial as the start and end points lie on the same square edges.

In the third case the traversal from ba to ca goes through two square vertices, but it is still a simple matter of comparing the square edges on which they lie:

1. ba lies on y=1, ca does not, so next vertex is (1,1)
2. (1,1) lies on x=1, ca does not, so next vertex is (1,0)
3. (1,0) lies on y=0, as does ca, so next vertex is ca.

Answer 2

Given the large number of triangles I would recommend scanline algorithm: sort all the points 1st by X and 2nd by Y, then proceed in X direction with a "scan line" that keeps a heap of Y-sorted intersections of all lines with that line. This approach has been widely used for Boolean operations on large collections of polygons: operations such as AND, OR, XOR, INSIDE, OUTSIDE, etc. all take O(n*log(n)).

It should be fairly straightforward to augment Boolean AND operation, implemented with the scanline algorithm to find the areas you need. The complexity will remain O(n*log(n)) on the number of triangles. The algorithm would also apply to intersections with arbitrary collections of arbitrary polygons, in case you would need to extend to that.

On the 2nd thought, if you don't need anything other than the triangle areas, you could do that in O(n), and scanline may be an overkill.

Answer 3

I came to this question late, but I think I've come up with a more fully flushed out solution along the lines of ryanm's answer. I'll give an outline of for others trying to do this problem at least somewhat efficiently.

First you have two trivial cases to check:

1) Triangle lies entirely within the square

2) Square lies entirely within the triangle (Just check if all corners are inside the triangle)

If neither is true, then things get interesting.

First, use either the Cohen-Sutherland or Liang-Barsky algorithm to clip each edge of the triangle to the square. (The linked article contains a nice bit of code that you can essentially just copy-paste if you're using C).

Given a triangle edge, these algorithms will output either a clipped edge or a flag denoting that the edge lies entirely outside the square. If all edges lie outsize the square, then the triangle and the square are disjoint.

Otherwise, we know that the endpoints of the clipped edges constitute at least some of the vertices of the polygon representing the intersection.

We can avoid a tedious case-wise treatment by making a simple observation. All other vertices of the intersection polygon, if any, will be corners of the square that lie inside the triangle.

Simply put, the vertices of the intersection polygon will be the (unique) endpoints of the clipped triangle edges in addition to the corners of the square inside the triangle.

We'll assume that we want to order these vertices in a counter-clockwise fashion. Since the intersection polygon will always be convex, we can compute its centroid (the mean over all vertex positions) which will lie inside the polygon.

Then to each vertex, we can assign an angle using the atan2 function where the inputs are the y- and x- coordinates of the vector obtained by subtracting the centroid from the position of the vertex (i.e. the vector from the centroid to the vertex).

Finally, the vertices can be sorted in ascending order based on the values of the assigned angles, which constitutes a counter-clockwise ordering. Successive pairs of vertices correspond to the polygon edges.