diff --git a/CerebNet/data_loader/dataset.py b/CerebNet/data_loader/dataset.py index 2781232e6..a9e3fcf45 100644 --- a/CerebNet/data_loader/dataset.py +++ b/CerebNet/data_loader/dataset.py @@ -242,6 +242,8 @@ def __init__( slice_thickness: int, primary_slice: str | None = None, ): + from numpy.linalg import inv + self.slice_thickness = slice_thickness self.transforms = Compose([ToTensor()]) self.img_org = img_org @@ -252,8 +254,6 @@ def __init__( # binarize the cerebellum from brain_seg cereb_aseg_mask = utils.get_aseg_cereb_mask(np.asarray(brain_seg.dataobj)) - from numpy.linalg import inv - affine = inv(brain_seg.affine) @ img_org.affine # print(brain_seg.affine, img_org.affine) @@ -264,9 +264,7 @@ def __init__( ) from scipy.ndimage import affine_transform - cereb_aseg = affine_transform( - cereb_aseg_mask.astype(np.float32), affine, output_shape=img_org.shape - ) + cereb_aseg = affine_transform(cereb_aseg_mask.astype(np.float32), affine, output_shape=img_org.shape) cereb_aseg_mask = cereb_aseg > 0.5 bbox = self.locate_mask_bbox(cereb_aseg_mask) @@ -274,17 +272,11 @@ def __init__( # create the roi from cereb_aseg (where labels after interpolation > 0.05 --> membership rounded to 1 decimal) self.roi: LocalizerROI = { "source_shape": img_org.shape, - "offsets": bounding_volume_offset( - bbox, patch_size, image_shape=cereb_aseg_mask.shape - ), + "offsets": bounding_volume_offset(bbox, patch_size, image_shape=cereb_aseg_mask.shape), "target_shape": patch_size, } # crop the region of interest - img = crop_transform( - self.img_org_data, - offsets=self.roi["offsets"], - target_shape=self.roi["target_shape"], - ) + img = crop_transform(self.img_org_data, offsets=self.roi["offsets"], target_shape=self.roi["target_shape"]) # reorient the data to lia img_lia, self.back_to_native = to_target_orientation(img, self.img_org.affine, target_orientation="LIA") @@ -298,13 +290,9 @@ def __init__( } for plane, data_i in data.items(): # data is transformed to 'plane'-direction in axis 2 - thick_slices = get_thick_slices( - data_i, self.slice_thickness - ) # [H, W, n_slices, C] + thick_slices = get_thick_slices(data_i, self.slice_thickness) # [H, W, n_slices, C] # it seems x and y are flipped with respect to expectations here - self.images_per_plane[plane] = np.transpose( - thick_slices, (2, 0, 1, 3) - ) # [n_slices, H, W, C] + self.images_per_plane[plane] = np.transpose(thick_slices, (2, 0, 1, 3)) # [n_slices, H, W, C] def locate_mask_bbox(self, mask: npt.NDArray[bool]): """Find the largest connected component of the mask. @@ -329,9 +317,7 @@ def get_bounding_offsets(self) -> LocalizerROI: def set_plane(self, plane: Plane): """Set the active plane.""" if plane not in self.images_per_plane.keys(): - raise ValueError( - f"Invalid plane name, must be in {tuple(self.images_per_plane.keys())}" - ) + raise ValueError(f"Invalid plane name, must be in {tuple(self.images_per_plane.keys())}") self._plane = plane @property diff --git a/CerebNet/inference.py b/CerebNet/inference.py index fdff34949..25449a0a5 100644 --- a/CerebNet/inference.py +++ b/CerebNet/inference.py @@ -56,6 +56,7 @@ class Inference: cerebnet_labels: Mapper[str, int] cereb_name2fs_id: Mapper[str, int] freesurfer_name2id: Mapper[str, int] + cereb_id2fs_id: Mapper[int, int] def __init__( self, @@ -129,9 +130,7 @@ def __init__( self.viewagg_device = _viewagg_device - def prep_lut( - file: Path, *args, **kwargs, - ) -> Future[TSVLookupTable | JsonColorLookupTable]: + def prep_lut(file: Path, *args, **kwargs) -> Future[TSVLookupTable | JsonColorLookupTable]: _cls = TSVLookupTable cls = {".json": JsonColorLookupTable, ".txt": _cls, ".tsv": _cls} return self.pool.submit(cls[file.suffix], file, *args, **kwargs) @@ -153,12 +152,8 @@ def lut_path(module: str, file: str) -> Path: self.cerebnet_labels = _cerebnet_mapper.result().labelname2id() self.freesurfer_name2id = fs_color_map.result().labelname2id() - cereb_name2fs_name: Mapper[str, str] = ( - cereb2freesurfer_mapper.result().labelname2id() - ) - cerebsag_name2cereb_name: Mapper[str, str] = ( - sagittal_cereb2cereb_mapper.result().labelname2id() - ) + cereb_name2fs_name: Mapper[str, str] = cereb2freesurfer_mapper.result().labelname2id() + cerebsag_name2cereb_name: Mapper[str, str] = sagittal_cereb2cereb_mapper.result().labelname2id() cereb_id2name = self.cerebnet_labels.__reversed__() self.cereb_name2fs_id = cereb_name2fs_name.chain(self.freesurfer_name2id) @@ -204,13 +199,9 @@ def __load_model(cfg: "yacs.config.CfgNode", plane: Plane) -> torch.nn.Module: return dict(zip(PLANES, self.pool.map(_load_model_func, PLANES), strict=False)) @torch.no_grad() - def _predict_single_subject( - self, subject_dataset: SubjectDataset - ) -> dict[Plane, list[torch.Tensor]]: - """Predict the classes based on a SubjectDataset.""" - img_loader = DataLoader( - subject_dataset, batch_size=self.batch_size, shuffle=False - ) + def _predict_single_subject(self, subject_dataset: SubjectDataset) -> dict[Plane, list[torch.Tensor]]: + """Predict the classes based on a SubjectDataset (operates fully in LIA).""" + img_loader = DataLoader(subject_dataset, batch_size=self.batch_size, shuffle=False) prediction_logits = {} try: for plane in PLANES: @@ -220,8 +211,7 @@ def _predict_single_subject( from CerebNet.data_loader.data_utils import slice_lia2ras, slice_ras2lia for img in img_loader: - # CerebNet is trained on RAS+ conventions, so we need to map between - # lia (FastSurfer) and RAS+ + # CerebNet is trained on RAS+ conventions, so we need to map between lia (FastSurfer) and RAS+ # map LIA 2 RAS img = slice_lia2ras(plane, img, thick_slices=True) batch = img.to(self.device) @@ -367,9 +357,7 @@ def _save_cerebnet_seg( if cerebnet_seg.shape != orig.shape: raise RuntimeError("Cereb segmentation shape inconsistent with Orig shape!") logger.info(f"Saving CerebNet cerebellum segmentation at {filename}") - return self.pool.submit( - save_image, orig.header, orig.affine, cerebnet_seg, filename, dtype=np.int16 - ) + return self.pool.submit(save_image, orig.header, orig.affine, cerebnet_seg, filename, dtype=np.int16) def _get_subject_dataset( self, subject: SubjectDirectory @@ -381,7 +369,7 @@ def _get_subject_dataset( from FastSurferCNN.data_loader.data_utils import load_image, load_maybe_conform - norm_file, norm_data, norm = None, None, None + norm_file, norm_data, _norm = None, None, None if subject.has_attribute("cereb_statsfile") : if not subject.can_resolve_attribute("cereb_statsfile"): from FastSurferCNN.utils.parser_defaults import ALL_FLAGS @@ -401,7 +389,7 @@ def _get_subject_dataset( norm_file = subject.filename_by_attribute("norm_name") # finally, load the bias field file - norm = self.pool.submit(load_maybe_conform, norm_file, norm_file, **self._conform_kwargs) + _norm = self.pool.submit(load_maybe_conform, norm_file, norm_file, **self._conform_kwargs) # localization if not subject.fileexists_by_attribute("asegdkt_segfile"): @@ -409,26 +397,23 @@ def _get_subject_dataset( f"The aseg.DKT-segmentation file '{subject.asegdkt_segfile}' did not " f"exist, please run FastSurferVINN first." ) - seg = self.pool.submit( - load_image, subject.filename_by_attribute("asegdkt_segfile") - ) + _seg = self.pool.submit(load_image, subject.filename_by_attribute("asegdkt_segfile")) # create conformed image - conf_img = self.pool.submit( + _conf_img = self.pool.submit( load_maybe_conform, subject.filename_by_attribute("conf_name"), subject.filename_by_attribute("orig_name"), **self._conform_kwargs, ) - seg, seg_data = seg.result() - conf_file, conf_img, conf_data = conf_img.result() + seg, seg_data = _seg.result() + conf_file, conf_img, conf_data = _conf_img.result() if np.allclose(conf_img.header.get_zooms(), 1.0, atol=0.01): logger.warning( "CerebNet does not support images that are not conformed to 1.0mm. We detected a voxel sizes of " f"{tuple(conf_img.header.get_zooms())} in {conf_file}!" ) - subject_dataset = SubjectDataset( img_org=conf_img, brain_seg=seg, @@ -437,8 +422,8 @@ def _get_subject_dataset( # obsolete: primary_slice=self.cfg.DATA.PRIMARY_SLICE_DIR, ) subject_dataset.transforms = ToTensorTest() - if norm is not None: - norm_file, _, norm_data = norm.result() + if _norm is not None: + norm_file, _, norm_data = _norm.result() return norm_data, norm_file, subject_dataset def run(self, subject_dirs: SubjectList): @@ -454,16 +439,13 @@ def run(self, subject_dirs: SubjectList): from FastSurferCNN.utils.common import iterate iter_subjects = iterate(self.pool, self._get_subject_dataset, subject_dirs) futures = [] - for idx, (subject, (norm, norm_file, subject_dataset)) in tqdm( - enumerate(iter_subjects), total=len(subject_dirs), desc="Subject", - ): + for idx, (subject, _data) in tqdm(enumerate(iter_subjects), total=len(subject_dirs), desc="Subject"): + norm, norm_file, subject_dataset = _data try: # predict CerebNet, returns logits (input and output are LIA) preds = self._predict_single_subject(subject_dataset) # create the folder for the output file, if it does not exist - _mkdir = self.pool.submit( - subject.segfile.parent.mkdir, exist_ok=True, parents=True, - ) + _mkdir = self.pool.submit(subject.segfile.parent.mkdir, exist_ok=True, parents=True) # postprocess logits (move axes, map sagittal to all classes, still LIA) preds_per_plane = self._post_process_preds(preds) @@ -487,11 +469,7 @@ def run(self, subject_dirs: SubjectList): # this is None, but synchronizes the creation of the directory _ = _mkdir.result() futures.append( - self._save_cerebnet_seg( - full_cereb_seg, - subject.segfile, - subject_dataset.get_nibabel_img(), - ) + self._save_cerebnet_seg(full_cereb_seg, subject.segfile, subject_dataset.get_nibabel_img()) ) if subject.has_attribute("cereb_statsfile"): @@ -521,10 +499,9 @@ def run(self, subject_dirs: SubjectList): ) ) - logger.info( - f"Subject {idx + 1}/{len(subject_dirs)} with id '{subject.id}' processed in " - f"{pred_time - start_time :.2f} sec." - ) + duration = pred_time - start_time + num = len(subject_dirs) + logger.info(f"Subject {idx + 1}/{num} with id '{subject.id}' processed in {duration:.2f} sec.") except Exception as e: logger.exception(e) return "\n".join(map(str, e.args)) diff --git a/CerebNet/run_prediction.py b/CerebNet/run_prediction.py index 7df291991..d8d432f8b 100644 --- a/CerebNet/run_prediction.py +++ b/CerebNet/run_prediction.py @@ -45,16 +45,11 @@ def setup_options(): # Training settings parser = argparse.ArgumentParser(description="Segmentation") - # 1. Directory information (where to read from, where to write from and to incl. - # search-tag) - parser = parser_defaults.add_arguments( - parser, ["in_dir", "tag", "csv_file", "sd", "sid", "remove_suffix"] - ) + # 1. Directory information (where to read from, where to write from and to incl. search-tag) + parser = parser_defaults.add_arguments(parser, ["in_dir", "tag", "csv_file", "sd", "sid", "remove_suffix"]) # 2. Options for the MRI volumes - parser = parser_defaults.add_arguments( - parser, ["t1", "conformed_name", "norm_name", "asegdkt_segfile"] - ) + parser = parser_defaults.add_arguments(parser, ["t1", "conformed_name", "norm_name", "asegdkt_segfile"]) parser.add_argument( "--cereb_segfile", dest="cereb_segfile", @@ -136,10 +131,9 @@ def main(args: argparse.Namespace) -> int | str: Returns ------- - int - Returns 0 upon successful execution to indicate success. - str - A message indicating the failure reason in case of an exception. + int, str + Returns 0 upon successful execution to indicate success or + a message indicating the failure reason in case of an exception. References ---------- diff --git a/Docker/README.md b/Docker/README.md index 7b24729c6..5a68f002f 100644 --- a/Docker/README.md +++ b/Docker/README.md @@ -190,13 +190,6 @@ To build a docker image with attestation and provenance, i.e. Software Bill Of M [[worker.containerd.gcpolicy]] all = true keepBytes = 1024000000 - # settings to push to a "local", registry with self-signed certificates - # see for example https://tech.paulcz.net/2016/01/secure-docker-with-tls/ https://github.com/paulczar/omgwtfssl - [registry."host:5000"] - ca=["/path/to/registry/ssl/ca.pem"] - [[registry."landau.dzne.ds:5000".keypair]] - key="/path/to/registry/ssl/key.pem" - cert="/path/to/registry/ssl/cert.pem" ``` 3. Attestation files are not supported by the standard docker image storage driver. Therefore, images cannot be tested locally. There are two solutions to this limitation. @@ -231,7 +224,7 @@ rocms=("rocm$rocm") # end of config # code -git clone --branch stable --single-branch gtihub.com/Deep-MI/FastSurfer $build_dir +git clone --branch stable --single-branch github.com/Deep-MI/FastSurfer $build_dir cd $build_dir all_tags=("latest" "gpu-latest" "cuda-v$version" "rocm-v$version" "cpu-latest") # build all distinct images diff --git a/FastSurferCNN/data_loader/conform.py b/FastSurferCNN/data_loader/conform.py index b057ca437..e403dfc28 100644 --- a/FastSurferCNN/data_loader/conform.py +++ b/FastSurferCNN/data_loader/conform.py @@ -1412,16 +1412,14 @@ def crop_transform( _target_shape = image.shape[:-len_off] + tuple( i - 2 * o for i, o in zip(image.shape[-len_off:], offsets, strict=False) ) - elif len_off != len(target_shape): - raise ValueError( - "Incompatible offset and target_shape dimensionality (at least once)." - ) - else: + elif len_off == len(target_shape): _target_shape = tuple( i if t == -1 else t for i, t in zip(image.shape[-len_off:], target_shape, strict=False) ) _target_shape = image.shape[:-len_off] + _target_shape + else: + raise ValueError("Incompatible offset and target_shape dimensionality (at least once).") if len_off > image.ndim: raise RuntimeError("shape of offsets is larger than dim of image allows.") diff --git a/FastSurferCNN/data_loader/data_utils.py b/FastSurferCNN/data_loader/data_utils.py index 67dc9524a..7c6a8b1db 100644 --- a/FastSurferCNN/data_loader/data_utils.py +++ b/FastSurferCNN/data_loader/data_utils.py @@ -206,31 +206,22 @@ def load_maybe_conform( dst_file = file else: # the image is not conformed to 1mm, do this now. - - fileext = [ - ext for ext in SUPPORTED_OUTPUT_FILE_FORMATS - if file.name.endswith("." + ext) - ] + fileext = [ext for ext in SUPPORTED_OUTPUT_FILE_FORMATS if file.name.endswith("." + ext)] if len(fileext) != 1: raise RuntimeError( - f"Invalid file extension of conf_name: {file}, must be one of " - f"{SUPPORTED_OUTPUT_FILE_FORMATS}." + f"Invalid file extension of conf_name: {file}, must be one of {SUPPORTED_OUTPUT_FILE_FORMATS}." ) file_no_fileext = str(file)[:-len(fileext[0]) - 1] - if (vox_size := conform_kwargs.get("vox_size", 1.0)) == "min": - vox_suffix = ".min" - else: - vox_suffix = f".{str(vox_size).replace('.', '')}mm" + vox_size = conform_kwargs.get("vox_size", 1.0) + vox_suffix = ".min" if vox_size == "min" else f".{str(vox_size).replace('.', '')}mm" if not file_no_fileext.endswith(vox_suffix): file_no_fileext += vox_suffix - # if the orig file is neither absolute nor in the subject path, use the - # conformed file + # if the orig file is neither absolute nor in the subject path, use the conformed file src_file = alt_file if alt_file.is_file() else file if not alt_file.is_file(): LOGGER.warning( - f"No valid alternative file (e.g. orig, here: {alt_file}) was given to " - f"interpolate from, so we might lose quality due to multiple chained " - f"interpolations." + f"No valid alternative file (e.g. orig, here: {alt_file}) was given to interpolate from, so we might " + f"lose quality due to multiple chained interpolations. " ) dst_file = Path(file_no_fileext + "." + fileext[0]) @@ -272,13 +263,8 @@ def save_image( Image array type; if provided, the image object is explicitly set to match this type. """ save_as = Path(save_as) - assert ( - save_as.suffix[1:] in SUPPORTED_OUTPUT_FILE_FORMATS or - save_as.suffixes[-2:] == [".nii", ".gz"] - ), ( - f"Output filename does not contain a supported file format " - f"{SUPPORTED_OUTPUT_FILE_FORMATS}!" - ) + valid_ext = save_as.suffix[1:] in SUPPORTED_OUTPUT_FILE_FORMATS or save_as.suffixes[-2:] == [".nii", ".gz"] + assert valid_ext, f"Output filename does not contain a supported file format {SUPPORTED_OUTPUT_FILE_FORMATS}!" mgh_img = None if save_as.suffix == ".mgz": diff --git a/FastSurferCNN/data_loader/dataset.py b/FastSurferCNN/data_loader/dataset.py index ce7a54a76..e40487b54 100644 --- a/FastSurferCNN/data_loader/dataset.py +++ b/FastSurferCNN/data_loader/dataset.py @@ -14,6 +14,7 @@ # IMPORTS import time +from collections.abc import Callable from typing import Optional import h5py @@ -40,7 +41,7 @@ def __init__( orig_data: npt.NDArray, orig_zoom: npt.NDArray, cfg: yacs.config.CfgNode, - transforms: Optional = None + transforms: Callable[[npt.NDArray[float]], npt.NDArray[float]] | None = None, ): """ Construct object. @@ -50,15 +51,14 @@ def __init__( orig_data : npt.NDArray Original Data. orig_zoom : npt.NDArray - Original zoomfactors. + Original zoom factors. cfg : yacs.config.CfgNode Configuration Node. - transforms : Optional - Transformer for the image. Defaults to None. + transforms : callable[[npt.NDArray[float]], npt.NDArray[float]], optional + Transforms for the image, defaults to no transformation. """ - assert ( - orig_data.max() > 0.8 - ), f"Multi Dataset - orig fail, max removed {orig_data.max()}" + orig_max = orig_data.max() + assert orig_max > 0.8, f"Multi Dataset - orig fail, max removed {orig_max}" self.plane = cfg.DATA.PLANE self.slice_thickness = cfg.MODEL.NUM_CHANNELS // 2 self.base_res = cfg.MODEL.BASE_RES diff --git a/FastSurferCNN/inference.py b/FastSurferCNN/inference.py index 9ba0306ad..fa0a812b7 100644 --- a/FastSurferCNN/inference.py +++ b/FastSurferCNN/inference.py @@ -98,7 +98,7 @@ def __init__( ckpt : str String or os.PathLike object containing the name to the checkpoint file (Default value = ""). lut : str, np.ndarray, DataFrame, optional - Lookup table for mapping (Default value = None). + Lookup table for mapping. """ # Set random seed from configs. np.random.seed(cfg.RNG_SEED) @@ -154,9 +154,7 @@ def setup_model(self, cfg=None, device: torch.device = None): device = self.default_device # Set up model - self._model_not_init = build_model( - self.cfg - ) # ~ model = FastSurferCNN(params_network) + self._model_not_init = build_model(self.cfg) # ~ model = FastSurferCNN(params_network) self._model_not_init.to(device) self.device = None @@ -181,9 +179,7 @@ def to(self, device: torch.device | None = None): The desired device of the parameters and buffers in this module (Default value = None). """ if self.model_parallel: - raise RuntimeError( - "Moving the model to other devices is not supported for multi-device models." - ) + raise RuntimeError("Moving the model to other devices is not supported for multi-device models.") _device = self.default_device if device is None else device self.device = _device self.model.to(device=_device) @@ -203,19 +199,17 @@ def load_checkpoint(self, ckpt: str | os.PathLike): # If device is None, the model has never been loaded (still in random initial configuration) if self.device is None: self.device = self.default_device + load_device = self.device # workaround for mps (directly loading to map_location=mps results in zeros) - device = self.device if self.device.type == "mps": - self.model.to("cpu") - device = "cpu" - else: - # make sure the model is, where it is supposed to be - self.model.to(self.device) + load_device = "cpu" + # make sure the model is, where it is supposed to be + self.model.to(load_device) # WARNING: weights_only=False can cause unsafe code execution, but here the # checkpoint can be considered to be from a safe source - model_state = torch.load(ckpt, map_location=device, weights_only=False) + model_state = torch.load(ckpt, map_location=load_device, weights_only=False) self.model.load_state_dict(model_state["model_state"]) # workaround for mps (move the model back to mps) @@ -335,7 +329,7 @@ def eval( Validation loader. out_scale : Optional Output scale (Default value = None). - out : Optional[torch.Tensor] + out : torch.Tensor, optional Previous prediction tensor (Default value = None). Returns @@ -366,14 +360,10 @@ def eval( log_batch_idx = None with logging_redirect_tqdm(): try: - for batch_idx, batch in tqdm( - enumerate(val_loader), total=len(val_loader), unit="batch" - ): + for batch_idx, batch in tqdm(enumerate(val_loader), total=len(val_loader), unit="batch"): log_batch_idx = batch_idx # move data to the model device - images, scale_factors = batch["image"].to(self.device), batch[ - "scale_factor" - ].to(self.device) + images, scale_factors = batch["image"].to(self.device), batch["scale_factor"].to(self.device) # predict the current batch, outputs logits pred = self.model(images, scale_factors, out_scale) @@ -382,14 +372,10 @@ def eval( # check if we need a special mapping (e.g. as for sagittal) if self.get_plane() == "sagittal": - pred = map_prediction_sagittal2full( - pred, num_classes=self.get_num_classes(), lut=self.lut - ) + pred = map_prediction_sagittal2full(pred, num_classes=self.get_num_classes(), lut=self.lut) # permute the prediction into the out slice order - pred = pred.permute(*self.permute_order[plane]).to( - out.device - ) # the to-operation is implicit + pred = pred.permute(*self.permute_order[plane]).to(out.device) # the to-operation is implicit # cut prediction to the image size pred = pred[pred_ii] @@ -400,14 +386,10 @@ def eval( start_index = end_index except: - logger.exception( - f"Exception in batch {log_batch_idx} of {plane} inference." - ) + logger.exception(f"Exception in batch {log_batch_idx + 1} of {plane} inference.") raise else: - logger.info( - f"Inference on {batch_idx + 1} batches for {plane} successful" - ) + logger.info(f"Inference on {log_batch_idx + 1} batches for {plane} successful") return out @@ -420,7 +402,7 @@ def run( orig_zoom: npt.NDArray, out: torch.Tensor | None = None, out_res: int | None = None, - batch_size: int = None, + batch_size: int | None = None, ) -> torch.Tensor: """ Run the loaded model on the data (T1) from orig_data and @@ -440,12 +422,12 @@ def run( Updated output tensor (Default = None). out_res : Optional[int] Output resolution (Default value = None). - batch_size : int - Batch size (Default = None). + batch_size : int, optional + Batch size. Returns ------- - toch.Tensor + torch.Tensor Prediction probability tensor. """ # Set up DataLoader @@ -467,8 +449,7 @@ def run( out = self.eval(init_pred, test_data_loader, out=out, out_scale=out_res) time_delta = time.time() - start logger.info( - f"{self.cfg.DATA.PLANE.capitalize()} inference on {img_filename} finished in " - f"{time_delta:0.4f} seconds" + f"{self.cfg.DATA.PLANE.capitalize()} inference on {img_filename} finished in {time_delta:0.4f} seconds" ) return out diff --git a/FastSurferCNN/run_prediction.py b/FastSurferCNN/run_prediction.py index ce74ddfdb..48bac2331 100644 --- a/FastSurferCNN/run_prediction.py +++ b/FastSurferCNN/run_prediction.py @@ -13,8 +13,7 @@ # limitations under the License. """ -This is the FastSurfer/run_prediction.py script, the backbone for whole brain -segmentation. +This is the FastSurfer/run_prediction.py script, the backbone for whole brain segmentation. Usage: @@ -76,9 +75,7 @@ def set_up_cfgs( batch_size: int = 1, ) -> yacs.config.CfgNode: """ - Set up configuration. - - Sets up configurations with given arguments inside the yaml file. + Set up configuration with given arguments inside the yaml file. Parameters ---------- @@ -239,13 +236,8 @@ def __init__( if self.device.type == "cpu" and viewagg_device in ("auto", "cpu"): self.viewagg_device = self.device else: - # check, if GPU is big enough to run view agg on it - # (this currently takes the memory of the passed device) - self.viewagg_device = find_device( - viewagg_device, - flag_name="viewagg_device", - min_memory=4 * (2**30), - ) + # check, if GPU is big enough to run view agg on it (this currently takes the memory of the passed device) + self.viewagg_device = find_device(viewagg_device, flag_name="viewagg_device", min_memory=4 * (2**30)) LOGGER.info(f"Running view aggregation on {self.viewagg_device}") @@ -253,38 +245,28 @@ def __init__( self.lut = du.read_classes_from_lut(lut) except FileNotFoundError as err: raise ValueError( - f"Could not find the ColorLUT in {lut}, please make sure the " - f"--lut argument is valid." + f"Could not find the ColorLUT in {lut}, please make sure the --lut argument is valid." ) from err self.labels = self.lut["ID"].values self.torch_labels = torch.from_numpy(self.lut["ID"].values) self.names = ["SubjectName", "Average", "Subcortical", "Cortical"] - self.cfg_fin, cfg_cor, cfg_sag, cfg_ax = args2cfg( - cfg_ax, cfg_cor, cfg_sag, batch_size=batch_size, - ) + self.cfg_fin, cfg_cor, cfg_sag, cfg_ax = args2cfg(cfg_ax, cfg_cor, cfg_sag, batch_size=batch_size) # the order in this dictionary dictates the order in the view aggregation self.view_ops = { "coronal": {"cfg": cfg_cor, "ckpt": ckpt_cor}, "sagittal": {"cfg": cfg_sag, "ckpt": ckpt_sag}, "axial": {"cfg": cfg_ax, "ckpt": ckpt_ax}, } - self.num_classes = max( - view["cfg"].MODEL.NUM_CLASSES for view in self.view_ops.values() - ) + self.num_classes = max(view["cfg"].MODEL.NUM_CLASSES for view in self.view_ops.values()) self.models = {} for plane, view in self.view_ops.items(): if all(view[key] is not None for key in ("cfg", "ckpt")): - self.models[plane] = Inference( - view["cfg"], ckpt=view["ckpt"], device=self.device, lut=self.lut, - ) + self.models[plane] = Inference(view["cfg"], ckpt=view["ckpt"], device=self.device, lut=self.lut) try: self.vox_size = _vox_size(vox_size) except argparse.ArgumentParser: - raise ValueError( - f"Invalid value for vox_size, must be between 0 and 1 or 'min', was " - f"{vox_size}." - ) from None + raise ValueError(f"Invalid value for vox_size, must be between 0 and 1 or 'min', was {vox_size}.") from None self.conform_to_1mm_threshold = conform_to_1mm_threshold @property @@ -294,10 +276,7 @@ def pool(self) -> Executor: specified in __init__). """ if not hasattr(self, "_pool"): - if not self._async_io: - self._pool = SerialExecutor() - else: - self._pool = ThreadPoolExecutor(self._threads) + self._pool = ThreadPoolExecutor(self._threads) if self._async_io else SerialExecutor() return self._pool def __del__(self): @@ -342,7 +321,7 @@ def conform_and_save_orig( if (self.orientation is None or self.orientation == "native") and \ not is_conform(orig, **self.__conform_kwargs(verbose=False, dtype=None)): raise RuntimeError( - f"To store images in native image space, the input image {subject.orig_name} must have isometric " + f"To store images in native image space, the input image {subject.orig_name} must have isotropic " f"voxels." ) LOGGER.info("Conforming image...") @@ -354,10 +333,7 @@ def conform_and_save_orig( self.async_save_img(subject.conf_name, orig_data, orig, dtype=np.uint8) LOGGER.info(f"Saving conformed image to {subject.conf_name}...") else: - raise RuntimeError( - "Cannot resolve the name to the conformed image, please specify an " - "absolute path." - ) + raise RuntimeError("Cannot resolve the name to the conformed image, please specify an absolute path.") return orig, orig_data @@ -426,8 +402,7 @@ def get_prediction( # map to freesurfer label space pred_classes = du.map_label2aparc_aseg(pred_classes, self.labels) # return numpy array - # TODO: split_cortex_labels requires a numpy ndarray input, maybe we can also - # use Mapper here + # TODO: split_cortex_labels requires a numpy ndarray input, maybe we can also use Mapper here pred_classes = du.split_cortex_labels(pred_classes.cpu().numpy()) return pred_classes @@ -450,16 +425,12 @@ def save_img( orig : nib.analyze.SpatialImage Original Image. dtype : type, optional - Data type to use for saving the image. If None, the original data type is - used (Default value = None). + Data type to use for saving the image. If None, the original data type is used. """ save_as = Path(save_as) # Create output directory if it does not already exist. if not save_as.parent.exists(): - LOGGER.info( - f"Output image directory {save_as.parent} does not exist. " - f"Creating it now..." - ) + LOGGER.info(f"Output image directory {save_as.parent} does not exist. Creating it now...") save_as.parent.mkdir(parents=True) np_data = data if isinstance(data, np.ndarray) else data.cpu().numpy() @@ -469,9 +440,7 @@ def save_img( else: _header = orig.header du.save_image(_header, orig.affine, np_data, save_as, dtype=dtype) - LOGGER.info( - f"Successfully saved image {'asynchronously ' if self._async_io else ''} as {save_as}." - ) + LOGGER.info(f"Successfully saved image {'asynchronously ' if self._async_io else ''}as {save_as}.") def async_save_img( self, @@ -488,19 +457,17 @@ def async_save_img( ---------- save_as : str, Path Filename to give the image. - data : Union[np.ndarray, torch.Tensor] + data : np.ndarray, torch.Tensor Image data. orig : nib.analyze.SpatialImage Original Image. dtype : type, optional - Data type to use for saving the image. If None, the original data type is - used. + Data type to use for saving the image. If None, the original data type is used. Returns ------- Future[None] - A Future object to synchronize (and catch/handle exceptions in the save_img - method). + A Future object to synchronize (and catch/handle exceptions in the save_img method). """ return self.pool.submit(self.save_img, save_as, data, orig, dtype) @@ -566,39 +533,22 @@ def make_parser(): # 1. Options for input directories and filenames parser = parser_defaults.add_arguments( - parser, ["t1", "sid", "in_dir", "tag", "csv_file", "lut", "remove_suffix"] + parser, + ["t1", "sid", "in_dir", "tag", "csv_file", "lut", "remove_suffix"], ) # 2. Options for output parser = parser_defaults.add_arguments( parser, - [ - "asegdkt_segfile", - "conformed_name", - "brainmask_name", - "aseg_name", - "sd", - "seg_log", - "qc_log", - ], + ["asegdkt_segfile", "conformed_name", "brainmask_name", "aseg_name", "sd", "seg_log", "qc_log"], ) # 3. Checkpoint to load files: dict[Plane, str | Path] = {k: "default" for k in PLANES} - parser = parser_defaults.add_plane_flags( - parser, - "checkpoint", - files, - CHECKPOINT_PATHS_FILE - ) + parser = parser_defaults.add_plane_flags(parser, "checkpoint", files, CHECKPOINT_PATHS_FILE) # 4. CFG-file with default options for network - parser = parser_defaults.add_plane_flags( - parser, - "config", - files, - CHECKPOINT_PATHS_FILE - ) + parser = parser_defaults.add_plane_flags(parser, "config", files, CHECKPOINT_PATHS_FILE) # 5. technical parameters image_flags = ["vox_size", "conform_to_1mm_threshold", "orientation", "image_size", "device"] @@ -715,19 +665,13 @@ def main( try: # The orig_t1_file is only used to populate verbose messages here pred_data = eval.get_prediction(subject.orig_name, data_array, orig_img.header.get_zooms(), orig_img.affine) - futures.append( - eval.async_save_img( - subject.segfile, pred_data, orig_img, dtype=np.int16 - ) - ) + futures.append(eval.async_save_img(subject.segfile, pred_data, orig_img, dtype=np.int16)) # Create aseg and brainmask - # There is a funny edge case in legacy FastSurfer 2.0, where the behavior is - # not well-defined, if orig_name is an absolute path, but out_dir is not - # set. Then, we would create a sub-folder in the folder of orig_name using - # the subject_id (passed by --sid or extracted from the orig_name) and use - # that as the subject folder. + # There is a funny edge case in legacy FastSurfer 2.0, where the behavior is not well-defined, if orig_name + # is an absolute path, but out_dir is not set. Then, we would create a sub-folder in the folder of orig_name + # using the subject_id (passed by --sid or extracted from the orig_name) and use that as the subject folder. bm = None store_brainmask = subject.can_resolve_filename(brainmask_name) store_aseg = subject.can_resolve_filename(aseg_name) @@ -737,17 +681,13 @@ def main( if store_brainmask: # get mask mask_name = subject.filename_in_subject_folder(brainmask_name) - futures.append( - eval.async_save_img(mask_name, bm, orig_img, dtype=np.uint8) - ) + futures.append(eval.async_save_img(mask_name, bm, orig_img, dtype=np.uint8)) else: - LOGGER.info( - "Not saving the brainmask, because we could not figure out where " - "to store it. Please specify a subject id with {sid[flag]}, or an " - "absolute brainmask path with {brainmask_name[flag]}.".format( - **subjects.flags, - ) + message = ( + "Not saving the brainmask, because we could not figure out where to store it. Please specify a " + "subject id with {sid[flag]}, or an absolute brainmask path with {brainmask_name[flag]}." ) + LOGGER.info(message.format(**subjects.flags)) if store_aseg: # reduce aparc to aseg and mask regions @@ -757,26 +697,19 @@ def main( aseg = rta.flip_wm_islands(aseg) aseg_name = subject.filename_in_subject_folder(aseg_name) # Change datatype to np.uint8, else mri_cc will fail! - futures.append( - eval.async_save_img(aseg_name, aseg, orig_img, dtype=np.uint8) - ) + futures.append(eval.async_save_img(aseg_name, aseg, orig_img, dtype=np.uint8)) else: - LOGGER.info( - "Not saving the aseg file, because we could not figure out where " - "to store it. Please specify a subject id with {sid[flag]}, or an " - "absolute aseg path with {aseg_name[flag]}.".format( - **subjects.flags, - ) + message = ( + "Not saving the aseg file, because we could not figure out where to store it. Please specify a " + "subject id with {sid[flag]}, or an absolute aseg path with {aseg_name[flag]}." ) + LOGGER.info(message.format(**subjects.flags)) # Run QC check LOGGER.info("Running volume-based QC check on segmentation...") seg_voxvol = np.prod(orig_img.header.get_zooms()) if not check_volume(pred_data, seg_voxvol): - LOGGER.warning( - "Total segmentation volume is too small. Segmentation may be " - "corrupted." - ) + LOGGER.warning("Total segmentation volume is too small. Segmentation may be corrupted.") if qc_file_handle is not None: qc_file_handle.write(subject.id + "\n") qc_file_handle.flush() diff --git a/FastSurferCNN/utils/common.py b/FastSurferCNN/utils/common.py index 747d55715..2e39f7840 100644 --- a/FastSurferCNN/utils/common.py +++ b/FastSurferCNN/utils/common.py @@ -997,9 +997,7 @@ def __getitem__(self, item: int | str) -> SubjectDirectory: """ if isinstance(item, int): if item < 0 or item >= self._num_subjects: - raise IndexError( - f"The index {item} is out of bounds for the subject list." - ) + raise IndexError(f"The index {item} is out of bounds for the subject list.") # subject is always an absolute path (or relative to the working directory) # ... of the input file diff --git a/FastSurferCNN/utils/logging.py b/FastSurferCNN/utils/logging.py index 195e8bf52..587d4af74 100644 --- a/FastSurferCNN/utils/logging.py +++ b/FastSurferCNN/utils/logging.py @@ -13,8 +13,10 @@ # limitations under the License. # IMPORTS +import logging as _logging from logging import CRITICAL, DEBUG, ERROR, INFO, WARNING, FileHandler, Logger, StreamHandler, basicConfig, getLogger from logging import getLogger as get_logger +from os import environ as _environ from pathlib import Path as _Path from sys import stdout as _stdout @@ -39,4 +41,8 @@ def setup_logging(log_file_path: _Path | str): handlers.append(FileHandler(filename=log_file_path, mode="a")) - basicConfig(level=INFO, format=_FORMAT, handlers=handlers) + log_level = _environ.get("FASTSURFER_LOG_LEVEL", "INFO").upper() + if log_level not in ("INFO", "DEBUG", "WARNING", "WARN", "ERROR", "CRITICAL", "FATAL"): + raise RuntimeError(f"Invalid log level: {log_level}") + + basicConfig(level=getattr(_logging, log_level), format=_FORMAT, handlers=handlers) diff --git a/FastSurferCNN/utils/parser_defaults.py b/FastSurferCNN/utils/parser_defaults.py index 08e904b0a..6b8b977db 100644 --- a/FastSurferCNN/utils/parser_defaults.py +++ b/FastSurferCNN/utils/parser_defaults.py @@ -21,9 +21,9 @@ >>> allows_root = args.root # instead of the default dest args.allow_root Values can also be extracted by ->>> print(ALL_FLAGS["allow_root"](dict, dest="root") ->>> # {'flag': '--allow_root', 'flags': ('--allow_root',), 'action': 'store_true', ->>> # 'dest': 'root', 'help': 'Allow execution as root user.'} +>>> print(ALL_FLAGS["seg_log"](dict, dest="log_file") +>>> # {'flag': '--seg_log', 'flags': ('--seg_log',), 'type': str, default='', 'dest': 'log_file', +>>> # 'help': 'Absolute path to file in which run logs will be saved. If not set, logs will not be saved.'} """ import argparse @@ -274,12 +274,7 @@ class SubjectDirectoryConfig: fieldname="search_tag", ), "csv_file": __arg("--csv_file", dc=SubjectDirectoryConfig), - "batch_size": __arg( - "--batch_size", - type=int, - default=1, - help="Batch size for inference. Default=1" - ), + "batch_size": __arg("--batch_size", type=int, default=1, help="Batch size for inference. Default=1"), "sd": __arg("--sd", dc=SubjectDirectoryConfig, fieldname="out_dir"), "qc_log": __arg( "--qc_log", diff --git a/Tutorial/Complete_FastSurfer_Tutorial.ipynb b/Tutorial/Complete_FastSurfer_Tutorial.ipynb index d2d1249a6..94c8cd00f 100644 --- a/Tutorial/Complete_FastSurfer_Tutorial.ipynb +++ b/Tutorial/Complete_FastSurfer_Tutorial.ipynb @@ -2296,10 +2296,10 @@ " (smallest per-direction voxel size) in the T1w\n", " image:\n", " If the minimal voxel size is bigger than 0.98mm,\n", - " the image is conformed to 1mm isometric.\n", + " the image is conformed to 1mm isotropic.\n", " If the minimal voxel size is smaller or equal to\n", " 0.98mm, the T1w image will be conformed to\n", - " isometric voxels of that voxel size.\n", + " isotropic voxels of that voxel size.\n", " The voxel size (whether set manually or derived)\n", " determines whether the surfaces are processed with\n", " highres options (below 1mm) or not.\n", diff --git a/doc/overview/FLAGS.md b/doc/overview/FLAGS.md index 8b1529afc..3f06d74dd 100644 --- a/doc/overview/FLAGS.md +++ b/doc/overview/FLAGS.md @@ -47,8 +47,8 @@ In the following, we give an overview of the most important options. You can vie * `--threads`, `--threads_seg` and `--threads_surf`: Target number of threads for all modules, segmentation, and surface pipeline. The default (`1`) tells FastSurfer to only use one core. Note, that the default value may change in the future for better performance on multi-core architectures. If threads for surface reconstruction is greater than 1, both hemispheres are processed in parallel with half the threads allocated to each hemisphere. * `--vox_size`: Forces processing at a specific voxel size. If a number between 0.7 and 1 is specified (below is experimental) the T1w image is conformed to that isotropic voxel size and processed. If "min" is specified (default), the voxel size is read from the size of the minimal voxel size (smallest per-direction voxel size) in the T1w image: - If the minimal voxel size is bigger than 0.98mm, the image is conformed to 1mm isometric. - If the minimal voxel size is smaller or equal to 0.98mm, the T1w image will be conformed to isometric voxels of that voxel size. + If the minimal voxel size is bigger than 0.98mm, the image is conformed to 1mm isotropic. + If the minimal voxel size is smaller or equal to 0.98mm, the T1w image will be conformed to isotropic voxels of that voxel size. The voxel size (whether set manually or derived) determines whether the surfaces are processed with highres options (below 1mm) or not. * `--py`: Command for python, used in both pipelines. Default: python3.10 * `--conformed_name`: Name of the file in which the conformed input image will be saved. Default location: \$SUBJECTS_DIR/\$sid/mri/orig.mgz diff --git a/recon_surf/long_prepare_template.sh b/recon_surf/long_prepare_template.sh index 0b23e710e..6efce5971 100755 --- a/recon_surf/long_prepare_template.sh +++ b/recon_surf/long_prepare_template.sh @@ -108,10 +108,10 @@ FLAGS: (smallest per-direction voxel size) in the T1w image: If the minimal voxel size is bigger than 0.98mm, - the image is conformed to 1mm isometric. + the image is conformed to 1mm isotropic. If the minimal voxel size is smaller or equal to 0.98mm, the T1w image will be conformed to - isometric voxels of that voxel size. + isotropic voxels of that voxel size. The voxel size (whether set manually or derived) determines whether the surfaces are processed with highres options (below 1mm) or not. diff --git a/recon_surf/recon-surf.sh b/recon_surf/recon-surf.sh index 18a091892..c3a925e0f 100755 --- a/recon_surf/recon-surf.sh +++ b/recon_surf/recon-surf.sh @@ -222,14 +222,14 @@ echo "T1 $t1" echo "asegdkt_segfile $asegdkt_segfile" echo "" -if [ -z "$SUBJECTS_DIR" ] +if [[ -z "$SUBJECTS_DIR" ]] then echo "ERROR: \$SUBJECTS_DIR not set. Either set it via the shell prior to" echo " running recon-surf.sh or supply it via the --sd flag." exit 1 fi -if [ -z "$FREESURFER_HOME" ] +if [[ -z "$FREESURFER_HOME" ]] then echo "ERROR: Did not find \$FREESURFER_HOME. A working version of FreeSurfer $FS_VERSION_SUPPORT" echo " is needed to run recon-surf locally." @@ -241,48 +241,40 @@ fi # needed in FS72 due to a bug in recon-all --fill using FREESURFER instead of FREESURFER_HOME export FREESURFER=$FREESURFER_HOME -if [ "$check_version" == "1" ] +if [[ "$check_version" == "1" ]] && grep -q -v "${FS_VERSION_SUPPORT}" "$FREESURFER_HOME/build-stamp.txt" then - if grep -q -v "${FS_VERSION_SUPPORT}" "$FREESURFER_HOME/build-stamp.txt" - then - echo "ERROR: You are trying to run recon-surf with FreeSurfer version $(cat "$FREESURFER_HOME/build-stamp.txt")." - echo " We are currently supporting only FreeSurfer $FS_VERSION_SUPPORT." - echo " Therefore, make sure to export and source the correct FreeSurfer version" - echo " before running recon-surf.sh: " - echo " export FREESURFER_HOME=/path/to/your/local/fs$FS_VERSION_SUPPORT" - echo " source \$FREESURFER_HOME/SetUpFreeSurfer.sh" - exit 1 - fi + echo "ERROR: You are trying to run recon-surf with FreeSurfer version $(cat "$FREESURFER_HOME/build-stamp.txt")." + echo " We are currently supporting only FreeSurfer $FS_VERSION_SUPPORT." + echo " Therefore, make sure to export and source the correct FreeSurfer version" + echo " before running recon-surf.sh: " + echo " export FREESURFER_HOME=/path/to/your/local/fs$FS_VERSION_SUPPORT" + echo " source \$FREESURFER_HOME/SetUpFreeSurfer.sh" + exit 1 fi -if [ -z "$PYTHONUNBUFFERED" ] -then - export PYTHONUNBUFFERED=0 -fi +if [[ -z "$PYTHONUNBUFFERED" ]] ; then export PYTHONUNBUFFERED=0 ; fi if [[ "$long" == "1" ]] && [[ "$base" == "1" ]] then echo "ERROR: You specified both --long and --base. You need to setup and then run base template first," echo "before you can run any longitudinal time points." - exit 1; + exit 1 fi -if [[ "$base" == "1" ]] +if [[ "$base" == "1" ]] && [[ ! -f "$SUBJECTS_DIR/$subject/base-tps.fastsurfer" ]] then - if [ ! -f "$SUBJECTS_DIR/$subject/base-tps.fastsurfer" ] ; then - echo "ERROR: $subject is either not found in SUBJECTS_DIR" - echo "or it is not a longitudinal template directory (base)," - echo "which needs to contain base-tps.fastsurfer file. Please ensure that" - echo "the base (template) has been created with long_prepare_template.sh." - exit 1 - fi + echo "ERROR: $subject is either not found in SUBJECTS_DIR" + echo "or it is not a longitudinal template directory (base)," + echo "which needs to contain base-tps.fastsurfer file. Please ensure that" + echo "the base (template) has been created with long_prepare_template.sh." + exit 1 fi basedir="" -if [ "$long" == "1" ] +if [[ "$long" == "1" ]] then basedir="$SUBJECTS_DIR/$baseid" - if [ ! -f "$basedir/base-tps.fastsurfer" ] ; then + if [[ ! -f "$basedir/base-tps.fastsurfer" ]] ; then echo "ERROR: $baseid is either not found in \$SUBJECTS_DIR or it is not a longitudinal" echo " template directory, which needs to contain base-tps.fastsurfer file. Please" echo " ensure that the base (template) has been created when running with --long flag." @@ -295,7 +287,7 @@ then fi fi -if [ -z "$t1" ] || [ ! -f "$t1" ] +if [[ -z "$t1" ]] || [[ ! -f "$t1" ]] then echo "ERROR: T1 image ($t1) could not be found. Must supply an existing T1 input" echo " (conformed, full head) via --t1 (absolute path and name)." @@ -303,19 +295,19 @@ then exit 1 fi -if [ -z "$subject" ] +if [[ -z "$subject" ]] then echo "ERROR: must supply subject name via --sid" exit 1 fi -if [ -z "$asegdkt_segfile" ] +if [[ -z "$asegdkt_segfile" ]] then # Set to default asegdkt_segfile="${SUBJECTS_DIR}/${subject}/mri/aparc.DKTatlas+aseg.deep.mgz" fi -if [ ! -f "$asegdkt_segfile" ] +if [[ ! -f "$asegdkt_segfile" ]] then # No segmentation found, exit with error echo "ERROR: Segmentation ($asegdkt_segfile) could not be found! " @@ -337,11 +329,9 @@ export OMP_NUM_THREADS=$threads export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$threads # define the fsthreads variable for the joint section -if [ "$threads" -gt "1" ] ; then fsthreads="-threads $threads -itkthreads $threads" -else fsthreads="" -fi +if [[ "$threads" -gt "1" ]] ; then fsthreads="-threads $threads -itkthreads $threads" ; else fsthreads="" ; fi -if [ "$(echo -n "${SUBJECTS_DIR}/${subject}" | wc -m)" -gt 185 ] +if [[ "$(echo -n "${SUBJECTS_DIR}/${subject}" | wc -m)" -gt 185 ]] then echo "ERROR: Subject directory path is very long." echo " This is known to cause errors due to some commands run by freesurfer versions built for Ubuntu." @@ -350,7 +340,7 @@ then fi # Check if running on an existing subject directory -if [ -f "$SUBJECTS_DIR/$subject/mri/wm.mgz" ] || [ -f "$SUBJECTS_DIR/$subject/mri/aparc.DKTatlas+aseg.orig.mgz" ] +if [[ -f "$SUBJECTS_DIR/$subject/mri/wm.mgz" ]] || [[ -f "$SUBJECTS_DIR/$subject/mri/aparc.DKTatlas+aseg.orig.mgz" ]] then on_existing_run="true" if [[ "$edits" == "true" ]] @@ -368,13 +358,6 @@ fi # collect info StartTime=$(date) tSecStart=$(date '+%s') -# unused -# year=$(date +%Y) -# month=$(date +%m) -# day=$(date +%d) -# hour=$(date +%H) -# min=$(date +%M) - # Setup dirs mkdir -p "$SUBJECTS_DIR/$subject/scripts" @@ -396,9 +379,9 @@ fi # Set up log file DoneFile="$SUBJECTS_DIR/$subject/scripts/recon-surf.done" -if [ "$DoneFile" != /dev/null ] ; then rm -f "$DoneFile" ; fi +if [[ "$DoneFile" != /dev/null ]] ; then rm -f "$DoneFile" ; fi LF="$SUBJECTS_DIR/$subject/scripts/recon-surf.log" -if [ "$LF" != /dev/null ] && [[ "$edits" != "true" ]]; then rm -f "$LF" ; fi +if [[ "$LF" != /dev/null ]] && [[ "$edits" != "true" ]]; then rm -f "$LF" ; fi echo "Log file for recon-surf.sh" >> "$LF" { # all output tee -a "$LF" date 2>&1 @@ -415,18 +398,16 @@ echo "Log file for recon-surf.sh" >> "$LF" echo "Running on top of an existing subject directory with edits=$edits!" fi echo " " - if [ "$base" == "1" ] ; then - echo " " + if [[ "$base" == "1" ]] ; then echo "================== BASE - Longitudinal Template Creation =========================" echo " " - elif [ "$long" == "1" ] ; then - echo " " + elif [[ "$long" == "1" ]] ; then echo "================== LONG - Longitudinal Timpe Point Creation ======================" echo "long: using template directory (base) $baseid" echo " " fi # Print parallelization parameters - if [ "$DoParallel" == "1" ] + if [[ "$DoParallel" == "1" ]] then echo " RUNNING both hemis in PARALLEL" else @@ -442,7 +423,6 @@ echo "Log file for recon-surf.sh" >> "$LF" cmd="$python $FASTSURFER_HOME/FastSurferCNN/quick_qc.py --asegdkt_segfile $asegdkt_segfile" RunIt "$cmd" "$LF" -echo "" | tee -a "$LF" ########################################## START ######################################################## @@ -504,14 +484,14 @@ popd > /dev/null || ( echo "Could not change to subject_dir" ; exit 1 ) # ============================= MASK & ASEG_noCC ======================================== -if [ "$long" == "1" ] ; then +if [[ "$long" == "1" ]] ; then # for long we copy mask from base cmda=(cp "$basedir/mri/mask.mgz" "$mask") run_it "$LF" "${cmda[@]}" fi aseg_nocc="aseg.auto_noCCseg.mgz" -if [ ! -f "$mask" ] || [ ! -f "$mdir/$aseg_nocc" ] ; then +if [[ ! -f "$mask" ]] || [[ ! -f "$mdir/$aseg_nocc" ]] ; then # independently of the existence of manedit files, generate the baseline files. # Mask or aseg.auto_noCCseg not found; create them from aparc.DKTatlas+aseg { @@ -528,11 +508,11 @@ if [ ! -f "$mask" ] || [ ! -f "$mdir/$aseg_nocc" ] ; then cmda=($python "$FASTSURFER_HOME/FastSurferCNN/reduce_to_aseg.py" -i "$mdir/aparc.DKTatlas+aseg.orig.mgz" -o "$mdir/$aseg_nocc" --fixwm) - if [ "$base" == "1" ] && [ ! -f "$mask" ] ; then + if [[ "$base" == "1" ]] && [[ ! -f "$mask" ]] ; then # for base we build union of mapped masks beforehand so it should be available echo "ERROR: $mask missing, but base run requires $mask!" | tee -a "$LF" exit 1 - elif [ "$long" != "1" ] && [ "$base" != 1 ] ; then + elif [[ "$long" != "1" ]] && [[ "$base" != 1 ]] ; then # cross-sectional processing, add outmask to cmd (not for or base long stream) cmda+=(--outmask "$mask") fi @@ -550,7 +530,7 @@ fi # ============================= NU BIAS CORRECTION ======================================= -if [ ! -f "$mdir/orig_nu.mgz" ] ; then +if [[ ! -f "$mdir/orig_nu.mgz" ]] ; then # only run the bias field correction, if the bias field corrected does not exist already { echo " " @@ -608,7 +588,7 @@ fi # create norm by masking nu (supports manedit-ed mask) cmda=(mri_mask "$mdir/nu.mgz" "$mask" "$mdir/norm.mgz") run_it "$LF" "${cmda[@]}" -if [ "$get_t1" == "1" ] +if [[ "$get_t1" == "1" ]] then # create T1.mgz from nu (!! here we could also try passing aseg?) # T1.mgz was needed by some 3rd party downstream tools such as fmriprep, so we provide it @@ -662,7 +642,7 @@ RunIt "$cmd" "$LF" echo " " } | tee -a "$LF" -if [ "$long" == "1" ] ; then +if [[ "$long" == "1" ]] ; then # in long we can skip fill as surfaces come from base # it would be great to also skip WM, but it is needed in place_surface to clip bright # maybe later add code to copy edits from base in maskbfs and wm segmentation, currently not supported! @@ -689,7 +669,7 @@ export OMP_NUM_THREADS=$threads_hemi export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$threads_hemi # define the fsthreads variable for the joint section -if [ "$threads_hemi" -gt "1" ] ; then fsthreads="-threads $threads_hemi -itkthreads $threads_hemi" +if [[ "$threads_hemi" -gt "1" ]] ; then fsthreads="-threads $threads_hemi -itkthreads $threads_hemi" else fsthreads="" fi @@ -709,66 +689,62 @@ for hemi in lh rh ; do # ============================= TESSELLATE - SMOOTH ===================================================== # In Long stream we skip these - if [ "$long" == "0" ] ; then - - { - echo "echo \" \"" - echo "echo \"================== Creating surfaces $hemi - orig.nofix ==================\"" - echo "echo \" \"" - } | tee -a "$CMDF" - - if [ "$fstess" == "1" ] + if [[ "$long" == "0" ]] then - cmd="recon-all -subject $subject -hemi $hemi -tessellate -smooth1 -no-isrunning $hiresflag $fsthreads" - RunIt "$cmd" "$LF" "$CMDF" - else - # instead of mri_tesselate lego land use marching cube - if [ $hemi == "lh" ] ; then - hemivalue=255 + { + echo "echo \" \"" + echo "echo \"================== Creating surfaces $hemi - orig.nofix ==================\"" + echo "echo \" \"" + } | tee -a "$CMDF" + + if [[ "$fstess" == "1" ]] + then + cmd="recon-all -subject $subject -hemi $hemi -tessellate -smooth1 -no-isrunning $hiresflag $fsthreads" + RunIt "$cmd" "$LF" "$CMDF" else - hemivalue=127 - fi + # instead of mri_tesselate lego land use marching cube + if [[ $hemi == "lh" ]] ; then hemivalue=255 ; else hemivalue=127 ; fi - # extract initial surface "?h.orig.nofix" - cmd="mri_pretess $mdir/filled.mgz $hemivalue $mdir/brain.mgz $mdir/filled-pretess$hemivalue.mgz" - RunIt "$cmd" "$LF" "$CMDF" + # extract initial surface "?h.orig.nofix" + cmd="mri_pretess $mdir/filled.mgz $hemivalue $mdir/brain.mgz $mdir/filled-pretess$hemivalue.mgz" + RunIt "$cmd" "$LF" "$CMDF" - # Marching cube does not return filename and wrong volume info! - outmesh=$sdir/$hemi.orig.nofix$hires_surface_suffix - cmd="mri_mc $mdir/filled-pretess$hemivalue.mgz $hemivalue $outmesh" - RunIt "$cmd" "$LF" "$CMDF" + # Marching cube does not return filename and wrong volume info! + outmesh=$sdir/$hemi.orig.nofix$hires_surface_suffix + cmd="mri_mc $mdir/filled-pretess$hemivalue.mgz $hemivalue $outmesh" + RunIt "$cmd" "$LF" "$CMDF" - # Rewrite surface orig.nofix to fix vertex locs bug (scannerRAS instead of surfaceRAS set with mc) - #cmd="$python ${binpath}rewrite_mc_surface.py --input $outmesh --output $outmesh --filename_pretess $mdir/filled-pretess$hemivalue.mgz" - #RunIt "$cmd" "$LF" "$CMDF" + # Rewrite surface orig.nofix to fix vertex locs bug (scannerRAS instead of surfaceRAS set with mc) + #cmd="$python ${binpath}rewrite_mc_surface.py --input $outmesh --output $outmesh --filename_pretess $mdir/filled-pretess$hemivalue.mgz" + #RunIt "$cmd" "$LF" "$CMDF" - # Check if the surfaceRAS was correctly set and exit otherwise (sanity check in case nibabel changes their default header behaviour) - { - cmd="mris_info $outmesh | tr -s ' ' | grep -q 'vertex locs : surfaceRAS'" - echo "echo \"$cmd\"" - echo "$timecmd $cmd" - } | tee -a "$CMDF" - echo "if [ \${PIPESTATUS[1]} -ne 0 ] ; then echo \"Incorrect header information detected in $outmesh: vertex locs is not set to surfaceRAS. Exiting... \" ; exit 1 ; fi" >> "$CMDF" + # Check if the surfaceRAS was correctly set and exit otherwise (sanity check in case nibabel changes their default header behaviour) + { + cmd="mris_info $outmesh | tr -s ' ' | grep -q 'vertex locs : surfaceRAS'" + echo "echo \"$cmd\"" + echo "$timecmd $cmd" + } | tee -a "$CMDF" + echo "if [ \${PIPESTATUS[1]} -ne 0 ] ; then echo \"Incorrect header information detected in $outmesh: vertex locs is not set to surfaceRAS. Exiting... \" ; exit 1 ; fi" >> "$CMDF" - # Reduce to largest component (usually there should only be one) - cmd="mris_extract_main_component $outmesh $outmesh" - RunIt "$cmd" "$LF" "$CMDF" - - # for hires decimate mesh - if [ -n "$hiresflag" ] ; then - DecimationFaceArea="0.5" - # Reduce the number of faces such that the average face area is - # DecimationFaceArea. If the average face area is already more - # than DecimationFaceArea, then the surface is not changed. - # set cmd = (mris_decimate -a $DecimationFaceArea ../surf/$hemi.orig.nofix.predec ../surf/$hemi.orig.nofix) - cmd="mris_remesh --desired-face-area $DecimationFaceArea --input $outmesh --output $sdir/$hemi.orig.nofix" + # Reduce to largest component (usually there should only be one) + cmd="mris_extract_main_component $outmesh $outmesh" + RunIt "$cmd" "$LF" "$CMDF" + + # for hires decimate mesh + if [[ -n "$hiresflag" ]] + then + DecimationFaceArea="0.5" + # Reduce the number of faces such that the average face area is DecimationFaceArea. If the average face + # area is already more than DecimationFaceArea, then the surface is not changed. + # set cmd = (mris_decimate -a $DecimationFaceArea ../surf/$hemi.orig.nofix.predec ../surf/$hemi.orig.nofix) + cmd="mris_remesh --desired-face-area $DecimationFaceArea --input $outmesh --output $sdir/$hemi.orig.nofix" + RunIt "$cmd" "$LF" "$CMDF" + fi + # -smooth1 (explicitly state 10 iteration (default) but may change in future) + cmd="mris_smooth -n 10 -nw -seed 1234 $sdir/$hemi.orig.nofix $sdir/$hemi.smoothwm.nofix" RunIt "$cmd" "$LF" "$CMDF" fi - # -smooth1 (explicitly state 10 iteration (default) but may change in future) - cmd="mris_smooth -n 10 -nw -seed 1234 $sdir/$hemi.orig.nofix $sdir/$hemi.smoothwm.nofix" - RunIt "$cmd" "$LF" "$CMDF" - fi else # LONG @@ -786,39 +762,41 @@ for hemi in lh rh ; do # ============================= INFLATE1 - QSPHERE ===================================================== # In Long stream we skip these - if [ "$long" == "0" ] ; then - - { - echo "echo \"\"" - echo "echo \"=================== Creating surfaces $hemi - qsphere ====================\"" - echo "echo \"\"" - } | tee -a "$CMDF" + if [[ "$long" == "0" ]] + then - #surface inflation (54sec both hemis) (needed for qsphere and for topo-fixer) - cmd="recon-all -subject $subject -hemi $hemi -inflate1 -no-isrunning $hiresflag $fsthreads" - RunIt "$cmd" "$LF" "$CMDF" + { + echo "echo \"\"" + echo "echo \"=================== Creating surfaces $hemi - qsphere ====================\"" + echo "echo \"\"" + } | tee -a "$CMDF" - if [ "$fsqsphere" == "1" ] - then - # quick spherical mapping (2min48sec) - cmd="recon-all -subject $subject -hemi $hemi -qsphere -no-isrunning $hiresflag $fsthreads" - RunIt "$cmd" "$LF" "$CMDF" - else - # instead of mris_sphere, directly project to sphere with spectral approach - # equivalent to -qsphere - # (23sec) - cmd="$python ${binpath}spherically_project_wrapper.py --hemi $hemi --sdir $sdir" - printf -v tmp %q "$python" - cmd="$cmd --subject $subject --threads=$threads_hemi --py ${tmp} --binpath ${binpath}" + #surface inflation (54sec both hemis) (needed for qsphere and for topo-fixer) + cmd="recon-all -subject $subject -hemi $hemi -inflate1 -no-isrunning $hiresflag $fsthreads" RunIt "$cmd" "$LF" "$CMDF" - fi + + if [ "$fsqsphere" == "1" ] + then + # quick spherical mapping (2min48sec) + cmd="recon-all -subject $subject -hemi $hemi -qsphere -no-isrunning $hiresflag $fsthreads" + RunIt "$cmd" "$LF" "$CMDF" + else + # instead of mris_sphere, directly project to sphere with spectral approach + # equivalent to -qsphere + # (23sec) + cmd="$python ${binpath}spherically_project_wrapper.py --hemi $hemi --sdir $sdir" + printf -v tmp %q "$python" + cmd="$cmd --subject $subject --threads=$threads_hemi --py ${tmp} --binpath ${binpath}" + RunIt "$cmd" "$LF" "$CMDF" + fi fi # not long # ============================= FIX - WHITEPREAPARC ================================================== # In Long stream we skip topo fix - if [ "$long" == "0" ] ; then + if [ "$long" == "0" ] + then # longitudinal base and cross-sectional { @@ -894,47 +872,49 @@ for hemi in lh rh ; do # ============================= SPHERE - SURFREG (optional) ============================================== # if we segment with FS or if surface registration is requested do it here: - if [ "$fsaparc" == "1" ] || [ "$fssurfreg" == "1" ] ; then + if [[ "$fsaparc" == "1" ]] || [[ "$fssurfreg" == "1" ]] + then { echo "echo \" \"" echo "echo \"============ Creating surfaces $hemi - FS sphere, surfreg ===============\"" echo "echo \" \"" } | tee -a "$CMDF" - if [ "$long" == "0" ] ; then + if [[ "$long" == "0" ]] + then - # SPHERE: Inflate to sphere with minimal metric distortion - cmd="recon-all -subject $subject -hemi $hemi -sphere $hiresflag -no-isrunning $fsthreads" - RunIt "$cmd" "$LF" "$CMDF" + # SPHERE: Inflate to sphere with minimal metric distortion + cmd="recon-all -subject $subject -hemi $hemi -sphere $hiresflag -no-isrunning $fsthreads" + RunIt "$cmd" "$LF" "$CMDF" - # SURFREG (sphere.reg) - # Surface registration for cross-subject correspondence (registration to fsaverage) - # (mr) FIX: sometimes FreeSurfer Sphere Reg. fails and moves pre and post central - # one gyrus too far posterior, FastSurferCNN's image-based segmentation does not - # seem to do this, so we initialize the spherical registration with the better - # cortical segmentation from FastSurferCNN, this replaces recon-all -surfreg - # 1. get alpha, beta, gamma for global alignment (rotation) based on aseg centers - # (note the former fix, initializing with pre-central label, is not working in FS7.2 - # as they broke the label initialization in mris_register) - cmd="$python ${binpath}/rotate_sphere.py \ - --srcsphere $sdir/${hemi}.sphere \ - --srcaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot \ - --trgsphere $FREESURFER_HOME/subjects/fsaverage/surf/${hemi}.sphere \ - --trgaparc $FREESURFER_HOME/subjects/fsaverage/label/${hemi}.aparc.annot \ - --out $sdir/${hemi}.angles.txt" - RunIt "$cmd" "$LF" "$CMDF" - # 2. use global rotation as initialization to non-linear registration: - cmd="mris_register -curv -norot -rotate \`cat $sdir/${hemi}.angles.txt\` \ - $sdir/${hemi}.sphere \ - $FREESURFER_HOME/average/${hemi}.folding.atlas.acfb40.noaparc.i12.2016-08-02.tif \ - $sdir/${hemi}.sphere.reg" - RunIt "$cmd" "$LF" "$CMDF" - # command to generate new aparc to check if registration was OK - # run only for debugging - # cmd="mris_ca_label -l $SUBJECTS_DIR/$subject/label/${hemi}.cortex.label \ - # -aseg $SUBJECTS_DIR/$subject/mri/aseg.presurf.mgz \ - # -seed 1234 $subject $hemi $SUBJECTS_DIR/$subject/surf/${hemi}.sphere.reg \ - # $SUBJECTS_DIR/$subject/label/${hemi}.aparc.DKTatlas-guided.annot" + # SURFREG (sphere.reg) + # Surface registration for cross-subject correspondence (registration to fsaverage) + # (mr) FIX: sometimes FreeSurfer Sphere Reg. fails and moves pre and post central + # one gyrus too far posterior, FastSurferCNN's image-based segmentation does not + # seem to do this, so we initialize the spherical registration with the better + # cortical segmentation from FastSurferCNN, this replaces recon-all -surfreg + # 1. get alpha, beta, gamma for global alignment (rotation) based on aseg centers + # (note the former fix, initializing with pre-central label, is not working in FS7.2 + # as they broke the label initialization in mris_register) + cmd="$python ${binpath}/rotate_sphere.py \ + --srcsphere $sdir/${hemi}.sphere \ + --srcaparc $ldir/$hemi.aparc.DKTatlas.mapped.annot \ + --trgsphere $FREESURFER_HOME/subjects/fsaverage/surf/${hemi}.sphere \ + --trgaparc $FREESURFER_HOME/subjects/fsaverage/label/${hemi}.aparc.annot \ + --out $sdir/${hemi}.angles.txt" + RunIt "$cmd" "$LF" "$CMDF" + # 2. use global rotation as initialization to non-linear registration: + cmd="mris_register -curv -norot -rotate \`cat $sdir/${hemi}.angles.txt\` \ + $sdir/${hemi}.sphere \ + $FREESURFER_HOME/average/${hemi}.folding.atlas.acfb40.noaparc.i12.2016-08-02.tif \ + $sdir/${hemi}.sphere.reg" + RunIt "$cmd" "$LF" "$CMDF" + # command to generate new aparc to check if registration was OK + # run only for debugging + # cmd="mris_ca_label -l $SUBJECTS_DIR/$subject/label/${hemi}.cortex.label \ + # -aseg $SUBJECTS_DIR/$subject/mri/aseg.presurf.mgz \ + # -seed 1234 $subject $hemi $SUBJECTS_DIR/$subject/surf/${hemi}.sphere.reg \ + # $SUBJECTS_DIR/$subject/label/${hemi}.aparc.DKTatlas-guided.annot" else # longitudinal @@ -970,7 +950,8 @@ for hemi in lh rh ; do # it is then used also below for surface placement. # we should consider, always computing it (when surfreg is available) -> test later what consequences this has #if [ "$fsaparc" == "1" ] || [ "$fssurfreg" == "1" ] ; then - if [ "$fsaparc" == "1" ] ; then + if [[ "$fsaparc" == "1" ]] + then { echo "echo \" \"" echo "echo \"============ Creating surfaces $hemi - FS aparc ===============\"" @@ -978,7 +959,8 @@ for hemi in lh rh ; do } | tee -a "$CMDF" longflag="" - if [ "$long" == "1" ] ; then + if [[ "$long" == "1" ]] + then # recon-all has different treatment for cortparc: # initialize with aparc.annot from base longflag="-long -R $basedir/label/${hemi}.aparc.annot" @@ -995,7 +977,8 @@ for hemi in lh rh ; do # first select what cortical parcellation to use to guide surface placement: aparc="" - if [ "$fsaparc" == "1" ] ; then + if [[ "$fsaparc" == "1" ]] + then { echo "echo \" \"" echo "echo \"============ Creating surfaces $hemi - white and pial using FS aparc ===============\"" @@ -1023,10 +1006,8 @@ for hemi in lh rh ; do --threads $threads_hemi --wm wm.mgz --invol brain.finalsurfs.mgz --$hemi --o ../surf/${hemi}.white \ --white --nsmooth 0 --rip-label ../label/${hemi}.cortex.label \ --rip-bg --rip-surf ../surf/${hemi}.white.preaparc --aparc $aparc" - if [ "$long" == "0" ] ; then # cross/regular/base - cmd="$cmd --i ../surf/$hemi.white.preaparc" - else # longitudinal processing ; also adds longmaxdist - cmd="$cmd --i ../surf/$hemi.orig_white --max-cbv-dist 3.5" + if [[ "$long" == "0" ]] ; then cmd="$cmd --i ../surf/$hemi.white.preaparc" # cross/regular/base + else cmd="$cmd --i ../surf/$hemi.orig_white --max-cbv-dist 3.5" # longitudinal processing ; also adds longmaxdist fi RunIt "$cmd" "$LF" "$CMDF" @@ -1037,8 +1018,7 @@ for hemi in lh rh ; do --pial --nsmooth 0 --rip-label ../label/${hemi}.cortex+hipamyg.label \ --pin-medial-wall ../label/${hemi}.cortex.label --aparc $aparc \ --repulse-surf ../surf/${hemi}.white --white-surf ../surf/${hemi}.white" - if [ "$long" == "0" ] ; then # cross/regular/base - cmd="$cmd --i ../surf/$hemi.white" + if [ "$long" == "0" ] ; then cmd="$cmd --i ../surf/$hemi.white" # cross/regular/base else # longitudinal processing ; also adds longmaxdist cmd="$cmd --i ../surf/$hemi.orig_pial --max-cbv-dist 3.5 --blend-surf .25 ../surf/$hemi.white" fi @@ -1074,10 +1054,8 @@ for hemi in lh rh ; do cmd="recon-all -subject $subject -hemi $hemi -curvstats -no-isrunning $hiresflag $fsthreads" RunIt "$cmd" "$LF" "$CMDF" - - - - if [ "$DoParallel" == "0" ] ; then + if [[ "$DoParallel" == "0" ]] + then { echo " " echo " RUNNING $hemi sequentially ... " @@ -1095,12 +1073,9 @@ export OMP_NUM_THREADS=$threads export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$threads # define the fsthreads variable for the joint section (again) -if [ "$threads" -gt "1" ] ; then fsthreads="-threads $threads -itkthreads $threads" -else fsthreads="" -fi +if [[ "$threads" -gt "1" ]] ; then fsthreads="-threads $threads -itkthreads $threads" ; else fsthreads="" ; fi - -if [ "$DoParallel" == 1 ] ; then +if [[ "$DoParallel" == 1 ]] ; then { echo "" echo " RUNNING HEMIs in PARALLEL !!! " @@ -1113,13 +1088,14 @@ fi # ============================= RIBBON =============================================== # Skip RIBBON in base -if [ "$base" != "1" ] ; then +if [[ "$base" != "1" ]] +then -{ - echo "" - echo "============================ Creating surfaces - ribbon ===========================" - echo "" -} | tee -a "$LF" + { + echo "" + echo "============================ Creating surfaces - ribbon ===========================" + echo "" + } | tee -a "$LF" # -cortribbon 4 minutes, ribbon is used in mris_anatomical stats to remove voxels from surface based volumes that should not be cortex # anatomical stats can run without ribbon, but will omit some surface based measures then # wmparc needs ribbon, probably other stuff (aparc to aseg etc). @@ -1131,10 +1107,11 @@ fi # skip in base # ============================= FSAPARC - parc23 surfcon hypo ... ========================================= -if [ "$fsaparc" == "1" ] ; then +if [[ "$fsaparc" == "1" ]] ; then # this per-hemi section does not get parallelized - for hemi in lh rh ; do + for hemi in lh rh + do { echo "" @@ -1144,7 +1121,8 @@ if [ "$fsaparc" == "1" ] ; then # Destrieux Atlas (recon-all -cortparc2): longflag="" - if [ "$long" == "1" ] ; then + if [[ "$long" == "1" ]] + then # recon-all has different treatment for cortparc: # initialize with destrieux annot from base longflag="-long -R $basedir/label/${hemi}.a2009s.annot" @@ -1156,11 +1134,8 @@ if [ "$fsaparc" == "1" ] ; then # DKT Atlas (recon-all -cortparc3): longflag="" - if [ "$long" == "1" ] ; then - # recon-all has different treatment for cortparc: - # initialize with destrieux annot from base - longflag="-long -R $basedir/label/${hemi}.DKTatlas.annot" - fi + # recon-all has different treatment for cortparc: initialize with destrieux annot from base + if [[ "$long" == "1" ]] ; then longflag="-long -R $basedir/label/${hemi}.DKTatlas.annot" ; fi CPAtlas="$FREESURFER_HOME/average/${hemi}.DKTaparc.atlas.acfb40.noaparc.i12.2016-08-02.gcs" annot="$ldir/${hemi}.aparc.DKTatlas.annot" cmd="mris_ca_label -l $ldir/${hemi}.cortex.label -aseg $mdir/aseg.presurf.mgz -seed 1234 $longflag $subject $hemi $sdir/${hemi}.sphere.reg $CPAtlas $annot" @@ -1169,7 +1144,8 @@ if [ "$fsaparc" == "1" ] ; then done # hemi loop # skip in base - if [ "$base" != "1" ] ; then + if [[ "$base" != "1" ]] + then cmd="recon-all -subject $subject -pctsurfcon -hyporelabel -apas2aseg -aparc2aseg -wmparc -parcstats -parcstats2 -parcstats3 $hiresflag $fsthreads" RunIt "$cmd" "$LF" # removed -balabels here and do that below independent of fsaparc flag @@ -1180,7 +1156,8 @@ fi # (FS-APARC) # Skip rest in case we have a base run, we are done here (probably we can skip stuff already in surface creation above) -if [ "$base" != "1" ] ; then +if [[ "$base" != "1" ]] +then # ============================= MAPPED SURF-STATS ========================================= @@ -1191,14 +1168,16 @@ if [ "$base" != "1" ] ; then } | tee -a "$LF" # 2x18sec create stats from mapped aparc - for hemi in lh rh ; do + for hemi in lh rh + do cmd="mris_anatomical_stats -th3 -mgz -cortex $ldir/$hemi.cortex.label -f $statsdir/$hemi.aparc.DKTatlas.mapped.stats -b -a $ldir/$hemi.aparc.DKTatlas.mapped.annot -c $ldir/aparc.annot.mapped.ctab $subject $hemi white" RunIt "$cmd" "$LF" done # ============================= FASTSURFER - surfcon hypo stats ========================================= - if [ "$fsaparc" == "0" ] ; then + if [[ "$fsaparc" == "0" ]] + then { echo "" echo "============= Creating surfaces - pctsurfcon, hypo, segstats ====================" @@ -1238,7 +1217,8 @@ if [ "$base" != "1" ] ; then # get stats for the aseg (note these are surface fine tuned, that may be good or bad, below we also do the stats for the input aseg (plus some processing) # cmd="recon-all -subject $subject -segstats $hiresflag $fsthreads" - if [[ "$segstats_legacy" == "true" ]] ; then + if [[ "$segstats_legacy" == "true" ]] + then cmda=($python "$FASTSURFER_HOME/FastSurferCNN/mri_brainvol_stats.py" --subject "$subject") run_it "$LF" "${cmda[@]}" @@ -1263,11 +1243,8 @@ if [ "$base" != "1" ] ; then "SubCortGray" "TotalGray" "SupraTentorial" "SupraTentorialNotVent" "Mask($mask)" "BrainSegVol-to-eTIV" "MaskVol-to-eTIV") - if [ "$long" == "0" ] ; then - # in long we do not have orig_nofix for surface hole computation as surfaces - # are inherited from base/template - cmda+=("lhSurfaceHoles" "rhSurfaceHoles" "SurfaceHoles") - fi + # in long we do not have orig_nofix for surface hole computation as surfaces are inherited from base/template + if [[ "$long" == "0" ]] ; then cmda+=("lhSurfaceHoles" "rhSurfaceHoles" "SurfaceHoles") ; fi cmda+=("EstimatedTotalIntraCranialVol") run_it "$LF" "${cmda[@]}" @@ -1353,7 +1330,8 @@ if [ "$base" != "1" ] ; then # ============================= FASTSURFER - SYMLINKS ========================================= # Create symlinks for downstream analysis (sub-segmentations, TRACULA, etc.) - if [ "$fsaparc" == "0" ] ; then + if [[ "$fsaparc" == "0" ]] + then # Symlink of aparc.DKTatlas+aseg.mapped.mgz pushd "$mdir" > /dev/null || (echo "Could not cd to $mdir" ; exit 1) softlink_or_copy "aparc.DKTatlas+aseg.mapped.mgz" "aparc.DKTatlas+aseg.mgz" "$LF" @@ -1374,7 +1352,8 @@ if [ "$base" != "1" ] ; then # ============================= BALABELS ========================================= # balabels need sphere.reg - if [ "$fssurfreg" == "1" ] ; then + if [[ "$fssurfreg" == "1" ]] + then # can be produced if surf registration exists #cmd="recon-all -subject $subject -balabels $hiresflag $fsthreads" #RunIt "$cmd" "$LF" @@ -1412,9 +1391,7 @@ tRunHours=$(LC_NUMERIC="en_US.UTF-8" printf %6.3f "$(bc -l <<< "($tSecEnd - $tSe # id -n sends an error message in docker (no user name), fall back to the USER environment variable or username=$(id -un 2>&1) if echo "$username" | grep -q "^id: " ; then - if [[ -n "$USER" ]] ; then username="$USER" - else username="$(id -u)" - fi + if [[ -n "$USER" ]] ; then username="$USER" ; else username="$(id -u)" ; fi fi echo "USER $username" echo "HOST $(hostname)" diff --git a/run_fastsurfer.sh b/run_fastsurfer.sh index 300c4536a..4a55bb162 100755 --- a/run_fastsurfer.sh +++ b/run_fastsurfer.sh @@ -130,10 +130,10 @@ FLAGS: (smallest per-direction voxel size) in the T1w image: If the minimal voxel size is bigger than 0.98mm, - the image is conformed to 1mm isometric. + the image is conformed to 1mm isotropic. If the minimal voxel size is smaller or equal to 0.98mm, the T1w image will be conformed to - isometric voxels of that voxel size. + isotropic voxels of that voxel size. The voxel size (whether set manually or derived) determines whether the surfaces are processed with highres options (below 1mm) or not. @@ -403,9 +403,7 @@ case $key in --py) python="$1" ; shift ;; -h|--help) usage ; exit ;; --version) - if [[ "$#" -lt 1 ]] || [[ "$1" =~ ^-- ]]; then - # no more args or next arg starts with -- - version_and_quit="1" + if [[ "$#" -lt 1 ]] || [[ "$1" =~ ^-- ]]; then version_and_quit="1" # no more args or next arg starts with -- else case "$(echo "$1" | tr '[:upper:]' '[:lower:]')" in all) version_and_quit="+checkpoints+git+pip" ;; @@ -443,10 +441,7 @@ case $key in esac shift # past value ;; - --no_cuda) - echo "WARNING: --no_cuda is deprecated and will be removed, use --device cpu." - device="cpu" - ;; + --no_cuda) echo "WARNING: --no_cuda is deprecated and will be removed, use --device cpu." ; device="cpu" ;; # asegdkt module options #============================================================= @@ -503,10 +498,8 @@ case $key in ############################################################## --base) base=1 ; run_cereb_module="0" ; run_hypvinn_module="0" ; surf_flags=("${surf_flags[@]}" "--base") ;; --long) long=1 ; baseid="$1" ; surf_flags=("${surf_flags[@]}" "--long" "$1") ; shift ;; - *) # unknown option - # if not empty arguments, error & exit - if [[ "$key" != "" ]] ; then echo "ERROR: Flag '$key' unrecognized." ; exit 1 ; fi - ;; + # unknown option ; if not empty arguments, error & exit + *) if [[ "$key" != "" ]] ; then echo "ERROR: Flag '$key' unrecognized." ; exit 1 ; fi ;; esac done set -- "${POSITIONAL[@]}" # restore positional parameters @@ -536,7 +529,7 @@ then version_cache_args=("${version_cache_args[@]}" --sections "$version_and_quit") fi $python "$FASTSURFER_HOME/FastSurferCNN/version.py" "${version_cache_args[@]}" - exit + exit 1 fi source "${reconsurfdir}/functions.sh" @@ -553,7 +546,8 @@ tmpLF=$(mktemp) # CHECKS -if [[ "$legacy_parallel_hemi" == 1 ]] ; then +if [[ "$legacy_parallel_hemi" == 1 ]] +then { echo "WARNING: The --parallel flag is obsolete and will be removed in FastSurfer 3." echo " Hemispheres are now automatically processed in parallel, if threads for surface " @@ -561,7 +555,8 @@ if [[ "$legacy_parallel_hemi" == 1 ]] ; then echo "IMPORTANT NOTE: The threads behavior has also changed, --threads used to define the" echo " number of threads per hemisphere, it now defines the number of threads in total!" } | tee -a "$tmpLF" - if [[ "$threads_surf" == 1 ]] ; then + if [[ "$threads_surf" == 1 ]] + then threads_surf=2 { echo "INFO: We have changed the requested number of threads from 1 to 2, to activate parallel" @@ -911,9 +906,8 @@ then if [[ "$edits" == "true" ]] then { - echo "INFO: $asegdkt_segfile_manedit (manedit file for ) detected," - echo " supersedes $asegdkt_segfile for creation of $aseg_segfile" - echo " and $mask_name!" + echo "INFO: $asegdkt_segfile_manedit (manedit file for ) detected, supersedes" + echo " $asegdkt_segfile for creation of $aseg_segfile and $mask_name!" } | tee -a "$seg_log" asegdkt_segfile="$asegdkt_segfile_manedit" cmd=($python "$fastsurfercnndir/reduce_to_aseg.py" -i "$asegdkt_segfile" -o "$aseg_segfile"