Initial commit

Mainly copied from spots-in-fibers.py of imcftodo-1089

2d_spots_in_fibers.py

+#@ 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