feat: begin setup of MedSAM2 benchmarking

.gitattributes

@@ -1,2 +1,2 @@
 # SCM syntax highlighting & preventing 3-way merges
-pixi.lock merge=binary linguist-language=YAML linguist-generated=true
+pixi.lock merge=binary linguist-language=YAML linguist-generated=true -diff

README.md

@@ -1,4 +1,4 @@
 Steps to add a model:
 
-1. Create a contract.py file in the models/<model_name> directory. That will get the input data ready to run the model.
+1. Create a prepare.py file in the models/<model_name> directory. That will get the input data ready to run the model.
 2. Create a .def file with the entry point for the model

docs/devnotes.md

@@ -1,47 +1,25 @@
 # Developer Notes
 
-## Purpose of This Section
-
-This section is for documenting technical decisions, challenges, and solutions encountered during your project. These notes are valuable for:
-
-- Future you (who will forget why certain decisions were made)
-- Collaborators who join the project later
-- People coming from your publication who want to reproduce your work
-- Anyone who might want to extend your research
-
-## What to Document
-
-### Design Decisions
-
+## Design Decisions  
 Document important decisions about your project's architecture, algorithms, or methodologies:
 
-``` markdown
-## Choice of RNA-Seq Analysis Pipeline
+### Choice of RECIST longest diameter calculation algorithm
+[2026-08-10] Katy S. compared the different algorithms developed by MedSAM2 and BHKLab teams to calculate the longest axial diameter of a 3D segmentation. MedSAM2's method in [`compute_recist_line`](../src/ab_testing/utils/conver_nifti_to_recist_npz.py) under function  utilizes OpenCV contours. BHKLab's method (by Kaitlyn K.) in [`get_recist_pts`](../src/ab_testing/utils/annotations.py) measures from the centroid of the slice based on a `transformation of the origin * half of the major axis length`. Based on some testing with the NSCLC-Radiomics dataset, the MedSAM2 method was visually assessed to be more accurate. [Exploration script](../sandbox/recist_comparison.py)
 
-[2025-04-25] We chose the kallisto over STAR pipeline for the following reasons:
-    1. The CCLE dataset is very large, and kallisto is faster for quantifying large datasets
-    2. GDSC used kallisto, so we can compare our results with theirs
-```
+### Multiple lesion handling  
+[2026-08-10] The nifti to recist npz setup right now expects an nnUnet directory structure for images where the label files can contain multiple lesions with each lesion having a different label number. I (Katy) typically run med-imagetools with the SEPARATE ROI strategy, but could use MERGE so that the data is set up this way. 
 
-### Technical Challenges
 
+## Technical Challenges  
 Record significant problems you encountered and how you solved them
 
-``` markdown
-## Sample Name Format Issue
 
-[2025-04-25] We encountered a problem with sample name formats between the CCLE and GDSC datasets.
-    The CCLE dataset uses "BRCA-XX-XXXX" format, while the GDSC dataset uses "BRCA-XX-XXXX-XX".
-    We had to write a script to remove the last two characters from the sample names in the GDSC dataset.
-```
-
-### Dependencies and Environment
+## Dependencies and Environment
 
 Document specific version requirements or compatibility issues:
 
 ``` markdown
-## Critical Version Dependencies
-
+### Critical Version Dependencies  
 [2025-04-25] SimpleITK 2.4.1 introduced a bug that flips images, so we froze version 2.4.0
 ```
 

pyproject.toml

