## 3.6 Triangle Meshes

The triangle is one of the most commonly used shapes in computer graphics; complex scenes may be modeled using millions of triangles to achieve great detail. (Figure 3.11 shows an image of a complex triangle mesh of over four million triangles.)

While a natural representation would be to have a Triangle shape implementation where each triangle stored the positions of its three vertices, a more memory-efficient representation is to separately store entire triangle meshes with an array of vertex positions where each individual triangle just stores three offsets into this array for its three vertices.

To see why this is the case, consider the celebrated Euler-Poincaré formula, which relates the number of vertices , edges , and faces on closed discrete meshes as

where is the genus of the mesh. The genus is usually a small number and can be interpreted as the number of “handles” in the mesh (analogous to a handle of a teacup). On a triangle mesh, the number of edges and vertices is furthermore related by the identity

This can be seen by dividing each edge into two parts associated with the two adjacent triangles. There are such half-edges, and all co-located pairs constitute the mesh edges. For large closed triangle meshes, the overall effect of the genus usually becomes negligible and we can combine the previous two equations (with ) to obtain

In other words, there are approximately twice as many faces as vertices. Since each face references three vertices, every vertex is (on average) referenced a total of six times. Thus, when vertices are shared, the total amortized storage required per triangle will be 12 bytes of memory for the offsets (at 4 bytes for three 32-bit integer offsets) plus half of the storage for one vertex—6 bytes, assuming three 4-byte floats are used to store the vertex position—for a total of 18 bytes per triangle. This is much better than the 36 bytes per triangle that storing the three positions directly would require. The relative storage savings are even better when there are per-vertex surface normals or texture coordinates in a mesh.

pbrt uses the TriangleMesh structure to store the shared information about a triangle mesh.

<<Triangle Declarations>>=
struct TriangleMesh { <<TriangleMesh Public Methods>>
TriangleMesh(const Transform &ObjectToWorld, int nTriangles, const int *vertexIndices, int nVertices, const Point3f *P, const Vector3f *S, const Normal3f *N, const Point2f *uv, const std::shared_ptr<Texture<Float>> &alphaMask);
<<TriangleMesh Data>>
const int nTriangles, nVertices; std::vector<int> vertexIndices; std::unique_ptr<Point3f[]> p; std::unique_ptr<Normal3f[]> n; std::unique_ptr<Vector3f[]> s; std::unique_ptr<Point2f[]> uv; std::shared_ptr<Texture<Float>> alphaMask;
};

The arguments to the TriangleMesh constructor are as follows:

• ObjectToWorld: The object-to-world transformation for the mesh.
• nTriangles: The total number of triangles in the mesh.
• vertexIndices: A pointer to an array of vertex indices. For the ith triangle, its three vertex positions are P[vertexIndices[3*i]], P[vertexIndices[3*i+1]], and P[vertexIndices[3*i+2]].
• nVertices: The total number of vertices in the mesh.
• P: An array of nVertices vertex positions.
• S: An optional array of tangent vectors, one per vertex in the mesh. These are used to compute shading tangents.
• N: An optional array of normal vectors, one per vertex in the mesh. If present, these are interpolated across triangle faces to compute shading normals.
• UV: An optional array of parametric values, one for each vertex.
• alphaMask: An optional alpha mask texture, which can be used to cut away parts of triangle surfaces.

Triangles have a dual role among the shapes in pbrt: not only are they frequently directly specified in scene description files, but other shapes often tessellate themselves into triangle meshes. For example, subdivision surfaces end up creating a mesh of triangles to approximate the smooth limit surface. Ray intersections are performed against these triangles, rather than directly against the subdivision surface (Section 3.8.2).

Due to this second role, it’s important that code that creates triangle meshes be able to specify the parameterization of the triangles. If a triangle was created by evaluating the position of a parametric surface at three particular coordinate values, for example, those values should be interpolated to compute the value at ray intersection points inside the triangle. Explicitly specified values are also useful for texture mapping, where an external program that created a triangle mesh may want to assign coordinates to the mesh so that a texture map assigns color to the mesh surface in the desired way.

The TriangleMesh constructor copies the relevant information and stores it in member variables. In particular, it makes its own copies of vertexIndices, P, N, S, and UV, allowing the caller to retain ownership of the data being passed in.

<<Triangle Method Definitions>>=
TriangleMesh::TriangleMesh(const Transform &ObjectToWorld, int nTriangles, const int *vertexIndices, int nVertices, const Point3f *P, const Vector3f *S, const Normal3f *N, const Point2f *UV, const std::shared_ptr<Texture<Float>> &alphaMask) : nTriangles(nTriangles), nVertices(nVertices), vertexIndices(vertexIndices, vertexIndices + 3 * nTriangles), alphaMask(alphaMask) { <<Transform mesh vertices to world space>>
p.reset(new Point3f[nVertices]); for (int i = 0; i < nVertices; ++i) p[i] = ObjectToWorld(P[i]);
<<Copy UV, N, and S vertex data, if present>>
if (UV) { uv.reset(new Point2f[nVertices]); memcpy(uv.get(), UV, nVertices * sizeof(Point2f)); } if (N) { n.reset(new Normal3f[nVertices]); for (int i = 0; i < nVertices; ++i) n[i] = ObjectToWorld(N[i]); } if (S) { s.reset(new Vector3f[nVertices]); for (int i = 0; i < nVertices; ++i) s[i] = ObjectToWorld(S[i]); }
}

