# Simple Ray Tracing in C# Part II (Triangles Intersection)

, 9 Apr 2013 GPL3
 Rate this:
Ray tracing and mapping in spheres
 Cube (8 vertices, 12 triangles) Ambient, diffuse and specular lights

## Introduction

Previously, we have seen ray tracing and mapping in spheres. Now, I am extending the algorithm to perform ray x triangles intersection calculation. This first approach is very basic, where the face normals are not interpolated and there is no texture mapping which will be shown in next articles.

## Background

The ray tracing technique consists in calculating each ray (rn) from the observer to the projection plane (screen) and from the light to the nearest intersection (if found) of rn to the objects within the viewer. After finding the intersection point, the ray can be reflected and/or refracted by the object depending on its material, generating another path to be computated.

Assuming the projection plane as the screen we have a rectangular area in integer bounding coordinates (x0, y0,x1,y1). A point (x,y) inside this area must be mapped to a real virtual screen with global bounding coordinates in floating point (rx0,rx1,ry0,ry1) obtaining (rx, ry) the point which will be used to define the ray in R3 assuming the screen iz Z=0 we have (rx, ry, 0).

### The Mapping Function

```public static double GetCoord(double i1, double i2, double w1, double w2,
double p)
{
return ((p - i1) / (i2 - i1)) * (w2 - w1) + w1;
}```

Since now we have two R3 points (Eye position and Screen pixel represented in virtual coordinates) we can obtain the ray representation:

3D lines equations can be represented in the form:

• x = px + t*vx
• y = py + t*vy
• z = pz + t*vz

Where (px,py,pz) are all the points lying in the 3D line, t is a scalar parameter, and (vx,vy,vz) is a direction vector. The above equation can be obtained by the definition where a line can be defined by 2 points, so given P1(x1,y1,z1) and P2(x2,y2,z2) we have:

• v = (x2-x1,y2-y1,z2-z1)

Replacing the found v and p, we have:

• x = x1 + t*(x2-x1)
• y = y1 + t*(y2-y1)
• z = z1 + t*(z2-z1)

So we can be sure that all (x,y,z) which satisfies the above equation belongs to the line defined by P1P2.

## Line x Sphere Intersection

Spheres can be represented in the form:

• r2 = (x-cx)2+(y-cy)2+(z-cz)2

...where:

• r is the sphere radius
• (cx,cy,cz) is the center of the sphere

So we can be sure that all x,y,x points lies on the sphere surface. Our objective now is to determine the intersection equation between a given line and a sphere it must be a set of (x,y,z) points which satisfies both equations. It is simple to imagine that a line intersecting a sphere can result 0 intersections, 1 intersection (if tangent) or at most 2 intersections.

### The Equations

• r2 = (x-cx)2+(y-cy)2+(z-cz)2
• x = x1 + t*(x2-x1)
• y = y1 + t*(y2-y1)
• z = z1 + t*(z2-z1)

Replacing x, y and z, we have:

• r2 = (x1 + t*(x2-x1)-cx)2 + (y1 + t*(y2-y1) -cy)2 + (z1+ t*(z2-z1)-cz)2

Let's create a variable for the vector:

• vx = x2 - x1
• vy = y2 - y1
• vz = z2 - z1

So now we have:

• r2 = (x1-cx+t*vx)2 + (y1-cy+t*vy)2 + (z1-cz+t*vz)2
• (x1-cx+t*vx)2 + (y1-cy+t*vy)2 + (z1-cz+t*vz)2 - r2 = 0

Let's replace (x1,y1,z1) with (px,py,pz) just to simplify... Now we have a perfect 2nd degree equation which can give us 0, 1 or 2 different solutions for 't' :

```double A = (vx * vx + vy * vy + vz * vz);
double B = 2.0 * (px * vx + py * vy + pz * vz - vx * cx - vy * cy - vz * cz);
double C = px * px - 2 * px * cx + cx * cx + py * py - 2 * py * cy + cy * cy +
pz * pz - 2 * pz * cz    + cz * cz - radius * radius;
double D = B * B - 4 * A * C;
double t = -1.0;
if (D >= 0)
{
double t1 = (-B - System.Math.Sqrt(D)) / (2.0 * A);
double t2 = (-B + System.Math.Sqrt(D)) / (2.0 * A);
if (t1 > t2)
t = t1;
else
t = t2;  // we choose the nearest t from the first point
}```

## Line x Triangle Intersection

### Plane Equation

• N dot (P - P3) = 0 (Eq1)

Where Line (P1-P2) Equation:

