Update LaVision processing

Author: constantinpape

flamingo_tools/extract_block_util.py

@@ -4,6 +4,7 @@
 import imageio.v3 as imageio
 import numpy as np
 import zarr
+from skimage.transform import rescale
 
 import flamingo_tools.s3_utils as s3_utils
 from flamingo_tools.file_utils import read_image_data
@@ -22,7 +23,8 @@ def extract_block(
     s3_credentials: Optional[str] = None,
     s3_bucket_name: Optional[str] = None,
     s3_service_endpoint: Optional[str] = None,
-):
+    scale_factor: Optional[Tuple[float, float, float]] = None,
+) -> None:
     """Extract block around coordinate from input data according to a given halo.
     Either from a local file or from an S3 bucket.
 
@@ -38,6 +40,7 @@ def extract_block(
         s3_bucket_name: S3 bucket name.
         s3_service_endpoint: S3 service endpoint.
         s3_credentials: File path to credentials for S3 bucket.
+        scale_factor: Optional factor for rescaling the extracted data.
     """
     coord_string = "-".join([str(int(round(c))).zfill(4) for c in coords])
 
@@ -93,9 +96,15 @@ def extract_block(
 
     data_ = read_image_data(input_path, input_key)
     data_roi = data_[roi]
+    if scale_factor is not None:
+        kwargs = {"preserve_range": True}
+        # Check if this is a segmentation.
+        if data_roi.dtype in (np.dtype("int32"), np.dtype("uint32"), np.dtype("int64"), np.dtype("uint64")):
+            kwargs.update({"order": 0, "anti_aliasing": False})
+        data_roi = rescale(data_roi, scale_factor, **kwargs).astype(data_roi.dtype)
 
     if tif:
         imageio.imwrite(output_file, data_roi, compression="zlib")
     else:
-        with zarr.open(output_file, mode="w") as f_out:
-            f_out.create_dataset(output_key, data=data_roi, compression="gzip")
+        f_out = zarr.open(output_file, mode="w")
+        f_out.create_dataset(output_key, data=data_roi, compression="gzip")

reproducibility/label_components/repro_label_components.py

@@ -40,6 +40,7 @@ def repro_label_components(
         component_list = dic["component_list"] if "component_list" in dic else default_component_list
 
         table_name = f"{cell_type.upper()}_{unet_version}"
+        # table_name = "PV_SGN_V2_DA"
         s3_path = os.path.join(f"{cochlea}", "tables", table_name, "default.tsv")
         tsv_path, fs = get_s3_path(s3_path, bucket_name=s3_bucket_name,
                                    service_endpoint=s3_service_endpoint, credential_file=s3_credentials)

reproducibility/templates_processing/REAMDE.md

@@ -13,7 +13,7 @@ For IHC segmentation run:
 - segment_unet_IHC_template.sbatch
 
 After this, run the following to add segmentation to MoBIE, create component labelings and upload to S3:
-- templates_transfer/mobie_segmentatio.sbatch
+- templates_transfer/mobie_segmentation_template.sbatch
 - templates_transfer/s3_seg_template.sh
 - repro_label_components.py
 - templates_transfer/s3_seg_template.sh

reproducibility/templates_processing/apply_unet_SGN_template.sbatch

@@ -4,12 +4,13 @@
 
 #SBATCH -p grete:shared         # the partition
 #SBATCH -G A100:1               # For requesting 1 A100 GPU.
-#SBATCH -c 1
-#SBATCH --mem 24G
+#SBATCH -c 4
+#SBATCH --mem 32G
 #SBATCH -a 0-9
 
 source ~/.bashrc
-micromamba activate micro-sam_gpu
+# micromamba activate micro-sam_gpu
+micromamba activate sam
 
 # Print out some info.
 echo "Submitting job with sbatch from directory: ${SLURM_SUBMIT_DIR}"
@@ -19,7 +20,8 @@ echo "Current node: ${SLURM_NODELIST}"
 
 # Run the script
 
-SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+# SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+SCRIPT_REPO=/user/pape41/u12086/Work/my_projects/flamingo-tools
 cd "$SCRIPT_REPO"/flamingo_tools/segmentation/ || exit
 
 export SCRIPT_DIR=$SCRIPT_REPO/scripts
@@ -37,9 +39,16 @@ export INPUT=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/"$COC
 
 export OUTPUT_FOLDER=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/"$COCHLEA"/"$SEG_NAME"
 
-export MODEL=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/trained_models/SGN/v2_cochlea_distance_unet_SGN_supervised_2025-05-27
+# The default v2 model
+# export MODEL=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/trained_models/SGN/v2_cochlea_distance_unet_SGN_supervised_2025-05-27
+
+# Domain adapted model for MLR99L
+export MODEL=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/trained_models/SGN/v2_domain_adaptation_mlr99l/best.pt
+
 export PREDICTION_INSTANCES=10
-export INPUT_KEY="setup$STAIN_CHANNEL/timepoint0/s0"
+
+# export INPUT_KEY="setup$STAIN_CHANNEL/timepoint0/s0"
+export INPUT_KEY="s0"
 
 echo "Input directory: ${INPUT}"
 echo "Output directory: ${OUTPUT_FOLDER}"

reproducibility/templates_processing/mean_std_SGN_template.sbatch

@@ -8,11 +8,13 @@
 #SBATCH --mem 128G
 
 source ~/.bashrc
-micromamba activate flamingo13
+# micromamba activate flamingo13
+micromamba activate sam
 
 # Run the script
 
-SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+# SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+SCRIPT_REPO=/user/pape41/u12086/Work/my_projects/flamingo-tools
 cd "$SCRIPT_REPO"/flamingo_tools/segmentation/ || exit
 
 export SCRIPT_DIR=$SCRIPT_REPO/scripts
@@ -26,15 +28,13 @@ STAIN_CHANNEL=$3
 # segmentation name, as it appears in MoBIE, e.g. SGN_v2 or Calb1_SGN_v2
 SEG_NAME=$4
 
-# TODO update model to
-# /mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/trained_models/SGN/v2_domain_adaptation_mlr99l
-# and download the PV channel again
-
 export INPUT=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/"$COCHLEA"/"$DATA"
 
 export OUTPUT_FOLDER=/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/"$COCHLEA"/"$SEG_NAME"
 export SEG_CLASS="sgn"
-export INPUT_KEY="setup$STAIN_CHANNEL/timepoint0/s0"
+
+# export INPUT_KEY="setup$STAIN_CHANNEL/timepoint0/s0"
+export INPUT_KEY="s0"
 
 if ! [[ -f $OUTPUT_FOLDER ]] ; then
 	mkdir -p "$OUTPUT_FOLDER"

reproducibility/templates_processing/segment_unet_SGN_template.sbatch

@@ -8,18 +8,17 @@
 #SBATCH --mem 400G
 
 source ~/.bashrc
-micromamba activate micro-sam_gpu
+# micromamba activate micro-sam_gpu
+micromamba activate sam
 
 # Print out some info.
 echo "Submitting job with sbatch from directory: ${SLURM_SUBMIT_DIR}"
 echo "Home directory: ${HOME}"
 echo "Working directory: $PWD"
 echo "Current node: ${SLURM_NODELIST}"
 
-# Run the script
-#python myprogram.py $SLURM_ARRAY_TASK_ID
-
-SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+# SCRIPT_REPO=/user/schilling40/u15000/flamingo-tools
+SCRIPT_REPO=/user/pape41/u12086/Work/my_projects/flamingo-tools
 cd "$SCRIPT_REPO"/flamingo_tools/segmentation/ || exit
 
 # name of cochlea, as it appears in MoBIE and the NHR

scripts/figures/segmentation_comparison.py

@@ -29,6 +29,7 @@ def _eval_seg(seg, eval_path):
 
 
 def sgn_comparison():
+    scale = (0.38,) * 2
     z = 10
 
     cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet"
@@ -57,22 +58,29 @@ def sgn_comparison():
             eval_im[seg_name] = eva
 
         v = napari.Viewer()
-        v.add_image(image)
+        v.add_image(image, scale=scale)
         for seg_name, bd in boundaries.items():
-            v.add_labels(bd, name=seg_name, colormap={1: "cyan"})
-            v.add_labels(eval_im[seg_name], name=f"{seg_name}_eval", colormap={1: "green", 2: "red"})
+            v.add_labels(bd, name=seg_name, colormap={1: "cyan"}, scale=scale)
+            v.add_labels(eval_im[seg_name], name=f"{seg_name}_eval", colormap={1: "green", 2: "red"}, scale=scale)
+
+        v.scale_bar.visible = True
+        v.scale_bar.unit = "μm"
+        v.scale_bar.font_size = 16
         v.title = Path(path).stem
+
         napari.run()
 
 
 def ihc_comparison():
+    scale = (0.38,) * 2
     z = 10
 
     cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet"
     val_sgn_dir = f"{cochlea_dir}/predictions/val_ihc"
     image_dir = f"{cochlea_dir}/AnnotatedImageCrops/F1ValidationIHCs"
 
     image_paths = sorted(glob(os.path.join(image_dir, "*.tif")))
+    scale = (0.38,) * 2
 
     for path in image_paths:
         image = imageio.imread(path)[z]
@@ -96,11 +104,16 @@ def ihc_comparison():
             eval_im[seg_name] = eva
 
         v = napari.Viewer()
-        v.add_image(image)
+        v.add_image(image, scale=scale)
         for seg_name, bd in boundaries.items():
-            v.add_labels(bd, name=seg_name, colormap={1: "cyan"})
-            v.add_labels(eval_im[seg_name], name=f"{seg_name}_eval", colormap={1: "green", 2: "red"})
+            v.add_labels(bd, name=seg_name, colormap={1: "cyan"}, scale=scale)
+            v.add_labels(eval_im[seg_name], name=f"{seg_name}_eval", colormap={1: "green", 2: "red"}, scale=scale)
+
+        v.scale_bar.visible = True
+        v.scale_bar.unit = "μm"
+        v.scale_bar.font_size = 16
         v.title = Path(path).stem
+
         napari.run()
 
 

scripts/figures/supp_fig2.py

@@ -0,0 +1,70 @@
+import json
+import os
+
+from glob import glob
+from pathlib import Path
+
+import numpy as np
+
+
+# FIXME something is off with cellpose-sam runtimes
+def runtimes_sgn():
+    for_comparison = ["distance_unet", "micro-sam", "cellpose3", "cellpose-sam", "stardist"]
+
+    cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet"
+    val_sgn_dir = f"{cochlea_dir}/predictions/val_sgn"
+    image_dir = f"{cochlea_dir}/AnnotatedImageCrops/F1ValidationSGNs/for_consensus_annotation"
+
+    image_paths = sorted(glob(os.path.join(image_dir, "*.tif")))
+
+    runtimes = {name: [] for name in for_comparison}
+
+    for path in image_paths:
+        eval_fname = Path(path).stem + "_dic.json"
+        for seg_name in for_comparison:
+            eval_path = os.path.join(val_sgn_dir, seg_name, eval_fname)
+            with open(eval_path, "r") as f:
+                result = json.load(f)
+            rt = result["time"]
+            runtimes[seg_name].append(rt)
+
+    for name, rts in runtimes.items():
+        print(name, ":", np.mean(rts), "+-", np.std(rts))
+
+
+def runtimes_ihc():
+    for_comparison = ["distance_unet_v3", "micro-sam", "cellpose3", "cellpose-sam"]
+
+    cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet"
+    val_sgn_dir = f"{cochlea_dir}/predictions/val_ihc"
+    image_dir = f"{cochlea_dir}/AnnotatedImageCrops/F1ValidationIHCs"
+
+    image_paths = sorted(glob(os.path.join(image_dir, "*.tif")))
+
+    runtimes = {name: [] for name in for_comparison}
+
+    for path in image_paths:
+        eval_fname = Path(path).stem + "_dic.json"
+        for seg_name in for_comparison:
+            eval_path = os.path.join(val_sgn_dir, seg_name, eval_fname)
+            if not os.path.exists(eval_path):
+                continue
+            with open(eval_path, "r") as f:
+                result = json.load(f)
+            rt = result["time"]
+            runtimes[seg_name].append(rt)
+
+    for name, rts in runtimes.items():
+        print(name, ":", np.mean(rts), "+-", np.std(rts))
+
+
+def main():
+    print("SGNs:")
+    runtimes_sgn()
+    print()
+    print("IHCs:")
+    runtimes_ihc()
+
+
+if __name__ == "__main__":
+    main()

scripts/la-vision/analyze_la_vision.py

@@ -0,0 +1 @@
+import

scripts/la-vision/import_marmoset_lsm.py

@@ -0,0 +1,36 @@
+import os
+
+import imageio.v3 as imageio
+from mobie import add_image
+
+INPUT_ROOT = "/mnt/ceph-hdd/cold/nim00007/cochlea-lightsheet/keppeler-et-al/marmoset"
+MOBIE_ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/mobie_project/cochlea-lightsheet"
+DS_NAME = "LaVision-Mar05"
+RESOLUTION = (3.0, 1.887779, 1.887779)
+
+
+# Marmoset_cochlea_05_LSFM_ch1_raw.tif
+def add_marmoset_05():
+    channel_names = ("PV", "7-AAD", "MYO")
+
+    scale_factors = 4 * [[2, 2, 2]]
+    chunks = (96, 96, 96)
+
+    for channel_id, channel_name in enumerate(channel_names, 1):
+        input_path = os.path.join(INPUT_ROOT, f"Marmoset_cochlea_05_LSFM_ch{channel_id}_raw.tif")
+        print("Load image data ...")
+        input_data = imageio.imread(input_path)
+        print(input_data.shape)
+        add_image(
+            input_path=input_data, input_key=None, root=MOBIE_ROOT,
+            dataset_name=DS_NAME, image_name=channel_name, resolution=RESOLUTION,
+            scale_factors=scale_factors, chunks=chunks, unit="micrometer", use_memmap=False,
+        )
+
+
+def main():
+    add_marmoset_05()
+
+
+if __name__ == "__main__":
+    main()