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.)

Figure 3.11: Ganesha Model. This triangle mesh contains over four million individual triangles. It was created from a real statue using a 3D scanner that uses structured light to determine shapes of objects.

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 upper V , edges upper E , and faces upper F on closed discrete meshes as

upper V minus upper E plus upper F equals 2 left-parenthesis 1 minus g right-parenthesis comma

where g element-of bold upper N Superscript 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

upper E equals three-halves upper F period

This can be seen by dividing each edge into two parts associated with the two adjacent triangles. There are 3 upper F such half-edges, and all co-located pairs constitute the upper E 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 g equals 0 ) to obtain

upper F almost-equals 2 upper V period

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 left-parenthesis u comma v right-parenthesis 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 left-parenthesis u comma v right-parenthesis coordinate values, for example, those left-parenthesis u comma v right-parenthesis values should be interpolated to compute the left-parenthesis u comma v right-parenthesis value at ray intersection points inside the triangle. Explicitly specified left-parenthesis u comma v right-parenthesis values are also useful for texture mapping, where an external program that created a triangle mesh may want to assign left-parenthesis u comma v right-parenthesis 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[3]) const { if (mesh->uv) { uv[0] = mesh->uv[v[0]]; uv[1] = mesh->uv[v[1]]; uv[2] = mesh->uv[v[2]]; } else { uv[0] = Point2f(0, 0); uv[1] = Point2f(1, 0); uv[2] = 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[0]]; const Point3f &p1 = mesh->p[v[1]]; const Point3f &p2 = mesh->p[v[2]];
return Union(Bounds3f((*WorldToObject)(p0), (*WorldToObject)(p1)), (*WorldToObject)(p2)); }

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

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[0]]; const Point3f &p1 = mesh->p[v[1]]; const Point3f &p2 = mesh->p[v[2]];
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[0]]; const Point3f &p1 = mesh->p[v[1]]; const Point3f &p2 = mesh->p[v[2]];
<<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 t 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 t 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 t is conservatively greater than zero>> 
<<Compute delta Subscript z term for triangle t error bounds>> 
Float maxZt = MaxComponent(Abs(Vector3f(p0t.z, p1t.z, p2t.z))); Float deltaZ = gamma(3) * maxZt;
<<Compute delta Subscript x and delta Subscript y terms for triangle t 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 delta Subscript e term for triangle t error bounds>> 
Float deltaE = 2 * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
<<Compute delta Subscript t term for triangle t 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[3]; GetUVs(uv); <<Compute deltas for triangle partial derivatives>> 
Vector2f duv02 = uv[0] - uv[2], duv12 = uv[1] - uv[2]; Vector3f dp02 = p0 - p2, dp12 = p1 - p2;
Float determinant = duv02[0] * duv12[1] - duv02[1] * duv12[0]; 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[1] * dp02 - duv02[1] * dp12) * invdet; dpdv = (-duv12[0] * dp02 + duv02[0] * 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 left-parenthesis u comma v right-parenthesis parametric coordinates and hit point>> 
Point3f pHit = b0 * p0 + b1 * p1 + b2 * p2; Point2f uvHit = b0 * uv[0] + b1 * uv[1] + b2 * uv[2];
<<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[0]] + b1 * mesh->n[v[1]] + b2 * mesh->n[v[2]]); else ns = isect->n;
<<Compute shading tangent ss for triangle>> 
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v[0]] + b1 * mesh->s[v[1]] + b2 * mesh->s[v[2]]); else ss = Normalize(isect->dpdu);
<<Compute shading bitangent ts for triangle and adjust ss>> 
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);
<<Compute partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v for triangle shading geometry>> 
Normal3f dndu, dndv; if (mesh->n) { <<Compute deltas for triangle partial derivatives of normal>> 
Vector2f duv02 = uv[0] - uv[2]; Vector2f duv12 = uv[1] - uv[2]; Normal3f dn1 = mesh->n[v[0]] - mesh->n[v[2]]; Normal3f dn2 = mesh->n[v[1]] - mesh->n[v[2]];
Float determinant = duv02[0] * duv12[1] - duv02[1] * duv12[0]; if (determinant == 0) dndu = dndv = Normal3f(0, 0, 0); else { Float invDet = 1 / determinant; dndu = ( duv12[1] * dn1 - duv02[1] * dn2) * invDet; dndv = (-duv12[0] * dn1 + duv02[0] * dn2) * invDet; } } else dndu = dndv = Normal3f(0,0,0);
isect->SetShadingGeometry(ss, ts, dndu, dndv, true);
} <<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 left-parenthesis 0 comma 0 comma 0 right-parenthesis in the transformed coordinate system and such that its direction is along the plus z 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 x and y 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 t 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 t 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 t is conservatively greater than zero>> 
<<Compute delta Subscript z term for triangle t error bounds>> 
Float maxZt = MaxComponent(Abs(Vector3f(p0t.z, p1t.z, p2t.z))); Float deltaZ = gamma(3) * maxZt;
<<Compute delta Subscript x and delta Subscript y terms for triangle t 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 delta Subscript e term for triangle t error bounds>> 
Float deltaE = 2 * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
<<Compute delta Subscript t term for triangle t 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 bold upper T , a coordinate permutation bold upper P , and a shear bold upper S . Rather than computing explicit transformation matrices for each of these and then computing an aggregate transformation matrix bold upper M equals bold upper S bold upper P bold upper T 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:

bold upper T equals Start 4 By 4 Matrix 1st Row 1st Column 1 2nd Column 0 3rd Column 0 4th Column minus normal o Subscript x Baseline 2nd Row 1st Column 0 2nd Column 1 3rd Column 0 4th Column minus normal o Subscript y Baseline 3rd Row 1st Column 0 2nd Column 0 3rd Column 1 4th Column minus normal o Subscript z Baseline 4th Row 1st Column 0 2nd Column 0 3rd Column 0 4th Column 1 EndMatrix period

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 z dimension is the one where the absolute value of the ray’s direction is largest. The x and y dimensions are arbitrarily assigned to the other two dimensions. This step ensures that if, for example, the original ray’s z direction is zero, then a dimension with non-zero magnitude is mapped to plus z .

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

bold upper P equals Start 4 By 4 Matrix 1st Row 1st Column 0 2nd Column 1 3rd Column 0 4th Column 0 2nd Row 1st Column 0 2nd Column 0 3rd Column 1 4th Column 0 3rd Row 1st Column 1 2nd Column 0 3rd Column 0 4th Column 0 4th Row 1st Column 0 2nd Column 0 3rd Column 0 4th Column 1 EndMatrix period

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 plus z axis:

bold upper S equals Start 4 By 4 Matrix 1st Row 1st Column 1 2nd Column 0 3rd Column minus bold d Subscript x Baseline slash bold d Subscript z Baseline 4th Column 0 2nd Row 1st Column 0 2nd Column 1 3rd Column minus bold d Subscript y Baseline slash bold d Subscript z Baseline 4th Column 0 3rd Row 1st Column 0 2nd Column 0 3rd Column 1 slash bold d Subscript z Baseline 4th Column 0 4th Row 1st Column 0 2nd Column 0 3rd Column 0 4th Column 1 EndMatrix period

To see how this transformation works, consider its operation on the ray direction vector left-bracket bold d Subscript x Baseline bold d Subscript y Baseline bold d Subscript z Baseline 0 right-bracket Superscript upper T .

For now, only the x and y dimensions are sheared; we can wait and shear the z 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 plus z 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 x , y coordinates left-parenthesis 0 comma 0 right-parenthesis are inside the x y projection of the triangle (Figure 3.12).

Figure 3.12: In the ray–triangle intersection coordinate system, the ray starts at the origin and goes along the plus z axis. The intersection test can be performed by considering only the x y projection of the ray and the triangle vertices, which in turn reduces to determining if the 2D point left-parenthesis 0 comma 0 right-parenthesis is within the triangle.

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 bold a and bold b , the area is

bold a Subscript x Baseline bold b Subscript y Baseline minus bold b Subscript x Baseline bold a Subscript y Baseline period

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 normal p 0 , normal p 1 , and normal p 2 is

