WIP

Script is working but segmentation results aren't perfect, maybe keeping the pre processing step is needed. To continue

WIP
8c384c68 · Laurent Guerard · 7003049f · 8c384c68
Commit 8c384c68 authored 4 months ago by Laurent Guerard
--- a/1_identify_fibers.py
+++ b/1_identify_fibers.py
@@ -5,8 +5,7 @@
 # TODO: are the imports RoiManager and ResultsTable needed when using the services?
 from ij import IJ, WindowManager as wm
 from ij.plugin import Duplicator, RoiEnlarger, RoiScaler
-from trainableSegmentation import WekaSegmentation
-from de.biovoxxel.toolbox import Extended_Particle_Analyzer
 from ij.measure import ResultsTable
 # Bio-formats imports
@@ -16,29 +15,47 @@ from loci.plugins.in import ImporterOptions
 # python imports
 import time
 import os
+import sys
+# TrackMate imports
+from fiji.plugin.trackmate import Logger, Model, SelectionModel, Settings, TrackMate
+from fiji.plugin.trackmate.action import LabelImgExporter
+from fiji.plugin.trackmate.cellpose import CellposeDetectorFactory
+from fiji.plugin.trackmate.cellpose.CellposeSettings import PretrainedModel
+from fiji.plugin.trackmate.features import FeatureFilter
+from fiji.plugin.trackmate.providers import (
+    SpotAnalyzerProvider,
+    SpotMorphologyAnalyzerProvider,
+)
+from fiji.plugin.trackmate.tracking.jaqaman import SparseLAPTrackerFactory
+from java.lang import Double
+from ch.epfl.biop.ij2command import Labels2Rois
 #@ String (visibility=MESSAGE, value="<html><b> Welcome to Myosoft - identify fibers! </b></html>") msg1
-#@ File (label="Select directory with classifiers", style="directory") classifiers_dir
 #@ File (label="Select directory for output", style="directory") output_dir
-#@ File (label="Select image file", description="select your image")  path_to_image
+#@ File (label="Select image file", description="select your image") path_to_image
+#@ File(label="Cellpose environment folder", style="directory", description="Folder with the cellpose env") cellpose_dir
 #@ Boolean (label="close image after processing", description="tick this box when using batch mode", value=False) close_raw
 #@ String (visibility=MESSAGE, value="<html><b> Morphometric Gates </b></html>") msg2
 #@ Integer (label="Min Area [um²]", value=10) minAr
 #@ Integer (label="Max Area [um²]", value=6000) maxAr
-#@ Float (label="Min Circularity", value=0.5) minCir
+#@ Double (label="Min Circularity", value=0.5) minCir
-#@ Float (label="Max Circularity", value=1) maxCir
+#@ Double (label="Max Circularity", value=1) maxCir
-#@ Float (label="Min solidity", value=0.0) minSol
+#@ Double (label="Min solidity", value=0.0) minSol
-#@ Float (label="Max solidity", value=1) maxSol
+#@ Double (label="Max solidity", value=1) maxSol
 #@ Integer (label="Min perimeter [um]", value=5) minPer
 #@ Integer (label="Max perimeter [um]", value=300) maxPer
 #@ Integer (label="Min min ferret [um]", value=0.1) minMinFer
 #@ Integer (label="Max min ferret [um]", value=100) maxMinFer
 #@ Integer (label="Min ferret AR", value=0) minFAR
 #@ Integer (label="Max ferret AR", value=8) maxFAR
-#@ Float (label="Min roundess", value=0.2) minRnd
+#@ Double (label="Min roundess", value=0.2) minRnd
-#@ Float (label="Max roundess", value=1) maxRnd
+#@ Double (label="Max roundess", value=1) maxRnd
 #@ String (visibility=MESSAGE, value="<html><b> Expand ROIS to match fibers </b></html>") msg3
-#@ Float (label="ROI expansion [microns]", value=1) enlarge
+#@ Double (label="ROI expansion [microns]", value=1) enlarge
 #@ String (visibility=MESSAGE, value="<html><b> channel positions in the hyperstack </b></html>") msg5
 #@ Integer (label="Membrane staining channel number", style="slider", min=1, max=5, value=1) membrane_channel
 #@ Integer (label="Fiber staining (MHC) channel number (0=skip)", style="slider", min=0, max=5, value=3) fiber_channel
@@ -163,45 +180,190 @@ def get_threshold_from_method(imp, channel, method):
    return lower_thr, upper_thr
