Given a triangle in the xy plane (specified by its vertices) and a point, find a computationally simple...

I have received a LARGE number of replies to my request for an algorithm,
and would like to give my thanks to all who participated.  Here is a
brief summary of what I received.

1.  General comments.

Many of the solutions violated the "computationally simple" requirement
by involving considerable overhead.  Some computed several square roots or
trig functions, although most of the people who did this made note of this
fact and pointed out exactly how many (how few) calls needed to be made.
Also, many people didn't notice that the algorithm would fail if one of
the lines were vertical or horizontal, resulting in zero division or an
undefined trig value.  Special handling of this, of course, would have made
the code more complex.  Some people even made statements of the form,
". . .and then you just solve the system of equations and you're done."
Naturally, solving a system of equations is not an instantaneous task.
(This reminded me of an old mathematician joke -- the punch line being,
"That reduces the problem to one I've already solved.")

It was also interesting to see that, although there were several
different approaches to the problem, a large number of them converged
to the same set of equations at the end!  This should not be surprising
in itself, but what *was* surprising was the large number of places to
start working from.

Finally, I'd like to point out that I got code in all sorts of
languages ranging from C and Algol to Pascal and FORTRAN.  (No BASIC
or LISP.)  :-)

2.  Methods suggested

The most popular suggestion was to use areas.  The reasoning went as
follows:  given the triangle ABC and point P, compute the areas of the
triangles ABC, ABP, ACP, and BCP.  If the point is inside of ABC, then
the line segments AP, BP, and CP partition it into three subtriangles,
and so the quantity a(ABP) + a(ACP) + a(BCP) should equal a(ABC).  If
it doesn't (in particular, the sum is greater than a(ABC)), the point
is not in the interior of the triangle.

There were two approaches involving angles.  The first (and more
popular) suggested computing the angles APB, BPC, and APC.  If the
point was in the interior of ABC, then APB + BPC + APC = 360 degrees.
The second method, sort of an inverse to this, claimed that the sum of
these angles would be less than 180 degrees if the point were outside the
triangle.

Another popular suggestion involved the use of a particular theorem.
(One person identified this theorem as the Jordan Curve Theorem.)  Simply
stated, a ray emanating from a point crosses the boundary of a closed region
an odd number of times if it is inside the region, or an even number of times
if it is outside.  Coincidentally, all of the people suggesting this method
chose as their ray the one pointing in the same direction as the positive
x-axis.  Obviously, the meat of the computation here lay in determining
whether or not this ray crossed one of the sides of the triangle.

There were several suggestions made by exactly one person (or one
group of collaborators).  These included creative ideas involving dot
products, cross products, determinants, and contour integration.  One
person came up with a scheme in which you systematically insert the point
P into each side of the triangle, making quadrilaterals, and note the
effects.

3.  The chosen solution

I selected this solution as the "best" because of two reasons:  (1)
it was computationally simpler than MOST of the above-mentioned solutions;
(2) although one or two of the solutions may actually be faster than this one,
this solution is (in my opinion) intuitively simpler.  (Enough to justify
the speed difference.)

This solution also took different forms -- vector equations, parametric
equations, determinants, and basic systems of equations.  They all reduced
to the same set of computations.

The basic premise was that, if we define the points A, B, C, and P
with 4 vectors V1, V2, V3, and X, the triangle can be defined by the
equations

        a1*V1 + a2*V2 + a3*V3 = X
        a1    + a2    + a3    = 1

Our equations can be simplified by translating one of the points (say, V3)
to the origin.  (Our situation is now triangle OAB, where we can define
OA = V1 and OB = V2.  O is the origin, of course.)

Points that lie on OA can be written a*V1 for 0 <= a <= 1.  Points on
OB can be written b*V2 for 0 <= b <= 1.  Points on AB are (a*V1 + b*V2)
for a >= 0, b >= 0, a+b = 1.  Points inside the triangle, therefore, are
characterized by

        X = a*V1 + b*V2
        a > 0
        b > 0
        a+b < 1

The vector equation produces two linear equations which, when solved, yields:

        a = (P.x*V2.y - P.y*V2.x) / (V1.x*V2.y - V1.y*V2.x)
        b = (P.x*V1.y - P.y*V1.x) / (V1.y*V2.x - V1.x*V2.y)

(This can be computed faster by computing D = V1.x*V2.y-V1.y*V2.x and
using D and -D for the denominators for a and b, respectively.)

What if the denominator is zero?  Then the points A, B, and C are
collinear, which is impossible if we're given a real triangle.

Complete algorithm:
        translate A, B, and P by (C.x,C.y)
        compute D
        compute a and b
        return (a>0 and b>0 and a+b<1)

[ Editor: "translate by (C.x,C.y)" means subtract C from each point, so
that C moves to the origin. ]

4.  Acknowledgements

Thanks go to:
        Mike Meyer <mwm%opal.Berkeley.EDU@BERKELEY.EDU>
        genrad!teddy!panda!tipsun!jjc@EDDIE.MIT.EDU (Jean J. Cote)
        g-rh@CCA.CCA.COM (Richard Harter)
        oddjob!matt@anl-mcs.ARPA (Matt Crawford)
        Paul Dworkin <mit-eddie!mit-amt!paul@harvard.UUCP>
        postpischil%being.DEC@decwrl.DEC.COM
        Glenn P. Downing <downing@MCC.COM>
        utah-gr!thomas@utah-cs.UUCP (Spencer W. Thomas)
        ags@h.cc.purdue.edu (Dave Seaman)
        pmontgom@sdcrdcf.UUCP (Peter Montgomery)
        Doug Gwyn (VLD/VMB) <gwyn@BRL.ARPA>
        topaz!rutgers!princeton!magic!mie@caip.rutgers.edu (Martin I. Eiger)
        cmcl2!lanl!a.LANL.ARPA!lhl@topaz.UUCP (Lewis Lowe)
        topaz!nike!pyramid!weitek!sci!daver@topaz.rutgers.edu (Dave Rickel)
        sequent!mntgfx!franka@seismo.UUCP (Frank A. Adrian)
plus any others I may have missed.

-- Bill Bulko <bulko@sally.UTEXAS.EDU>

In article <5760@decwrl.DEC.COM>, postpischil@being.dec.com.UUCP writes:

    So, let's define a function f(a,b,c) such that f(a,b,c) =
    x[b]y[a]+x[a]y[c]+x[c]y[b]-(x[a]y[b]+x[c]y[a]+x[b]y[c]).  Then the
    desired predicate, given points v1, v2, and v3 for the triangle and p
    for the point being tested, is:

        (f(v1,v2,c) > 0 and f(v2,v3,c) > 0 and f(v3,v1,c) > 0) or
        (f(v1,v2,c) < 0 and f(v2,v3,c) < 0 and f(v3,v1,c) < 0)

    -- edp

[ Editor: in this predicate the third argument c stands for the test
point p. ]

If the points are ordered either clockwise or counterclockwise, the
signs of the cross product z coordinates of the vectors from the point
in question to adjacent ordered vertices will be the same if the point
is in the triangle.  Treating the triangle as if it were in an
arbitrary z plane = z, define f(a,b,c) as

 f(a,b,c) = (x[a]-x[c])*(y[b]-y[c]) - (x[b]-x[c])*(y[a]-y[c]).

Using the same evaluation criteria as edp, this function call
now requires four less multiplications than edp's function
each time it is called.

-- Jon