one-half left-parenthesis left-parenthesis normal p 1 Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis left-parenthesis normal p 2 Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis minus left-parenthesis normal p 2 Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis right-parenthesis period

Figure 3.13 visualizes this geometrically.

Figure 3.13: The area of a triangle with two edges given by vectors v 1 and v 2 is one-half of the area of the parallelogram shown here. The parallelogram area is given by the length of the cross product of bold v 1 and bold v 2 .

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

e left-parenthesis normal p Subscript Baseline right-parenthesis equals left-parenthesis normal p 1 Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis left-parenthesis normal p Subscript Baseline Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis minus left-parenthesis normal p Subscript Baseline Subscript x Baseline minus normal p 0 Subscript x Baseline right-parenthesis left-parenthesis normal p 1 Subscript y Baseline minus normal p 0 Subscript y Baseline right-parenthesis period
(3.1)

(See Figure 3.14.)

Figure 3.14: The edge function e left-parenthesis normal p Subscript Baseline right-parenthesis characterizes points with respect to an oriented line between two points normal p 0 and normal p 1 . The value of the edge function is positive for points normal p Subscript to the right of the line, zero for points on the line, and negative for points to the left of the line. The ray–triangle intersection algorithm uses an edge function that is twice the signed area of the triangle formed by the three points.

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 normal p Subscript is left-parenthesis 0 comma 0 right-parenthesis . This simplifies the edge function expressions. For example, for the edge e 0 from normal p 1 to normal p 2 , we have:

StartLayout 1st Row 1st Column e 0 left-parenthesis normal p Subscript Baseline right-parenthesis 2nd Column equals left-parenthesis normal p 2 Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis left-parenthesis normal p Subscript Baseline Subscript y Baseline minus normal p 1 Subscript y Baseline right-parenthesis minus left-parenthesis normal p Subscript Baseline Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis left-parenthesis normal p 2 Subscript y Baseline minus normal p 1 Subscript y Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-parenthesis normal p 2 Subscript x Baseline minus normal p 1 Subscript x Baseline right-parenthesis left-parenthesis minus normal p 1 Subscript y Baseline right-parenthesis minus left-parenthesis minus normal p 1 Subscript x Baseline right-parenthesis left-parenthesis normal p 2 Subscript y Baseline minus normal p 1 Subscript y Baseline right-parenthesis 3rd Row 1st Column Blank 2nd Column equals normal p 1 Subscript x Baseline normal p 2 Subscript y Baseline minus normal p 2 Subscript x Baseline normal p 1 Subscript y Baseline period EndLayout
(3.2)

In the following, we’ll use the indexing scheme that the edge function e Subscript i corresponds to the directed edge from vertex normal p Subscript Baseline Subscript left-parenthesis i plus 1 right-parenthesis normal m normal o normal d 3 and normal p Subscript Baseline Subscript left-parenthesis i plus 2 right-parenthesis normal m normal o normal d 3 .

<<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 left-parenthesis 0 comma 0 right-parenthesis 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 plus z axis, the z coordinate value of the intersection point is equal to the intersection’s parametric t value. To compute this z value, we first need to go ahead and apply the shear transformation to the z coordinates of the triangle vertices. Given these z 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:

b Subscript i Baseline equals StartFraction e Subscript i Baseline Over e 0 plus e 1 plus e 2 EndFraction period

Thus, the b Subscript i sum to one.

The interpolated z value is given by

z equals b 0 z 0 plus b 1 z 1 plus b 2 z 2 comma

where z Subscript i 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 b Subscript i in cases where the final t value is out of the range of valid t values, the implementation here first computes t by interpolating z Subscript i with e Subscript i (in other words, not yet performing the division by d equals e 0 plus e 1 plus e 2 ). If the sign of d and the sign of the interpolated t value are different, then the final t value will certainly be negative and thus not a valid intersection.

Along similar lines, the check t less-than t Subscript normal m normal a normal x can be equivalent performed in two ways:

