diff --git a/README.md b/README.md
index fea5655..2e13841 100644
--- a/README.md
+++ b/README.md
@@ -6,6 +6,7 @@
 [![Anaconda-Server Badge](https://anaconda.org/conda-forge/puma/badges/version.svg)](https://anaconda.org/conda-forge/puma)
 [![PyPI version](https://badge.fury.io/py/pumapy.svg)](https://badge.fury.io/py/pumapy)
 ![pumapy Tests](https://github.com/nasa/puma/actions/workflows/test-pumapy.yml/badge.svg)
+[![Gitter](https://badges.gitter.im/puma-nasa/community.svg)](https://gitter.im/puma-nasa/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge)
 
 -----
 
@@ -14,7 +15,7 @@
 The stable releases of PuMA can be found on the [official NASA Github repository](https://github.com/nasa/puma),
 whereas the latest development can be found on the
 [development Gitlab repository](https://gitlab.com/jcfergus/puma-dev). Access the
-[PuMA documentation](https://puma-nasa.readthedocs.io) for detailed use of each function or to get started with the tutorial.
+[PuMA documentation](https://puma-nasa.readthedocs.io) for detailed use of each function.
 
 The Porous Microstructure Analysis (PuMA) software has been developed to
 compute effective material properties and perform material response simulations on
@@ -25,7 +26,8 @@ Version 3 includes modules to compute simple morphological properties such as po
 volume fractions, pore diameter, and specific surface area. Additional capabilities include
 the determination of effective thermal and electrical conductivity (both radiative and solid conduction - 
 including the ability to simulate local anisotropy for the latter); effective diffusivity and 
-tortuosity from the continuum to the rarefied regime; techniques to determine the local material orientation, as well as the mechanical properties (elasticity coefficient), and the permeability of a material.
+tortuosity from the continuum to the rarefied regime; techniques to determine the local material orientation, 
+as well as the mechanical properties (elasticity coefficient), and the permeability of a material.
 
 Some examples of microstructures that have been run in the past are shown in the pictures below,
 together with PuMA's software architecture schematic.
diff --git a/python/pumapy/__init__.py b/python/pumapy/__init__.py
index d164a52..f37ae86 100644
--- a/python/pumapy/__init__.py
+++ b/python/pumapy/__init__.py
@@ -31,7 +31,7 @@
 # - git push nasa main
 # - git tag -a v$(python setup.py --version) -m 'INPUT DESCRIPTION'
 # - gh release create v$(python setup.py --version) --target main
-__version__ = "3.1.6"
+__version__ = "3.1.7"
 
 
 # utilities
diff --git a/python/pumapy/io/output.py b/python/pumapy/io/output.py
index b11318e..a87a95e 100644
--- a/python/pumapy/io/output.py
+++ b/python/pumapy/io/output.py
@@ -15,7 +15,7 @@ def export_vti(filename, dict_data, voxel_length=None):
         :param filename: filepath and name
         :type filename: string
         :param dict_data: dictionary setup as {"name1": data1, "name2": data2 ...} containing either Workspaces or ndarrays
-        :type dict_data: dict
+        :type dict_data: dict or Workspace or np.ndarray
         :param voxel_length: with voxel length to give to Numpy arrays (if any)
         :type voxel_length: float
         :return: True if successful, False otherwise.
diff --git a/python/pumapy/utilities/isosurface.py b/python/pumapy/utilities/isosurface.py
index 0724a42..3fbc57b 100644
--- a/python/pumapy/utilities/isosurface.py
+++ b/python/pumapy/utilities/isosurface.py
@@ -1,7 +1,6 @@
 import numpy as np
 from skimage import measure
 import scipy.ndimage as ndimage
-import pyvista as pv
 from pumapy.utilities.workspace import Workspace
 from pumapy.utilities.generic_checks import check_ws_cutoff
 
@@ -46,6 +45,7 @@ def copy(self, other):
         self.values = np.copy(other.values)
 
     def create_mesh(self):
+        import pyvista as pv  # lazily import pyvista to solve Dragonfly's crash
         f = np.zeros((self.faces.shape[0], 4), dtype=np.uint32)
         f[:, 0] = 3
         f[:, 1:] = self.faces
@@ -83,7 +83,7 @@ def compute(self):
         import warnings
         warnings.simplefilter(action='ignore', category=FutureWarning)
         self.tri_mesh.verts, self.tri_mesh.faces, self.tri_mesh.normals, self.tri_mesh.values \
-            = measure.marching_cubes_lewiner(self.matrix, self.iso_value)
+            = measure.marching_cubes(self.matrix, self.iso_value)
 
     def flip_matrix(self):
         self.matrix = np.copy(self.workspace.matrix).astype(np.float32)
diff --git a/python/pumapy/visualization/render.py b/python/pumapy/visualization/render.py
index 3334301..7852b4a 100644
--- a/python/pumapy/visualization/render.py
+++ b/python/pumapy/visualization/render.py
@@ -1,7 +1,6 @@
 from pumapy.utilities.workspace import Workspace
 from pumapy.utilities.isosurface import generate_isosurface
 from pumapy.utilities.logger import print_warning
-import pyvista as pv
 import numpy as np
 import string
 import random
@@ -15,7 +14,7 @@ def render_volume(workspace, cutoff=None, solid_color=None, style='surface', ori
         :param workspace: domain
         :type workspace: Workspace or np.ndarray
         :param cutoff: specifying the values to render
-        :type cutoff: (int, int) or None
+        :type cutoff: (int, int) or (float, float) or None
         :param solid_color: if set to None (default), the material is colored by the matrix's values.
                             Otherwise, a solid color can be specified (e.g. for white (1., 1., 1.))
         :type solid_color: (float, float, float) or None
@@ -120,7 +119,7 @@ def render_orientation(workspace, scale_factor=1., solid_color=None, style='surf
         :type scale_factor: float
         :param solid_color: a solid color for the arrows. Deafult is None, which colors the vectors by their magnitude.
                             To color white, input set solid_color=(1., 1., 1.)
-        :type solid_color: None or (float, float, float)
+        :type solid_color: None or (float, float, float) or str
         :param style: specifying the representation style ('surface', 'edges', 'wireframe', 'points')
         :type style: string
         :param origin: origin of the data as
@@ -276,6 +275,7 @@ def render_contour_multiphase(workspace, cutoffs, solid_colors=None, style='surf
     """
 
     if add_to_plot is None:
+        import pyvista as pv  # lazily import pyvista to solve Dragonfly's crash
         p = pv.Plotter(notebook=notebook)
     else:
         p = add_to_plot
@@ -317,6 +317,8 @@ def __init__(self, existing_plot=None, filter_type=None, workspace=None, cutoff=
         self.voxel_length = 1
         self.scale_factor = scale_factor
 
+        import pyvista as pv  # lazily import pyvista to solve Dragonfly's crash
+
         if isinstance(workspace, Workspace):
             self.voxel_length = workspace.voxel_length
 
@@ -405,6 +407,7 @@ def build_plotter(self):
                 tmp[:, i] = self.array[:, :, :, i].ravel(order='F')
             self.grid.cell_data["scalars"] = np.linalg.norm(self.array, axis=3).ravel(order='F')
             self.grid.cell_data["vectors"] = tmp
+            import pyvista as pv  # lazily import pyvista to solve Dragonfly's crash
             self.filter = self.grid.glyph(orient="vectors", scale="scalars", factor=self.scale_factor, geom=pv.Arrow())
 
         elif self.filter_type == "warp":
diff --git a/setup.py b/setup.py
index 036fa74..b1ee4de 100644
--- a/setup.py
+++ b/setup.py
@@ -62,7 +62,7 @@ def run(self):
     cmdclass={'clean': CleanCommand},
     install_requires=[  # TexGen also required, but it can be installed as add-on
         "numpy",
-        "scikit-image",
+        "scikit-image >=0.17",  # in order to have marching_cubes in the API
         "scipy",
         "matplotlib",
         "pyevtk",