Expand skinseg test to show vtk 9.2.2 issue (OpenwaterHealth#237)

Author: ebrahimebrahim

tests/test_skinseg.py

@@ -37,6 +37,20 @@ def add_ball(
 
     volume[dist_sq <= radius**2] = value
 
+def apply_affine_to_polydata(affine:np.ndarray, polydata:vtk.vtkPolyData) -> vtk.vtkPolyData:
+    """Apply an affine transform to a vtkPolyData"""
+    affine_vtkmat = vtk.vtkMatrix4x4()
+    for i in range(4):
+        for j in range(4):
+            affine_vtkmat.SetElement(i, j, affine[i, j])
+    affine_vtktransform = vtk.vtkTransform()
+    affine_vtktransform.SetMatrix(affine_vtkmat)
+    transform_filter = vtk.vtkTransformPolyDataFilter()
+    transform_filter.SetTransform(affine_vtktransform)
+    transform_filter.SetInputData(polydata)
+    transform_filter.Update()
+    return transform_filter.GetOutput()
+
 def test_take_largest_connected_component():
     vol_array = np.zeros((20,20,20))
     add_ball(vol_array, (5,5,5), 4) # ball of radius 4 at (5,5,5)
@@ -76,11 +90,18 @@ def test_create_closed_surface_from_labelmap():
     sphere_radius = 7
     sphere_center = np.array([10,10,10])
     add_ball(labelmap, tuple(sphere_center), sphere_radius)
-    labelmap_vtk = vtk_img_from_array_and_affine(labelmap, affine = np.eye(4))
+    rng = np.random.default_rng(6548)
+    affine = np.eye(4)
+    affine[:3,:3] = expm((lambda A: (A - A.T)/2)(rng.normal(size=(3,3)))) # generate a random orthogonal matrix
+    affine[:3,3] = rng.random(3) # generate a random origin
+    labelmap_vtk = vtk_img_from_array_and_affine(labelmap, affine = affine)
 
     # run the algorithm to be tested
     surface = create_closed_surface_from_labelmap(labelmap_vtk, decimation_factor=0.1)
 
+    # the mesh is in "physical space" mapped to by `affine`, transform it back to the ijk space of the original `labelmap`
+    surface = apply_affine_to_polydata(np.linalg.inv(affine), surface)
+
     # verify that the points on the generated mesh are not too far off being at distance 7 from the ball center
     points = surface.GetPoints()
     for i in range(points.GetNumberOfPoints()):