diff --git a/src/polymesh/algorithms/normal_estimation.hh b/src/polymesh/algorithms/normal_estimation.hh
index 21464054b302f5ccdfc8244d8b3e73a216b29b32..a1f3bd19b68f2abb37d4fb1f2f0f38663002ff04 100644
--- a/src/polymesh/algorithms/normal_estimation.hh
+++ b/src/polymesh/algorithms/normal_estimation.hh
@@ -19,35 +19,41 @@ template <class Vec3, class IsHardEdgeF>
 
     auto const hard_edges = m.edges().map([&](edge_handle e) { return e.is_boundary() ? true : is_hard_edge(e); });
 
-    return m.halfedges().map([&](halfedge_handle h) {
-        Vec3 n = face_normals[h.face()];
+    return m.halfedges().map(
+        [&](halfedge_handle h)
+        {
+            if (h.is_boundary())
+                return Vec3{};
 
-        auto h0 = h;
-        auto h1 = h.next().opposite();
+            POLYMESH_ASSERT(h.face().is_valid());
+            Vec3 n = face_normals[h.face()];
 
-        // iterate over h0
-        auto hh = h0;
-        while (hh != h1 && !hard_edges[hh])
-        {
-            hh = hh.opposite().prev();
-            n += face_normals[hh.face()];
-        }
+            auto h0 = h;
+            auto h1 = h.next().opposite();
 
-        // iterate over h1 if not reached around
-        if (hh != h1)
-        {
-            hh = h1;
-            while (hh != h0 && !hard_edges[hh])
+            // iterate over h0
+            auto hh = h0;
+            while (hh != h1 && !hard_edges[hh])
             {
+                hh = hh.opposite().prev();
                 n += face_normals[hh.face()];
-                hh = hh.next().opposite();
             }
-        }
 
-        // normalize
-        auto l = std::sqrt(n.x * n.x + n.y * n.y + n.z * n.z);
-        return n / (l + 1e-30f);
-    });
+            // iterate over h1 if not reached around
+            if (hh != h1)
+            {
+                hh = h1;
+                while (hh != h0 && !hard_edges[hh])
+                {
+                    n += face_normals[hh.face()];
+                    hh = hh.next().opposite();
+                }
+            }
+
+            // normalize
+            auto l = std::sqrt(n.x * n.x + n.y * n.y + n.z * n.z);
+            return n / (l + 1e-30f);
+        });
 }
 
 /// same as normal_estimation(face_attribute<Vec3>, ...)