"""
This module contains some operators used to filter the various connected components
of a set according to their number of holes.
"""
# Programs realised by Serge BEUCHER, (C) 2014
# Importing mamba and mambaComposed.
from mamba import *
from mambaComposed import *
def holesSieving(imIn, imOut, n):
"""
Puts in imOut all the connected components of imIn which contains exactly
n holes.
"""
imWrk0 = imageMb(imIn, 32)
imWrk1 = imageMb(imIn)
imOut.reset()
# Initial image labelling
nParticles = label(imIn, imWrk0)
vLabel = computeRange(imWrk0)[1]
for i in range(vLabel):
# each particle is extracted and its connectivity number
# is measured
if ((i+1) % 256)<>0:
threshold(imWrk0, imWrk1, i+1, i+1)
nc = computeConnectivityNumber(imWrk1)
# if the number of holes is equal to n, the particle is
# added to the output image
if (n==(1 - nc)):
logic(imWrk1, imOut, imOut, "sup")
def homotopyTreeBuild(imIn, imOut):
"""
Builds the homotopy tree of the initial binary set imIn. The result is stored
in image imOut. Each grey level corresponds to a hierarchy level of the homotopy
tree. Extracting each level i of embedding of particles and holes of the initial
set consists simply in a [i, i] thresholding of imOut.
"""
imWrk0 = imageMb(imIn)
imWrk1 = imageMb(imIn)
imWrk2 = imageMb(imIn, 8)
copy(imIn, imWrk0)
v = computeVolume(imIn)
imOut.reset()
i = 0
# Loop performed until no connected component remains.
while v!=0:
i = i+1
# Background extraction.
closeHoles(imWrk0, imWrk1)
negate(imWrk1, imWrk1)
# Connected components adjacent to the background are rebuilt.
dilate(imWrk1, imWrk1)
build(imWrk0, imWrk1)
# These particles are at level i in the homotopy tree.
# They are removed from the current set, giving access to the next level (in imWrk0).
diff(imWrk0, imWrk1, imWrk0)
v = computeVolume(imWrk0)
# The level-i particles are labelled with i and added to the label image.
convertByMask(imWrk1, imWrk2, 0, i)
add(imOut, imWrk2, imOut)
def closeOneHole(imIn, imOut):
"""
This procedure allows to close one and only one hole in each connected component
of the binary image imIn. When a connected component has no hole, it remains unchanged
in the final image stored in imOut.
This algorithm requires a single level of homotopy in the initial image.
"""
imWrk0 = imageMb(imIn, 32)
imWrk1 = imageMb(imIn, 32)
imWrk2 = imageMb(imIn, 32)
imWrk3 = imageMb(imIn)
imWrk4 = imageMb(imIn)
# The holes are filled in imWrk3 and extracted in imWrk4.
closeHoles(imIn, imWrk3)
diff(imWrk3, imIn, imWrk4)
# The holes are labelled.
nb = label(imWrk4, imWrk0)
# The image with filled holes is converted in 32-bit.
convertByMask(imWrk3, imWrk1, 0, computeMaxRange(imWrk1)[1])
# Geodesic reconstruction of the label image.
copy(imWrk0, imWrk2)
build(imWrk1, imWrk2)
# The holes with same label as the built image are extracted...
generateSupMask(imWrk0, imWrk2, imWrk4, False)
logic(imWrk4, imWrk3, imWrk3, "inf")
# ... and filled.
logic(imIn, imWrk3, imOut, "sup")
def holesFiltering(imIn, imOut, n):
"""
Removes from image imIn all the particles which contain up to n holes and
puts the result in imOut. imOut contains particles with more than n holes.
The homotopy tree of imIn is supposed to be simple (one level only).
n must be positive or equal to zero..
"""
imWrk1 = imageMb(imIn)
imWrk2 = imageMb(imIn)
copy(imIn, imWrk1)
i = n
while i>0:
closeOneHole(imWrk1, imWrk1)
i -= 1
closeHoles(imWrk1, imWrk2)
diff(imWrk2, imWrk1, imWrk2)
dilate(imWrk2, imWrk2)
build(imIn, imWrk2)
copy(imWrk2, imOut)
# The measure labelling procedure is defined.
def measureLabelling(imIn, imMeasure, imOut):
"""
Labelling each particle of the binary image 'imIn' with the number of pixels
in image 'imMeasure' contained in each particle. The result is put is the 32-bit
image 'imOut'.
"""
# Working images.
imWrk1 = imageMb(imIn, 32)
imWrk2 = imageMb(imIn)
imWrk3 = imageMb(imIn, 8)
imWrk4 = imageMb(imIn, 8)
imWrk5 = imageMb(imIn, 8)
imWrk6 = imageMb(imIn, 32)
# Output image is emptied.
imOut.reset()
# Labelling the initial image.
nbParticles = label(imIn, imWrk1)
# Defining output LUTs.
outLuts = [[0 for i in range(256)] for i in range(4)]
# Convering the imMeasure image to 8-bit.
convert(imMeasure, imWrk4)
while nbParticles > 0:
# particles with labels between 1 and 255 are extracted.
threshold(imWrk1, imWrk2, 0, 255)
convert(imWrk2, imWrk3)
copyBytePlane(imWrk1, 0, imWrk5)
logic(imWrk3, imWrk5, imWrk3, "inf")
# The points contained in each particle are labelled.
logic(imWrk3, imWrk4, imWrk5, "inf")
# The histogram is computed.
histo = getHistogram(imWrk5)
# The same operation is performed for the 255 particles.
for i in range(1, 256):
# The number of points in each particle is obtained from the histogram.
value = histo[i]
j = 3
# This value is splitted in powers of 256 and stored in the four
# output LUTs.
while j >= 0:
n = 2 ** (8 * j)
outLuts[j][i] = value / n
value = value % n
j -= 1
# each LUT is used to label each byte plane of a temporary image with the
# corresponding value.
for i in range(4):
lookup(imWrk3, imWrk5, outLuts[i])
copyBytePlane(imWrk5, i, imWrk6)
# The intermediary result is accumulated in the final image.
logic(imOut, imWrk6, imOut, "sup")
# 256 (in order to skip multiple of 256 values) is subtracted from the
# initial labelled image in order to process the next 255 particles.
floorSubConst(imWrk1, 256, imWrk1)
nbParticles -= 255
# Labelling each particle with its number of holes (up to 1).
def holesLabelling(imIn, imOut):
"""
Labels each particle in imIn with a value equal to its number of holes +1.
The result is put in 32-bit image imOut.
"""
# Working images.
imWrk1 = imageMb(imIn)
imWrk2 = imageMb(imIn, 32)
# Initializing the label image with 2.
convertByMask(imIn, imOut, 0, 2)
# Determining the 2nd configurations in the connectivity number calculation.
# and adding their number to the label image.
hitOrMiss(imIn, imWrk1, 2, 5)
measureLabelling(imIn, imWrk1, imWrk2)
add(imOut, imWrk2, imOut)
# Determining the 1st configurations and subtracting them to get the
# number of holes + 1.
hitOrMiss(imIn, imWrk1, 66, 1)
measureLabelling(imIn, imWrk1, imWrk2)
sub(imOut, imWrk2, imOut)