• P and P3 are on the plane
• N is the plane Normal vector
• P = P1 + v (P2 - P1) (Eq2)

So the intersection is given by:

• Eq1 = Eq2 or
• N dot (P1 + v (P2 - P1)) = N dot P3

Now the intersection between line and plane is found, so we need to verify if the given intersection lies inside the triangle:

```public bool InternalSide(double p1x, double p1y, double p1z,
double p2x, double p2y, double p2z,
double ax, double ay, double az,
double bx, double by, double bz)
{
double cp1x = 0, cp1y = 0, cp1z = 0, cp2x = 0, cp2y = 0,
cp2z = 0;
tAlgebra.Cross3(bx - ax, by - ay, bz - az, p1x - ax,
p1y - ay, p1z - az, ref cp1x,
ref cp1y, ref cp1z);
tAlgebra.Cross3(bx - ax, by - ay, bz - az, p2x - ax, p2y - ay,
p2z - az, ref cp2x, ref cp2y, ref cp2z);

if (tAlgebra.Dot3(cp1x, cp1y, cp1z, cp2x, cp2y, cp2z) >= 0)
return true;
else
return false;
}

public bool PointInTriangle(double px, double py, double pz)
{
if (InternalSide(px, py, pz,
tp1x, tp1y, tp1z,
tp2x, tp2y, tp2z,
tp3x, tp3y, tp3z) &&
InternalSide(px, py, pz,
tp2x, tp2y, tp2z,
tp1x, tp1y, tp1z,
tp3x, tp3y, tp3z) &&
InternalSide(px, py, pz,
tp3x, tp3y, tp3z,
tp1x, tp1y, tp1z,
tp2x, tp2y, tp2z))
return true;
else
return false;
}```

## The Code

### Cube.obj file

• v = vertex
• f = face
```#
# cube.obj
#

v -1.000000 -1.000000  1.000000
v -1.000000  1.000000  1.000000
v  1.000000  1.000000  1.000000
v  1.000000 -1.000000  1.000000
v -1.000000 -1.000000 -1.000000
v -1.000000  1.000000 -1.000000
v  1.000000  1.000000 -1.000000
v  1.000000 -1.000000 -1.000000

f 5 8 7
f 7 6 5
f 1 4 3
f 3 2 1
f 5 1 2
f 2 6 5
f 8 4 3
f 3 7 8
f 8 4 1
f 1 5 8
f 7 3 2
f 2 6 7```

### Algebra Class

```public class tAlgebra
{
public tAlgebra()
{
}

//i' = i - (2 * n * dot(i, n))
public static void Reflect(double ix, double iy, double iz,
double nx, double ny, double nz,
ref double iix, ref double iiy, ref double iiz)
{
iix = ix - (2.0 * nx * Dot3(nx, ny, nz, ix, iy, iz));
iiy = iy - (2.0 * ny * Dot3(nx, ny, nz, ix, iy, iz));
iiz = iz - (2.0 * nz * Dot3(nx, ny, nz, ix, iy, iz));
Normalize(ref iix, ref iiy, ref iiz);
}
public static void Cross3(double ux, double uy, double uz,
double wx,
double wy, double wz,
ref double vx, ref double vy,
ref double vz)
{
// u x w
vx = wz * uy - wy * uz;
vy = wx * uz - wz * ux;
vz = wy * ux - wx * uy;
}
public static double Dot3(double x1, double y1, double z1,
double x2, double y2, double z2)
{
return (x1 * x2) + (y1 * y2) + (z1 * z2);
}
public static double GetCosAngleV1V2(double v1x, double v1y,
double v1z, double v2x, double v2y, double v2z)
{
/* incident angle
// inters pt (i)
double ix, iy, iz;
ix = px+t*vx;
iy = py+t*vy;
iz = pz+t*vz;

// normal at i
double nx, ny, nz;
nx = ix - cx;
ny = iy - cy;
nz = iz - cz;
*/
double x, y, z;
x = v1x; y = v1y; z = v1z;
Normalize(ref x, ref y, ref z);
v1x = x; v1y = y; v1z = z;

x = v2x; y = v2y; z = v2z;
Normalize(ref x, ref y, ref z);
v2x = x; v2y = y; v2z = z;

// cos(t) = (v.w) / (|v|.|w|)
double n = (v1x * v2x + v1y * v2y + v1z * v2z) ;
double d = (modv(v1x, v1y, v1z) * modv(v2x, v2y, v2z));

if(Math.Abs(d)<1.0E-10) return 0;
return n/d ;
}
public static double modv(double vx, double vy, double vz)
{
return System.Math.Sqrt(vx * vx + vy * vy + vz * vz);
}
public static double GetCoord(double i1, double i2, double w1,
double w2, double p)
{
return ((p - i1) / (i2 - i1)) * (w2 - w1) + w1;
}

public static void Normalize(ref double vx, ref double vy,
ref double vz)
{
double mod_v = tAlgebra.modv(vx, vy, vz);
if (Math.Abs(mod_v) < 1.0E-10) return;
vx = vx / mod_v;
vy = vy / mod_v;
vz = vz / mod_v;
}
}```

