Archived

This topic is now archived and is closed to further replies.

quasar3d

raytracer

Recommended Posts

I made a raytracer, and as you can see there are some errors on the edges of the polygons. They are obviously floating point errors, but what can I do against them? I'm quite sure anybody who has programmed a raytracer did get this problem, so does anybody know how to get rid of these errors? My Site

Share this post


Link to post
Share on other sites
I´ve never done a raytracer myself but to avoid floating point errors there´s two way in general:

1) switch to longer floating point structures (single -> double). Thats the easiest way but it only pushes your problem away from the ".".

2) Always try to add numbers that a roughly equal in size. Say you are tracing a ray x(t) from the origin into any direction by adding dx to x each step: x(t+1) = x + dx.
Now let the values x can become be much larger than dx : x(tn) >> dx. Your floating point number only has a certain number of decimals so fractions of your dx might get lost. Assume your floating point number had 5 decimals (the real number of decimals depends on your floating point type; check your manual about it):

CORRECT COMPUTED
x(tn) = 142.34 142.34
+ dx = 0.0254 000.02
------------------------------
x(tn+1) = 142.3654 142.36

This might not look much but in fact you´re loosing 25% auf your information in dx during each step.

To solve this you might ...
a) use x(tn) = x(0) + n*dx (multiplications should work fine, but I´m not 100% sure about this one).
b) introduce a (perhaps even multiple but I guess You won´t need that) mid-level vector dx2 and some x2. You now have x(t) = x2 + dx2. Increase dx2 by dx each step. When it has reached a certain level add dx2 to x2 and clear the decimals that have been added. This way You don´t loose precious information in dx.


By taking a 2nd look at your picture there is something else that might be the reason for the "false" pixels:
1) You want dx to be small for good results.
2) Check your detection algorithm. Do you just say "hooray, I hit a plane this step" or do you check if maybe you hit more than one plane and then try to determine which is the correct one?

Like I said I never tried writing a raytracer. Also I dunno about your skill-level with coding so maybe all I said was quite obvious 2 you. Hope that helped you. If not: Try to be a bit more specific.

Share this post


Link to post
Share on other sites
use antialiasing, for example adaptive supersampling.

also, use backface culling for closed polyhedra, as in the image.

[edited by - Michbeck on July 7, 2003 8:15:05 AM]

Share this post


Link to post
Share on other sites
First of all, I am assuming you aren't using ray marching, in which case Atheist's advice is entirely useless. Antialiasing may blur out the artifacts but they will not go away. Lastly, it is not always wise to backface cull polygons in a raytracer -- you'll learn that the hard way when you go to implement refraction or shadowing.

The real problem looks like you are doing an == test on a floating point number, which is a Big No No Do Not Ever Do This Period. Floating points are rarely precisely equal to anything, so instead of (f == 0) use (f < some_very_small_constant). The very small constant is generally referred to as Epsilon. I usually use a #define or float const with a value of 0.001 for floats, or 0.00001 for doubles. You may need to tweak the epsilon a bit to get good results.

Proper epsilon usage will also cure your back face problem.


[edit] Just in case it isn't immediately obvious, if you need to test a == b, then use fabs(a - b) < epsilon, and so on.

[edited by - ApochPiQ on July 7, 2003 11:50:18 AM]

Share this post


Link to post
Share on other sites
What I do is just shooting a ray through every pixel of the image. My raytrace function returns the intersection point of the ray and the triangle. I do this for each triangle, and than use the colour of the intersection point that's the closest to the camera. The problem is that the ray at the edges intersect with more than one polygons, and the intersection points are very close, so there occur errors.
I am not doing any == test, only <

This is my raytrace function for a triangle:

bool TRIANGLE::RayTrace(const VECTOR3D &start,const VECTOR3D &end,VERTEX &intersection) const
{
float dist1 = DotProduct(start - verts[0],normal);
float dist2 = DotProduct(end - verts[0],normal);

if(!fsigncmp(dist1,dist2))
{
VECTOR3D onplane = end - start;
onplane *= (dist1 / (dist1 - dist2));
onplane += start;

int nInside = 0;
float dists[3];
for(int i = 0;i<3;i++)
{
if((dists[i] = DotProduct(onplane - verts[i],edgenormals[i])) >= 0)
nInside++;
}
if(nInside == 3)
{
// calculate the intersection point

intersection = verts[0] * (dists[1] * reciedgedistances[0]);
intersection += verts[1] * (dists[2] * reciedgedistances[1]);
intersection += verts[2] * (dists[0] * reciedgedistances[2]);
return true;
}
}
return false;
}


and the function for triangle soup:


bool RayTrace(const VECTOR3D &start,const VECTOR3D &end,VERTEX &intersection) const
{
bool hit = false;

VERTEX curintersection;
for(int i = 0;i<triangles.size();i++)
{
if(triangles[i].RayTrace(start,end,curintersection) &&
(curintersection.Len() < intersection.Len() || !hit))
{
intersection = curintersection;
hit = true;
}
}
return hit;
}


I don't want to use doubles, because my entire code base uses floats, and I don't want to rewrite everything, and it doesn't solve the problem. I also know it's possible with floats.

I know I can reduce it with antialiasing and backface culling, but they don't solve the problem, they only make it less noticable, and I want to use my raytracer for a lot of other things (like lightmaps etc), so I want it to be perfect.



My Site

[edited by - Quasar3D on July 7, 2003 12:35:12 PM]

Share this post


Link to post
Share on other sites
quote:
Original post by quasar3d
I don''t want to use doubles, because my entire code base uses floats, and I don''t want to rewrite everything


4 words: Search And Replace All.

Okay, I dont have the time to look at all your code, but I doubt the problem comes from simply using floats. First, Ive done raytracers with floats and it worked well...

Second, I mean check the dark blue pixel on the middle diagonal edge, its not just on the wrong surface, it got a totally wrong color.

I think you should look elsewhere for your problem.
Sorry if Im not of much help

Share this post


Link to post
Share on other sites
Maybe the problem is something else, but I don''t know what could be wrong. I have checked everything a billion times, and it should be correct.
that dark blue pixel on the edge is a pixel of the polygon at the back, and the lightblue pixel on the edge is on the polygon at the right side.

My Site

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
Try to get all triangle intersectors you can find on the web, and replace your with each one, and pick the one with the least errors

for starters, here''s three (sorry, forgot code tag :/ )

/* Ray-Triangle Intersection Test Routines */
/* Different optimizations of my and Ben Trumbore''s */
/* code from journals of graphics tools (JGT) */
/* http://www.acm.org/jgt/ */
/* by Tomas Moller, May 2000 */

#include <math.h>

#define EPSILON 0.000001
#define CROSS(dest,v1,v2) \
dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
#define SUB(dest,v1,v2) \
dest[0]=v1[0]-v2[0]; \
dest[1]=v1[1]-v2[1]; \
dest[2]=v1[2]-v2[2];

/* the original jgt code */
int intersect_triangle(double orig[3], double dir[3],
double vert0[3], double vert1[3], double vert2[3],
double *t, double *u, double *v)
{
double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
double det,inv_det;

/* find vectors for two edges sharing vert0 */
SUB(edge1, vert1, vert0);
SUB(edge2, vert2, vert0);

/* begin calculating determinant - also used to calculate U parameter */
CROSS(pvec, dir, edge2);

/* if determinant is near zero, ray lies in plane of triangle */
det = DOT(edge1, pvec);

if (det > -EPSILON && det < EPSILON)
return 0;
inv_det = 1.0 / det;

/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);

/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec) * inv_det;
if (*u < 0.0 || *u > 1.0)
return 0;

/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) * inv_det;
if (*v < 0.0 || *u + *v > 1.0)
return 0;

/* calculate t, ray intersects triangle */
*t = DOT(edge2, qvec) * inv_det;

return 1;
}


