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

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.

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]

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]

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

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

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;
}

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]