<<TriangleMesh Data>>=
const int nTriangles, nVertices; std::vector<int> vertexIndices; std::unique_ptr<Point3f[]> p; std::unique_ptr<Normal3f[]> n; std::unique_ptr<Vector3f[]> s; std::unique_ptr<Point2f[]> uv; std::shared_ptr<Texture<Float>> alphaMask;

Unlike the other shapes that leave the shape description in object space and then transform incoming rays from world space to object space, triangle meshes transform the shape into world space and thus save the work of transforming incoming rays into object space and the work of transforming the intersection’s geometric representation out to world space. This is a good idea because this operation can be performed once at start-up, avoiding transforming rays many times during rendering. Using this approach with quadrics is more complicated, although possible—see Exercise 3.9.7 at the end of the chapter for more information.

<<Transform mesh vertices to world space>>=
p.reset(new Point3f[nVertices]); for (int i = 0; i < nVertices; ++i) p[i] = ObjectToWorld(P[i]);

The fragment <<Copy uv, N, and S vertex data, if present>> just allocates the appropriate amount of space and copies the appropriate values. Normals and tangent vectors, if present, are also transformed to world space. This fragment’s implementation isn’t included here.

### 3.6.1 Triangle

The Triangle class actually implements the Shape interface. It represents a single triangle.

<<Triangle Declarations>>+=
class Triangle : public Shape { public: <<Triangle Public Methods>>
Triangle(const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, const std::shared_ptr<TriangleMesh> &mesh, int triNumber) : Shape(ObjectToWorld, WorldToObject, reverseOrientation), mesh(mesh) { v = &mesh->vertexIndices[3 * triNumber]; } Bounds3f ObjectBound() const; Bounds3f WorldBound() const; bool Intersect(const Ray &ray, Float *tHit, SurfaceInteraction *isect, bool testAlphaTexture) const; bool IntersectP(const Ray &ray, bool testAlphaTexture) const; Float Area() const; Interaction Sample(const Point2f &u) const;
private: <<Triangle Private Methods>>
void GetUVs(Point2f uv) const { if (mesh->uv) { uv = mesh->uv[v]; uv = mesh->uv[v]; uv = mesh->uv[v]; } else { uv = Point2f(0, 0); uv = Point2f(1, 0); uv = Point2f(1, 1); } }
<<Triangle Private Data>>
std::shared_ptr<TriangleMesh> mesh; const int *v;
};

Triangle doesn’t store much data—just a pointer to the parent TriangleMesh that it came from and a pointer to its three vertex indices in the mesh.

<<Triangle Public Methods>>=
Triangle(const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, const std::shared_ptr<TriangleMesh> &mesh, int triNumber) : Shape(ObjectToWorld, WorldToObject, reverseOrientation), mesh(mesh) { v = &mesh->vertexIndices[3 * triNumber]; }

Note that the implementation stores a pointer to the first vertex index, instead of storing three pointers to the vertices themselves. This reduces the amount of storage required for each Triangle at a cost of another level of indirection.

<<Triangle Private Data>>=
std::shared_ptr<TriangleMesh> mesh; const int *v;

Because a number of other shape representations in pbrt convert themselves into triangle meshes, the utility function CreateTriangleMesh() takes care of creating an underlying TriangleMesh as well as a Triangle for each triangle in the mesh. It returns a vector of triangle shapes.

<<Triangle Method Definitions>>+=
std::vector<std::shared_ptr<Shape>> CreateTriangleMesh( const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, int nTriangles, const int *vertexIndices, int nVertices, const Point3f *p, const Vector3f *s, const Normal3f *n, const Point2f *uv, const std::shared_ptr<Texture<Float>> &alphaMask) { std::shared_ptr<TriangleMesh> mesh = std::make_shared<TriangleMesh>( *ObjectToWorld, nTriangles, vertexIndices, nVertices, p, s, n, uv, alphaMask); std::vector<std::shared_ptr<Shape>> tris; for (int i = 0; i < nTriangles; ++i) tris.push_back(std::make_shared<Triangle>(ObjectToWorld, WorldToObject, reverseOrientation, mesh, i)); return tris; }

The object space bound of a triangle is easily found by computing a bounding box that encompasses its three vertices. Because the vertex positions p are transformed to world space in the constructor, the implementation here has to transform them back to object space before computing their bound.

<<Triangle Method Definitions>>+=
Bounds3f Triangle::ObjectBound() const { <<Get triangle vertices in p0, p1, and p2>>
const Point3f &p0 = mesh->p[v]; const Point3f &p1 = mesh->p[v]; const Point3f &p2 = mesh->p[v];
return Union(Bounds3f((*WorldToObject)(p0), (*WorldToObject)(p1)), (*WorldToObject)(p2)); }

<<Get triangle vertices in p0, p1, and p2>>=
const Point3f &p0 = mesh->p[v]; const Point3f &p1 = mesh->p[v]; const Point3f &p2 = mesh->p[v];

The Triangle shape is one of the shapes that can compute a better world space bound than can be found by transforming its object space bounding box to world space. Its world space bound can be directly computed from the world space vertices.

