diff --git a/2d_spots_in_fibers.py b/2d_spots_in_fibers.py new file mode 100644 index 0000000000000000000000000000000000000000..3757c8e19cf1daad0f47a44e30e94fba71558b2e --- /dev/null +++ b/2d_spots_in_fibers.py @@ -0,0 +1,768 @@ +#@ OpService ops +#@ CommandService command +#@ File (label="select image", description="select your input image", style=file) path_to_image +#@ Integer (label="select image file series", description="leave 1 if not needed", value=1, min=1, max=20, stepSize=1, persist=false, style=slider) series_number +#@ Integer (label="Fiber Segmentation", description="select the Fiber Segmentation RoiSet .zip file") fiber_segmentation_roiset +#@ Integer (label="DAPI channel", description="select the DAPI channel", min=1, max=20, stepSize=1, persist=true, style=slider) dapi_channel_number +#@ Integer (label="DAPI threshold", description="0 = Auto", min=0) threshold +#@ String (label="processing channels", description="comma separated list of channels, e.g. 2,3,6,7", value="1,2,3,4") processing_channels_string +#@ String (label="Spot detection quality threshold", description="comma separated list of values, one per channel", value="80,80,80,80" ) quality_thresholds_string + + +# trackmate imports +from fiji.plugin.trackmate import Settings +from fiji.plugin.trackmate import TrackMate +from fiji.plugin.trackmate.cellpose import CellposeDetectorFactory +from fiji.plugin.trackmate.cellpose.CellposeSettings import PretrainedModel +from fiji.plugin.trackmate.detection import DogDetectorFactory +from fiji.plugin.trackmate.action import LabelImgExporter +# ij imports +from ij import IJ +from ij.plugin import ZProjector, ImageCalculator +from ij.plugin.frame import RoiManager +from ij.measure import ResultsTable +# Bio-Formats imports +from loci.plugins import BF +from loci.common import Region +from loci.plugins.in import ImporterOptions +from loci.formats import ImageReader, MetadataTools, Memoizer +# MorpholibJ imports +from inra.ijpb.binary import BinaryImages +# CLIJ imports +from net.haesleinhuepf.clij2 import CLIJ2 +# imglib2 imports +from net.imglib2.img.display.imagej import ImageJFunctions +from net.imglib2.roi import Regions +from net.imglib2.roi.labeling import LabelRegions +from net.imglib2.algorithm.labeling import ConnectedComponents +# BIOP imports +from ch.epfl.biop.ij2command import Labels2Rois +# python imports +import os + + +def get_parentdir_filename_ext_from_path(path): + """ returns parent directory, filename, and file extension from a given path + + Parameters + ---------- + path : String + the full path to a file + + Returns + ------- + String + the parent directory, the filename without extension and the extension split after the dot + """ + parent_dir = os.path.dirname(path) + base = os.path.basename(path) + (filename, ext) = os.path.splitext(base) + + return (parent_dir, filename, ext) + + +def tile_image(image_width, image_height, tile_width, tile_height, overlap): + """ Calculate tiles of an image of given its dimensions. Calculates tile start coordinates and tile dimensions including overlap, + which can directly be used in Bio-Formats "crop" function. Also calculates tile crop start coordinates and dimensions + to remove the overlap and the crops start coordinates to re-assemble the original size image. Useful to e.g. segment a large + already stitched image tile by tile. + + Parameters + ---------- + image_width : int + the input images width in pixels + image_height : int + the input images height in pixels + tile_width : int + the desired tile width in pixels + tile_height : int + the desired tile height in pixels + overlap : int + the desired tile overlap in pixels + + Returns + ------- + list + all_tile_regions: tile coordinates and dimensions including overlap + all_crop_regions: tile coordinates and dimensions to remove the overlap again + all_tile_start_no_overlap: start coordinates of the tiles without overlap to re-assemble the original size image + """ + # TODO: this function likely needs refactoring to be more clear and readable + number_of_columns = image_width // tile_width # is a floored division + column_remainder_px = image_width % tile_width + number_of_rows = image_height // tile_height + row_remainder_px = image_height % tile_height + all_tile_regions = [] + all_crop_regions = [] + all_tile_start_no_overlap = [] + + for current_row in range(number_of_rows): + if current_row == number_of_rows -1 : + tile_height_boarder_addition = row_remainder_px # increase tile height in the last row + else: + tile_height_boarder_addition = 0 + + for current_column in range (number_of_columns): + if current_column == number_of_columns - 1 : + tile_width_boarder_addition = column_remainder_px # increase tile width in the last column + else: + tile_width_boarder_addition = 0 + + tile_start_x_no_overlap = current_column * tile_width + tile_start_y_no_overlap = current_row * tile_height + + if current_column > 0: + tile_start_x = current_column * tile_width - overlap + crop_start_x = overlap + else: + tile_start_x = 0 + crop_start_x = 0 + + if current_row > 0: + tile_start_y = current_row * tile_height - overlap + crop_start_y = overlap + else: + tile_start_y = 0 + crop_start_y = 0 + + if current_column == 0 or current_column == number_of_columns - 1: + current_tile_width = tile_width + overlap + tile_width_boarder_addition + else: + current_tile_width = tile_width + 2 * overlap + + if current_row == 0 or current_row == number_of_rows - 1: + current_tile_height = tile_height + overlap + tile_height_boarder_addition + else: + current_tile_height = tile_height + 2 * overlap + + tile_region = [tile_start_x, tile_start_y, current_tile_width, current_tile_height] + crop_region = [crop_start_x, crop_start_y, current_tile_width - overlap, current_tile_height - overlap] + tile_start_no_overlap = [tile_start_x_no_overlap, tile_start_y_no_overlap] + + all_tile_regions.append(tile_region) + all_crop_regions.append(crop_region) + all_tile_start_no_overlap.append(tile_start_no_overlap) + + return all_tile_regions, all_crop_regions, all_tile_start_no_overlap + + +def create_empty_image(bit_depth, image_width, image_height): + """creates an empty image of given dimensions. The image is 2D, filled with black pixels of value 0. + + Parameters + ---------- + bit_depth : int + the image bit dept + image_width : int + image width in px + image_height : int + image height in px + + Returns + ------- + ImagePlus + the empty image + """ + + canvas = IJ.createImage("mosaic", str(bit_depth) + "-bit black", image_width, image_height, 1) + + return canvas + + +def write_bf_memoryfile(path_to_file): + """Write a BF memo-file so subsequent access to the same file is faster. + + The Bio-Formats memo-file is written next to the image file (i.e. in the + same folder as the given file). + + Parameters + ---------- + string + The full path to the image file. + """ + reader = Memoizer(ImageReader()) + reader.setId(path_to_file) + reader.close() + + +def get_ome_metadata(path_to_image, series_number): + """Get some metadata fields using BF from an image on disk from a given series. + Gets number of channels, the image width, image height and bit depth. + + Parameters + ---------- + path_to_image : str + full path to the image + series_number : int + the Bio-Formats series number + + Returns + ------- + int + n_channels : the number of channels + int + image_width : the image width in px + int + image_height : the image height in px + int + bit_depth : the image bit depth + """ + series_number -= 1 # python starts counting from 0 + reader = ImageReader() + ome_meta = MetadataTools.createOMEXMLMetadata() + reader.setMetadataStore(ome_meta) + reader.setId(path_to_image) + reader.setSeries(series_number) + n_channels = reader.getSizeC() + image_width = reader.getSizeX() + image_height = reader.getSizeY() + bit_depth = reader.getBitsPerPixel() + + reader.close() + + return n_channels, image_width, image_height, bit_depth + + +def BFopen_image(path_to_image, channel_number, series_number, region=None, z_slice_number=None): + """Use Bio-Formats to open a given image or subset. + + Parameters + ---------- + path_to_image : str + full path to the image file on disk + channel_number : int + the number of the channel to open + series_number : int + the Bio-Formats series number to open + region : list, optional + Bio-Formats crop region, by default None. format = [start_x, start_y, width in px, height in px] + z_slice_number : int, optional + the number of the z-slice to open, by default None + + Returns + ------- + ImagePlus + the requested image + """ + series_number -= 1 # python starts counting from 0 + channel_number -= 1 # python starts counting from 0 + options = ImporterOptions() + options.setId(path_to_image) + options.setColorMode(ImporterOptions.COLOR_MODE_GRAYSCALE) + options.setSplitChannels(False) + options.setSplitFocalPlanes(False) + options.setSplitTimepoints(False) + options.setSeriesOn(series_number, True) + + if region is not None: + options.setCrop(True) + options.setCropRegion(series_number, Region(region[0], region[1], region[2], region[3])) # x,y,w,h + + options.setSpecifyRanges(True) + options.setCBegin(series_number, channel_number) + options.setCEnd(series_number, channel_number) + options.setCStep(series_number, 1) + + if z_slice_number is not None: + z_slice_number -= 1 # python starts counting from 0 + options.setZBegin(series_number, z_slice_number) + options.setZEnd(series_number, z_slice_number) + options.setZStep(series_number, 1) + + imps = BF.openImagePlus(options) + + return imps[0] + + +def run_trackmate_dog_spot_detector(imp, quality_threshold): + """Run TrackMates DoG detector with a given quality threshold + on a target image + + Parameters + ---------- + imp : ImagePlus + the input image + quality_threshold : float or double + the spot detection quality threshold + + Returns + ------- + ImagePlus + a label image of the detected spots. + """ + + cal = imp.getCalibration() + settings = Settings(imp) + + settings.detectorFactory = DogDetectorFactory() + settings.detectorSettings['DO_SUBPIXEL_LOCALIZATION'] = True + settings.detectorSettings['RADIUS'] = 0.15 # type = double + settings.detectorSettings['TARGET_CHANNEL'] = 0 + settings.detectorSettings['THRESHOLD'] = quality_threshold # type = double + settings.detectorSettings['DO_MEDIAN_FILTERING'] = False + + trackmate = TrackMate(settings) + trackmate.computeSpotFeatures(False) + trackmate.computeTrackFeatures(False) + + if not trackmate.checkInput(): + print( str( trackmate.getErrorMessage() ) ) + + if not trackmate.process(): + print( str( trackmate.getErrorMessage() ) ) + + exportSpotsAsDots = True + exportTracksOnly = False + label_imp = LabelImgExporter.createLabelImagePlus(trackmate, exportSpotsAsDots, exportTracksOnly, False) + label_imp.setCalibration(cal) + imp.close() + + return label_imp + + +def run_trackmate_cellpose_detector(imp): + """Run TrackMates CellPose detector on a target image + + Parameters + ---------- + imp : ImagePlus + the input image + + Returns + ------- + ImagePlus + a label image of the detected regions. + """ + + cal = imp.getCalibration() + settings = Settings(imp) + + # TODO: expose more parameters + settings.detectorFactory = CellposeDetectorFactory() + settings.detectorSettings['TARGET_CHANNEL'] = 0 + settings.detectorSettings['OPTIONAL_CHANNEL_2'] = 0 + settings.detectorSettings['CELLPOSE_PYTHON_FILEPATH'] = 'S:/cellpose_env/python.exe' + settings.detectorSettings['CELLPOSE_MODEL_FILEPATH'] = '' + settings.detectorSettings['CELLPOSE_MODEL'] = PretrainedModel.CYTO + settings.detectorSettings['CELL_DIAMETER'] = 30.0 + settings.detectorSettings['USE_GPU'] = True + settings.detectorSettings['SIMPLIFY_CONTOURS'] = True + + trackmate = TrackMate(settings) + trackmate.computeSpotFeatures(False) + trackmate.computeTrackFeatures(False) + + if not trackmate.checkInput(): + print( str( trackmate.getErrorMessage() ) ) + + if not trackmate.process(): + print( str( trackmate.getErrorMessage() ) ) + + exportSpotsAsDots = False + exportTracksOnly = False + label_imp = LabelImgExporter.createLabelImagePlus(trackmate, exportSpotsAsDots, exportTracksOnly, False) + label_imp.setCalibration(cal) + imp.close() + + return label_imp + + +def get_threshold_from_method(imp, method): + """returns the threshold value of chosen IJ AutoThreshold method from a stack (hyperstack?) + + Parameters + ---------- + imp : ImagePlus (maybe also accepts ij2 dataset, img?) + the image from which to get the threshold value + method : string + the AutoThreshold method to use, e.g. "otsu", "yen", ... + + Returns + ------- + int + """ + histogram = ops.run("image.histogram", imp) + threshold_value = ops.run("threshold.%s" % method, histogram) + threshold_value = int(round(threshold_value.get())) + + return threshold_value + + +def convert_to_binary(imp, threshold=1, scale_binary=True): + """Convert grayscale image to a mask / binary image + + Parameters + ---------- + imp : ImagePlus + the input image + threshold : int, optional + the threshold above which pixels are considered foreground, by default 1 + scale_binary : bool, optional + scale the mask image from 0/255 to 0/1, by default True + + Returns + ------- + ImagePlus + the mask / binary image + """ + mask_imp = imp.duplicate() + bitdepth = mask_imp.getBitDepth() + if bitdepth > 16.0: + IJ.run("Conversions...", " ") + IJ.run(mask_imp, "16-bit", "") + IJ.run("Conversions...", "scale") + + IJ.setRawThreshold(mask_imp, threshold, (2**bitdepth)-1, None) + IJ.run(mask_imp, "Convert to Mask", "") + if scale_binary == True: + IJ.run(mask_imp, "Divide...", "value=255") + IJ.run(mask_imp, "Enhance Contrast", "saturated=0.35") + mask_imp.changes = False + + return mask_imp + + +def erode_labelimp_gpu(label_imp, radius, relabel_islands=True): + """Erode a label image on the GPU. + + Parameters + ---------- + label_imp : ImagePlus + the input image + radius : float + radius used for erosion + relabel_islands : bool, optional + if True, split objects will get new labels each , by default True + + Returns + ------- + ImagePlus + the eroded label image + """ + + clij2 = CLIJ2.getInstance() + labels_input = clij2.push(label_imp) + labels_destination = clij2.create(labels_input) + radius = 1.0 + clij2.erodeLabels(labels_input, labels_destination, radius, relabel_islands) + eroded_labelimp = clij2.pull(labels_destination) + clij2.release(labels_input) + clij2.release(labels_destination) + label_imp.close() + + return eroded_labelimp + + +def remove_tileoverlap(imp, crop_region): + """Crop an image to specified dimensions. + The operation is destructive, the original image is lost. + + Parameters + ---------- + imp : ImagePlus + the input image + crop_region : list or array of int + the crop region in the format of [start_x, start_y, crop_width, crop_height] + + Returns + ------- + ImagePlus + the cropped image + """ + + start_x = crop_region[0] + start_y = crop_region[1] + crop_width = crop_region[2] + crop_height = crop_region[3] + imp.setRoi(start_x, start_y, crop_width, crop_height) + cropped_imp = imp.crop() + imp.close() + + return cropped_imp + + +def copy_and_paste_image(source, target, start_coordinates_xy): + """Copy the source image and paste it into the target image + + Parameters + ---------- + source : ImagePlus + the source image + target : ImagePlus + the target image + start_coordinates_xy : list or array of int + the x and y coordinates in the target image where the source image will be pasted + """ + source.copy() + start_x = start_coordinates_xy[0] + start_y = start_coordinates_xy[1] + target.paste(start_x, start_y) + source.close() + + +def merge_touching_labels(label_imp, bit_depth): + """merge touching labels in a label image so they have the same label ID. + + Parameters + ---------- + label_imp : ImagePlus + the label image with touching labels of different label ID + bit_depth : int + the image bit depth + + Returns + ------- + ImagePlus + the label image with merged touching labels + """ + + binary_imp = BinaryImages.binarize(label_imp) + connectivity = 8 + merged_label_imp = BinaryImages.componentsLabeling(binary_imp, connectivity, bit_depth) + label_imp.close() + binary_imp.close() + + return merged_label_imp + + +def save_image_as_IJtif(imp, filename, suffix, target): + """Save the image as ImageJ-Tiff + + Parameters + ---------- + imp : ImagePlus + the input Image to be saved + filename : str + the image filename without extension + suffix : int, float or str + a suffix that will be added to the filename as a string + target : str + the target directory in with to save the image + """ + + savename = filename + "_" + str(suffix) + ".tif" + savepath = os.path.join(target, savename) + IJ.log("now saving: " + str(savepath)) + IJ.saveAs(imp, "Tiff", savepath) + + +def measure_label_size_shape_2d(label_imp): + rm = RoiManager.getInstance() + if not rm: + rm = RoiManager() + rm.reset() + rt = ResultsTable() + + command.run( Labels2Rois , False , 'imp' , label_imp , 'rm', rm).get() + IJ.run("Set Measurements...", "area centroid perimeter shape feret's redirect=None decimal=3") + rm.runCommand(label_imp,"Deselect") + rm.runCommand(label_imp,"Measure") + rm.reset() + IJ.renameResults("Results", "Size and shape features") + results_table = rt.getResultsTable("Size and shape features") + + + return results_table + + +def convert_labelimage_to_imglib2regions(imp): + """Convert a ImagePlus label image to imglib2 regions + + Parameters + ---------- + imp : ImagePlus + the input label image + + Returns + ------- + Imglib2 LabelRegions + the imglib2 label regions + """ + + # fist convert ImagePlus to img as suggested by curtis: + # from https://forum.image.sc/t/how-to-wrap-any-kind-of-imageplus-to-an-imglib2-img-floattype/178/6 + wrapImg = ImageJFunctions.wrap(imp) + # workaround to convert a label image into an imglib2 ImgLabeling (which contains not only an image) + img_labeling = ops.labeling().cca(wrapImg, ConnectedComponents.StructuringElement.EIGHT_CONNECTED) + # img = img_labeling.getIndexImg() # the img label image + # ImageJFunctions.show(img) + regions = LabelRegions(img_labeling) + + return regions + + +def convert_labelimage_to_binary(label_imp, scale_binary=True): + """Convert a label image to a mask / binary image + + Parameters + ---------- + label_imp : ImagePlus + the input label image + scale_binary : bool, optional + scale the mask image from 0/255 to 0/1, by default True + + Returns + ------- + ImagePlus + the mask / binary image + """ + binary_imp = BinaryImages.binarize(label_imp) + if scale_binary == True: + IJ.run(binary_imp, "Divide...", "value=255") + + return binary_imp + + +def measure_intensity_sum(label_imp, target_imp): + """Measure the sum intensity for each label in a label image on another target image. + Uses imglib2 and ops for this and works with both 2D and 3D label images. + + Parameters + ---------- + label_imp : ImagePlus + the input label image + target_imp : ImagePlus + the target image in which to measure + + Returns + ------- + array + label id and corresponding sum intensity [label_id, sum_intensity] + """ + label_id = [] + sum_intensity = [] + target_img = ImageJFunctions.wrap(target_imp) # convert ImagePlus to img to use ops/region measurements on it + regions = convert_labelimage_to_imglib2regions(label_imp) + + for region in regions: + label_id.append(region.getLabel() + 1) # region.getlabel() starts counting from 0. So the label with int 1 has the index 0. + input = Regions.sample(region, target_img) # returns an iterableInterval + sum = ops.stats().sum(input).getRealDouble() + sum_intensity.append(sum) + # other measurements see https://forum.image.sc/t/can-i-get-to-measurements-using-imagej-ops/4448/5 + result = [label_id, sum_intensity] + + return result + + +def add_results_to_resultstable(results_table, column, values): + """Add values to the ResultsTable starting from row 0 of a given column. + + Parameters + ---------- + results_table : ij.measure.ResultsTable + a reference of the IJ-ResultsTable + column : string + the column in which to add the values + values : list(int, double or float) + array with values to be added + """ + for index, value in enumerate(values): + results_table.setValue(column, index, value) + + +def save_results_table(results_table, filename, suffix, target): + """Save an IJ-ResultsTable as txt to a given location + + Parameters + ---------- + results_table : IJ-ResultsTable + an IJ-ResultsTable + filename : str + filename to include in the tables name + suffix : str, int or float + a suffix to be included in the tables name + target : str + the target directory in which to save the table + """ + + savename = filename + "_" + str(suffix) + ".txt" + savepath = os.path.join(target, savename) + results_table.save(savepath) + + +def save_labelimage_as_ijroiset(label_imp, filename, suffix, target): + """Save all labels i a label Images as IJ-ROIset to target location. + Only 2D, requires PTBIOP update site. + + Parameters + ---------- + label_imp : ImagePlus + the input label image + filename : str + the filename to include in the ROIset name + suffix : str, int or float + a suffix to be included in the ROIset name + target : str + the target directory in which to save the table + """ + rm = RoiManager.getInstance() + if not rm: + rm = RoiManager() + rm.reset() + command.run( Labels2Rois , False , 'imp' , label_imp , 'rm', rm).get() + savename = filename + "_" + str(suffix) + ".zip" + savepath = os.path.join(target, savename) + rm.runCommand("Save", savepath) + rm.reset() + + +def close_images(list_of_imps): + """Close given ImagePlus images + + Parameters + ---------- + list_of_imps : list of ImagePlus + list of ImagePlus images to be closed + """ + for imp in list_of_imps: + imp.changes = False + imp.close() + + +path_to_image = str(path_to_image).replace("\\","/") +parent_dir, filename, ext = get_parentdir_filename_ext_from_path(path_to_image) +write_bf_memoryfile(path_to_image) +n_channels, image_width, image_height, bit_depth = get_ome_metadata(path_to_image, series_number) +processing_channels_string = processing_channels_string.replace(" ", "") +processing_channels = processing_channels_string.split(",") +quality_thresholds_string = quality_thresholds_string.replace(" ", "") +quality_thresholds = quality_thresholds_string.split(",") + +# TODO: Get the fiber segmentation and convert to labelimage + +# threshold DAPI channel and convert to binary +dapi_channel = BFopen_image(path_to_image, dapi_channel_number, series_number, z_slice_number=dapi_channel_zslice) +if threshold <= 0: + threshold = get_threshold_from_method(dapi_channel, "otsu") +dapi_binary = convert_to_binary(dapi_channel, threshold) +close_images([dapi_channel]) + +# detect spots and count them per fiber +results_table = ResultsTable() +for index, channel in enumerate(processing_channels): + channel = int(channel) + quality_threshold = float(quality_thresholds[index]) + spots_channel = BFopen_image(path_to_image, channel, series_number) # can be a stack + spots_label_imp = run_trackmate_dog_spot_detector(spots_channel, quality_threshold) # spot detection is 3D + save_image_as_IJtif(spots_label_imp, filename, "spots_ch" + str(channel), parent_dir) + save_labelimage_as_ijroiset(spots_label_imp, filename, "spots_ch" + str(channel), parent_dir) + spots_binary_imp = convert_labelimage_to_binary(spots_label_imp) # can be a stack + spots_binary_imp_sum_proj = ZProjector.run(spots_binary_imp, "sum") # if 2 (n) spots perfecly overlap in z, the pixel value will be 2 (n) + dapi_positive_spots_binary = ImageCalculator.run(spots_binary_imp_sum_proj, dapi_binary, "Multiply create") + dapi_negative_spots_binary = ImageCalculator.run(spots_binary_imp_sum_proj, dapi_positive_spots_binary, "Subtract create") + + n_spots_per_fiber = measure_intensity_sum(fibers_label_imp, spots_binary_imp_sum_proj) + n_dapi_positive_spots_per_fiber = measure_intensity_sum(fibers_label_imp, dapi_positive_spots_binary) + n_dapi_negative_spots_per_fiber = measure_intensity_sum(fibers_label_imp, dapi_negative_spots_binary) + + # TODO: the column fiber label ID needs to be added only one time. Maybe check before if column exists. + add_results_to_resultstable(results_table, "fiber label ID", n_spots_per_fiber[0]) + add_results_to_resultstable(results_table, "n spots channel " + str(channel), n_spots_per_fiber[1]) + add_results_to_resultstable(results_table, "n spots dapi positive channel " + str(channel), n_dapi_positive_spots_per_fiber[1]) + add_results_to_resultstable(results_table, "n spots dapi negative channel " + str(channel), n_dapi_negative_spots_per_fiber[1]) + close_images([spots_channel, spots_label_imp, spots_binary_imp, spots_binary_imp_sum_proj, dapi_binary, dapi_positive_spots_binary, dapi_negative_spots_binary]) + +results_table.show("Spots") +save_results_table(results_table, filename, "Spots", parent_dir) +IJ.log("DONE") \ No newline at end of file