diff --git a/lib/rendering/pbr/light/MeshLight.cc b/lib/rendering/pbr/light/MeshLight.cc
index 0629866..57bb021 100644
--- a/lib/rendering/pbr/light/MeshLight.cc
+++ b/lib/rendering/pbr/light/MeshLight.cc
@@ -805,7 +805,11 @@ MeshLight::setMesh(geom::internal::Primitive* prim)
             }
         }
         mArea += area;
-        faces[f].mInvArea = 1 / area;
+        // Guard against degenerate (zero-area) faces storing inf in mInvArea.
+        // Ideally, zero-area (which are also zero energy) faces will not be chosen.
+        // However, if both branches have zero energy, we choose one branch randomly
+        // (see the edge case in drawSampleRecurse). This guard prevents inf from propagating.
+        faces[f].mInvArea = (area > 0.0f) ? (1.0f / area) : 0.0f;
 
         // get energy
         // TODO: more fine detail on texture sampling? Take average instead of max?
@@ -1633,8 +1637,9 @@ MeshLight::intersect(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::V
 }
 
 bool
-MeshLight::sample(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::Vec3f *n, float time, const scene_rdl2::math::Vec3f& r,
-    scene_rdl2::math::Vec3f &wi, LightIntersection &isect, float rayDirFootprint) const
+MeshLight::sample(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::Vec3f *n, float time, 
+                  const scene_rdl2::math::Vec3f& r, scene_rdl2::math::Vec3f &wi, LightIntersection &isect, 
+                  float rayDirFootprint) const
 {
     MNRY_ASSERT(mOn);
 
@@ -1654,7 +1659,10 @@ MeshLight::sample(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::Vec3
         faceIndex = drawSampleRecurse(transformedP, nullptr, 0, r3, pdf);
     }
 
-    MNRY_ASSERT(scene_rdl2::math::isfinite(pdf) && pdf >= 0.0f);
+    MNRY_ASSERT(scene_rdl2::math::isfinite(pdf) && pdf > 0.0f);
+    if (!scene_rdl2::math::isfinite(pdf) || pdf <= 0.0f) {
+        return false;
+    }
     MNRY_ASSERT(mBVH[faceIndex].isLeaf());
     const Face& face = *mBVH[faceIndex].mFace;
 
@@ -1668,6 +1676,20 @@ MeshLight::sample(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::Vec3
         scene_rdl2::math::Vec3f p2 = getFaceVertex(face, 1, time);
         scene_rdl2::math::Vec3f p3 = getFaceVertex(face, 2, time);
 
+        if (mDeformationMb) {
+            const scene_rdl2::math::Vec3f crossProd = cross(p2 - p1, p3 - p1);
+            const float crossLen = length(crossProd);
+            // Handle degenerate, zero-area triangles
+            if (crossLen < scene_rdl2::math::sEpsilon) {
+                return false;
+            }
+            // We cannot directly use face.mNormal here because it is the normal at
+            // time = centroidTime. Therefore we compute it here.
+            normal = crossProd / crossLen;
+        } else {
+            normal = face.mNormal;
+        }
+
         // the random numbers r1 and r2 can be used as uv coordinates
         float u = r1;
         float v = r2;
@@ -1680,14 +1702,6 @@ MeshLight::sample(const scene_rdl2::math::Vec3f &p, const scene_rdl2::math::Vec3
         }
         hit = w*p1 + u*p2 + v*p3;
 
-        if (mDeformationMb) {
-            // We cannot directly use face.mNormal here because it is the normal at
-            // time = centroidTime. Therefore we compute it here.
-            normal = normalize(cross(p2 - p1, p3 - p1));
-        } else {
-            normal = face.mNormal;
-        }
-
         if (mMapShader) {
             // get barycentric coordinate for hit
             uv.x = u;
diff --git a/lib/rendering/pbr/light/MeshLight.ispc b/lib/rendering/pbr/light/MeshLight.ispc
index bfda566..77af541 100644
--- a/lib/rendering/pbr/light/MeshLight.ispc
+++ b/lib/rendering/pbr/light/MeshLight.ispc
@@ -397,7 +397,10 @@ MeshLight_sample(const uniform Light * uniform li,
         faceIndex = MeshLight_drawSampleRecurse(light, transformedP, nullptr, 0, r3, pdf);
     }
 
-    MNRY_ASSERT(isfinite(pdf) && pdf >= 0.0f);
+    MNRY_ASSERT(isfinite(pdf) && pdf > 0.0f);
+    if (!isfinite(pdf) || pdf <= 0.0f) {
+        return false;
+    }
     MNRY_ASSERT(Node_isLeaf(&light->mBVHPtr[faceIndex]));
 
     const uniform Face * varying face = light->mBVHPtr[faceIndex].mFace;
@@ -412,6 +415,20 @@ MeshLight_sample(const uniform Light * uniform li,
         Vec3f p2 = MeshLight_getFaceVertex(light, face, 1, time);
         Vec3f p3 = MeshLight_getFaceVertex(light, face, 2, time);
 
+        if (light->mDeformationMb) {
+            const Vec3f crossProd = cross(p2 - p1, p3 - p1);
+            const float crossLen = length(crossProd);
+            // Handle degenerate, zero-area triangles
+            if (crossLen < sEpsilon) {
+                return false;
+            }
+            // We cannot directly use face.mNormal here because it is the normal at
+            // time = centroidTime. Therefore we compute it here.
+            normal = crossProd / crossLen;
+        } else {
+            normal = face->mNormal;
+        }
+
         // the random numbers r1 and r2 can be used as uv coordinates
         float u = r1;
         float v = r2;
@@ -424,14 +441,6 @@ MeshLight_sample(const uniform Light * uniform li,
         }
         hit = w*p1 + u*p2 + v*p3;
 
-        if (light->mDeformationMb) {
-            // We cannot directly use face.mNormal here because it is the normal at
-            // time = centroidTime. Therefore we compute it here.
-            normal = normalize(cross(p2 - p1, p3 - p1));
-        } else {
-            normal = face->mNormal;
-        }
-
         if (light->mMapShader) {
             // get barycentric coordinate for hit
             uv.x = u;