<<Triangle Method Definitions>>+=
Bounds3f Triangle::WorldBound() const { <<Get triangle vertices in p0, p1, and p2>>
const Point3f &p0 = mesh->p[v]; const Point3f &p1 = mesh->p[v]; const Point3f &p2 = mesh->p[v];
return Union(Bounds3f(p0, p1), p2); }

### 3.6.2 Triangle Intersection

The structure of the triangle shape’s Intersect() method follows the form of earlier intersection test methods: a geometric test is applied to determine if there is an intersection and, if so, further information is computed about the intersection to return in the given SurfaceInteraction.

<<Triangle Method Definitions>>+=
bool Triangle::Intersect(const Ray &ray, Float *tHit, SurfaceInteraction *isect, bool testAlphaTexture) const { <<Get triangle vertices in p0, p1, and p2>>
const Point3f &p0 = mesh->p[v]; const Point3f &p1 = mesh->p[v]; const Point3f &p2 = mesh->p[v];
<<Perform ray–triangle intersection test>>
<<Transform triangle vertices to ray coordinate space>>
<<Translate vertices based on ray origin>>
Point3f p0t = p0 - Vector3f(ray.o); Point3f p1t = p1 - Vector3f(ray.o); Point3f p2t = p2 - Vector3f(ray.o);
<<Permute components of triangle vertices and ray direction>>
int kz = MaxDimension(Abs(ray.d)); int kx = kz + 1; if (kx == 3) kx = 0; int ky = kx + 1; if (ky == 3) ky = 0; Vector3f d = Permute(ray.d, kx, ky, kz); p0t = Permute(p0t, kx, ky, kz); p1t = Permute(p1t, kx, ky, kz); p2t = Permute(p2t, kx, ky, kz);
<<Apply shear transformation to translated vertex positions>>
Float Sx = -d.x / d.z; Float Sy = -d.y / d.z; Float Sz = 1.f / d.z; p0t.x += Sx * p0t.z; p0t.y += Sy * p0t.z; p1t.x += Sx * p1t.z; p1t.y += Sy * p1t.z; p2t.x += Sx * p2t.z; p2t.y += Sy * p2t.z;
<<Compute edge function coefficients e0, e1, and e2>>
Float e0 = p1t.x * p2t.y - p1t.y * p2t.x; Float e1 = p2t.x * p0t.y - p2t.y * p0t.x; Float e2 = p0t.x * p1t.y - p0t.y * p1t.x;
<<Fall back to double-precision test at triangle edges>>
if (sizeof(Float) == sizeof(float) && (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f)) { double p2txp1ty = (double)p2t.x * (double)p1t.y; double p2typ1tx = (double)p2t.y * (double)p1t.x; e0 = (float)(p2typ1tx - p2txp1ty); double p0txp2ty = (double)p0t.x * (double)p2t.y; double p0typ2tx = (double)p0t.y * (double)p2t.x; e1 = (float)(p0typ2tx - p0txp2ty); double p1txp0ty = (double)p1t.x * (double)p0t.y; double p1typ0tx = (double)p1t.y * (double)p0t.x; e2 = (float)(p1typ0tx - p1txp0ty); }
<<Perform triangle edge and determinant tests>>
if ((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0)) return false; Float det = e0 + e1 + e2; if (det == 0) return false;
<<Compute scaled hit distance to triangle and test against ray range>>
p0t.z *= Sz; p1t.z *= Sz; p2t.z *= Sz; Float tScaled = e0 * p0t.z + e1 * p1t.z + e2 * p2t.z; if (det < 0 && (tScaled >= 0 || tScaled < ray.tMax * det)) return false; else if (det > 0 && (tScaled <= 0 || tScaled > ray.tMax * det)) return false;
<<Compute barycentric coordinates and value for triangle intersection>>
Float invDet = 1 / det; Float b0 = e0 * invDet; Float b1 = e1 * invDet; Float b2 = e2 * invDet; Float t = tScaled * invDet;
<<Ensure that computed triangle is conservatively greater than zero>>
<<Compute term for triangle error bounds>>
Float maxZt = MaxComponent(Abs(Vector3f(p0t.z, p1t.z, p2t.z))); Float deltaZ = gamma(3) * maxZt;
<<Compute and terms for triangle error bounds>>
Float maxXt = MaxComponent(Abs(Vector3f(p0t.x, p1t.x, p2t.x))); Float maxYt = MaxComponent(Abs(Vector3f(p0t.y, p1t.y, p2t.y))); Float deltaX = gamma(5) * (maxXt + maxZt); Float deltaY = gamma(5) * (maxYt + maxZt);
<<Compute term for triangle error bounds>>
Float deltaE = 2 * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
<<Compute term for triangle error bounds and check t>>
Float maxE = MaxComponent(Abs(Vector3f(e0, e1, e2))); Float deltaT = 3 * (gamma(3) * maxE * maxZt + deltaE * maxZt + deltaZ * maxE) * std::abs(invDet); if (t <= deltaT) return false;
<<Compute triangle partial derivatives>>
Vector3f dpdu, dpdv; Point2f uv; GetUVs(uv); <<Compute deltas for triangle partial derivatives>>
Vector2f duv02 = uv - uv, duv12 = uv - uv; Vector3f dp02 = p0 - p2, dp12 = p1 - p2;
Float determinant = duv02 * duv12 - duv02 * duv12; if (determinant == 0) { <<Handle zero determinant for triangle partial derivative matrix>>
CoordinateSystem(Normalize(Cross(p2 - p0, p1 - p0)), &dpdu, &dpdv);
} else { Float invdet = 1 / determinant; dpdu = ( duv12 * dp02 - duv02 * dp12) * invdet; dpdv = (-duv12 * dp02 + duv02 * dp12) * invdet; }
<<Compute error bounds for triangle intersection>>
Float xAbsSum = (std::abs(b0 * p0.x) + std::abs(b1 * p1.x) + std::abs(b2 * p2.x)); Float yAbsSum = (std::abs(b0 * p0.y) + std::abs(b1 * p1.y) + std::abs(b2 * p2.y)); Float zAbsSum = (std::abs(b0 * p0.z) + std::abs(b1 * p1.z) + std::abs(b2 * p2.z)); Vector3f pError = gamma(7) * Vector3f(xAbsSum, yAbsSum, zAbsSum);
<<Interpolate parametric coordinates and hit point>>
Point3f pHit = b0 * p0 + b1 * p1 + b2 * p2; Point2f uvHit = b0 * uv + b1 * uv + b2 * uv;
<<Test intersection against alpha texture, if present>>
if (testAlphaTexture && mesh->alphaMask) { SurfaceInteraction isectLocal(pHit, Vector3f(0,0,0), uvHit, Vector3f(0,0,0), dpdu, dpdv, Normal3f(0,0,0), Normal3f(0,0,0), ray.time, this); if (mesh->alphaMask->Evaluate(isectLocal) == 0) return false; }
<<Fill in SurfaceInteraction from triangle hit>>
*isect = SurfaceInteraction(pHit, pError, uvHit, -ray.d, dpdu, dpdv, Normal3f(0, 0, 0), Normal3f(0, 0, 0), ray.time, this); <<Override surface normal in isect for triangle>>
isect->n = isect->shading.n = Normal3f(Normalize(Cross(dp02, dp12)));
if (mesh->n || mesh->s) { <<Initialize Triangle shading geometry>>
<<Compute shading normal ns for triangle>>
Normal3f ns; if (mesh->n) ns = Normalize(b0 * mesh->n[v] + b1 * mesh->n[v] + b2 * mesh->n[v]); else ns = isect->n;
<<Compute shading tangent ss for triangle>>
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v] + b1 * mesh->s[v] + b2 * mesh->s[v]); else ss = Normalize(isect->dpdu);
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);
<<Compute and for triangle shading geometry>>
Normal3f dndu, dndv; if (mesh->n) { <<Compute deltas for triangle partial derivatives of normal>>
Vector2f duv02 = uv - uv; Vector2f duv12 = uv - uv; Normal3f dn1 = mesh->n[v] - mesh->n[v]; Normal3f dn2 = mesh->n[v] - mesh->n[v];
Float determinant = duv02 * duv12 - duv02 * duv12; if (determinant == 0) dndu = dndv = Normal3f(0, 0, 0); else { Float invDet = 1 / determinant; dndu = ( duv12 * dn1 - duv02 * dn2) * invDet; dndv = (-duv12 * dn1 + duv02 * dn2) * invDet; } } else dndu = dndv = Normal3f(0,0,0);
} <<Ensure correct orientation of the geometric normal>>
if (mesh->n) isect->n = Faceforward(isect->n, isect->shading.n); else if (reverseOrientation ^ transformSwapsHandedness) isect->n = isect->shading.n = -isect->n;
*tHit = t; return true; }

