DeepMoma Ground truth segment generation - nimwegenLab/moma GitHub Wiki
Generation of ground truth segments
Introduction
This page describes how the ground truth segments are generated, which are used during training of U-Net.
The ground truth segment (GT) is a binary mask for each input image with same size as the image, where a pixel value of 1 indicates that the pixel belongs to the cell (foreground pixel) and 0 does not belong to a cell (background pixel).
GT data is generated from the fluorescence (FL) color channel that accompanies the phase contrast images (if available for a given experiment).
The code for GT generation is written in Python and relies on the Scikit-Image library (https://scikit-image.org/).
GT segment generation
Segment generation is based on local thresholding and watershedding of the FL images and consists of the following steps, which are performed for each FL image in the TIFF stack. It consists of the following steps. A summary of the code is provided in the snippet below:
- Image smoothing/noise filtering:
Images are preprocessed for noise reduction using (1) Gaussian smoothing and (2) non-local means noise reduction.
Gaussian smoothing results in slightly larger segment masks, which match cell contours in the PhC image better (evaluated qualitatively by eye).
I use the non-local means method for noise-reduction to segment data from experiments, which where not deliberately recorded for GT generation and which have lower signal-to-noise ratio, because of lower exposure times. Non-local means is an edge preserving noise-filter. - Correcting for intensity inhomogeneity:
To handle strong intensity imbalances between adjacent cells and cases with high intensity at the channel exit (i.e. from cells that cluster outside the GL exit) I perform the following operations to make the background intensity more uniform: - Obtaining watershed markers:
Watershedding requires marker pixels with a unique integer value (label value) per cell (note: there can be multiple marker pixels for the same cell, which all have the same label value; these will yield a single region during watershedding).
I use two different methods to obtain wattershedding markers depending on whether (1) we are generating them from previous GT masks or (2) generating new GT segments without prior GT data. I describe both methods in the section below. - Watershedding:
We obtain cell segment using the markers by first generating a watershedding mask using local thresholding.
The preprocessed image is then watershed using the mask and markers. We set `watershed_line=True`, which puts a set of background pixels (with value 0) between adjacent/touching segments. This is needed so that we have separate binary masks for each cell during training. - Finally, the resulting segments are processed to remove holes that can occur in segments in noisy data. This is done separately for each segment of each image
Code summary as single Python function (this code is actually part the Python class that does the segmentation):
def segment_fl_image(input_fl_image, watershed_markers):
# preprocess image
image = skimage.filters.gaussian(input_fl_image, sigma=1.5)
sigma_orig = skimage.restoration.estimate_sigma(input_fl_image)
image = skimage.restoration.denoise_nl_means(image, h=sigma_orig * 5., fast_mode=True, patch_size=5, patch_distance=3, multichannel=False)
image = image - skimage.filters.gaussian(image, sigma=20)
# watershed image
local_threshholds = skimage.filters.threshold_local(image, block_size=15)
watershed_mask = image > local_threshholds
cell_labels = skimage.segmentation.watershed(-image, markers=watershed_markers, mask=watershed_mask, watershed_line=True) # we have to invert image intensity for watershedding
# filter holes in cell labels
cell_labels_final = np.zeros_like(cell_labels)
regions = skimage.measure.regionprops(cell_labels)
for region in regions:
label_mask = cell_labels == region.label
new_label_mask = skimage.morphology.remove_small_holes(label_mask)
cell_labels_final[new_label_mask] = region.label
return cell_labels_final
Generation of watershed markers:
There are two different situations in which I generate watershed markers: (1) When starting from previous GT data and (2) without pre-existing GT data.
-
Using pre-existing GT data:
In this case, we start from a pre-existing image with cells labels (called `previous_gt_labels`), where the pixel-values of each segment in the image has a unique label-value (integer value). I skeletonize the segments and use the resulting image with skeletonized pixels as watershed markers:skeletonized_segments = np.int16(skeletonize(np.uint8(preexisting_labels), method='lee')) markers = skeletonized_segments * preexisting_labels # go from binary skeleton values to skeletonized label values as markersThe reason for skeletonizing the previous segments (as opposed to e.g. keeping a single label-pixel at the centroid position) is that skeletonizing ensures that I will have pixels inside the target-region of the watershed mask. This is not garanteed by keeping only a single pixel at the centroid position, which can fall outside the target region of the watershed mask (e.g. for “banana”-shaped segments):
[EXAMPLE IMAGES]
-
Without pre-existing GT data:
A code-snippet how this is done, is provided below.
In method `get_peak_markers` we calculate a 1D intensity profile by tracing the “ridge” of maximal cell intensities along direction of the GL. The algorithm moves vertically through each pixel-row and shifts the column-position in the direction of the maximal intensity using a look-ahead of 3 rows (meaning it considers 3 pixel values along the vertical and diagonal directions). The shift in x-position is either -1, 0 and +1 pixels for each iteration, so that we step at maximum one pixel horizontally in the direction of the maximal intensity. This yields a contiguous trace of pixel values and positions, which approximately follows the maximal intensity in each row (barring sudden jumps in position of the maximal value):IMAGE, INTENSITY PROFILE
We determine the peak values and their index from the obtained intensity profile. We then check for each pair of peaks, if any profile-value in-between falls below a threshold value. If so, we consider the two peaks as belonging to separate cells. The threshold value is set as a (user-defined) fraction of the higher value of the two peaks being considered. We use the higher value of the two peak intensities, because for two cells with dissimilar intensities, the intensity-minimum between the two peaks will be strongly influenced (shift to higher values) by the brighter cell.
This is an excerpt from the Python class that does the segmentation, where the intensity profile and watershed markers are generated:
intensity_profile, profile_coordinates = self.get_intensity_profile(fl_image, self.max_center_offset) markers = self.get_peak_markers(fl_image, intensity_profile, profile_coordinates, self.cell_division_threshold_ratio) def get_intensity_profile(self, image, max_center_offset): """ This function moves along the maximal intensity of each row while using a horizontal step-size of one pixel. This effectively traces the maximum intensity of fluorescent cells and traverses the the low intensity regions between them. From this we obtain an intensity profile where the maxima correspond to the maximum intensity of cells and the minima correspond to the maximal intensity between cells. :param image: :param horizontal_image_center: :param max_center_offset: :return: """ horizontal_image_center = np.argmax(np.mean(image, axis=0)) coords = [] coord = [0, horizontal_image_center] index_range = range(image.shape[0]) image_padded = np.pad(image, ((0, 3), (0, 0)), mode='reflect') # pad at bottom of image; for the look ahead of detecting the next direction for ind in index_range: coords.append(coord) diagonal_left_coords = np.array([[coord[0] + 1, coord[1] - 1], [coord[0] + 2, coord[1] - 2], [coord[0] + 3, coord[1] - 3]]) center_coords = np.array([coord[0] + 1, coord[1](/nimwegenLab/moma/wiki/coord[0]-+-1,-coord[1), [coord[0] + 2, coord[1]], [coord[0] + 3, coord[1]]]) diagonal_right_coords = np.array([[coord[0] + 1, coord[1] + 1], [coord[0] + 2, coord[1] + 2], [coord[0] + 3, coord[1] + 3]]) left_max = np.max(image_padded[diagonal_left_coords[:, 0], diagonal_left_coords[:, 1]]) center_max = np.max(image_padded[center_coords[:, 0], center_coords[:, 1]]) right_max = np.max(image_padded[diagonal_right_coords[:, 0], diagonal_right_coords[:, 1]]) if left_max > center_max and left_max > right_max: y_shift = -1 if right_max > center_max and right_max > left_max: y_shift = 1 if center_max >= left_max and center_max >= right_max: y_shift = 0 if np.abs((coord[1] + y_shift) - horizontal_image_center) > max_center_offset: y_shift = 0 next_coord = [ind + 1, coord[1] + y_shift] coord = next_coord coords_arr = np.array(coords) return image[coords_arr[:, 0], coords_arr[:, 1]], np.roll(coords_arr, 1, axis=1) def get_peak_markers(self, fl_image, intensity_profile, profile_coordinates, cell_division_threshold_ratio): """ This method determines which of the peaks in intensity_profile belong to separate cells. It determines this by comparing neighboring pairs of peaks and testing if the intermediate intensity values fall below threshold, which is a user-defined ratio of the peak value of higher of the two. We use the higher peak of the two, because for adjacent cells with dissimilar brightness, the brighter cell (with higher peak), will determine the minimum intensity between the peaks. :param cell_division_threshold_ratio: :param intensity_profile: :param peak_inds: :param peak_vals: :return: """ padding_length = 1 # We have 0-pad the intensity profile added at the front and back, because 'find_peaks' will not detect maxima directly at the start/end of the profile (because technically they are not maxima). intensity_profile_padded = np.pad(intensity_profile, (padding_length, padding_length)) peaks = find_peaks(intensity_profile_padded, height=self.min_intensity, distance=self.peak_inter_distance, prominence=self.peak_prominence) peak_inds = peaks[0] - padding_length peak_vals = intensity_profile[peak_inds] peak_coords = profile_coordinates[peak_inds] output_peak_inds = [] current_peak_val = peak_vals[0] current_peak_ind = peak_inds[0] peak_labels = [] label_val = 1 for ind in range(len(peak_vals) - 1): next_peak_val = peak_vals[ind + 1] next_peak_ind = peak_inds[ind + 1] max_peak_val = np.max([current_peak_val, next_peak_val]) if np.any(intensity_profile[current_peak_ind:next_peak_ind] < max_peak_val * (1 - cell_division_threshold_ratio)): # add current peak, since the next one will start a new cell output_peak_inds.append(current_peak_ind) current_peak_val = next_peak_val current_peak_ind = next_peak_ind peak_labels.append(label_val) label_val += 1 else: # set the highest of the two peaks as current peak; this way we will keep only the highest peak until the if-condition is met current_peak_val = current_peak_val if current_peak_val == max_peak_val else next_peak_val current_peak_ind = current_peak_ind if current_peak_val == max_peak_val else next_peak_ind peak_labels.append(label_val) next_peak_val = peak_vals[-1] next_peak_ind = peak_inds[-1] max_peak_val = np.max([current_peak_val, next_peak_val]) if np.any(intensity_profile[current_peak_ind:next_peak_ind] < max_peak_val * (1 - cell_division_threshold_ratio)): output_peak_inds.append(next_peak_ind) label_val += 1 peak_labels.append(label_val) else: peak_labels.append(label_val) markers = np.zeros_like(fl_image, dtype=np.int) for ind, label_val in enumerate(peak_labels): markers[peak_coords[ind][1], peak_coords[ind][0]] = label_val return markers
Filtering/Curation of GT data
Again I distinguish between generation of GT data working with pre-existing GT data and without.
-
When using pre-existing GT data the resulting data is filtered automatically by comparing each new GT image with the corresponding previous one. For each image we check that (1) the number of segments in the image is the same and that (2) the difference between the limits of the bounding-boxes of eachsegment is below a user-selected threshold. As threshold for the difference of bounding-box limits, I used a threshold of +/-2 pixels in the horizontal direction (width) and +/-5 pixels in the vertical limits (height) ([SEE RESULTS SECTION FOR DISCUSSION]).
-
Without pre-existing GT data user-curation is required. GT is generated with one Python script per experiment. In this Python script the user specifies the the experiment-path, the GL numbers to use, the image ROI of in the GL TIFF-stack and the frames to used or conversely, which frames to ignore. The script produces 4 TIFF files per input GL:
images_fluorescence.tiff
images_phase_contrast.tiff
labels.tiff
segmented_cell_growth_animation.tiffThe first 3 files is the GT data that is used during training. The last file is used for user curation during the generation of the GT data. It contains a collage of the PhC and FL images with overlayed GT along with information on the final frame and the source frame:
User curation then proceeds iteratively with the user (1) specifing in the Python script the GL and frames to use, (2) generating the data, and (3) afterwards specifying which frames should be discarded, if they contain incorrect segments.
# %% import os from deep_moma.dataset_generation import DatasetConfiguration from deep_moma.dataset_generation import DatasetGenerator from deep_moma.dataset_generation import generate_dataset_config_from_template_with_time_ranges_with_ignores, write_rois_to_dataset_configs from deep_moma.dataset_generation import get_mm_data_path from deep_moma.segmentation_methods import ThresholdSegmenterV004, FluorescenceImageGlCenterer # %% get environment dataset_script_path = os.path.dirname(os.path.realpath(__file__)) data_input_base_directory = get_mm_data_path(dataset_script_path) # %% input settings # dataset_path_template: should contain {pos_ind} for the position index and {gl_ind:02d} for the growthlane index dataset_path_template = "/scicore/home/nimwegen/micuby52/Documents/01_work/DeepLearning/02_redo_training_data/06_with_reservoirs/00_lis_20210303/20210302_process_positions__combined__output__20210317-220506/Pos{pos_ind}/Pos{pos_ind}_GL{gl_ind}/20210303_VNG40_12h_2hAB_10h_1_MMStack_Pos{pos_ind}_GL{gl_ind}.tiff" dataset_config = DatasetConfiguration( mm_data_path=os.path.join(data_input_base_directory, dataset_path_template), number_of_color_channels=3, channel_number_phase_contrast=1, channel_number_fluorescence=2, roi=(30, None, None, None), ) # %% set output path script_name = os.path.splitext(os.path.basename(__file__))[0] dataset_config.dataset_path = os.path.join(dataset_script_path, script_name) # %% select a subset of positions to analyze. Should be representative of different cases e.g emtpy space at bottom. dataset_config.start_time = 0 dataset_config.stop_time = None # setting 'stop_timestep = None' will use 'mom.get_max_time()' positions_and_growthlanes = None output_width = 64 x_start = 21 x_stop = x_start + output_width rois = { 0: { # position index 22: (90, None, x_start, x_stop), # GL_ID: (ROI_LIMIT_TOP, ROI_LIMIT_BOTTOM, ROI_LIMIT_LEFT, ROI_LIMIT_RIGHT) }, } frame_range = (1, 480, 2) # range: (START, STOP, STEPSIZE) positions_and_growthlanes_and_timesteps = { # format: POS_ID: {GL_ID: [RANGE_TUPLE_1, INDIVIDUAL_FRAME_INDEX, ...], GL_ID: [RANGE_TUPLE_1, ...], ...} } 0: { # POS_ID 22: [frame_range], # GL_ID: [RANGE_1, RANGE_2, ...] }, } timesteps_to_ignore = { # format: POS_ID: {GL_ID: [RANGE_TUPLE_1, INDIVIDUAL_FRAME_INDEX, ...], GL_ID: [RANGE_TUPLE_1, ...], ...} } 0: { 22: [(81, 120, 1)], }, } dataset_configs = generate_dataset_config_from_template_with_time_ranges_with_ignores(dataset_config, positions_and_growthlanes_and_timesteps, timesteps_to_ignore) dataset_configs = write_rois_to_dataset_configs(dataset_configs, rois) # %% segmentation configuration segmenter = ThresholdSegmenterV004(gauss_sigma=1.5, block_size=15, cell_division_threshold=0.25) # %% define GL centerer; an object that will re-center the image gl_centerer = FluorescenceImageGlCenterer(vertical_detection_region_start=90) # %% generate datasets generator = DatasetGenerator(dataset_configs, run_interactive=False, segmenter=segmenter, show_debug_images=False) generator.generate_datasets()