Basic Ray tracing

Started by
12 comments, last by RJSkywalker 12 years ago
Hey people, I am trying to implement a basic raytracer application. I m not using in graphics library. Just making use of the fundamental concepts.
I m having a problem while displaying my image. Its a simple utah teapot.

This is how I am implementing my ray tracing. I m not considering shadows and secondary rays yet until I get the teapot rendered.

1. Constructing the ray through the center of the pixel. Given: camera position, lookat, up vector, camera FOV. All of them are in world space.
- I m using the ray equation r = origin + t*direction
origin is my camera position. I find the 3 camera axes u,v,w:
w = lookat - camposition; u = crossproduct(w,camera->up); v = crossproduct(u,w); Normalize u, v, w
I m choosing distance to the image plane, d = -5
I m using my image plane as : left = -4, right = 4, top = 3, bottom = -3
using this i m calculating the pixelPositions: x = left + (right-left) * (i+0.5) / screenX; y = bottom + (top-bottom) * (j+0.5)/screenY
i,j is the i,j the pixel
So I get the ray direction = pixelScreenPosition - Origin
= xU + yV - d*W
I normalize this.
This completes this function

2. Before performing the intersection, i do the transformation. I m using the triangle primitives. All the vertices are in Model Space and even the Vertex Normals. So I convert them to world space.
3. Then I do the intersection of the ray with the triangle and return the closest point.
4. Then I perform the shading calculations.

I m not sure about the ray construction step. I could include the code. I hope the theory is right
Advertisement
Hello,

What is the problem you are experiencing? Can you tell us more on what the output is, show some code? Here is some information:

1. You should precalculate the inverse camera matrix (and actually, precalculate the focal plane's corners), but essentially this is correct. You can always generalize this later if you want to use more elaborate camera systems (like a spherical camera).

2. You should do that before you even start tracing anything, since this does not depend on any ray. So upon loading the teapot, you automatically transform its vertices in world space.

3. Correct, but if you want maximum speed you will need to use an acceleration structure (which can be difficult to implement), otherwise it will be quite slow.

4. Yes, and if you want to do shadow rays and more enhancements there is where you do reflections/refractions and possibly go back to step 3 with the new ray

The theory is essentially right and will give you a working raytracer, but that's only scratching the surface. Don't forget that one massive advantage of raytracers is the ability to display any analytical surface (as long as it's simple enough to be efficiently intersected with a ray), so don't just restrict yourself to triangles later on. Spheres, cylinders, planes come to mind. Of course most models are done purely in triangles.

I m not sure about the ray construction step. I could include the code. I hope the theory is right[/quote]
I can provide you with the code from my camera code from my path tracer, if you want to take a look. It's in Free Pascal but should be easily understandable.

“If I understand the standard right it is legal and safe to do this but the resulting value could be anything.”

I would really appreciate it if you could help me out and could include the code you mentioned. I even tried orthographic projection but it did not work. A bunch of random triangles get drawn instead of the teapot. I tried it the code in world space as well as in camera space but it was giving me the same result. Here is my Constructing the ray code: It is for orthographic projection

void ConstructRayThroughPixel(GzRender* render, int i, int j,GzRay *ray)
{

GzCoord u,v,w; //cam axes

for(int k = 0; k < 3; k++) //cam z-axis
w[k] = render->camera.lookat[k] - render->camera.position[k];

Normalize(&w);

crossProduct(w, render->camera.worldup, &u);
Normalize(&u);

crossProduct(u,w,&v);
Normalize(&v);

GzCoord rayDirection;

//for orthographic projection rayDir = -w
//calculate the coordinates of the pixel in world space
//Pixelscreen = Camera origin + Au + Bv
//plane coordinates in world space
float planeLeft = -4.0;
float planeRight = 4.0;
float planeTop = 3.0;
float planeBottom = -3.0;

float pixU, pixV;


pixU = planeLeft + (imagePlaneWidth * (i + 0.5) / (render->display->xres));
pixV = planeBottom + (imagePlaneHeight * (j + 0.5) / (render->display->yres));

float d = 1; //distance to the image plane from the z-axis (any value)

rayDirection[X] = d * w[X];
rayDirection[Y] = d * w[Y];
rayDirection[Z] = d * w[Z];

//normalize the direction
Normalize(&rayDirection);

GzCoord rOrigin;
//origin = e + Au + Bv
rOrigin[X] = render->camera.position[X] + (pixU * u[X]) + (pixV * v[X]);
rOrigin[Y] = render->camera.position[Y] + (pixU * u[Y]) + (pixV * v[Y]);
rOrigin[Z] = render->camera.position[Z] + (pixU * u[Z]) + (pixV * v[Z]);

//ray equation is r = rOrigin + t * direction

memcpy(&ray->origin, &rOrigin, sizeof(GzCoord));
memcpy(&ray->direction, &rayDirection, sizeof(GzCoord));
}
Ok, here is the code (this is for a perspective projection, not orthographic):

This part creates the view matrix. The matrix is 4x4, depending on how your memory is aligned you might have to switch around the entries a bit.
// This function creates an orthogonal view matrix, using a position point and a target vector.
function LookAt(Position, Target: TVector): TMatrix;
Var
Up, XAxis, YAxis, ZAxis: TVector;
begin
// Create the Up vector.
Up := Vector(0, 1, 0);

// Create the rotation axis basis.
zAxis := NormalizeVector(Target - Position);
xAxis := NormalizeVector(CrossVector(Up, zAxis));
yAxis := NormalizeVector(CrossVector(zAxis, xAxis));

// Build the matrix from those vectors.
Result[0][0] := xAxis.X;
Result[0][1] := xAxis.Y;
Result[0][2] := xAxis.Z;
Result[0][3] := -DotVector(xAxis, Position);
Result[1][0] := yAxis.X;
Result[1][1] := yAxis.Y;
Result[1][2] := yAxis.Z;
Result[1][3] := -DotVector(yAxis, Position);
Result[2][0] := zAxis.X;
Result[2][1] := zAxis.Y;
Result[2][2] := zAxis.Z;
Result[2][3] := -DotVector(zAxis, Position);
Result[3][0] := 0;
Result[3][1] := 0;
Result[3][2] := 0;
Result[3][3] := 1;
end;


This code calculates the positions of the four corners of the "virtual screen" in world space, where HorizontalFOV is your field of view in radians (probably something between pi/4 and pi/2) - the tangent function is a correcting factor, but anyway:
// Find the corner points (all the others can be interpolated with a bilinear interpolation).
FOV := Tan(HorizontalFOV * 0.5);
Corners[0] := MulMatrix(ViewMatrix, Vector(-(+FOV) * AspectRatio, (-FOV), 1));
Corners[1] := MulMatrix(ViewMatrix, Vector(-(-FOV) * AspectRatio, (-FOV), 1));
Corners[2] := MulMatrix(ViewMatrix, Vector(-(-FOV) * AspectRatio, (+FOV), 1));
Corners[3] := MulMatrix(ViewMatrix, Vector(-(+FOV) * AspectRatio, (+FOV), 1));

MulMatrix is your standard matrix multiplication (you may ignore the matrix's last row or column as it's just (0, 0, 0, 1))

And finally, this one actually traces the ray (the code above is preprocessing which is done only once when the camera moves):
// Returns the camera ray going through the pixel at U, V.
function TCamera.Trace(U, V: Double): TRay;
begin
// Interpolate the focal plane point and trace the ray to it.
Result.Origin := CameraPosition;
Result.Direction := NormalizeVector(LerpVector(LerpVector(Corners[0], Corners[1], U), LerpVector(Corners[3], Corners[2], U), V));
end;

It's basically a bilinear interpolation if you look closely. Note that U and V are assumed to be between 0 and 1 here, so if yours are between -1 and 1, just add 1 to each coordinate and halve them. LerpVector is a vector linear interpolation, e.g. lerp(a, b, t) = a + (b - a)t.

This code should give you a proper camera ray construction (in world space) which can then be intersected against your teapot.

Note that from your description "a bunch of random triangles get drawn" I don't think it's actually the camera code that's failing. Can you show your triangle intersection code? I suspect this is where the problem lies, because with camera bugs in general it's all or nothing.

“If I understand the standard right it is legal and safe to do this but the resulting value could be anything.”

Yes, i think the problem is in the intersection code. I changed it and now I can see a much clearer but still vague silhouette of the teapot. I'll post my code and the changes I made. I'll also post the image. One thing I did not understand in your corner points finding code. What does (-(+FOV) mean? I know that you are assigning the vector but I did not understand the notation of (-(+) FOV)? I have 2 versions of Intersect function. I m posting them both. My aim is to find the closest intersection point and its a triangle-ray intersection



void FindIntersection(GzRay ray, vector<SceneMesh> pScene, GzCoord* intersectionPoint,int& primitiveIndex)
{
GzCoord ClosestPoint;
int ClosestPrimitiveIdx = -1; //closest primitive not found
float closestDistance = INT_MAX; //sentinel value
ClosestPoint[X] = INT_MAX;
ClosestPoint[Y] = INT_MAX;
ClosestPoint[Z] = INT_MAX;
for(int i = 0; i < pScene.size(); i++) //loop through all the primitives (for now all the primitives are triangles that make up the teapot)
{
float tValue;
bool intersects = Intersect(ray,pScene.vertexList,tValue);

if(intersects == true)
{
//assign closest point
if(tValue < closestDistance)
{
closestDistance = tValue;
ClosestPrimitiveIdx = i;
//point is
ClosestPoint[X] = ray.origin[X] + (closestDistance * ray.direction[X]);
ClosestPoint[Y] = ray.origin[Y] + (closestDistance * ray.direction[Y]);
ClosestPoint[Z] = ray.origin[Z] + (closestDistance * ray.direction[Z]);
}
}
}
memcpy(intersectionPoint,ClosestPoint,sizeof(GzCoord));
primitiveIndex = ClosestPrimitiveIdx;
}



Here is the 1st version of Intersect function


//1st find the surface normal to the plane containing the triangle vertices
GzCoord normal;
GzCoord tempPoint;
GetDirectionNormals(pVertexList[0],pVertexList[1],pVertexList[2],&normal);
//equation of the plane containing the triangle : ax+by+cz = d
float dValue = dotProduct(normal,pVertexList[0]); // technically its -d but just check both cases
//1st check if ray.direction . N > 0, then no intersection - early elimination
float dotDirNorm = dotProduct(normal,ray.direction);

if(!(dotDirNorm < 0.0))
{
//no intersection - so early elimination
tValue = INT_MAX;
return false;
}

//now find t = (d - rayOrigin.N)/direction.N
tValue = (dValue - dotProduct(ray.origin,normal)); //don't divide yet
if(!(tValue <= 0.0)) //we won't intersection only at the front side
{
//no frontal intersection
tValue = INT_MAX;
return false;
}
tValue /= dotDirNorm;
//now that we have t, we can find the point from the ray equation
tempPoint[X] = ray.origin[X] + (tValue * ray.direction[X]);
tempPoint[Y] = ray.origin[Y] + (tValue * ray.direction[Y]);
tempPoint[Z] = ray.origin[Z] + (tValue * ray.direction[Z]);
//check if the point lies inside or outside the triangle using barycentric coordinates
if(BarycentricInterpolationTest(tempPoint,pVertexList,normal) == true)
{
return true;
}
else //point probably lies outside
{
tValue = INT_MAX;
return false;
}
}


Here is the barycentric interpolation test function:

bool BarycentricInterpolationTest(GzCoord point, GzCoord* pVertexList, GzCoord normal)
{
//we have to discard either x,y or z depending on which has higher abs value in the normal
GzCoord u,v; //edges
if(fabs(normal[X]) > fabs(normal[Y]))
{
if(fabs(normal[X]) > fabs(normal[Z]))
{
//discard x
u[0] = point[Y] - pVertexList[0][Y];
u[1] = pVertexList[1][Y] - pVertexList[0][Y];
u[2] = pVertexList[2][Y] - pVertexList[0][Y];
v[0] = point[Z] - pVertexList[0][Z];
v[1] = pVertexList[1][Z] - pVertexList[0][Z];
v[2] = pVertexList[2][Z] - pVertexList[0][Z];
}
else
{
//discard z
u[0] = point[X] - pVertexList[0][X];
u[1] = pVertexList[1][X] - pVertexList[0][X];
u[2] = pVertexList[2][X] - pVertexList[0][X];
v[0] = point[Y] - pVertexList[0][Y];
v[1] = pVertexList[1][Y] - pVertexList[0][Y];
v[2] = pVertexList[2][Y] - pVertexList[0][Y];
}
}
else
{
if(fabs(normal[Y]) > fabs(normal[Z]))
{
//discard Y
u[0] = point[X] - pVertexList[0][X];
u[1] = pVertexList[1][X] - pVertexList[0][X];
u[2] = pVertexList[2][X] - pVertexList[0][X];
v[0] = point[Z] - pVertexList[0][Z];
v[1] = pVertexList[1][Z] - pVertexList[0][Z];
v[2] = pVertexList[2][Z] - pVertexList[0][Z];
}
else
{
//discard Z
u[0] = point[X] - pVertexList[0][X];
u[1] = pVertexList[1][X] - pVertexList[0][X];
u[2] = pVertexList[2][X] - pVertexList[0][X];
v[0] = point[Y] - pVertexList[0][Y];
v[1] = pVertexList[1][Y] - pVertexList[0][Y];
v[2] = pVertexList[2][Y] - pVertexList[0][Y];
}
}
//now before computing the coordinates compute the denominator of the equation
//if it is zero return false
float denom = (u[Y] * v[Z]) - (v[Y] * u[Z]);

if(denom == 0)
{
return false;
}
//else compute the barycentric coordinates
denom = 1.0f / denom;
float alpha = ((u[X] * v[Z]) - (v[Z] * u[Z])) * denom;
if(!(alpha >= 0.0f)) //actually we might have to check if <= 1.0f
return false;

float beta = ((u[Y] * v[X]) - (v[Y] * u[X])) * denom;
if(!(beta >= 0.0f))
return false;
float gamma = 1.0f - alpha - beta;
if(!(gamma >= 0.0f))
return false;
return true;
}


Here is the 2nd and more recent version of Intersect function:

bool Intersect(GzRay ray, GzCoord* pVertexList,float& tValue)
{
//1st find the surface normal to the plane containing the triangle vertices
GzCoord normal;
GzCoord tempPoint;
GetDirectionNormals(pVertexList[0],pVertexList[1],pVertexList[2],&normal);
//equation of the plane containing the triangle : ax+by+cz = d
float dValue = dotProduct(normal,pVertexList[0]); // technically its -d but just check both cases
//1st check if ray.direction . N > 0, then no intersection - early elimination
float dotDirNorm = dotProduct(normal,ray.direction);
if(dotDirNorm != 0)
{
tValue = -(dotProduct(ray.origin, normal) + dValue) / dotDirNorm;

if(tValue > 0)
return true;
}
return false;
}


And I m also posting the images
[sharedmedia=gallery:images:2097]
[sharedmedia=gallery:images:2098]

Output2 image file I was getting using the older version of Intersect. and output1 image file I was getting using the updated version of Intersect.
Oh the weird sign stuff was just me writing down each corner mechanically so I didn't mess up the signs when I was writing this code, I didn't bother changing it after. Regular sign rules apply, i.e. -(-FOV) = FOV, +(-FOV) = -FOV, etc..

Ok. Firstly you can simplify your main intersection loop by just keeping the closest distance. Only calculate the final intersection point at the very end, when you have found the closest intersection for all triangles..Then the main intersection loop looks like this, in pseudocode:

float closest = -1
int closest_index = -1

for all primitives
{
float distance = intersection(ray, this_primitive)
if (closest < 0) or ((distance > 0) and (distance < closest))
{
closest = distance;
closest_index = this_primitive_index;
}
}

if (closest < 0)
{
// there is no intersection
}
else
{
// the intersection is ray.origin + closest * ray.direction
// and the intersected primitive is closest_index
}


Also, your triangle-ray intersection code is pretty big, there is a lot of potential for error there. Here is a version from my raytracer - it's not the fastest but it does work, it still uses barycentric coordinates but by using vectors the code is a lot shorter:

// point1, point2 and point3 are the three vertices of the triangle
// L1 is one side of the triangle (i.e. point2 - point1)
// L2 is another side of the triangle (i.e. point3 - point1)

float distance = ray.origin - point1;

vector S1 = cross(ray.direction, L2);
d = 1 / dot(S1, L1);
float u = dot(distance, S1) * d; // first barycentric coordinate

// if u < 0 or u > 1, then the ray doesn't intersect the triangle by definition of barycentric coordinates, so return a distance of -1
if (u < 0) or (u > 1) return -1;

vector S2 = cross(distance, L1);
float v = dot(ray.direction, S2) * d; // second barycentric coordinate

// if (v < 0) or (u + v > 1) then the ray doesn't intersect the triangle, by definition of barycentric coordinates
if (v < 0) or (u + v > 1) return -1;

// then the distance of the ray to the triangle is ...
return dot(L2, S2) * d:

The above code is a bit obscure I know but it does work (you can also work it out geometrically). You will need to precompute a couple of vectors for each triangle upon loading the model (L1 and L2) and store it along with the normal. Note that for now, the way you calculate them does not matter, but it will later on when triangle winding direction becomes important.

“If I understand the standard right it is legal and safe to do this but the resulting value could be anything.”

Also, your second image strikes me as being very similar to an vertex index bug. Is your model using vertex indices? It could just be you are loading your model wrong (because of how indices work, being off even a single index will garble the model and give a "spider web"-like triangle distribution)

Note debugging a raytracer isn't easy, but when it works it all fits together perfectly (from my experience).

“If I understand the standard right it is legal and safe to do this but the resulting value could be anything.”

Yes it is using vertex indices. I m loading the primitives from a .asc file that stores the vertices, and vertex normals of the triangles in Object space. I tried various intersection equations. After applying your technique, I was getting the teapot kind of right. This is what I m getting. I think now I may have to alter the camera position?

[sharedmedia=gallery:images:2099]
This looks a lot better. There are still some bugs but you are getting something. You probably want to move the camera a bit farther away from the teapot so that you can see what's going on better. I assume your teapot is centered on the origin so you can just increase the camera's distance from the teapot while keeping the target as the origin.

It seems all the triangles are actually there now, simply shaded a bit oddly. What is your shading code?

“If I understand the standard right it is legal and safe to do this but the resulting value could be anything.”

This topic is closed to new replies.

Advertisement