bugfixes for hires images, README updates, error handling

Author: ClePol

CorpusCallosum/README.md

@@ -136,7 +136,7 @@ The pipeline produces the following outputs in the specified output directory:
 - `surf/callosum.thickness.w`: FreeSurfer overlay file containing thickness values
 - `surf/callosum_mesh.vtk`: VTK format mesh file for 3D visualization 
 
-### Template Files (when --save_template is used)
+**Template Files (when --save_template is used):**
 
 - `contours.txt`: Corpus callosum contour coordinates for visualization
 - `thickness_values.txt`: Thickness measurements at each contour point

CorpusCallosum/cc_visualization.py

@@ -128,7 +128,7 @@ def main(
             contour_idx=len(cc_mesh.contours) // 2, save_path=str(output_dir / "midslice_2d.png")
         )
 
-        cc_mesh.to_fs_coordinates()
+        cc_mesh.to_fs_coordinates(vox_size=[resolution, resolution, resolution])
         cc_mesh.write_vtk(str(output_dir / "cc_mesh.vtk"))
         cc_mesh.write_fssurf(str(output_dir / "cc_mesh.fssurf"))
         cc_mesh.write_overlay(str(output_dir / "cc_mesh_overlay.curv"))

CorpusCallosum/fastsurfer_cc.py

