POV-Ray : Newsgroups : povray.programming : Explanation of point inside a triangle code Server Time
29 Mar 2024 11:21:38 EDT (-0400)
  Explanation of point inside a triangle code (Message 1 to 5 of 5)  
From: queercore
Subject: Explanation of point inside a triangle code
Date: 6 Apr 2009 11:30:00
Message: <web.49da1f19de1ed1545807d1b10@news.povray.org>
Hi,

I am an undergraduate student wiritng a thesis report on pov-ray version 3.6.0.
I have been studying pov-ray source code but since there is no description or
explanation about many parts of the code, I have beem facing lots of
difficulties. I have been trying to figure out the point inside a triangle
test. I understood the plane intersection part, which is simple but cant figure
out the point inside the triangle test. I googled to find out what algorithm was
implemented but didnt find any similarity with any algorithm I have searched.
So, can someone please explain the test to me.

Here is the code ray triangle intersection...

static int Intersect_Triangle(RAY *Ray, TRIANGLE *Triangle, DBL *Depth)
{
  DBL NormalDotOrigin, NormalDotDirection;
  DBL s, t;

  Increase_Counter(stats[Ray_Triangle_Tests]);

  if (Test_Flag(Triangle, DEGENERATE_FLAG))
  {
    return(false);
  }

  VDot(NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);

  if (fabs(NormalDotDirection) < EPSILON)
  {
    return(false);
  }

  VDot(NormalDotOrigin, Triangle->Normal_Vector, Ray->Initial);

  *Depth = -(Triangle->Distance + NormalDotOrigin) / NormalDotDirection;

  if ((*Depth < DEPTH_TOLERANCE) || (*Depth > Max_Distance))
  {
    return(false);
  }

  switch (Triangle->Dominant_Axis)
  {
    case X:

      s = Ray->Initial[Y] + *Depth * Ray->Direction[Y];
      t = Ray->Initial[Z] + *Depth * Ray->Direction[Z];

      if ((Triangle->P2[Y] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
          (Triangle->P2[Z] - t) * (Triangle->P2[Y] - Triangle->P1[Y]))
      {
        return(false);
      }

      if ((Triangle->P3[Y] - s) * (Triangle->P3[Z] - Triangle->P2[Z]) <
          (Triangle->P3[Z] - t) * (Triangle->P3[Y] - Triangle->P2[Y]))
      {
        return(false);
      }

      if ((Triangle->P1[Y] - s) * (Triangle->P1[Z] - Triangle->P3[Z]) <
          (Triangle->P1[Z] - t) * (Triangle->P1[Y] - Triangle->P3[Y]))
      {
        return(false);
      }

      Increase_Counter(stats[Ray_Triangle_Tests_Succeeded]);

      return(true);

    case Y:

      s = Ray->Initial[X] + *Depth * Ray->Direction[X];
      t = Ray->Initial[Z] + *Depth * Ray->Direction[Z];

      if ((Triangle->P2[X] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
          (Triangle->P2[Z] - t) * (Triangle->P2[X] - Triangle->P1[X]))
      {
        return(false);
      }

      if ((Triangle->P3[X] - s) * (Triangle->P3[Z] - Triangle->P2[Z]) <
          (Triangle->P3[Z] - t) * (Triangle->P3[X] - Triangle->P2[X]))
      {
        return(false);
      }

      if ((Triangle->P1[X] - s) * (Triangle->P1[Z] - Triangle->P3[Z]) <
          (Triangle->P1[Z] - t) * (Triangle->P1[X] - Triangle->P3[X]))
      {
        return(false);
      }

      Increase_Counter(stats[Ray_Triangle_Tests_Succeeded]);

      return(true);

    case Z:

      s = Ray->Initial[X] + *Depth * Ray->Direction[X];
      t = Ray->Initial[Y] + *Depth * Ray->Direction[Y];

      if ((Triangle->P2[X] - s) * (Triangle->P2[Y] - Triangle->P1[Y]) <
          (Triangle->P2[Y] - t) * (Triangle->P2[X] - Triangle->P1[X]))
      {
        return(false);
      }

      if ((Triangle->P3[X] - s) * (Triangle->P3[Y] - Triangle->P2[Y]) <
          (Triangle->P3[Y] - t) * (Triangle->P3[X] - Triangle->P2[X]))
      {
        return(false);
      }

      if ((Triangle->P1[X] - s) * (Triangle->P1[Y] - Triangle->P3[Y]) <
          (Triangle->P1[Y] - t) * (Triangle->P1[X] - Triangle->P3[X]))
      {
        return(false);
      }

      Increase_Counter(stats[Ray_Triangle_Tests_Succeeded]);

      return(true);
  }

  return(false);
}



Thanks in advance.


Post a reply to this message

From: clipka
Subject: Re: Explanation of point inside a triangle code
Date: 6 Apr 2009 12:35:01
Message: <web.49da2ed2b10008142b82b7c80@news.povray.org>
"queercore" <nomail@nomail> wrote:
> Here is the code ray triangle intersection...


> static int Intersect_Triangle(RAY *Ray, TRIANGLE *Triangle, DBL *Depth)
> {
>   DBL NormalDotOrigin, NormalDotDirection;
>   DBL s, t;
>
>   Increase_Counter(stats[Ray_Triangle_Tests]);
>
>   if (Test_Flag(Triangle, DEGENERATE_FLAG))
>   {
>     return(false);
>   }
>
>   VDot(NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
>
>   if (fabs(NormalDotDirection) < EPSILON)
>   {
>     return(false);
>   }
>
>   VDot(NormalDotOrigin, Triangle->Normal_Vector, Ray->Initial);
>
>   *Depth = -(Triangle->Distance + NormalDotOrigin) / NormalDotDirection;
>
>   if ((*Depth < DEPTH_TOLERANCE) || (*Depth > Max_Distance))
>   {
>     return(false);
>   }

The above code mainly computes *Depth, which tells us the distance along the ray
at which it intersects the plane of the triangle. Aside from that, it's just
framework stuff and sanity checks.


Next thing we do is test whether the resulting intersection point is inside or
outside the triangle. As the problem is invariant to linear transformations
(except for pathological cases), and projection onto an arbitrary plane *is* a
linear transformation (being pathological to this problem only if the plane is
perpendicular to the triangle's plane), we take the liberty to project the
whole smash onto a plane of our liking.

Obviously, the XY, XZ and YZ planes lend themselves to this purpose, as
projecting onto these is brain-dead easy; to avoid pathological or
near-pathological cases, we choose between these three according to the most
dominant axis of the triangle's surface normal:

>   switch (Triangle->Dominant_Axis)
>   {
>
>     case X:

Determine (projection of) intersection point (s;t):

>       s = Ray->Initial[Y] + *Depth * Ray->Direction[Y];
>       t = Ray->Initial[Z] + *Depth * Ray->Direction[Z];

Use straightforward 2D linear math to check whether we're "offside" any of the
three lines along the (projected) triangle's edges:

>       if ((Triangle->P2[Y] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
>           (Triangle->P2[Z] - t) * (Triangle->P2[Y] - Triangle->P1[Y]))
>       {
>         return(false);
>       }
>
>       if ((Triangle->P3[Y] - s) * (Triangle->P3[Z] - Triangle->P2[Z]) <
>           (Triangle->P3[Z] - t) * (Triangle->P3[Y] - Triangle->P2[Y]))
>       {
>         return(false);
>       }
>
>       if ((Triangle->P1[Y] - s) * (Triangle->P1[Z] - Triangle->P3[Z]) <
>           (Triangle->P1[Z] - t) * (Triangle->P1[Y] - Triangle->P3[Y]))
>       {
>         return(false);
>       }

That's it; it may appear disappointingly cheap, but it does the job ;). The
remainder is just framework stuff again, and the same algorithm for the other
two planes of our choice.

>       Increase_Counter(stats[Ray_Triangle_Tests_Succeeded]);
>
>       return(true);
>
>     case Y:
>       [...]
>
>     case Z:
>       [...]
>
>   }
>
>   return(false);
> }


Post a reply to this message

From: queercore
Subject: Re: Explanation of point inside a triangle code
Date: 6 Apr 2009 14:20:00
Message: <web.49da468fb1000814440d46cd0@news.povray.org>
"clipka" <nomail@nomail> wrote:
> "queercore" <nomail@nomail> wrote:
> > Here is the code ray triangle intersection...
>
>
> > static int Intersect_Triangle(RAY *Ray, TRIANGLE *Triangle, DBL *Depth)
> > {
> >   DBL NormalDotOrigin, NormalDotDirection;
> >   DBL s, t;
> >
> >   Increase_Counter(stats[Ray_Triangle_Tests]);
> >
> >   if (Test_Flag(Triangle, DEGENERATE_FLAG))
> >   {
> >     return(false);
> >   }
> >
> >   VDot(NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
> >
> >   if (fabs(NormalDotDirection) < EPSILON)
> >   {
> >     return(false);
> >   }
> >
> >   VDot(NormalDotOrigin, Triangle->Normal_Vector, Ray->Initial);
> >
> >   *Depth = -(Triangle->Distance + NormalDotOrigin) / NormalDotDirection;
> >
> >   if ((*Depth < DEPTH_TOLERANCE) || (*Depth > Max_Distance))
> >   {
> >     return(false);
> >   }
>
> The above code mainly computes *Depth, which tells us the distance along the ray
> at which it intersects the plane of the triangle. Aside from that, it's just
> framework stuff and sanity checks.
>
>
> Next thing we do is test whether the resulting intersection point is inside or
> outside the triangle. As the problem is invariant to linear transformations
> (except for pathological cases), and projection onto an arbitrary plane *is* a
> linear transformation (being pathological to this problem only if the plane is
> perpendicular to the triangle's plane), we take the liberty to project the
> whole smash onto a plane of our liking.
>
> Obviously, the XY, XZ and YZ planes lend themselves to this purpose, as
> projecting onto these is brain-dead easy; to avoid pathological or
> near-pathological cases, we choose between these three according to the most
> dominant axis of the triangle's surface normal:
>
> >   switch (Triangle->Dominant_Axis)
> >   {
> >
> >     case X:
>
> Determine (projection of) intersection point (s;t):
>
> >       s = Ray->Initial[Y] + *Depth * Ray->Direction[Y];
> >       t = Ray->Initial[Z] + *Depth * Ray->Direction[Z];
>
> Use straightforward 2D linear math to check whether we're "offside" any of the
> three lines along the (projected) triangle's edges:
>
> >       if ((Triangle->P2[Y] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
> >           (Triangle->P2[Z] - t) * (Triangle->P2[Y] - Triangle->P1[Y]))
> >       {
> >         return(false);
> >       }
> >
> >       if ((Triangle->P3[Y] - s) * (Triangle->P3[Z] - Triangle->P2[Z]) <
> >           (Triangle->P3[Z] - t) * (Triangle->P3[Y] - Triangle->P2[Y]))
> >       {
> >         return(false);
> >       }
> >
> >       if ((Triangle->P1[Y] - s) * (Triangle->P1[Z] - Triangle->P3[Z]) <
> >           (Triangle->P1[Z] - t) * (Triangle->P1[Y] - Triangle->P3[Y]))
> >       {
> >         return(false);
> >       }
>
> That's it; it may appear disappointingly cheap, but it does the job ;). The
> remainder is just framework stuff again, and the same algorithm for the other
> two planes of our choice.
>
> >       Increase_Counter(stats[Ray_Triangle_Tests_Succeeded]);
> >
> >       return(true);
> >
> >     case Y:
> >       [...]
> >
> >     case Z:
> >       [...]
> >
> >   }
> >
> >   return(false);
> > }




Thank you very much for your reply. I understood the projection part but can you
please tell me what is the 2D linear math involved in finding whether the point
is inside or offside of the triangle edges.


Post a reply to this message

From: clipka
Subject: Re: Explanation of point inside a triangle code
Date: 6 Apr 2009 15:15:00
Message: <web.49da53e8b10008142b82b7c80@news.povray.org>
"queercore" <nomail@nomail> wrote:
> Thank you very much for your reply. I understood the projection part but can you
> please tell me what is the 2D linear math involved in finding whether the point
> is inside or offside of the triangle edges.

Okay, let's look at this in detail:

       if ((Triangle->P2[Y] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
           (Triangle->P2[Z] - t) * (Triangle->P2[Y] - Triangle->P1[Y]))
       {
         return(false);
       }

It becomes clearer if we flip the signs:

       if ((s - Triangle->P2[Y]) * (Triangle->P1[Z] - Triangle->P2[Z]) <
           (t - Triangle->P2[Z]) * (Triangle->P1[Y] - Triangle->P2[Y]))
       {
         return(false);
       }

Note that all the differences in there actually constitute a translation by -P2,
so we get a statement of the form:

    if (A[Y] * B[Z] < A[Z] * B[Y]) ...

which is a more robust way of comparing A[Y]/A[Z] with B[Y]/B[Z], which gets
unnecessarily nasty in certain cases (and, as it seems, would be more
problematic regarding the signs, and choice of > vs. <).

So in essence the idea is to compare the "direction coefficients" (don't know if
that's the proper English word) of the line equations of (A) the line through P2
and the intersection point, and (B) the line through P2 and P1.

You can also think of it in terms of a vector dot product:

    if (A[Y] * B[Z] < A[Z] * B[Y]) ...

is equal to:

    if (A[Y] * B[Z] + A[Z] * (-B[Y]) < 0) ...

which happens to be the dot product of A and a vector perpendicular to B
(compared with zero).

Interpreted this way, the test boils down to checking whether the angle between
A, and that particular perpendicular to B, is greater or smaller than 90
degrees.


Post a reply to this message

From: queercore
Subject: Re: Explanation of point inside a triangle code
Date: 6 Apr 2009 16:30:00
Message: <web.49da6553b1000814c66344120@news.povray.org>
> Okay, let's look at this in detail:
>
>        if ((Triangle->P2[Y] - s) * (Triangle->P2[Z] - Triangle->P1[Z]) <
>            (Triangle->P2[Z] - t) * (Triangle->P2[Y] - Triangle->P1[Y]))
>        {
>          return(false);
>        }
>
> It becomes clearer if we flip the signs:
>
>        if ((s - Triangle->P2[Y]) * (Triangle->P1[Z] - Triangle->P2[Z]) <
>            (t - Triangle->P2[Z]) * (Triangle->P1[Y] - Triangle->P2[Y]))
>        {
>          return(false);
>        }
>
> Note that all the differences in there actually constitute a translation by -P2,
> so we get a statement of the form:
>
>     if (A[Y] * B[Z] < A[Z] * B[Y]) ...
>
> which is a more robust way of comparing A[Y]/A[Z] with B[Y]/B[Z], which gets
> unnecessarily nasty in certain cases (and, as it seems, would be more
> problematic regarding the signs, and choice of > vs. <).
>
> So in essence the idea is to compare the "direction coefficients" (don't know if
> that's the proper English word) of the line equations of (A) the line through P2
> and the intersection point, and (B) the line through P2 and P1.
>
> You can also think of it in terms of a vector dot product:
>
>     if (A[Y] * B[Z] < A[Z] * B[Y]) ...
>
> is equal to:
>
>     if (A[Y] * B[Z] + A[Z] * (-B[Y]) < 0) ...
>
> which happens to be the dot product of A and a vector perpendicular to B
> (compared with zero).
>
> Interpreted this way, the test boils down to checking whether the angle between
> A, and that particular perpendicular to B, is greater or smaller than 90
> degrees.


Thanks a lot. I get it now.


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.