Cc subsegment ordering

CorpusCallosum/fastsurfer_cc.py

@@ -834,7 +834,13 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4:
 
     logger.info(f"Processing slices with selection mode: {slice_selection}")
     try:
-        slice_results, slice_io_futures, cc_contours, num_failed_slices = recon_cc_surf_measures_multi(
+        (
+            slice_results,
+            slice_io_futures,
+            cc_contours,
+            num_failed_slices,
+            result_slice_indices,
+        ) = recon_cc_surf_measures_multi(
             segmentation=cc_fn_seg_labels,
             slice_selection=slice_selection,
             orig_header=orig.header,
@@ -859,11 +865,16 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4:
         slice_io_futures = []
         cc_contours = []
         num_failed_slices = cc_fn_seg_labels.shape[0] if slice_selection == "all" else 1
+        if slice_selection == "middle":
+            result_slice_indices = [cc_fn_seg_labels.shape[0] // 2]
+        elif slice_selection == "all":
+            result_slice_indices = list(range(cc_fn_seg_labels.shape[0]))
+        else:
+            result_slice_indices = [int(slice_selection)]
 
     # Filter out None results for further processing
     valid_slice_results = [r for r in slice_results if r is not None]
     valid_cc_contours = [c for c in cc_contours if c is not None]
-    num_failed_slices += len(cc_contours) - len(valid_cc_contours)
     outer_contours = []
     cc_volume_contour = None
 
@@ -878,14 +889,35 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4:
             )
         outer_contours = [slice_result["split_contours"][0] for slice_result in valid_slice_results]
 
-        if len(outer_contours) > 1 and not check_area_changes(outer_contours):
+        if num_failed_slices == 0 and len(outer_contours) > 1 and not check_area_changes(outer_contours):
             logger.warning(
                 "Large area changes detected between consecutive slices, this is likely due to a segmentation error."
             )
 
-    # Get middle slice result if available
-    middle_slice_idx = len(slice_results) // 2
-    middle_slice_result = slice_results[middle_slice_idx] if middle_slice_idx < len(slice_results) else None
+    def _select_middle_valid_slice():
+        """Select the center slice result, or the nearest valid result."""
+        if not slice_results:
+            return None, None, None
+        target_idx = len(slice_results) // 2
+        if slice_results[target_idx] is not None:
+            return target_idx, result_slice_indices[target_idx], slice_results[target_idx]
+
+        valid_results = [
+            (abs(i - target_idx), i, result)
+            for i, result in enumerate(slice_results)
+            if result is not None
+        ]
+        if not valid_results:
+            return None, None, None
+
+        _, nearest_idx, nearest_result = min(valid_results)
+        logger.warning(
+            f"Middle slice morphometry failed; using nearest valid slice result at position {nearest_idx + 1} "
+            f"of {len(slice_results)}."
+        )
+        return nearest_idx, result_slice_indices[nearest_idx], nearest_result
+
+    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)
     if sd.has_attribute("cc_orig_segfile"):
@@ -944,7 +976,7 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4:
     voxel_volume = np.prod(vox_size)
     additional_metrics["voxel_volume"] = voxel_volume
 
-    if len(valid_cc_contours) > 1:
+    if num_failed_slices == 0 and len(valid_cc_contours) > 1:
         logger.info(
             f"CC voxel count: {cc_num_voxel} at {voxel_volume:.2f} mm^3 voxel volume => "
             f"{cc_num_voxel * voxel_volume:.2f} mm^3"
@@ -980,6 +1012,13 @@ def _orig2midslab_vox2vox(additional_context: int = 0) -> AffineMatrix4x4:
     additional_metrics["num_thickness_points"] = num_thickness_points
     additional_metrics["subdivision_method"] = subdivision_method
     additional_metrics["subdivision_ratios"] = subdivisions
+    additional_metrics["selected_morphometry_slice"] = selected_slice_idx
+    additional_metrics["selected_morphometry_slice_position"] = selected_slice_position
+    additional_metrics["selected_morphometry_slice_is_fallback"] = bool(
+        selected_slice_position is not None
+        and slice_results
+        and selected_slice_position != len(slice_results) // 2
+    )
     additional_metrics["contour_smoothing"] = contour_smoothing
     additional_metrics["slice_selection"] = slice_selection
     additional_metrics["midline_refine_shift_vox"] = float(midline_shift_vox)

CorpusCallosum/shape/postprocessing.py

@@ -97,7 +97,13 @@ def recon_cc_surf_measures_multi(
     subdivision_method: SubdivisionMethod,
     contour_smoothing: int,
     subject_dir: SubjectDirectory,
-) -> tuple[list[CCMeasuresDict], list[concurrent.futures.Future], list[CCContour]]:
+) -> tuple[
+    list[CCMeasuresDict | None],
+    list[concurrent.futures.Future],
+    list[CCContour],
+    int,
+    list[int],
+]:
     """Surface reconstruction and metrics computation of corpus callosum slices based on selection mode.
 
     Parameters
@@ -134,16 +140,19 @@ def recon_cc_surf_measures_multi(
 
     Returns
     -------
-    list of CCMeasuresDict
-        List of slice processing results.
+    list of CCMeasuresDict or None
+        List of slice processing results in requested slice order. Failed slices
+        are represented by None.
     list of concurrent.futures.Future
         List of background IO processes.
     list of CCContour
         List of CC contours.
     int
         Number of failed slices.
+    list of int
+        Slice indices requested for reconstruction, in result order.
     """
-    slice_cc_measures: list[CCMeasuresDict] = []
+    slice_cc_measures: list[CCMeasuresDict | None] = []
     io_futures = []
 
     if subdivision_method == "angular" and not np.allclose(np.diff(subdivisions), np.diff(subdivisions)[0]):
@@ -172,7 +181,7 @@ def recon_cc_surf_measures_multi(
         num_slices = segmentation.shape[0]
         start_slice = 0
         end_slice = segmentation.shape[0]
-        slices_to_recon = range(start_slice, end_slice)
+        slices_to_recon = list(range(start_slice, end_slice))
     else:  # specific slice number
         num_slices = 1
         slices_to_recon = [int(slice_selection)]
@@ -184,43 +193,40 @@ def _gen_slice2slab_vox2vox(_slice_idx: int) -> AffineMatrix4x4:
     fsavg_midslab_vox2ras = fsavg_vox2ras @ np.linalg.inv(fsavg2midslab_vox2vox)
     per_slice_vox2ras = fsavg_midslab_vox2ras @ np.stack(list(map(_gen_slice2slab_vox2vox, slices_to_recon)), axis=0)
 
-    per_slice_recon = process_executor().map(_each_slice, slices_to_recon, per_slice_vox2ras, chunksize=1)
-    cc_contours = []
+    slice_cc_measures = [None] * num_slices
+    cc_contours_by_slice: list[CCContour | None] = [None] * num_slices
 
     run = thread_executor().submit
     wants_output = subject_dir.has_attribute
     output_path = subject_dir.filename_by_attribute
 
-    def _zip_failed(it_idx, it_affine, it_result):
-        """Zip slice indices, affines, and results, logging errors for failed slices."""
-        _sentinel = object()
-        for idx, affine in zip(it_idx, it_affine, strict=True):
-            try:
-                result = next(it_result, _sentinel)
-            except Exception as e:
-                logger.error(f"Processing slice {idx} failed: {e}")
-                logger.exception(e)
-                yield idx, affine, (None, None)
-                continue
-            if result is _sentinel:
-                logger.error("Number of items in idx and affine did not match results")
-                return
-            yield idx, affine, result
-
-    slice_iterator = _zip_failed(slices_to_recon, per_slice_vox2ras, iter(per_slice_recon))
-    for i, (slice_idx, this_slice_vox2ras, _results) in enumerate(slice_iterator):
-        progress = f" ({i+1} of {num_slices})" if num_slices > 1 else ""
-        # unpack values from _results
-        cc_measures: CCMeasuresDict | None = _results[0]
-        _contour: CCContour | None = _results[1]
+    executor = process_executor()
+    recon_futures = {
+        executor.submit(_each_slice, slice_idx, this_slice_vox2ras): (
+            i, slice_idx, this_slice_vox2ras,
+        )
+        for i, (slice_idx, this_slice_vox2ras) in enumerate(
+            zip(slices_to_recon, per_slice_vox2ras, strict=True)
+        )
+    }
+    for completion_idx, future in enumerate(concurrent.futures.as_completed(recon_futures), start=1):
+        i, slice_idx, this_slice_vox2ras = recon_futures[future]
+        progress = f" ({completion_idx} of {num_slices})" if num_slices > 1 else ""
+        try:
+            cc_measures, _contour = future.result()
+        except Exception as e:
+            logger.error(f"Processing slice {slice_idx} failed: {e}")
+            logger.exception(e)
+            cc_measures = None
+            _contour = None
 
         if cc_measures is None or _contour is None:
             logger.warning(f"Calculating CC measurements for slice {slice_idx+1}{progress} failed")
             continue
 
         logger.info(f"Calculating CC measurements for slice {slice_idx+1}{progress}")
-        cc_contours.append(_contour)
-        slice_cc_measures.append(cc_measures)
+        cc_contours_by_slice[i] = _contour
+        slice_cc_measures[i] = cc_measures
         is_debug = logger.getEffectiveLevel() <= logging.DEBUG
         is_midslice = i == num_slices // 2
         if wants_output("cc_qc_image") and (is_debug or is_midslice):
@@ -249,17 +255,21 @@ def _zip_failed(it_idx, it_affine, it_result):
                 )
             )
 
