12,446,146 members (25,613 online)
alternative version

89.4K views
90 bookmarked
Posted

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

, 15 Jul 2016 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 Brazil
Senior Analyst

## You may also be interested in...

 Pro Pro

 First Prev Next
 vote of 5 Beginner Luck24-Jan-16 20:51 Beginner Luck 24-Jan-16 20:51
 good shiftwik11-Jun-14 19:03 shiftwik 11-Jun-14 19:03
 Links to other parts Pawel Gielmuda9-Apr-13 2:22 Pawel Gielmuda 9-Apr-13 2:22
 Re: Links to other parts Yks19-Jul-16 5:38 Yks 19-Jul-16 5:38
 Latest source code and library andalmeida6-Jan-12 5:14 andalmeida 6-Jan-12 5:14
 My vote of 4 Abinash Bishoyi28-Jul-11 0:06 Abinash Bishoyi 28-Jul-11 0:06
 Black region elafali26-Jul-11 0:08 elafali 26-Jul-11 0:08
 Re: Black region andalmeida26-Jul-11 2:03 andalmeida 26-Jul-11 2:03
 Re: Black region andalmeida26-Jul-11 10:05 andalmeida 26-Jul-11 10:05
 Re: Black region elafali28-Jul-11 7:46 elafali 28-Jul-11 7:46
 Re: Black region andalmeida28-Jul-11 21:25 andalmeida 28-Jul-11 21:25
 Re: Black region elafali29-Jul-11 2:07 elafali 29-Jul-11 2:07
 Re: Black region elafali9-Aug-11 15:58 elafali 9-Aug-11 15:58
 Re: Black region andalmeida9-Aug-11 21:19 andalmeida 9-Aug-11 21:19
 Hi, I'm real busy and didn't find time to look into the issue. By your problem description you have issues in the triangles adjacent lines, so for sure it is a problem with some range comparison. You can define epsilons like 1.0E-6 for example to make the comparisons in the functions like PointInTriangle and also where you do have code like A [>,<] B replace to | A - B | [>,<] epsilon. This should make difference. Another thing is that PointInTriangle returns false you can never reach a point outside the image map so it must be within [0-size-1](x,y). I'm sorry but I do not have time to review the code at the moment. If you find any improvement feel free to write me. Best Regards, Anderson Anderson J. Almeida Systems Analyst TIHunter http://www.tihunter.com
 It is simple to imagine that a line intersecting a sphere can result 0 intersections, 1 intersection (if tangent) or at most 2 intersections BooGhost10-Feb-10 5:15 BooGhost 10-Feb-10 5:15
 Re: It is simple to imagine that a line intersecting a sphere can result 0 intersections, 1 intersection (if tangent) or at most 2 intersections andalmeida10-Feb-10 5:20 andalmeida 10-Feb-10 5:20
 how to compile this code using Visual C# express Roozbeh Rock24-Sep-09 5:51 Roozbeh Rock 24-Sep-09 5:51
 How to generate images as your sample? PDU Cuong10-Mar-09 19:00 PDU Cuong 10-Mar-09 19:00
 Re: One Question andalmeida10-Aug-07 2:51 andalmeida 10-Aug-07 2:51
 Re: One Question..are you interested in some software developments? Nevi6-Feb-08 3:56 Nevi 6-Feb-08 3:56
 Great series, but.... Christian Graus3-Aug-07 13:17 Christian Graus 3-Aug-07 13:17
 Re: Great series, but.... [modified] andalmeida4-Aug-07 11:18 andalmeida 4-Aug-07 11:18
 Re: Great series, but.... andalmeida5-Aug-07 14:45 andalmeida 5-Aug-07 14:45
 Re: Great series, but.... Sander van Driel6-Aug-07 13:48 Sander van Driel 6-Aug-07 13:48
 Re: Great series, but.... Lukasz Swiatkowski8-Aug-07 11:10 Lukasz Swiatkowski 8-Aug-07 11:10
 Re: Great series, but.... Sander van Driel10-Aug-07 6:31 Sander van Driel 10-Aug-07 6:31
 Re: Great series, but.... Drew Noakes11-Aug-07 0:29 Drew Noakes 11-Aug-07 0:29
 Re: Great series, but.... colinmeier18-Feb-15 6:35 colinmeier 18-Feb-15 6:35
 Last Visit: 31-Dec-99 18:00     Last Update: 24-Aug-16 8:51 Refresh 1