Make vtk contour take an affine

dipy/viz/actor.py

@@ -201,6 +201,127 @@ def copy(self):
     return image_actor
 
 
+def surface_actor(data, affine=None,
+                  color=np.array([1, 0, 0]), opacity=1):
+    """Take a binary roi and generate a surface actor of the specified color
+    and opacity
+
+    Parameters
+    ----------
+    data : array, shape (X, Y, Z)
+        an ROI file that will be binarized and displayed
+    affine : array, shape (4, 4)
+        Grid to space (usually RAS 1mm) transformation matrix. Default is None.
+        If None then the identity matrix is used.
+    color : (1, 3) ndarray
+        RGB values in [0,1].
+    opacity : float
+        Opacity of surface between 0 and 1.
+
+    Returns
+    -------
+    contour_assembly : vtkAssembly
+        ROI surface object displayed in space
+        coordinates as calculated by the affine parameter.
+
+    """
+
+    if data.ndim != 3:
+        raise ValueError('Only 3D arrays are currently supported.')
+    else:
+        nb_components = 1
+
+    data = (data>0)*1
+    vol = np.interp(data, xp=[data.min(), data.max()], fp=[0, 255])
+    vol = vol.astype('uint8')
+
+    im = vtk.vtkImageData()
+    if major_version <= 5:
+        im.SetScalarTypeToUnsignedChar()
+    di, dj, dk = vol.shape[:3]
+    im.SetDimensions(di, dj, dk)
+    voxsz = (1., 1., 1.)
+    # im.SetOrigin(0,0,0)
+    im.SetSpacing(voxsz[2], voxsz[0], voxsz[1])
+    if major_version <= 5:
+        im.AllocateScalars()
+        im.SetNumberOfScalarComponents(nb_components)
+    else:
+        im.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, nb_components)
+
+    # copy data
+    vol = np.swapaxes(vol, 0, 2)
+    vol = np.ascontiguousarray(vol)
+
+    if nb_components == 1:
+        vol = vol.ravel()
+    else:
+        vol = np.reshape(vol, [np.prod(vol.shape[:3]), vol.shape[3]])
+
+    uchar_array = numpy_support.numpy_to_vtk(vol, deep=0)
+    im.GetPointData().SetScalars(uchar_array)
+
+    if affine is None:
+        affine = np.eye(4)
+
+    # Set the transform (identity if none given)
+    transform = vtk.vtkTransform()
+    transform_matrix = vtk.vtkMatrix4x4()
+    transform_matrix.DeepCopy((
+        affine[0][0], affine[0][1], affine[0][2], affine[0][3],
+        affine[1][0], affine[1][1], affine[1][2], affine[1][3],
+        affine[2][0], affine[2][1], affine[2][2], affine[2][3],
+        affine[3][0], affine[3][1], affine[3][2], affine[3][3]))
+    transform.SetMatrix(transform_matrix)
+    transform.Inverse()
+
+    # Set the reslicing
+    image_resliced = vtk.vtkImageReslice()
+    set_input(image_resliced, im)
+    image_resliced.SetResliceTransform(transform)
+    image_resliced.AutoCropOutputOn()
+
+    # Adding this will allow to support anisotropic voxels
+    # and also gives the opportunity to slice per voxel coordinates
+
+    rzs = affine[:3, :3]
+    zooms = np.sqrt(np.sum(rzs * rzs, axis=0))
+    image_resliced.SetOutputSpacing(*zooms)
+
+    image_resliced.SetInterpolationModeToLinear()
+    image_resliced.Update()
+
+    # code from fvtk contour function
+
+    skin_extractor = vtk.vtkContourFilter()
+    if major_version <= 5:
+        skin_extractor.SetInput(image_resliced)
+    else:
+        skin_extractor.SetInputData(image_resliced.GetOutput())
+
+    skin_extractor.SetValue(0,1)
+
+    skin_normals = vtk.vtkPolyDataNormals()
+    skin_normals.SetInputConnection(skin_extractor.GetOutputPort())
+    skin_normals.SetFeatureAngle(60.0)
+
+    skin_mapper = vtk.vtkPolyDataMapper()
+    skin_mapper.SetInputConnection(skin_normals.GetOutputPort())
+    skin_mapper.ScalarVisibilityOff()
+
+    skin_actor = vtk.vtkActor()
+
+    skin_actor.SetMapper(skin_mapper)
+    skin_actor.GetProperty().SetOpacity(opacity)
+
+    skin_actor.GetProperty().SetColor(color[0], color[1], color[2])
+
+    del skin_mapper
+    del skin_extractor
+
+    return skin_actor
+
+
 def streamtube(lines, colors=None, opacity=1, linewidth=0.1, tube_sides=9,
                lod=True, lod_points=10 ** 4, lod_points_size=3,
                spline_subdiv=None, lookup_colormap=None):

