FIX: Add test for normals

Author: larsoner

surfer/tests/test_utils.py

@@ -10,12 +10,36 @@
 data_path = pjoin(subj_dir, subject_id)
 
 
+def _slow_compute_normals(rr, tris):
+    """Efficiently compute vertex normals for triangulated surface"""
+    # first, compute triangle normals
+    r1 = rr[tris[:, 0], :]
+    r2 = rr[tris[:, 1], :]
+    r3 = rr[tris[:, 2], :]
+    tri_nn = np.cross((r2 - r1), (r3 - r1))
+
+    #   Triangle normals and areas
+    size = np.sqrt(np.sum(tri_nn * tri_nn, axis=1))
+    zidx = np.where(size == 0)[0]
+    size[zidx] = 1.0  # prevent ugly divide-by-zero
+    tri_nn /= size[:, np.newaxis]
+
+    # accumulate the normals
+    nn = np.zeros((len(rr), 3))
+    for p, verts in enumerate(tris):
+        nn[verts] += tri_nn[p, :]
+    size = np.sqrt(np.sum(nn * nn, axis=1))
+    size[size == 0] = 1.0  # prevent ugly divide-by-zero
+    nn /= size[:, np.newaxis]
+    return nn
+
+
 @utils.requires_fsaverage
 def test_surface():
     """Test IO for Surface class"""
     for subjects_dir in [None, subj_dir]:
         surface = utils.Surface('fsaverage', 'lh', 'inflated',
-                             subjects_dir=subjects_dir)
+                                subjects_dir=subjects_dir)
         surface.load_geometry()
         surface.load_label('BA1')
         surface.load_curvature()
@@ -25,3 +49,7 @@ def test_surface():
         surface.apply_xfm(xfm)
         x_ = surface.x
         assert_array_almost_equal(x + 2, x_)
+
+        # normals
+        nn = _slow_compute_normals(surface.coords, surface.faces)
+        assert_array_almost_equal(nn, surface.nn)

surfer/utils.py

@@ -171,7 +171,6 @@ def _fast_cross_3d(x, y):
 def _compute_normals(rr, tris):
     """Efficiently compute vertex normals for triangulated surface"""
     # first, compute triangle normals
-    t0 = time.time()
     r1 = rr[tris[:, 0], :]
     r2 = rr[tris[:, 1], :]
     r3 = rr[tris[:, 2], :]
@@ -197,6 +196,9 @@ def _compute_normals(rr, tris):
         vals = np.r_[np.zeros((1, 3)), np.cumsum(tri_nn[reord, :], 0)]
         idx = np.cumsum(np.r_[0, counts])
         nn += vals[idx[1:], :] - vals[idx[:-1], :]
+    size = np.sqrt(np.sum(nn * nn, axis=1))
+    size[size == 0] = 1.0  # prevent ugly divide-by-zero
+    nn /= size[:, np.newaxis]
     return nn