StartLayout 1st Row 1st Column sigma-summation Underscript i Endscripts e Subscript i Baseline z Subscript i Baseline less-than t Subscript normal m normal a normal x Baseline left-parenthesis e 0 plus e 1 plus e 2 right-parenthesis 2nd Column If e 0 plus e 1 plus e 2 greater-than 0 2nd Row 1st Column sigma-summation Underscript i Endscripts e Subscript i Baseline z Subscript i Baseline greater-than t Subscript normal m normal a normal x Baseline left-parenthesis e 0 plus e 1 plus e 2 right-parenthesis 2nd Column otherwise period EndLayout

<<Compute scaled hit distance to triangle and test against ray t 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 t value for the intersection.

<<Compute barycentric coordinates and t 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 partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v using the parametric left-parenthesis u comma v right-parenthesis 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

normal p Subscript normal o Baseline plus u StartFraction partial-differential normal p Over partial-differential u EndFraction plus v StartFraction partial-differential normal p Over partial-differential v EndFraction comma

for some normal p Subscript normal o , where u and v range over the parametric coordinates of the triangle. We also know the three vertex positions normal p Subscript i , i equals 0 comma 1 comma 2 , and the texture coordinates left-parenthesis u Subscript i Baseline comma v Subscript i Baseline right-parenthesis at each vertex. From this it follows that the partial derivatives of normal p Subscript must satisfy

normal p Subscript i Baseline equals normal p Subscript normal o Baseline plus u Subscript i Baseline StartFraction partial-differential normal p Over partial-differential u EndFraction plus v Subscript i Baseline StartFraction partial-differential normal p Over partial-differential v EndFraction period

In other words, there is a unique affine mapping from the 2D left-parenthesis u comma v right-parenthesis 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 partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v , we start by computing the differences normal p 0 minus normal p 2 and normal p 1 minus normal p 2 , giving the matrix equation

Start 2 By 2 Matrix 1st Row 1st Column u 0 minus u 2 2nd Column v 0 minus v 2 2nd Row 1st Column u 1 minus u 2 2nd Column v 1 minus v 2 EndMatrix StartBinomialOrMatrix partial-differential normal p slash partial-differential u Choose partial-differential normal p slash partial-differential v EndBinomialOrMatrix equals StartBinomialOrMatrix normal p 0 minus normal p 2 Choose normal p 1 minus normal p 2 EndBinomialOrMatrix period

Thus,

StartBinomialOrMatrix partial-differential normal p slash partial-differential u Choose partial-differential normal p slash partial-differential v EndBinomialOrMatrix equals Start 2 By 2 Matrix 1st Row 1st Column u 0 minus u 2 2nd Column v 0 minus v 2 2nd Row 1st Column u 1 minus u 2 2nd Column v 1 minus v 2 EndMatrix Superscript negative 1 Baseline StartBinomialOrMatrix normal p 0 minus normal p 2 Choose normal p 1 minus normal p 2 EndBinomialOrMatrix period

Inverting a 2 times 2 matrix is straightforward. The inverse of the left-parenthesis u comma v right-parenthesis differences matrix is

StartFraction 1 Over left-parenthesis u 0 minus u 2 right-parenthesis left-parenthesis v 1 minus v 2 right-parenthesis minus left-parenthesis v 0 minus v 2 right-parenthesis left-parenthesis u 1 minus u 2 right-parenthesis EndFraction Start 2 By 2 Matrix 1st Row 1st Column v 1 minus v 2 2nd Column minus left-parenthesis v 0 minus v 2 right-parenthesis 2nd Row 1st Column minus left-parenthesis u 1 minus u 2 right-parenthesis 2nd Column u 0 minus u 2 EndMatrix period

<<Compute triangle partial derivatives>>= 
Vector3f dpdu, dpdv; Point2f uv[3]; GetUVs(uv); <<Compute deltas for triangle partial derivatives>> 
Vector2f duv02 = uv[0] - uv[2], duv12 = uv[1] - uv[2]; Vector3f dp02 = p0 - p2, dp12 = p1 - p2;
Float determinant = duv02[0] * duv12[1] - duv02[1] * duv12[0]; 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[1] * dp02 - duv02[1] * dp12) * invdet; dpdv = (-duv12[0] * dp02 + duv02[0] * dp12) * invdet; }

