Source code for histomicstk.filters.shape.clog

import numpy as np


[docs] def clog(im_input, im_mask, sigma_min, sigma_max): """Constrained Laplacian of Gaussian filter. Takes as input a grayscale nuclear image and binary mask of cell nuclei, and uses the distance transform of the nuclear mask to constrain the LoG filter response of the image for nuclear seeding. Returns a LoG filter image of type float. Local maxima are used for seeding cells. Parameters ---------- im_input : array_like A hematoxylin intensity image obtained from ColorDeconvolution. Objects are assumed to be dark with a light background. im_mask : array_like A binary image where nuclei pixels have value 1/True, and non-nuclear pixels have value 0/False. sigma_min : double Minimum sigma value for the scale space. For blob detection, set this equal to minimum-blob-radius / sqrt(2). sigma_max : double Maximum sigma value for the scale space. For blob detection, set this equal to maximum-blob-radius / sqrt(2). Returns ------- im_log_max : array_like An intensity image containing the maximal LoG filter response across all scales for each pixel im_sigma_max : array_like An intensity image containing the sigma value corresponding to the maximal LoG response at each pixel. The nuclei/blob radius value for a given sigma can be estimated to be equal to sigma * sqrt(2). References ---------- .. [#] Y. Al-Kofahi et al "Improved Automatic Detection and Segmentation of Cell Nuclei in Histopathology Images" in IEEE Transactions on Biomedical Engineering,vol.57,no.4,pp.847-52, 2010. """ from scipy.ndimage import distance_transform_edt, gaussian_laplace # convert intensity image type to float im_input = im_input.astype(float) # generate distance map im_dmap = distance_transform_edt(im_mask) # compute max sigma at each pixel as 2 times the distance to background im_sigma_ubound = 2.0 * im_dmap # clip max sigma values to specified range im_sigma_ubound = np.clip(im_sigma_ubound, sigma_min, sigma_max) # initialize log filter response array MIN_FLOAT = np.finfo(im_input.dtype).min im_log_max = np.zeros_like(im_input) im_log_max[:, :] = MIN_FLOAT im_sigma_max = np.zeros_like(im_input) # Compute maximal LoG filter response across the scale space sigma_start = int(np.floor(sigma_min)) sigma_end = int(np.ceil(sigma_max)) sigma_list = np.linspace(sigma_start, sigma_end, int(sigma_end - sigma_start + 1)) for sigma in sigma_list: # generate normalized filter response im_log_cur = sigma**2 * gaussian_laplace(im_input, sigma, mode='mirror') # constrain LoG response im_log_cur[im_sigma_ubound < sigma] = MIN_FLOAT # update maxima max_update_pixels = np.where(im_log_cur > im_log_max) if len(max_update_pixels[0]) > 0: im_log_max[max_update_pixels] = im_log_cur[max_update_pixels] im_sigma_max[max_update_pixels] = sigma # replace min floats im_log_max[im_log_max == MIN_FLOAT] = 0 return im_log_max, im_sigma_max