Implemented openCV segmentation

src/cli.py

@@ -32,12 +32,14 @@ def process(directory, check_segmentation, save_overlay, save_plot, verbose):
 
 @cli.command()
 @click.argument('directory', type=click.Path(exists=True))
+@click.option('-c', '--compare', is_flag=True, help="Compare scikit-image and openCV segmentation implementation")
 @click.option('-o', '--save-overlay', is_flag=True, help="Save segmented overlay to disk")
 @click.option('-v', '--verbose', count=True, help="Increase verbosity level")
-def segment(directory, save_overlay, verbose):
+def segment(directory, compare, save_overlay, verbose):
     '''Segment an image or directory of images and saves extracted droplets to disk'''
 
     from .data.segment_droplets import segment_droplets_to_file
+    from .data.compare_segmentation import segmentation_compare
 
     # Setup logging
     if verbose == 1:
@@ -46,6 +48,10 @@ def segment(directory, save_overlay, verbose):
         logging.basicConfig(level=logging.DEBUG, format='%(levelname)s - %(message)s')
     else:
         logging.basicConfig(level=logging.WARNING, format='%(levelname)s - %(message)s')
+
+    if compare:
+        segmentation_compare(directory)
+        return
 
     logging.info("Extracting droplets from: %s", directory)
     segment_droplets_to_file(directory, save_overlay=save_overlay)

src/crystal_processing/process_image_folder.py

@@ -43,7 +43,7 @@ def process_image(image_path, crop_box, model, save_overlay=False):
     cropped = crop(image, crop_box)
 
     # Segment image
-    (labeled, _, _) = segment(cropped)
+    (labeled, _) = segment(cropped)
 
     # Extract individual droplets
     drop_images, regProps = extract_indiv_droplets(cropped, labeled)
@@ -130,10 +130,10 @@ def process_image_folder(directory, crop_box=None, show_plot=False, save_overlay
         idx_80 = int(len(image_list) * 0.8)
         image_80 = crop(open_grey_scale_image(image_list[idx_80]), crop_box)
         logging.info("Segmentation check requested. Segmenting image %s", image_list[idx_80])
-        labeled, _, _ = segment(image_80)
+        labeled, _ = segment(image_80)
 
         from skimage.color import label2rgb
-        overlay_image = label2rgb(labeled, image=image_80, bg_label=0)
+        overlay_image = label2rgb(labeled, image=image_80, bg_label=1)
 
         import matplotlib
         matplotlib.use('Qt4Agg', force=True)
@@ -156,7 +156,6 @@ def process_image_folder(directory, crop_box=None, show_plot=False, save_overlay
                 result = input("Please press 'y' for yes or 'n' for no\n")
 
         plt.close()
-        plt.pause(0.1)
 
     # Compute the number of batches necessary
     num_images = len(image_list)

src/data/compare_segmentation.py

@@ -0,0 +1,88 @@
+'''
+Compare different segmentation implementations.
+'''
+
+import os
+import logging
+from time import time
+
+from skimage.color import label2rgb
+from skimage.exposure import equalize_adapthist
+import matplotlib
+matplotlib.use('Qt4Agg', force=True)
+import matplotlib.pyplot as plt
+
+from src.data.utils import open_grey_scale_image, select_rectangle, crop
+from src.data.segment_droplets import segment, segment_skimage
+
+def segmentation_compare(image_path, crop_box=None):
+    '''
+    Compares droplet segmentation algorithm implemented with scikit-image and CV2
+
+    Parameters
+    ----------
+    img: numpy.ndarray
+        Array representing the greyscale values (0-255) of an image cropped to show only the droplets region
+    '''
+
+    if not os.path.isfile(image_path):
+        raise ValueError("Must provide path to a single image.")
+
+    img = open_grey_scale_image(image_path)
+
+    # Get the crop box from the first image if not provided
+    logging.info('Getting crop box from image {}'.format(os.path.basename(image_path)))
+    if not crop_box:
+        crop_box = select_rectangle(img)
+
+    # Crop image
+    cropped = crop(img, crop_box)
+
+    # Segment image scikit-image
+    logging.info("Starting scikit-image segmentation")
+    t_0_ski = time()
+    (labeled_ski, num_reg_ski) = segment_skimage(cropped)
+    time_ski = time() - t_0_ski
+    logging.info("Finished scikit-image segmentation")
+
+    # Segment image openCV
+    logging.info("Starting OpenCV segmentation")
+    t_0_cv = time()
+    (labeled_cv, num_reg_cv) = segment(cropped)
+    time_cv = time() - t_0_cv
+    logging.info("Finished OpenCV segmentation")
+
+    # Print report
+    print(f"Scikit-image - Droplets: {num_reg_ski} - Time: {time_ski}s")
+    print(f"OpenCV       - Droplets: {num_reg_cv} - Time: {time_cv}s")
+
+    # Show segmented images
+
+    img_contrast = equalize_adapthist(cropped, clip_limit=0.06)
+
+    plt.ion()
+    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 15))
+    ax = axes.ravel()
+
+    ax[0].set_title('Skimage - Contour')
+    ax[0].imshow(img_contrast, cmap=plt.cm.gray, interpolation='nearest')
+    ax[0].contour(labeled_ski, [0.5], linewidths=1, colors='r')
+    ax[1].set_title('Skimage - Overlay')
+    ax[1].imshow(label2rgb(labeled_ski, image=img_contrast, bg_label=0), cmap=plt.cm.gray, interpolation='nearest')
+    ax[2].set_title('OpenCV - Contour')
+    ax[2].imshow(img_contrast, cmap=plt.cm.gray, interpolation='nearest')
+    ax[2].contour(labeled_cv, [0.5], linewidths=1, colors='r')
+    ax[3].set_title('OpenCV - Overlay')
+    ax[3].imshow(label2rgb(labeled_cv, image=img_contrast, bg_label=1), cmap=plt.cm.gray, interpolation='nearest')
+
+    for a in ax:
+        a.set_axis_off()
+
+    plt.tight_layout()
+    plt.show(block=False)
+
+    input("Press any key to exit...")
+
+    plt.close()
+
+    return