+    cc_contours = [contour for contour in cc_contours_by_slice if contour is not None]
+
     if wants_output("save_template_dir"):
         template_dir = output_path("save_template_dir")
         # ensure directory exists
         template_dir.mkdir(parents=True, exist_ok=True)
         logger.info("Saving template files (contours.txt, thickness_values.txt, "
                     f"thickness_measurement_points.txt) to {template_dir}")
-        for j in range(len(cc_contours)):
-            io_futures.append(run(cc_contours[j].save_contour, template_dir / f"contour_{j}.txt"))
-            io_futures.append(run(cc_contours[j].save_thickness_values, template_dir / f"thickness_values_{j}.txt"))
+        for j, contour in enumerate(cc_contours_by_slice):
+            if contour is None:
+                continue
+            io_futures.append(run(contour.save_contour, template_dir / f"contour_{j}.txt"))
+            io_futures.append(run(contour.save_thickness_values, template_dir / f"thickness_values_{j}.txt"))
 
-    num_failed_slices = num_slices - len(cc_contours)
+    num_failed_slices = slice_cc_measures.count(None)
 
     mesh_outputs = ("html", "mesh", "thickness_overlay", "surf", "thickness_image")
     if len(cc_contours) > 1 and any(wants_output(f"cc_{n}") for n in mesh_outputs):
