ColonyCounter/colonycounter.py

from typing import List, Generator, Optional, Tuple
from fileutils import FILE, dslice, search, DIRIN, DIROUT, DIRCACHE
import window
from math import pi
import cv2
import pandas as pd
import numpy as np
from scipy import ndimage as ndi
from skimage.filters import threshold_otsu, threshold_li
from skimage import feature, measure, restoration, segmentation
from pythreshold.global_th.entropy.kapur import kapur_threshold

PIXELS_TO_UM = 3.2  # 1 pixel equals 3.2 um

def rball(img, *argv, **args):
    return img - restoration.rolling_ball(img, *argv, **args)


def watershed_new(img, blur: int, min_thresh: int, dist: int, min_size: int, min_roundness: float, min_mean_brightness: int, ignore: List[Tuple[int, int]]):
    MAX_VAL = 255 if img.dtype == "uint8" else 65535
    BRIGHT_MUL = 65535/255 if img.dtype == "uint8" else 1

    c_img = cv2.cvtColor(img,cv2.COLOR_GRAY2RGB)

    if blur > 0:
        b_img = cv2.GaussianBlur(img, (blur, blur), 0)
        thresh_val = kapur_threshold(b_img)
        thresh_val = thresh_val if thresh_val >= min_thresh else min_thresh
        thresh = cv2.threshold(cv2.bitwise_not(b_img), MAX_VAL-thresh_val, MAX_VAL, cv2.THRESH_BINARY_INV)[1]
    else:
        thresh_val = kapur_threshold(img)
        thresh_val = thresh_val if thresh_val >= min_thresh else min_thresh
        thresh = cv2.threshold(cv2.bitwise_not(img), MAX_VAL-thresh_val, MAX_VAL, cv2.THRESH_BINARY_INV)[1]
    thresh = segmentation.clear_border(thresh)
    thresh = ndi.binary_fill_holes(thresh)

    distance = ndi.distance_transform_edt(thresh)
    if dist < 1:
        localMax = distance > distance.min()
    else:
        localMax = np.zeros_like(img, dtype=np.bool)
        localMax[tuple(feature.peak_local_max(distance, exclude_border=False, min_distance=dist, labels=thresh).T)] = True


    markers = ndi.label(localMax, structure=np.ones((3, 3)))[0]
    labels = segmentation.watershed(-img, markers, mask=thresh)

    grain_num = 1
    meta = {
        "outlines": np.full_like(img, MAX_VAL, dtype=img.dtype),
        "data": [[ "Grain","Area","Mean","Min","Max","Circ", "Eccentricity", "IntDen","RawIntDen","AR","Round","Solidity" ]]
    }

    for prop, label in list(zip(measure.regionprops(labels, intensity_image=img), np.unique(labels)[1:]))[:50]:
        bbox = prop["bbox"]

        if prop["area"]*PIXELS_TO_UM**2 < min_size \
           or 1-prop["eccentricity"] < min_roundness \
           or prop["mean_intensity"]*BRIGHT_MUL < min_mean_brightness*257 \
           or grain_num > 50 \
           or any([bbox[1] <= x <= bbox[3] and bbox[0] <= y <= bbox[2]
                   for (x, y) in ignore]):
            continue

        # Create mask of image size and mark labeled area
        mask = np.zeros(img.shape, dtype=img.dtype)
        mask[labels == label] = MAX_VAL
        # cnts = measure.find_contours(mask.copy())
        # print(cnts)
        # c_img[cnts] = (65535,65535)
        # measure.drawContours()

        # c_img = segmentation.mark_boundaries(c_img, mask, (0,255,255))
        # Create and draw contour around labeled area
        if cv2.__version__.startswith("3."):
            cnts = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL,
                                    cv2.CHAIN_APPROX_SIMPLE)[1]
        else:
            cnts = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL,
                                    cv2.CHAIN_APPROX_SIMPLE)[0]

        c_img = cv2.drawContours(c_img, cnts, -1, (0,MAX_VAL,MAX_VAL), 1)
        meta["outlines"] = cv2.drawContours(meta["outlines"], cnts, -1, (0,0), 1)
        c = max(cnts, key=cv2.contourArea)
        (x,y), r = cv2.minEnclosingCircle(c)
        x = int(float(x) + float(r)/2)
        y = int(float(y) + float(r)/2)
        cv2.putText(c_img, str(grain_num), (x, y), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 0, 255), 1)
        cv2.putText(meta["outlines"], str(grain_num), (x, y), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 0), 1)
        meta["data"].append([
            grain_num,
            round(prop["area"]*PIXELS_TO_UM**2, 4),
            round(prop["mean_intensity"]*BRIGHT_MUL, 4),
            round(prop["min_intensity"]*BRIGHT_MUL, 4),
            round(prop["max_intensity"]*BRIGHT_MUL, 4),
            (1 if (4 * pi * prop.area) / (prop.perimeter_crofton * prop.perimeter_crofton) >=1 else round((4 * pi * prop.area) / (prop.perimeter_crofton * prop.perimeter_crofton), 4)),
            round(1-prop["eccentricity"], 4), # 4*pi*(prop["area"])/prop["perimeter"]**2,
            round((prop["area"]*PIXELS_TO_UM**2)*prop["mean_intensity"]*BRIGHT_MUL, 4),
            round(prop["intensity_image"].sum()*BRIGHT_MUL, 4),
            (None if prop["minor_axis_length"] == 0 else round(prop["major_axis_length"]/prop["minor_axis_length"], 4)),
            (None if prop["major_axis_length"] == 0 else round(4*prop["area"]*PIXELS_TO_UM**2/(pi*prop["major_axis_length"]**2)/10, 4)),
            round(prop["solidity"], 4)
        ])
        grain_num += 1
    return c_img, meta

def main():
    """
    tif metadata documentation:
    threshold-low  -> lower bound bit range (0)
    threshold-high -> upper bound bit range (65535)
    https://docs.opencv.org/3.0-beta/modules/core/doc/operations_on_arrays.html#cv2.normalize
    """
    files = [FILE(p) for p in search([DIRIN]) if p.endswith(".tif")]
    if files != []:
        files.sort()
        normalize = ["normalized", (cv2.normalize, [], {"dst": None, "alpha": 0, "beta": 65535, "norm_type": cv2.NORM_MINMAX}, True)]
        to8bit = ["8bit", (cv2.convertScaleAbs, [] , {"alpha": 255.0/65535.0}, True)]
        rollingball = ["rball", (rball, [], {"radius": 50.0}, True)]

        wat = ["wshed", (watershed_new, [], {"blur": 0, "dist": 0, "min_thresh": 20, "min_size": 250, "min_roundness": 0.2, "min_mean_brightness": 50, "ignore": []}, False)]
        [f.opqueue.__init__([normalize, to8bit, rollingball, wat]) for f in files]
        [f.apply(dslice(f.opqueue,None,-1)) for f in files]

        window.init(files)
        window.reset()
        return True
    print(f"ERR: No files in {DIRIN}, please add some.")
    return False

if __name__ == '__main__':
    if main():
        window.run()