/* code rewritten to do tests on the sign of the determinant */
/* the division is at the end in the code */
int intersect_triangle1(double orig[3], double dir[3],
double vert0[3], double vert1[3], double vert2[3],
double *t, double *u, double *v)
{
double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
double det,inv_det;

/* find vectors for two edges sharing vert0 */
SUB(edge1, vert1, vert0);
SUB(edge2, vert2, vert0);

/* begin calculating determinant - also used to calculate U parameter */
CROSS(pvec, dir, edge2);

/* if determinant is near zero, ray lies in plane of triangle */
det = DOT(edge1, pvec);

if (det > EPSILON)
{
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);

/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u < 0.0 || *u > det)
return 0;

/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec);
if (*v < 0.0 || *u + *v > det)
return 0;

}
else if(det < -EPSILON)
{
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);

/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
/* printf("*u=%f\n",(float)*u); */
/* printf("det=%f\n",det); */
if (*u > 0.0 || *u < det)
return 0;

/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) ;
if (*v > 0.0 || *u + *v < det)
return 0;
}
else return 0; /* ray is parallell to the plane of the triangle */


inv_det = 1.0 / det;

/* calculate t, ray intersects triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;

return 1;
}

/* code rewritten to do tests on the sign of the determinant */
/* the division is before the test of the sign of the det */
int intersect_triangle2(double orig[3], double dir[3],
double vert0[3], double vert1[3], double vert2[3],
double *t, double *u, double *v)
{
double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
double det,inv_det;

/* find vectors for two edges sharing vert0 */
SUB(edge1, vert1, vert0);
SUB(edge2, vert2, vert0);

/* begin calculating determinant - also used to calculate U parameter */
CROSS(pvec, dir, edge2);

/* if determinant is near zero, ray lies in plane of triangle */
det = DOT(edge1, pvec);

/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
inv_det = 1.0 / det;

if (det > EPSILON)
{
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u < 0.0 || *u > det)
return 0;

/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec);
if (*v < 0.0 || *u + *v > det)
return 0;

}
else if(det < -EPSILON)
{
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u > 0.0 || *u < det)
return 0;

/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) ;
if (*v > 0.0 || *u + *v < det)
return 0;
}
else return 0; /* ray is parallell to the plane of the triangle */

/* calculate t, ray intersects triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;

return 1;
}

/* code rewritten to do tests on the sign of the determinant */
/* the division is before the test of the sign of the det */
/* and one CROSS has been moved out from the if-else if-else */
int intersect_triangle3(double orig[3], double dir[3],
double vert0[3], double vert1[3], double vert2[3],
double *t, double *u, double *v)
{
double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
double det,inv_det;

/* find vectors for two edges sharing vert0 */
SUB(edge1, vert1, vert0);
SUB(edge2, vert2, vert0);

/* begin calculating determinant - also used to calculate U parameter */
CROSS(pvec, dir, edge2);

/* if determinant is near zero, ray lies in plane of triangle */
det = DOT(edge1, pvec);

/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
inv_det = 1.0 / det;

CROSS(qvec, tvec, edge1);

if (det > EPSILON)
{
*u = DOT(tvec, pvec);
if (*u < 0.0 || *u > det)
return 0;

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec);
if (*v < 0.0 || *u + *v > det)
return 0;

}
else if(det < -EPSILON)
{
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u > 0.0 || *u < det)
return 0;

/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) ;
if (*v > 0.0 || *u + *v < det)
return 0;
}
else return 0; /* ray is parallell to the plane of the triangle */

*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;

return 1;
}


Share this post


Link to post
Share on other sites
I solved it. The problem was in this part:

// calculate the intersection point
intersection = verts[0] * (dists[1] * reciedgedistances[0]);
intersection += verts[1] * (dists[2] * reciedgedistances[1]);
intersection += verts[2] * (dists[0] * reciedgedistances[2]);

now I do this.

// calculate the intersection point
intersection = verts[0] * (dists[1] * reciedgedistances[0]);
intersection += verts[1] * (dists[2] * reciedgedistances[1]);
intersection += verts[2] * (dists[0] * reciedgedistances[2]);
intersection.x = onplane.x;
intersection.y = onplane.y;
intersection.z = onplane.z;

