Merge branch 'convexFillHoles' into 'master'

Convex fill holes

See merge request !24

tools/label/label.py

@@ -1858,4 +1858,90 @@ def erodeLabels(lab, erosion=1, boundingBoxes=None, centresOfMass=None, numberOf
                                        labelDilate=-erosion, 
                                        numberOfThreads=numberOfThreads)
     return erodeImage
+
+def convexFillHoles(lab, boundingBoxes=None, centresOfMass=None):
+    """
+    This function fills the holes computing the convex volume around each label.
+
+    Parameters
+    -----------
+        lab : 3D numpy array
+            Labelled image
+            
+        boundingBoxes : lab.max()x6 array of ints, optional
+            Bounding boxes in format returned by ``boundingBoxes``.
+            If not defined (Default = None), it is recomputed by running ``boundingBoxes``
+
+        centresOfMass : lab.max()x3 array of floats, optional
+            Centres of mass in format returned by ``centresOfMass``.
+            If not defined (Default = None), it is recomputed by running ``centresOfMass``
+
+    Returns
+    --------
+        labOut : 3D numpy array
+            New labelled image.
+            
+    Note
+    ----
+        The function works nicely for convex particles. For non-convex particles, it will alter the shape.
+
+    """
+
+    # Check for boundingBoxes
+    if boundingBoxes is None:
+        boundingBoxes = spam.label.boundingBoxes(lab)
+    # Check for centresOfMass
+    if centresOfMass is None:
+        centresOfMass = spam.label.centresOfMass(lab)
+    # Create output label image
+    labOut = numpy.zeros_like(lab, dtype=spam.label.labelType)
+    # Get number of labels
+    numberOfLabels = lab.max()
+    # Create progressbar
+    widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
+    pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfLabels)
+    pbar.start()
+    for i in range(1,numberOfLabels+1,1):
+        # Get label
+        getLabelReturn = spam.label.getLabel(lab,
+                                             i,
+                                             labelDilate=0,
+                                             margin=3,
+                                             boundingBoxes=boundingBoxes,
+                                             centresOfMass=centresOfMass,
+                                             maskOtherLabels=False)
+        # Get subvolume
+        subVol = getLabelReturn['subvol']
+        # Transform to binary
+        subVolBinMask = (subVol > 0).astype(int)
+        # Mask out all the other labels
+        subVolBinMaskLabel = numpy.where(subVol == i, 1, 0).astype(int)
+        # Mask only the current label - save all the other labels
+        subVolMaskOtherLabel = subVolBinMask - subVolBinMaskLabel
+        # Fill holes with convex volume
+        points = numpy.transpose(numpy.where(subVolBinMaskLabel))
+        hull = scipy.spatial.ConvexHull(points)
+        deln = scipy.spatial.Delaunay(points[hull.vertices])
+        idx = numpy.stack(numpy.indices(subVol.shape), axis = -1)
+        out_idx = numpy.nonzero(deln.find_simplex(idx) + 1)
+        hullIm = numpy.zeros(subVol.shape)
+        hullIm[out_idx] = 1
+        hullIm = hullIm > 0
+        # Identify added voxels
+        subVolAdded = hullIm - subVolBinMaskLabel
+        # Identify the wrong voxels - they are inside other labels
+        subVolWrongAdded = subVolAdded * subVolMaskOtherLabel
+        # Remove wrong filling areas
+        subVolCorrect = (hullIm - subVolWrongAdded) > 0
+        # Get slice
+        grainSlice = (slice(getLabelReturn['slice'][0].start, getLabelReturn['slice'][0].stop),
+                        slice(getLabelReturn['slice'][1].start, getLabelReturn['slice'][1].stop),
+                        slice(getLabelReturn['slice'][2].start, getLabelReturn['slice'][2].stop))
+        # Add it to the output file
+        labOut[grainSlice][subVolCorrect] = i
+        # Update the progressbar
+        widgets[0] = progressbar.FormatLabel("{}/{} ".format(i, numberOfLabels))
+        pbar.update(i)
+        
+    return labOut
     

tools/tests/test_label.py

@@ -745,5 +745,42 @@ class TestFunctionLabel(unittest.TestCase):
         newCOM = spam.label.centresOfMass(imLab2)[-1]
         self.assertEqual(0, numpy.sum(iniCOM - newCOM))
         self.assertGreater(1.0, newVol/iniVol)
+        
+    def test_convexFillHoles(self):
+        
+        # Generate two spheres
+        grain1 = spam.label.label.Spheroid(15, 15, numpy.asarray([0,1,0])).digitize()
+        grain1 = numpy.pad(grain1, (5), 'constant', constant_values=(0))
+        grain1 = numpy.where(grain1, 1, 0)
+        grain2 = spam.label.label.Spheroid(15, 15, numpy.asarray([0,1,0])).digitize()
+        grain2 = numpy.pad(grain2, (5), 'constant', constant_values=(0))
+        grain2 = numpy.where(grain2, 2, 0)
+        grainIm = numpy.concatenate((grain1,grain2))
+        grainIm = numpy.zeros(grainIm.shape)
+        grainIm[:grain1.shape[0]-1,:,:] = grain1[:grain1.shape[0]-1,:,:]
+        grainIm[grain2.shape[0]-5:-5,:,:] =  grainIm[grain2.shape[0]-5:-5,:,:] + grain2[:,:,:]
+        grainIm = grainIm.astype(int)
+        # Compute initial volume & COM
+        COM = spam.label.centresOfMass(grainIm)
+        iniVol = spam.label.volumes(grainIm).astype(int)
+        # Dry run of the function
+        res1 = spam.label.convexFillHoles(grainIm)
+        # Compute new volume
+        res1Vol = spam.label.volumes(res1).astype(int)
+        # Test#1 -> Check that the volumes are the same
+        self.assertEqual(res1Vol[1], iniVol[1])
+        self.assertEqual(res1Vol[2], iniVol[2])
+        # Add a marker in the center of one of the labels
+        grainIm[COM[1][0].astype(int), COM[1][1].astype(int), COM[1][2].astype(int)] = 2
+        # Run again
+        res2 = spam.label.convexFillHoles(grainIm)
+        # Compute final volume
+        endVol = spam.label.volumes(res2).astype(int)
+        # Compare volumes
+        normVol1 = (endVol[1] - iniVol[1]) / iniVol[1]
+        normVol2 = (endVol[2] - iniVol[2]) / iniVol[2]
+        # Test#2 -> Check that Vol1 is the same and Vol2 increased
+        self.assertAlmostEqual(normVol1, 0, places = 3)
+        self.assertGreater(normVol2, 0)
 if __name__ == '__main__':
     unittest.main()