import numpy as np
import scipy as sp
from openpnm.network import Cubic
from openpnm import topotools
from openpnm.utils import logging
logger = logging.getLogger(__name__)


class CubicTemplate(Cubic):
    r"""
    Simple cubic lattice with arbitrary domain shape specified by a template
    image

    The class creates a standard Cubic network the same shape as the provided
    image, then trims pores from the network that are not in the mask.

    Parameters
    ----------
    template : array_like
        The array (image) describing the desired shape of the domain.  All
        locations in the image that are marked as ``True`` are kept while the
        rest of trimmed to yeild the shape.

    spacing : array_like, optional
        The spacing between pore centers in each direction. If not given, then
        [1, 1, 1] is assumed.

    name : string
        An optional name for the object to help identify it.  If not given,
        one will be generated.

    project : OpenPNM Project object, optional
        Each OpenPNM object must be part of a *Project*.  If none is supplied
        then one will be created and this Network will be automatically
        assigned to it.  To create a *Project* use ``openpnm.Project()``.

    Notes
    -----
    The other arguments are the same as ``Cubic`` except that ``shape`` is
    inferred from the ``template`` image.

    See Also
    --------
    The following methods in ``topotools`` can help generate template images:

    template_cylinder_annulus
    template_sphere_shell

    Examples
    --------
    >>> import openpnm as op
    >>> im = op.topotools.template_cylinder_annulus(15, 10, 5)
    >>> pn = op.network.CubicTemplate(template=im)

    And it can be plotted for quick visualization using:

    >>> fig = op.topotools.plot_connections(network=pn)

    .. image:: /../docs/static/images/cubic_template_network.png
        :align: center

    For larger networks and more control over presentation use `Paraview
    <http://www.paraview.org>`_.

    """
    def __init__(self, template, spacing=[1, 1, 1], **kwargs):

        template = np.atleast_3d(template)
        if 'shape' in kwargs:
            del kwargs['shape']
            logger.warning('shape argument ignored, inferred from template')
        super().__init__(shape=template.shape, spacing=spacing, **kwargs)

        coords = np.unravel_index(range(template.size), template.shape)
        self['pore.template_coords'] = np.vstack(coords).T
        self['pore.template_indices'] = self.Ps
        self['pore.drop'] = template.flatten() == 0
        topotools.trim(network=self, pores=self.pores('drop'))
        del self['pore.drop']
        # remove labels pertaining to surface pores, then redo post-trim
        self.clear(mode='labels')
        self['pore.internal'] = True
        self['throat.internal'] = True
        topotools.find_surface_pores(self)

@@ -16,27 +16,8 @@
  np%0A
-import scipy as sp%0A
 from
@@ -2207,21 +2207,23 @@
 arning('
+%22
 shape
+%22
  argumen
@@ -2524,28 +2524,43 @@
     
+topotools.trim(network=
 self
-%5B'
+, 
 pore
-.drop'%5D = 
+s=
 temp
@@ -2577,13 +2577,12 @@
 en()
- 
 ==
- 
 0
+)
 %0A   
@@ -2590,306 +2590,410 @@
     
-topotools.trim(network=self, pores=self.pores('drop'))%0A        del self%5B'pore.drop'%5D%0A        # remove labels pertaining to surface pores, then redo post-trim%0A        self.clear(mode='labels')%0A       
+# Add %22internal_surface%22 label to %22fake%22 surface pores!%0A        ndims = topotools.dimensionality(self).sum()%0A        max_neighbors = 6 if ndims == 3 else 4%0A        num_neighbors = np.diff(self.get_adjacency_matrix(fmt=%22csr%22).indptr)%0A        mask_surface =
  self%5B
-'
+%22
 pore.
-internal'%5D = True%0A        self%5B'throat.internal'%5D = True%0A        topotools.find
+surface%22%5D%0A        mask_internal_surface = (num_neighbors %3C max_neighbors) & ~mask_surface%0A        self.set_label(%22pore.internal
 _surface
 _por
@@ -2992,17 +2992,36 @@
 face
-_
+%22, 
 pores
-(self
+=mask_internal_surface
 )%0A

Enhance CubicTemplate to not call find_surface_pores (takes ages to run) + add "internal_surface" label, indicating "surface" pores that are not on the exterior boundary of the network.