From 831048d7750731297538ca4f01071428184646a9 Mon Sep 17 00:00:00 2001 From: ClePol Date: Wed, 22 Jul 2026 14:49:38 +0200 Subject: [PATCH] added CC editing --- CorpusCallosum/README.md | 79 +++++- CorpusCallosum/fastsurfer_cc.py | 234 +++++++++++++----- CorpusCallosum/registration/landmarks.py | 102 ++++++++ .../midsagittal_plane_alignment.py | 43 +++- CorpusCallosum/utils/editing.py | 106 ++++++++ CorpusCallosum/utils/mapping_helpers.py | 66 ++++- doc/overview/EDITING.md | 72 +++++- doc/overview/OUTPUT_FILES.md | 21 +- doc/overview/modules/CC.md | 2 + doc/scripts/fastsurfer_cc.rst | 133 +++++++++- run_fastsurfer.sh | 42 +++- 11 files changed, 813 insertions(+), 87 deletions(-) create mode 100644 CorpusCallosum/registration/landmarks.py create mode 100644 CorpusCallosum/utils/editing.py diff --git a/CorpusCallosum/README.md b/CorpusCallosum/README.md index cc7507153..576f6e77f 100644 --- a/CorpusCallosum/README.md +++ b/CorpusCallosum/README.md @@ -10,8 +10,83 @@ The documentation is split into three files, please refer to: Quickstart ---------- +The expert command expects an existing FastSurfer subject containing at least `mri/orig.mgz` and the FastSurfer +segmentation. It reads those inputs and writes all CC outputs into the subject directory. + +### Native + +```bash +python3 CorpusCallosum/fastsurfer_cc.py \ + --sd /data/subjects \ + --sid sub001 \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png +``` + +Paths that are not absolute are resolved relative to `/data/subjects/sub001`. The upright volume is useful for +checking the midplane and is the anatomical reference for manual CC editing. + +### Docker + +The FastSurfer Docker image normally starts `run_fastsurfer.sh`. Expert commands therefore override the entrypoint +with FastSurfer's environment-setup wrapper and explicitly invoke `fastsurfer_cc.py`: + +```bash +SUBJECTS_DIR=/data/fastsurfer +SID=sub001 + +docker run --gpus all --rm \ + --user "$(id -u):$(id -g)" \ + --volume "$SUBJECTS_DIR:/output" \ + --entrypoint /fastsurfer/tools/Docker/entrypoint.sh \ + deepmi/fastsurfer:latest \ + python3 /fastsurfer/CorpusCallosum/fastsurfer_cc.py \ + --sd /output \ + --sid "$SID" \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png +``` + +Use an image tag containing the FastSurfer-CC features documented here. For CPU execution, use the matching +`deepmi/fastsurfer:cpu-latest` image and omit `--gpus all`. The bind-mounted output directory must be writable by the +user passed with `--user`. + +### Singularity or Apptainer + +Build or download an image as described in the [Singularity documentation](../doc/overview/SINGULARITY.md), then bind +the subject directory and invoke the expert script directly: + ```bash -python3 fastsurfer_cc.py --sd /path/to/fastsurfer/output --sid test-case --verbose +SUBJECTS_DIR=/data/fastsurfer +SID=sub001 +FASTSURFER_SIF=/containers/fastsurfer-gpu.sif + +singularity exec --nv --no-mount home,cwd -e \ + --bind "$SUBJECTS_DIR:/output" \ + "$FASTSURFER_SIF" \ + python3 /fastsurfer/CorpusCallosum/fastsurfer_cc.py \ + --sd /output \ + --sid "$SID" \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png ``` -Gives all standard outputs. The corpus callosum morphometry can be found at `stats/callosum.CC.midslice.json` including 100 thickness measurements and the areas of sub-segments. +The same command works with Apptainer by replacing `singularity` with `apptainer`. For CPU execution, omit `--nv` and +use an appropriate CPU image. `--no-mount home,cwd -e` prevents host Python packages and environment variables from +leaking into the container. + +These commands generate all standard CC outputs. Morphometry is written to `stats/callosum.CC.midslice.json`, +including 100 thickness measurements and the areas of subsegments. + +The expert interface also supports supplied 3D AC/PC voxel coordinates and manual corrections of the upright CC +segmentation. Copy `mri/callosum.CC.upright.mgz` to `mri/callosum.CC.upright.manedit.mgz`, edit label 192 using +`mri/upright_volume.mgz` as the reference, and pass the edited path to `--segmentation_manedit`. A manual segmentation +containing any fornix label-250 voxels is self-contained and authoritative for the entire fornix. If it has no label +250, the automatic upright and original-space CC segmentations from a previous run are required and supply the +fornix. Choose supplied AC/PC coordinates before creating the edit and reuse the exact same coordinates and midplane +method for every edit rerun; otherwise regenerate the upright reference and recreate or rebase the correction. See the +[advanced documentation](../doc/scripts/fastsurfer_cc.rst) for complete native, Docker, and Singularity/Apptainer +commands. diff --git a/CorpusCallosum/fastsurfer_cc.py b/CorpusCallosum/fastsurfer_cc.py index a704b69dd..216712b5f 100644 --- a/CorpusCallosum/fastsurfer_cc.py +++ b/CorpusCallosum/fastsurfer_cc.py @@ -33,6 +33,7 @@ CC_LABEL, DEFAULT_INPUT_PATHS, DEFAULT_OUTPUT_PATHS, + FORNIX_LABEL, THIRD_VENTRICLE_LABEL, ) from CorpusCallosum.data.read_write import MGHHeaderDict, convert_numpy_to_json_serializable @@ -47,10 +48,12 @@ offset_affine, recon_cc_surf_measures_multi, ) +from CorpusCallosum.utils.editing import add_file_suffix, load_manual_upright_segmentation, validate_landmarks_in_image from CorpusCallosum.utils.mapping_helpers import ( apply_transform_to_pt, apply_transform_to_volume, calc_mapping_to_standard_space, + map_hard_segmentation_to_orig, map_softlabels_to_orig, ) from CorpusCallosum.utils.types import SliceSelection, SubdivisionMethod @@ -204,6 +207,13 @@ def _slice_selection(a: str) -> SliceSelection: help="Output path for corpus callosum and fornix segmentation output.", default=Path(DEFAULT_OUTPUT_PATHS["segmentation"]), ) + advanced.add_argument( + "--segmentation_manedit", + type=path_or_none, + help="Edited upright CC segmentation. If it contains fornix label 250, it is self-contained; otherwise the " + "fornix is retained from the automatic segmentation output of a previous run.", + default=None, + ) advanced.add_argument( "--segmentation_in_orig", type=path_or_none, @@ -253,6 +263,22 @@ def _slice_selection(a: str) -> SliceSelection: "'fsaverage_symmetry': fsaverage alignment + LR label-symmetry shift refinement (default); " "'fsaverage_distance_map': fsaverage alignment + distance-map plane-fitting refinement.", ) + advanced.add_argument( + "--ac_coords", + type=float, + nargs=3, + metavar=("X", "Y", "Z"), + help="Optional AC point as three floating-point voxel coordinates in orig.mgz space. Requires --pc_coords.", + default=None, + ) + advanced.add_argument( + "--pc_coords", + type=float, + nargs=3, + metavar=("X", "Y", "Z"), + help="Optional PC point as three floating-point voxel coordinates in orig.mgz space. Requires --ac_coords.", + default=None, + ) advanced.add_argument( "--qc_image", type=path_or_none, @@ -332,6 +358,9 @@ def options_parse() -> argparse.Namespace: parser = make_parser() args = parser.parse_args() + if (args.ac_coords is None) != (args.pc_coords is None): + parser.error("--ac_coords and --pc_coords must be supplied together.") + # Reconstruct subject_dir from sd and sid (but sd might be stored as out_dir by parser_defaults) sd_value = getattr(args, "out_dir", None) if sd_value and hasattr(args, "sid") and args.sid: @@ -379,6 +408,7 @@ def options_parse() -> argparse.Namespace: all_paths = ( "segmentation", + "segmentation_manedit", "segmentation_in_orig", "cc_measures", "upright_lta", @@ -539,6 +569,7 @@ def main( device: str | torch.device = "auto", upright_volume: str | Path | None = None, segmentation: str | Path | None = None, + segmentation_manedit: str | Path | None = None, cc_measures: str | Path | None = None, cc_mid_measures: str | Path | None = None, upright_lta: str | Path | None = None, @@ -554,6 +585,8 @@ def main( softlabels_cc: str | Path | None = None, softlabels_fn: str | Path | None = None, softlabels_background: str | Path | None = None, + ac_coords: Iterable[float] | None = None, + pc_coords: Iterable[float] | None = None, ) -> Literal[0] | str: """Main pipeline function for corpus callosum analysis. @@ -587,6 +620,9 @@ def main( Path to save upright volume. segmentation : str or Path, optional Path to save segmentation. + segmentation_manedit : str or Path, optional + Path to an edited upright CC segmentation. A file containing fornix label 250 is self-contained; otherwise + the fornix is retained from the automatic segmentation outputs of a previous run. cc_measures : str or Path, optional Path to save post-processing results. cc_mid_measures : str or Path, optional @@ -615,6 +651,10 @@ def main( Path to save fornix soft labels. softlabels_background : str or Path, optional Path to save background soft labels. + ac_coords : tuple of float, optional + Supplied AC voxel coordinates in original image space; requires ``pc_coords``. + pc_coords : tuple of float, optional + Supplied PC voxel coordinates in original image space; requires ``ac_coords``. Returns ------- @@ -664,6 +704,7 @@ def main( save_template_dir=save_template_dir, upright_volume=upright_volume, cc_segmentation=segmentation, + cc_segmentation_manedit=segmentation_manedit, cc_measures=cc_measures, cc_mid_measures=cc_mid_measures, upright_lta=upright_lta, @@ -689,13 +730,20 @@ def main( #### setup variables futures = [] - # load models + manual_edit = sd.has_attribute("cc_segmentation_manedit") + supplied_landmarks = ac_coords is not None or pc_coords is not None + + # load only models needed by the selected path device = find_device(device) logger.info(f"Using device: {device}") logger.info("Loading models") - _model_localization = thread_executor().submit(localization_inference.load_model, device=device) - _model_segmentation = thread_executor().submit(segmentation_inference.load_model, device=device) + _model_localization = ( + None if supplied_landmarks else thread_executor().submit(localization_inference.load_model, device=device) + ) + _model_segmentation = ( + None if manual_edit else thread_executor().submit(segmentation_inference.load_model, device=device) + ) _aseg_fut = thread_executor().submit(nib.load, sd.filename_by_attribute("aseg_name")) orig = cast(nibabelImage, nib.load(sd.conf_name)) @@ -711,6 +759,12 @@ def main( logger.error(error_message) return error_message + try: + ac_coords_orig, pc_coords_orig = validate_landmarks_in_image(ac_coords, pc_coords, orig.shape) + except ValueError as e: + logger.error(str(e)) + return str(e) + # Analysis-width slab around the midplane (guaranteed to be aligned RAS by as_closest_canonical). vox_size_ras: tuple[float, float, float] = nib.as_closest_canonical(orig).header.get_zooms() vox_size = vox_size_ras[0], vox_size_ras[2], vox_size_ras[1] # convert from RAS to LIA @@ -732,7 +786,17 @@ def main( logger.error(error_message) return error_message - midplane = find_midplane_transform(orig=orig, aseg_img=aseg_img, midplane_method=midplane_method) + try: + midplane = find_midplane_transform( + orig=orig, + aseg_img=aseg_img, + midplane_method=midplane_method, + ac_coords_orig=ac_coords_orig, + pc_coords_orig=pc_coords_orig, + ) + except ValueError as e: + logger.error(str(e)) + return str(e) orig2fsavg_vox2vox = midplane.orig2fsavg_vox2vox fsavg_vox2ras = midplane.fsavg_vox2ras _fsavg_header_dict = midplane.fsavg_header_dict @@ -771,18 +835,25 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4: fsavg2midslab_vox2vox = offset_affine([slices_to_analyze // 2, 0, 0]) @ fsavg2midslice_vox2vox fsaverage_midslab_vox2ras: AffineMatrix4x4 = fsavg_vox2ras @ np.linalg.inv(fsavg2midslab_vox2vox) - #### do localization and segmentation inference - logger.info("Starting AC/PC localization") - target_shape: tuple[int, int, int] = (slices_to_analyze, fsavg_header["dims"][1], fsavg_header["dims"][2]) - # predict ac and pc coordinates in upright AS space - ac_coords_vox, pc_coords_vox = localize_ac_pc( - np.asarray(orig.dataobj), - aseg_img, - _orig2midslab_vox2vox(additional_context=2), - _model_localization.result(), - target_shape, - ) - logger.info("Starting corpus callosum segmentation") + #### resolve landmarks and segmentation + target_shape: Shape3d = (slices_to_analyze, fsavg_header["dims"][1], fsavg_header["dims"][2]) + if ac_coords_orig is not None and pc_coords_orig is not None: + logger.info("Using supplied AC/PC coordinates from orig.mgz voxel space") + ac_coords_vox = apply_transform_to_pt(ac_coords_orig, orig2midslice_vox2vox)[1:] + pc_coords_vox = apply_transform_to_pt(pc_coords_orig, orig2midslice_vox2vox)[1:] + landmark_source = "supplied" + else: + logger.info("Starting AC/PC localization") + assert _model_localization is not None + ac_coords_vox, pc_coords_vox = localize_ac_pc( + np.asarray(orig.dataobj), + aseg_img, + _orig2midslab_vox2vox(additional_context=2), + _model_localization.result(), + target_shape, + ) + landmark_source = "model" + num_context = 8 # 8 extra in x-direction for context slices target_shape: Shape3d = (slices_to_analyze + num_context, fsavg_header["dims"][1], fsavg_header["dims"][2]) midslices: Image3d = affine_transform( @@ -796,42 +867,60 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4: cval=0, prefilter=True, # unclear, why we are using a smoothing filter here ) - cc_fn_seg_labels, cc_fn_softlabels = segment_cc( - midslices, - ac_coords_vox, - pc_coords_vox, - aseg_img, - _model_segmentation.result(), - ) + cc_fn_softlabels: Image4d | None + manual_has_fornix = False + if manual_edit: + logger.info(f"Using manual upright CC segmentation {sd.filename_by_attribute('cc_segmentation_manedit')}") + try: + cc_fn_seg_labels, manual_has_fornix = load_manual_upright_segmentation( + automatic_path=( + sd.filename_by_attribute("cc_segmentation") if sd.has_attribute("cc_segmentation") else None + ), + manual_path=sd.filename_by_attribute("cc_segmentation_manedit"), + expected_shape=(slices_to_analyze, fsavg_header["dims"][1], fsavg_header["dims"][2]), + expected_affine=fsaverage_midslab_vox2ras, + ) + except ValueError as e: + logger.error(str(e)) + return str(e) + cc_fn_softlabels = None + segmentation_source = "manual" + else: + logger.info("Starting corpus callosum segmentation") + assert _model_segmentation is not None + cc_fn_seg_labels, cc_fn_softlabels = segment_cc( + midslices, + ac_coords_vox, + pc_coords_vox, + aseg_img, + _model_segmentation.result(), + ) + segmentation_source = "model" + + # save automatic segmentation softlabels + for i, (attr, name) in enumerate((("background",) * 2, ("cc", "Corpus Callosum"), ("fn", "Fornix"))): + if sd.has_attribute(f"cc_softlabels_{attr}"): + logger.info(f"Saving {name} softlabels to {sd.filename_by_attribute(f'cc_softlabels_{attr}')}") + futures.append( + thread_executor().submit( + nib.save, + nib.MGHImage(cc_fn_softlabels[..., i], fsaverage_midslab_vox2ras, orig.header), + sd.filename_by_attribute(f"cc_softlabels_{attr}"), + ) + ) - # save segmentation softlabels - for i, (attr, name) in enumerate((("background",) * 2, ("cc", "Corpus Callosum"), ("fn", "Fornix"))): - if sd.has_attribute(f"cc_softlabels_{attr}"): - logger.info(f"Saving {name} softlabels to {sd.filename_by_attribute(f'cc_softlabels_{attr}')}") + if sd.has_attribute("cc_segmentation"): + _cc_seg_path = sd.filename_by_attribute("cc_segmentation") + _cc_seg_path.parent.mkdir(exist_ok=True, parents=True) + logger.info(f"Saving CC segmentation to {_cc_seg_path}") futures.append( thread_executor().submit( nib.save, - nib.MGHImage(cc_fn_softlabels[..., i], fsaverage_midslab_vox2ras, orig.header), - sd.filename_by_attribute(f"cc_softlabels_{attr}"), + nib.MGHImage(cc_fn_seg_labels, fsaverage_midslab_vox2ras, orig.header), + _cc_seg_path, ) ) - # Create a temporary segmentation image with proper affine for enhanced postprocessing - # Process slices based on selection mode - - # save segmentation labels - if sd.has_attribute("cc_segmentation"): - _cc_seg_path = sd.filename_by_attribute("cc_segmentation") - _cc_seg_path.parent.mkdir(exist_ok=True, parents=True) - logger.info(f"Saving CC segmentation to {_cc_seg_path}") - futures.append( - thread_executor().submit( - nib.save, - nib.MGHImage(cc_fn_seg_labels, fsaverage_midslab_vox2ras, orig.header), - _cc_seg_path, - ) - ) - logger.info(f"Processing slices with selection mode: {slice_selection}") try: ( @@ -919,7 +1008,7 @@ def _select_middle_valid_slice(): selected_slice_position, selected_slice_idx, middle_slice_result = _select_middle_valid_slice() - # map soft labels to original space (in parallel because this takes a while, and we only do it to save the labels) + # Map automatic soft labels or the edited hard CC mask to original space. if sd.has_attribute("cc_orig_segfile"): if middle_slice_result is not None and len(middle_slice_result["split_contours"]) <= 5: cc_subseg_midslice = make_subdivision_mask( @@ -935,17 +1024,45 @@ def _select_middle_valid_slice(): cc_subseg_midslice = None # if num_threads is not large enough (>1), this might be blocking ; serial_executor runs the function in submit executor = thread_executor() if get_num_threads() > 2 else serial_executor() - futures.append( - executor.submit( - map_softlabels_to_orig, - cc_fn_softlabels=cc_fn_softlabels, - orig=orig, - orig_space_segmentation_path=sd.filename_by_attribute("cc_orig_segfile"), - orig2slab_vox2vox=_orig2midslab_vox2vox(), - cc_subseg_midslice=cc_subseg_midslice, - orig2midslice_vox2vox=orig2midslice_vox2vox, + if manual_edit: + automatic_orig_path = sd.filename_by_attribute("cc_orig_segfile") + if not manual_has_fornix and not automatic_orig_path.is_file(): + error_message = ( + f"Manual upright segmentation does not contain fornix label {FORNIX_LABEL}, and automatic " + f"original-space segmentation {automatic_orig_path} is missing. Supply the fornix in the manual " + "segmentation or run FastSurfer-CC once without edits." + ) + logger.error(error_message) + return error_message + if cc_subseg_midslice is None: + error_message = "Cannot map the edited CC to original space without a valid middle-slice subdivision." + logger.error(error_message) + return error_message + futures.append( + executor.submit( + map_hard_segmentation_to_orig, + edited_segmentation=cc_fn_seg_labels, + reference_image=orig, + orig2slab_vox2vox=_orig2midslab_vox2vox(), + cc_subseg_midslice=cc_subseg_midslice, + orig2midslice_vox2vox=orig2midslice_vox2vox, + output_path=add_file_suffix(automatic_orig_path, "manedit"), + automatic_orig_segmentation=None if manual_has_fornix else nib.load(automatic_orig_path), + ) + ) + else: + assert cc_fn_softlabels is not None + futures.append( + executor.submit( + map_softlabels_to_orig, + cc_fn_softlabels=cc_fn_softlabels, + orig=orig, + orig_space_segmentation_path=sd.filename_by_attribute("cc_orig_segfile"), + orig2slab_vox2vox=_orig2midslab_vox2vox(), + cc_subseg_midslice=cc_subseg_midslice, + orig2midslice_vox2vox=orig2midslice_vox2vox, + ) ) - ) metrics: tuple[CCMeasures] = get_args(CCMeasures) @@ -1023,6 +1140,8 @@ def _select_middle_valid_slice(): additional_metrics["slice_selection"] = slice_selection additional_metrics["midline_refine_shift_vox"] = float(midline_shift_vox) additional_metrics["midline_refine_diagnostics"] = midline_shift_diagnostics + additional_metrics["landmark_source"] = landmark_source + additional_metrics["segmentation_source"] = segmentation_source # QC checks if len(outer_contours) > 1 and cc_volume_contour is not None: @@ -1179,6 +1298,7 @@ def save_cc_measures_json(cc_mid_measure_file: Path, metrics: dict[str, object]) device=options.device, upright_volume=options.upright_volume, segmentation=options.segmentation, + segmentation_manedit=options.segmentation_manedit, cc_measures=options.cc_measures, cc_mid_measures=options.cc_mid_measures, upright_lta=options.upright_lta, @@ -1194,5 +1314,7 @@ def save_cc_measures_json(cc_mid_measure_file: Path, metrics: dict[str, object]) softlabels_cc=options.softlabels_cc, softlabels_fn=options.softlabels_fn, softlabels_background=options.softlabels_background, + ac_coords=options.ac_coords, + pc_coords=options.pc_coords, ) ) diff --git a/CorpusCallosum/registration/landmarks.py b/CorpusCallosum/registration/landmarks.py new file mode 100644 index 000000000..346c3293b --- /dev/null +++ b/CorpusCallosum/registration/landmarks.py @@ -0,0 +1,102 @@ +import nibabel as nib +import numpy as np + +from CorpusCallosum.data.read_write import convert_numpy_to_json_serializable +from FastSurferCNN.utils import AffineMatrix4x4, logging + +logger = logging.get_logger(__name__) + + +def _rotation_from_vectors(source: np.ndarray, target: np.ndarray) -> np.ndarray: + """Compute a 3x3 rotation matrix mapping one direction onto another.""" + source = source / np.linalg.norm(source) + target = target / np.linalg.norm(target) + cross = np.cross(source, target) + dot = float(np.clip(np.dot(source, target), -1.0, 1.0)) + cross_norm = np.linalg.norm(cross) + if cross_norm < 1e-8: + return np.eye(3, dtype=float) + kx, ky, kz = cross / cross_norm + skew = np.array( + [ + [0.0, -kz, ky], + [kz, 0.0, -kx], + [-ky, kx, 0.0], + ] + ) + angle = np.arccos(dot) + return np.eye(3, dtype=float) + np.sin(angle) * skew + (1 - np.cos(angle)) * (skew @ skew) + + +def _affine_from_rotation_and_center(rotation: np.ndarray, center: np.ndarray) -> AffineMatrix4x4: + """Construct an affine applying a rotation about a voxel-space center.""" + affine = np.eye(4, dtype=float) + affine[:3, :3] = rotation + affine[:3, 3] = center - rotation @ center + return affine + + +def adjust_midplane_to_landmarks( + orig2fsavg_vox2vox: AffineMatrix4x4, + ac_coords_orig: np.ndarray, + pc_coords_orig: np.ndarray, + base_middle_vox: float, + warning_tilt_deg: float = 15.0, +) -> tuple[AffineMatrix4x4, dict[str, object]]: + """Minimally adjust a midsagittal transform to contain supplied AC/PC points.""" + ac_orig = np.asarray(ac_coords_orig, dtype=float) + pc_orig = np.asarray(pc_coords_orig, dtype=float) + if ac_orig.shape != (3,) or pc_orig.shape != (3,): + raise ValueError("AC and PC coordinates must each contain exactly three values.") + if not np.all(np.isfinite((ac_orig, pc_orig))): + raise ValueError("AC and PC coordinates must be finite.") + + ac_fsavg, pc_fsavg = nib.affines.apply_affine(orig2fsavg_vox2vox, (ac_orig, pc_orig)) + acpc_direction = pc_fsavg - ac_fsavg + acpc_length = float(np.linalg.norm(acpc_direction)) + if acpc_length < 1e-6: + raise ValueError("AC and PC coordinates must be distinct.") + acpc_unit = acpc_direction / acpc_length + + initial_normal = np.array([1.0, 0.0, 0.0], dtype=float) + adjusted_normal = initial_normal - np.dot(initial_normal, acpc_unit) * acpc_unit + adjusted_normal_norm = float(np.linalg.norm(adjusted_normal)) + if adjusted_normal_norm < 1e-6: + raise ValueError( + "The AC-PC line is parallel to the midsagittal plane normal; " + "a stable sagittal plane containing both points cannot be constructed." + ) + adjusted_normal /= adjusted_normal_norm + if np.dot(adjusted_normal, initial_normal) < 0: + adjusted_normal *= -1 + + tilt_deg = float( + np.degrees(np.arccos(np.clip(np.dot(adjusted_normal, initial_normal), -1.0, 1.0))) + ) + if tilt_deg > warning_tilt_deg: + logger.warning( + "Supplied AC/PC points require a large midsagittal plane adjustment (%.2f degrees). " + "Check that the coordinates are in orig.mgz voxel space.", + tilt_deg, + ) + + midpoint = (ac_fsavg + pc_fsavg) / 2.0 + rotation = _rotation_from_vectors(adjusted_normal, initial_normal) + rotate_affine = _affine_from_rotation_and_center(rotation, midpoint) + rotated_midpoint = nib.affines.apply_affine(rotate_affine, midpoint) + translate_affine = np.eye(4, dtype=float) + translate_affine[0, 3] = base_middle_vox - rotated_midpoint[0] + update_affine = translate_affine @ rotate_affine + updated_vox2vox = update_affine @ orig2fsavg_vox2vox + + ac_adjusted, pc_adjusted = nib.affines.apply_affine(updated_vox2vox, (ac_orig, pc_orig)) + residuals = [float(ac_adjusted[0] - base_middle_vox), float(pc_adjusted[0] - base_middle_vox)] + diagnostics: dict[str, object] = { + "landmark_source": "supplied", + "ac_coords_orig_vox": ac_orig.tolist(), + "pc_coords_orig_vox": pc_orig.tolist(), + "plane_tilt_deg": tilt_deg, + "off_plane_residuals_vox": residuals, + "update_affine": convert_numpy_to_json_serializable(update_affine), + } + return updated_vox2vox, diagnostics diff --git a/CorpusCallosum/registration/midsagittal_plane_alignment.py b/CorpusCallosum/registration/midsagittal_plane_alignment.py index 5039452e1..5eecb2fba 100644 --- a/CorpusCallosum/registration/midsagittal_plane_alignment.py +++ b/CorpusCallosum/registration/midsagittal_plane_alignment.py @@ -11,6 +11,7 @@ from CorpusCallosum.data.constants import FSAVERAGE_MIDDLE, FSAVERAGE_REGISTRATION_LABELS, FSAVERAGE_TARGET_PATH from CorpusCallosum.data.read_write import MGHHeaderDict, convert_numpy_to_json_serializable +from CorpusCallosum.registration.landmarks import adjust_midplane_to_landmarks from CorpusCallosum.shape.postprocessing import offset_affine from FastSurferCNN.utils import AffineMatrix4x4, Shape3d, logging, nibabelImage from FastSurferCNN.utils.brainvolstats import hemi_masks_from_aseg @@ -795,6 +796,8 @@ def find_midplane_transform( orig: nibabelImage, aseg_img: nibabelImage, midplane_method: str, + ac_coords_orig: np.ndarray | None = None, + pc_coords_orig: np.ndarray | None = None, ) -> MidplaneTransformResult: """Resolve the fsaverage midplane transform for a selected refinement method. @@ -817,8 +820,32 @@ def find_midplane_transform( vox_size = vox_size_ras[0], vox_size_ras[2], vox_size_ras[1] aseg_data = np.asarray(aseg_img.dataobj) - if midplane_method == "center": + if (ac_coords_orig is None) != (pc_coords_orig is None): + raise ValueError("AC and PC coordinates must be supplied together.") + + def _adjust_with_landmarks(result: MidplaneTransformResult) -> MidplaneTransformResult: + if ac_coords_orig is None or pc_coords_orig is None: + return result + adjusted_transform, landmark_diagnostics = adjust_midplane_to_landmarks( + result.orig2fsavg_vox2vox, + ac_coords_orig, + pc_coords_orig, + result.base_middle_vox, + ) + diagnostics = dict(result.midline_shift_diagnostics) + diagnostics["landmark_adjustment"] = landmark_diagnostics return MidplaneTransformResult( + orig2fsavg_vox2vox=adjusted_transform, + fsavg_vox2ras=result.fsavg_vox2ras, + fsavg_header_dict=result.fsavg_header_dict, + fsavg_shape=result.fsavg_shape, + base_middle_vox=result.base_middle_vox, + midline_shift_vox=result.midline_shift_vox, + midline_shift_diagnostics=diagnostics, + ) + + if midplane_method == "center": + return _adjust_with_landmarks(MidplaneTransformResult( orig2fsavg_vox2vox=np.eye(4), fsavg_vox2ras=orig.affine, fsavg_header_dict={"dims": list(orig.shape[:3])}, @@ -826,7 +853,7 @@ def find_midplane_transform( base_middle_vox=orig.shape[0] / 2.0, midline_shift_vox=0.0, midline_shift_diagnostics={}, - ) + )) orig2fsavg_vox2vox, _, fsavg_vox2ras, fsavg_header_dict = register_centroids_to_fsavg(aseg_img) fsavg_shape = tuple(fsavg_header_dict["dims"]) @@ -839,7 +866,7 @@ def find_midplane_transform( fsavg_shape=fsavg_shape, base_middle_vox=base_middle_vox, ) - return MidplaneTransformResult( + return _adjust_with_landmarks(MidplaneTransformResult( orig2fsavg_vox2vox=dm_result.updated_vox2vox, fsavg_vox2ras=fsavg_vox2ras, fsavg_header_dict=fsavg_header_dict, @@ -847,7 +874,7 @@ def find_midplane_transform( base_middle_vox=base_middle_vox, midline_shift_vox=dm_result.center_shift_vox, midline_shift_diagnostics=dm_result.diagnostics, - ) + )) if midplane_method == "fsaverage_symmetry": orig2fsavg_vox2vox, lr_shift, _, diagnostics = refine_midline_lr_shift( @@ -856,7 +883,7 @@ def find_midplane_transform( fsavg_shape=fsavg_shape, base_middle_vox=base_middle_vox, ) - return MidplaneTransformResult( + return _adjust_with_landmarks(MidplaneTransformResult( orig2fsavg_vox2vox=orig2fsavg_vox2vox, fsavg_vox2ras=fsavg_vox2ras, fsavg_header_dict=fsavg_header_dict, @@ -864,10 +891,10 @@ def find_midplane_transform( base_middle_vox=base_middle_vox, midline_shift_vox=float(lr_shift), midline_shift_diagnostics=diagnostics, - ) + )) if midplane_method == "fsaverage": - return MidplaneTransformResult( + return _adjust_with_landmarks(MidplaneTransformResult( orig2fsavg_vox2vox=orig2fsavg_vox2vox, fsavg_vox2ras=fsavg_vox2ras, fsavg_header_dict=fsavg_header_dict, @@ -875,6 +902,6 @@ def find_midplane_transform( base_middle_vox=base_middle_vox, midline_shift_vox=0.0, midline_shift_diagnostics={}, - ) + )) raise ValueError(f"Unsupported midplane_method: {midplane_method!r}") diff --git a/CorpusCallosum/utils/editing.py b/CorpusCallosum/utils/editing.py new file mode 100644 index 000000000..c687554c4 --- /dev/null +++ b/CorpusCallosum/utils/editing.py @@ -0,0 +1,106 @@ +from collections.abc import Iterable +from pathlib import Path + +import nibabel as nib +import numpy as np + +from CorpusCallosum.data.constants import CC_LABEL, FORNIX_LABEL +from FastSurferCNN.utils import AffineMatrix4x4, Image3d, Shape3d + + +def add_file_suffix(path: str | Path, suffix: str) -> Path: + """Insert a suffix before a file extension, including compound NIfTI extensions.""" + path = Path(path) + extension = ".nii.gz" if path.name.endswith(".nii.gz") else path.suffix + stem = path.name[:-len(extension)] if extension else path.name + return path.with_name(f"{stem}.{suffix}{extension}") + + +def validate_landmarks_in_image( + ac_coords: Iterable[float] | None, + pc_coords: Iterable[float] | None, + image_shape: tuple[int, ...], +) -> tuple[np.ndarray | None, np.ndarray | None]: + """Validate optional paired landmarks in original-image voxel coordinates.""" + if (ac_coords is None) != (pc_coords is None): + raise ValueError("AC and PC coordinates must be supplied together.") + if ac_coords is None or pc_coords is None: + return None, None + ac = np.asarray(ac_coords, dtype=float) + pc = np.asarray(pc_coords, dtype=float) + if ac.shape != (3,) or pc.shape != (3,): + raise ValueError("AC and PC coordinates must each contain exactly three values.") + if not np.all(np.isfinite((ac, pc))): + raise ValueError("AC and PC coordinates must be finite.") + upper = np.asarray(image_shape[:3], dtype=float) - 1 + if np.any(ac < 0) or np.any(pc < 0) or np.any(ac > upper) or np.any(pc > upper): + raise ValueError( + f"AC and PC coordinates must lie inside orig.mgz voxel bounds [0, {upper.tolist()}]." + ) + if np.linalg.norm(pc - ac) < 1e-6: + raise ValueError("AC and PC coordinates must be distinct.") + return ac, pc + + +def load_manual_upright_segmentation( + automatic_path: Path | None, + manual_path: Path, + expected_shape: Shape3d, + expected_affine: AffineMatrix4x4, +) -> tuple[Image3d, bool]: + """Load an upright edit, using its fornix or falling back to the automatic result.""" + if not manual_path.is_file(): + raise ValueError(f"Manual upright segmentation {manual_path} does not exist.") + + manual_img = nib.load(manual_path) + if manual_img.shape[:3] != expected_shape: + raise ValueError( + f"The manual upright segmentation has shape {manual_img.shape[:3]}, expected {expected_shape}. " + "Use the same midplane method and supplied AC/PC coordinates as the original run." + ) + if not np.allclose(manual_img.affine, expected_affine): + raise ValueError( + "The manual upright segmentation affine does not match the current midsagittal slab. " + "Use the same midplane method and supplied AC/PC coordinates as the original run." + ) + + manual = np.asarray(manual_img.dataobj) + unknown_labels = set(np.unique(manual).astype(int).tolist()) - {0, CC_LABEL, FORNIX_LABEL} + if unknown_labels: + raise ValueError( + f"Manual upright segmentation contains unsupported labels {sorted(unknown_labels)}; " + f"only 0, {CC_LABEL}, and {FORNIX_LABEL} are allowed." + ) + manual_cc = manual == CC_LABEL + if not np.any(manual_cc): + raise ValueError("Manual upright segmentation does not contain any corpus callosum voxels (label 192).") + + manual_has_fornix = np.any(manual == FORNIX_LABEL) + if manual_has_fornix: + return manual.astype(np.uint8), True + + if automatic_path is None or not automatic_path.is_file(): + automatic_location = automatic_path if automatic_path is not None else "not configured" + raise ValueError( + f"Manual upright segmentation does not contain fornix label {FORNIX_LABEL}, and automatic upright " + f"segmentation ({automatic_location}) is missing. Supply the fornix in the manual segmentation or run " + "FastSurfer-CC once without edits." + ) + + automatic_img = nib.load(automatic_path) + if automatic_img.shape[:3] != expected_shape: + raise ValueError( + f"The automatic upright segmentation has shape {automatic_img.shape[:3]}, expected {expected_shape}. " + "Use the same midplane method and supplied AC/PC coordinates as the original run." + ) + if not np.allclose(automatic_img.affine, expected_affine): + raise ValueError( + "The automatic upright segmentation affine does not match the current midsagittal slab. " + "Use the same midplane method and supplied AC/PC coordinates as the original run." + ) + + automatic = np.asarray(automatic_img.dataobj) + edited = np.zeros(expected_shape, dtype=np.uint8) + edited[automatic == FORNIX_LABEL] = FORNIX_LABEL + edited[manual_cc] = CC_LABEL + return edited, False diff --git a/CorpusCallosum/utils/mapping_helpers.py b/CorpusCallosum/utils/mapping_helpers.py index c043cc542..197de433e 100644 --- a/CorpusCallosum/utils/mapping_helpers.py +++ b/CorpusCallosum/utils/mapping_helpers.py @@ -6,7 +6,7 @@ import SimpleITK as sitk from scipy.ndimage import affine_transform -from CorpusCallosum.data.constants import CC_LABEL, FORNIX_LABEL +from CorpusCallosum.data.constants import CC_LABEL, FORNIX_LABEL, SUBSEGMENT_LABELS from CorpusCallosum.utils.types import Polygon3dType from FastSurferCNN.utils import ( AffineMatrix3x3, @@ -396,3 +396,67 @@ def _map_softlabel_to_orig(data: Image3d, fill: int) -> Image3d: orig_space_segmentation_path, ) return seg_orig_space + + +def map_hard_segmentation_to_orig( + edited_segmentation: Image3d, + reference_image: nibabelImage, + orig2slab_vox2vox: AffineMatrix4x4, + cc_subseg_midslice: Image2d, + orig2midslice_vox2vox: AffineMatrix4x4, + output_path: str | Path, + automatic_orig_segmentation: nibabelImage | None = None, +) -> np.ndarray[Shape3d, np.dtype[np.int_]]: + """Map an edited upright segmentation to original space. + + Subdivision labels are regenerated from the edited middle-slice CC contour. + The fornix is copied from ``automatic_orig_segmentation`` when provided; + otherwise label ``FORNIX_LABEL`` is mapped from ``edited_segmentation``. + """ + orig_shape = reference_image.shape[:3] + edited_cc_orig = affine_transform( + np.equal(edited_segmentation, CC_LABEL).astype(np.uint8), + orig2slab_vox2vox, + output_shape=orig_shape, + order=0, + mode="constant", + cval=0, + prefilter=False, + ).astype(bool) + cc_subseg_orig = affine_transform( + cc_subseg_midslice[None], + orig2midslice_vox2vox, + output_shape=orig_shape, + order=0, + mode="nearest", + prefilter=False, + ) + + output = np.zeros(orig_shape, dtype=np.uint8) + if automatic_orig_segmentation is None: + edited_fornix_orig = affine_transform( + np.equal(edited_segmentation, FORNIX_LABEL).astype(np.uint8), + orig2slab_vox2vox, + output_shape=orig_shape, + order=0, + mode="constant", + cval=0, + prefilter=False, + ).astype(bool) + output[edited_fornix_orig] = FORNIX_LABEL + else: + automatic_orig = np.asarray(automatic_orig_segmentation.dataobj) + output[automatic_orig == FORNIX_LABEL] = FORNIX_LABEL + valid_subsegments = np.isin(cc_subseg_orig, SUBSEGMENT_LABELS) + if np.any(np.logical_and(edited_cc_orig, ~valid_subsegments)): + raise ValueError("Could not assign CC subdivision labels to all edited CC voxels.") + output[edited_cc_orig] = cc_subseg_orig[edited_cc_orig].astype(np.uint8) + + output_path = Path(output_path) + output_path.parent.mkdir(exist_ok=True, parents=True) + logger.info(f"Saving edited segmentation in original space to {output_path}") + nib.save( + nib.MGHImage(output, reference_image.affine, reference_image.header), + output_path, + ) + return output diff --git a/doc/overview/EDITING.md b/doc/overview/EDITING.md index 25244bf78..06998e654 100644 --- a/doc/overview/EDITING.md +++ b/doc/overview/EDITING.md @@ -37,9 +37,10 @@ Possible Edits FastSurfer supports the following edits: 1. [Bias field corrected inputs](#bias-field-correction) (for improved image quality, not really an edit) 2. [asegdkt_segfile](#asegdkt_segfile): `/mri/aparc.DKTatlas+aseg.deep.mgz` via `/mri/aparc.DKTatlas+aseg.deep.manedit.mgz` -3. [Talairach registration](#talairach-registration): `/mri/transforms/talairach.xfm` (overwrites automatic results from `/mri/transforms/talairach.auto.xfm`) -4. [White matter segmentation](#white-matter-segmentation): `/mri/wm.mgz` and `/mri/filled.mgz` -5. [Pial placement](#pial-surface-placement): `/mri/brain.finalsurfs.mgz` via `/mri/brain.finalsurfs.manedit.mgz` +3. [Corpus callosum segmentation](#corpus-callosum-segmentation): `/mri/callosum.CC.upright.mgz` via `/mri/callosum.CC.upright.manedit.mgz` +4. [Talairach registration](#talairach-registration): `/mri/transforms/talairach.xfm` (overwrites automatic results from `/mri/transforms/talairach.auto.xfm`) +5. [White matter segmentation](#white-matter-segmentation): `/mri/wm.mgz` and `/mri/filled.mgz` +6. [Pial placement](#pial-surface-placement): `/mri/brain.finalsurfs.mgz` via `/mri/brain.finalsurfs.manedit.mgz` Note, as FastSurfer's surface pipeline is derived from FreeSurfer, some editing options and corresponding naming schemes are inherited from FreeSurfer. @@ -131,6 +132,71 @@ In specific, such errors are inspected in `/mri/aparc.DKTatlas+aseg 2. Open `/mri/aparc.DKTatlas+aseg.deep.manedit.mgz` (for example using Freeview) and resolve all errors/quality issues. 3. [Re-run FastSurfer](#general-process) to propagate the changes into other results. Among others, this updates `/mri/aseg.auto_noCCseg.mgz` and `/mri/mask.mgz`. +Corpus callosum segmentation +---------------------------- + +### When to use +The corpus callosum mask is over- or under-segmented, or the error affects the CC contours, thickness, subdivisions, +surface, QC images, or volumetric inpainting. + +### What to do +1. Run FastSurfer with `--qc_snap`. This creates `mri/upright_volume.mgz`, the intensity reference in the same space as + the editable CC segmentation. For example: + + ```bash + ./run_fastsurfer.sh \ + --t1 /data/input/sub001_T1w.nii.gz \ + --sd /data/subjects \ + --sid sub001 \ + --seg_only \ + --qc_snap + ``` + +2. Copy `/mri/callosum.CC.upright.mgz` to + `/mri/callosum.CC.upright.manedit.mgz`: + + ```bash + SUBJECT_DIR=/data/subjects/sub001 + cp "$SUBJECT_DIR/mri/callosum.CC.upright.mgz" \ + "$SUBJECT_DIR/mri/callosum.CC.upright.manedit.mgz" + ``` + +3. Edit label 192 in the manedit file using `mri/upright_volume.mgz` as the reference image. For example: + + ```bash + freeview \ + -v "$SUBJECT_DIR/mri/upright_volume.mgz" \ + -v "$SUBJECT_DIR/mri/callosum.CC.upright.manedit.mgz":colormap=lut:opacity=0.5 + ``` + + A manual file containing any fornix label-250 voxels is self-contained and authoritative for the entire fornix. To + retain the automatic fornix instead, remove all label-250 voxels from the manual file; the automatic upright and + original-space CC segmentations must then be present. + +4. Re-run the original FastSurfer command with `--edits` and otherwise identical options: + + ```bash + ./run_fastsurfer.sh \ + --t1 /data/input/sub001_T1w.nii.gz \ + --sd /data/subjects \ + --sid sub001 \ + --seg_only \ + --qc_snap \ + --edits + ``` + +FastSurfer recomputes CC contours, morphometry, subdivisions, surfaces, statistics, and QC outputs. It creates +`mri/callosum.CC.orig.manedit.mgz` and uses that derived volume for downstream CC inpainting. The automatic +`callosum.CC.upright.mgz` and `callosum.CC.orig.mgz` files remain unchanged. Do not edit the original-space file +directly. If an upright manedit exists but `--edits` is omitted, FastSurfer exits with an error rather than silently +ignoring the correction. + +If the upright segmentation was created with supplied AC/PC coordinates through `fastsurfer_cc.py`, pass exactly the +same coordinates on the edit rerun so the edited slab geometry remains valid. + +For direct expert usage, including Docker and Singularity/Apptainer commands, see the +[FastSurfer-CC expert documentation](../scripts/fastsurfer_cc.rst). + Talairach registration ---------------------- diff --git a/doc/overview/OUTPUT_FILES.md b/doc/overview/OUTPUT_FILES.md index e7c32a4f8..a599fb7f7 100644 --- a/doc/overview/OUTPUT_FILES.md +++ b/doc/overview/OUTPUT_FILES.md @@ -23,20 +23,21 @@ The Corpus Callosum module outputs the files in the table shown below. It create |:----------------|--------------------------------|--------|--------------------------------------------------------------------------------------------------------------| | mri | callosum.CC.upright.mgz | cc | corpus callosum segmentation in upright space | | mri | callosum.CC.orig.mgz | cc | corpus callosum segmentation in conformed image orientation | +| mri | callosum.CC.orig.manedit.mgz | cc | corrected CC and supplied or retained fornix mapped to conformed space for downstream inpainting | | mri | callosum.CC.soft.mgz | cc | corpus callosum soft labels (in upright space) | | mri | fornix.CC.soft.mgz | cc | fornix soft labels (in upright space) | | mri | background.CC.soft.mgz | cc | background soft labels (in upright space) | -| mri | upright_volume.mgz | cc | conformed image mapped to upright space (only with fastsurfer_cc.py `--upright_volume`) | -| mri/transforms | cc_up.lta | cc | transform from conformed to upright space | +| mri | upright_volume.mgz | cc | conformed image mapped to upright space (with `--qc_snap` or fastsurfer_cc.py `--upright_volume`) | +| mri/transforms | cc_up.lta | cc | transform from conformed to upright space | | mri/transforms | orient_volume.lta | cc | transform to standardized space | -| stats | callosum.CC.midslice.json | cc | measurements from the mid-sagittal slice (landmarks, area, thickness, etc.) | -| stats | callosum.CC.all_slices.json | cc | comprehensive per-slice analysis | -| qc_snapshots | callosum.png | cc | debug visualization of CC contours, AC, PC and thickness (only with run_fastsurfer.sh `--qc_snap`) | -| qc_snapshots | callosum.thickness.png | cc | 3D thickness visualization (only with run_fastsurfer.sh `--qc_snap`) | -| qc_snapshots | corpus_callosum.html | cc | interactive 3D mesh visualization (only with run_fastsurfer.sh `--qc_snap`) | -| surf | callosum.surf | cc | 3D Corpus Callosum mesh in FreeSurfer surface format (open with freeview) | -| surf | callosum.thickness.w | cc | FreeSurfer overlay file containing thickness values (open with callosum.surf in freeview) | -| surf | callosum.vtk | cc | VTK format mesh file for 3D visualization | +| stats | callosum.CC.midslice.json | cc | measurements from the mid-sagittal slice (landmarks, area, thickness, etc.) | +| stats | callosum.CC.all_slices.json | cc | comprehensive per-slice analysis | +| qc_snapshots | callosum.png | cc | debug visualization of CC contours, AC, PC and thickness (only with run_fastsurfer.sh `--qc_snap`) | +| qc_snapshots | callosum.thickness.png | cc | 3D thickness visualization (only with run_fastsurfer.sh `--qc_snap`) | +| qc_snapshots | corpus_callosum.html | cc | interactive 3D mesh visualization (only with run_fastsurfer.sh `--qc_snap`) | +| surf | callosum.surf | cc | 3D Corpus Callosum mesh in FreeSurfer surface format (open with freeview) | +| surf | callosum.thickness.w | cc | FreeSurfer overlay file containing thickness values (open with callosum.surf in freeview) | +| surf | callosum.vtk | cc | VTK format mesh file for 3D visualization | CerebNet module --------------- diff --git a/doc/overview/modules/CC.md b/doc/overview/modules/CC.md index b08f4f5dc..b4a8d0ae7 100644 --- a/doc/overview/modules/CC.md +++ b/doc/overview/modules/CC.md @@ -62,6 +62,8 @@ All anatomical landmarks are given image voxel coordinates (LIA orientation) - `pc_center_oriented_volume`: PC coordinates in standardized space (orient_volume.lta) - `ac_center_upright`: AC coordinates in upright space (cc_up.lta) - `pc_center_upright`: PC coordinates in upright space (cc_up.lta) +- `landmark_source`: Whether AC/PC were obtained from the model or supplied by the user +- `segmentation_source`: Whether morphometry used the automatic or manual upright CC segmentation ### `stats/callosum.CC.all_slices.json` (Multi-Slice Analysis) This file contains comprehensive per-slice analysis when using `--slice_selection all`: diff --git a/doc/scripts/fastsurfer_cc.rst b/doc/scripts/fastsurfer_cc.rst index a30a01d09..17bc246f4 100644 --- a/doc/scripts/fastsurfer_cc.rst +++ b/doc/scripts/fastsurfer_cc.rst @@ -23,7 +23,6 @@ The following section provides a detailed overview of the command-line interface :prog: fastsurfer_cc.py - Midplane extraction ------------------- The ``--midplane_method`` flag controls how the corpus callosum pipeline refines the midsagittal plane before segmentation. This is implemented in ``CorpusCallosum/registration/midsagittal_plane_alignment.py``. @@ -38,6 +37,138 @@ Available modes are: The refinement after fsaverage alignment is intentionally conservative, expected to only make adjustments for unusual anatomies, or significant asymmetry. +Supplying AC/PC landmarks +~~~~~~~~~~~~~~~~~~~~~~~~~ +The expert interface accepts paired ``--ac_coords X Y Z`` and ``--pc_coords X Y Z`` arguments. Values are floating-point +voxel coordinates in ``orig.mgz`` space. When supplied, the landmark network is skipped and the points are used for +segmentation conditioning, morphometry, QC, orientation transforms, and output measurements. + +The selected midsagittal plane is minimally rotated and translated so that it contains both 3D landmarks exactly. A +plane adjustment larger than 15 degrees emits a warning because this commonly indicates a coordinate-space error. + +For example: + +.. code-block:: bash + + python3 CorpusCallosum/fastsurfer_cc.py \ + --sd /data/subjects \ + --sid sub001 \ + --ac_coords 127.4 126.8 128.1 \ + --pc_coords 127.9 103.6 128.7 \ + --upright_volume mri/upright_volume.mgz + +For Docker or Singularity/Apptainer, append the same two options to the corresponding expert command. Coordinates are +interpreted in the voxel space of the subject's input ``mri/orig.mgz``, not scanner RAS coordinates. Both points are +required, must be finite and distinct, and must lie inside the image. A left-right AC-PC line is rejected because it +does not define a stable sagittal plane. These options are available through the expert interface and are not exposed +by ``run_fastsurfer.sh``. + +Manual CC edits +--------------- +.. warning:: + Choose and record any supplied AC/PC coordinates before creating a manual CC correction. The automatic run and + every edit rerun must use exactly the same ``--ac_coords``, ``--pc_coords``, and ``--midplane_method`` values. If + these values change, regenerate ``mri/upright_volume.mgz`` and ``mri/callosum.CC.upright.mgz`` and recreate or + explicitly rebase the manual correction on the new upright image. Do not apply an edit created on a different + upright plane. Because upright files use a standardized affine, FastSurfer-CC cannot reliably detect every stale + edit from file geometry alone. + +To reprocess a manual CC correction, first copy the automatic upright segmentation: + +.. code-block:: bash + + SUBJECT_DIR=/data/subjects/sub001 + cp "$SUBJECT_DIR/mri/callosum.CC.upright.mgz" \ + "$SUBJECT_DIR/mri/callosum.CC.upright.manedit.mgz" + +Edit label 192 in ``mri/callosum.CC.upright.manedit.mgz`` using ``mri/upright_volume.mgz`` as the anatomical reference. +The manual file may also contain fornix label 250. Then rerun the expert command with the manual input: + +.. code-block:: bash + + python3 CorpusCallosum/fastsurfer_cc.py \ + --sd /data/subjects \ + --sid sub001 \ + --segmentation_manedit mri/callosum.CC.upright.manedit.mgz \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png + +If the manual segmentation contains any voxels with fornix label 250, FastSurfer-CC treats its CC and fornix as +self-contained input. If label 250 is absent, FastSurfer-CC instead retains the fornix from the automatic upright and +original-space CC segmentations, which must therefore exist from a previous run. It recomputes all CC-derived results +and writes ``mri/callosum.CC.orig.manedit.mgz``. The top-level ``run_fastsurfer.sh --edits`` workflow also uses this +derived file for downstream inpainting. Any label-250 voxels make the manual file authoritative for the entire +fornix; remove all label-250 voxels from the manual file to use the automatic fornix instead. Do not manually edit the +original-space file. + +.. note:: + A direct expert invocation produces the CC outputs, including ``callosum.CC.orig.manedit.mgz``, but does not paint + them into the broader aseg outputs. Use the top-level ``run_fastsurfer.sh --edits`` workflow when that downstream + integration is required. The top-level workflow does not currently accept supplied AC/PC coordinates. + +The complete Docker edit rerun is: + +.. code-block:: bash + + SUBJECTS_DIR=/data/fastsurfer + SID=sub001 + + docker run --gpus all --rm \ + --user "$(id -u):$(id -g)" \ + --volume "$SUBJECTS_DIR:/output" \ + --entrypoint /fastsurfer/tools/Docker/entrypoint.sh \ + deepmi/fastsurfer:latest \ + python3 /fastsurfer/CorpusCallosum/fastsurfer_cc.py \ + --sd /output \ + --sid "$SID" \ + --segmentation_manedit mri/callosum.CC.upright.manedit.mgz \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png + +The example above resolves the relative manual path inside ``/output/$SID``. A manual segmentation stored elsewhere +on the host can instead be mounted read-only and passed by its absolute container path: + +.. code-block:: bash + + MANUAL_EDIT=/data/annotations/sub001_cc_manedit.mgz + + docker run --gpus all --rm \ + --user "$(id -u):$(id -g)" \ + --volume "$SUBJECTS_DIR:/output" \ + --volume "$MANUAL_EDIT:/manual-edit.mgz:ro" \ + --entrypoint /fastsurfer/tools/Docker/entrypoint.sh \ + deepmi/fastsurfer:latest \ + python3 /fastsurfer/CorpusCallosum/fastsurfer_cc.py \ + --sd /output \ + --sid "$SID" \ + --segmentation_manedit /manual-edit.mgz \ + --qc_image qc_snapshots/callosum.png + +The complete Singularity edit rerun is: + +.. code-block:: bash + + SUBJECTS_DIR=/data/fastsurfer + SID=sub001 + FASTSURFER_SIF=/containers/fastsurfer-gpu.sif + + singularity exec --nv --no-mount home,cwd -e \ + --bind "$SUBJECTS_DIR:/output" \ + "$FASTSURFER_SIF" \ + python3 /fastsurfer/CorpusCallosum/fastsurfer_cc.py \ + --sd /output \ + --sid "$SID" \ + --segmentation_manedit mri/callosum.CC.upright.manedit.mgz \ + --upright_volume mri/upright_volume.mgz \ + --qc_image qc_snapshots/callosum.png \ + --thickness_image qc_snapshots/callosum.thickness.png + +An automatic CC run is only required when the manual segmentation does not contain fornix label 250. If the +segmentation used as the editing reference was generated with supplied AC/PC landmarks, append the exact same +``--ac_coords X Y Z --pc_coords X Y Z`` arguments to the edit command. + Quality Control diff --git a/run_fastsurfer.sh b/run_fastsurfer.sh index 45a36262e..feffe2b08 100755 --- a/run_fastsurfer.sh +++ b/run_fastsurfer.sh @@ -152,7 +152,10 @@ FLAGS: with the segmentation pipeline. --edits Enables manual edits by replacing select intermediate/ result files by manedit substitutes (*.manedit.). - Segmentation: and . + Segmentation edits (default paths): + mri/aparc.DKTatlas+aseg.deep.manedit.mgz + mri/mask.manedit.mgz + mri/callosum.CC.upright.manedit.mgz Surface: Disables check for existing recon-surf.sh run; edits of mri/wm.mgz and brain.finalsurfs.mgz as well as FreeSurfer-style WM control points. @@ -228,8 +231,9 @@ SEGMENTATION PIPELINE: CORPUS CALLOSUM MODULE: --no_cc Skip the segmentation and analysis of the corpus callosum. - --qc_snap Create QC snapshots in \$SUBJECTS_DIR/\$sid/qc_snapshots - to simplify the QC process. + --qc_snap Create quality control images in \$SUBJECTS_DIR/\$sid/qc_snapshots + to simplify the QC process. Also creates additional volumes + in mri/ for QC. HYPOTHALAMUS MODULE (HypVINN): --no_hypothal Skip the hypothalamus segmentation. @@ -525,7 +529,7 @@ case $key in --qc_snap) hypvinn_flags+=(--qc_snap) ; cc_flags+=(--qc_image "qc_snapshots/callosum.png" --thickness_image "qc_snapshots/callosum.thickness.png" - --cc_html "qc_snapshots/corpus_callosum.html") + --cc_html "qc_snapshots/corpus_callosum.html" --upright_volume "mri/upright_volume.mgz") ;; ############################################################## @@ -1291,8 +1295,31 @@ then asegdkt_withcc_segfile="$(add_file_suffix "$asegdkt_segfile" "withCC")" asegdkt_withcc_vinn_statsfile="$(add_file_suffix "$asegdkt_vinn_statsfile" "withCC")" aseg_auto_statsfile="$(dirname "$aseg_vinn_statsfile")/aseg.auto.mgz" - # note: callosum manedit currently only affects inpainting and not internal FastSurferCC processing (surfaces etc) + callosum_upright_seg="$subject_dir/mri/callosum.CC.upright.mgz" + callosum_upright_seg_manedit="$(add_file_suffix "$callosum_upright_seg" "manedit")" callosum_seg_manedit="$(add_file_suffix "$callosum_seg" "manedit")" + if [[ -f "$callosum_seg_manedit" ]] && [[ ! -f "$callosum_upright_seg_manedit" ]] + then + { + echo "ERROR: Legacy original-space CC edit $callosum_seg_manedit detected without" + echo " $callosum_upright_seg_manedit. Edit the upright segmentation instead; the" + echo " original-space manedit file is generated automatically during the edit rerun." + } | tee -a "$seg_log" + exit 1 + fi + if [[ -f "$callosum_upright_seg_manedit" ]] + then + if [[ "$edits" == "true" ]] + then + cc_flags+=(--segmentation_manedit "$callosum_upright_seg_manedit") + else + { + echo "ERROR: $callosum_upright_seg_manedit (manedit file for the upright CC segmentation) detected," + echo " but --edits was not passed. Delete the manedit file or add --edits." + } | tee -a "$seg_log" + exit 1 + fi + fi # generate callosum segmentation, mesh, shape and downstream measure files cmd=($python "$CorpusCallosumDir/fastsurfer_cc.py" --sd "$sd" --sid "$subject" "--threads" "$threads_seg" "--conformed_name" "$conformed_name" "--aseg_name" "$aseg_segfile" @@ -1304,7 +1331,10 @@ then echo "ERROR: FastSurferCC corpus callosum analysis failed!" | tee -a "$seg_log" exit "$exit_code" fi - if [[ "$edits" == "true" ]] && [[ -f "$callosum_seg_manedit" ]] ; then callosum_seg="$callosum_seg_manedit" ; fi + if [[ "$edits" == "true" ]] && [[ -f "$callosum_upright_seg_manedit" ]] + then + callosum_seg="$callosum_seg_manedit" + fi { # add CC into aparc.DKTatlas+aseg.deep.mgz and aseg.auto.mgz as mri_cc did before. cmd=($python "$CorpusCallosumDir/paint_cc_into_pred.py" -in_cc "$callosum_seg" -in_pred "$asegdkt_segfile"