Debug

.gitignore

@@ -169,4 +169,5 @@ dicom_select
 examples
 slicerio_data
 *.nrrd
-nnUNet_results
+nnUNet_results
+CLAUDE*.md

TPTBox/core/dicom/dicom_extract.py

@@ -31,7 +31,7 @@
 
 import string
 
-from TPTBox.core.dicom.dicom2nii_utils import get_json_from_dicom, load_json, save_json, test_name_conflict
+from TPTBox.core.dicom.dicom2nii_utils import get_json_from_dicom, load_json, save_json, secure_save_json, test_name_conflict
 
 logger = Print_Logger()
 
@@ -552,9 +552,108 @@ def _from_dicom_to_nii(
 
     if add_grid:
         _add_grid_info_to_json(nii_path, json_file_name)
+        # Multi-echo Philips DIXON (magnitude/phase) arrives as a 4-D NIfTI.
+        # Split it into per-echo 3-D files with `-eco<i>` appended to `part`.
+        if json_bids.get("part") in ("magnitude", "phase"):
+            _split_multi_echo_dixon(Path(nii_path), Path(json_file_name), dcm_data_l)
     return nii_path if add_grid else None
 
 
+def _split_multi_echo_dixon(nii_path: Path, json_path: Path, dcm_data_l) -> list[Path] | None:
+    """Split a 4-D multi-echo DIXON (magnitude/phase) NIfTI into per-echo 3-D files.
+
+    A no-op unless the NIfTI on disk is 4-D. Per-echo outputs are named
+    ``..._part-<part>-eco<i>_<mod>.nii.gz`` for ``i = 0..N-1`` (ascending
+    echo). Per-echo TEs are taken from the source DICOM headers (grouped by
+    ``EchoNumbers`` in ascending order) and written into each per-echo
+    sidecar JSON as ``EchoTime`` (with ``EchoNumbers`` set to the 1-based
+    echo index). The 4-D file and its sidecar are removed on success.
+
+    Args:
+        nii_path: Path to the just-written 4-D NIfTI.
+        json_path: Path to its sidecar JSON.
+        dcm_data_l: The source DICOM datasets that produced *nii_path*.
+
+    Returns:
+        List of per-echo NIfTI paths on success, or ``None`` when the input
+        was not 4-D / the split could not be performed.
+    """
+    if not Path(nii_path).exists():
+        return None
+    nii = NII.load(nii_path, False)
+    if nii.get_num_dims() != 4:
+        return None
+    n_echo = nii.shape[-1]
+
+    # Ascending (EchoNumbers -> EchoTime) from the source DICOMs. Classic
+    # Philips single-frame series carry both tags on every file.
+    te_by_echo: dict[int, float] = {}
+    if isinstance(dcm_data_l, list):
+        for d in dcm_data_l:
+            try:
+                en = int(getattr(d, "EchoNumbers", 0) or 0)
+                te = float(getattr(d, "EchoTime", 0.0) or 0.0)
+            except Exception:  # noqa: BLE001
+                continue
+            if en > 0:
+                te_by_echo.setdefault(en, te)
+    tes: list[float | None] = [te_by_echo[k] for k in sorted(te_by_echo)] if te_by_echo else []
+    if len(tes) != n_echo:
+        if len(tes) != 0:
+            logger.on_warning(
+                f"Multi-echo DIXON split: {n_echo} volumes but {len(tes)} unique EchoNumbers "
+                f"in DICOM; falling back to axis order without per-echo TE."
+            )
+        tes = [None] * n_echo  # type: ignore[list-item]
+
+    parent_json = load_json(json_path) if Path(json_path).exists() else {}
+    frames = nii.split_4D_image_to_3D()
+    out_paths: list[Path] = []
+    for i, (frame, te) in enumerate(zip(frames, tes)):
+        new_nii = _with_echo_suffix(nii_path, i)
+        new_json = _with_echo_suffix(json_path, i)
+        frame.save(new_nii)
+        j = dict(parent_json)
+        j.pop("grid", None)
+        if te is not None:
+            j["EchoTime"] = te
+        j["EchoNumbers"] = i + 1
+        secure_save_json(new_json, j, indent=4)
+        _add_grid_info_to_json(new_nii, new_json)
+        out_paths.append(new_nii)
+
+    # Only delete the 4-D originals once every per-echo file is on disk.
+    if all(p.exists() for p in out_paths):
+        for p in (Path(nii_path), Path(json_path)):
+            try:
+                p.unlink()
+            except FileNotFoundError:
+                pass
+    return out_paths
+
+
+def _with_echo_suffix(p: Path, eco_index: int) -> Path:
+    """Insert ``-eco<i>`` into the ``part`` BIDS entity of a NIfTI/JSON filename.
+
+    ``sub-X_..._part-magnitude_dixon.nii.gz`` → ``..._part-magnitude-eco0_dixon.nii.gz``.
+    Falls back to appending before the extension when no ``_part-`` entity is present.
+    """
+    name = p.name
+    import re
+
+    m = re.search(r"_part-([^_]+)_", name)
+    if m:
+        old = m.group(0)
+        new = f"_part-{m.group(1)}-eco{eco_index}_"
+        name = name.replace(old, new, 1)
+    else:
+        for ext in (".nii.gz", ".json"):
+            if name.endswith(ext):
+                name = f"{name[: -len(ext)]}-eco{eco_index}{ext}"
+                break
+    return p.with_name(name)
+
+
 def _add_grid_info_to_json(nii_path: Path | str, simp_json: Path | str, force_update: bool = False, add: bool = True) -> dict:
     """Append grid metadata (shape, spacing, orientation, affine) to a sidecar JSON file.
 
@@ -596,6 +695,80 @@ def _add_grid_info_to_json(nii_path: Path | str, simp_json: Path | str, force_up
     return json_dict
 
 
+_EXTRACT_CACHE_DIR = ".extract_cache"
+
+
+def _folder_fingerprint(folder: Path) -> tuple[int, str] | None:
+    """Cheap folder identity: (file count, sha1 of sorted (relpath, size) list).
+
+    No DICOM headers are read — only ``os.stat`` on each file. Returns
+    ``None`` if the folder can't be scanned. Used by
+    :func:`extract_dicom_folder`'s fast-skip path.
+    """
+    import hashlib
+    import json as _json
+
+    try:
+        entries: list[tuple[str, int]] = []
+        for p in sorted(folder.rglob("*")):
+            if p.is_file():
+                try:
+                    size = p.stat().st_size
+                except OSError:
+                    continue
+                entries.append((str(p.relative_to(folder)), size))
+        payload = _json.dumps(entries, separators=(",", ":")).encode()
+        return len(entries), hashlib.sha1(payload).hexdigest()  # noqa: S324
+    except OSError:
+        return None
+
+
+def _extract_marker_path(source_folder: Path, dataset_path_out: Path) -> Path:
+    """Path of the fast-skip marker for a source folder, kept under the OUTPUT dataset.
+
+    Using the dataset root instead of the source tree means we never write
+    into the input DICOMs (which may be read-only or on removable media).
+    """
+    import hashlib
+
+    key = hashlib.sha1(str(source_folder.resolve()).encode()).hexdigest()  # noqa: S324
+    return dataset_path_out / _EXTRACT_CACHE_DIR / f"{key}.json"
+
+
+def _is_already_extracted(source_folder: Path, dataset_path_out: Path) -> bool:
+    """True when the source folder was extracted before and its file list is unchanged."""
+    import json as _json
+
+    marker = _extract_marker_path(source_folder, dataset_path_out)
+    if not marker.is_file():
+        return False
+    try:
+        prev = _json.loads(marker.read_text())
+    except (OSError, ValueError):
+        return False
+    fp = _folder_fingerprint(source_folder)
+    if fp is None:
+        return False
+    count, h = fp
+    return prev.get("count") == count and prev.get("hash") == h
+
+
+def _write_extract_marker(source_folder: Path, dataset_path_out: Path) -> None:
+    """Record the current file-list fingerprint for the source folder."""
+    import json as _json
+
+    fp = _folder_fingerprint(source_folder)
+    if fp is None:
+        return
+    count, h = fp
+    marker = _extract_marker_path(source_folder, dataset_path_out)
+    try:
+        marker.parent.mkdir(parents=True, exist_ok=True)
+        marker.write_text(_json.dumps({"source": str(source_folder), "count": count, "hash": h}))
+    except OSError:
+        pass  # writing the marker is best-effort; missing it just disables the fast-path
+
+
 def _find_all_files(dcm_dirs: Path | list[Path], verbose=False):
     """Recursively find all DICOM directories or files in the given paths.
 
@@ -811,6 +984,8 @@ def extract_dicom_folder(
     skip_localizer=True,
     parent: str = "rawdata",
     censor_list: list | None = None,
+    skip_already_extracted: bool = True,
+    force_rescan: bool = False,
 ) -> dict:
     """Extract DICOM files from a directory or list of directories, convert them to NIfTI format, and store the output.
 
@@ -835,6 +1010,13 @@ def extract_dicom_folder(
         parent (str, optional): Parent folder inside ``dataset_path_out`` under which subjects are written
             (typically ``"rawdata"``). Defaults to ``"rawdata"``.
         censor_list (list | None, optional): List of series keys to skip entirely. Defaults to an empty list.
+        skip_already_extracted (bool, optional): If True, use a per-source-folder marker under
+            ``dataset_path_out/.extract_cache/`` to skip folders whose file listing (relative
+            paths + sizes) is unchanged since the last successful extraction. No DICOM headers
+            are read for a skipped folder — great for re-runs that only need to pick up newly
+            added subjects. Defaults to True.
+        force_rescan (bool, optional): If True, bypass the fast-skip marker and re-read every
+            DICOM. Defaults to False.
 
     Returns:
         dict: A dictionary with keys representing DICOM series and values as paths to the generated NIfTI files.
@@ -856,6 +1038,18 @@ def extract_dicom_folder(
 
         if str(dicom_path).endswith(".pkl"):
             continue
+        # Fast-skip: identical file listing since last successful extraction
+        # → no DICOM headers read for this folder. Zips are excluded because
+        # their inner file list isn't visible without unpacking.
+        if (
+            skip_already_extracted
+            and not force_rescan
+            and not str(dicom_path).endswith(".zip")
+            and Path(dicom_path).is_dir()
+            and _is_already_extracted(Path(dicom_path), Path(dataset_path_out))
+        ):
+            logger.print(f"Skip {dicom_path} (already extracted; fingerprint matches)", verbose=verbose)
+            continue
         temp_dir = None
         try:
             if str(dicom_path).endswith(".zip"):
@@ -904,6 +1098,15 @@ def process_series(key, files, parts):
                     except Exception:
                         logger.print_error()
 
+            # Record the fingerprint only when the whole folder went through
+            # without an exception AND the source is a real directory (not a
+            # zip mount that's about to disappear). Errors above are caught
+            # per-series so this fires even if individual series were skipped
+            # (e.g. localizers) — but not if `_read_dicom_files` itself raised
+            # (that path lands in the outer `finally` without reaching here).
+            if skip_already_extracted and temp_dir is None and Path(dicom_path).is_dir():
+                _write_extract_marker(Path(dicom_path), Path(dataset_path_out))
+
         finally:
             if temp_dir is not None:
                 shutil.rmtree(temp_dir)

TPTBox/core/dicom/dicom_header_to_keys.py

@@ -39,6 +39,9 @@
     "phase": "phase",
     "mag": "mag",
     "sub": "subtraction",
+    "m_ffe": "magnitude",
+    "p_ffe": "phase",
+    "magnitude": "magnitude",
 }
 dixon_mapping = {**dixon_mapping, **{v: v for v in dixon_mapping.values()}}
 map_series_description_to_file_format_default = {

TPTBox/core/nii_wrapper.py

@@ -454,7 +454,8 @@ def nii(self) -> Nifti1Image:
                 nii2 = Nifti1Image(arr_for_nib, self.affine, self.header)
                 nii2.set_data_dtype(safe_dtype)
                 nii = Nifti1Image(arr_for_nib, nii2.affine, nii2.header)  # type: ignore
-            if all(a is None for a in self.header.get_slope_inter()):
+            c = self.get_c_val()
+            if all(a is None for a in self.header.get_slope_inter()) and not np.isnan(c):
                 nii.header.set_slope_inter(1,self.get_c_val()) # type: ignore
             #if self.header is not None:
             #    self.header.set_sform(self.affine, code=1)

TPTBox/stitching/stitching.py

@@ -118,6 +118,35 @@ def argmin(lst: list) -> int:
     return lst.index(min(lst))
 
 
+def _occupancy_bbox(occ: np.ndarray) -> tuple[np.ndarray, np.ndarray] | None:
+    """Axis-aligned bounding box (min, max inclusive) of the nonzero region of `occ`.
+
+    Uses three 1-D `np.any` reductions rather than `np.where`, so cost is O(N)
+    with cache-friendly access. Returns ``None`` for an all-zero occupancy.
+    """
+    mask = occ > 0
+    axes = mask.ndim
+    lo = np.empty(axes, dtype=np.int64)
+    hi = np.empty(axes, dtype=np.int64)
+    for ax in range(axes):
+        collapsed = np.any(mask, axis=tuple(i for i in range(axes) if i != ax))
+        idx = np.flatnonzero(collapsed)
+        if idx.size == 0:
+            return None
+        lo[ax] = idx[0]
+        hi[ax] = idx[-1]
+    return lo, hi
+
+
+def _aabb_overlaps(a: tuple[np.ndarray, np.ndarray] | None, b: tuple[np.ndarray, np.ndarray] | None) -> bool:
+    """Whether two axis-aligned bounding boxes touch or intersect. ``None`` means empty."""
+    if a is None or b is None:
+        return False
+    a_lo, a_hi = a
+    b_lo, b_hi = b
+    return bool(np.all(a_lo <= b_hi) and np.all(b_lo <= a_hi))
+
+
 def get_max_affine_and_shape(
     points: np.ndarray,
     affines: list[np.ndarray],
@@ -431,6 +460,21 @@ def buffer_reference(path: str | Path, bias_field: bool, crop: bool = False) ->
 }
 
 
+def _auto_output_dtype(niis: list[nib.nifti1.Nifti1Image]) -> type:
+    """Pick the smallest lossless output dtype from the inputs' on-disk dtypes.
+
+    Reads only the NIfTI headers — no pixel data is loaded. If every input is
+    integer-typed, returns the widest integer dtype that covers them all;
+    otherwise falls back to float32. This lets the stitcher store magnitude
+    MR outputs as uint16 (or int16) when the source already fit in 16 bits,
+    halving the on-disk and downstream RAM footprint versus float64.
+    """
+    dtypes = [np.dtype(nii.get_data_dtype()) for nii in niis]
+    if all(np.issubdtype(d, np.integer) for d in dtypes):
+        return max(dtypes, key=lambda d: d.itemsize).type  # e.g. np.uint16
+    return np.float32
+
+
 def main(  # noqa: C901
     images: list[str] | list[Path] | list[nib.nifti1.Nifti1Image],
     output: str | None,
@@ -572,6 +616,11 @@ def main(  # noqa: C901
             dtype2 = np.uint64
         dtype = dtype2
     else:
+        # Auto-detect the output dtype from the input headers when the caller
+        # asked for "auto". Blending math still runs in float; only the final
+        # cast at save time uses the picked dtype (e.g. uint16 for magnitude MR).
+        if isinstance(dtype, str) and dtype == "auto":
+            dtype = _auto_output_dtype(niis)
         dtype2 = float
     nii_out = get_max_affine_and_shape(corners_current, affines, min_spacing=min_spacing, dtype=dtype2, verbose=verbose)
     target_list = []
@@ -582,8 +631,10 @@ def main(  # noqa: C901
         print(f"{i:2}/{len(niis):2} resampled", end="\r") if verbose else None
         nii_new = nip.resample_from_to(nii, nii_out, 0 if is_segmentation else 3, mode="constant", cval=min_value)
         arr_new = get_array(nii_new)
+        if not is_segmentation and np.issubdtype(arr_new.dtype, np.floating):
+            np.nan_to_num(arr_new, copy=False, nan=min_value, posinf=min_value, neginf=min_value)
         target_list.append(arr_new)
-        b = nib.Nifti1Image(get_array(nii) * 0 + 1, affine=nii.affine)  # type: ignore
+        b = nib.Nifti1Image(np.ones(nii.shape, dtype=np.float32), affine=nii.affine)  # type: ignore
         b = nip.resample_from_to(b, nii_new, 0, cval=0, mode="constant")
         if is_segmentation:
             x = arr_new > 0
@@ -593,10 +644,18 @@ def main(  # noqa: C901
             occupancy_list.append(get_array(b).astype(np.float32))
 
     print("\n### ramp stitching ###") if verbose else None
+    # Per-chunk axis-aligned bounding box in target-space voxel coords.
+    # Precomputing once avoids the O(N_chunks^2) full-volume copy + multiply
+    # inside the ramp loop for pairs that don't touch. `_occupancy_bbox`
+    # returns None for an empty occupancy — treated as "no overlap possible".
+    bboxes = [_occupancy_bbox(occ) for occ in occupancy_list]
     # ramp stitching
     combinations = list(itertools.combinations(range(len(target_list)), 2))
     for idx, item in enumerate(combinations, 1):
         print(f"{idx:2}/{len(combinations):2} ramp stitching", end="\r") if verbose else None
+        # Skip disjoint pairs before touching the full-volume arrays.
+        if not _aabb_overlaps(bboxes[item[0]], bboxes[item[1]]):
+            continue
         # TODO fix intersection with more than two occupancies
         arr_1_full = occupancy_list[item[0]]
         arr_2_full = occupancy_list[item[1]]