choose_seeds.py

# This function was written as part of the analysis pipeline in 
# Lu et al 2018, IEEE ISBI: https://arxiv.org/abs/1803.07534
#
# The purpose of this function was to convert the predicted mask
# of cilia into an empirical probability density function that we
# would then sample "seeds" from to extract patches of ciliary motion
# (with high probability; hence, the empirical PDF) to feed to the
# downstream Conv-LSTM for classification. The reason we didn't
# just uniformly sample from the masked areas of cilia was to try
# and minimize uncertainty in the mask edges (hence, the distance
# transform to more heavily weight the interiors of the masks) while
# also acknowledging that different types of motion occur at the 
# cilia base vs the tips.

import numpy as np
import scipy.ndimage as spnd

def choose_seeds(mask, overlap, patchsize):
    """
    Randomly chooses a number of seeds for patch extraction.
    """
    # Create a distance transform of the mask, and restrict the results only
    # to those pixels within the mask. We'll use this as a sort of "spatial"
    # probability distribution for randomly choosing patch seeds.
    dist = spnd.distance_transform_edt(mask)

    # Need to first zero out the borders.
    width = int(patchsize / 2)
    dist[:width, :] = dist[-width:, :] = dist[:, :width] = dist[:, -width:] = 0

    # Now, proceed with sampling.
    cilia_indices = np.where(dist > 0)
    weights = dist[cilia_indices]
    indices = np.column_stack([cilia_indices[0], cilia_indices[1]])
    weights /= weights.sum()

    # How many patches do we actually extract?
    num_pixels = weights.shape[0]
    dense_patch_num = int((num_pixels / (patchsize ** 2))) # * (1 - overlap))
    if dense_patch_num == 0:
        return None
    ### TODO: More intelligently determine the number of patches.

    seeds = np.random.choice(num_pixels, size = dense_patch_num, replace = False, p = weights)
    ### TODO: Alter the probability density after each sampling to discourage
    # additional patches from being sampled right around the previous one.

    selected = indices[seeds]
    return selected