@@ -0,0 +1,42 @@
+[project]
+authors = [{name = "Joshua Siraj", email = "joshua.siraj@ryerson.ca"}, {name = "Katy Scott", email = "bhklab.katyscott@gmail.com"}]
+dependencies = [
+    "simpleitk>=2.5.6,<3",
+    "pydantic>=2.13.4,<3",
+    "damply>=0.26.1,<0.27",
+    "med-imagetools>=2.16.0,<3", 
+    "connected-components-3d>=4.0.0,<5", 
+    "opencv-python>=5.0.0.93,<6", 
+    "nibabel>=5.4.2,<6", 
+    "readii",
+    "click>=8.4.2,<9",
+]
+name = "ab-testing"
+requires-python = ">= 3.11"
+version = "0.1.0"
+
+[build-system]
+build-backend = "hatchling.build"
+requires = ["hatchling"]
+
+[tool.pixi.workspace]
+channels = ["conda-forge"]
+platforms = ["linux-64", "osx-64", "osx-arm64"]
+
+[tool.pixi.pypi-dependencies]
+ab_testing = { path = ".", editable = true }
+
+[tool.pixi.activation]
+# convenient variables which can be used in scripts
+env.CONFIG = "${PIXI_PROJECT_ROOT}/config"
+env.METADATA = "${PIXI_PROJECT_ROOT}/metadata"
+env.LOGS = "${PIXI_PROJECT_ROOT}/logs"
+env.RAWDATA = "${PIXI_PROJECT_ROOT}/data/rawdata"
+env.PROCDATA = "${PIXI_PROJECT_ROOT}/data/procdata"
+env.RESULTS = "${PIXI_PROJECT_ROOT}/data/results"
+env.SCRIPTS = "${PIXI_PROJECT_ROOT}/workflow/scripts"
+
+[tool.pixi.tasks]
+
+[tool.pixi.environments.qc.dependencies]
+ruff = ">=0.16.2,<0.17"

ruff.toml

@@ -73,8 +73,8 @@ ignore = [
     "D100", # Allows modules without docstrings
     # https://docs.astral.sh/ruff/formatter/#conflicting-lint-rules
     "COM812", # Avoids conflict with Ruff formatter trailing comma rules
-    "D206"
-    # Allows docstrings to be formatted flexibly
+    "D206", # Allows docstrings to be formatted flexibly
+    "N813", # Allows for Camelcase to lowercase imports
 ]
 
 # https://docs.astral.sh/ruff/settings/#isort

sandbox/recist_comparison.py