### Object Classes

```public class tObject
{
public tObject()
{
}
public double ambientR, ambientG, ambientB;
public double diffuseR, diffuseG, diffuseB;
public double specularR, specularG, specularB;
public double shininess;
}
public class tTriangle : tObject
{
public tTriangle()
{
}

public void Init()
{
getNormal(ref tnormalX, ref tnormalY, ref tnormalZ);
}
public bool SameSide(double p1x, double p1y, double p1z,
double p2x, double p2y, double p2z,
double ax, double ay, double az,
double bx, double by, double bz)
{
double cp1x = 0, cp1y = 0, cp1z = 0, cp2x = 0, cp2y = 0,
cp2z = 0;
tAlgebra.Cross3(bx - ax, by - ay, bz - az, p1x - ax,
p1y - ay, p1z - az, ref cp1x, ref cp1y, ref cp1z);
tAlgebra.Cross3(bx - ax, by - ay, bz - az, p2x - ax,
p2y - ay, p2z - az, ref cp2x, ref cp2y, ref cp2z);
if (tAlgebra.Dot3(cp1x, cp1y, cp1z, cp2x, cp2y, cp2z) >= 0)
return true;
else
return false;
}
public bool PointInTriangle(double px, double py, double pz)
{
if (SameSide(px, py, pz, tp1x, tp1y, tp1z,
tp2x, tp2y, tp2z, tp3x, tp3y, tp3z) &&
SameSide(px, py, pz, tp2x, tp2y, tp2z,
tp1x, tp1y, tp1z, tp3x, tp3y, tp3z) &&
SameSide(px, py, pz, tp3x, tp3y, tp3z,
tp1x, tp1y, tp1z, tp2x, tp2y, tp2z))
return true;
else
return false;
}

// ray p1, ray p2
public double GetInterSect(double p1x, double p1y, double p1z,
double p2x, double p2y, double p2z)
{
double v1x = tp3x - p1x; double v1y = tp3y - p1y;
double v1z = tp3z - p1z;
double v2x = p2x - p1x;  double v2y = p2y - p1y;
double v2z = p2z - p1z;
double dot1 = tAlgebra.Dot3(tnormalX, tnormalY, tnormalZ,
v1x, v1y, v1z);
double dot2 = tAlgebra.Dot3(tnormalX, tnormalY, tnormalZ,
v2x, v2y, v2z);
if (Math.Abs(dot2) < 1.0E-6)
return -1; // division by 0 means parallel
double u = dot1 / dot2;
// point in triangle?
if(!PointInTriangle(p1x+u*(p2x-p1x),
p1y+u*(p2y-p1y), p1z+u*(p2z-p1z)))
return -1;
return u;
}
public void getNormal(ref double vx, ref double vy, ref double vz)
{
double ux = tp3x - tp1x, uy = tp3y - tp1y, uz = tp3z - tp1z;
double wx = tp2x - tp1x, wy = tp2y - tp1y, wz = tp2z - tp1z;

// u x w
vx = wz * uy - wy * uz;
vy = wx * uz - wz * ux;
vz = wy * ux - wx * uy;
}

public double tp1x, tp1y, tp1z;
public double tp2x, tp2y, tp2z;
public double tp3x, tp3y, tp3z;
public double tnormalX, tnormalY, tnormalZ;
}```

### The Program

