diff --git a/CorpusCallosum/cc_visualization.py b/CorpusCallosum/cc_visualization.py index a793ba44..1f2f3beb 100644 --- a/CorpusCallosum/cc_visualization.py +++ b/CorpusCallosum/cc_visualization.py @@ -3,12 +3,15 @@ from pathlib import Path from typing import Literal +import nibabel as nib import numpy as np from CorpusCallosum.data.fsaverage_cc_template import load_fsaverage_cc_template from CorpusCallosum.shape.contour import CCContour from CorpusCallosum.shape.mesh import CCMesh +from FastSurferCNN.utils import AffineMatrix4x4 from FastSurferCNN.utils.logging import get_logger, setup_logging +from FastSurferCNN.utils.lta import read_lta logger = get_logger(__name__) @@ -39,7 +42,7 @@ def make_parser() -> argparse.ArgumentParser: "cc_mesh.vtk - VTK mesh file format " "cc_mesh.fssurf - FreeSurfer surface file " "cc_mesh_overlay.curv - FreeSurfer curvature overlay file " - "cc_mesh_snap.png - Screenshot/snapshot of the 3D mesh (requires whippersnappy>=1.3.1)", + "cc_mesh_snap.png - Screenshot/snapshot of the 3D mesh (requires whippersnappy>=2.1)", metavar="OUTPUT_DIR" ) parser.add_argument( @@ -213,6 +216,17 @@ def main( # 3D visualization cc_mesh = CCMesh.from_contours(contours, smooth=0) + if Path(output_dir / "mri" / "upright.mgz").exists(): + header = nib.load(output_dir / "mri" / "upright.mgz").header + # we need to get the upright image header, which is the same as cc_up.lta applied to orig. + elif Path(template_dir / "mri/orig.mgz").exists() and Path(template_dir / "mri/transforms/cc_up.lta").exists(): + image = nib.load(template_dir / "mri" / "orig.mgz") + lta_mat: AffineMatrix4x4 = read_lta(template_dir / "mri/transforms/cc_up.lta")["lta"] + image.affine = lta_mat @ image.affine + header = image.header + else: + header = None + plot_kwargs = dict( colormap=colormap, color_range=color_range, @@ -225,15 +239,16 @@ def main( logger.info(f"Writing vtk file to {output_dir / 'cc_mesh.vtk'}") cc_mesh.write_vtk(str(output_dir / "cc_mesh.vtk")) logger.info(f"Writing freesurfer surface file to {output_dir / 'cc_mesh.fssurf'}") - cc_mesh.write_fssurf(str(output_dir / "cc_mesh.fssurf")) + cc_mesh.write_fssurf(str(output_dir / "cc_mesh.fssurf"), image=header) logger.info(f"Writing freesurfer overlay file to {output_dir / 'cc_mesh_overlay.curv'}") cc_mesh.write_morph_data(str(output_dir / "cc_mesh_overlay.curv")) try: - cc_mesh.snap_cc_picture(str(output_dir / "cc_mesh_snap.png")) + cc_mesh.snap_cc_picture(str(output_dir / "cc_mesh_snap.png"), ref_header=header) logger.info(f"Writing 3D snapshot image to {output_dir / 'cc_mesh_snap.png'}") - except RuntimeError: - logger.warning("The cc_visualization script requires whippersnappy>=1.3.1 to makes screenshots, install with " - "`pip install whippersnappy>=1.3.1` !") + except Exception: + logger.warning("The cc_visualization script requires whippersnappy>=2.1 to makes screenshots, install with " + "`pip install whippersnappy>=2.1` !") + raise return 0 if __name__ == "__main__": diff --git a/CorpusCallosum/paint_cc_into_pred.py b/CorpusCallosum/paint_cc_into_pred.py index 39d0e61a..91272ed4 100644 --- a/CorpusCallosum/paint_cc_into_pred.py +++ b/CorpusCallosum/paint_cc_into_pred.py @@ -19,23 +19,20 @@ import sys from functools import partial from pathlib import Path -from typing import TypeVar, cast import nibabel as nib import numpy as np -from numpy import typing as npt from scipy import ndimage -import FastSurferCNN.utils.logging as logging from CorpusCallosum.data.constants import FORNIX_LABEL, SUBSEGMENT_LABELS from FastSurferCNN.data_loader.conform import is_conform +from FastSurferCNN.data_loader.data_utils import load_image from FastSurferCNN.reduce_to_aseg import reduce_to_aseg_and_save +from FastSurferCNN.utils import Mask2d, Mask3d, Shape3d, logging from FastSurferCNN.utils.arg_types import path_or_none from FastSurferCNN.utils.brainvolstats import mask_in_array from FastSurferCNN.utils.parallel import thread_executor -_T = TypeVar("_T", bound=np.number) - logger = logging.get_logger(__name__) HELPTEXT = """ @@ -55,7 +52,8 @@ Original Author: Leonie Henschel Date: Jul-10-2020 - +Modified by: Clemens Pollak, David Kügler +Date: Dec-2025 """ @@ -110,26 +108,23 @@ def make_parser() -> argparse.ArgumentParser: return parser -def paint_in_cc(pred: npt.NDArray[np.int_], - aseg_cc: npt.NDArray[np.int_]) -> npt.NDArray[np.int_]: +def paint_in_cc( + pred: np.ndarray[Shape3d, np.dtype[int]], + aseg_cc: np.ndarray[Shape3d, np.dtype[int]], +) -> np.ndarray[Shape3d, np.dtype[int]]: """Paint corpus callosum segmentation into aseg+dkt segmentation map. Parameters ---------- - pred : npt.NDArray[np.int_] + pred : np.ndarray Deep-learning segmentation map. - aseg_cc : npt.NDArray[np.int_] + aseg_cc : np.ndarray Aseg segmentation with CC. Returns ------- - npt.NDArray[np.int_] + np.ndarray Segmentation map with added CC. - - Notes - ----- - This function modifies the original array and does not create a copy. - The CC labels (251-255) from aseg_cc are copied into pred. """ cc_mask = mask_in_array(aseg_cc, SUBSEGMENT_LABELS) @@ -142,14 +137,14 @@ def paint_in_cc(pred: npt.NDArray[np.int_], logger.info(f"Painting CC: {np.sum(cc_mask)} voxels (replacing {num_wm_replaced} WM, " f"{num_background_replaced} background, {num_other_replaced} other)") - pred[cc_mask] = aseg_cc[cc_mask] - return pred - -def _fill_gaps_in_direction( - corrected_pred: npt.NDArray[np.int_], - potential_fill: npt.NDArray[np.bool_], - source_binary: npt.NDArray[np.bool_], - target_binary: npt.NDArray[np.bool_], + out = np.where(cc_mask, aseg_cc, pred) + return out + +def _fill_gaps_in_direction_( + corrected_pred: np.ndarray[Shape3d, np.dtype[int]], + potential_fill: Mask2d, + source_binary: Mask2d, + target_binary: Mask2d, x_slice: int, direction: str, max_gap_voxels: int, @@ -159,13 +154,13 @@ def _fill_gaps_in_direction( Parameters ---------- - corrected_pred : npt.NDArray[np.int_] + corrected_pred : np.ndarray The segmentation array to modify in place. - potential_fill : npt.NDArray[np.bool_] + potential_fill : np.ndarray 2D mask of potential fill regions for this slice. - source_binary : npt.NDArray[np.bool_] + source_binary : np.ndarray 2D binary mask of source structure (e.g., CC). - target_binary : npt.NDArray[np.bool_] + target_binary : np.ndarray 2D binary mask of target structure (e.g., ventricle). x_slice : int The x-coordinate of the current slice. @@ -254,10 +249,10 @@ def _fill_gaps_in_direction( return voxels_filled -def _fill_gaps_between_structures( - corrected_pred: npt.NDArray[np.int_], - source_mask: npt.NDArray[np.bool_], - target_mask: npt.NDArray[np.bool_], +def _fill_gaps_between_structures_( + corrected_pred: np.ndarray[Shape3d, np.dtype[int]], + source_mask: Mask3d, + target_mask: Mask3d, voxel_size: tuple[float, float, float], close_gap_size_mm: float, fillable_labels: set[int], @@ -267,11 +262,11 @@ def _fill_gaps_between_structures( Parameters ---------- - corrected_pred : npt.NDArray[np.int_] + corrected_pred : np.ndarray The segmentation array to modify in place. - source_mask : npt.NDArray[np.bool_] + source_mask : np.ndarray 3D binary mask of source structure (e.g., CC). - target_mask : npt.NDArray[np.bool_] + target_mask : np.ndarray 3D binary mask of target structure (e.g., ventricle or background). voxel_size : tuple[float, float, float] Voxel size in mm. @@ -315,13 +310,13 @@ def _fill_gaps_between_structures( potential_fill = (source_dilated & target_dilated) & ~(source_binary | target_binary) # Fill gaps in inferior-superior direction - voxels_filled += _fill_gaps_in_direction( + voxels_filled += _fill_gaps_in_direction_( corrected_pred, potential_fill, source_binary, target_binary, x, 'inferior-superior', max_gap_vox_inferior_superior, fillable_labels ) # Fill gaps in anterior-posterior direction - voxels_filled += _fill_gaps_in_direction( + voxels_filled += _fill_gaps_in_direction_( corrected_pred, potential_fill, source_binary, target_binary, x, 'anterior-posterior', max_gap_vox_anterior_posterior, fillable_labels ) @@ -333,11 +328,11 @@ def _fill_gaps_between_structures( def correct_wm_ventricles( - aseg_cc: npt.NDArray[np.int_], - fornix_mask: npt.NDArray[np.bool_], + aseg_cc: np.ndarray[Shape3d, np.dtype[int]], + fornix_mask: Mask3d, voxel_size: tuple[float, float, float], close_gap_size_mm: float = 3.0 -) -> npt.NDArray[np.int_]: +) -> np.ndarray[Shape3d, np.dtype[int]]: """Fill small gaps between corpus callosum, ventricles, and background. This function performs two gap-filling operations: @@ -349,9 +344,9 @@ def correct_wm_ventricles( Parameters ---------- - aseg_cc : npt.NDArray[np.int_] + aseg_cc : np.ndarray Aseg segmentation with CC already painted in. - fornix_mask : npt.NDArray[np.bool_] + fornix_mask : np.ndarray Mask of the fornix. Not currently used (kept for interface compatibility). voxel_size : tuple[float, float, float] Voxel size of the aseg image in mm. @@ -360,7 +355,7 @@ def correct_wm_ventricles( Returns ------- - npt.NDArray[np.int_] + np.ndarray Corrected segmentation map with filled gaps. """ # Create a copy to avoid modifying the original @@ -374,26 +369,28 @@ def correct_wm_ventricles( # Get background mask background_mask = aseg_cc == 0 - + print(np.unique(corrected_pred)) + # 1. Fill gaps between CC and ventricles (replace WM and background with ventricle labels) - _fill_gaps_between_structures( + _fill_gaps_between_structures_( corrected_pred, cc_mask, ventricle_mask, voxel_size, close_gap_size_mm, fillable_labels={0, 2, 41}, # background and WM description="between CC and ventricles (WM/background → ventricle)" ) - + print(np.unique(corrected_pred)) + # 2. Fill WM gaps between CC and background (replace WM with background) - _fill_gaps_between_structures( + _fill_gaps_between_structures_( corrected_pred, cc_mask, background_mask, voxel_size, close_gap_size_mm, fillable_labels={2, 41}, # only WM description="between CC and background (WM → background)" ) + print(np.unique(corrected_pred)) return corrected_pred if __name__ == "__main__": - from FastSurferCNN.utils import nibabelImage # Command Line options are error checking done here options = argument_parse() @@ -401,10 +398,9 @@ def correct_wm_ventricles( logging.setup_logging() logger.info(f"Reading inputs: {options.input_cc} {options.input_pred}...") - cc_seg_image = cast(nibabelImage, nib.load(options.input_cc)) - cc_seg_data = np.asanyarray(cc_seg_image.dataobj) - aseg_image = cast(nibabelImage, nib.load(options.input_pred)) - aseg_data = np.asanyarray(aseg_image.dataobj) + + tmap = thread_executor().map + (cc_seg_image, cc_seg_data), (aseg_image, aseg_data) = tmap(load_image, (options.input_cc, options.input_pred)) def _is_conform(img, dtype, verbose): return is_conform(img, vox_size=None, img_size=None, verbose=verbose, dtype=dtype) @@ -433,8 +429,8 @@ def _is_conform(img, dtype, verbose): initial_wm = np.sum((aseg_data == 2) | (aseg_data == 41)) initial_ventricles = np.sum((aseg_data == 4) | (aseg_data == 43)) - # Paint CC into prediction (modifies aseg_data in place) - paint_in_cc(aseg_data, cc_seg_data) + # Paint CC into prediction + aseg_data = paint_in_cc(aseg_data, cc_seg_data) # Apply ventricle gap filling corrections fornix_mask = cc_seg_data == FORNIX_LABEL diff --git a/CorpusCallosum/shape/mesh.py b/CorpusCallosum/shape/mesh.py index 8dcdfa53..d978b1d0 100644 --- a/CorpusCallosum/shape/mesh.py +++ b/CorpusCallosum/shape/mesh.py @@ -12,7 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -import tempfile from pathlib import Path from typing import TypeVar @@ -27,8 +26,8 @@ import FastSurferCNN.utils.logging as logging from CorpusCallosum.shape.contour import CCContour from CorpusCallosum.shape.thickness import make_mesh_from_contour -from FastSurferCNN.utils import AffineMatrix4x4, nibabelImage -from FastSurferCNN.utils.common import suppress_stdout, update_docstring +from FastSurferCNN.utils import AffineMatrix4x4, nibabelHeader, nibabelImage +from FastSurferCNN.utils.common import update_docstring try: from pyrr import Matrix44 @@ -478,11 +477,12 @@ def __create_cc_viewmat() -> "Matrix44": - -8 degrees around z-axis 3. Adds a small translation for better centering """ + from whippersnappy import ViewType, get_view_matrix if not HAS_PYRR: raise ImportError("Pyrr not installed, install pyrr with `pip install pyrr`.") - viewLeft = np.array([[0, 0, -1, 0], [-1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]]) # left w top up // right + viewLeft = get_view_matrix(ViewType.LEFT) # left w top up // right transl = Matrix44.from_translation((0, 0, 0.4)) viewmat = transl * viewLeft @@ -503,9 +503,7 @@ def __create_cc_viewmat() -> "Matrix44": def snap_cc_picture( self, output_path: Path | str, - fssurf_file: Path | str | None = None, - overlay_file: Path | str | None = None, - ref_image: Path | str | nibabelImage | None = None, + ref_header: Path | str | nibabelHeader | None = None, ) -> None: """Snap a picture of the corpus callosum mesh. @@ -513,34 +511,19 @@ def snap_cc_picture( ---------- output_path : Path, str Path where to save the snapshot image. - fssurf_file : Path, str, optional - Path to a FreeSurfer surface file to use for the snapshot. - If None, the mesh is saved to a temporary file. - overlay_file : Path, str, optional - Path to a FreeSurfer overlay file to use for the snapshot. - If None, the mesh is saved to a temporary file. - ref_image : Path, str, nibabelImage, optional - Path to reference image to use for tkr creation. If None, ignores the file for saving. + ref_header : Path, str, nibabelHeader, optional + Path to reference image header to use for tkr creation. If None, ignores the file for saving. Raises ------ - Warning + ImportError + If whippersnappy is not installed or if the version is too old. + ValueError If the mesh has no faces and cannot create a snapshot. - - Notes - ----- - The function: - 1. Creates temporary files for mesh and overlay data if needed. - 2. Uses whippersnappy to create a snapshot with: - - Custom view matrix for standard orientation. - - Ambient lighting and colorbar settings. - - Thickness overlay if available. - 3. Cleans up temporary files after use. """ + from packaging.version import parse try: - # Dummy import of OpenCL to ensure it's available for whippersnappy - import OpenGL.GL # noqa: F401 - from whippersnappy.core import snap1 + import whippersnappy except ImportError as e: # whippersnappy not installed raise ImportError( @@ -548,63 +531,50 @@ def snap_cc_picture( f"Please install {e.name}!", name=e.name, path=e.path ) from None - except Exception as e: - # Catch all other types of errors, - raise RuntimeError( - "Could not import OpenGL or whippersnappy. The snap_cc_picture method of CCMesh requires OpenGL and " - "whippersnappy to render the QC thickness image. On headless servers, this also requires a virtual " - "framebuffer like xvfb.", - ) from e - self.__make_parent_folder(output_path) + from nibabel.affines import apply_affine + + if parse(whippersnappy.__version__) < parse("2.1.0"): + raise ImportError( + f"The snap_cc_picture method of CCMesh requires whippersnappy>=2.1, but version " + f"{whippersnappy.__version__} is installed. Please upgrade whippersnappy to version 2.1 or higher!", + name="whippersnappy", path=None + ) # Skip snapshot if there are no faces if len(self.t) == 0: - logger.warning("Cannot create snapshot - no faces in mesh") - return + raise ValueError("Cannot create snapshot - no faces in mesh") - # create temp file - if fssurf_file: - fssurf_file = Path(fssurf_file) - else: - fssurf_file = tempfile.NamedTemporaryFile(suffix=".fssurf", delete=True).name - - ref_image_arg = str(ref_image) if isinstance(ref_image, (Path, str)) else ref_image - self.write_fssurf(fssurf_file, image=ref_image_arg) + self.__make_parent_folder(output_path) - if overlay_file: - overlay_file = Path(overlay_file) + if ref_header is None: + v = self.v else: - overlay_file = tempfile.NamedTemporaryFile(suffix=".w", delete=True).name - # Write thickness values in FreeSurfer '*.w' overlay format - self.write_morph_data(overlay_file) - - try: - with suppress_stdout(): - snap1( - fssurf_file, - overlaypath=overlay_file, - view=None, - viewmat=self.__create_cc_viewmat(), - width=3 * 500, - height=3 * 300, - outpath=output_path, - ambient=0.6, - colorbar_scale=0.5, - colorbar_y=0.88, - colorbar_x=0.19, - brain_scale=2.1, - fthresh=0, - caption="Corpus Callosum thickness (mm)", - caption_y=0.85, - caption_x=0.17, - caption_scale=0.5, - ) - except Exception as e: - raise e from None - - if fssurf_file and hasattr(fssurf_file, "close"): - fssurf_file.close() - if overlay_file and hasattr(overlay_file, "close"): - overlay_file.close() + from nibabel.freesurfer.mghformat import MGHHeader + + # if header is a file, load its header from the file + if isinstance(ref_header, (str, Path)): + ref_header = nib.load(ref_header).header + # if header is not already an MGHHeader, convert it to MGHHeader, so we have the get_vox2ras_tkr function + mgh_header = ref_header if isinstance(ref_header, MGHHeader) else MGHHeader.from_header(ref_header) + v = apply_affine(mgh_header.get_vox2ras_tkr(), self.v) + + whippersnappy.snap1( + mesh=(v, self.t), + overlay=self.mesh_vertex_colors, + view=self.__create_cc_viewmat(), + width=3 * 500, + height=3 * 300, + outpath=str(output_path), + ambient=0.6, + colorbar_scale=0.5, + colorbar_y=0.88, + colorbar_x=0.19, + scale=2.1, + fthresh=0, + caption="Corpus Callosum thickness (mm)", + caption_y=0.85, + caption_x=0.17, + caption_scale=0.5, + ) def smooth_(self, iterations: int = 1) -> None: """Smooth the mesh while preserving the z-coordinates. @@ -674,7 +644,7 @@ def to_vox_coordinates( return new_object @update_docstring(parent_doc=TriaMesh.write_fssurf.__doc__) - def write_fssurf(self, filename: Path | str, image: str | nibabelImage | None = None) -> None: + def write_fssurf(self, filename: Path | str, image: str | nibabelImage | nibabelHeader | None = None) -> None: """{parent_doc} Also creates parent directory if needed before writing the file. """ diff --git a/CorpusCallosum/shape/postprocessing.py b/CorpusCallosum/shape/postprocessing.py index e1a3df6c..b99abf21 100644 --- a/CorpusCallosum/shape/postprocessing.py +++ b/CorpusCallosum/shape/postprocessing.py @@ -272,49 +272,32 @@ def _zip_failed(it_idx, it_affine, it_result): logger.info(f"Saving vtk file to {vtk_file_path}") io_futures.append(run(cc_mesh.write_vtk, vtk_file_path)) - if wants_output("cc_thickness_overlay") and not wants_output("cc_thickness_image"): + if wants_output("cc_thickness_overlay"): overlay_file_path = output_path("cc_thickness_overlay") logger.info(f"Saving overlay file to {overlay_file_path}") io_futures.append(run(cc_mesh.write_morph_data, overlay_file_path)) if any(wants_output(f"cc_{n}") for n in ("thickness_image", "surf")): - import nibabel as nib - up_data: Image3d[np.uint8] = np.empty(upright_header["dims"][:3], dtype=upright_header.get_data_dtype()) - upright_img = nib.MGHImage(up_data, fsavg_vox2ras, upright_header) # the mesh is generated in upright coordinates, so we need to also transform to orig coordinates # Mesh is fsavg_midplane (RAS); we need to transform to voxel coordinates # fsavg ras is also on the midslice, so this is fine and we multiply in the IA and SP offsets - cc_mesh = cc_mesh.to_vox_coordinates(mesh_ras2vox=np.linalg.inv(fsavg_vox2ras @ orig2fsavg_vox2vox)) - cc_surf_generated = False + cc_mesh_orig = cc_mesh.to_vox_coordinates(mesh_ras2vox=np.linalg.inv(fsavg_vox2ras @ orig2fsavg_vox2vox)) + if wants_output("cc_surf"): + surf_file_path = output_path("cc_surf") + logger.info(f"Saving surf file to {surf_file_path}") + io_futures.append(run(cc_mesh_orig.write_fssurf, surf_file_path, image=upright_header)) + if wants_output("cc_thickness_image"): - # this will also write overlay and surface thickness_image_path = output_path("cc_thickness_image") logger.info(f"Saving thickness image to {thickness_image_path}") - kwargs = { - "fssurf_file": output_path("cc_surf") if wants_output("cc_surf") else None, - "overlay_file": output_path("cc_thickness_overlay") - if wants_output("cc_thickness_overlay") else None, - "ref_image": upright_img, - } try: - cc_mesh.snap_cc_picture(thickness_image_path, **kwargs) - cc_surf_generated = True - except (ImportError, ModuleNotFoundError) as e: - logger.error( - "The thickness image was not generated because whippersnappy, glfw or OpenGL are not installed." - ) - logger.exception(e) + cc_mesh_orig.snap_cc_picture(thickness_image_path, ref_header=upright_header) except Exception as e: logger.error( - "The thickness image was not generated (see below). On headless Linux systems or if the " - "x-server cannot/should not be accessed due to other reasons, xvfb-run may be used to provide " - "a virtual framebuffer for offscreen rendering." + "Generation of the thickness image failed (see below). Please ensure that whippersnappy and " + "(for headless rendering) EGL libraries (libegl1) are available." ) logger.exception(e) - if not cc_surf_generated and wants_output("cc_surf"): - surf_file_path = output_path("cc_surf") - logger.info(f"Saving surf file to {surf_file_path}") - io_futures.append(run(cc_mesh.write_fssurf, str(surf_file_path), image=upright_img)) if not slice_cc_measures: logger.error("Error: No valid slices were found for postprocessing") @@ -609,14 +592,13 @@ def make_subdivision_mask( # Use only as many labels as needed based on the number of subdivisions # Number of regions = number of division lines + 1 num_labels_needed = len(subdivision_lines) + 1 - cc_labels_posterior_to_anterior = SUBSEGMENT_LABELS[:num_labels_needed] + cc_labels_anterior_to_posterior = SUBSEGMENT_LABELS[:num_labels_needed][::-1] # Initialize with first segment label - subdivision_mask = np.full(slice_shape, cc_labels_posterior_to_anterior[0], dtype=np.int32) - + subdivision_mask = np.full(slice_shape, cc_labels_anterior_to_posterior[0], dtype=np.int32) # Process each subdivision line, subdivision_lines has for each division line the two points that are on the # contour and divide the subsegments - for label, segment_points in zip(cc_labels_posterior_to_anterior[1:], subdivision_lines, strict=True): + for label, segment_points in zip(cc_labels_anterior_to_posterior[1:], subdivision_lines, strict=True): # line_start and line_end are the intersection points of the CC subsegmentation boundary and the contour line line_start, line_end = segment_points @@ -634,14 +616,14 @@ def make_subdivision_mask( from FastSurferCNN.utils.plotting import backend with backend("qtagg"): plt.figure(figsize=(10, 8)) - plt.imshow(subdivision_mask, cmap='tab10') + plkwargs = {f"v{op}": getattr(np, op)(cc_labels_anterior_to_posterior) for op in ("min", "max")} + plt.imshow(subdivision_mask, cmap='tab10', **plkwargs) plt.colorbar(label='Subdivision') plt.title('CC Subdivision Mask') plt.xlabel('X') plt.ylabel('Y') plt.tight_layout() plt.show() - return subdivision_mask diff --git a/FastSurferCNN/reduce_to_aseg.py b/FastSurferCNN/reduce_to_aseg.py index bbf2c705..59b785ac 100644 --- a/FastSurferCNN/reduce_to_aseg.py +++ b/FastSurferCNN/reduce_to_aseg.py @@ -119,11 +119,9 @@ def reduce_to_aseg(data_inseg: np.ndarray[ShapeType, _TDType]) -> np.ndarray[Sha The reduced segmentation. """ LOGGER.info("Reducing to aseg ...") - # replace 2000... with 42 - data_inseg[data_inseg >= 2000] = 42 - # replace 1000... with 3 - data_inseg[data_inseg >= 1000] = 3 - return data_inseg + cortical_fill = np.full_like(data_inseg, 3) + cortical_fill[data_inseg >= 2000] = 42 + return np.where(data_inseg >= 1000, cortical_fill, data_inseg) def create_mask(aseg_data: np.ndarray[ShapeType, _TDType], dnum: int, enum: int) \ diff --git a/HypVINN/run_prediction.py b/HypVINN/run_prediction.py index 14486c28..95fbb2c8 100644 --- a/HypVINN/run_prediction.py +++ b/HypVINN/run_prediction.py @@ -31,7 +31,7 @@ load_checkpoint_config_defaults, ) from FastSurferCNN.utils.common import update_docstring -from FastSurferCNN.utils.parallel import get_num_threads, thread_executor +from FastSurferCNN.utils.parallel import get_num_threads, set_num_threads, thread_executor from HypVINN.config.hypvinn_files import HYPVINN_MASK_NAME, HYPVINN_SEG_NAME from HypVINN.data_loader.data_utils import hypo_map_label2subseg, rescale_image from HypVINN.inference import Inference @@ -155,6 +155,7 @@ def main( hypo_segfile: str = HYPVINN_SEG_NAME, hypo_maskfile: str = HYPVINN_MASK_NAME, qc_snapshots: bool = False, + threads: int | None = None, reg_mode: Literal["coreg", "robust", "none"] = "coreg", batch_size: int = 1, async_io: bool = False, @@ -192,6 +193,8 @@ def main( The name of the hypothalamus mask file. qc_snapshots : bool, default=False Whether to create QC snapshots. + threads : int, optional + If not None, updates the FastSurfer global setting in `FastSurfer.utils.parallel`. reg_mode : "coreg", "robust", "none", default="coreg" The registration mode to use. batch_size : int, default=1 @@ -208,6 +211,9 @@ def main( int, str 0, if successful, an error message describing the cause for the failure otherwise. """ + if threads is not None and threads > 1: + set_num_threads(threads) + from concurrent.futures import Future prep_tasks: dict[str, Future] = {} diff --git a/pyproject.toml b/pyproject.toml index 5b5bed33..4226a7a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -60,7 +60,7 @@ dependencies = [ [project.optional-dependencies] qc = [ - 'whippersnappy>=1.3.1', + 'whippersnappy>=2.1', ] doc = [ 'fastsurfer[qc]', diff --git a/requirements.mac.txt b/requirements.mac.txt index fc713c6a..acfd1534 100644 --- a/requirements.mac.txt +++ b/requirements.mac.txt @@ -19,5 +19,5 @@ yacs>=0.1.8 monai>=1.4.0 meshpy>=2025.1.1 pyrr>=0.10.3 -whippersnappy>=1.3.1 +whippersnappy>=2.1 pip>=25.0 \ No newline at end of file diff --git a/run_fastsurfer.sh b/run_fastsurfer.sh index 3bafea86..471790a3 100755 --- a/run_fastsurfer.sh +++ b/run_fastsurfer.sh @@ -764,32 +764,39 @@ then fi fi -maybe_xvfb=() -# check if we are running on a headless system (CC QC needs a (virtual) display that support OpenGL) +# Check if --thickness_image is in cc_flags and whippersnappy version is >= 2.1 if [[ "$run_seg_pipeline" == "true" ]] && [[ "$run_cc_module" == "true" ]] && \ - [[ "${cc_flags[*]}" =~ --thickness_image ]] + [[ "${cc_flags[*]}" == *"--thickness_image"* ]] then - # if we have xvfb-run, we can use it to provide a virtual display - if [[ -n "$(which xvfb-run)" ]] ; then maybe_xvfb=("xvfb-run" "-a") ; fi - - # try loading opengl, if this is successful we are fine - py_opengltest="import sys ; import glfw ; import whippersnappy.core ; sys.exit(1-glfw.init())" - opengl_error_message="$("${maybe_xvfb[@]}" $python -c "$py_opengltest" 2>&1 > /dev/null)" - exit_code="$?" - if [[ "$exit_code" != "0" ]] + # Check if whippersnappy is installed and version is >= 2.1 + whippersnappy_check=$($python -c " +try: + import whippersnappy as wspy + from packaging.version import parse + print('OK' if parse(wspy.__version__) >= parse('2.1') else ('OLD_VERSION:' + wspy.__version__)) +except ImportError: + print('NOT_INSTALLED') +except Exception as e: + print('ERROR:' + str(e)) +" 2>&1) + + if [[ "$whippersnappy_check" != "OK" ]] then - # if we cannot import OpenGL or whippersnappy, its an environment installation issue - if [[ "$opengl_error_message" =~ "ModuleNotFoundError" ]] || [[ "$opengl_error_message" =~ "ImportError" ]] + if [[ "$whippersnappy_check" == "NOT_INSTALLED" ]] then - echo "WARNING: The --qc_snap option of the corpus callosum module requires the Python packages PyOpenGL, glfw and" - echo " whippersnappy to be installed, but python could not import those three. Please install them and their" - echo " dependencies via 'pip install pyopengl glfw whippersnappy'." + echo "ERROR: The --qc_snap flag requires the 'whippersnappy' package (version >= 2.1) to generate the qc" + echo " thickness image, but whippersnappy is not installed in your Python environment." + elif [[ "$whippersnappy_check" == OLD_VERSION:* ]] + then + installed_version="${whippersnappy_check#OLD_VERSION:}" + echo "ERROR: The --qc_snap flag requires whippersnappy version >= 2.1 to generate the qc thickness image," + echo " but you only have version $installed_version installed." else - echo "WARNING: The --qc_snap option of the corpus callosum module requires OpenGL support, but we could not" - echo " create OpenGL handles. For Linux headless systems, you may install xvfb-run to provide a virtual display." + echo "ERROR: Failed to check whippersnappy installation: $whippersnappy_check" fi - echo " FastSurfer will not fail due to the unavailability of OpenGL, but some QC snapshots (rendered thickness" - echo " image) will not be created." + echo " Please install or upgrade whippersnappy with one of the following commands:" + echo " pip install 'whippersnappy>=2.1'" + exit 1 fi fi @@ -1181,10 +1188,9 @@ then # note: callosum manedit currently only affects inpainting and not internal FastSurferCC processing (surfaces etc) callosum_seg_manedit="$(add_file_suffix "$callosum_seg" "manedit")" # generate callosum segmentation, mesh, shape and downstream measure files - cmd=("${maybe_xvfb[@]}" $python "$CorpusCallosumDir/fastsurfer_cc.py" --sd "$sd" --sid "$subject" - "--threads" "$threads_seg" "--conformed_name" "$conformed_name" "--aseg_name" "$asegdkt_segfile" + cmd=($python "$CorpusCallosumDir/fastsurfer_cc.py" --sd "$sd" --sid "$subject" + "--threads" "$threads_seg" "--conformed_name" "$conformed_name" "--aseg_name" "$aseg_segfile" "--segmentation_in_orig" "$callosum_seg" "${cc_flags[@]}") - # if we are trying to create the thickness image in a headless setting, wrap call in xvfb-run echo_quoted "${cmd[@]}" | tee -a "$seg_log" "${wrap[@]}" "${cmd[@]}" 2>&1 | tee -a "$seg_log" exit_code=${PIPESTATUS[0]} diff --git a/tools/Docker/Dockerfile b/tools/Docker/Dockerfile index b0295c90..06f81f1b 100644 --- a/tools/Docker/Dockerfile +++ b/tools/Docker/Dockerfile @@ -272,11 +272,11 @@ RUN <