@@ -0,0 +1,64 @@
+# Script to compare the calculation of a RECIST longest diameter from a 3D tumour segmentation mask.
+# Compared method from MedSAM2 implementation and BHKLab implementation by Kaitlyn Kobayashi. 
+
+import cv2
+import matplotlib as mpl
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
+import numpy as np
+from damply import dirs
+from readii.image_processing import displayCTSegOverlay
+
+from ab_testing.utils.annotations import get_line_from_recist, get_recist_pts
+from ab_testing.utils.conver_nifti_to_recist_npz import (
+                 compute_recist_line,
+                 generate_recist,
+                 read_case,
+)
+
+cool_cmap = mpl.colormaps['cool']
+binary_magenta = mcolors.ListedColormap([(0,0,0,0.1), cool_cmap(400)])
+
+sample_scan = "LUNG1-012_0011/CT_55636694"
+sample_mask = "LUNG1-012_0011/RTSTRUCT_38424583"
+
+sample_id = sample_scan.split("/", maxsplit=1)[0]
+
+scan_path = dirs.RAWDATA / "TCIA_NSCLC-Radiomics/images/mit_NSCLC-Radiomics" / f"{sample_scan}" / "CT.nii.gz"
+mask_path = dirs.RAWDATA / "TCIA_NSCLC-Radiomics/images/mit_NSCLC-Radiomics" / f"{sample_mask}" / "GTV.nii.gz"
+
+scan, mask, spacing, direction, origin, reader = read_case(scan_path, mask_path)
+
+key_slice = int(np.argmax(np.sum(mask, axis=(1, 2))))
+
+r_fig, r_ax = plt.subplots(1,3, figsize=(21,8.5), layout='tight')
+plt.suptitle(f'RECIST line comparison for NSCLC-Radiomics {sample_id}')
+
+# MedSAM2 method
+p1, p2 = compute_recist_line(mask[key_slice])
+m_recist_line, short_ids = generate_recist(mask, spacing)
+
+displayCTSegOverlay(scan, mask, sliceIdx=key_slice, dispMin=-1350, dispMax=150, ax=r_ax[0], alpha=0.3)
+r_ax[0].imshow(m_recist_line[key_slice], cmap=binary_magenta, alpha=0.7)
+r_ax[0].set_title(f"MedSAM2 \n Coordinates: p1={p1}, p2={p2}")
+
+# Kaitlyn's method
+recist_pts = get_recist_pts(mask[key_slice])
+k_recist_line = get_line_from_recist(recist_pts, key_slice, mask.shape)
+f_k_pts = f"Coordinates: p1=[{recist_pts[2]:.2f} {recist_pts[3]:.2f}], p2=[{recist_pts[0]:.2f} {recist_pts[1]:.2f}]"
+
+displayCTSegOverlay(scan, mask, sliceIdx=key_slice, dispMin=-1350, dispMax=150, ax=r_ax[1], alpha=0.3)
+r_ax[1].imshow(k_recist_line[key_slice], cmap=binary_magenta, alpha=0.7)
+r_ax[1].set_title(f"BHKLab \n {f_k_pts}")
+
+# Draw line like MedSAM method with Kaitlyn's coordinates
+k_thick_line = np.zeros_like(mask, dtype=np.uint16)
+cv2.line(k_thick_line[key_slice], (int(recist_pts[0]), int(recist_pts[1])), (int(recist_pts[2]), int(recist_pts[3])),
+                 color=1, thickness=2)
+f_kt_pts = f"Coordinates: p1=[{int(recist_pts[2])} {int(recist_pts[3])}], p2=[{int(recist_pts[0])} {int(recist_pts[1])}]"
+
+displayCTSegOverlay(scan, mask, sliceIdx=key_slice, dispMin=-1350, dispMax=150, ax=r_ax[2], alpha=0.3)
+r_ax[2].imshow(k_thick_line[key_slice], cmap=binary_magenta, alpha=0.7)
+r_ax[2].set_title(f"BHKLab drawn like MedSAM2 \n {f_kt_pts}")
+
+r_fig.savefig(dirs.RESULTS / f"recist_comparison_{sample_id}.png")

contract.py → src/ab_testing/contract.py

@@ -1,6 +1,7 @@
 from pydantic import BaseModel
-from input import SUPPORTED_INPUT_TYPES
-from output import SUPPORTED_OUTPUT_TYPES
+
+from ab_testing.input import SUPPORTED_INPUT_TYPES
+from ab_testing.output import SUPPORTED_OUTPUT_TYPES
 
 
 class ModelContract(BaseModel):

input.py → src/ab_testing/input.py

@@ -1,4 +1,3 @@
-from tarfile import SUPPORTED_TYPES
 from pydantic import BaseModel
 
 

src/ab_testing/models/medsam2/prepare.py