```int scz = 200;

System.IO.FileStream file = new System.IO.FileStream(
Server.MapPath("./data/cube.obj"),
System.IO.FileMode.Open,
file.Close();

string[] lines = filebody.Split('\n');
ArrayList pointArrayX = new ArrayList();
ArrayList pointArrayY = new ArrayList();
ArrayList pointArrayZ = new ArrayList();

System.Collections.ArrayList obj3dArrayList;
obj3dArrayList = new System.Collections.ArrayList();

double rx = 0.5000;
double ry = 0.5000;
double rz = 0.0000;

foreach (string line in lines)
{
if (line.Length < 1) continue;
string auxline = line;
if (auxline.IndexOf("mtllib") >= 0)
{
auxline = auxline.Replace("mtllib", "");
GetMaterial(auxline);
}
else
if (auxline[0] == 'v')
{
auxline = auxline.Replace("v ", "");
string[] points = auxline.Split(' ');
double x = double.Parse(points[0]);
double y = double.Parse(points[1]);
double z = double.Parse(points[2]);
}
else
if (auxline[0] == 'f')
{
auxline = auxline.Replace("f ", "");
string[] vertices = auxline.Split(' ');
int a = int.Parse(vertices[0])-1;
int b = int.Parse(vertices[1])-1;
int c = int.Parse(vertices[2])-1;

tTriangle tri1 = new tTriangle();
tri1.tp1x = (double) pointArrayX[a];
tri1.tp1y = (double) pointArrayY[a];
tri1.tp1z = (double) pointArrayZ[a];

tri1.tp2x = (double) pointArrayX;
tri1.tp2y = (double) pointArrayY;
tri1.tp2z = (double) pointArrayZ;
tri1.tp3x = (double) pointArrayX[c];
tri1.tp3y = (double) pointArrayY[c];
tri1.tp3z = (double) pointArrayZ[c];

Algebra.RotX(rx, ref tri1.tp1y, ref tri1.tp1z);
Algebra.RotX(rx, ref tri1.tp2y, ref tri1.tp2z);
Algebra.RotX(rx, ref tri1.tp3y, ref tri1.tp3z);
Algebra.RotY(ry, ref tri1.tp1x, ref tri1.tp1z);
Algebra.RotY(ry, ref tri1.tp2x, ref tri1.tp2z);
Algebra.RotY(ry, ref tri1.tp3x, ref tri1.tp3z);
Algebra.RotZ(rz, ref tri1.tp1x, ref tri1.tp1y);
Algebra.RotZ(rz, ref tri1.tp2x, ref tri1.tp2y);
Algebra.RotZ(rz, ref tri1.tp3x, ref tri1.tp3y);

// ambient properties for the material
tri1.ambientR = 0.3;
tri1.ambientG = 0.4;
tri1.ambientB = 0.5;
// specular properties for the material
tri1.specularR = 0.911;
tri1.specularG = 0.911;
tri1.specularB = 0.911;
tri1.shininess = 10.6;
tri1.diffuseR = 0.3775;
tri1.diffuseG = 0.3775;
tri1.diffuseB = 0.5775;
tri1.Init();
}
}

Bitmap newBitmap = new Bitmap(scz, scz,
PixelFormat.Format32bppArgb);
Graphics g = Graphics.FromImage(newBitmap);

Color clrBackground = Color.Black;

Rectangle rect = new Rectangle(0, 0, scz, scz);
g.FillRectangle(new SolidBrush(clrBackground),
rect);

Graphics graphics = g;

double px = 0.0;
double py = 0.0;
double pz = 600.0;

double lpx = 1.5;
double lpy = 1.7;
double lpz = 500.0;

double lp2x = lpx;
double lp2y = lpy;
double lp2z = lpz-0.1;

double lvx = lp2x - lpx;
double lvy = lp2y - lpy;
double lvz = lp2z - lpz;

tAlgebra.Normalize(ref lvx, ref lvy, ref lvz);

// virtual mapping size
double fMax = 2.5;

for (int i = rect.Left; i <= rect.Right; i++)
{
double x = tAlgebra.GetCoord(rect.Left,
rect.Right, -fMax, fMax, i);

for (int j = rect.Top; j <= rect.Bottom; j++)
{
double y = tAlgebra.GetCoord(rect.Top,
rect.Bottom, fMax, -fMax, j);

double t = 1.0E10;

double vx = x - px, vy = y - py, vz = -pz;

double mod_v = tAlgebra.modv(vx, vy, vz);
vx = vx / mod_v;
vy = vy / mod_v;
vz = vz / mod_v;

tTriangle triangleHit = null;

for (int k = 0; k < (int)obj3dArrayList.Count; k++)
{
tTriangle triN = (tTriangle)obj3dArrayList[k];
double taux = triN.GetInterSect(px, py, pz, x, y, 0.0);
if (taux < 0) continue;

if (taux > 0 && taux < t)
{
t = taux;
triangleHit = triN;
}
}

Color color = Color.FromArgb(10, 20, 10);

if (triangleHit != null)
{
double intersx = px + t * vx,
intersy = py + t * vy, intersz = pz + t * vz;
double l2px = intersx - lpx,
l2py = intersy - lpy, l2pz = intersz - lpz;
tAlgebra.Normalize(ref l2px, ref l2py, ref l2pz);

double vNormalX = triangleHit.tnormalX,
vNormalY = triangleHit.tnormalY,
vNormalZ = triangleHit.tnormalZ;
tAlgebra.Normalize(ref vNormalX, ref vNormalY,
ref vNormalZ);

double cost = tAlgebra.GetCosAngleV1V2(l2px, l2py, l2pz,
vNormalX, vNormalY, vNormalZ);
if (cost < 0) cost = 0;

double vReflX = 0, vReflY = 0, vReflZ = 0;
double vEye2IntersX = intersx - px,
vEye2IntersY = intersy - py,
vEye2IntersZ = intersz - pz;
tAlgebra.Reflect(l2px, l2py, l2pz, vNormalX, vNormalY,
vNormalZ, ref vReflX, ref vReflY, ref vReflZ);

tAlgebra.Normalize(ref vReflX, ref  vReflY, ref vReflZ);
tAlgebra.Normalize(ref vEye2IntersX, ref vEye2IntersY,
ref vEye2IntersZ);
double cosf = tAlgebra.GetCosAngleV1V2(vReflX, vReflY,
vReflZ, vEye2IntersX, vEye2IntersY,
vEye2IntersZ);

if (cosf < 0)
cosf = 0;

double result1 = cost * 255.0;
double result2 = Math.Pow(cosf, triangleHit.shininess) *
255.0;

double rgbR = (triangleHit.ambientR * 255.0) +
(triangleHit.diffuseR * result1) +
(triangleHit.specularR * result2);
double rgbG = (triangleHit.ambientG * 255.0) +
(triangleHit.diffuseG * result1) +
(triangleHit.specularG * result2);
double rgbB = (triangleHit.ambientB * 255.0) +
(triangleHit.diffuseB * result1) +
(triangleHit.specularB * result2);

rgbR = Math.Min(rgbR, 255);
rgbG = Math.Min(rgbG, 255);
rgbB = Math.Min(rgbB, 255);
rgbR = Math.Max(rgbR, 0);
rgbG = Math.Max(rgbG, 0);
rgbB = Math.Max(rgbB, 0);

color = Color.FromArgb((int)rgbR, (int)rgbG, (int)rgbB);
}

Brush brs = new SolidBrush(color);
graphics.FillRectangle(brs, i, j, 1, 1);
brs.Dispose();

}// for pixels lines
}// for pixels columns
///////////////////////////////////////

MemoryStream tempStream = new MemoryStream();
newBitmap.Save(tempStream, ImageFormat.Png);
Response.ClearContent();
Response.ContentType = "image/png";
Response.BinaryWrite(tempStream.ToArray());
Response.Flush();```

