Convex hull algorithm to determine points inside of the triangle

Question

I have a bichromatic set S with red and green points in plane general position (no three points of S are collinear). The red triangle is made such that all of its three vertices are red points of S. Green point p of S is said to be enclosed by a red triangle If p lies in the interior of that triangle.

I am supposed to make an algorithm that finds all the green points inside the triangles. I tried to make my own algorithm after seeing a bit from online resources such as: http://totologic.blogspot.com/2014/01/accurate-point-in-triangle-test.html

And here is my "solution":

green_inside = all_the_greens
for each triangle T in get_all_triangles()
  for each greenie G in green_inside
    if G in T
      remove G from green_inside

greenies_in_triangles = all_the_greens - green_inside

I do not know whether this algorithm works. And its runtime is n^3. Also, I am not sure If I can just say get_all_triangles. This algorithm that I am supposed to create should be the most efficient one I can come up with, so please any other solutions or suggestions are welcome!

Answer 1

Let's start with the image shown in the question. It shows red dots which are not part of triangles. And the text gives them no specific meaning. So in my answer, I ignore them.

This transforms the question to:

Given a set of (green) points and a set of triangles, report all points which are inside any of the triangles.

I hope that is what you are looking for.

According to Murphys law, inside every small problem lurks a bunch of bigger problems, struggling to get out. Here, the first bigger problem is, how to test if a point is inside a triangle. Wiki, google etc. show that there are various solutions to the problem. For the given question, I consider the Barycentric coordinate system approach to be a good (the best?) candidate, given we have to retest for each point and this approach allows computing a set of parameters for each triangle, which then is reused in each point test.

The second question is if we can do better than the naive approach in terms of number of tests performed, given K green points and N triangles.

• naive: #Tests = N * K (test each point against all triangles). (I cannot see, how O(N^3) comes into that, but then the author of the question does not say what N is).

• improved: Try to find a way to only test points against a subset of triangles. The subset being triangles which "are likely": #Tests = K * N' where N' hopefully is smaller than or equal to N. K does not change here - after all we need to look at each point.

So, I propose a scan line algorithm as a (potentially) improved solution: Look at the image in the question and imagine a vertical line progressing from left to right over the picture, stepping to the next green dot, whenever it moves.
Now, depending on where that vertical scan line is at a given moment, only triangles which intersect with the scan line are candidates for hit tests of the point at which the scan line is positioned. This is our subset of triangles we want to test the point against.

To implement the scan line algorithm, we need to do some reordering first.

• sort all green points by x coordinate in ascending order. (O (K log K)).
• sort all triangle vertices of each triangle in x-ascending order.
• sort all triangles by the x coordinate of their leftmost points. (O (N log N)).

    Then, for each green point P:
        * drop all of the remaining triangles which are entirely to the left of P.
        * select all of the remaining triangles, which start but do not end to the left of P. 
          (open triangles)
        * test P against each of the open triangles. 
          P is contained if it is contained in any of those open triangles.

While it grew a bit longer than I expected, here is my prototype implementation in everyones favorite language Haskell:

import Data.List
-- While colors are part of the question, we have no real use for them in the solution.
-- data Color = Red | Green deriving (Show,Eq)
-- data Vertex = Vertex { loc :: Point, col :: Color } deriving (Show,Eq)

data Point = Point { x :: Float, y :: Float } deriving (Show,Eq)
data Triangle = Triangle [Point] deriving (Show,Eq)
data BaryParms = BaryParms 
    { dt :: Float   -- the DeTerminant.
    , dy23 :: Float
    , dx32 :: Float
    , dy31 :: Float
    , dx13 :: Float
    , x3 :: Float
    , y3 :: Float
    } deriving (Show,Eq)

-- Get the list of points from a triangle. Of course that list is always of length 3.    
triPoints :: Triangle -> [Point]    
triPoints (Triangle ps) = ps

-- We do some sorting (of triangles and points) and this is our 
-- ordering (left to right on the x-coordinates of a point.)
leftToRight :: Point -> Point -> Ordering
leftToRight p1 p2 
    | x p1 < x p2 = LT
    | x p1 == x p2 = EQ
    | otherwise = GT

