1_identify_fibers.py

@@ -67,6 +67,7 @@ from ij.plugin import Duplicator, ImageCalculator, RoiEnlarger
 from imcflibs import pathtools
 from imcflibs.imagej import bioformats as bf
 from imcflibs.imagej import misc
+from loci.formats import ImageReader, MetadataTools
 
 # ─── Functions ────────────────────────────────────────────────────────────────
 
@@ -625,6 +626,62 @@ def setup_defined_ij(rm, rt):
     IJ.log("\\Clear")
 
 
+def get_series_info_from_ome_metadata(path_to_file, skip_labels=False):
+    """Get the Bio-Formates series count from a file on disk.
+
+    Useful to access a specific image in a container format like .czi, .nd2, .lif...
+
+    Parameters
+    ----------
+    path_to_file : str
+        The full path to the image file.
+
+    Returns
+    -------
+    int
+        The number of Bio-Formats series detected in the image file metadata.
+    """
+
+    if not skip_labels:
+        reader = ImageReader()
+        reader.setFlattenedResolutions(False)
+        ome_meta = MetadataTools.createOMEXMLMetadata()
+        reader.setMetadataStore(ome_meta)
+        reader.setId(path_to_file)
+        series_count = reader.getSeriesCount()
+
+        reader.close()
+        return series_count, range(series_count)
+
+    else:
+        reader = ImageReader()
+        # reader.setFlattenedResolutions(True)
+        ome_meta = MetadataTools.createOMEXMLMetadata()
+        reader.setMetadataStore(ome_meta)
+        reader.setId(path_to_file)
+        series_count = reader.getSeriesCount()
+
+        series_ids = []
+        series_names = []
+        x = 0
+        y = 0
+        for i in range(series_count):
+            reader.setSeries(i)
+
+            if reader.getSizeX() > x and reader.getSizeY() > y:
+                name = ome_meta.getImageName(i)
+
+                if name not in ["label image", "macro image"]:
+                    series_ids.append(i)
+                    series_names.append(name)
+
+            x = reader.getSizeX()
+            y = reader.getSizeY()
+
+        print(series_names)
+        return len(series_ids), series_ids
+
+
 # ─── Main Code ────────────────────────────────────────────────────────────────
 
 if __name__ == "__main__":
@@ -643,141 +700,201 @@ if __name__ == "__main__":
         # 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"
+        _, series_index = get_series_info_from_ome_metadata(
+            file_info["full"], skip_labels=True
         )
-        print("output_dir: ", str(output_dir))
+        print(series_index)
+        for list_index, serie_index in enumerate(series_index):
+            misc.progressbar(list_index + 1, len(series_index), 2, "Opening serie : ")
+            raw = bf.import_image(file_info["full"], series_number=serie_index)[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 settings once for the first file and first series
+            if index == 0 and serie_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 -- ")
+
+            # 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
+            )
 
-        if not os.path.exists(str(output_dir)):
-            os.makedirs(str(output_dir))
+            if (membrane.getWidth() * membrane.getHeight()) > 100000:
+                misc.timed_log("Image is too large, resizing to speed up processing")
+                membrane = membrane.resize(
+                    membrane.getWidth() / 2,
+                    membrane.getHeight() / 2,
+                    "none",
+                )
+            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],
+            )
 
-        # 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]."
+            imp_result = imp_result.resize(
+                raw.getWidth(),
+                raw.getHeight(),
+                "none",
             )
-        # 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 -- ")
-
-        # 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()
+            IJ.saveAs(
+                imp_result,
+                "Tiff",
+                os.path.join(
+                    output_dir,
+                    raw_image_title + "_" + str(serie_index) + "_all_fibers_binary",
+                ),
+            )
 
-        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")
-        )
+            command.run(Labels2CompositeRois, True, "rm", rm, "imp", imp_result).get()
 
-        # 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(
+            enlarge_all_rois(enlarge_radius, rm, raw_image_calibration.pixelWidth)
+            renumber_rois(rm)
+            save_all_rois(
                 rm,
-                positive_fibers,
                 os.path.join(
-                    output_dir, raw_image_title + "_mhc_positive_fiber_rois.zip"
+                    output_dir,
+                    raw_image_title + "_" + str(serie_index) + "_all_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")
+            # 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
+                        + "_"
+                        + str(serie_index)
+                        + "_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)
 
-        # print(rt.size())
+            rt = ResultsTable.getResultsTable("Results")
 
-        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()
+            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
+                    + "_"
+                    + str(serie_index)
+                    + "_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
+                + "_"
+                + str(serie_index)
+                + "_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
+                    + "_"
+                    + str(serie_index)
+                    + "_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))