-
Laurent Guerard authoredLaurent Guerard authored
1_identify_fibers.py 24.12 KiB
# @ String (visibility=MESSAGE, value="<html><b> Welcome to Myosoft - identify fibers! </b></html>") msg1
# @ File (label="Select folder with your images", description="select folder with your images", style="directory") src_dir
# @ String(label="Extension for the images to look for", value="czi") filename_filter
# @ File (label="Select directory for output", style="directory") output_dir
# @ 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
# @ Double (label="Min Circularity", value=0.5) minCir
# @ Double (label="Max Circularity", value=1) maxCir
# @ Integer (label="Min perimeter [um]", value=5) minPer
# @ Integer (label="Max perimeter [um]", value=300) maxPer
# @ String (visibility=MESSAGE, value="<html><b> Expand ROIS to match fibers </b></html>") msg3
# @ Double (label="ROI expansion [microns]", value=1) enlarge_radius
# @ 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
# @ Integer (label="minimum fiber intensity (0=auto)", description="0 = automatic threshold detection", value=0) min_fiber_intensity
# @ CommandService command
# @ RoiManager rm
# @ ResultsTable rt
# this is a python rewrite of the original ijm published at
# https://github.com/Hyojung-Choo/Myosoft/blob/Myosoft-hub/Scripts/central%20nuclei%20counter.ijm
# ─── Requirements ─────────────────────────────────────────────────────────────
# List of update sites needed for the code
# * TrackMate-Cellpose
# * IMCF
# * PTBIOP
# * CLIJ-CLIJ2
# ─── Imports ──────────────────────────────────────────────────────────────────
# IJ imports
# TODO: are the imports RoiManager and ResultsTable needed when using the services?
import os
import sys
# Bio-formats imports
# from loci.plugins import BF
# from loci.plugins.in import ImporterOptions
# python imports
import time
from ch.epfl.biop.ij2command import Labels2CompositeRois
# TrackMate imports
from fiji.plugin.trackmate import Logger, Model, 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 ij import IJ
from ij import WindowManager as wm
from ij.measure import ResultsTable
from ij.plugin import Duplicator, ImageCalculator, RoiEnlarger
from imcflibs import pathtools
from imcflibs.imagej import bioformats as bf
from imcflibs.imagej import misc
# ─── Functions ────────────────────────────────────────────────────────────────
def fix_ij_options():
"""Put IJ into a defined state."""
# disable inverting LUT
IJ.run("Appearance...", " menu=0 16-bit=Automatic")
# set foreground color to be white, background black
IJ.run("Colors...", "foreground=white background=black selection=red")
# black BG for binary images and pad edges when eroding
IJ.run("Options...", "black pad")
# set saving format to .txt files
IJ.run("Input/Output...", "file=.txt save_column save_row")
# ============= DON’T MOVE UPWARDS =============
# set "Black Background" in "Binary Options"
IJ.run("Options...", "black")
# scale when converting = checked
IJ.run("Conversions...", "scale")
def fix_ij_dirs(path):
"""use forward slashes in directory paths.
Parameters
----------
path : string
a directory path obtained from dialogue or script parameter
Returns
-------
string
a more robust path with forward slashes as separators
"""
fixed_path = str(path).replace("\\", "/")
# fixed_path = fixed_path + "/"
return fixed_path
def fix_BF_czi_imagetitle(imp):
"""Fix the title of an image read using the bio-formats importer.
The title is modified to remove the ".czi" extension and replace
spaces with underscores.
Parameters
----------
imp : ij.ImagePlus
The image to be processed.
Returns
-------
string
The modified title of the image.
"""
image_title = os.path.basename(imp.getShortTitle())
# remove the ".czi" extension
image_title = image_title.replace(".czi", "")
# replace spaces with underscores
image_title = image_title.replace(" ", "_")
# remove any double underscores
image_title = image_title.replace("_-_", "")
# remove any double underscores
image_title = image_title.replace("__", "_")
# remove any "#" characters
image_title = image_title.replace("#", "Series")
return image_title
def do_background_correction(imp, gaussian_radius=20):
"""Perform background correction on an image.
This is done by applying a Gaussian blur to the image and then dividing the
original image by the blurred image.
Parameters
----------
imp : ij.ImagePlus
The image to be corrected.
gaussian_radius : int
The radius of the Gaussian filter to be used. Default value is 20.
Returns
-------
ij.ImagePlus
The background-corrected image.
"""
imp_bgd = imp.duplicate()
IJ.run(
imp_bgd,
"Gaussian Blur...",
"sigma=" + str(gaussian_radius) + " scaled",
)
return ImageCalculator.run(imp, imp_bgd, "Divide create 32-bit")
def get_threshold_from_method(imp, channel, method):
"""Get the value of automated threshold method.
Returns the threshold value of chosen IJ AutoThreshold method in desired channel.
Parameters
----------
imp : ij.ImagePlus
the imp from which to get the threshold value
channel : integer
the channel in which to get the treshold
method : string
the AutoThreshold method to use
Returns
-------
list
the upper and the lower threshold (integer values)
"""
imp.setC(channel) # starts at 1
ip = imp.getProcessor()
ip.setAutoThreshold(method + " dark")
lower_thr = ip.getMinThreshold()
upper_thr = ip.getMaxThreshold()
ip.resetThreshold()
return lower_thr, upper_thr
def run_tm(
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],
area_thresh=[0, 0],
crop_roi=None,
use_gpu=True,
):
"""
Function to run TrackMate on open data, applying filters to spots.
Parameters
----------
implus : ij.ImagePlus
ImagePlus on which to run the function
channel_seg : int
Channel of interest
cellpose_env : str
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 None
intensity_thresh : float, optional
Threshold for intensity filtering, by default None
circularity_thresh : float, optional
Threshold for circularity filtering, by default None
perimeter_thresh : float, optional
Threshold for perimeter filtering, by default None
area_thresh : float, optional
Threshold for area filtering, by default None
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
-------
ij.ImagePlus
Label image with the segmented objects
"""
# Get image dimensions and calibration
dims = implus.getDimensions()
cal = implus.getCalibration()
# 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)
# 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
----------
settings : Settings
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.
"""
filter = FeatureFilter(filter_name, filter_value[0], True)
settings.addSpotFilter(filter)
filter = FeatureFilter(filter_name, filter_value[1], False)
settings.addSpotFilter(filter)
return settings
def delete_channel(imp, channel_number):
"""Delete a channel from target imp.
Parameters
----------
imp : ij.ImagePlus
the imp from which to delete target channel
channel_number : integer
the channel number to be deleted. starts at 0.
"""
imp.setC(channel_number)
IJ.run(imp, "Delete Slice", "delete=channel")
def measure_in_all_rois(imp, channel, rm):
"""Gives measurements for all ROIs in ROIManager.
Measures in all ROIS on a given channel of imp all parameters that are set in IJ "Set Measurements".
Parameters
----------
imp : ij.ImagePlus
the imp to measure on
channel : integer
the channel to measure in. starts at 1.
rm : RoiManager
a reference of the IJ-RoiManager
"""
imp.setC(channel)
rm.runCommand(imp, "Deselect")
rm.runCommand(imp, "Measure")
def change_all_roi_color(rm, color):
"""Cchange the color of all ROIs in the RoiManager.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
color : string
the desired color. e.g. "green", "red", "yellow", "magenta" ...
"""
number_of_rois = rm.getCount()
for roi in range(number_of_rois):
rm.select(roi)
rm.runCommand("Set Color", color)
def change_subset_roi_color(rm, selected_rois, color):
"""Change the color of selected ROIs in the RoiManager.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
selected_rois : array
ROIs in the RoiManager to change
color : string
the desired color. e.g. "green", "red", "yellow", "magenta" ...
"""
rm.runCommand("Deselect")
rm.setSelectedIndexes(selected_rois)
rm.runCommand("Set Color", color)
rm.runCommand("Deselect")
def show_all_rois_on_image(rm, imp):
"""Shows all ROIs in the ROiManager on imp.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
imp : ij.ImagePlus
the imp on which to show the ROIs
"""
imp.show()
rm.runCommand(imp, "Show All")
def save_all_rois(rm, target):
"""Save all ROIs in the RoiManager as zip to target path.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
target : string
the path in to store the ROIs. e.g. /my-images/resulting_rois.zip
"""
rm.runCommand("Save", target)
def save_selected_rois(rm, selected_rois, target):
"""Save selected ROIs in the RoiManager as zip to target path.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
selected_rois : array
ROIs in the RoiManager to save
target : string
the path in to store the ROIs. e.g. /my-images/resulting_rois_subset.zip
"""
rm.runCommand("Deselect")
rm.setSelectedIndexes(selected_rois)
rm.runCommand("save selected", target)
rm.runCommand("Deselect")
def enlarge_all_rois(amount_in_um, rm, pixel_size_in_um):
"""Enlarges all ROIs in the RoiManager by x scaled units.
Parameters
----------
amount_in_um : float
the value by which to enlarge in scaled units, e.g 3.5
rm : RoiManager
a reference of the IJ-RoiManager
pixel_size_in_um : float
the pixel size, e.g. 0.65 px/um
"""
amount_px = amount_in_um / pixel_size_in_um
all_rois = rm.getRoisAsArray()
rm.reset()
for roi in all_rois:
enlarged_roi = RoiEnlarger.enlarge(roi, amount_px)
rm.addRoi(enlarged_roi)
def select_positive_fibers(imp, channel, rm, min_intensity):
"""Select ROIs in ROIManager based on intensity in specific channel.
For all ROIs in the RoiManager, select ROIs based on intensity measurement in given channel of imp.
See https://imagej.nih.gov/ij/developer/api/ij/process/ImageStatistics.html
Parameters
----------
imp : ij.ImagePlus
the imp on which to measure
channel : integer
the channel on which to measure. starts at 1
rm : RoiManager
a reference of the IJ-RoiManager
min_intensity : integer
the selection criterion (here: intensity threshold)
Returns
-------
array
a selection of ROIs which passed the selection criterion (are above the threshold)
"""
imp.setC(channel)
all_rois = rm.getRoisAsArray()
selected_rois = []
for i, roi in enumerate(all_rois):
imp.setRoi(roi)
stats = imp.getStatistics()
if stats.mean > min_intensity:
selected_rois.append(i)
return selected_rois
def preset_results_column(rt, column, value):
"""Pre-set values in selected column from the ResultsTable.
Pre-set all rows in given column of the IJ-ResultsTable with desired value.
Parameters
----------
rt : ResultsTable
a reference of the IJ-ResultsTable
column : string
the desired column. will be created if it does not yet exist
value : string or float or integer
the value to be set
"""
for i in range(rt.size()):
rt.setValue(column, i, value)
rt.show("Results")
def add_results(rt, column, row, value):
"""Adds a value in desired rows of a given column.
Parameters
----------
rt : ResultsTable
a reference of the IJ-ResultsTable
column : string
the column in which to add the values
row : array
the row numbers in which too add the values.
value : string or float or integer
the value to be set
"""
for i in range( len( row ) ):
rt.setValue(column, row[i], value)
rt.show("Results")
def enhance_contrast(imp):
"""Use "Auto" Contrast & Brightness settings in each channel of imp.
Parameters
----------
imp : ij.ImagePlus
the imp on which to change C&B
"""
for channel in range(imp.getDimensions()[2]):
imp.setC(channel + 1) # IJ channels start at 1
IJ.run(imp, "Enhance Contrast", "saturated=0.35")
def renumber_rois(rm):
"""Rename all ROIs in the RoiManager according to their number.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
"""
number_of_rois = rm.getCount()
for roi in range(number_of_rois):
rm.rename(roi, str(roi + 1))
def setup_defined_ij(rm, rt):
"""Set up a clean and defined Fiji user environment.
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
rt : ResultsTable
a reference of the IJ-ResultsTable
"""
fix_ij_options()
rm.runCommand('reset')
rt.reset()
IJ.log("\\Clear")
execution_start_time = time.time()
setup_defined_ij(rm, rt)
print rt.size()
# open image using Bio-Formats
path_to_image = fix_ij_dirs(path_to_image)
raw = open_image_with_BF(path_to_image)
# get image info
raw_image_calibration = raw.getCalibration()
raw_image_title = fix_BF_czi_imagetitle(raw)
print("raw image title: ", str(raw_image_title))
# take care of paths and directories
output_dir = fix_ij_dirs(output_dir) + "/" + str(raw_image_title) + "/1_identify_fibers"
print("output_dir: ", str(output_dir))
if not os.path.exists( str(output_dir) ):
os.makedirs( str(output_dir) )
# update the log for the user
IJ.log( "Now working on " + str(raw_image_title) )
if raw_image_calibration.scaled() == False:
IJ.log("Your image is not spatially calibrated! Size measurements are only possible in [px].")
IJ.log( " -- settings used -- ")
IJ.log( "area = " + str(minAr) + "-" + str(maxAr) )
IJ.log( "perimeter = " + str(minPer) + "-" + str(maxPer) )
IJ.log( "circularity = " + str(minCir) + "-" + str(maxCir) )
IJ.log( "roundness = " + str(minRnd) + "-" + str(maxRnd) )
IJ.log( "solidity = " + str(minSol) + "-" + str(maxSol) )
IJ.log( "feret_ar = " + str(minFAR) + "-" + str(maxFAR) )
IJ.log( "min_feret = " + str(minMinFer) + "-" + str(maxMinFer) )
IJ.log( "ROI expansion [microns] = " + str(enlarge) )
IJ.log( "Membrane channel = " + str(membrane_channel) )
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)# imp, firstC, lastC, firstZ,
# lastZ, firstT, lastT
membrane = Duplicator().run(raw, membrane_channel, membrane_channel, 1, 1, 1, 1)
preprocess_membrane_channel(membrane)
imp_result = run_tm(membrane, 1, cellpose_dir.getPath(), PretrainedModel.CYTO2, 30.0, area_thresh=[minAr, maxAr], circularity_thresh=[minCir, maxCir],
perimeter_thresh=[minPer, maxPer],
# feret_thresh=[minMinFer, maxMinFer],
)
IJ.saveAs(imp_result, "Tiff", output_dir + "/" + raw_image_title + "_all_fibers_binary")
sys.exit()
# modify rois
rm.hide()
raw.hide()
enlarge_all_rois( enlarge, rm, raw_image_calibration.pixelWidth )
renumber_rois(rm)
save_all_rois( rm, output_dir + "/" + raw_image_title + "_all_fiber_rois.zip" )
# check for positive fibers
if fiber_channel > 0:
if min_fiber_intensity == 0:
min_fiber_intensity = get_threshold_from_method(raw, fiber_channel, "Mean")[0]
IJ.log( "automatic intensity threshold detection: True" )
IJ.log( "fiber intensity threshold: " + str(min_fiber_intensity) )
change_all_roi_color(rm, "blue")
positive_fibers = select_positive_fibers( raw, fiber_channel, rm, min_fiber_intensity )
change_subset_roi_color(rm, positive_fibers, "magenta")
save_selected_rois( rm, positive_fibers, output_dir + "/" + raw_image_title + "_mhc_positive_fiber_rois.zip")
# measure size & shape, save
IJ.run("Set Measurements...", "area perimeter shape feret's redirect=None decimal=4")
IJ.run("Clear Results", "")
measure_in_all_rois( raw, membrane_channel, rm )
rt = ResultsTable.getResultsTable("Results")
print rt.size()
IJ.log(" -- settings used -- ")
IJ.log("area = " + str(minAr) + "-" + str(maxAr))
IJ.log("perimeter = " + str(minPer) + "-" + str(maxPer))
IJ.log("circularity = " + str(minCir) + "-" + str(maxCir))
IJ.log("ROI expansion [microns] = " + str(enlarge_radius))
IJ.log("Membrane channel = " + str(membrane_channel))
IJ.log("MHC positive fiber channel = " + str(fiber_channel))
# IJ.log("sub-tiling = " + str(tiling_factor))
IJ.log(" -- settings used -- ")
print "saved the all_fibers_results.csv"
# dress up the original image, save a overlay-png, present original to the user
rm.show()
raw.show()
show_all_rois_on_image( rm, raw )
raw.setDisplayMode(IJ.COMPOSITE)
enhance_contrast( raw )
IJ.run("From ROI Manager", "") # ROIs -> overlays so they show up in the saved png
qc_duplicate = raw.duplicate()
IJ.saveAs(qc_duplicate, "PNG", output_dir + "/" + raw_image_title + "_all_fibers")
qc_duplicate.close()
wm.toFront( raw.getWindow() )
IJ.run("Remove Overlay", "")
raw.setDisplayMode(IJ.GRAYSCALE)
show_all_rois_on_image( rm, raw )
total_execution_time_min = (time.time() - execution_start_time) / 60.0
IJ.log("total time in minutes: " + str(total_execution_time_min))
IJ.log( "~~ all done ~~" )
IJ.selectWindow("Log")
IJ.saveAs("Text", str(output_dir + "/" + raw_image_title + "_all_fibers_Log"))
if close_raw == True:
raw.close()