Merge to generic binary geodesic reconstruction

.coveragerc

@@ -11,6 +11,9 @@ exclude_lines =
     elif mpi:
     def create2Patche*
     def create8Patche*
+    raise InputError
+    if verbose:
+
 
 [run]
 omit =
@@ -20,6 +23,7 @@ omit =
      */plotting/*
      *optionsParser.py
      */visual/*
+     tools/tests/errors.py
 
 source =
      spam

.gitignore

@@ -8,6 +8,7 @@
 .ipynb_checkpoints/
 .~notebook.ipynb
 .coverage
+.coverage.*
 coverage
 
 build/

examples/filters/plot_geodesic_reconstruction.py

+# -*- coding: utf-8 -*-
+"""
+Geodesic reconstruction
+========================
+
+This example shows how to use the geodesic reconstruction
+"""
+
+# sphinx_gallery_thumbnail_number = 1
+import matplotlib.pyplot as plt
+import numpy
+import spam.excursions
+import spam.filters
+
+############################################
+# Generate a morphology
+############################################
+# Generation and excrusion of a realisation of a correlated random field
+covariance = {'type': 'stable', 'alpha': 2.0, 'correlation_length': 0.01}
+realisation = spam.excursions.simulateRandomField(nNodes=500, covariance=covariance, nRea=1, dim=2, vtkFile="morpho")
+excursion = realisation > 0.1
+
+############################################
+# Geodesic reconstruction
+############################################
+# Direct definition of the marker (cube in the center)
+marker = numpy.zeros(excursion.shape, dtype=bool)
+marker[230:270, 230:270] = excursion[230:270, 230:270]
+recA = spam.filters.binaryGeodesicReconstruction(excursion, marker)
+
+##############################################
+# Definition of the marker with a list of plans
+marker = numpy.zeros(excursion.shape, dtype=bool)
+marker = [0, 0]  # direction 0 distance 0 (top side)
+recB = spam.filters.binaryGeodesicReconstruction(excursion, marker)
+marker = [1, 0, 0, -1]  # direction 1 distance 0 (left side) and direction 0 distance max (bottom side)
+recC = spam.filters.binaryGeodesicReconstruction(excursion, marker)
+
+######################
+# Plots
+fig, axs = plt.subplots(2, 2)
+plt.setp(axs, xticks=[], yticks=[])
+axs[0, 0].imshow(excursion)
+axs[0, 1].imshow(recA)
+axs[1, 0].imshow(recB)
+axs[1, 1].imshow(recC)
+axs[0, 0].set_title("Morphology")
+axs[0, 1].set_title("Geodesic reconstruction from a square in the center")
+axs[1, 0].set_title("Geodesic reconstruction from the top side")
+axs[1, 1].set_title("Geodesic reconstruction from the left and bottom side")
+
+plt.show()
+
+
+############################################
+# Notes
+############################################
+# The geodesic reconstruction is implemented for 2D and 3D binary morphologies

tools/errors.py

+class Error(Exception):
+    """Base class for exceptions in this module."""
+    pass
+
+class InputError(Error):
+    """Exception raised for errors in the input.
+
+    Attributes:
+        variable -- variable name which caused the error
+        explanation -- explanation of the error
+    """
+
+    def __init__(self, variable, explanation=None):
+        self.variable = variable
+        self.explanation = explanation
+        self.message = f"input error on variable {self.variable}"
+        if self.explanation is not None:
+            self.message += f' ({self.explanation})'
+
+        super().__init__(self.message)

tools/filters/morphologicalOperations.py

@@ -181,55 +181,20 @@ def binaryMorphologicalGradient(im, sub=False):
     return (numpy.logical_xor(binaryDilation(im, sub=sub), im)).astype('<u1')
 
 
-# def _binary_reconstruction_from_edges(im, dmax=0):
-#     """
-#     Calculate the morphological reconstruction of an image (binary) with the edges of a CUBE as a marker!
-#     PARAMETERS:
-#     - im (3D numpy.array): The input image (binary)
-#     - dmax (int): the maximum geodesic distance. If zero, the reconstruction is complete.
-#     RETURNS:
-#     - (3D numpy.array): The reconstructed image
-#     TODO:
-#     - Consider mask for other geometry than cube
-#     #- Consider different marker than the edges (for example, two opposite faces for percolation test)
-#     HISTORY:
-#     2016-04-21 (Sun Yue): First version of the function
-#     """
-#     # Step 1: compute marker
-#     size = im.shape[0]
-#     bord = numpy.zeros((size, size, size), dtype=bool)
-#     bord[0, :, :] = im[0, :, :]
-#     bord[-1, :, :] = im[-1, :, :]
-#     bord[:, 0, :] = im[:, 0, :]
-#     bord[:, -1, :] = im[:, -1, :]
-#     bord[:, :, 0] = im[:, :, 0]
-#     bord[:, :, -1] = im[:, :, -1]
-#
-#     # Step 2: first dilation and intersection
-#     temp1 = (binaryDilation(bord)) & (im)
-#     temp2 = temp1
-#     temp1 = (binaryDilation(temp2)) & (im)
-#     distance = 1
-#
-#     if dmax == 0:
-#         dmax = 1e99  # perform complete reconstruction
-#     while ((not(numpy.array_equal(temp1, temp2))) & (distance < dmax)):
-#         temp2 = temp1
-#         temp1 = (binaryDilation(temp2)) & (im)
-#         distance += 1
-#         print('distance =', distance)
-#
-#     return temp1  # send the reconstructed image
-
-
-def binaryReconstructionFromEdges(im, dmax=None, verbose=False):
+def binaryGeodesicReconstruction(im, marker, dmax=None, verbose=False):
     """
-    Calculate the morphological reconstruction of an binary image with the edges of a cube as a marker
+    Calculate the geodesic reconstruction of a binary image with a given marker
 
     Parameters
     -----------
-        im: numpy array
-            The input image
+        im: 3D numpy array
+            The input binary image
+        marker: 3D numpy array or list
+            If numpy array: direct input of the marker (must be the size of im)
+            If list: description of the plans of the image considered as the marker
+                example: [1, 0] plan defined by all voxels at x1=0
+                         [0, -1] plan defined by all voxels at x0=x0_max
+                         [0, 0, 2, 5] plans defined by all voxels at x0=0 and x2=5
         dmax: int, default=None
             The maximum geodesic distance. If None, the reconstruction is complete.
         verbose: bool, default=False
@@ -240,95 +205,93 @@ def binaryReconstructionFromEdges(im, dmax=None, verbose=False):
         numpy.array
             The reconstructed image
     """
-    size = im.shape[0]
-    bord = numpy.zeros((size, size, size), dtype=bool)
-    bord[0, :, :] = im[0, :, :]
-    bord[-1, :, :] = im[-1, :, :]
-    bord[:, 0, :] = im[:, 0, :]
-    bord[:, -1, :] = im[:, -1, :]
-    bord[:, :, 0] = im[:, :, 0]
-    bord[:, :, -1] = im[:, :, -1]
-
-    # Step 2: first dilation and intersection
-    rec1 = (binaryDilation(bord)) & (im)
-    rec2 = rec1
-    rec1 = (binaryDilation(rec2)) & (im)
-    distance = 1
-
-    if dmax is None:
-        dmax = numpy.inf  # perform complete reconstruction
-    while ((not(numpy.array_equal(rec1, rec2))) & (distance < dmax)):
-        rec2 = rec1
-        rec1 = (binaryDilation(rec2)) & (im)
-        distance += 1
-        if verbose:
-            print('Distance = {}'.format(distance))
-
-    return rec1  # send the reconstructed image
-
-
-def reconstructionFromOppositeFaces(im, dmax=None, verbose=False):
-    """
-    Calculate the morphological reconstruction of an binary image with two opposite faces of a cube as a marker
+    from spam.errors import InputError
+
+    # Step 1: Define marker
+    if isinstance(marker, list):
+        # marker based on list of plans
+        if len(marker) % 2:
+            raise InputError("marker", explanation="len(marker) must be a multiple of 2")
+
+        plans = marker[:]
+        marker = numpy.zeros(im.shape, dtype=bool)
+        for i in range(len(plans) // 2):
+            direction = plans[2 * i]
+            distance = plans[2 * i + 1]
+            if len(im.shape) == 2:
+                if direction == 0:
+                    marker[distance, :] = im[distance, :]
+                elif direction == 1:
+                    marker[:, distance] = im[:, distance]
+                else:
+                    raise InputError("marker", explanation=f"Wrong marker plan direction {direction}")
+
+            elif len(im.shape) == 3:
+                if direction == 0:
+                    marker[distance, :, :] = im[distance, :, :]
+                elif direction == 1:
+                    marker[:, distance, :] = im[:, distance, :]
+                elif direction == 2:
+                    marker[:, :, distance] = im[:, :, distance]
+                else:
+                    raise InputError("marker", explanation=f"Wrong marker plan direction {direction}")
+
+            else:
+                raise InputError("marker", explanation=f"Image dimension should be 2 or 3, not {len(im.shape)}")
 
-    Parameters
-    -----------
-        im: numpy array
-            The input image
-        dmax: int, default=None
-            The maximum geodesic distance. If None, the reconstruction is complete.
-        verbose: bool, default=False
-            Verbose mode
+            if verbose:
+                print(f"binaryGeodesicReconstruction: marker -> set plan in direction {direction} at distance {distance}")
 
-    Returns
-    --------
-        numpy.array
-            The reconstructed image
-    """
-    # Step 1: compute marker
-    size = im.shape[0]
-    bord = numpy.zeros((size, size, size), dtype=bool)
-    bord[0, :, :] = im[0, :, :]
-    bord[-1, :, :] = im[-1, :, :]
+    elif isinstance(marker, numpy.ndarray):
+        # direct input of the marker
+        if im.shape != marker.shape:
+            raise InputError("marker", explanation="im and marker must have same shape")
+    else:
+        raise InputError("marker", explanation="must be a numpy array or a list")
 
     # Step 2: first dilation and intersection
-    rec1 = (binaryDilation(bord)) & (im)
-    rec2 = rec1
-    rec1 = (binaryDilation(rec2)) & (im)
-    distance = 1
-
-    if dmax is None:
-        dmax = numpy.inf  # perform complete reconstruction
-    while ((not(numpy.array_equal(rec1, rec2))) & (distance < dmax)):
-        rec2 = rec1
-        rec1 = (binaryDilation(rec2)) & (im)
-        distance += 1
+    r1 = binaryDilation(marker) & im
+    r2 = r1
+    r1 = binaryDilation(r2) & im
+    d = 1
+    dmax = numpy.inf if dmax is None else dmax
+    if verbose:
+        print(f'binaryGeodesicReconstruction: geodesic distance = {d} (sum = {r1.sum()} & {r2.sum()})')
+
+    # binary dilation until:
+    #   geodesic distance reach dmax
+    #   reconstuction complete
+    while not numpy.array_equal(r1, r2) and d < dmax:
+        r2 = r1
+        r1 = binaryDilation(r2) & im
+        d += 1
         if verbose:
-            print('Distance =', distance)
+            print(f'binaryGeodesicReconstruction: geodesic distance = {d} (sum = {r1.sum()} & {r2.sum()})')
+
+    return r1  # send the reconstructed image
 
-    return rec1  # send the reconstructed image
 
 def directionalErosion(bwIm, vect, a, c, numberOfThreads=1, verbose = False):
     """
     This functions performs direction erosion over the binarized image using
     an ellipsoidal structuring element over a range of directions. It is highly
-    recommended that the structuring element is slightly smaller than the 
+    recommended that the structuring element is slightly smaller than the
     expected particle (50% smaller in each axis is a fair guess)
-    
+
     Parameters
     -----------
         bwIm : 3D numpy array
             Binarized image to perform the erosion
-            
+
         vect : list of n elements, each element correspond to a 1X3 array of floats
             List of directional vectors for the structuring element
-        
+
         a : int or float
             Length of the secondary semi-axis of the structuring element in px
-            
+
         c : int or float
             Lenght of the principal semi-axis of the structuring element in px
-            
+
         numberOfThreads : integer, optional
             Number of Threads for multiprocessing.
             Default = 1
@@ -345,18 +308,18 @@ def directionalErosion(bwIm, vect, a, c, numberOfThreads=1, verbose = False):
     Note
     -----
         Taken from https://sbrisard.github.io/posts/20150930-orientation_correlations_among_rice_grains-06.html
-        
+
     """
-    
-    
+
+
     #Check if the directional vector input is a list
     if isinstance(vect,list) == False:
         print("spam.contacts.directionalErosion: The directional vector must be a list")
         return
-        
+
     numberOfJobs = len( vect )
     imEroded = numpy.zeros(bwIm.shape)
-    
+
     def worker( workerNumber, qJobs, qResults ):
         while True:
             job = qJobs.get()
@@ -364,35 +327,35 @@ def directionalErosion(bwIm, vect, a, c, numberOfThreads=1, verbose = False):
             if job == "STOP":
                 qResults.put("STOP")
                 break
-            
+
             maxDim = numpy.max([a,c])
             spheroid = spam.kalisphera.makeBlurryNoisySpheroid([maxDim,maxDim,maxDim],
                                                                [numpy.floor(maxDim/2), numpy.floor(maxDim/2), numpy.floor(maxDim/2)],
-                                                               [a,c], 
-                                                               vect[job], 
-                                                               background=0, 
+                                                               [a,c],
+                                                               vect[job],
+                                                               background=0,
                                                                foreground=1)
             imEroded_i = scipy.ndimage.binary_erosion(bwIm, structure = spheroid)
-            
+
             qResults.put( [imEroded_i])
-            
-            
+
+
     qJobs = multiprocessing.Queue()
     qResults = multiprocessing.Queue()
-    
+
     # print "Master: Adding jobs to queues"
     for x in range(numberOfJobs):
         # qJobs.put( contactList[x,0] )
         qJobs.put(x)
-    
+
     for i in range(numberOfThreads):
         qJobs.put("STOP")
-        
+
     # print "Master: Launching workers"
     for i in range(numberOfThreads):
         p = multiprocessing.Process(target=worker, args=(i, qJobs, qResults, ))
         p.start()
-        
+
     if verbose:
         pbar = progressbar.ProgressBar(maxval=numberOfJobs).start()
     finishedThreads = 0
@@ -407,7 +370,7 @@ def directionalErosion(bwIm, vect, a, c, numberOfThreads=1, verbose = False):
 
         else:
             # print "Master: got {}".format( result )
-            
+
             imEroded = imEroded + result[0]
             finishedJobs += 1
             if verbose:
@@ -423,13 +386,13 @@ def morphologicalReconstruction(im, selem = None):
     by dilation followed by a reconstruction by erosion. The ouput image presents less
     variability in the greyvalues inside each phase, without modifying the original
     shape of the objects of the image.
-    - 
-    
+    -
+
     Parameters
     -----------
         im : 3D numpy array
             Greyscale image to perform the reconstuction
-            
+
         selem : structuring element, optional
             Structuring element
             Default = None
@@ -437,8 +400,8 @@ def morphologicalReconstruction(im, selem = None):
     Returns
     --------
         imReconstructed : 3D boolean array
-            Greyscale image after the reconstuction   
-        
+            Greyscale image after the reconstuction
+
     """
     #Create the seed for the dilation
     imErosion = skimage.morphology.erosion(im,selem = selem)
@@ -448,5 +411,5 @@ def morphologicalReconstruction(im, selem = None):
     imDilation = skimage.morphology.dilation(imOpening, selem = selem)
     #Erode the seed until it reaches the mask
     imReconstructed = skimage.morphology.reconstruction(seed = imDilation, mask = imOpening, method='erosion')
-    
+
     return imReconstructed

tools/tests/test_filters.py

-# -*- coding: utf-8 -*-
+# -*- coding: utf-8 -*-
 from __future__ import print_function
 
 import unittest
@@ -60,131 +60,134 @@ class testAll(unittest.TestCase):
         self.assertEqual(sum, 0)
 
     def test_DirectionalErosion(self):
-        #Generate random ellipsoid grain
+        # Generate random ellipsoid grain
         theta = numpy.radians(random.randrange(30,300,1))
         phi = numpy.radians(random.randrange(0,90,1))
         vect = numpy.array([numpy.cos(phi), numpy.sin(phi)*numpy.sin(theta), numpy.sin(phi)*numpy.cos(theta)])
-        if vect[0]<0: vect[:]=-1*vect[:] 
+        if vect[0]<0: vect[:]=-1*vect[:]
         a = random.randrange(1,10,1)
         c = random.randrange(1,10,1)
-        print(a,c)
         if numpy.abs(a-c)<2:a = a*2
         maxDim = numpy.max([a+10,c+10])
         imTest = spam.kalisphera.makeBlurryNoisySpheroid([maxDim,maxDim,maxDim],
                                                            [numpy.floor(maxDim/2), numpy.floor(maxDim/2), numpy.floor(maxDim/2)],
-                                                           [a,c], 
-                                                           vect, 
-                                                           background=0, 
+                                                           [a,c],
+                                                           vect,
+                                                           background=0,
                                                            foreground=1)
-        #Check that a non-list entry of vect generate error
+        # Check that a non-list entry of vect generate error
         imEroded = morph.directionalErosion(imTest, vect, a, c)
         self.assertIs(imEroded, None)
-        #Perform watershed to label the grains
-        labIm = ws.watershed(imTest) 
-        #Perform directional erosion
+        # Perform watershed to label the grains
+        labIm = ws.watershed(imTest)
+        # Perform directional erosion
         imEroded = morph.directionalErosion(imTest, [vect], a, c)
         # Label the markers
         markers = spam.label.watershed(imEroded)
         # Compute COM
         COM = spam.label.centresOfMass(markers)
-        #Check that the center of mass of the marker lies inside the labelled grain        
+        # Check that the center of mass of the marker lies inside the labelled grain
         self.assertEqual(labIm[int(COM[1][0]),int(COM[1][1]),int(COM[1][2])], numpy.unique(labIm)[1])
-        
+
     def test_greyDilation(self):
-        #Generate single sphere
+        # Generate single sphere
         im1 = numpy.zeros((20, 20, 20), dtype="<f8")
         spam.kalisphera.makeSphere(im1, [10, 10, 10], 5)
         im1 = im1*255
-        #Run
+        # Run
         res = spam.filters.greyDilation(im1)
-        #Test that it work without problem
+        # Test that it work without problem
         self.assertIsNotNone(res)
-        #Test that the result make sense for the operation 
+        # Test that the result make sense for the operation
         self.assertGreater(numpy.sum(res), numpy.sum(im1))
-        
+
     def test_greyErosion(self):
-        #Generate single sphere
+        # Generate single sphere
         im1 = numpy.zeros((20, 20, 20), dtype="<f8")
         spam.kalisphera.makeSphere(im1, [10, 10, 10], 5)
         im1 = im1*255
-        #Run
+        # Run
         res = spam.filters.greyErosion(im1)
-        #Test that it work without problem
+        # Test that it work without problem
         self.assertIsNotNone(res)
-        #Test that the result make sense for the operation 
+        # Test that the result make sense for the operation
         self.assertGreater(numpy.sum(im1), numpy.sum(res))
-        
+
     def test_greyMorphologicalGradient(self):
-        #Generate single sphere
+        # Generate single sphere
         im1 = numpy.zeros((20, 20, 20), dtype="<f8")
         spam.kalisphera.makeSphere(im1, [10, 10, 10], 5)
         im1 = im1*255
-        #Run
+        # Run
         res = spam.filters.greyMorphologicalGradient(im1)
-        #Test that it work without problem
+        # Test that it work without problem
         self.assertIsNotNone(res)
-        
+
     def test_binaryDilation(self):
-        #Generate single sphere
+        # Generate single sphere
         im1 = numpy.zeros((20, 20, 20), dtype="<f8")
         spam.kalisphera.makeSphere(im1, [10, 10, 10], 5)
         im1 = im1*255
-        #Binarise
+        # Binarise
         im1Bin = im1 > 255/2
-        #Run
+        # Run
         res = spam.filters.binaryDilation(im1Bin)
-        #Test that it work without problem
+        # Test that it work without problem
         self.assertIsNotNone(res)
-        #Test that the result make sense for the operation 
+        # Test that the result make sense for the operation
         self.assertGreater(numpy.sum(res), numpy.sum(im1Bin))
-        
+
     def test_binaryErosion(self):
-        #Generate single sphere
+        # Generate single sphere
         im1 = numpy.zeros((20, 20, 20), dtype="<f8")
         spam.kalisphera.makeSphere(im1, [10, 10, 10], 5)
         im1 = im1*255
-        #Binarise
+        # Binarise
         im1Bin = im1 > 255/2
-        #Run
+        # Run
         res = spam.filters.binaryErosion(im1Bin)
-        #Test that it work without problem
+        # Test that it work without problem
         self.assertIsNotNone(res)
-        #Test that the result make sense for the operation 
+        # Test that the result make sense for the operation
         self.assertGreater(numpy.sum(im1Bin), numpy.sum(res))
-        
-    def test_binaryReconstructionFromEdges(self):
-        #Generate image
-        im = numpy.ones((20,20,20))
-        im[5:15, 5:15, 5:15] = 0
-        im[10,10,10] = 1
-        im = im.astype('uint8')
-        #Run
-        res = spam.filters.morphologicalOperations.binaryReconstructionFromEdges(im, verbose = True)
-        #Test that it work without problem
-        self.assertIsNotNone(res)
-        #Test that it makes sense
-        self.assertGreater(numpy.sum(im), numpy.sum(res))
-        
-    def test_reconstructionFromOppositeFaces(self):
-        #Generate image
-        im = numpy.ones((20,20,20))
-        im[5:15, 5:15, 5:15] = 0
-        im[10,10,10] = 1
-        im = im.astype('uint8')
-        #Run
-        res = spam.filters.morphologicalOperations.reconstructionFromOppositeFaces(im, verbose = True)