The first version calculates the intersection point the same way as it calculates the vertex colour, for it calculates the barycentric coordinates, and than combines the coordinates of three triangle vertices.
the first version would indead work with a perfect floating point format, but it's obviously very sensitive for errors.



There are still some errors. Do you think they can be solved without too much trouble, or should I just live with it.

My Site

[edited by - Quasar3D on July 7, 2003 3:56:01 PM]

Share this post


Link to post
Share on other sites
Rats, I don''t get to blame evil floating point inaccuracies

The errors you have left seem to be unavoidable. In my own tracer, when I have closely overlapping edges like that, they occasionally stick out an extra pixel. I haven''t found a way around it so far. Perhaps it can be fixed by using doubles, but they break my memory manager so I can''t go test it out for you, as I don''t have time to redo all the memory code.

You may also get better results by increasing the image resolution.

Share this post


Link to post
Share on other sites
Biggest problems seem to be at the polygon edges, so probably the error is in the point-in-poly testing routine.

Some people have a weird tendency of optimizing point-in-poly so much that it can''t handle roundoff error anymore So remember to make it work correct at first, then make it work faster (but still correct).

- Mikko Kauppila

Share this post


Link to post
Share on other sites
Today I came up with an idea on how to solve it. I haven't implemented it, and I'm also not sure if I'll ever do it, because I need to change the entire architecture of my mesh data etc for it.

How I do it now:
For the point in polygon test, I first project the point onto the plane. Than I test if that projected point is in front of the three planes on the edges of the polygon. One tangent vector of these planes is the normal of the polygon, and the other tangent is the edge of the polygon. So the normal of the plane for the first edge is:

edgenormal[0] = CrossProduct(vert[1] - vert[0],normal);

Due to the float errors the point can be on the back side of both planes on the same edge.

My solution
What I was thinking was to use one plane for an edge that's shared by two polygon. To do this I use the average of the normals both polygons sharing that edge in the statement above.

No matter how inaccurate your floating points are, one point will never lay both on the front side and the backside of one and the same plane So the problem is solved. Well, almost, since you still have to project on the 2 different planes of the polygons

My Site


[edited by - Quasar3D on July 9, 2003 9:07:03 PM]

Share this post


Link to post
Share on other sites
quote:
Original post by quasar3d
No matter how inaccurate your floating points are, one point will never lay both on the front side and the backside of one and the same plane



True, but it could still lay on the "wrong" side due to the inaccuracy of floats.

Somehow I feel that backface culling is the best solution to this, although that means that more care has to be taken when doing refraction.

Share this post


Link to post
Share on other sites
Yes but if it's on the wrong side, it will be on the wrong side too for the same point on the neighbour polygon, and for that pol the point needs to be on the other side to.

The edge at the eft can indeed be solved with backface culling, but that sort of t-junction in the middle edge can not. That's what this is for

My Site

[edited by - Quasar3D on July 10, 2003 6:48:25 AM]

Share this post


Link to post
Share on other sites
Backface culling in a raytracer is evil. You actually need MORE precision when tracing shadow/reflection/refraction rays, because they are usually originating a lot closer to the surface, so every tiny miniscule error is magnified. Besides, backface culling, as I mentioned, wreaks total havoc with recursive rays. Remember, raytracers are not like rasterizers -- they do not simply project geometry onto a view plane. They actually have to deal with a full 3D data set, not just a projection of polygons. Because of this, removing pieces of data will only cause more problems. First-hit testing in a raytracer is the functional equivalent of backface culling in a rasterizer. Each belongs in one type of renderer, and it should not be carried over to the other, or nasty problems will result. And yes, I speak from experience on this one.


Are you using quads or pairs of triangles? If you are using orthogonal planar quads (i.e. rectangles) then any errors will be resulting from floating point inaccuracies, because a rectangle point-in-poly test is just an if(x < max && x > min) series of tests. However, point-in-arbitrary-quad or point-in-triangle tests are, as mentioned earlier, often optimized into inaccuracy, so be sure you are using a very well-known and well-tested polygon collision test. I use Plucker coordinates to great effect in my raytracer.

Share this post


Link to post
Share on other sites