-- Shortcut function, since we use it a few times.
-- Left To Right Sort.    
ltrSort = sortBy leftToRight

-- rearrange the points of a triangle in left to right manner.
ltrTriangle :: Triangle -> Triangle
ltrTriangle (Triangle points) = Triangle $ ltrSort points

-- compute the bary coordinates for a point with respect to 
-- the bary parameters of a given triangle.
toBaryCoords :: BaryParms -> Point -> (Float,Float,Float)
toBaryCoords bp p =
    (l1,l2,l3)
    where
        l1 = (dy23 bp * (x p - x3 bp) + (dx32 bp * (y p - y3 bp))) / dt bp
        l2 = (dy31 bp * (x p - x3 bp) + (dx13 bp * (y p - y3 bp))) / dt bp
        l3 = 1 - l1 - l2

-- precompute some values we need for testing points against a triangle.        
triToBaryParms :: Triangle -> BaryParms
triToBaryParms (Triangle [r1,r2,r3]) =
    BaryParms 
    { 
        dt = det (x r1 - x r3) (x r2 - x r3) (y r1 - y r3) (y r2 - y r3), 
        dy23 = y r2 - y r3,
        dx32 = x r3 - x r2,
        dy31 = y r3 - y r1,
        dx13 = x r1 - x r3,
        x3 = x r3,
        y3 = y r3
    }
    where
        det x11 x12 x21 x22 = x11 * x22 - x12 * x21

-- test if a (green) point is contained in a given triangle.
triangleContainsPoint :: Point -> BaryParms -> Bool
triangleContainsPoint p parms  =    
    l1 >= 0.0 && l2 >= 0.0 && l3 >= 0.0
    where
        (l1,l2,l3) = toBaryCoords parms p

-- our improved (?) scan line algorithm.
reportGreensInTriangles :: [Point] -> [Triangle] -> [Point]
reportGreensInTriangles greens triangles =  
    scan sgreens stris []
    where 
        -- sort the green points left to right.
        sgreens :: [Point]
        sgreens = ltrSort greens
        -- sort the triangles left to right by their leftmost vertex.
        -- and compute the barycentric coordinate parameters.
        stris :: [(Triangle,BaryParms)]
        stris = 
            fmap (\t -> (t,triToBaryParms t))
            . sortBy (\t1 t2 -> leftToRight (head . triPoints $ t1) (head . triPoints $ t2)) 
            . fmap ltrTriangle
                $ triangles
        -- recursive iteration over each (green) point. It's tail recursive. 
        scan :: [Point] -> [(Triangle,BaryParms)] -> [Point] -> [Point]
        scan [] _ found = found  -- no more points to test
        scan (g:gs) remainingTris found =
            scan gs remainingTris' found'
            where
                found' = 
                    if any (triangleContainsPoint g . snd) open 
                    then g : found
                    else found
                inScope (Triangle [r1,_,r3],_) = x r1 <= x g && x r3 >= x g
                -- the list of open triangles
                open = takeWhile inScope remainingTris
                -- triangles not discarded yet.
                remainingTris' = dropWhile (\ (Triangle [_, _, r3]  ,_) -> x r3 < x g ) remainingTris

Update

Now, all that is missing is to create an outer function, which lazily creates batches of triangles from the red dots and tests the green dots against that list with the function given above. Any green dot reported by the function above (which is contained in a triangle in that particular call) need not be tested ever again. The outer function creates batches of triangles and calls the function above until either:

• The set of remaining points is empty
• There are no more triangles (worst case).

The set of remaining points in the latter situation is the set of points not contained in any triangle.

Thus, it is easily possible to control the space complexity, since the triangles are being lazily generated in batches.

Answer 2

If all possible triangles (all possible three choices of red points), You can do it simpler:

Create the convex hull (CH) of the red points.
Check each of the green points is in the CH or not.
    If it is, there is a red triangle that contains that point.
    If it is not, there is no red triangle that contains the point.

Creating the convex hull of N red points is in O(N log(N)), and checking a point is inside a convex polygon can be done in O(log(N)). Hence, if the number of green points is M, the total complexity of the above algorithm is in O((M+N)log(N)).