17 replies to this topic

[broli86]

I think there is definetely something wrong with my code because I get some wierd results( too few intersections). I am using a two step procedure :
1. Find whether ray itnersects triangle plane.
2. Find if the intersection point lies inside the plane.

Heres my code:

// r is ray, i is intersection point, n is normal at intersection point, t is the //triangle intersected


bool ray_triangle_intersection(ray **r, triangle t, vector **i, vector **n)

{

    vector temp;

    // First I calculate the distance using famous formula

    // t =  -(O-A) dot n / (d dot n)  

    vector_sub(&((*r)->origin), &t.a, &temp);  // Finding (O-A)

    double dotnum;

    dotnum = vector_dot(&temp, &t.n); // finding (O-A) dot n

    double dotdenom;

    dotdenom = vector_dot(&((*r)->d), &t.n); // finding  d dot n 

    if(fabs(dotdenom) <= TOLERANCE) // if it is zero then the ray is parallel to 

      return FALSE;                          // plane of triangle

     

     double dist  = -dotnum/dotdenom; // find distance

     if(dist < 0.0)  // if distance is -ve then point is behind the origin of ray

       return FALSE;

     

      (*r)->t = dist;  

      **n = t.n;

      (*i)->x = (*r)->origin.x + dist * (*r)->d.x;   // Calculating intersection 

      (*i)->y = (*r)->origin.y + dist * (*r)->d.y;   // point

      (*i)->z = (*r)->origin.z + dist * (*r)->d.z;

      vector temp1, temp2, temp3;

      double dot00, dot01, dot02;


// condition for intersection point to be inside the triangle is:

//  ((b-a) cross (i-a) ) dot n >= 0

// and ((c-b) cross(i-b)) dot n >= 0

// and ((a-c) cross (i-c)) dot n >=0



      vector_cross(vector_sub(&t.b, &t.a, &temp1),  vector_sub(*i, &t.a, &temp2), &temp3); //  (b-a) cross (i - a) 


      dot00 = vector_dot(&temp3, &t.n); // ((b-a) cross (i-a) ) dot n


      vector_cross(vector_sub(&t.c, &t.b, &temp1),  vector_sub(*i, &t.b, &temp2), &temp3); // (c-b) cross (i-b)


      dot01 = vector_dot(&temp3, &t.n); // ((c-b) cross (i-b) ) dot n 


      vector_cross(vector_sub(&t.a, &t.c, &temp1),  vector_sub(*i, &t.c, &temp2), &temp3); // (a-c) cross (i-c)


      dot02 = vector_dot(&temp3, &t.n); // ((a-c) cross (i-c) ) dot n


      if(dot00 >= 0.0 && dot01 >= 0.0 && dot02 >= 0.0) // condition

      {

		//printf("ip: %f %f %f  closest distance: %f\n", i->x, i->y, i->z, r->t);

        return TRUE;

      }

	    else

	      return FALSE;    

}

[anubis]

http://www.devmaster...le_intersection

[anubis]

Also I think people here are not your debugger, so the least you can do is describe what you think might be wrong and what you did so far to analyze the error. There is little point in us going through your code checking things you already checked.

[roel]

I ain't gonna check your code either, if you simply need a working function, get one from the web, if you want to create your own for fun, debug it for fun. But why are some your arguments pointerpointers?

[Rofar]

I hardly find that looking over that code to see if there is an obvious mistake would equate to debugging the code. All ray-triangle intersections tests do the same thing. It's just a simple matter of looking at the code to see if there is something blatantly incorrect.

Unfortunately, I can't do that without my ray-tri intersect code in front of me to use for a reference (and I'm not where I can get to it right now) or I would give it a look over. I'll try to remember to check it later.

[broli86]

In that article on devmaster for ray triangle intersection using barycentric cordinates, I do not understand this particular statement:

"choose the dominant axis of the triangles normal for the projection"

let's say Z is the dominant axis so does it mean that we project the points onto XY plane ??

[Reedbeta]

Yes.

[broli86]

Reedbeta said:

Yes.

Wouldn't that also mean that based on which plane the points are projected, we have to solve the equations for u and v accordingly ??

For eg. if we project along XY plane , then we have an equation in x and y z is 0. but if we peoject on YZ plane, then we should have an equation in y and z and x is zero ??

I think the article solves the problem considering the point is projected on XY plane.

[anubis]

Quote

Wouldn't that also mean that based on which plane the points are projected, we have to solve the equations for u and v accordingly ??