pbrt’s ray–triangle intersection test is based on first computing an affine transformation that transforms the ray such that its origin is at in the transformed coordinate system and such that its direction is along the axis. Triangle vertices are also transformed into this coordinate system before the intersection test is performed. In the following, we’ll see that applying this coordinate system transformation simplifies the intersection test logic since, for example, the and coordinates of any intersection point must be zero. Later, in Section 3.9.3, we’ll see that this transformation makes it possible to have a watertight ray–triangle intersection algorithm, such that intersections with tricky rays like those that hit the triangle right on the edge are never incorrectly reported as misses.

<<Perform ray–triangle intersection test>>=
<<Transform triangle vertices to ray coordinate space>>
<<Translate vertices based on ray origin>>
Point3f p0t = p0 - Vector3f(ray.o); Point3f p1t = p1 - Vector3f(ray.o); Point3f p2t = p2 - Vector3f(ray.o);
<<Permute components of triangle vertices and ray direction>>
int kz = MaxDimension(Abs(ray.d)); int kx = kz + 1; if (kx == 3) kx = 0; int ky = kx + 1; if (ky == 3) ky = 0; Vector3f d = Permute(ray.d, kx, ky, kz); p0t = Permute(p0t, kx, ky, kz); p1t = Permute(p1t, kx, ky, kz); p2t = Permute(p2t, kx, ky, kz);
<<Apply shear transformation to translated vertex positions>>
Float Sx = -d.x / d.z; Float Sy = -d.y / d.z; Float Sz = 1.f / d.z; p0t.x += Sx * p0t.z; p0t.y += Sy * p0t.z; p1t.x += Sx * p1t.z; p1t.y += Sy * p1t.z; p2t.x += Sx * p2t.z; p2t.y += Sy * p2t.z;
<<Compute edge function coefficients e0, e1, and e2>>
Float e0 = p1t.x * p2t.y - p1t.y * p2t.x; Float e1 = p2t.x * p0t.y - p2t.y * p0t.x; Float e2 = p0t.x * p1t.y - p0t.y * p1t.x;
<<Fall back to double-precision test at triangle edges>>
if (sizeof(Float) == sizeof(float) && (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f)) { double p2txp1ty = (double)p2t.x * (double)p1t.y; double p2typ1tx = (double)p2t.y * (double)p1t.x; e0 = (float)(p2typ1tx - p2txp1ty); double p0txp2ty = (double)p0t.x * (double)p2t.y; double p0typ2tx = (double)p0t.y * (double)p2t.x; e1 = (float)(p0typ2tx - p0txp2ty); double p1txp0ty = (double)p1t.x * (double)p0t.y; double p1typ0tx = (double)p1t.y * (double)p0t.x; e2 = (float)(p1typ0tx - p1txp0ty); }
<<Perform triangle edge and determinant tests>>
if ((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0)) return false; Float det = e0 + e1 + e2; if (det == 0) return false;
<<Compute scaled hit distance to triangle and test against ray range>>
p0t.z *= Sz; p1t.z *= Sz; p2t.z *= Sz; Float tScaled = e0 * p0t.z + e1 * p1t.z + e2 * p2t.z; if (det < 0 && (tScaled >= 0 || tScaled < ray.tMax * det)) return false; else if (det > 0 && (tScaled <= 0 || tScaled > ray.tMax * det)) return false;
<<Compute barycentric coordinates and value for triangle intersection>>
Float invDet = 1 / det; Float b0 = e0 * invDet; Float b1 = e1 * invDet; Float b2 = e2 * invDet; Float t = tScaled * invDet;
<<Ensure that computed triangle is conservatively greater than zero>>
<<Compute term for triangle error bounds>>
Float maxZt = MaxComponent(Abs(Vector3f(p0t.z, p1t.z, p2t.z))); Float deltaZ = gamma(3) * maxZt;
<<Compute and terms for triangle error bounds>>
Float maxXt = MaxComponent(Abs(Vector3f(p0t.x, p1t.x, p2t.x))); Float maxYt = MaxComponent(Abs(Vector3f(p0t.y, p1t.y, p2t.y))); Float deltaX = gamma(5) * (maxXt + maxZt); Float deltaY = gamma(5) * (maxYt + maxZt);
<<Compute term for triangle error bounds>>
Float deltaE = 2 * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
<<Compute term for triangle error bounds and check t>>
Float maxE = MaxComponent(Abs(Vector3f(e0, e1, e2))); Float deltaT = 3 * (gamma(3) * maxE * maxZt + deltaE * maxZt + deltaZ * maxE) * std::abs(invDet); if (t <= deltaT) return false;