src/data/segment_droplets.py

@@ -3,10 +3,10 @@
 '''
 
 import os
+import warnings
 
 import numpy as np
 
-
 from skimage import io, exposure
 from skimage.color import label2rgb
 from skimage.exposure import equalize_adapthist
@@ -18,14 +18,15 @@
 import cv2
 from tqdm import tqdm
 
+
 from scipy import ndimage as ndi
 
-from src.data.utils import select_rectangle, open_grey_scale_image, crop
+from src.data.utils import select_rectangle, open_grey_scale_image, crop, clear_border
 
 
-def segment(img, exp_clip_limit=0.06, closing_disk_radius=4, rm_holes_area=8192, minima_minDist=100, mask_val=0.1):
+def segment_skimage(img, exp_clip_limit=0.06, closing_disk_radius=4, rm_holes_area=8192, minima_minDist=100, mask_val=0.1):
     '''
-    Segments droplets in an image using a watershed algorithm.
+    Segments droplets in an image using a watershed algorithm. Scikit-image implementation.
 
     Parameters
     ----------
@@ -44,14 +45,13 @@ def segment(img, exp_clip_limit=0.06, closing_disk_radius=4, rm_holes_area=8192,
 
     Returns
     -------
-    (labeled: numpy.ndarray, num_maxima: int, num_regions: int)
+    (labeled: numpy.ndarray, num_regions: int)
         labeled: labeled array of the same shape as input image where each region is assigned a disctinct integer label.
-        num_maxima: Number of maxima detected from the distance transform
         num_regions: number of labeled regions
     '''
 
     # Adaptive equalization
-    img_adapteq = equalize_adapthist(img, clip_limit = exp_clip_limit)
+    img_adapteq = equalize_adapthist(img, clip_limit=exp_clip_limit)
 
     # Minimum threshold
     threshold = threshold_otsu(img_adapteq)
@@ -64,7 +64,6 @@ def segment(img, exp_clip_limit=0.06, closing_disk_radius=4, rm_holes_area=8192,
 
     # Calculate the distance to the dark background
     distance = ndi.distance_transform_edt(rm_holes_closed)
-    #distance = cv2.distanceTransform(rm_holes_closed.astype('uint8'),cv2.DIST_L2,3) # TODO: test cv2 implementation for speed and acuraccy
 
     # Increase contrast of the the distance image
     cont_stretch = exposure.rescale_intensity(distance, in_range='image')
@@ -85,7 +84,73 @@ def segment(img, exp_clip_limit=0.06, closing_disk_radius=4, rm_holes_area=8192,
     # Label the segments of the image
     labeled, num_regions = ndi.label(segmented)
 
-    return (labeled, num_maxima, num_regions)
+    return (labeled, num_regions)
+
+def segment(img, exp_clip_limit=15):
+    '''
+    Segments droplets in an image using a watershed algorithm. OpenCV implementation.
+
+    Parameters
+    ----------
+    img: numpy.ndarray
+        Array representing the greyscale values (0-255) of an image cropped to show only the droplets region
+    exp_clip_limit: float [0-1], optional
+        clip_limit parameter for adaptive equalisation
+    closing_disk_radius: int, optional
+        diamater of selection disk for the closing function
+    rm_holes_area: int, optional
+        maximum area of holes to remove
+    minima_minDist: int, optional
+        minimum distance between peaks in local minima determination
+    mask_val: float, optional
+        Masking value (0-1) for the distance plot to remove small regions. Default 0.2
+
+    Returns
+    -------
+    (labeled: numpy.ndarray, num_regions: int)
+        labeled: labeled array of the same shape as input image where each region is assigned a disctinct integer label.
+        num_regions: number of labeled regions
+    '''
+
+    # Adaptive Equalization
+    clahe = cv2.createCLAHE(clipLimit=exp_clip_limit, tileGridSize=(8,8))
+    img_adapteq = clahe.apply(img)
+
+    # Thresholding (OTSU)
+    blur = cv2.GaussianBlur(img_adapteq, (5,5), 0)
+    _, binary = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
+
+    # Remove small dark regions
+    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(4,4))
+    closed = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations = 2)
+    fill_holes = ndi.morphology.binary_fill_holes(closed, structure=np.ones((3, 3))).astype('uint8')
+
+    # Sure background area
+    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3, 3))
+    sure_bg = np.uint8(cv2.dilate(fill_holes, kernel, iterations=1))
+
+    # Sure foreground area
+    dist_transform_fg = cv2.distanceTransform(fill_holes, cv2.DIST_L2, 5)
+    _, sure_fg = cv2.threshold(dist_transform_fg, 0.25*dist_transform_fg.max(), 255, 0)
+    clear_border(sure_fg)
+    sure_fg = np.uint8(sure_fg)
+
+    # Unknown region
+    unknown = cv2.subtract(sure_bg, sure_fg)
+
+    # Marker labelling
+    _, markers = cv2.connectedComponents(sure_fg)
+    # Add one to all labels so that sure background is not 0, but 1
+    markers = markers+1
+    # Now, mark the region of unknown with zero
+    markers[unknown > 0] = 0
+
+    # Run the watershed algorithm
+    three_channels = cv2.cvtColor(fill_holes, cv2.COLOR_GRAY2BGR)
+    segmented = cv2.watershed(three_channels.astype('uint8'), markers)
+
+    return (segmented, segmented.max()-1)
+
 
 def extract_indiv_droplets(img, labeled, border = 25, ecc_cutoff = 0.8):
     '''
@@ -155,13 +220,15 @@ def segment_droplets_to_file(image_filename, crop_box=None, save_overlay=False):
         cropped = crop(image, crop_box)
 
         # Segment image
-        (labeled, num_maxima, num_regions) = segment(cropped)
+        (labeled, num_regions) = segment(cropped)
 
         # Save the overlay image if requested
         if save_overlay:
             image_overlay = label2rgb(labeled, image=cropped, bg_label=0)
             filename = image_file.split('.')[0] + '_segmented.jpg'
-            io.imsave(filename, image_overlay)
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                io.imsave(filename, image_overlay)
 
         # Extract individual droplets
         drop_images, _ = extract_indiv_droplets(cropped, labeled)
@@ -175,4 +242,6 @@ def segment_droplets_to_file(image_filename, crop_box=None, save_overlay=False):
         # Save all the images in the output directory
         for (i, img) in enumerate(drop_images):
             name = out_directory + image_file.split('.')[0].split('/')[-1] + '_drop_' + str(i) + '.jpg'
-            io.imsave(name, img, check_contrast=False)
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                io.imsave(name, img, check_contrast=False)

src/data/utils.py

@@ -11,6 +11,7 @@
 from matplotlib.widgets import RectangleSelector
 from matplotlib import pyplot as plt
 import matplotlib as mpl
+import numpy as np
 import logging
 
 
@@ -93,4 +94,20 @@ def open_grey_scale_image(path):
 def crop(img, crop_box):
     '''Returns a cropped image for crop_box = (minRow, maxRow, minCol, maxCol)'''
     (minRow, minCol, maxRow, maxCol) = crop_box
-    return img[minRow:maxRow, minCol:maxCol]
+    return img[minRow:maxRow, minCol:maxCol]
+
+def clear_border(image):
+    '''Removes connected (white) items from the border of an image.'''
+    h, w = image.shape
+    mask = np.zeros((h + 2, w + 2), np.uint8)
+    for i in range(h-1): # Iterate on the lines
+        if image[i, 0] == 255:
+            cv2.floodFill(image, mask, (0, i), 0)
+        if image[i, w-1] == 255:
+            cv2.floodFill(image, mask, (w-1, i), 0)
+    for i in range(w-1): # Iterate on the columns
+        if image[0, i] == 255:
+            cv2.floodFill(image, mask, (i, 0), 0)
+        if image[h-1, i] == 255:
+            cv2.floodFill(image, mask, (i, h-1), 0)
+    return image