diff --git a/CorpusCallosum/fastsurfer_cc.py b/CorpusCallosum/fastsurfer_cc.py index 409974f2..a704b69d 100644 --- a/CorpusCallosum/fastsurfer_cc.py +++ b/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) diff --git a/CorpusCallosum/shape/postprocessing.py b/CorpusCallosum/shape/postprocessing.py index ac9f05f0..5915f36a 100644 --- a/CorpusCallosum/shape/postprocessing.py +++ b/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: diff --git a/CorpusCallosum/shape/thickness.py b/CorpusCallosum/shape/thickness.py index e86eb74e..2624d3f9 100644 --- a/CorpusCallosum/shape/thickness.py +++ b/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)