## Share

Engineer IBM
Brazil
Senior Analyst

Founder of TIHunter Vagas de TI

 First PrevNext
 good shiftwik 11-Jun-14 20:03
 Links to other parts Pawel Gielmuda 9-Apr-13 3:22
 Latest source code and library andalmeida 6-Jan-12 6:14
 My vote of 4 Abinash Bishoyi 28-Jul-11 1:06
 Black region elafali 26-Jul-11 1:08
 Re: Black region andalmeida 26-Jul-11 3:03
 Re: Black region andalmeida 26-Jul-11 11:05
 Re: Black region elafali 28-Jul-11 8:46
 Re: Black region andalmeida 28-Jul-11 22:25
 Re: Black region elafali 29-Jul-11 3:07
 Re: Black region elafali 9-Aug-11 16:58
 Re: Black region andalmeida 9-Aug-11 22:19
 how to compile this code using Visual C# express Roozbeh Rock 24-Sep-09 6:51
 How to generate images as your sample? PDU Cuong 10-Mar-09 20:00
 Re: One Question andalmeida 10-Aug-07 3:51
 Great series, but.... Christian Graus 3-Aug-07 14:17
 Re: Great series, but.... [modified] andalmeida 4-Aug-07 12:18
 Re: Great series, but.... andalmeida 5-Aug-07 15:45
 Re: Great series, but.... Sander van Driel 6-Aug-07 14:48
 Re: Great series, but.... Lukasz Swiatkowski 8-Aug-07 12:10
 Re: Great series, but.... Sander van Driel 10-Aug-07 7:31
 Both; in this context it is the same.
 Last Visit: 31-Dec-99 19:00     Last Update: 20-Dec-14 1:19 Refresh 12 Next »