Source code for ClearMap.ImageProcessing.GreyReconstruction

"""
GreyReconstruction
==================

This morphological reconstruction routine was adapted from 
`CellProfiler <http://www.cellprofiler.org>`_.

Authors
-------
Original author: Lee Kamentsky, Massachusetts Institute of Technology
Modifed by Christoph Kirst for ClearMap integration.
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'


import numpy as np

from skimage.filters._rank_order import rank_order


import ClearMap.ImageProcessing.Filter.StructureElement as se

import ClearMap.Utils.Timer as tmr
import ClearMap.Utils.HierarchicalDict as hdict


###############################################################################
### Grey reconstruction
###############################################################################


[docs]
def reconstruct(seed, mask = None, method = 'dilation', selem = None, offset = None):
  """Performs a morphological reconstruction of an image.
  
  Arguments
  ---------
  seed : array
    Seed image to be dilated or eroded.
  mask : array
    Maximum (dilation) / minimum (erosion) allowed
  method : {'dilation'|'erosion'}
    The method to use.
  selem : array
    Structuring element.
  offset : array or None
    The offset of the structuring element, None is centered.

  Returns
  -------
  reconstructed : array
    Result of morphological reconstruction.
  
  Note
  ----
  Reconstruction uses a seed image, which specifies the values
  to dilate and a mask image that gives the maximum allowed dilated value at
  each pixel.
  
  The algorithm is taken from [1]_. Applications for greyscale 
  reconstruction are discussed in [2]_ and [3]_.
  
  Effectively operates on 2d images.
  
  Reference:
  
  .. [1] Robinson, "Efficient morphological reconstruction: a downhill
         filter", Pattern Recognition Letters 25 (2004) 1759-1767.
  .. [2] Vincent, L., "Morphological Grayscale Reconstruction in Image
         Analysis: Applications and Efficient Algorithms", IEEE Transactions
         on Image Processing (1993)
  .. [3] Soille, P., "Morphological Image Analysis: Principles and
         Applications", Chapter 6, 2nd edition (2003), ISBN 3540429883.
  """
  if mask is None:
    mask = seed.copy();
  
  if seed.shape != mask.shape:
    raise ValueError('Seed shape % and mask shape %r do not match' % (seed.shape, mask.shape))
  
  if method == 'dilation' and np.any(seed > mask):
    raise ValueError("Intensity of seed image must be less than that "
                     "of the mask image for reconstruction by dilation.")
  elif method == 'erosion' and np.any(seed < mask):
    raise ValueError("Intensity of seed image must be greater than that "
                     "of the mask image for reconstruction by erosion.")
  try:
      from skimage.morphology._greyreconstruct import reconstruction_loop
  except ImportError:
      raise ImportError("_greyreconstruct extension not available.")

  if selem is None:
      selem = np.ones([3] * seed.ndim, dtype=bool)
  else:
      selem = selem.copy()

  if offset is None:
    if not all([d % 2 == 1 for d in selem.shape]):
        ValueError("Footprint dimensions must all be odd")
    offset = np.array([d // 2 for d in selem.shape])
  
  # Cross out the center of the structural_element
  selem[tuple(slice(d, d + 1) for d in offset)] = False

  # Make padding for edges of reconstructed image so we can ignore boundaries
  padding = (np.array(selem.shape) / 2).astype(int)
  dims = np.zeros(seed.ndim + 1, dtype=int)
  dims[1:] = np.array(seed.shape) + 2 * padding
  dims[0] = 2
  inside_slices = tuple(slice(p, -p) for p in padding)
  # Set padded region to minimum image intensity and mask along first axis so
  # we can interleave image and mask pixels when sorting.
  if method == 'dilation':
      pad_value = np.min(seed)
  elif method == 'erosion':
      pad_value = np.max(seed)
  images = np.ones(dims, dtype = seed.dtype) * pad_value
  images[(0,) + inside_slices] = seed
  images[(1,) + inside_slices] = mask

  # Create a list of strides across the array to get the neighbors within
  # a flattened array
  value_stride = np.array(images.strides[1:]) / images.dtype.itemsize
  image_stride = images.strides[0] // images.dtype.itemsize
  selem_mgrid = np.mgrid[[slice(-o, d - o)
                          for d, o in zip(selem.shape, offset)]]
  selem_offsets = selem_mgrid[:, selem].transpose()
  nb_strides = np.array([np.sum(value_stride * selem_offset)
                         for selem_offset in selem_offsets], np.int32)

  images = images.flatten()

  # Erosion goes smallest to largest; dilation goes largest to smallest.
  index_sorted = np.argsort(images).astype(np.int32)
  if method == 'dilation':
      index_sorted = index_sorted[::-1]

  # Make a linked list of pixels sorted by value. -1 is the list terminator.
  prev = -np.ones(len(images), np.int32)
  next = -np.ones(len(images), np.int32)
  prev[index_sorted[1:]] = index_sorted[:-1]
  next[index_sorted[:-1]] = index_sorted[1:]

  # Cython inner-loop compares the rank of pixel values.
  if method == 'dilation':
      value_rank, value_map = rank_order(images)
  elif method == 'erosion':
      value_rank, value_map = rank_order(-images)
      value_map = -value_map

  start = index_sorted[0]
  reconstruction_loop(value_rank, prev, next, nb_strides, start, image_stride)

  # Reshape reconstructed image to original image shape and remove padding.
  rec_img = value_map[value_rank[:image_stride]]
  rec_img.shape = np.array(seed.shape) + 2 * padding
  
  return rec_img[inside_slices]




[docs]
def grey_reconstruct(source, mask = None, sink = None, method = None, shape = 3, verbose = False):
  """Calculates the grey reconstruction of the image 
  
  Arguments
  ---------
  
  source : array
    The source image data.
  method : 'dilation' or 'erosion' or None
    The mehtjod to use, if None return original image.
  shape : in or tuple
    Shape of the strucuturing element for the grey reconstruction.
  verbose : boo;
    If True, print progress info.

  Returns
  -------
  reconstructed: array
    Grey reconstructed image.
    
  Note
  ----
  The reconstruction is done slice by slice along the z-axis.
  """
  
  if verbose:
    timer = tmr.Timer();
    hdict.pprint(head='Grey reconstruction', method=method, shape=shape)
  
  if method is None:
    return source;
  
  if sink is None:
    sink = np.empty(source.shape, dtype=source.dtype);
  
  # background subtraction in each slice
  selem = se.structure_element(form='Disk', shape=shape, ndim=2).astype('uint8');
  for z in range(source.shape[2]):
    #img[:,:,z] = img[:,:,z] - grey_opening(img[:,:,z], structure = structureElement('Disk', (30,30)));
    #img[:,:,z] = img[:,:,z] - morph.grey_opening(img[:,:,z], structure = self.structureELement('Disk', (150,150)));
    sink[:,:,z] = source[:,:,z] - reconstruct(source[:,:,z], mask=mask[:,:,z], method=method, selem=selem)
  
  if verbose:
    timer.print_elapsed_time('Grey reconstruction');
  
  return sink 




###############################################################################
### Test
###############################################################################

def _test():
  import numpy as np
  import ClearMap.ImageProcessing.GreyReconstruction as gr
  import ClearMap.Visualization.Plot3d as p3d

  x = np.random.rand(*(200,200,10));
  r = gr.grey_reconstruct(x, mask=0.5 * x, method='erosion', shape=3);
  
  p3d.plot([x,r])