There are three steps to computing the transformation from world space to the ray–triangle intersection coordinate space: a translation , a coordinate permutation , and a shear . Rather than computing explicit transformation matrices for each of these and then computing an aggregate transformation matrix to transform vertices to the coordinate space, the following implementation applies each step of the transformation directly, which ends up being a more efficient approach.

<<Transform triangle vertices to ray coordinate space>>=
<<Translate vertices based on ray origin>>
Point3f p0t = p0 - Vector3f(ray.o); Point3f p1t = p1 - Vector3f(ray.o); Point3f p2t = p2 - Vector3f(ray.o);
<<Permute components of triangle vertices and ray direction>>
int kz = MaxDimension(Abs(ray.d)); int kx = kz + 1; if (kx == 3) kx = 0; int ky = kx + 1; if (ky == 3) ky = 0; Vector3f d = Permute(ray.d, kx, ky, kz); p0t = Permute(p0t, kx, ky, kz); p1t = Permute(p1t, kx, ky, kz); p2t = Permute(p2t, kx, ky, kz);
<<Apply shear transformation to translated vertex positions>>
Float Sx = -d.x / d.z; Float Sy = -d.y / d.z; Float Sz = 1.f / d.z; p0t.x += Sx * p0t.z; p0t.y += Sy * p0t.z; p1t.x += Sx * p1t.z; p1t.y += Sy * p1t.z; p2t.x += Sx * p2t.z; p2t.y += Sy * p2t.z;

The translation that places the ray origin at the origin of the coordinate system is:

This transformation doesn’t need to be explicitly applied to the ray origin, but we will apply it to the three triangle vertices.

<<Translate vertices based on ray origin>>=
Point3f p0t = p0 - Vector3f(ray.o); Point3f p1t = p1 - Vector3f(ray.o); Point3f p2t = p2 - Vector3f(ray.o);

Next, the three dimensions of the space are permuted so that the dimension is the one where the absolute value of the ray’s direction is largest. The and dimensions are arbitrarily assigned to the other two dimensions. This step ensures that if, for example, the original ray’s direction is zero, then a dimension with non-zero magnitude is mapped to .

For example, if the ray’s direction had the largest magnitude in , the permutation would be:

As before, it’s easiest to just permute the dimensions of the ray direction and the translated triangle vertices directly.

