Ray-triangle intersection -- how to tell if some pseudrandomly placed vertex is inside the mesh

Shawn Halayka · 2021-04-21T17:14:14

I am having difficulty with ray-triangle intersection. My full code is at: https://github.com/sjhalayka/ray_triangle_intersection The sample fractal.stl file is included with the code. When I run the code, using the example fractal.stl, I detect 2 intersections, where there should only be 1 intersection. Any ideas? The pertinent code is: /* bool rayTriangleIntersect( const vec3& orig, const vec3& dir, const vec3& v0, const vec3& v1, const vec3& v2, float& t) { // compute plane's normal vec3 v0v1 = v1 - v0; vec3 v0v2 = v2 - v0; // no need to normalize vec3 N = cross(v0v1, v0v2); // N // Step 1: finding P // check if ray and plane are parallel ? float NdotRayDirection = dot(N, dir); if (fabs(NdotRayDirection) < 0.00001) // almost 0 return false; // they are parallel so they don't intersect ! // compute d parameter using equation 2 float d = dot(N, v0); // compute t (equation 3) t = (dot(N, orig) + d) / NdotRayDirection; // check if the triangle is in behind the ray if (t < 0) return false; // the triangle is behind // compute the intersection point using equation 1 vec3 P = orig + t * dir; // Step 2: inside-outside test vec3 C; // vector perpendicular to triangle's plane // edge 0 vec3 edge0 = v1 - v0; vec3 vp0 = P - v0; C = cross(edge0, vp0); if (dot(N, C) < 0) return false; // P is on the right side // edge 1 vec3 edge1 = v2 - v1; vec3 vp1 = P - v1; C = cross(edge1, vp1); if (dot(N, C) < 0) return false; // P is on the right side // edge 2 vec3 edge2 = v0 - v2; vec3 vp2 = P - v2; C = cross(edge2, vp2); if (dot(N, C) < 0) return false; // P is on the right side; return true; // this ray hits the triangle } */ bool RayIntersectsTriangle(const vec3 rayOrigin, const vec3 rayVector, const vec3 v0, const vec3 v1, const vec3 v2, vec3& outIntersectionPoint) { const float EPSILON = 0.0000001f; vec3 vertex0 = v0;// inTriangle->vertex0; vec3 vertex1 = v1;// inTriangle->vertex1; vec3 vertex2 = v2;// inTriangle->vertex2; vec3 edge1, edge2, h, s, q; float a, f, u, v; edge1 = vertex1 - vertex0; edge2 = vertex2 - vertex0; h = cross(rayVector, edge2); a = dot(edge1, h); if (a > -EPSILON && a < EPSILON) return false; // This ray is parallel to this triangle. f = 1.0f / a; s = rayOrigin - vertex0; u = f * dot(s, h); if (u < 0.0f || u > 1.0f) return false; q = cross(s, edge1); v = f * dot(rayVector, q); if (v < 0.0f || u + v > 1.0f) return false; // At this stage we can compute t to find out where the intersection point is on the line. float t = f * dot(edge2, q); if (t > EPSILON) // ray intersection { outIntersectionPoint = rayOrigin + rayVector * t; return true; } else // This means that there is a line intersection but not a ray intersection. return false; } int main(void) { if (false == read_triangles_from_binary_stereo_lithography_file(triangles, "fractal.stl")) { cout << "Error: Could not properly read file fractal.stl" << endl; return 2; } // get_vertices_and_normals_from_triangles(triangles, face_normals, vertices, vertex_normals); vec3 pos = vec3(0, 0, 0); vec3 look_at_ray = vec3(-1, 0, 0); size_t intersection_count = 0; for (size_t i = 0; i < triangles.size(); i++) { vec3 v0 = vec3(triangles[i].vertex[0].x, triangles[i].vertex[0].y, triangles[i].vertex[0].z); vec3 v1 = vec3(triangles[i].vertex[1].x, triangles[i].vertex[1].y, triangles[i].vertex[1].z); vec3 v2 = vec3(triangles[i].vertex[2].x, triangles[i].vertex[2].y, triangles[i].vertex[2].z); vec3 outIntersectionPoint; bool found_intersection = RayIntersectsTriangle(pos, look_at_ray, v0, v1, v2, outIntersectionPoint); if (found_intersection) { intersection_count++; } } cout << intersection_count << endl; return 0; }

Math and Physics Programming C++ OpenGL

Started by taby April 15, 2021 07:41 PM

23 comments, last by JoeJ 3 years, 9 months ago

taby

Author

1,557

April 21, 2021 04:23 PM

P.S. Have you seen the stuff on trajectories?

http://paulbourke.net/fractals/trajectories/

https://vixra.org/abs/1807.0418

I'm currently working with a guy from Brazil, and we're rewriting the document, to make it way better. The guy's hard to please, but it's worth the trouble in the end. I'll let you know about the final version of the rewrite, if you're interested.

JoeJ

4,410

April 21, 2021 04:31 PM

taby said:
are you also talking about sets like Z = sin(Z) + C sin(Z)?

I don't know these. Is the sin applied to all 4 numbers of the quat here?
Looks very interesting, also the one from the other thread which looked like some crab animal.

But i assume you can not magnify a patch of surface by a factor of 10000 and new details pop up? This does not really work for Mandelbulb either. Projection trickery like twisting seems global and creates no kinds of attractive self repetitions hidden in smaller scales.
Though, in practice we may not even want such details because they generate noise, but it is expected from watching Mandelbrot zooms where precision is the only limit.

taby

Author

1,557

April 21, 2021 04:45 PM

I learned about the quaternion sin function from:

https://theworld.com/~sweetser/quaternions/intro/tools/tools.html

All I know is that if there is no small-scale for the quaternion Julia set, then the surface is differential, and it is of dimension ~2.0 for sufficiently small grid steps – like 2-sphere would. I assume that the quaternion Mandelbrot set has a non-differentiable surface, and is of dimension (2.0, 3.0).

Sound about right?