dipy/viz/tests/test_fvtk_actors.py

@@ -23,10 +23,13 @@
 if actor.have_vtk:
     if actor.major_version == 5 and use_xvfb:
         skip_slicer = True
+        skip_surface = True
     else:
         skip_slicer = False
+        skip_surface = False
 else:
     skip_slicer = False
+    skip_surface = False
 
 
 @npt.dec.skipif(skip_slicer)
@@ -156,6 +159,113 @@ def test_slicer():
     npt.assert_array_equal([1, 3, 2] * np.array(data.shape),
                            np.array(slicer.shape))
 
+@npt.dec.skipif(skip_surface)
+@npt.dec.skipif(not run_test)
+@xvfb_it
+def test_surface():
+
+    #Render volume
+    renderer = window.renderer()
+    data = np.zeros((50, 50, 50))
+    data[20:30,25,25]=1.
+    data[25, 20:30, 25] = 1.
+    affine = np.eye(4)
+    surface = actor.surface_actor(data, affine,
+                                  colors=np.array([1, 0, 1]),
+                                  opacity=.5)
+    renderer.add(surface)
+
+    renderer.reset_camera()
+    renderer.reset_clipping_range()
+    #window.show(renderer)
+
+    #Test binarization
+    renderer2 = window.renderer()
+    data2 = np.zeros((50, 50, 50))
+    data2[20:30, 25, 25] = 1.
+    data2[35:40, 25, 25] = 1.
+    affine = np.eye(4)
+    surface2 = actor.surface_actor(data2, affine,
+                                  colors=np.array([0, 1, 1]),
+                                  opacity=.5)
+    renderer2.add(surface2)
+
+    renderer2.reset_camera()
+    renderer2.reset_clipping_range()
+    #window.show(renderer2)
+
+    arr = window.snapshot(renderer, 'test_surface.png', offscreen=False)
+    arr2 = window.snapshot(renderer2, 'test_surface2.png', offscreen=False)
+
+    report = window.analyze_snapshot(arr, find_objects=True)
+    report2 = window.analyze_snapshot(arr2, find_objects=True)
+
+    npt.assert_equal(report.objects, 1)
+    npt.assert_equal(report2.objects, 2)
+
+    print(report)
+
+
+    #test on real streamlines using tracking example
+    from dipy.data import read_stanford_labels
+    from dipy.reconst.shm import CsaOdfModel
+    from dipy.data import default_sphere
+    from dipy.direction import peaks_from_model
+    from dipy.tracking.local import ThresholdTissueClassifier
+    from dipy.tracking import utils
+    from dipy.tracking.local import LocalTracking
+    from dipy.viz import fvtk
+    from dipy.viz.colormap import line_colors
+
+
+    hardi_img, gtab, labels_img = read_stanford_labels()
+    data = hardi_img.get_data()
+    labels = labels_img.get_data()
+    affine = hardi_img.get_affine()
+
+    white_matter = (labels == 1) | (labels == 2)
+
+    csa_model = CsaOdfModel(gtab, sh_order=6)
+    csa_peaks = peaks_from_model(csa_model, data, default_sphere,
+                                 relative_peak_threshold=.8,
+                                 min_separation_angle=45,
+                                 mask=white_matter)
+
+    classifier = ThresholdTissueClassifier(csa_peaks.gfa, .25)
+
+    seed_mask = labels == 2
+    seeds = utils.seeds_from_mask(seed_mask, density=[1, 1, 1], affine=affine)
+
+    # Initialization of LocalTracking. The computation happens in the next step.
+    streamlines = LocalTracking(csa_peaks, classifier, seeds, affine,
+                                step_size=2)
+
+    # Compute streamlines and store as a list.
+    streamlines = list(streamlines)
+
+    # Prepare the display objects.
+    streamlines_actor = fvtk.line(streamlines, line_colors(streamlines))
+    seedroi_actor = actor.surface_actor(seed_mask, affine, [0, 1, 1], 0.5)
+
+    # Create the 3d display.
+    r = fvtk.ren()
+    r2 = fvtk.ren()
+    fvtk.add(r, streamlines_actor)
+    arr3 = window.snapshot(r, 'test_surface3.png', offscreen=False)
+    report3 = window.analyze_snapshot(arr3, find_objects=True)
+    fvtk.add(r2, streamlines_actor)
+    fvtk.add(r2, seedroi_actor)
+    arr4 = window.snapshot(r2, 'test_surface4.png', offscreen=False)
+    report4 = window.analyze_snapshot(arr4, find_objects=True)
+
+    #assert that the seed ROI rendering isn't far away from the streamlines (affine error)
+    npt.assert_equal(report3.objects,report4.objects)
+    #fvtk.show(r)
+    #fvtk.show(r2)
+
+
+
+
 
 @npt.dec.skipif(not run_test)
 @xvfb_it