-def apply_weka_model(model_path, imp, tiles_per_dim):
+def run_tm(
-    """apply a pretrained WEKA model to an ImagePlus
+    implus,
+    channel_seg,
+    cellpose_env,
+    seg_model,
+    diam_seg,
+    channel_sec=0,
+    quality_thresh=[0,0],
+    intensity_thresh=[0,0],
+    circularity_thresh=[0,0],
+    perimeter_thresh=[0,0],
+    feret_thresh=[0,0],
+    area_thresh=[0,0],
+    crop_roi=None,
+    use_gpu=True,
+):
+    """
+    Function to run TrackMate on open data
    Parameters
    ----------
-    model_path : string
+    implus : ImagePlus
-        path to the model file
+        ImagePlus on which to run the function
-    imp : ImagePlus
+    channel_seg : int
-        ImagePlus to apply the model to
+        Channel of interest
-    tiles_per_dim : integer
+    cellpose_env : str
-        tiles the imp to save RAM
+        Path to the cellpose environment
+    seg_model : PretrainedModel
+        Model to use for the segmentation
+    diam_seg : float
+        Diameter to use for segmentation
+    channel_sec : int, optional
+        Secondary channel to use for segmentation, by default 0
+    quality_thresh : float, optional
+        Threshold for quality filtering, by default 0
+    intensity_thresh : float, optional
+        Threshold for intensity filtering, by default 0
+    circularity_thresh : float, optional
+        Threshold for circularity filtering, by default 0
+    perimeter_thresh : float, optional
+        Threshold for perimeter filtering, by default 0
+    feret_thresh : float, optional
+        Threshold for Feret filtering, by default 0
+    area_thresh : float, optional
+        Threshold for area filtering, by default 0
+    crop_roi : ROI, optional
+        ROI to crop on the image, by default None
+    use_gpu : bool, optional
+        Boolean to use GPU or not, by default True
    Returns
    -------
    ImagePlus
-        the result of the WEKA segmentation. One channel per class.
+        Label image with the segmented objects
    """
-    segmentator = WekaSegmentation()
-    segmentator.loadClassifier( model_path )
-    result = segmentator.applyClassifier( imp, [tiles_per_dim, tiles_per_dim], 0, True ) #ImagePlus imp, int[x,y,z] tilesPerDim, int numThreads (0=all), boolean probabilityMaps
-    return result
+    # Get image dimensions and calibration
-def process_weka_result(imp):
+    dims = implus.getDimensions()
-    """apply myosoft pre-processing steps for the imp after WEKA classification to prepare it
+    cal = implus.getCalibration()
-    for ROI detection with the extended particle analyzer
+    # If the image has more than one slice, adjust the dimensions
+    if implus.getNSlices() > 1:
+        implus.setDimensions(dims[2], dims[4], dims[3])
+    # Set ROI if provided
+    if crop_roi is not None:
+        implus.setRoi(crop_roi)
+    # Initialize TrackMate model
+    model = Model()
+    model.setLogger(Logger.IJTOOLBAR_LOGGER)
+    # Prepare settings for TrackMate
+    settings = Settings(implus)
+    settings.detectorFactory = CellposeDetectorFactory()
+    # Configure detector settings
+    settings.detectorSettings["TARGET_CHANNEL"] = channel_seg
+    settings.detectorSettings["OPTIONAL_CHANNEL_2"] = channel_sec
+    settings.detectorSettings["CELLPOSE_PYTHON_FILEPATH"] = os.path.join(
+        cellpose_env, "python.exe"
+    )
+    settings.detectorSettings["CELLPOSE_MODEL_FILEPATH"] = os.path.join(
+        os.environ["USERPROFILE"], ".cellpose", "models"
+    )
+    settings.detectorSettings["CELLPOSE_MODEL"] = seg_model
+    settings.detectorSettings["CELL_DIAMETER"] = diam_seg
+    settings.detectorSettings["USE_GPU"] = use_gpu
+    settings.detectorSettings["SIMPLIFY_CONTOURS"] = True
+    settings.initialSpotFilterValue = -1.0
+    # Add spot analyzers
+    spotAnalyzerProvider = SpotAnalyzerProvider(1)
+    spotMorphologyProvider = SpotMorphologyAnalyzerProvider(1)
+    for key in spotAnalyzerProvider.getKeys():
+        settings.addSpotAnalyzerFactory(spotAnalyzerProvider.getFactory(key))
+    for key in spotMorphologyProvider.getKeys():
+        settings.addSpotAnalyzerFactory(spotMorphologyProvider.getFactory(key))
+    # Apply spot filters based on thresholds
+    if any(quality_thresh):
+        settings = set_trackmate_filter(settings, "QUALITY", quality_thresh)
+    if any(intensity_thresh):
+        settings = set_trackmate_filter(settings, "MEAN_INTENSITY_CH" + str(channel_seg), intensity_thresh)
+    if any(circularity_thresh):
+        settings = set_trackmate_filter(settings, "CIRCULARITY", circularity_thresh)
+    if any(area_thresh):
+        settings = set_trackmate_filter(settings, "AREA", area_thresh)
+    if any(perimeter_thresh):
+        settings = set_trackmate_filter(settings, "PERIMETER", perimeter_thresh)
+    if any(feret_thresh):
+        settings = set_trackmate_filter(settings, "FERET", feret_thresh)
+    print(settings)
+    # Configure tracker
+    settings.trackerFactory = SparseLAPTrackerFactory()
+    settings.trackerSettings = settings.trackerFactory.getDefaultSettings()
+    # settings.addTrackAnalyzer(TrackDurationAnalyzer())
+    settings.trackerSettings["LINKING_MAX_DISTANCE"] = 3.0
+    settings.trackerSettings["GAP_CLOSING_MAX_DISTANCE"] = 3.0
+    settings.trackerSettings["MAX_FRAME_GAP"] = 2
+    # Initialize TrackMate with model and settings
+    trackmate = TrackMate(model, settings)
+    trackmate.computeSpotFeatures(True)
+    trackmate.computeTrackFeatures(False)
+    # Check input validity
+    if not trackmate.checkInput():
+        sys.exit(str(trackmate.getErrorMessage()))
+        return
+    # Process the data
+    if not trackmate.process():
+        if "[SparseLAPTracker] The spot collection is empty." in str(
+            trackmate.getErrorMessage()
+        ):
+            return IJ.createImage(
+                "Untitled",
+                "8-bit black",
+                implus.getWidth(),
+                implus.getHeight(),
+                implus.getNFrames(),
+            )
+        else:
+            sys.exit(str(trackmate.getErrorMessage()))
+            return
+    # Export the label image
+    # sm = SelectionModel(model)
+    exportSpotsAsDots = False
+    exportTracksOnly = False
+    label_imp = LabelImgExporter.createLabelImagePlus(
+        trackmate, exportSpotsAsDots, exportTracksOnly, False
+    )
+    label_imp.setDimensions(1, dims[3], dims[4])
+    label_imp.setCalibration(cal)
+    implus.setDimensions(dims[2], dims[3], dims[4])
+    return label_imp
+def set_trackmate_filter(settings, filter_name, filter_value):
+    """Sets a TrackMate spot filter with specified filter name and values.
    Parameters
    ----------
-    imp : ImagePlus
+    settings : Settings
-        a single channel (= desired class) of the WEKA classification result imp
+        TrackMate settings object to which the filter will be added.
+    filter_name : str
+        The name of the filter to be applied.
+    filter_value : list
+        A list containing two values for the filter. The first value is
+        applied as an above-threshold filter, and the second as a below-threshold filter.
    """
-    IJ.run(imp, "8-bit", "")
+    filter = FeatureFilter(filter_name, filter_value[0], True)
-    IJ.run(imp, "Median...", "radius=3")
+    settings.addSpotFilter(filter)
-    IJ.run(imp, "Gaussian Blur...", "sigma=2")
+    filter = FeatureFilter(filter_name, filter_value[1], False)
-    IJ.run(imp, "Auto Threshold", "method=MaxEntropy")
+    settings.addSpotFilter(filter)
-    IJ.run(imp, "Invert", "")
+    return settings
 def delete_channel(imp, channel_number):
    """delete a channel from target imp
@@ -514,10 +676,6 @@ print("output_dir: ", str(output_dir))
 if not os.path.exists( str(output_dir) ):
    os.makedirs( str(output_dir) )
-classifiers_dir = fix_ij_dirs(classifiers_dir)
-primary_model = classifiers_dir + "/" + "primary.model"
-secondary_model = classifiers_dir + "/" + "secondary_central_nuclei.model"
 # update the log for the user
 IJ.log( "Now working on " + str(raw_image_title) )
 if raw_image_calibration.scaled() == False:
@@ -536,19 +694,14 @@ IJ.log( "MHC positive fiber channel = " + str(fiber_channel) )
 IJ.log( "sub-tiling = " + str(tiling_factor) )
 IJ.log( " -- settings used -- ")
-# image (pre)processing and segmentation (-> ROIs)
+# image (pre)processing and segmentation (-> ROIs)# imp, firstC, lastC, firstZ, lastZ, firstT, lastT
-membrane = Duplicator().run(raw, membrane_channel, membrane_channel, 1, 1, 1, 1) # imp, firstC, lastC, firstZ, lastZ, firstT, lastT
+imp_result = run_tm(raw, membrane_channel, cellpose_dir.getPath(), PretrainedModel.CYTO2, 30.0, area_thresh=[minAr, maxAr], circularity_thresh=[minCir, maxCir],
-preprocess_membrane_channel(membrane)
+        perimeter_thresh=[minPer, maxPer],
-weka_result1 = apply_weka_model(primary_model, membrane, tiling_factor )
+        # feret_thresh=[minMinFer, maxMinFer],
-delete_channel(weka_result1, 1)
+    )
-weka_result2 = apply_weka_model(secondary_model, weka_result1, tiling_factor )
+IJ.saveAs(imp_result, "Tiff", output_dir + "/" + raw_image_title + "_all_fibers_binary")
-delete_channel(weka_result2, 1)
-weka_result2.setCalibration(raw_image_calibration)
+sys.exit()
-process_weka_result(weka_result2)
-IJ.saveAs(weka_result2, "Tiff", output_dir + "/" + raw_image_title + "_all_fibers_binary")
-eda_parameters = [minAr, maxAr, minPer, maxPer, minCir, maxCir, minRnd, maxRnd, minSol, maxSol, minFAR, maxFAR, minMinFer, maxMinFer]
-raw.show() # EPA will not work if no image is shown
-run_extended_particle_analyzer(weka_result2, eda_parameters)
 # modify rois
 rm.hide()