text: import numpy as np
  import openpnm as op
  from tqdm import tqdm
  
  
  __all__ = [
      'plot_connections',
      'plot_coordinates',
      'plot_networkx',
      'plot_vpython',
      'plot_tutorial',
      'plot_network_jupyter',
      'generate_voxel_image',
  ]
  
  
  def plot_connections(network,
                       throats=None,
                       ax=None,
                       size_by=None,
                       color_by=None,
                       cmap='jet',
                       color='b',
                       alpha=1.0,
                       linestyle='solid',
                       linewidth=1,
                       **kwargs):  # pragma: no cover
      r"""
      Produce a 3D plot of the network topology.
  
      This shows how throats connect for quick visualization without having
      to export data to veiw in Paraview.
  
      Parameters
      ----------
      network : GenericNetwork
          The network whose topological connections to plot
      throats : array_like (optional)
          The list of throats to plot if only a sub-sample is desired.  This is
          useful for inspecting a small region of the network.  If no throats are
          specified then all throats are shown.
      fig : Matplotlib figure handle and line property arguments (optional)
          If a ``fig`` is supplied, then the topology will be overlaid on this
          plot.  This makes it possible to combine coordinates and connections,
          and to color throats differently for instance.
      size_by : array_like (optional)
          An ndarray of throat values (e.g. alg['throat.rate']).  These
          values are used to scale the ``linewidth``, so if the lines are too
          thin, then increase ``linewidth``.
      color_by : str or array_like (optional)
          An ndarray of throat values (e.g. alg['throat.rate']).
      cmap : str or cmap object (optional)
          The matplotlib colormap to use if specfying a throat property
          for ``color_by``
      color : str, optional (optional)
          A matplotlib named color (e.g. 'r' for red).
      alpha : float (optional)
          The transparency of the lines, with 1 being solid and 0 being invisible
      linestyle : str (optional)
          Can be one of {'solid', 'dashed', 'dashdot', 'dotted'}.  Default is
          'solid'.
      linewidth : float (optional)
          Controls the thickness of drawn lines.  Is used to scale the thickness
          if ``size_by`` is given. Default is 1. If a value is provided for
          ``size_by`` then they are used to scale the ``linewidth``.
      **kwargs : dict
          All other keyword arguments are passed on to the ``Line3DCollection``
          class of matplotlib, so check their documentation for additional
          formatting options.
  
      Returns
      -------
      lc : LineCollection or Line3DCollection
          Matplotlib object containing the lines representing the throats.
  
      Notes
      -----
      To create a single plot containing both pore coordinates and throats,
      consider creating an empty figure and then pass the ``ax`` object as
      an argument to ``plot_connections`` and ``plot_coordinates``.
      Otherwise, each call to either of these methods creates a new figure.
  
      See Also
      --------
      plot_coordinates
  
      Examples
      --------
      >>> import openpnm as op
      >>> import matplotlib as mpl
      >>> import matplotlib.pyplot as plt
      >>> mpl.use('Agg')
      >>> pn = op.network.Cubic(shape=[10, 10, 3])
      >>> pn.add_boundary_pores()
      >>> Ts = pn.throats('*boundary', mode='not')        # find internal throats
      >>> fig, ax = plt.subplots()                        # create empty figure
      >>> _ = op.topotools.plot_connections(network=pn,
      ...                                   throats=Ts)   # plot internal throats
      >>> Ts = pn.throats('*boundary')                    # find boundary throats
      >>> _ = op.topotools.plot_connections(network=pn,
      ...                                   throats=Ts,
      ...                                   ax=ax,
      ...                                   color='r')    # plot boundary throats in red
  
      """
      import matplotlib.pyplot as plt
      from matplotlib import cm
      from matplotlib import colors as mcolors
      from mpl_toolkits.mplot3d import Axes3D
      from matplotlib.collections import LineCollection
      from mpl_toolkits.mplot3d.art3d import Line3DCollection
      from openpnm.topotools import dimensionality
  
      Ts = network.Ts if throats is None else network._parse_indices(throats)
      dim = dimensionality(network)
      ThreeD = True if dim.sum() == 3 else False
      # Add a dummy axis for 1D networks
      if dim.sum() == 1:
          dim[np.argwhere(~dim)[0]] = True
  
      if "fig" in kwargs.keys():
          raise Exception("'fig' argument is deprecated, use 'ax' instead.")
      if ax is None:
          fig, ax = plt.subplots()
      else:
          # The next line is necessary if ax was created using plt.subplots()
          fig, ax = ax.get_figure(), ax.get_figure().gca()
      if ThreeD and ax.name != '3d':
          fig.delaxes(ax)
          ax = fig.add_subplot(111, projection='3d')
  
      # Collect coordinates
      Ps = np.unique(network['throat.conns'][Ts])
      X, Y, Z = network['pore.coords'][Ps].T
      xyz = network["pore.coords"][:, dim]
      P1, P2 = network["throat.conns"][Ts].T
      throat_pos = np.column_stack((xyz[P1], xyz[P2])).reshape((Ts.size, 2, dim.sum()))
  
      # Deal with optional style related arguments
      if 'c' in kwargs.keys():
          color = kwargs.pop('c')
      color = mcolors.to_rgb(color) + tuple([alpha])
      # Override colors with color_by if given
      if color_by is not None:
          color = cm.get_cmap(name=cmap)(color_by / color_by.max())
          color[:, 3] = alpha
      if size_by is not None:
          linewidth = size_by / size_by.max() * linewidth
  
      if ThreeD:
          lc = Line3DCollection(throat_pos, colors=color, cmap=cmap,
                                linestyles=linestyle, linewidths=linewidth,
                                antialiaseds=np.ones_like(network.Ts), **kwargs)
      else:
          lc = LineCollection(throat_pos, colors=color, cmap=cmap,
                              linestyles=linestyle, linewidths=linewidth,
                              antialiaseds=np.ones_like(network.Ts), **kwargs)
      ax.add_collection(lc)
  
      _scale_axes(ax=ax, X=X, Y=Y, Z=Z)
      _label_axes(ax=ax, X=X, Y=Y, Z=Z)
      fig.tight_layout()
  
      return lc
  
  
  def plot_coordinates(network,
                       pores=None,
                       ax=None,
                       size_by=None,
                       color_by=None,
                       cmap='jet',
                       color='r',
                       alpha=1.0,
                       marker='o',
                       markersize=10,
                       **kwargs):  # pragma: no cover
      r"""
      Produce a 3D plot showing specified pore coordinates as markers.
  
      Parameters
      ----------
      network : GenericNetwork
          The network whose topological connections to plot.
      pores : array_like (optional)
          The list of pores to plot if only a sub-sample is desired. This is
          useful for inspecting a small region of the network. If no pores
          are specified then all are shown.
      ax : Matplotlib axis handle
          If ``ax`` is supplied, then the coordinates will be overlaid.
          This enables the plotting of multiple different sets of pores as
          well as throat connections from ``plot_connections``.
      size_by : str or array_like
          An ndarray of pore values (e.g. alg['pore.concentration']). These
          values are normalized by scaled by ``markersize``.
      color_by : str or array_like
          An ndarray of pore values (e.g. alg['pore.concentration']).
      cmap : str or cmap object
          The matplotlib colormap to use if specfying a pore property
          for ``color_by``
      color : str
          A matplotlib named color (e.g. 'r' for red).
      alpha : float
          The transparency of the lines, with 1 being solid and 0 being invisible
      marker : 's'
          The marker to use.  The default is a circle.  Options are explained
          `here <https://matplotlib.org/3.2.1/api/markers_api.html>`_
      markersize : scalar
          Controls size of marker, default is 1.0.  This value is used to scale
          the ``size_by`` argument if given.
      **kwargs
          All other keyword arguments are passed on to the ``scatter``
          function of matplotlib, so check their documentation for additional
          formatting options.
  
      Returns
      -------
      pc : PathCollection
          Matplotlib object containing the markers representing the pores.
  
      Notes
      -----
      To create a single plot containing both pore coordinates and throats,
      consider creating an empty figure and then pass the ``ax`` object as
      an argument to ``plot_connections`` and ``plot_coordinates``.
      Otherwise, each call to either of these methods creates a new figure.
  
      See Also
      --------
      plot_connections
  
      Examples
      --------
      >>> import openpnm as op
      >>> import matplotlib as mpl
      >>> import matplotlib.pyplot as plt
      >>> mpl.use('Agg')
      >>> pn = op.network.Cubic(shape=[10, 10, 3])
      >>> pn.add_boundary_pores()
      >>> Ps = pn.pores('internal')                       # find internal pores
      >>> fig, ax = plt.subplots()                        # create empty figure
      >>> _ = op.topotools.plot_coordinates(network=pn,
      ...                                   pores=Ps,
      ...                                   color='b',
      ...                                   ax=ax)        # plot internal pores
      >>> Ps = pn.pores('*boundary')                      # find boundary pores
      >>> _ = op.topotools.plot_coordinates(network=pn,
      ...                                   pores=Ps,
      ...                                   color='r',
      ...                                   ax=ax)        # plot boundary pores in red
  
      """
      import matplotlib.pyplot as plt
      from matplotlib import cm
      from mpl_toolkits.mplot3d import Axes3D
      from openpnm.topotools import dimensionality
  
      Ps = network.Ps if pores is None else network._parse_indices(pores)
  
      dim = dimensionality(network)
      ThreeD = True if dim.sum() == 3 else False
      # Add a dummy axis for 1D networks
      if dim.sum() == 1:
          dim[np.argwhere(~dim)[0]] = True
      # Add 2 dummy axes for 0D networks (1 pore only)
      if dim.sum() == 0:
          dim[[0, 1]] = True
  
      if "fig" in kwargs.keys():
          raise Exception("'fig' argument is deprecated, use 'ax' instead.")
      if ax is None:
          fig, ax = plt.subplots()
      else:
          # The next line is necessary if ax was created using plt.subplots()
          fig, ax = ax.get_figure(), ax.get_figure().gca()
      if ThreeD and ax.name != '3d':
          fig.delaxes(ax)
          ax = fig.add_subplot(111, projection='3d')
  
      # Collect specified coordinates
      X, Y, Z = network['pore.coords'][Ps].T
      # The bounding box for fig is the entire ntwork (to fix the problem with
      # overwriting figures' axes lim)
      Xl, Yl, Zl = network['pore.coords'].T
  
      # Parse formatting kwargs
      if 'c' in kwargs.keys():
          color = kwargs.pop('c')
      if 's' in kwargs.keys():
          markersize = kwargs.pop('s')
      if color_by is not None:
          color = cm.get_cmap(name=cmap)(color_by / color_by.max())
      if size_by is not None:
          markersize = size_by / size_by.max() * markersize
  
      if ThreeD:
          sc = ax.scatter(X, Y, Z, c=color, s=markersize, marker=marker, alpha=alpha, **kwargs)
          _scale_axes(ax=ax, X=Xl, Y=Yl, Z=Zl)
      else:
          _X, _Y = np.column_stack((X, Y, Z))[:, dim].T
          sc = ax.scatter(_X, _Y, c=color, s=markersize, marker=marker, alpha=alpha, **kwargs)
          _scale_axes(ax=ax, X=Xl, Y=Yl, Z=np.zeros_like(Yl))
  
      _label_axes(ax=ax, X=Xl, Y=Yl, Z=Zl)
      fig.tight_layout()
  
      return sc
  
  
  def _label_axes(ax, X, Y, Z):
      labels = ["X", "Y", "Z"]
      dim = np.zeros(3, dtype=bool)
      for i, arr in enumerate([X, Y, Z]):
          if np.unique(arr).size > 1:
              dim[i] = True
      # Add a dummy axis for 1D networks
      if dim.sum() == 1:
          dim[np.argwhere(~dim)[0]] = True
      # Add 2 dummy axes for 0D networks (1 pore only)
      if dim.sum() == 0:
          dim[[0, 1]] = True
      dim_idx = np.argwhere(dim).squeeze()
      ax.set_xlabel(labels[dim_idx[0]])
      ax.set_ylabel(labels[dim_idx[1]])
      if hasattr(ax, "set_zlim"):
          ax.set_zlabel("Z")
  
  
  def _scale_axes(ax, X, Y, Z):
      max_range = np.ptp([X, Y, Z]).max() / 2
      mid_x = (X.max() + X.min()) * 0.5
      mid_y = (Y.max() + Y.min()) * 0.5
      mid_z = (Z.max() + Z.min()) * 0.5
      ax.set_xlim(mid_x - max_range, mid_x + max_range)
      ax.set_ylim(mid_y - max_range, mid_y + max_range)
      if hasattr(ax, "set_zlim"):
          ax.set_zlim(mid_z - max_range, mid_z + max_range)
      else:
          ax.axis("equal")
  
  
  def plot_networkx(network,
                    plot_throats=True,
                    labels=None,
                    colors=None,
                    scale=1,
                    ax=None,
                    alpha=1.0):  # pragma: no cover
      r"""
      Creates a pretty 2d plot for 2d OpenPNM networks.
  
      Parameters
      ----------
      network : GenericNetwork
      plot_throats : bool, optional
          Plots throats as well as pores, if True.
      labels : list, optional
          List of OpenPNM labels
      colors : list, optional
          List of corresponding colors to the given `labels`.
      scale : float, optional
          Scale factor for size of pores.
      ax : matplotlib.Axes, optional
          Matplotlib axes object
      alpha: float, optional
          Transparency value, 1 is opaque and 0 is transparent
  
      """
      import matplotlib.pyplot as plt
      from matplotlib.collections import PathCollection
      from networkx import Graph, draw_networkx_nodes, draw_networkx_edges
      from openpnm.topotools import dimensionality
  
      dims = dimensionality(network)
      if dims.sum() > 2:
          raise Exception("NetworkX plotting only works for 2D networks.")
      temp = network['pore.coords'].T[dims].squeeze()
      if dims.sum() == 1:
          x = temp
          y = np.zeros_like(x)
      if dims.sum() == 2:
          x, y = temp
      try:
          node_size = scale * network['pore.diameter']
      except KeyError:
          node_size = np.ones_like(x) * scale * 0.5
      G = Graph()
      pos = {network.Ps[i]: [x[i], y[i]] for i in range(network.Np)}
      if not np.isfinite(node_size).all():
          node_size[~np.isfinite(node_size)] = np.nanmin(node_size)
      node_color = np.array(['k'] * len(network.Ps))
  
      if labels:
          if not isinstance(labels, list):
              labels = [labels]
          if not isinstance(colors, list):
              colors = [colors]
          if len(labels) != len(colors):
              raise Exception('len(colors) must be equal to len(labels)!')
          for label, color in zip(labels, colors):
              node_color[network.pores(label)] = color
  
      if ax is None:
          fig, ax = plt.subplots()
      ax.set_aspect('equal', adjustable='datalim')
      offset = node_size.max() * 0.5
      ax.set_xlim((x.min() - offset, x.max() + offset))
      ax.set_ylim((y.min() - offset, y.max() + offset))
      ax.axis("off")
  
      # Keep track of already plotted nodes
      temp = [id(item) for item in ax.collections if isinstance(item, PathCollection)]
  
      # Plot pores
      gplot = draw_networkx_nodes(G, ax=ax, pos=pos, nodelist=network.Ps.tolist(),
                                  alpha=alpha, node_color=node_color, edgecolors=node_color,
                                  node_size=node_size)
      # (Optionally) Plot throats
      if plot_throats:
          draw_networkx_edges(G, pos=pos, edge_color='k', alpha=alpha,
                              edgelist=network['throat.conns'].tolist(), ax=ax)
  
      spi = 2700  # 1250 was obtained by trial and error
      figwidth, figheight = ax.get_figure().get_size_inches()
      figsize_ratio = figheight / figwidth
      data_ratio = ax.get_data_ratio()
      corr = min(figsize_ratio / data_ratio, 1)
      xrange = np.ptp(ax.get_xlim())
      markersize = np.atleast_1d((corr*figwidth)**2 / xrange**2 * node_size**2 * spi)
      for item in ax.collections:
          if isinstance(item, PathCollection) and id(item) not in temp:
              item.set_sizes(markersize)
  
      return gplot
  
  
  def plot_vpython(network,
                   Psize='pore.diameter',
                   Tsize='throat.diameter',
                   Pcolor=None,
                   Tcolor=None,
                   cmap='jet',
                   **kwargs):  # pragma: no cover
      r"""
      Quickly visualize a network in 3D using VPython.
  
      Parameters
      ----------
      network : GenericNetwork
          The network to visualize.
      Psize : str (default = 'pore.diameter')
          The dictionary key pointing to the pore property by which sphere
          diameters should be scaled
      Tsize : str (default = 'throat.diameter')
          The dictionary key pointing to the throat property by which cylinder
          diameters should be scaled
      Pcolor : str
          The dictionary key pointing to the pore property which will control
          the sphere colors.  The default is None, which results in a bright
          red for all pores.
      Tcolor : str
          The dictionary key pointing to the throat property which will control
          the cylinder colors.  The default is None, which results in a unform
          pale blue for all throats.
      cmap : str or Matplotlib colormap object (default is 'jet')
          The color map to use when converting pore and throat properties to
          RGB colors.  Can either be a string indicating which color map to
          fetch from matplotlib.cmap, or an actual cmap object.
      kwargs : dict
          Any additional kwargs that are received are passed to the VPython
          ``canvas`` object.  Default options are:
  
          *'height' = 500* - Height of canvas
  
          *'width' = 800* - Width of canvas
  
          *'background' = [0, 0, 0]* - Sets the background color of canvas
  
          *'ambient' = [0.2, 0.2, 0.3]* - Sets the brightness of lighting
  
      Returns
      -------
      canvas : VPython Canvas object
          The canvas object containing the generated scene. The object has
          several useful methods.
  
      Notes
      -----
      **Important**
  
      a) This does not work in Spyder. It should only be called from a
      Jupyter Notebook.
  
      b) This is only meant for relatively small networks. For proper
      visualization use Paraview.
  
      """
      import matplotlib.pyplot as plt
  
      try:
          from vpython import canvas, vec, sphere, cylinder
      except ModuleNotFoundError:
          raise Exception('VPython must be installed to use this function')
  
      if isinstance(cmap, str):
          cmap = getattr(plt.cm, cmap)
  
      if Pcolor is None:
          Pcolor = [vec(230/255, 57/255, 0/255)]*network.Np
      else:
          a = cmap(network[Pcolor]/network[Pcolor].max())
          Pcolor = [vec(row[0], row[1], row[2]) for row in a]
  
      if Tcolor is None:
          Tcolor = [vec(51/255, 153/255, 255/255)]*network.Nt
      else:
          a = cmap(network[Tcolor]/network[Tcolor].max())
          Tcolor = [vec(row[0], row[1], row[2]) for row in a]
  
      # Set default values for canvas properties
      if 'background' not in kwargs.keys():
          kwargs['background'] = vec(1.0, 1.0, 1.0)
      if 'height' not in kwargs.keys():
          kwargs['height'] = 500
      if 'width' not in kwargs.keys():
          kwargs['width'] = 800
      # Parse any given values for canvas properties
      for item in kwargs.keys():
          try:
              kwargs[item] = vec(*kwargs[item])
          except TypeError:
              pass
      scene = canvas(title=network.name, **kwargs)
  
      for p in network.Ps:
          r = network[Psize][p]/2
          xyz = network['pore.coords'][p]
          c = Pcolor[p]
          sphere(pos=vec(*xyz), radius=r, color=c,
                 shininess=.5)
  
      for t in network.Ts:
          head = network['throat.endpoints.head'][t]
          tail = network['throat.endpoints.tail'][t]
          v = tail - head
          r = network[Tsize][t]
          L = np.sqrt(np.sum((head-tail)**2))
          c = Tcolor[t]
          cylinder(pos=vec(*head), axis=vec(*v), opacity=1, size=vec(L, r, r),
                   color=c)
  
      return scene
  
  
  def plot_tutorial(network,
                    font_size=12,
                    line_width=2,
                    node_color='b',
                    edge_color='r',
                    node_size=500):  # pragma: no cover
      r"""
      Generate a network plot suitable for tutorials and explanations.
  
      Parameters
      ----------
      network : GenericNetwork
          The network to plot, should be 2D, since the z-coordinate will be
          ignored.
      font_size : int
          Size of font to use for labels.
      line_width : int
          Thickness of edge lines and node borders.
      node_color : str
          Color of node border.
      edge_color : str
          Color of edge lines.
      node_size : int
          Size of node circle.
  
      Returns
      -------
      g : NetworkX plot object
  
      """
      import networkx as nx
      import matplotlib.pyplot as plt
      from openpnm.io import to_networkx
  
      G = to_networkx(network=network)
      pos = {i: network['pore.coords'][i, 0:2] for i in network.Ps}
      labels = {i: i for i in network.Ps}
      edge_labels = {tuple(network['throat.conns'][i, :]): i for i in network.Ts}
  
      gplot = nx.draw_networkx_nodes(G, pos,
                                     node_size=node_size,
                                     node_color='w',
                                     edgecolors=node_color,
                                     linewidths=line_width)
      nx.draw_networkx_edges(
          G, pos, width=line_width, edge_color=edge_color)
      nx.draw_networkx_labels(
          G, pos, labels=labels, font_size=font_size, font_color='k')
      nx.draw_networkx_edge_labels(
          G, pos, edge_labels=edge_labels, font_size=font_size, font_color='k')
  
      # Prettify the figure (margins, etc.)
      plt.axis('off')
      ax = plt.gca()
      ax.margins(0.1, 0.1)
      ax.set_aspect("equal")
      fig = plt.gcf()
      fig.tight_layout()
      dims = op.topotools.dimensionality(network)
      xy_range = network.coords.ptp(axis=0)[dims]
      aspect_ratio = xy_range[0] / xy_range[1]
      fig.set_size_inches(5, 5 / aspect_ratio)
  
      return gplot
  
  
  def plot_network_jupyter(network,
                           node_color=0,
                           edge_color=0,
                           node_size=1,
                           node_scale=20,
                           edge_scale=5,
                           colormap='viridis'):
      r"""
      Visualize a network in 3D using Plotly.
  
      The pores and throats are scaled and colored by their properties.
      The final figure can be rotated and zoomed.
  
      Parameters
      ----------
      network : GenericNetwork
          The network to visualize
      node_color : ndarray
          An array of values used for coloring the pores. If not given, the
          lowest value of the employed colormap is assigned to all markers.
      edge_color : ndarray
          An array of values used for coloring the throats. If not given, the
          lowest value of the employed colormap is assigned to all lines.
      node_size : ndarray
          An array of values controlling the size of the markers.  If not given
          all markers will be the same size
      node_scale : scalar
          A scaler to resize the markers
      edge_scale : scalar
          A scaler to the line thickness
      colormap : str
          The colormap to use
  
      Returns
      -------
      fig : Plotly graph object
          The graph object containing the generated plots. The object has
          several useful methods.
  
      Notes
      -----
      **Important**
  
      a) This does not work in Spyder. It should only be called from a
      Jupyter Notebook.
  
      b) This is only meant for relatively small networks. For proper
      visualization use Paraview.
  
      """
      try:
          import plotly.graph_objects as go
      except ImportError:
          raise Exception('Plotly is not installed.'
                          'Please install Plotly using "pip install plotly"')
  
      # Get xyz coords for points
      x_nodes, y_nodes, z_nodes = network.coords.T
  
      node_size = np.ones(network.Np)*node_size
      node_color = np.ones(network.Np)*node_color
      edge_color = np.ones(network.Nt)*edge_color
  
      node_labels = [str(i)+ ': ' + str(x) for i, x in enumerate(zip(node_size, node_color))]
      edge_labels = [str(i)+ ': ' + str(x) for i, x in enumerate(edge_color)]
  
      # Create edges and nodes coordinates
      N = network.Nt*3
  
      x_edges = np.zeros(N)
      x_edges[np.arange(0, N, 3)] = network.coords[network.conns[:, 0]][:, 0]
      x_edges[np.arange(1, N, 3)] = network.coords[network.conns[:, 1]][:, 0]
      x_edges[np.arange(2, N, 3)] = np.nan
  
      y_edges = np.zeros(network.Nt*3)
      y_edges[np.arange(0, N, 3)] = network.coords[network.conns[:, 0]][:, 1]
      y_edges[np.arange(1, N, 3)] = network.coords[network.conns[:, 1]][:, 1]
      y_edges[np.arange(2, N, 3)] = np.nan
  
      z_edges = np.zeros(network.Nt*3)
      z_edges[np.arange(0, N, 3)] = network.coords[network.conns[:, 0]][:, 2]
      z_edges[np.arange(1, N, 3)] = network.coords[network.conns[:, 1]][:, 2]
      z_edges[np.arange(2, N, 3)] = np.nan
  
      # Create plotly's Scatter3d object for pores and throats
      trace_edges = go.Scatter3d(x=x_edges,
                                 y=y_edges,
                                 z=z_edges,
                                 mode='lines',
                                 line=dict(color=edge_color,
                                           width=edge_scale,
                                           colorscale=colormap),
                                 text=edge_labels, hoverinfo='text')
  
      trace_nodes = go.Scatter3d(x=x_nodes,
                                 y=y_nodes,
                                 z=z_nodes,
                                 mode='markers',
                                 marker=dict(symbol='circle',
                                             size=node_size*node_scale,
                                             color=node_color,
                                             colorscale=colormap,
                                             line=dict(color='black', width=0.5)),
                                 text=node_labels, hoverinfo='text')
  
      axis = dict(showbackground=False,
                  showline=False,
                  zeroline=False,
                  showgrid=False,
                  showticklabels=False,
                  title='')
  
      layout = go.Layout(width=650,
                         height=625,
                         showlegend=False,
                         scene=dict(xaxis=dict(axis),
                                    yaxis=dict(axis),
                                    zaxis=dict(axis),),
                         margin=dict(t=100),
                         hovermode='closest')
  
      data = [trace_edges, trace_nodes]
      fig = go.Figure(data=data, layout=layout)
      return fig
  
  
  def _generate_voxel_image(network, pore_shape, throat_shape, max_dim=200):
      r"""
      Generates a 3d numpy array from an OpenPNM network
  
      Parameters
      ----------
      network : OpenPNM GenericNetwork
          Network from which voxel image is to be generated
      pore_shape : str
          Shape of pores in the network, valid choices are "sphere", "cube"
      throat_shape : str
          Shape of throats in the network, valid choices are "cylinder", "cuboid"
      max_dim : int
          Number of voxels in the largest dimension of the network
  
      Returns
      -------
      im : ndarray
          Voxelated image corresponding to the given pore network model
  
      Notes
      -----
      (1) The generated voxel image is labeled with 0s, 1s and 2s signifying
      solid phase, pores, and throats respectively.
  
      """
      from skimage.morphology import cube, ball
      from porespy.tools import overlay, insert_cylinder
      xyz = network["pore.coords"]
      cn = network["throat.conns"]
  
      # Distance bounding box from the network by a fixed amount
      delta = network["pore.diameter"].mean() / 2
      if isinstance(network, op.network.Cubic):
          try:
              delta = op.topotools.get_spacing(network).mean() / 2
          except AttributeError:
              delta = network.spacing.mean() / 2
  
      # Shift everything to avoid out-of-bounds
      extra_clearance = int(max_dim * 0.05)
  
      # Transform points to satisfy origin at (0, 0, 0)
      xyz0 = xyz.min(axis=0) - delta
      xyz += -xyz0
      res = (xyz.ptp(axis=0).max() + 2 * delta) / max_dim
      shape = np.rint((xyz.max(axis=0) + delta) / res).astype(int) + 2 * extra_clearance
  
      # Transforming from real coords to matrix coords
      xyz = np.rint(xyz / res).astype(int) + extra_clearance
      pore_radi = np.rint(network["pore.diameter"] * 0.5 / res).astype(int)
      throat_radi = np.rint(network["throat.diameter"] * 0.5 / res).astype(int)
  
      im_pores = np.zeros(shape, dtype=np.uint8)
      im_throats = np.zeros_like(im_pores)
  
      if pore_shape == "cube":
          pore_elem = cube
          rp = pore_radi * 2 + 1  # +1 since num_voxel must be odd
          rp_max = int(2 * round(delta / res)) + 1
      if pore_shape == "sphere":
          pore_elem = ball
          rp = pore_radi
          rp_max = int(round(delta / res))
      if throat_shape == "cuboid":
          raise Exception("Not yet implemented, try 'cylinder'.")
  
      # Generating voxels for pores
      for i, pore in enumerate(tqdm(network.Ps)):
          elem = pore_elem(rp[i])
          try:
              im_pores = overlay(im1=im_pores, im2=elem, c=xyz[i])
          except ValueError:
              elem = pore_elem(rp_max)
              im_pores = overlay(im1=im_pores, im2=elem, c=xyz[i])
      # Get rid of pore overlaps
      im_pores[im_pores > 0] = 1
  
      # Generating voxels for throats
      for i, throat in enumerate(tqdm(network.Ts)):
          try:
              im_throats = insert_cylinder(
                  im_throats, r=throat_radi[i], xyz0=xyz[cn[i, 0]], xyz1=xyz[cn[i, 1]])
          except ValueError:
              im_throats = insert_cylinder(
                  im_throats, r=rp_max, xyz0=xyz[cn[i, 0]], xyz1=xyz[cn[i, 1]])
      # Get rid of throat overlaps
      im_throats[im_throats > 0] = 1
  
      # Subtract pore-throat overlap from throats
      im_throats = (im_throats.astype(bool) * ~im_pores.astype(bool)).astype(np.uint8)
      im = im_pores * 1 + im_throats * 2
  
      return im[extra_clearance:-extra_clearance,
                extra_clearance:-extra_clearance,
                extra_clearance:-extra_clearance]
  
  
  def generate_voxel_image(network, pore_shape="sphere", throat_shape="cylinder",
                           max_dim=None, rtol=0.1):
      r"""
      Generate a voxel image from an GenericNetwork
  
      Parameters
      ----------
      network : OpenPNM GenericNetwork
          Network from which voxel image is to be generated
      pore_shape : str
          Shape of pores in the network, valid choices are "sphere", "cube"
      throat_shape : str
          Shape of throats in the network, valid choices are "cylinder", "cuboid"
      max_dim : int
          Number of voxels in the largest dimension of the network
      rtol : float
          Stopping criteria for finding the smallest voxel image such that
          further increasing the number of voxels in each dimension by 25% would
          improve the predicted porosity of the image by less that ``rtol``
  
      Returns
      -------
      im : ndarray
          Voxelated image corresponding to the given pore network model
  
      Notes
      -----
      (1) The generated voxelated image is labeled with 0s, 1s and 2s signifying
      solid phase, pores, and throats respectively.
  
      (2) If max_dim is not provided, the method calculates it such that the
      further increasing it doesn't change porosity by much.
  
      """
      # If max_dim is provided, generate voxel image using max_dim
      if max_dim is not None:
          return _generate_voxel_image(
              network, pore_shape, throat_shape, max_dim=max_dim)
      max_dim = 200
      # If max_dim is not provided, find best max_dim that predicts porosity
      err = 100
      eps_old = 200
      while err > rtol:
          im = _generate_voxel_image(
              network, pore_shape, throat_shape, max_dim=max_dim)
          eps = im.astype(bool).sum() / np.prod(im.shape)
          err = abs(1 - eps / eps_old)
          eps_old = eps
          max_dim = int(max_dim * 1.25)
      return im
  