# Library of SPAM morphological functions.
# Copyright (C) 2020 SPAM Contributors
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>.
import multiprocessing
import numpy
import spam.helpers
try:
multiprocessing.set_start_method("fork")
except RuntimeError:
pass
import progressbar
import scipy.ndimage
import skimage.morphology
import spam.label
# operations on greyscale images
# Global number of processes
nProcessesDefault = multiprocessing.cpu_count()
[docs]
def greyDilation(im, nBytes=1):
"""
This function applies a dilation on a grey scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
nBytes: int, default=1
Number of bytes used to substitute the values on the border.
Returns
--------
numpy array
The dilated image
"""
# Step 1: check type and dimension
dim = len(im.shape)
# Step 2: Determine substitution value
sub = 2 ** (8 * nBytes) - 1
# Step 3: apply dilation # x y z
outputIm = im # +0 0 0
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 0, sub=sub)) # +1 0 0
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 0, sub=sub)) # -1 0 0
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 1, sub=sub)) # +0 1 0
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 1, sub=sub)) # +0 -1 0
if dim == 3:
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, 1, 2, sub=sub)) # 0 0 1
outputIm = numpy.maximum(outputIm, spam.helpers.singleShift(im, -1, 2, sub=sub)) # 0 0 -1
return outputIm
[docs]
def greyErosion(im, nBytes=1):
"""
This function applies a erosion on a grey scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
nBytes: int, default=1
Number of bytes used to substitute the values on the border.
Returns
--------
numpy array
The eroded image
"""
# Step 1: check type and dimension
dim = len(im.shape)
# Step 2: Determine substitution value
sub = 2 ** (8 * nBytes) - 1
# Step 1: apply erosion # x y z
outputIm = im # 0 0 0
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 0, sub=sub)) # 1 0 0
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 0, sub=sub)) # -1 0 0
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 1, sub=sub)) # 0 1 0
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 1, sub=sub)) # 0 -1 0
if dim == 3:
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, 1, 2, sub=sub)) # 0 0 1
outputIm = numpy.minimum(outputIm, spam.helpers.singleShift(im, -1, 2, sub=sub)) # 0 0 -1
return outputIm
[docs]
def greyMorphologicalGradient(im, nBytes=1):
"""
This function applies a morphological gradient on a grey scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
nBytes: int, default=1
Number of bytes used to substitute the values on the border.
Returns
--------
numpy array
The morphologycal gradient of the image
"""
return greyDilation(im, nBytes=nBytes) - im
# operation on binary images
[docs]
def binaryDilation(im, sub=False):
"""
This function applies a dilation on a binary scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
sub: bool, default=False
Subtitute value.
Returns
--------
numpy array
The dilated image
"""
# Step 0: import as bool
im = im.astype(bool)
# Step 1: check type and dimension
dim = len(im.shape)
# Step 1: apply dilation # x y z
outputIm = im # 0 0 0
outputIm = outputIm | spam.helpers.singleShift(im, 1, 0, sub=sub) # 1 0 0
outputIm = outputIm | spam.helpers.singleShift(im, -1, 0, sub=sub) # -1 0 0
outputIm = outputIm | spam.helpers.singleShift(im, 1, 1, sub=sub) # 0 1 0
outputIm = outputIm | spam.helpers.singleShift(im, -1, 1, sub=sub) # 0 -1 0
if dim == 3:
outputIm = outputIm | spam.helpers.singleShift(im, 1, 2, sub=sub) # 0 0 1
outputIm = outputIm | spam.helpers.singleShift(im, -1, 2, sub=sub) # 0 0 -1
return numpy.array(outputIm).astype("<u1")
[docs]
def binaryErosion(im, sub=False):
"""
This function applies a erosion on a binary scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
sub: bool, default=False
Substitute value.
Returns
--------
numpy array
The eroded image
"""
# Step 1: apply erosion with dilation --> erosion = ! dilation( ! image )
return numpy.logical_not(binaryDilation(numpy.logical_not(im), sub=sub)).astype("<u1")
[docs]
def binaryMorphologicalGradient(im, sub=False):
"""
This function applies a morphological gradient on a binary scale image
Parameters
-----------
im: numpy array
The input image (greyscale)
nBytes: int, default=False
Number of bytes used to substitute the values on the border.
Returns
--------
numpy array
The morphologycal gradient of the image
"""
return (numpy.logical_xor(binaryDilation(im, sub=sub), im)).astype("<u1")
[docs]
def binaryGeodesicReconstruction(im, marker, dmax=None, verbose=False):
"""
Calculate the geodesic reconstruction of a binary image with a given marker
Parameters
-----------
im: numpy.array
The input binary image
marker: 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
| ``[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
Verbose mode
Returns
--------
numpy.array
The reconstructed image
"""
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)}")
if verbose:
print(f"binaryGeodesicReconstruction: marker -> set plan in direction {direction} at distance {distance}")
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
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(f"binaryGeodesicReconstruction: geodesic distance = {d} (sum = {r1.sum()} & {r2.sum()})")
return r1 # send the reconstructed image
[docs]
def directionalErosion(bwIm, vect, a, c, nProcesses=nProcessesDefault, 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
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
nProcesses : integer (optional, default = nProcessesDefault)
Number of processes for multiprocessing
Default = number of CPUs in the system
verbose : boolean, optional (Default = False)
True for printing the evolution of the process
False for not printing the evolution of process
Returns
--------
imEroded : 3D boolean array
Booean array with the result of the erosion
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) is False:
print("spam.contacts.directionalErosion: The directional vector must be a list")
return
numberOfJobs = len(vect)
imEroded = numpy.zeros(bwIm.shape)
# Function for directionalErosion
global funDirectionalErosion
def funDirectionalErosion(job):
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, foreground=1
)
imEroded_i = scipy.ndimage.binary_erosion(bwIm, structure=spheroid)
return imEroded_i
if verbose:
widgets = [progressbar.FormatLabel(""), " ", progressbar.Bar(), " ", progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfJobs)
pbar.start()
finishedNodes = 0
# Run multiprocessing
with multiprocessing.Pool(processes=nProcesses) as pool:
for returns in pool.imap_unordered(funDirectionalErosion, range(0, numberOfJobs)):
if verbose:
finishedNodes += 1
widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedNodes, numberOfJobs))
pbar.update(finishedNodes)
imEroded = imEroded + returns
pool.close()
pool.join()
if verbose:
pbar.finish()
return imEroded
[docs]
def morphologicalReconstruction(im, selem=skimage.morphology.ball(1)):
"""
This functions performs a morphological reconstruction (greyscale opening followed by greyscale closing).
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
Returns
--------
imReconstructed : 3D boolean array
Greyscale image after the reconstuction
"""
# Perform the opening
imOpen = scipy.ndimage.grey_opening(im, footprint=selem)
# Perform the closing
imReconstructed = (scipy.ndimage.grey_closing(imOpen, footprint=selem)).astype(numpy.float32)
return imReconstructed
