skotinator 122 Report post Posted March 24, 2008 Hi, I've spent the last few hours trying to get the seemingly-simple ray-triangle barycentric intersection routine going in my raytracer with absolutely no luck - and getting some very whacky results. Stepping through the code seems to yield some strange values for 'distance', which should contain the distance along the ray from the origin that the plane intersects at. Given that i am testing against a cube defined from z=24 to 30 and some distances are > 1000.0f im not sure whats going on - especially seeing as more or less the same distance routine is used elsewhere fine to calculate distance to planes. Would really appreciate some fresh eyes on this, its driving me insane :) float OBJTriangleProxy::testIntersection (Ray &r, float maxDist) { Vector3d O = r.getOrigin(); Vector3d D = r.getDirection(); Vector3d A = getVert (0); Vector3d B = getVert (1); Vector3d C = getVert (2); //Calculate average normal Vector3d N = getNormal(0) + getNormal(1) + getNormal(2); N.normalise(); //Calculate distance from ray to the triangle plane float distance = ((O - A).dot (N)) / (D.dot (N)); //Ray is behind the plane the triangle's plane if (distance < 0) return -1; //Work out dominant axis for projection of triangle points into 2d int dominant = 0; if (abs(N.getX()) > abs(N.getY())) { if (abs(N.getX() > abs(N.getZ()))) dominant = 0; //x dominant else dominant = 2; //z dominant } else { if (abs(N.getY() > abs(N.getZ()))) dominant = 1; //y dominant else dominant = 2; // z dominant } //Find point where the ray intersects the plane Vector3d P = O + D * distance; //Project points onto appropriate plane B = C - A; C = B - A; switch (dominant) { case 0: //Project onto x B = Vector3d (B.getY(), A.getZ(), 0); C = Vector3d (C.getY(), C.getZ(), 0); P = Vector3d (P.getY(), P.getZ(), 0); break; case 1: //Project onto Y B = Vector3d (B.getX(), A.getZ(), 0); C = Vector3d (C.getX(), C.getZ(), 0); P = Vector3d (P.getX(), P.getZ(), 0); break; case 2: //Project onto Z B = Vector3d (B.getX(), A.getY(), 0); C = Vector3d (C.getX(), C.getY(), 0); P = Vector3d (P.getX(), P.getY(), 0); break; } //Calculate 2d barycentric components float u = ((P.getY() * C.getX()) - (P.getX()*B.getY())) / ((B.getY() * C.getX()) - (B.getX() - C.getY())); float v = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY())); //Point in triangle? if (u >= 0 && v >= 0 && u + v <= 1) return distance; return -1; } 0 Share this post Link to post Share on other sites
Pragma 395 Report post Posted March 24, 2008 Don't average the vertex normals. Use the true geometric normal of the triangle.N = normalize(cross(B-A, C-A));Also make sure that these lines really do what you expect them toB = C - A;C = B - A; 0 Share this post Link to post Share on other sites
haegarr 7373 Report post Posted March 24, 2008 Presumbly not directly related to the problem, but: Check that the absolute value of D.dot(N) is greater than some epsilon before computing the distance. If you don't do so and the ray is nearly parallel to the plane, then a division by zero will happen. 0 Share this post Link to post Share on other sites
haegarr 7373 Report post Posted March 24, 2008 Oops, herein are 2 inconsistencies I've not seen during my first visit:float u = ((P.getY() * C.getX()) - (P.getX()*B.getY())) / ((B.getY() * C.getX()) - (B.getX() - C.getY()));float v = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY())); Shouldn't the last part of u be like B.getX() * C.getY() or so? Furthurmore, P.getX()*B.getY() appearing in both formulas seems me curious, too. 0 Share this post Link to post Share on other sites
skotinator 122 Report post Posted March 25, 2008 Thanks for the input!I've fixed the normal calculation to use cross product (although this didn't matter as the test mesh im using uses the same normals for each vertex on a face anyway). Also, I forgot to multiply the distance by negative 1, and the calculation of the second triangle edge just before projection onto an axis-plane was incorrect (thanks pragma!). Still getting wierd results, although starting to look a bit more normal.The plane/ray intersection is working correctly now so im sure my problem is in the projection/baycentric part. //Project points onto appropriate plane Vector3d B2 = B; B = C - A; C = B2 - A; switch (dominant) { case 0: //Project onto x B = Vector3d (B.getY(), B.getZ(), 0); C = Vector3d (C.getY(), C.getZ(), 0); P = Vector3d (P.getY(), P.getZ(), 0); break; case 1: //Project onto Y B = Vector3d (B.getX(), B.getZ(), 0); C = Vector3d (C.getX(), C.getZ(), 0); P = Vector3d (P.getX(), P.getZ(), 0); break; case 2: //Project onto Z B = Vector3d (B.getX(), B.getY(), 0); C = Vector3d (C.getX(), C.getY(), 0); P = Vector3d (P.getX(), P.getY(), 0); break; } //Calculate 2d barycentric components float u = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY())); float v = ((C.getY() * P.getX()) - (C.getX()*P.getY())) / ((B.getX() * C.getY()) - (B.getY() * C.getX())); //Point in triangle? if (u >= 0 && v >= 0 && u + v <= 1) return distance;Will have more of a play with it tonight :)[Edited by - skotinator on March 27, 2008 12:48:46 AM] 0 Share this post Link to post Share on other sites
haegarr 7373 Report post Posted March 25, 2008 There is still a (C.getX() - B.getY())in the computation of u, where I expect a (C.getX() * B.getY())instead. 0 Share this post Link to post Share on other sites
skotinator 122 Report post Posted March 25, 2008 I can't believe i missed that!Well it is printing out correctly proportioned triangles now, but they are not correctly positioned.for instancev 30.0 0.0 20.0v 15.0 10.0 20.0v 10.0 0.0 20.0Is the right shape, but has the right most edge (which should be at x=30) at x=0. If i change tov 50.0 0.0 20.0v 15.0 10.0 20.0v 10.0 0.0 20.0the triangle is deformed as you'd expect, but the extra length shifts the left hand edge, not the right - which remains fixed at x=0.Generally, for the XY plane the triangle is rendered into, an extension of X or Y of a certain amount results in the triangle stretching in the opposite direction to accomodate - while keeping the most centered point on the origin (ie x=0 y=0).I think my barycentric equations are to blame 0 Share this post Link to post Share on other sites
haegarr 7373 Report post Posted March 26, 2008 Okay, lets have a deeper look.Your plane equation looks like( P - A ) . n = 0and your ray equation looks likeP(k) := O + k * d(where O is not a Null 0 but an Ooh ;) ).That yields in the( O + k * d - A ) . n = 0<=>k = ( ( A - O ) . n ) / ( d . n )If I made no mistake, then there is already a difference to your implementation: Your distance is negated!Your calculation of the dominant and intersection point seems me correct.Then you compute the difference vectors. For any reason you chosed to "exchange" the variables. I would not do so for clarity, but I've followed your decision here:b := C - Ac := B - AHowever, you've forgotten to relate P to A, too. IMHO here ap := P - Ais missed, isn't it?Your projection stuff seems me correct.Now you have to compute u,v from linear combinations. From the 2 possibilities, you've chosenu * c + v * b = por, as the equal equation systemu * c_{x} + v * b_{x} = p_{x}u * c_{y} + v * b_{y} = p_{y}Eliminating u,v, resp., and re-ordering terms, let me assume thatu = ( p_{x} * b_{y} - p_{y} * b_{x} ) / ( c_{x} * b_{y} - c_{y} * b_{x} )v = ( p_{x} * c_{y} - p_{y} * c_{x} ) / ( b_{x} * c_{y} - b_{y} * c_{x} )are the results. Your implementation isn't identical, but it's just an expansion of the fractions by -1/-1, so it doesn't really matter. 0 Share this post Link to post Share on other sites
skotinator 122 Report post Posted March 26, 2008 Haegar!Thankyou so much for your time explaining that through.I didn't entirely understand the purpose of calculating the line segments of the triangle edges relative to A, so I didn't think to do the same to the intersection point :)I was negating distance because I was using (O-A).N not (A-0).N. Guess thats how i saw the formula written to begin with :PAnyway, here is complete & working, hopefully clearly commented code for anyone else that happens to stumble across this looking for a full implementation of the barycentric intersection technique./**Function to return the distance distance along a ray an intersection with the triangle occurs atRay R The ray to compute intersection withfloat maxDist The maximum acceptable distance along the ray an intersection can occur atreturn float +ve the distance along the ray an intersection occurs at or -ve no intersection**/float OBJTriangleProxy::testIntersection (Ray &r, float maxDist){ //Get ray origin and direction to simplify equations later on Vector3d O = r.getOrigin(); Vector3d D = r.getDirection(); //Get triangle vertices Vector3d A = getVert (0); Vector3d B = getVert (1); Vector3d C = getVert (2); //Calculate average normal Vector3d N = (B - A).cross (C-A); N.normalise(); //Calculate distance from ray to the plane the triangle lays in float distance = ((A - O).dot (N)) / (D.dot (N)); //Ray is behind the plane the triangle's plane if (distance < 0) return -1; //Work out dominant axis for projection of triangle points into 2d int dominant = 0; if (abs(N.getX()) > abs(N.getY())) { if (abs(N.getX() > abs(N.getZ()))) dominant = 0; //x dominant else dominant = 2; //z dominant } else { if (abs(N.getY() > abs(N.getZ()))) dominant = 1; //y dominant else dominant = 2; // z dominant } //Find point where the ray intersects the plane Vector3d P = O + D * distance; //Obtain points relative to A B = B - A; C = C - A; P = P - A; //Project vectors onto dominant axis plane switch (dominant) { case 0: //Project onto x B = Vector3d (B.getY(), B.getZ(), 0); C = Vector3d (C.getY(), C.getZ(), 0); P = Vector3d (P.getY(), P.getZ(), 0); break; case 1: //Project onto Y B = Vector3d (B.getX(), B.getZ(), 0); C = Vector3d (C.getX(), C.getZ(), 0); P = Vector3d (P.getX(), P.getZ(), 0); break; case 2: //Project onto Z B = Vector3d (B.getX(), B.getY(), 0); C = Vector3d (C.getX(), C.getY(), 0); P = Vector3d (P.getX(), P.getY(), 0); break; } //Calculate 2d barycentric components float u = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY())); float v = ((C.getY() * P.getX()) - (C.getX()*P.getY())) / ((B.getX() * C.getY()) - (B.getY() * C.getX())); //Point in triangle? if (u >= 0 && v >= 0 && u + v <= 1) return distance; //No intersection return -1;}There is obviously massive potential for precalulation here that i've skipped over, but for the sake of gaining some understanding of the algorithm, it is better to keep it all in one routine to begin with.[Edited by - skotinator on March 28, 2008 12:53:17 AM] 0 Share this post Link to post Share on other sites