Partially implement dense and dilute phase mask generation and intensity value...

Partially implement dense and dilute phase mask generation and intensity value calculations. Not complete, not tested.

processMFLLPS.py

 import sys
 
-from skimage import io, filters, restoration, util, segmentation, feature, morphology, exposure, measure, draw, registration, transform
+from skimage import io, filters, restoration, util, segmentation, feature, morphology, exposure, measure, draw, \
+    registration, transform
 from scipy import ndimage as ndi
 from scipy import signal, misc
 import matplotlib.pyplot as plt
 import numpy as np
 import argparse
+import os
 
 
 def parse_args(args):
@@ -14,14 +16,17 @@ def parse_args(args):
     parser.add_argument("--input", help="Input raw image file(s).")
     parser.add_argument("--dark", help="Dark image for background subtraction.")
     parser.add_argument("--flatfield", help="Image for flatfield correction.")
-    parser.add_argument("--colorchannels", help="Which color channels to use in segmentation.",
+    parser.add_argument("--segchannels", help="Which color channels to use in droplet boundary segmentation.",
                         default=[0], nargs="+", type=int)
+    parser.add_argument("--densechannels", help="Which color channels to use for dense phase segmentation.",
+                        default=[1, 2], nargs="+", type=int)
     parser.add_argument("--writeintermediates", help="Write intermediate images during segmentation",
                         action="store_true")
     parser.add_argument("--zproject", "For 3D images, method to be used for z-projection.", default="median")
     parser.add_argument("--histequalsize", "Window size for local histogram equalization", default=64, type=int)
     parser.add_argument("--blur", "Sigma parameter for Gaussian blur during droplet segmentation.", default=2, type=int)
-    parser.add_argument("--rollingball", "Radius of rolling ball for background subtration during droplet segmentation.",
+    parser.add_argument("--rollingball",
+                        "Radius of rolling ball for background subtration during droplet segmentation.",
                         default=10, type=int)
     parser.add_argument("--sauvolapars", help="Window, k, and R parameters for Sauvola local thresholding.",
                         nargs=3, default=[15, 0.5, 128], type=int)
@@ -79,7 +84,7 @@ def segment_image_stack(image,
                         sauvola_pars=(15, 0.5, 128),
                         dist_trans_rel_thresh=0.7,
                         use_peak_local_max=False,
-                        peak_local_max_window=(20,20)):
+                        peak_local_max_window=(20, 20)):
     # Prepare a list for saving intermediate images
     intermediate_images = []
     # Handle multicolor images - get a single channel or sum across given channels
@@ -107,20 +112,20 @@ def segment_image_stack(image,
     if nz > 1:
         step += 1
         work_image = z_project(work_image, 0)
-        append_intermediate(intermediate_images, work_image, str(step)+"_"+"z_project", write_intermediates)
+        append_intermediate(intermediate_images, work_image, str(step) + "_" + "z_project", write_intermediates)
     # Rescale to (0, 1) and convert to 8-bit
     step += 1
     work_image = util.img_as_ubyte(exposure.rescale_intensity(work_image, out_range=(0., 1.)))
-    append_intermediate(intermediate_images, work_image, str(step)+"_"+"rescale", write_intermediates)
+    append_intermediate(intermediate_images, work_image, str(step) + "_" + "rescale", write_intermediates)
     # Local histogram equalization
     # https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_local_equalize.html#sphx-glr-auto-examples-color-exposure-plot-local-equalize-py
     step += 1
     work_image = filters.rank.equalize(work_image, footprint=morphology.disk(hist_equal_size))
-    append_intermediate(intermediate_images, work_image, str(step)+"_"+"hist_equal", write_intermediates)
+    append_intermediate(intermediate_images, work_image, str(step) + "_" + "hist_equal", write_intermediates)
     # Gaussian blur
     step += 1
     work_image = filters.gaussian(work_image, sigma=gaussian_blur_r, preserve_range=True)
-    append_intermediate(intermediate_images, work_image, str(step)+"_"+"gauss_blur", write_intermediates)
+    append_intermediate(intermediate_images, work_image, str(step) + "_" + "gauss_blur", write_intermediates)
     # Background subtraction
     step += 1
     work_image = work_image - restoration.rolling_ball(work_image, radius=rolling_ball_radius)
@@ -176,7 +181,51 @@ def segment_image_stack(image,
     labels = segmentation.watershed(-dist_transform, markers, watershed_line=True)
     append_intermediate(intermediate_images, labels, str(step) + "_" + "watershed", write_intermediates)
 
-    return labels, intermediate_images
+    return labels, intermediate_images, step
+
+
+def extract_droplets(intensity_img, label_img):
+    img_type = get_image_type(intensity_img)
+    nz, nr, nc, nchan = get_image_dims(intensity_img)
+    max_row = nr
+    max_col = nc
+    droplets = []
+    masks = []
+    label_index = []
+    label_img_no_edge = np.zeros(label_img.shape)
+    for region in measure.regionprops(label_img):
+        if region.bbox[0] > 0 and region.bbox[1] > 0 and region.bbox[2] < max_row and region.bbox[3] < max_col:
+            if img_type == 'grayscale':  # 2D image - no slices, grayscale
+                droplets.append(intensity_img[region.slice[0], region.slice[1]])
+            elif img_type == 'multicolor':  # 2D image - no slices, last dimension has color channels
+                droplets.append(intensity_img[region.slice[0], region.slice[1], :])
+            elif img_type == '3d_grayscale':  # 3D image - has z-slices, grayscale
+                droplets.append(intensity_img[:, region.slice[0], region.slice[1]])
+            elif img_type == '3d_multicolor':  # 3D color image - has z-slices, last dimension has color channels
+                droplets.append(intensity_img[:, region.slice[0], region.slice[1], :])
+            masks.append(region.image)
+            label_img_no_edge[tuple(region.coords.T)] = region.label
+            label_index.append(region.label)
+    return droplets, masks, label_index, label_img_no_edge
+
+
+def get_intensity_values(cropped_droplet_img, zplane, dense_mask, dil_mask, img_type, densechannels):
+    # Phase masks as list for iteration purposes
+    phase_masks = [dense_mask, dil_mask]
+    # Row 1: dense, Row 2: dilute, Columns: channels
+    mean_intensities = np.array((2, len(densechannels)))
+    sd_intensities = np.array((2, len(densechannels)))
+    for phase in range(0, 2):
+        for channel in range(0, len(densechannels)):
+            # Masked image
+            masked_image = cropped_droplet_img[zplane, :, :, densechannels[channel]][phase_masks[phase]]
+            mean_intensities[phase, channel] = np.mean(masked_image)
+            sd_intensities[phase, channel] = np.mean(masked_image)
+    dense_means = mean_intensities[0, :]
+    dil_means = mean_intensities[1, :]
+    dense_sds = sd_intensities[0, :]
+    dil_sds = sd_intensities[1, :]
+    return dense_means, dil_means, dense_sds, dil_sds
 
 
 def main(args):
@@ -184,21 +233,77 @@ def main(args):
     z_project_methods = {"median": np.median,
                          "max": np.max,
                          "mean": np.mean}
+    # Begin processing loop
+    # Iterate over droplet image stacks
     for img_file in args.input:
         img = io.imread(img_file)
-        droplet_mask, intermediate_images = segment_image_stack(img,
-                                                                channels=args.channels,
-                                                                write_intermediates=args.writeintermediates,
-                                                                z_project=z_project_methods[args.zproject],
-                                                                hist_equal_size=args.histequalsize,
-                                                                gaussian_blur_r=args.blur,
-                                                                rolling_ball_radius=args.rollingball,
-                                                                sauvola_pars=args.sauvolapars,
-                                                                dist_trans_rel_thresh=args.threshdisttrans,
-                                                                use_peak_local_max=args.peaklocalmax,
-                                                                peak_local_max_window=2*[args.windowpeaklocalmax])
+        img_type = get_image_type(img)
+        droplet_mask, intermediate_images, step = segment_image_stack(img,
+                                                                      channels=args.segchannels,
+                                                                      write_intermediates=args.writeintermediates,
+                                                                      z_project=z_project_methods[args.zproject],
+                                                                      hist_equal_size=args.histequalsize,
+                                                                      gaussian_blur_r=args.blur,
+                                                                      rolling_ball_radius=args.rollingball,
+                                                                      sauvola_pars=args.sauvolapars,
+                                                                      dist_trans_rel_thresh=args.threshdisttrans,
+                                                                      use_peak_local_max=args.peaklocalmax,
+                                                                      peak_local_max_window=2 * [
+                                                                          args.windowpeaklocalmax])
+        # Extract cropped droplets
+        # Create an image with relevant color channels summed
+        if img_type == 'grayscale' or img_type == '3d_grayscale':
+            # Don't need to do any summing - only one color
+            dense_phase_sum = img
+        elif img_type == 'multicolor':
+            dense_phase_sum = np.zeros((img.shape[0], img.shape[1]))
+            for channel in args.densechannels:
+                dense_phase_sum += img[:, :, channel]
+            # Rescale intensities
+            dense_phase_sum = util.img_as_uint(exposure.rescale_intensity(dense_phase_sum, out_range=(0., 1.)))
+        elif img_type == "3d_multicolor":
+            dense_phase_sum = np.zeros((img.shape[0], img.shape[1], img.shape[2]))
+            for channel in args.densechannels:
+                dense_phase_sum += img[:, :, :, channel]
+            # Rescale intensities
+            dense_phase_sum = util.img_as_uint(exposure.rescale_intensity(dense_phase_sum, out_range=(0., 1.)))
+        step += 1
+        droplets, masks, label_index, label_img_no_edge = extract_droplets(dense_phase_sum, droplet_mask)
+        append_intermediate(intermediate_images, label_img_no_edge, str(step) + "_" + "no_edge_droplets",
+                            args.write_intermediates)
+        # Save intermediates
+        if args.writeintermediates:
+            img_base_name = os.path.basename(img_file)
+            os.mkdir(img_base_name)
+            for img, name in intermediate_images:
+                img_name = img_base_name + "_" + name + ".tif"
+                io.imsave(img_name, img)
+        # Iterate over cropped droplet intensity images
+        dil_intensity_mean = []
+        dil_intensity_sd = []
+        dense_intensity_mean = []
+        dense_intensity_sd = []
+        for droplet, mask in zip(droplets, masks):
+            blob = feature.blob_log(droplet, 1/np.sqrt(3), 10/np.sqrt(3), overlap=0.1)
+            if blob.shape[0] == 1:
+                # Create the dense phase mask
+                rcoords, ccoords = draw.disk(center=tuple(blob.flat)[1:3], radius=blob.flat[-1], shape=mask.shape)
+                dense_mask = np.full(fill_value=False, shape=mask.shape, dtype='bool')
+                dense_mask[rcoords, ccoords] = True
+                # Create the dilute phase mask
+                rcoords, ccoords = draw.disk(center=tuple(blob.flat)[1:3], radius=4*blob.flat[-1], shape=mask.shape)
+                dilute_inner_mask = np.full(fill_value=False, shape=mask.shape, dtype='bool')
+                dilute_inner_mask[rcoords, ccoords] = True
+                dilute_mask = np.logical_xor(mask, dilute_inner_mask)
+                # Get intensity values
+                dense_means, dil_means, dense_sds, dil_sds = get_intensity_values(droplet, tuple(blob.flat)[0],
+                                                                                  dense_mask, dilute_mask,
+                                                                                  args.densechannels)
+                dense_intensity_mean.append(dense_means)
+                dil_intensity_mean.append(dil_means)
+                dense_intensity_sd.append(dense_sds)
+                dil_intensity_sd.append(dil_sds)
 
 
 if __name__ == "__main__":
     main(sys.argv[1:])
-