<<Permute components of triangle vertices and ray direction>>=
int kz = MaxDimension(Abs(ray.d)); int kx = kz + 1; if (kx == 3) kx = 0; int ky = kx + 1; if (ky == 3) ky = 0; Vector3f d = Permute(ray.d, kx, ky, kz); p0t = Permute(p0t, kx, ky, kz); p1t = Permute(p1t, kx, ky, kz); p2t = Permute(p2t, kx, ky, kz);

Finally, a shear transformation aligns the ray direction with the axis:

To see how this transformation works, consider its operation on the ray direction vector .

For now, only the and dimensions are sheared; we can wait and shear the dimension only if the ray actually intersects the triangle.

<<Apply shear transformation to translated vertex positions>>=
Float Sx = -d.x / d.z; Float Sy = -d.y / d.z; Float Sz = 1.f / d.z; p0t.x += Sx * p0t.z; p0t.y += Sy * p0t.z; p1t.x += Sx * p1t.z; p1t.y += Sy * p1t.z; p2t.x += Sx * p2t.z; p2t.y += Sy * p2t.z;

Note that the calculations for the coordinate permutation and the shear coefficients only depend on the given ray; they are independent of the triangle. In a high-performance ray tracer, we might want to compute these values once and store them in the Ray class, rather than recomputing them for each triangle the ray is intersected with.

With the triangle vertices transformed to this coordinate system, our task now is to find if the ray starting from the origin and traveling along the axis intersects the transformed triangle. Because of the way the coordinate system was constructed, this problem is equivalent to the 2D problem of determining if the , coordinates are inside the projection of the triangle (Figure 3.12).

To understand how the intersection algorithm works, first recall from Figure 2.5 that the length of the cross product of two vectors gives the area of the parallelogram that they define. In 2D, with vectors and , the area is

Half of this area is the area of the triangle that they define. Thus, we can see that in 2D, the area of a triangle with vertices , , and is

Figure 3.13 visualizes this geometrically.

We’ll use this expression of triangle area to define a signed edge function: given two triangle vertices and , then we can define the directed edge function as the function that gives twice the area of the triangle given by , , and a given third point :

(3.1)

(See Figure 3.14.)

The edge function gives a positive value for points to the left of the line, and negative value for points to the right. Thus, if a point has edge function values of the same sign for all three edges of a triangle, it must be on the same side of all three edges and thus must be inside the triangle.

Thanks to the coordinate system transformation, the point we’re testing is . This simplifies the edge function expressions. For example, for the edge from to , we have:

(3.2)

In the following, we’ll use the indexing scheme that the edge function corresponds to the directed edge from vertex and .

<<Compute edge function coefficients e0, e1, and e2>>=
Float e0 = p1t.x * p2t.y - p1t.y * p2t.x; Float e1 = p2t.x * p0t.y - p2t.y * p0t.x; Float e2 = p0t.x * p1t.y - p0t.y * p1t.x;

In the rare case that any of the edge function values is exactly zero, it’s not possible to be sure if the ray hits the triangle or not, and the edge equations are reevaluated using double-precision floating-point arithmetic. (Section 3.9.3 discusses the need for this step in more detail.) The fragment that implements this computation, <<Fall back to double-precision test at triangle edges>>, is just a reimplementation of <<Compute edge function coefficients e0, e1, and e2>> using doubles and so isn’t included here.

Given the values of the three edge functions, we have our first two opportunities to determine that there is no intersection. First, if the signs of the edge function values differ, then the point is not on the same side of all three edges and therefore is outside the triangle. Second, if the sum of the three edge function values is zero, then the ray is approaching the triangle edge-on, and we report no intersection. (For a closed triangle mesh, the ray will hit a neighboring triangle instead.)

<<Perform triangle edge and determinant tests>>=
if ((e0 < 0 || e1 < 0 || e2 < 0) && (e0 > 0 || e1 > 0 || e2 > 0)) return false; Float det = e0 + e1 + e2; if (det == 0) return false;

Because the ray starts at the origin, has unit length, and is along the axis, the coordinate value of the intersection point is equal to the intersection’s parametric value. To compute this value, we first need to go ahead and apply the shear transformation to the coordinates of the triangle vertices. Given these values, the barycentric coordinates of the intersection point in the triangle can be used to interpolate them across the triangle. They are given by dividing each edge function value by the sum of edge function values:

Thus, the sum to one.

The interpolated value is given by

where are the coordinates of the three vertices in the ray–triangle intersection coordinate system.

In order to save the cost of the floating-point division to compute in cases where the final value is out of the range of valid values, the implementation here first computes by interpolating with (in other words, not yet performing the division by ). If the sign of and the sign of the interpolated value are different, then the final value will certainly be negative and thus not a valid intersection.

Along similar lines, the check can be equivalent performed in two ways:

<<Compute scaled hit distance to triangle and test against ray range>>=
p0t.z *= Sz; p1t.z *= Sz; p2t.z *= Sz; Float tScaled = e0 * p0t.z + e1 * p1t.z + e2 * p2t.z; if (det < 0 && (tScaled >= 0 || tScaled < ray.tMax * det)) return false; else if (det > 0 && (tScaled <= 0 || tScaled > ray.tMax * det)) return false;

