def add_boundary_regions(regions=None, faces=['front', 'back', 'left',
                                              'right', 'top', 'bottom']):
    r"""
    Given an image partitioned into regions, pads specified faces with new
    regions

    Parameters
    ----------
    regions : ND-array
        An image of the pore space partitioned into regions and labeled
    faces : list of strings
        The faces of ``regions`` which should have boundaries added.  Options
        are:

        *'right'* - Adds boundaries to the x=0 face (``im[0, :, :]``)

        *'left'* - Adds boundaries to the x=X face (``im[-1, :, :]``)

        *'front'* - Adds boundaries to the y=0 face (``im[:, ), :]``)

        *'back'* - Adds boundaries to the x=0 face (``im[:, -1, :]``)

        *'bottom'* - Adds boundaries to the x=0 face (``im[:, :, 0]``)

        *'top'* - Adds boundaries to the x=0 face (``im[:, :, -1]``)

        The default is all faces.

    Returns
    -------
    image : ND-array
        A copy of ``regions`` with the specified boundaries added, so will be
        slightly larger in each direction where boundaries were added.

    """
    # -------------------------------------------------------------------------
    # Edge pad segmentation and distance transform
    if faces is not None:
        regions = sp.pad(regions, 1, 'edge')
        # ---------------------------------------------------------------------
        if regions.ndim == 3:
            # Remove boundary nodes interconnection
            regions[:, :, 0] = regions[:, :, 0] + regions.max()
            regions[:, :, -1] = regions[:, :, -1] + regions.max()
            regions[0, :, :] = regions[0, :, :] + regions.max()
            regions[-1, :, :] = regions[-1, :, :] + regions.max()
            regions[:, 0, :] = regions[:, 0, :] + regions.max()
            regions[:, -1, :] = regions[:, -1, :] + regions.max()
            regions[:, :, 0] = (~find_boundaries(regions[:, :, 0],
                                                 mode='outer')) * regions[:, :, 0]
            regions[:, :, -1] = (~find_boundaries(regions[:, :, -1],
                                                  mode='outer')) * regions[:, :, -1]
            regions[0, :, :] = (~find_boundaries(regions[0, :, :],
                                                 mode='outer')) * regions[0, :, :]
            regions[-1, :, :] = (~find_boundaries(regions[-1, :, :],
                                                  mode='outer')) * regions[-1, :, :]
            regions[:, 0, :] = (~find_boundaries(regions[:, 0, :],
                                                 mode='outer')) * regions[:, 0, :]
            regions[:, -1, :] = (~find_boundaries(regions[:, -1, :],
                                                  mode='outer')) * regions[:, -1, :]
            # -----------------------------------------------------------------
            regions = sp.pad(regions, 2, 'edge')

            # Remove unselected faces
            if 'front' not in faces:
                regions = regions[:, 3:, :]  # y
            if 'back' not in faces:
                regions = regions[:, :-3, :]
            if 'left' not in faces:
                regions = regions[3:, :, :]  # x
            if 'right' not in faces:
                regions = regions[:-3, :, :]
            if 'bottom' not in faces:
                regions = regions[:, :, 3:]  # z
            if 'top' not in faces:
                regions = regions[:, :, :-3]

        elif regions.ndim == 2:
            # Remove boundary nodes interconnection
            regions[0, :] = regions[0, :] + regions.max()
            regions[-1, :] = regions[-1, :] + regions.max()
            regions[:, 0] = regions[:, 0] + regions.max()
            regions[:, -1] = regions[:, -1] + regions.max()
            regions[0, :] = (~find_boundaries(regions[0, :],
                                              mode='outer')) * regions[0, :]
            regions[-1, :] = (~find_boundaries(regions[-1, :],
                                               mode='outer')) * regions[-1, :]
            regions[:, 0] = (~find_boundaries(regions[:, 0],
                                              mode='outer')) * regions[:, 0]
            regions[:, -1] = (~find_boundaries(regions[:, -1],
                                               mode='outer')) * regions[:, -1]
            # -----------------------------------------------------------------
            regions = sp.pad(regions, 2, 'edge')

            # Remove unselected faces
            if 'left' not in faces:
                regions = regions[3:, :]  # x
            if 'right' not in faces:
                regions = regions[:-3, :]
            if 'front' not in faces and 'bottom' not in faces:
                regions = regions[:, 3:]  # y
            if 'back' not in faces and 'top' not in faces:
                regions = regions[:, :-3]
        else:
            print('add_boundary_regions works only on 2D and 3D images')
        # ---------------------------------------------------------------------
        # Make labels contiguous
        regions = make_contiguous(regions)
    else:
        regions = regions

    return regions
r"""
    Given an image partitioned into regions, pads specified faces with new
    regions

    Parameters
    ----------
    regions : ND-array
        An image of the pore space partitioned into regions and labeled
    faces : list of strings
        The faces of ``regions`` which should have boundaries added.  Options
        are:

        *'right'* - Adds boundaries to the x=0 face (``im[0, :, :]``)

        *'left'* - Adds boundaries to the x=X face (``im[-1, :, :]``)

        *'front'* - Adds boundaries to the y=0 face (``im[:, ), :]``)

        *'back'* - Adds boundaries to the x=0 face (``im[:, -1, :]``)

        *'bottom'* - Adds boundaries to the x=0 face (``im[:, :, 0]``)

        *'top'* - Adds boundaries to the x=0 face (``im[:, :, -1]``)

        The default is all faces.

    Returns
    -------
    image : ND-array
        A copy of ``regions`` with the specified boundaries added, so will be
        slightly larger in each direction where boundaries were added.