doc/examples/viz_surface_actor.py

@@ -0,0 +1,98 @@
+"""
+==================================
+Visualization of 3D Surface Rendered ROI with Streamlines
+==================================
+
+Here is a simple tutorial following the probabilistic CSA Tracking Example in
+which we generate a dataset of streamlines from a corpus callosum ROI, and
+then display them with the seed ROI rendered in 3D with 50% transparency.
+
+"""
+
+from dipy.data import read_stanford_labels
+from dipy.reconst.shm import CsaOdfModel
+from dipy.data import default_sphere
+from dipy.direction import peaks_from_model
+from dipy.tracking.local import ThresholdTissueClassifier
+from dipy.tracking import utils
+from dipy.tracking.local import LocalTracking
+from dipy.viz import fvtk
+from dipy.viz.colormap import line_colors
+
+"""
+First, we need to generate some streamlines. For a more complete description of
+these steps, please refer to the CSA Probabilistic Tracking Tutorial.
+"""
+
+hardi_img, gtab, labels_img = read_stanford_labels()
+data = hardi_img.get_data()
+labels = labels_img.get_data()
+affine = hardi_img.get_affine()
+
+white_matter = (labels == 1) | (labels == 2)
+
+csa_model = CsaOdfModel(gtab, sh_order=6)
+csa_peaks = peaks_from_model(csa_model, data, default_sphere,
+                                 relative_peak_threshold=.8,
+                                 min_separation_angle=45,
+                                 mask=white_matter)
+
+classifier = ThresholdTissueClassifier(csa_peaks.gfa, .25)
+
+seed_mask = labels == 2
+seeds = utils.seeds_from_mask(seed_mask, density=[1, 1, 1], affine=affine)
+
+# Initialization of LocalTracking. The computation happens in the next step.
+streamlines = LocalTracking(csa_peaks, classifier, seeds, affine,
+                                step_size=2)
+
+# Compute streamlines and store as a list.
+streamlines = list(streamlines)
+
+"""
+We will create a streamline actor from the streamlines
+"""
+
+streamlines_actor = fvtk.line(streamlines, line_colors(streamlines))
+
+"""
+Next, we create a surface actor from the corpus callosum seed ROI. We
+provide the ROI data, the affine, the color in [R,G,B], and the opacity as
+a decimal between zero and one. Here, we set the color as blue/green with
+50% opacity.
+"""
+surface_opacity = 0.5
+surface_color = [0,1,1]
+
+seedroi_actor = fvtk.actor.surface_actor(seed_mask, affine,
+                                         surface_color, surface_opacity)
+
+"""
+Next, we initialize a ''Renderer'' object and add both of the actors
+to the rendering.
+"""
+
+ren = fvtk.ren()
+fvtk.add(ren, streamlines_actor)
+fvtk.add(ren, seedroi_actor)
+
+"""
+If you uncomment the following line, the rendering will pop up in an interactive
+window.
+"""
+
+#fvtk.show(ren)
+
+ren.zoom(1.5)
+ren.reset_clipping_range()
+
+window.record(ren, out_path='surface_actor_tutorial.png', size=(1200, 900),
+              reset_camera=False)
+
+"""
+.. figure:: surface_actor_tutorial.png
+   :align: center
+
+   **A top view of corpus callosum streamlines with the blue transparent
+   seed ROI in the center**.
+"""