We now know that there is a valid intersection and will go ahead and pay the cost of the floating-point division to compute actual barycentric coordinates as well as the actual value for the intersection.

<<Compute barycentric coordinates and value for triangle intersection>>=
Float invDet = 1 / det; Float b0 = e0 * invDet; Float b1 = e1 * invDet; Float b2 = e2 * invDet; Float t = tScaled * invDet;

In order to generate consistent tangent vectors over triangle meshes, it is necessary to compute the partial derivatives and using the parametric values at the triangle vertices, if provided. Although the partial derivatives are the same at all points on the triangle, the implementation here recomputes them each time an intersection is found. Although this results in redundant computation, the storage savings for large triangle meshes can be significant.

A triangle can be described by the set of points

for some , where and range over the parametric coordinates of the triangle. We also know the three vertex positions , , and the texture coordinates at each vertex. From this it follows that the partial derivatives of must satisfy

In other words, there is a unique affine mapping from the 2D space to points on the triangle (such a mapping exists even though the triangle is specified in 3D space because the triangle is planar). To compute expressions for and , we start by computing the differences and , giving the matrix equation

Thus,

Inverting a 22 matrix is straightforward. The inverse of the differences matrix is

<<Compute triangle partial derivatives>>=
Vector3f dpdu, dpdv; Point2f uv; GetUVs(uv); <<Compute deltas for triangle partial derivatives>>
Vector2f duv02 = uv - uv, duv12 = uv - uv; Vector3f dp02 = p0 - p2, dp12 = p1 - p2;
Float determinant = duv02 * duv12 - duv02 * duv12; if (determinant == 0) { <<Handle zero determinant for triangle partial derivative matrix>>
CoordinateSystem(Normalize(Cross(p2 - p0, p1 - p0)), &dpdu, &dpdv);
} else { Float invdet = 1 / determinant; dpdu = ( duv12 * dp02 - duv02 * dp12) * invdet; dpdv = (-duv12 * dp02 + duv02 * dp12) * invdet; }

<<Compute deltas for triangle partial derivatives>>=
Vector2f duv02 = uv - uv, duv12 = uv - uv; Vector3f dp02 = p0 - p2, dp12 = p1 - p2;

Finally, it is necessary to handle the case when the matrix is singular and therefore cannot be inverted. Note that this only happens when the user-supplied per-vertex parameterization values are degenerate. In this case, the Triangle just chooses an arbitrary coordinate system about the triangle’s surface normal, making sure that it is orthonormal:

<<Handle zero determinant for triangle partial derivative matrix>>=
CoordinateSystem(Normalize(Cross(p2 - p0, p1 - p0)), &dpdu, &dpdv);

To compute the intersection point and the parametric coordinates at the hit point, the barycentric interpolation formula is applied to the vertex positions and the coordinates at the vertices. As we’ll see in Section 3.9.4, this gives a more precise result for the intersection point than evaluating the parametric ray equation using t.

<<Interpolate parametric coordinates and hit point>>=
Point3f pHit = b0 * p0 + b1 * p1 + b2 * p2; Point2f uvHit = b0 * uv + b1 * uv + b2 * uv;

The utility routine GetUVs() returns the coordinates for the three vertices of the triangle, either from the Triangle, if it has them, or returning default values if explicit coordinates were not specified with the mesh.

<<Triangle Private Methods>>=
void GetUVs(Point2f uv) const { if (mesh->uv) { uv = mesh->uv[v]; uv = mesh->uv[v]; uv = mesh->uv[v]; } else { uv = Point2f(0, 0); uv = Point2f(1, 0); uv = Point2f(1, 1); } }

Before a successful intersection is reported, the intersection point is tested against an alpha mask texture, if one has been assigned to the shape. This texture can be thought of as a 1D function over the triangle’s surface, where at any point where its value is zero, the intersection is ignored, effectively treating that point on the triangle as not being present. (Chapter 10 defines the texture interface and implementations in more detail.) Alpha masks can be helpful for representing objects like leaves: a leaf can be modeled as a single triangle, with an alpha mask “cutting out” the edges so that a leaf shape remains. This functionality is less often useful for other shapes, so pbrt only supports it for triangles.

<<Test intersection against alpha texture, if present>>=
if (testAlphaTexture && mesh->alphaMask) { SurfaceInteraction isectLocal(pHit, Vector3f(0,0,0), uvHit, Vector3f(0,0,0), dpdu, dpdv, Normal3f(0,0,0), Normal3f(0,0,0), ray.time, this); if (mesh->alphaMask->Evaluate(isectLocal) == 0) return false; }

Now we certainly have a valid intersection and can update the values pointed to by the pointers passed to the intersection routine. Unlike other shapes’ implementations, the code that initializes the SurfaceInteraction structure here doesn’t need to transform the partial derivatives to world space, since the triangle’s vertices were already transformed to world space. Like the disk, the partial derivatives of the triangle’s normal are also both , since it is flat.