@@ -532,8 +532,8 @@ def main(
     # Validate subdivision fractions
     for i in subdivisions:
         if i < 0 or i > 1:
-            print("Error: Subdivision fractions must be between 0 and 1, but got: ", i)
-            exit(1)
+            logger.error(f"Error: Subdivision fractions must be between 0 and 1, but got: {i}")
+            raise ValueError(f"Subdivision fractions must be between 0 and 1, but got: {i}")
 
     #### setup variables
     IO_processes = []
@@ -551,9 +551,9 @@ def main(
             "center around the mid-sagittal plane)"
         )
 
-    if not is_conform(orig):
+    if not is_conform(orig, vox_size='min', img_size=None):
         logger.error("Error: MRI is not conformed, please run conform.py or mri_convert to conform the image.")
-        exit(1)
+        raise ValueError("MRI is not conformed, please run conform.py or mri_convert to conform the image.")
 
     # load models
     device = torch.device("cuda" if torch.cuda.is_available() and not cpu else "cpu")
@@ -651,7 +651,8 @@ def main(
         cc_html_path=cc_html_path,
         vtk_file_path=vtk_file_path,
         thickness_image_path=thickness_image_path,
-        vox_size=orig.header.get_zooms()[0],
+        vox_size=orig.header.get_zooms(),
+        image_size=orig.shape,
         verbose=verbose,
         save_template=save_template,
     )

CorpusCallosum/segmentation/segmentation_postprocessing.py

@@ -100,7 +100,11 @@ def get_cc_volume_contour(desired_width_mm: int, cc_contours: list[np.ndarray],
 
     # Calculate cross-sectional areas for each contour
     areas = []
+    
     for contour in cc_contours:
+        contour = contour.copy()
+        assert voxel_size[1] == voxel_size[2], "volume must be isotropic"
+        contour *= voxel_size[1]
         # Calculate area using the shoelace formula for polygon area
         if contour.shape[1] < 3:
             areas.append(0.0)
@@ -117,15 +121,25 @@ def get_cc_volume_contour(desired_width_mm: int, cc_contours: list[np.ndarray],
 
     # Calculate spacing between slices (left-right direction)
     lr_spacing = voxel_size[0]  # x-direction voxel size
+
+    measurement_points = np.arange(-voxel_size[0]*(areas.shape[0]//2), 
+                                    voxel_size[0]*((areas.shape[0]+1)//2), lr_spacing)
 
-    # interpolate areas at 0 and 5
-    areas_interpolated = np.interp(x=[0, 5], xp=np.arange(lr_spacing/2, 5, lr_spacing), fp=areas)
+    # interpolate areas at 0.25 and 5
+    areas_interpolated = np.interp(x=[-2.5, 2.5], 
+                                   xp=measurement_points, 
+                                   fp=areas)
 
 
+    # remove measurement points that are outside of the desired range
+    # not sure if this can happen, but let's be safe
+    outside_range = (measurement_points < -2.5) | (measurement_points > 2.5)
+    measurement_points = [-2.5] + measurement_points[~outside_range].tolist() + [2.5]
+    areas = [areas_interpolated[0]] + areas[~outside_range].tolist() + [areas_interpolated[1]]
+    
 
-    measurements = [0,0.5,1.5,2.5,3.5,4.5,5]
-    # can also use cumulative trapezoidal rule
-    return integrate.simpson([areas_interpolated[0]] + areas.tolist() + [areas_interpolated[1]], x=measurements)
+    # can also use trapezoidal rule
+    return integrate.simpson(areas, x=measurement_points)
 
 
 def get_largest_cc(seg_arr: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
@@ -192,7 +206,7 @@ def clean_cc_segmentation(seg_arr: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
 
     if 250 not in unique_labels:
         clean_seg[seg_arr == 250] = 250
-        mask [seg_arr == 250] = True
+        mask[seg_arr == 250] = True
     if 192 not in unique_labels:
         clean_seg[seg_arr == 192] = 192
         mask[seg_arr == 192] = True

CorpusCallosum/shape/cc_mesh.py

@@ -1317,14 +1317,19 @@ def load_thickness_values(self, input_path: str, original_thickness_vertices_pat
                     try:
                         new_thickness_values[slice_idx][vertex_indices] = loaded_thickness_values[slice_idx][
                             ~np.isnan(loaded_thickness_values[slice_idx])]
-                    except IndexError:
-                        print(
+                    except IndexError as err:
+                        logger.error(
                             f"Tried to load "
                             f"{loaded_thickness_values[slice_idx][~np.isnan(loaded_thickness_values[slice_idx])]} "
                             f"values, but template has {new_thickness_values[slice_idx][vertex_indices]} values, "
                             "supply a correct template to visualize the thickness values"
                         )
-                        sys.exit(1)
+                        raise ValueError(
+                            f"Tried to load "
+                            f"{loaded_thickness_values[slice_idx][~np.isnan(loaded_thickness_values[slice_idx])]} "
+                            f"values, but template has {new_thickness_values[slice_idx][vertex_indices]} values, "
+                            "supply a correct template to visualize the thickness values"
+                        ) from err
             self.thickness_values = new_thickness_values
 
     @staticmethod
@@ -1334,15 +1339,17 @@ def __make_parent_folder(filename: str):
         output_folder = Path(filename).parent
         output_folder.mkdir(parents=False, exist_ok=True)
 
-    def to_fs_coordinates(self):
+    def to_fs_coordinates(self, vox_size: tuple[int, int, int], image_size: tuple[int, int, int]):
         """Convert mesh coordinates to FreeSurfer coordinate system.
 
         Transforms the mesh vertices from the original coordinate system to the
         FreeSurfer coordinate system by reordering axes and applying appropriate offsets.
         """
-        self.v = self.v[:, [2, 0, 1]]
-        self.v[:, 1] -= 128
-        self.v[:, 2] += 128
+        self.v = self.v[:, [2, 0, 1]] # LIA to ALI?
+        self.v *= (vox_size[0] **2) ## ???
+        self.v[:, 1] -= image_size[1] * vox_size[1] // 2 # move 0 to center of image
+        self.v[:, 2] += image_size[2] * vox_size[2] // 2
+        
 
     def write_fssurf(self, filename):
         """Write the mesh to a FreeSurfer surface file.

CorpusCallosum/shape/cc_postprocessing.py

@@ -134,7 +134,6 @@ def process_slice(segmentation, slice_idx, ac_coords, pc_coords, affine, num_thi
             'angular', or 'eigenvector')
         contour_smoothing (float): 
             Gaussian sigma for contour smoothing
-        
     Returns:
         slice_data (dict | None): 
             Dictionary containing measurements if successful, including:
@@ -246,7 +245,7 @@ def process_slice(segmentation, slice_idx, ac_coords, pc_coords, affine, num_thi
 def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac_coords, pc_coords, 
                   num_thickness_points, subdivisions, subdivision_method, contour_smoothing, 
                   debug_image_path=None, one_debug_image=False,
-                  thickness_image_path=None, vox_size=None, 
+                  thickness_image_path=None, vox_size=None, image_size=None,
                   save_template=None, surf_file_path=None, overlay_file_path=None, cc_html_path=None, 
                   vtk_file_path=None, verbose=False):
     """Process corpus callosum slices based on selection mode.
@@ -281,7 +280,7 @@ def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac
     if slice_selection == "middle":
         cc_mesh = CC_Mesh(num_slices=1)
         cc_mesh.set_acpc_coords(ac_coords, pc_coords)
-        cc_mesh.set_resolution(vox_size) # contour is always scaled to 1 mm
+        cc_mesh.set_resolution(vox_size[0]) # contour is always scaled to 1 mm
 
         # Process only the middle slice
         slice_idx = segmentation.shape[0] // 2
@@ -309,12 +308,12 @@ def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac
             if verbose:
                 logger.info(f"Saving segmentation qc image to {debug_image_path}")
             IO_processes.append(create_visualization(subdivision_method, result, midslices, 
-                                                   debug_image_path, ac_coords, pc_coords, vox_size))
+                                                   debug_image_path, ac_coords, pc_coords, vox_size[0]))
     else:
         num_slices = segmentation.shape[0]
         cc_mesh = CC_Mesh(num_slices=num_slices)
         cc_mesh.set_acpc_coords(ac_coords, pc_coords)
-        cc_mesh.set_resolution(vox_size) # contour is always scaled to 1 mm
+        cc_mesh.set_resolution(vox_size[0]) # contour is always scaled to 1 mm
 
         # Process multiple slices or specific slice
         if slice_selection == "all":
@@ -367,7 +366,7 @@ def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac
                     IO_processes.append(create_visualization(subdivision_method, result, 
                                                             midslices[current_slice_in_volume:current_slice_in_volume+1], 
                                                             debug_output_path_slice, ac_coords, pc_coords, 
-                                                            vox_size, f' (Slice {slice_idx})'))
+                                                            vox_size[0], f' (Slice {slice_idx})'))
 
     if save_template is not None:
         # Convert to Path object and ensure directory exists
@@ -394,7 +393,7 @@ def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac
         #cc_mesh.write_vtk(str(output_dir / 'cc_mesh.vtk'))
 
 
-        cc_mesh.to_fs_coordinates()
+        cc_mesh.to_fs_coordinates(vox_size=vox_size, image_size=image_size)
 
         if overlay_file_path is not None:
             if verbose: 
@@ -417,7 +416,7 @@ def process_slices(segmentation, slice_selection, temp_seg_affine, midslices, ac
 
     if not slice_results:
         logger.error("Error: No valid slices were found for postprocessing")
-        exit(1)
+        raise ValueError("No valid slices were found for postprocessing")
 
     return slice_results, IO_processes