For eg. if we project along XY plane , then we have an equation in x and y z is 0. but if we peoject on YZ plane, then we should have an equation in y and z and x is zero ??

I think the article solves the problem considering the point is projected on XY plane.

The projection assumes that certain aspects of the triangle are preserved. One is that barycentric coordinates do not change. For numerical reasons it is convenient to pick the dominant axis for the projection.

[anubis]

Quote

I hardly find that looking over that code to see if there is an obvious mistake would equate to debugging the code. All ray-triangle intersections tests do the same thing. It's just a simple matter of looking at the code to see if there is something blatantly incorrect.

Still I think that it's better form if there is some kind of explanation posted along with the code, instead of just putting it there. I mean... it's not very hard to write up something along the lines of : "I think the ray-plane intersection is wrong, because it gives this, that and the other result, for a certain triangle, which is obviously wrong."

[broli86]

anubis said:

The projection assumes that certain aspects of the triangle are preserved. One is that barycentric coordinates do not change. For numerical reasons it is convenient to pick the dominant axis for the projection.

That was not my question. This is what the article has done to compute u,v from linear combinations

u * c + v * b = p
or, as the equal equation system
u * bx + v * cx = px
u * by + v * cy = py
Eliminating u,v, resp., and re-ordering terms we get the results.

note that the above equations are in x and y which means z = 0 which further implies that z was the dominant axis. So shouldn't I have cases for situations where dominant axis is x and y ?

eg if x is dominant axis then x = 0 so the equations are:

u * by + v * cy = py
u * bz + v * cz = pz

Similary if y is the dominant axis then y = 0 so the equations are:

u * bz + v * cz = pz
u * by + v * cy = py


I ask this because in some codes I saw codes using only the equations in terms of x and y (underlying assumption that z is the dominant axis) and solvign them to find u and v.

[anubis]

broli86 said:

That was not my question. This is what the article has done to compute u,v from linear combinations

u * c + v * b = p
or, as the equal equation system
u * bx + v * cx = px
u * by + v * cy = py
Eliminating u,v, resp., and re-ordering terms we get the results.

note that the above equations are in x and y which means z = 0 which further implies that z was the dominant axis. So shouldn't I have cases for situations where dominant axis is x and y ?

eg if x is dominant axis then x = 0 so the equations are:

u * by + v * cy = py
u * bz + v * cz = pz

Similary if y is the dominant axis then y = 0 so the equations are:

u * bz + v * cz = pz
u * by + v * cy = py


I ask this because in some codes I saw codes using only the equations in terms of x and y (underlying assumption that z is the dominant axis) and solvign them to find u and v.

That might just have been an unwise choice of names. It would certainly make more sense to label the projected points differently. I will fix that... thank you.

[broli86]

anubis said:

That might just have been an unwise choice of names. It would certainly make more sense to label the projected points differently. I will fix that... thank you.

I think you misunderstood my question. MY question did not have anything to do with the labelling. I simply asked:

Is the article simply talking (Where they solve the equations to find u and v, look at the middle of the page) about the situation where z = 0 (z is dominant axis)? If yes, then shouldn't we also find u and v for situations where x and y may be zero ? IT will consist of three if loops :

if(dominantaxis == x)
{
solve for u and v

}

if(dominantaxis == y)
{
solve for u and v

}

if(dominantaxis == z)
{
solve for u and v

}

[anubis]

Quote

I think you misunderstood my question. MY question did not have anything to do with the labelling. I simply asked:

Is the article simply talking (Where they solve the equations to find u and v, look at the middle of the page) about the situation where z = 0 (z is dominant axis)? If yes, then shouldn't we also find u and v for situations where x and y may be zero ?

I think now you misunderstood. With the new labeling the projected points are in R2 and x and y are the only coordinates. There is no third coordinate to take care of anymore and the x and y coordinates of the projected points stand for whatever coordinates were not eliminated during the projection.

So to answer your question. No, the article already considers all cases, with projections into all three planes.

[broli86]

anubis said:

I think now you misunderstood. With the new labeling the projected points are in R2 and x and y are the only coordinates. There is no third coordinate to take care of anymore and the x and y coordinates of the projected points stand for whatever coordinates were not eliminated during the projection.

So to answer your question. No, the article already considers all cases, with projections into all three planes.

then is this the proper procedure now:

1. Find the dominant axis(greatest of Nx, Ny and Nz).

2. Depending on dominant axis, make the appropriate component in A, B, C and P as zero.

3. Calculate b, c, p using the values in 2.

4. Using b, c, p , determine the value of u and v with the following 2 equations:

u * bx + v * cx = px
u * by + v * cy = py

[anubis]

broli86 said:

then is this the proper procedure now:

1. Find the dominant axis(greatest of Nx, Ny and Nz).

2. Depending on dominant axis, make the appropriate component in A, B, C and P as zero.

3. Calculate b, c, p using the values in 2.

4. Using b, c, p , determine the value of u and v with the following 2 equations:

u * bx + v * cx = px
u * by + v * cy = py

Yes... that is correct. In my code I create a lookup table : modulo = { 0, 1, 2, 0, 1 }, and then use that to lookup the correct coordinates (dominant + 1) and (dominant + 2), that I need to use as indices into the triangle coordinates (that may not actually be any faster today, than actually computing the right indices). It's also really worth doing the optimization that is mentioned at the end of the article.

[broli86]

anubis said:

Yes... that is correct. In my code I create a lookup table : modulo = { 0, 1, 2, 0, 1 }, and then use that to lookup the correct coordinates (dominant + 1) and (dominant + 2), that I need to use as indices into the triangle coordinates (that may not actually be any faster today, than actually computing the right indices). It's also really worth doing the optimization that is mentioned at the end of the article.

Yeah I used a similar notation :

0 - dominant axis is x axis.
1 - dominant axis is y axis.
2 - dominant axis is z axis.

[anubis]

I thought it might be nice to post a working implementation for reference.

First the data structure :

class TriangleIntersector {
  
    public:
    
        TriangleIntersector(const point3& a, const point3& b, const point3& c);
    
        bool intersect(ray& r, float& u, float &v);
  
    private:
    
        // precomputed data for the intersection test
        unsigned char dominant;
    
        float nu, nv, nd;
        float bnu, bnv, bnd, cnu, cnv, cnd; 
};

This part constructs the precomputed data :

TriangleIntersector::TriangleIntersector(const point3& pa, const point3& pb, const point3& pc)
{
    vector3 b = pc - pa;
    vector3 c = pb - pa;
    
    vector3 n = c.cross(b);
    if(fabs(n.x) > fabs(n.y)){
            
        dominant = (unsigned char)((fabs(n.x) > fabs(n.z)) ? 0 : 2);
    }
    else{
            
        dominant = (unsigned char)((fabs(n.y) > fabs(n.z)) ? 1 : 2);
    }
            
    int u_axis = (dominant + 1) % 3;
    int v_axis = (dominant + 2) % 3;
    
    float ooa = 1.0f / n.v[dominant];
            
    nu = n.v[u_axis] * ooa;
    nv = n.v[v_axis] * ooa;
    nd = n.dot(pa) * ooa;
            
    float ooi = 1.0f / (b.v[u_axis] * c.v[v_axis] - b.v[v_axis] * c.v[u_axis]);
            
    bnu =  b.v[u_axis] * ooi;
    bnv = -b.v[v_axis] * ooi;
    bnd = (b.v[v_axis] * pa.v[u_axis] - b.v[u_axis] * pa.v[v_axis]) * ooi;
            
    cnu =  c.v[v_axis] * ooi;
    cnv = -c.v[u_axis] * ooi;
    cnd = (c.v[u_axis] * pa.v[v_axis] - c.v[v_axis] * pa.v[u_axis]) * ooi;
}

And finally the intersection code :

bool TriangleIntersector::intersect(ray& ray, float& u, float& v)
{
    int u_axis = modulo[dominant + 1];
    int v_axis = modulo[dominant + 2];

    // compute plane distance
    float d        = 1.0f / (ray.direction.v[dominant] + nu * ray.direction.v[u_axis] + nv * ray.direction.v[v_axis]);
    float distance = (nd - ray.origin.v[dominant] - nu * ray.origin.v[u_axis] - nv * ray.origin.v[v_axis]) * d;
            
    if(!(distance < ray.max && distance > ray.min)) return false;
    
    // calculate hitpoint
    float pu = ray.origin.v[u_axis] + distance * ray.direction.v[u_axis];
    float pv = ray.origin.v[v_axis] + distance * ray.direction.v[v_axis];
            
    // solve barycentric coordinate equation in R2
    float a2 = pv * bnu + pu * bnv + bnd;
    if(a2 < 0.0f) return false;
            
    float a3 = pu * cnu + pv * cnv + cnd;
    if(a3 < 0.0f) return false;
            
    if(a2 + a3 > 1.0f){

        return false;
    }
        
    ray.max = distance;
    u            = a2;
    v            = a3;
        
    return true;
}