@@ -0,0 +1,141 @@
+import logging
+from pathlib import Path
+
+import click
+from damply import dirs
+
+from ab_testing.utils.conver_nifti_to_recist_npz import (
+    MIN_LABEL_SLICES,
+    MIN_RECIST_MM,
+    run_dataset,
+)
+
+logfile = dirs.LOGS / "models" / "prepare_medsam2.log"
+logfile.parent.mkdir(parents=True, exist_ok=True)
+
+logging.basicConfig(
+	level=logging.INFO, 
+    format='%(asctime)s - %(levelname)s - %(message)s',
+    filename=logfile
+)
+
+logger = logging.getLogger(__name__)
+
+ANATOMICAL_BUCKETS = [
+    "lung",
+    "soft_tissue",   # includes heart
+    "bone",
+    "brain",
+]
+
+@click.command(no_args_is_help=True)
+@click.argument(
+    "dataset",
+    type=str,
+)
+@click.argument(
+    "mit_directory",
+    type=click.Path(
+        file_okay=False, dir_okay=True, writable=False, path_type=Path, resolve_path=True, exists=True
+    ),
+)
+@click.argument(
+    "output_directory",
+    type=click.Path(
+        file_okay=False, dir_okay=True, writable=True, path_type=Path, resolve_path=True
+    ),
+)
+@click.option(
+    "--anat-window",
+    type=click.Choice(ANATOMICAL_BUCKETS),
+    required=True,
+    help=f"Anatomical windowing to apply for this dataset. Must be one of {ANATOMICAL_BUCKETS}."
+)
+@click.option(
+    "--min-slices",
+    type=int,
+    default=MIN_LABEL_SLICES,
+    show_default=True,
+    help="Minimum labelled slice count required for npz conversion."
+)
+@click.option(
+    "--min-recist-mm",
+    type=float,
+    default=MIN_RECIST_MM,
+    show_default=True,
+    help="Minimum RECIST diameter measurement required for npz conversion in millimetres (mm)."
+)
+@click.option(
+    "--workers",
+    type=int,
+    default=24,
+    help="Number of parallel processes to use."
+)
+@click.option(
+    "--overwrite",
+    is_flag=True,
+    help="Whether to overwrite existing files."
+)
+def prepare_data(
+    dataset: str, 
+    mit_directory: str | Path, 
+    output_directory: str | Path,
+    anat_window: str,
+    *,
+    min_slices: int = MIN_LABEL_SLICES,
+    min_recist_mm: float = MIN_RECIST_MM,
+    workers: int = 24,
+    overwrite: bool = False,
+    ) -> int:
+    """Process a med-imagetools nifti dataset by pair of image-masks and save into npzs.
+
+    Converted npz contains:
+        * imgs: The image, windowed based on the anatomical window specified,
+        * gts: The corresponding 3D segmentation mask, 
+        * spacing: Original image/mask spacing, from SimpleITK
+        * direction: Original image/mask direction, from SimpleITK
+        * origin: Original image/mask origin, from SimpleTIK 
+        * reader: What library was used to load in the image/mask, will be sitk or nibabel_orthofix
+        * recist: Reverse-engineered RECIST longest diameter annotation derived from gts
+
+    Based on min_slices and min_recist_mm, samples will be pruned and not saved as npz if they fall below the threshold.
+
+    Args:
+        dataset: str, dataset name
+        mit_directory: str | Path, pathway to the image and mask directory
+        output_directory: str | Path, output directory npz files will be saved to
+        anat_window: anatomical windowing to apply for this dataset
+        min_slices: int, minimum number of axial slices with label for lesion inclusion (default 5)
+        min_recist_mm: float, RECIST diameter threshold in millimetres for lesion inclusion (default 10)
+        workers: int, number of parallel processes to use (default 24)
+        overwrite: bool, whether to reconvert every case (default False, skip existing npz)
+
+    Returns:
+        0 if successful, 1 if errors occurred.
+    """
+    written, skipped, filtered, errors, dropped_lesions, _z_kept, _z_total = run_dataset(
+        ds = dataset,
+        root = Path(mit_directory),
+        out_root = Path(output_directory),
+        workers = workers,
+        with_recist = True,
+        overwrite = overwrite,
+        min_slices = min_slices,
+        min_recist_mm = min_recist_mm,
+        prune_filtered = True,
+        anat_window = anat_window,
+        pair_builder = 'mit',
+    )
+
+    logger.info(f"\nTOTAL written={written} skipped={skipped} "
+                f"filtered={filtered} errors={errors} "
+                f"small_lesions_dropped={dropped_lesions}")
+
+    return 1 if errors else 0
+
+
+
+
+
+if __name__ == "__main__":
+    prepare_data()