<<Fill in SurfaceInteraction from triangle hit>>=
*isect = SurfaceInteraction(pHit, pError, uvHit, -ray.d, dpdu, dpdv, Normal3f(0, 0, 0), Normal3f(0, 0, 0), ray.time, this); <<Override surface normal in isect for triangle>>
isect->n = isect->shading.n = Normal3f(Normalize(Cross(dp02, dp12)));
if (mesh->n || mesh->s) { <<Initialize Triangle shading geometry>>
<<Compute shading normal ns for triangle>>
Normal3f ns; if (mesh->n) ns = Normalize(b0 * mesh->n[v] + b1 * mesh->n[v] + b2 * mesh->n[v]); else ns = isect->n;
<<Compute shading tangent ss for triangle>>
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v] + b1 * mesh->s[v] + b2 * mesh->s[v]); else ss = Normalize(isect->dpdu);
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);
<<Compute and for triangle shading geometry>>
Normal3f dndu, dndv; if (mesh->n) { <<Compute deltas for triangle partial derivatives of normal>>
Vector2f duv02 = uv - uv; Vector2f duv12 = uv - uv; Normal3f dn1 = mesh->n[v] - mesh->n[v]; Normal3f dn2 = mesh->n[v] - mesh->n[v];
Float determinant = duv02 * duv12 - duv02 * duv12; if (determinant == 0) dndu = dndv = Normal3f(0, 0, 0); else { Float invDet = 1 / determinant; dndu = ( duv12 * dn1 - duv02 * dn2) * invDet; dndv = (-duv12 * dn1 + duv02 * dn2) * invDet; } } else dndu = dndv = Normal3f(0,0,0);
} <<Ensure correct orientation of the geometric normal>>
if (mesh->n) isect->n = Faceforward(isect->n, isect->shading.n); else if (reverseOrientation ^ transformSwapsHandedness) isect->n = isect->shading.n = -isect->n;

The SurfaceInteraction constructor initializes the geometric normal n as the normalized cross product of dpdu and dpdv. This works well for most shapes, but in the case of triangle meshes it is preferable to rely on an initialization that does not depend on the underlying texture coordinates: it is fairly common to encounter meshes with bad parameterizations that do not preserve the orientation of the mesh, in which case the geometric normal would have an incorrect orientation.

We therefore initialize the geometric normal using the normalized cross product of the edge vectors dp02 and dp12, which results in the same normal up to a potential sign difference that depends on the exact order of triangle vertices (also known as the triangle’s winding order). 3D modeling packages generally try to ensure that triangles in a mesh have consistent winding orders, which makes this approach more robust.

<<Override surface normal in isect for triangle>>=
isect->n = isect->shading.n = Normal3f(Normalize(Cross(dp02, dp12)));

When interpolated normals are available, then we consider those to be the most authoritative source of orientation information. In this case, we flip the orientation of isect->n if the angle between it and the interpolated normal is greater than 90 degrees.

<<Ensure correct orientation of the geometric normal>>=
if (mesh->n) isect->n = Faceforward(isect->n, isect->shading.n); else if (reverseOrientation ^ transformSwapsHandedness) isect->n = isect->shading.n = -isect->n;

With Triangles, the user can provide normal vectors and tangent vectors at the vertices of the mesh that are interpolated to give normals and tangents at points on the faces of triangles. Shading geometry with interpolated normals can make otherwise faceted triangle meshes appear to be smoother than they geometrically are. If either shading normals or shading tangents have been provided, they are used to initialize the shading geometry in the SurfaceInteraction.

<<Compute shading normal ns for triangle>>
Normal3f ns; if (mesh->n) ns = Normalize(b0 * mesh->n[v] + b1 * mesh->n[v] + b2 * mesh->n[v]); else ns = isect->n;
<<Compute shading tangent ss for triangle>>
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v] + b1 * mesh->s[v] + b2 * mesh->s[v]); else ss = Normalize(isect->dpdu);
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);
<<Compute and for triangle shading geometry>>
Normal3f dndu, dndv; if (mesh->n) { <<Compute deltas for triangle partial derivatives of normal>>
Vector2f duv02 = uv - uv; Vector2f duv12 = uv - uv; Normal3f dn1 = mesh->n[v] - mesh->n[v]; Normal3f dn2 = mesh->n[v] - mesh->n[v];
Float determinant = duv02 * duv12 - duv02 * duv12; if (determinant == 0) dndu = dndv = Normal3f(0, 0, 0); else { Float invDet = 1 / determinant; dndu = ( duv12 * dn1 - duv02 * dn2) * invDet; dndv = (-duv12 * dn1 + duv02 * dn2) * invDet; } } else dndu = dndv = Normal3f(0,0,0);

Given the barycentric coordinates of the intersection point, it’s straightforward to compute the shading normal by interpolating among the appropriate vertex normals, if present.

<<Compute shading normal ns for triangle>>=
Normal3f ns; if (mesh->n) ns = Normalize(b0 * mesh->n[v] + b1 * mesh->n[v] + b2 * mesh->n[v]); else ns = isect->n;

The shading tangent is computed similarly.

<<Compute shading tangent ss for triangle>>=
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v] + b1 * mesh->s[v] + b2 * mesh->s[v]); else ss = Normalize(isect->dpdu);

The bitangent vector ts is found using the cross product of ns and ss, giving a vector orthogonal to the two of them. Next, ss is overwritten with the cross product of ts and ns; this ensures that the cross product of ss and ts gives ns. Thus, if per-vertex and values are provided and if the interpolated and values aren’t perfectly orthogonal, will be preserved and will be modified so that the coordinate system is orthogonal.