Newer
Older
# @ 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
# from loci.plugins import BF
# from loci.plugins.in import ImporterOptions
from ch.epfl.biop.ij2command import Labels2CompositeRois
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 ────────────────────────────────────────────────────────────────
# 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):
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 + "/"
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.
image_title = os.path.basename(imp.getShortTitle())
# remove the ".czi" extension
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.
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.
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],
Function to run TrackMate on open data, applying filters to spots.
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 circularity 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
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
# 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)
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
# 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.
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
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".
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()
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")
Parameters
----------
rm : RoiManager
a reference of the IJ-RoiManager
the imp on which to show the ROIs
"""
imp.show()
"""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
----------
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)
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
"""
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
"""
rt.setValue(column, row[i], value)
rt.show("Results")
def enhance_contrast(imp):
"""Use "Auto" Contrast & Brightness settings in each channel of imp.
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))
"""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()
# ─── Main Code ────────────────────────────────────────────────────────────────
if __name__ == "__main__":
execution_start_time = time.time()
IJ.log("\\Clear")
misc.timed_log("Script starting")
setup_defined_ij(rm, rt)
file_list = pathtools.listdir_matching(
src_dir.getPath(), filename_filter, fullpath=True
)
out_dir_info = pathtools.parse_path(output_dir)
for index, file in enumerate(file_list):
# open image using Bio-Formats
file_info = pathtools.parse_path(file)
misc.progressbar(index + 1, len(file_list), 1, "Opening : ")
raw = bf.import_image(file_info["full"])[0]
# 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 = os.path.join(
out_dir_info["full"], 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
misc.timed_log("Now working on " + str(raw_image_title))
if raw_image_calibration.scaled() is False:
IJ.log(
"Your image is not spatially calibrated! Size measurements are only possible in [px]."
)
# Only print it once since we'll use the same settings everytime
if index == 0:
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 -- ")
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
# 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_bgd_corrected = do_background_correction(membrane)
IJ.run("Conversions...", "scale")
IJ.run(imp_bgd_corrected, "16-bit", "")
imp_result = run_tm(
imp_bgd_corrected,
1,
cellpose_dir.getPath(),
PretrainedModel.CYTO2,
30.0,
area_thresh=[minAr, maxAr],
circularity_thresh=[minCir, maxCir],
perimeter_thresh=[minPer, maxPer],
)
IJ.saveAs(
imp_result,
"Tiff",
os.path.join(output_dir, raw_image_title + "_all_fibers_binary"),
)
command.run(Labels2CompositeRois, True, "rm", rm, "imp", imp_result).get()
enlarge_all_rois(enlarge_radius, rm, raw_image_calibration.pixelWidth)
renumber_rois(rm)
save_all_rois(
rm, os.path.join(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,
os.path.join(
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())
if fiber_channel > 0:
# print(rt.size())
preset_results_column(rt, "MHC Positive Fibers (magenta)", "NO")
# print(rt.size())
add_results(rt, "MHC Positive Fibers (magenta)", positive_fibers, "YES")
# print(rt.size())
rt.save(os.path.join(output_dir, raw_image_title + "_all_fibers_results.csv"))
# 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)
IJ.selectWindow("Log")
IJ.saveAs("Text", str(output_dir + "/" + raw_image_title + "_all_fibers_Log"))
membrane.close()
imp_bgd_corrected.close()
imp_result.close()
if close_raw == True:
raw.close()
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 ~~")