count_fluid_drop_areas.py

from nd2reader import ND2Reader
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import scipy.ndimage as ndi
import skimage.filters as F
import skimage.morphology as M
import skimage.measure as Ms
import skimage.exposure as E
from skimage.morphology.selem import disk
from skimage.morphology import remove_small_objects

import matplotlib as mpl
mpl.rcParams['figure.figsize'] = [10, 10]

def remove_big_objects(ar, max_size=64, connectivity=1, in_place=False):
    if in_place:
        out = ar
    else:
        out = ar.copy()

    if max_size == 0:  # shortcut for efficiency
        return out

    if out.dtype == bool:
        selem = ndi.generate_binary_structure(ar.ndim, connectivity)
        ccs = np.zeros_like(ar, dtype=np.int32)
        ndi.label(ar, selem, output=ccs)
    else:
        ccs = out

    try:
        component_sizes = np.bincount(ccs.ravel())
    except ValueError:
        raise ValueError("Negative value labels are not supported. Try "
                         "relabeling the input with `scipy.ndimage.label` or "
                         "`skimage.morphology.label`.")
    too_big = component_sizes > max_size
    too_big_mask = too_big[ccs]
    out[too_big_mask] = 0

    return out
  

base = Path("./2021-06-15 GBP1 mutation/2021-06-15 GBP1 mutation/")

def read_nd_img(path, index):
    with ND2Reader(path) as images:
        img = images[index]
    img = np.array(img)
    return img

img = read_nd_img(str(base/"C589A/4.nd2"), 0)
plt.imshow(img)

m = E.adjust_gamma(img, 2)
m = E.rescale_intensity(m, out_range=(0, 1))
plt.imshow(m)

th = 0.9
mask = m > th
mask = remove_small_objects(mask, min_size=5)
#mask = M.closing(mask, selem=disk(5))
max_size = 500
mask = remove_big_objects(mask, max_size)
#ccs = Ms.regionprops(Ms.label(mask))
#print([cc.area for cc in ccs])
plt.imshow(mask)

out_base = Path("./output")
if not out_base.exists():
    out_base.mkdir()
for dir_ in base.glob("*"):
    print(dir_)
    out_dir = out_base / dir_.name
    if not out_dir.exists():
        out_dir.mkdir()
    for path in dir_.glob("*.nd2"):
        print(path)
        img = read_nd_img(str(path), 0)
        m = E.adjust_gamma(img, 2)
        m = E.rescale_intensity(m, out_range=(0, 1))
        th = 0.9
        mask = m > th
        mask = remove_small_objects(mask, min_size=5)
        #mask = M.closing(mask, selem=disk(10))
        max_size = 500
        mask = remove_big_objects(mask, max_size)
        fig, axes = plt.subplots(figsize=(20,10), ncols=2)
        axes[0].imshow(img)
        axes[1].imshow(mask)
        fig.savefig(str(out_dir/(path.name+".jpg")))
        ccs = Ms.regionprops(Ms.label(mask))
        with open(out_dir/(path.name+".tsv"), 'w') as f:
            f.write("center\tarea\n")
            for cc in ccs:
                outline = "\t".join([str(i) for i in [cc.centroid, cc.area]]) + "\n"
                f.write(outline)