dicom

README.md

@@ -28,7 +28,7 @@ in a program, whilst mantaining access to the full range of VTK native classes.
 It includes a [large set of working examples](https://github.com/marcomusy/vtkplotter/tree/master/examples)
 for the all following functionalities:
 
-  - Import meshes from VTK format, STL, Wavefront OBJ, 3DS, XML, Neutral, GMSH, OFF, PCD (PointCloud), volumetric TIFF stacks, SLC, MHD, 2D images PNG, JPEG.
+  - Import meshes from VTK format, STL, Wavefront OBJ, 3DS, XML, Neutral, GMSH, OFF, PCD (PointCloud), volumetric TIFF stacks, DICOM, SLC, MHD, 2D images PNG, JPEG.
   - Export meshes as ASCII or binary to VTK, STL, OBJ, PLY formats.
   - Mesh analysis through the built-in methods of VTK package. Additional analysis tools like *Moving Least Squares*, mesh morphing.
   - Tools to visualize and edit meshes (cutting a mesh with another mesh, slicing, normalizing, moving vertex positions, etc..). Interactive cutter widget.
@@ -79,7 +79,7 @@ Voxel-data (_vti, slc, tiff_) files can also be visualized with options `-g` and
 e.g.:
 ```bash
 vtkplotter -g -c blue examples/data/embryo.slc  # (3D scan of a mouse embryo)
-vtkplotter --slicer   examples/data/embryo.slc    
+vtkplotter --slicer   examples/data/embryo.slc  # can read a DICOM directory  
 ```
 ![e2](https://user-images.githubusercontent.com/32848391/50738810-58af4380-11d8-11e9-8fc7-6c6959207224.jpg)
 

bin/vtkplotter

@@ -253,29 +253,39 @@ if args.ray_cast_mode or args.z_spacing or args.y_spacing:
     exit()
 
 ##########################################################
-# special case of SLC/TIFF volumes with --slicer option
+# special case of SLC/TIFF/DICOM volumes with --slicer option
 elif args.slicer:
 
     filename = args.files[0]
     if not os.path.isfile(filename):
-        printc("~times File not found:", filename, c=1)
-        exit()
+        if os.path.isdir(filename):
+            reader = vtk.vtkDICOMImageReader()
+            reader.SetDirectoryName(filename)
+        else:
+            printc("~times File not found:", filename, c=1)
+            exit()
+
     printc("...loading", filename, c="m", bold=0)
 
     if ".tif" in filename.lower():
         reader = vtk.vtkTIFFReader()
+        reader.SetFileName(filename)
 
     elif ".slc" in filename.lower():
         reader = vtk.vtkSLCReader()
+        reader.SetFileName(filename)
         if not reader.CanReadFile(filename):
             printc("~bomb Bad SLC file:", filename, c=1)
 
     elif ".vti" in filename.lower():
         reader = vtk.vtkXMLImageDataReader()
+        reader.SetFileName(filename)
+
+    elif ".dcm" in filename.lower() or ".dicom" in filename.lower():
+        reader = vtk.vtkDICOMImageReader()
+        dirname = os.path.dirname(filename)
+        reader.SetDirectoryName(dirname)
 
-    else:
-        printc("~sad Use --slicer option only with voxel data files [slc,tif,vti].", c=1)
-    reader.SetFileName(filename)
     reader.Update()
     img = reader.GetOutput()
 

docs/source/conf.py

@@ -36,7 +36,7 @@
 author = 'Marco Musy'
 
 # The short X.Y version
-version = '2019.1.2'
+version = '2019.1.3'
 
 
 # -- General configuration ---------------------------------------------------

examples/basic/distance2mesh.py

@@ -0,0 +1,16 @@
+"""
+Computes the (signed) distance 
+from one mesh to another.
+"""
+from vtkplotter import Sphere, Cube, show, Text
+
+s1 = Sphere()
+s2 = Cube(pos=[1,0,0], c='white', alpha=0.4)
+
+s1.distanceToMesh(s2, signed=True, negate=False)
+
+s1.addScalarBar(title='Signed\nDistance')
+
+#print(s1.scalars("Distance"))
+
+show(s1, s2, Text(__doc__))

examples/other/dolfin/markmesh.py

@@ -11,9 +11,9 @@ def inside(self, x, on_boundary):
         return on_boundary and abs(x[0]) < DOLFIN_EPS
 left = left()
 
-tcond = MeshFunction("double", mesh, 0)
-tcond.set_all(0.0)
-left.mark(tcond, 3.14)
+tcond = MeshFunction("size_t", mesh, 0)
+tcond.set_all(0)
+left.mark(tcond, 1)
 
 ##################################
 from vtkplotter.dolfin import plot

examples/other/spherical_harmonics1.py

@@ -75,4 +75,4 @@
         pts.append(p)
 
 act = Points(pts, c="r", r=8, alpha=0.5)
-vp.show(at=1, actors=act, interactive=1)
+vp.show(act, at=1, interactive=1)

examples/other/spherical_harmonics2.py

@@ -6,7 +6,6 @@
 
 try:
     import pyshtools
-
     print(__doc__)
 except:
     print("Please install pyshtools to run this example")
@@ -69,10 +68,10 @@ def morph(clm1, clm2, t, lmax):
 shape2 = vp.load(datadir+"shapes/icosahedron.vtk").normalize().lineWidth(1)
 
 agrid1, actorpts1 = makeGrid(shape1, N)
-vp.show(at=0, actors=[shape1, actorpts1])
+vp.show(shape1, actorpts1, at=0)
 
 agrid2, actorpts2 = makeGrid(shape2, N)
-vp.show(at=1, actors=[shape2, actorpts2])
+vp.show(shape2, actorpts2, at=1)
 vp.camera.Zoom(1.2)
 vp.interactive = False
 
@@ -84,8 +83,8 @@ def morph(clm1, clm2, t, lmax):
 for t in arange(0, 1, 0.005):
     act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
     act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)
-    vp.show(at=2, actors=act21, resetcam=0, legend="time: " + str(int(t * 100)))
-    vp.show(at=3, actors=act12)
+    vp.show(act21, at=2, resetcam=0)
+    vp.show(act12, at=3)
     vp.camera.Azimuth(2)
 
 vp.show(interactive=1)

examples/other/voronoi3d.py

@@ -0,0 +1,19 @@
+'''
+Voronoi in 3D with Voro++ library.
+'''
+from vtkplotter import voronoi3D, Points, show
+import numpy as np
+
+N = 2000
+nuclei = np.random.rand(N, 3) - (0.5,0.5,0.5)
+ncl = Points(nuclei).clean(0.1) # clean makes points evenly spaced
+nuclei = ncl.coordinates()
+
+actor = voronoi3D(nuclei, tol=.001)
+#print(len(actor.info['cells']), actor.info['volumes'])
+
+pts_inside = actor.insidePoints(nuclei)
+inpts = Points(pts_inside, r=50, c='r', alpha=0.2)
+
+show(actor, inpts)
+

examples/run_all.sh

@@ -42,6 +42,9 @@ python basic/buildpolydata.py
 echo Running basic/delaunay2d.py
 python basic/delaunay2d.py
 
+echo Running basic/distance2mesh.py
+python basic/distance2mesh.py
+
 echo Running basic/clustering.py
 python basic/clustering.py
 
@@ -250,6 +253,9 @@ python simulations/gyroscope2.py
 echo Running simulations/multiple_pendulum.py
 python simulations/multiple_pendulum.py
 
+echo Running simulations/pendulum.py
+python simulations/pendulum.py
+
 echo Running simulations/wave_equation.py
 python simulations/wave_equation.py
 

examples/simulations/README.md

@@ -25,6 +25,8 @@ python example.py  # on mac OSX try 'pythonw' instead
 |    |    |
 | [![mpend](https://user-images.githubusercontent.com/32848391/50738892-db380300-11d8-11e9-807c-fb320c7b7917.gif)](https://github.com/marcomusy/vtkplotter/blob/master/examples/simulations/multiple_pendulum.py)<br/> `multiple_pendulum.py` | Simulation of an elastic multiple pendulums with viscous friction. |
 |    |    |
+| [![pendulum](https://user-images.githubusercontent.com/32848391/55420020-51e56200-5576-11e9-8513-4a5d93913b17.png)](https://github.com/marcomusy/vtkplotter/blob/master/examples/simulations/pendulum.py)<br/> `pendulum.py`       | Visualize the phase space of a simple pendulum (from [3Blue1Brown](https://www.youtube.com/watch?v=p_di4Zn4wz4)). |
+|    |    |
 | [![ruth](https://user-images.githubusercontent.com/32848391/50738891-db380300-11d8-11e9-84c2-0f55be7228f1.gif)](https://github.com/marcomusy/vtkplotter/blob/master/examples/simulations/particle_simulator.py)<br/> `particle_simulator.py` | Simulates interacting charged particles in 3D space. |
 |    |    |
 | [![tunneling1](https://vtkplotter.embl.es/gifs/tunnelling2.gif)](https://github.com/marcomusy/vtkplotter/blob/master/examples/simulations/tunnelling1.py)<br/> `tunnelling1.py`                                                    | Quantum Tunnelling effect using 4th order Runge-Kutta method with arbitrary potential shape. <br>The animation shows the evolution of a particle of relatively well defined momentum (hence undefined position) in a box hitting a potential barrier. The wave function is forced to be zero at the box walls. |

examples/simulations/pendulum.py

@@ -0,0 +1,43 @@
+'''
+Visualize the phase space of a simple pendulum.
+x = starting angle theta
+y = starting angular velocity 
+'''
+# From https://www.youtube.com/watch?v=p_di4Zn4wz4
+# Overview of differential equations | Chapter 1
+# (by 3Blue1Brown)
+
+#Install with:
+# pip install vtkplotter
+
+from vtkplotter import *
+
+# Physical constants
+g = 9.81 # m/s^2
+L = 2    # m
+mu = 0.1 # friction 1/s
+delta_t = 0.01 # Some time step
+t_tot = 50 # seconds total time
+
+# Definition of ODE
+def get_theta_dot_dot(theta, theta_dot):
+    return -mu * theta_dot - (g/L) * sin(theta)
+
+lines = []
+for THETA_0 in arange(0, 3.1415, 0.2):
+    for THETA_DOT_0 in arange(4, 9, 1):
+
+        # Solution to the differential equation
+        theta = THETA_0
+        theta_dot = THETA_DOT_0
+        pts = []
+        for time in arange(0, t_tot, delta_t):
+            theta_dot_dot = get_theta_dot_dot(theta, theta_dot)
+            theta     += theta_dot * delta_t
+            theta_dot += theta_dot_dot * delta_t
+            pts.append([theta, theta_dot])
+
+        l = Line(pts).color(int(THETA_DOT_0))
+        lines.append(l)
+
+show(lines, Text(__doc__), axes=2, bg='white')

setup.py

@@ -2,7 +2,7 @@
 
 setup(
     name='vtkplotter',
-    version='2019.1.2', # change also in vtkplotter/__init__.py and docs/source/conf.py
+    version='2019.1.3', # change also in vtkplotter/__init__.py and docs/source/conf.py
     packages=['vtkplotter'],
     scripts=['bin/vtkplotter', 'bin/vtkconvert'],
     install_requires=['vtk'],
@@ -50,15 +50,15 @@
 
 # # check version number here and in vtkplotter/__init__.py
 
+# git status
+# git commit -a -m 'comment'
+# git push
+
 # git status
 # (sudo apt install twine)
 # (python -m pip install --user --upgrade twine)
 # python setup.py sdist bdist_wheel
 # twine upload dist/vtkplotter-?.?.?.tar.gz -r pypi
-
-# git status
-# git commit -a -m 'comment'
-# git push
 # make release
 
 # pip:

vtkplotter/__init__.py

@@ -46,7 +46,7 @@
 __email__ = "[email protected]"
 __status__ = "dev"
 __website__ = "https://github.com/marcomusy/vtkplotter"
-__version__ = "2019.1.2"  ### defined also above, in setup.py and docs/source/conf.py
+__version__ = "2019.1.3"  ### defined also above, in setup.py and docs/source/conf.py
 
 from vtkplotter.plotter import *
 from vtkplotter.analysis import *

vtkplotter/actors.py

@@ -1081,6 +1081,34 @@ def closestPoint(self, pt, N=1, radius=None, returnIds=False):
         else:
             return np.array(trgp)
 
+
+    def distanceToMesh(self, actor, signed=False, negate=False):
+        '''
+        Computes the (signed) distance from one mesh to another.
+        
+        Example: |distance2mesh.py|_
+        '''
+        poly1 = self.polydata()
+        poly2 = actor.polydata()
+        df = vtk.vtkDistancePolyDataFilter()
+        df.SetInputData(0, poly1)
+        df.SetInputData(1, poly2)
+        if signed:
+            df.SignedDistanceOn()
+        if negate:
+            df.NegateDistanceOn()
+        df.Update()
+
+        scals = df.GetOutput().GetPointData().GetScalars()
+        poly1.GetPointData().AddArray(scals)
+
+        poly1.GetPointData().SetActiveScalars(scals.GetName())
+        rng = scals.GetRange()
+        self.mapper.SetScalarRange(rng[0], rng[1])
+        self.mapper.ScalarVisibilityOn()
+        return self
+
+
     def clone(self, transformed=True):
         """
         Clone a ``Actor(vtkActor)`` and make an exact copy of it.
@@ -1173,7 +1201,8 @@ def mirror(self, axis="x"):
             pass
         else:
             colors.printc("~times Error in mirror(): mirror must be set to x, y, z or n.", c=1)
-            exit()
+            raise RuntimeError()
+
         if axis != "n":
             for j in range(polyCopy.GetNumberOfPoints()):
                 p = [0, 0, 0]
@@ -1797,7 +1826,7 @@ def subdivide(self, N=1, method=0):
             sdf = vtk.vtkButterflySubdivisionFilter()
         else:
             colors.printc("~times Error in subdivide: unknown method.", c="r")
-            exit(1)
+            exit()
         if method != 2:
             sdf.SetNumberOfSubdivisions(N)
         sdf.SetInputData(originalMesh)
@@ -2080,15 +2109,16 @@ def insidePoints(self, points, invert=False, tol=1e-05):
         """Return the sublist of points that are inside a polydata closed surface."""
         poly = self.polydata(True)
         # check if the stl file is closed
-        featureEdge = vtk.vtkFeatureEdges()
-        featureEdge.FeatureEdgesOff()
-        featureEdge.BoundaryEdgesOn()
-        featureEdge.NonManifoldEdgesOn()
-        featureEdge.SetInputData(poly)
-        featureEdge.Update()
-        openEdges = featureEdge.GetOutput().GetNumberOfCells()
-        if openEdges != 0:
-            colors.printc("~lightning Warning: polydata is not a closed surface", c=5)
+
+        #featureEdge = vtk.vtkFeatureEdges()
+        #featureEdge.FeatureEdgesOff()
+        #featureEdge.BoundaryEdgesOn()
+        #featureEdge.NonManifoldEdgesOn()
+        #featureEdge.SetInputData(poly)
+        #featureEdge.Update()
+        #openEdges = featureEdge.GetOutput().GetNumberOfCells()
+        #if openEdges != 0:
+        #    colors.printc("~lightning Warning: polydata is not a closed surface", c=5)
 
         vpoints = vtk.vtkPoints()
         vpoints.SetData(numpy_to_vtk(points, deep=True))