What is wrong with my code while trying to calculate if the point "point" is in the triangle or not?

Question

Can anybody see what is wrong with my code?

float intersect(Ray ray, Triangle triangle) {
    float scalar = (0 - dot(triangle.normal, ray.origin)) / dot(triangle.normal, ray.dir);
    vec3 point = ray.origin + scalar * ray.dir;

Berrycentric coordinates

    vec3 A = vec3(1, 0, 0);
    vec3 B = vec3(0, 1, 0);
    vec3 C = vec3(0, 0, 1);

Using the Berrycentric coordinates to calculate tby matrixtby = [A-point] ==> tby = inverse(matrix)[A-point]

    mat3 matrix = mat3(point, A - B, A - C);
    vec3 tby = transpose(matrix) * vec3(A - point);

    float t = tby.x;
    float beta = tby.y;
    float gamma = tby.z;
    if (beta + gamma < 1 && beta > 0 && gamma > 0 && t > 0)
        return scalar;
    return -1.0;
}

What I created so far TriangleStruct

    struct Triangle
    {
     vec3 p1;
     vec3 p2;
     vec3 p3;
     vec3 normal;
     Material material;
    };

The points

    vec3 p1 = vec3(-0.3,0.2,0.5);
    vec3 p2 = vec3(0.3,0.2,0.5);
    vec3 p3 = vec3(0.15,0.0,0.0);
    vec3 p4 = vec3(0.15,0.2,0.5);

The triangles

    {

Triangle 1

    scene.triangles[0].p1 = p1;
    scene.triangles[0].p2 = p4;
    scene.triangles[0].p3 = p3;
    scene.triangles[0].normal = normalize(cross((p4-p1), (p3-p1)));

Triangle 2

    scene.triangles[1].p1 = p3;
    scene.triangles[1].p2 = p2;
    scene.triangles[1].p3 = p1;
    scene.triangles[1].normal = normalize(cross((p2-p3), (p1-p3)));

Triangle 3

    scene.triangles[2].p1 =p3;
    scene.triangles[2].p2 = p4;
    scene.triangles[2].p3 = p2;
    scene.triangles[2].normal = normalize(cross((p4-p3), (p2-p3)));

...

    }

Answer 1

I recommend to rewrite your code like this:

bool PointInOrOn( vec3 P1, vec3 P2, vec3 A, vec3 B )
{
    vec3 CP1 = cross( B - A, P1 - A )
    vec3 CP2 = cross( B - A, P2 - A )
    return step(0.0, dot( CP1, CP2 ));
}

float intersect(Ray ray, Triangle triangle)
{
    vec3 D = normalize(ray.dir);         // skip normalize, if ray.dir is normalized
    vec3 N = normalize(triangle.normal); // skip normalize, if triangle.normal is normalized

    float d = dot(triangle.p1 - ray.origin, N) / dot(D, N)
    vec3  X = ray.origin + D * d;

    float isIn = PointInOrOn( X, triangle.p1, triangle.p2, triangle.p3 ) *
                 PointInOrOn( X, triangle.p2, triangle.p3, triangle.p1 ) *
                 PointInOrOn( X, triangle.p3, triangle.p1, triangle.p2 );

    if ( isIn > 0.01 )
        return d;
    return -1.0;
}

See the following explanation.

---

Intersection of a ray and a triangle primitive

The ray is defined by a point R0 and a normalized direction D.
The plane is defined by a triangle with the three points PA, PB, and PC.

The normal vector of the plane can be calculated by the cross product of 2 legs of the triangle:

N  =  normalize( cross(PC-PA, PB-PA)

The normal distance n of the point R0 to the plane is:

n  =  | R0 - PA | * cos(alpha)  =  dot(PA - R0, N)

It follows that the distance d of the intersection point X to the origin of the ray R0 is:

d  =  n / cos(beta)  =  n / dot(D, N)

The intersection point X is:

X  =  R0 + D * d  =  R0 + D * dot(PA - R0, N) / dot(D, N)

---

To find out, if a point is inside a triangle, has to be tested, if the line from a corner point to the intersection point is between the to legs which are connect to the corner point. The triangle is defined by the points A, B, C and the point to be tested is P:

bool PointInOrOn( P1, P2, A, B )
{
    CP1 = cross( B - A, P1 - A )
    CP2 = cross( B - A, P2 - A )
    return dot( CP1, CP2 ) >= 0
}

bool PointInOrOnTriangle( P, A, B, C )
{
    return PointInOrOn( P, A, B, C ) &&
           PointInOrOn( P, B, C, A ) &&
           PointInOrOn( P, C, A, B )
} 

Answer 2

From what I understood you try to check whether ray goes through triangle. The very first error that I see is that point:

vec3 point = ray.origin + scalar * ray.dir;

has a meaningless definition. You calculate intersection of the ray with a parallel plane to the triangle that goes through the origin. Unless this plane for some miraculous reason coincides with the triangle's plane, all forthcoming calculations with this point is devout of any meaning.

To fix this you need to define scalar like that:

float scalar = dot(triangle.normal,triangle.x) - dot(triangle.normal,ray.origin); 
scalar /= dot(triangle.normal,ray.dir);

where triangle.x is any point of the triangle.

After you fix this issue it will be possible to discuss whether other parts of your code make sense.

Also please provide more information about the details of the code, i.e., relevant pieces of the implementation.

---

Now about how to check whether ray intersects the triangle. The method of intersecting ray with triangle's plane and then checking whether the point is inside the triangle in 2-dimansional setting - is not very good in terms of stability. So avoid it.

A more simple and direct method is to compute the difference vectors:

vec3 dvec1 = tringle.p1 - ray.origin;
vec3 dvec2 = tringle.p2 - ray.origin;
vec3 dvec3 = tringle.p3 - ray.origin;

and then check whether ray.dir can be expressed as a linear sum of dvec1, dvec2, and dvec3 with positive coefficients. This can be achieved by computing the inverse of the matrix mat3(dvec1,dvec2,dvec3) and multiplying it by ray.dir (this way you obtain the coefficients that are needed to express ray.dir linear sum of dvec1, dvec2, and dvec3).

However, the matrix method isn't perfectly stable due divisions. This can be improved further by implementing a logically equivalent code just without divisions.

vec3 dvec12 = cross_product(dvec1, dvec2);
if(dot(dvec12,dvec3)*dot(dvec12,ray.dir) < 0.) return false;

vec3 dvec23 = cross_product(dvec2, dvec3);
if(dot(dvec23,dvec1)*dot(dvec23,ray.dir) < 0.) return false;

vec3 dvec31 = cross_product(dvec3, dvec1);
if(dot(dvec31,dvec2)*dot(dvec31,ray.dir) < 0.) return false;

return true;

In the first if we check if ray.dir is on the same half plane as dvec3 with respect to the plane spanned by {dvec1,dvec2}. If it isn't then ray won't intersect the triangle.

Then we repeat this check for other vectors.