@@ -269,7 +279,7 @@ def _zip_failed(it_idx, it_affine, it_result):
                 f"{num_failed_slices} of {num_slices} requested slices. "
                 "Surface generation requires a complete slice set to preserve the physical slab spacing.",
             )
-            return slice_cc_measures, io_futures, cc_contours, num_failed_slices
+            return slice_cc_measures, io_futures, cc_contours, num_failed_slices, slices_to_recon
 
         _cc_contours = thread_executor().map(_resample_thickness, cc_contours)
         # Surface vertices represent the continuous analysis slab. The
@@ -320,7 +330,7 @@ def _zip_failed(it_idx, it_affine, it_result):
             logger.error("Error: No valid slices were found for postprocessing")
             raise ValueError("No valid slices were found for postprocessing")
 
-    return slice_cc_measures, io_futures, cc_contours, num_failed_slices
+    return slice_cc_measures, io_futures, cc_contours, num_failed_slices, slices_to_recon
 
 
 def _resample_thickness(contour: CCContour) -> CCContour:

CorpusCallosum/shape/thickness.py

@@ -333,6 +333,15 @@ def cc_thickness(
             # If midline computation fails, we cannot proceed with thickness calculation
             raise RuntimeError(f"Corpus callosum midline computation failed: {e}") from e
 
+        # LaPy orders open level paths by internal graph traversal, not by
+        # anatomical endpoint semantics. Downstream subdivision expects the
+        # midline to run from anterior to posterior.
+        anterior_endpoint = contour_2d[anterior_endpoint_idx]
+        start_dist_to_anterior = np.linalg.norm(midline_equidistant_asz[0, :2] - anterior_endpoint)
+        end_dist_to_anterior = np.linalg.norm(midline_equidistant_asz[-1, :2] - anterior_endpoint)
+        if start_dist_to_anterior > end_dist_to_anterior:
+            midline_equidistant_asz = midline_equidistant_asz[::-1]
+
         midline_equidistant_contour_space: np.ndarray = midline_equidistant_asz[:, :2]
 
         gf = compute_rotated_f(tria_asz, vfunc)