<<Compute deltas for triangle partial derivatives>>= 
Vector2f duv02 = uv[0] - uv[2], duv12 = uv[1] - uv[2]; 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 left-parenthesis u comma v right-parenthesis parametric coordinates at the hit point, the barycentric interpolation formula is applied to the vertex positions and the left-parenthesis u comma v right-parenthesis 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 left-parenthesis u comma v right-parenthesis parametric coordinates and hit point>>= 
Point3f pHit = b0 * p0 + b1 * p1 + b2 * p2; Point2f uvHit = b0 * uv[0] + b1 * uv[1] + b2 * uv[2];

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

<<Triangle Private Methods>>= 
void GetUVs(Point2f uv[3]) const { if (mesh->uv) { uv[0] = mesh->uv[v[0]]; uv[1] = mesh->uv[v[1]]; uv[2] = mesh->uv[v[2]]; } else { uv[0] = Point2f(0, 0); uv[1] = Point2f(1, 0); uv[2] = 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 left-parenthesis 0 comma 0 comma 0 right-parenthesis , 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[0]] + b1 * mesh->n[v[1]] + b2 * mesh->n[v[2]]); else ns = isect->n;
<<Compute shading tangent ss for triangle>> 
Vector3f ss; if (mesh->s) ss = Normalize(b0 * mesh->s[v[0]] + b1 * mesh->s[v[1]] + b2 * mesh->s[v[2]]); else ss = Normalize(isect->dpdu);
<<Compute shading bitangent ts for triangle and adjust ss>> 
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);
<<Compute partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v for triangle shading geometry>> 
Normal3f dndu, dndv; if (mesh->n) { <<Compute deltas for triangle partial derivatives of normal>> 
Vector2f duv02 = uv[0] - uv[2]; Vector2f duv12 = uv[1] - uv[2]; Normal3f dn1 = mesh->n[v[0]] - mesh->n[v[2]]; Normal3f dn2 = mesh->n[v[1]] - mesh->n[v[2]];
Float determinant = duv02[0] * duv12[1] - duv02[1] * duv12[0]; if (determinant == 0) dndu = dndv = Normal3f(0, 0, 0); else { Float invDet = 1 / determinant; dndu = ( duv12[1] * dn1 - duv02[1] * dn2) * invDet; dndv = (-duv12[0] * dn1 + duv02[0] * dn2) * invDet; } } else dndu = dndv = Normal3f(0,0,0);
isect->SetShadingGeometry(ss, ts, dndu, dndv, true);
} <<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;

3.6.3 Shading Geometry

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.

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

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[0]] + b1 * mesh->n[v[1]] + b2 * mesh->n[v[2]]); 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[0]] + b1 * mesh->s[v[1]] + b2 * mesh->s[v[2]]); 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 bold n Subscript and bold s values are provided and if the interpolated bold n Subscript and bold s values aren’t perfectly orthogonal, bold n Subscript will be preserved and bold s will be modified so that the coordinate system is orthogonal.

<<Compute shading bitangent ts for triangle and adjust ss>>= 
Vector3f ts = Cross(ns, ss); if (ts.LengthSquared() > 0.f) { ts = Normalize(ts); ss = Cross(ts, ns); } else CoordinateSystem((Vector3f)ns, &ss, &ts);

The code to compute the partial derivatives partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v of the shading normal is almost identical to the code to compute the partial derivatives partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v . Therefore, it has been elided from the text here.

3.6.4 Surface Area

Using the fact that the area of a parallelogram is given by the length of the cross product of the two vectors along its sides, the Area() method computes the triangle area as half the area of the parallelogram formed by two of its edge vectors (Figure 3.13).

<<Triangle Method Definitions>>+=  
Float Triangle::Area() const { <<Get triangle vertices in p0, p1, and p2>> 
const Point3f &p0 = mesh->p[v[0]]; const Point3f &p1 = mesh->p[v[1]]; const Point3f &p2 = mesh->p[v[2]];
return 0.5 * Cross(p1 - p0, p2 - p0).Length(); }