[BUG] Wrong results during the atlas_registration step

[huashengmao]

Describe the bug
For some patients, the MRI images produce snowflake-like artifacts during the atlas_registration step, but no error is raised. Both on Linux and Windows

To Reproduce
Steps to reproduce the behavior:

1. Create anaconda envs python=3.10.19
2. Install 'pytorch==2.4.0 torchvision==0.19.0 torchaudio==2.4.0 pytorch-cuda=12.4'
3. Install 'brainles-preprocessing=0.6.8'
4. Install 'brainles_preprocessing[itk-elastix]'
5. Install 'mkl=2023'
6. edit brainles_preprocessing\registration\elastix\elastix.py from

        parameter_object = itk.ParameterObject.New()
        default_rigid_parameter_map = parameter_object.GetDefaultParameterMap("rigid")
        parameter_object.AddParameterMap(default_rigid_parameter_map)

to

        parameter_object = itk.ParameterObject.New()
        translation_map = parameter_object.GetDefaultParameterMap('translation')
        parameter_object.AddParameterMap(translation_map)
        default_rigid_parameter_map = parameter_object.GetDefaultParameterMap('rigid')
        parameter_object.AddParameterMap(default_rigid_parameter_map)

to avoid errors'Too many samples map outside moving image buffer'

7. Run the code below to preprocess T1, T1c, Flair MRI

from pathlib import Path
from brainles_preprocessing.modality import Modality, CenterModality
from brainles_preprocessing.normalization.percentile_normalizer import (
    PercentileNormalizer,
)
from brainles_preprocessing.preprocessor import AtlasCentricPreprocessor
from brainles_preprocessing.registration import ElastixRegistrator

patient_folder = Path(patient_path)
output_folder = Path(output_path)

percentile_normalizer = PercentileNormalizer(
    lower_percentile=0.1,
    upper_percentile=99.9,
    lower_limit=0,
    upper_limit=1,
)

center = CenterModality(
    modality_name="t1c",  # required
    input_path=patient_folder / "t1ce.nii.gz",  # required
    normalizer=percentile_normalizer,  # optional: required for normalized_* outputs
    raw_bet_output_path=output_folder / "raw_bet/t1c_bet_raw.nii.gz",  # brain extracted
    raw_skull_output_path=output_folder / "raw_skull/t1c_skull_raw.nii.gz",  # with skull
    raw_defaced_output_path=output_folder / "raw_defaced/t1c_defaced_raw.nii.gz",  # defaced
    normalized_bet_output_path=output_folder / "normalized_bet/t1c_bet_normalized.nii.gz",
    normalized_skull_output_path=output_folder / "normalized_skull/t1c_skull_normalized.nii.gz",
    normalized_defaced_output_path=output_folder / "normalized_defaced/t1c_defaced_normalized.nii.gz",
    bet_mask_output_path=output_folder / "masks/t1c_bet_mask.nii.gz",
    defacing_mask_output_path=output_folder / "masks/t1c_defacing_mask.nii.gz",
    atlas_correction=True,  # default: True
    n4_bias_correction=False,  # default: False
)

moving_modalities = [
    Modality(
        modality_name="flair",  # required
        input_path=patient_folder / "t2_flair.nii.gz",  # required
        normalizer=percentile_normalizer,  # optional: required for normalized_* outputs
        raw_bet_output_path=output_folder / "raw_bet/flair_bet_raw.nii.gz",
        raw_skull_output_path=output_folder / "raw_skull/flair_skull_raw.nii.gz",
        raw_defaced_output_path=output_folder / "raw_defaced/flair_defaced_raw.nii.gz",
        normalized_bet_output_path=output_folder / "normalized_bet/flair_bet_normalized.nii.gz",
        normalized_skull_output_path=output_folder / "normalized_skull/flair_skull_normalized.nii.gz",
        normalized_defaced_output_path=output_folder / "normalized_defaced/flair_defaced_normalized.nii.gz",
        atlas_correction=True,  # default: True
        n4_bias_correction=False,  # default: False
    ),
    Modality(
        modality_name="t1",  # required
        input_path=patient_folder / "t1.nii.gz",  # required
        normalizer=percentile_normalizer,  # optional: required for normalized_* outputs
        raw_bet_output_path=output_folder / "raw_bet/t1_bet_raw.nii.gz",
        raw_skull_output_path=output_folder / "raw_skull/t1_skull_raw.nii.gz",
        raw_defaced_output_path=output_folder / "raw_defaced/t1_defaced_raw.nii.gz",
        normalized_bet_output_path=output_folder / "normalized_bet/t1_bet_normalized.nii.gz",
        normalized_skull_output_path=output_folder / "normalized_skull/t1_skull_normalized.nii.gz",
        normalized_defaced_output_path=output_folder / "normalized_defaced/t1_defaced_normalized.nii.gz",
        atlas_correction=True,  # default: True
        n4_bias_correction=False,  # default: False
    )
]

preprocessor = AtlasCentricPreprocessor(
    center_modality=center,
    moving_modalities=moving_modalities,
    registrator=ElastixRegistrator(),
    temp_folder=Path('temp'),
    use_gpu=True,
    limit_cuda_visible_devices='1'
)

preprocessor.run(
    save_dir_coregistration=output_folder / "coregistration",
    save_dir_atlas_registration=output_folder / "atlas_registration",
    save_dir_atlas_correction=output_folder / "atlas_correction",
    save_dir_n4_bias_correction=output_folder / "n4_bias_correction",
    save_dir_brain_extraction=output_folder / "brain_extraction",
    save_dir_defacing=output_folder / "defacing",
)

this code should register center modality to atlas,but it produce wrong results

Screenshots
this is an example of mri that produce snowflake-like artifacts during the atlas_registration step. The first row shows the original T1c, and the FLAIR and T1 images after registration to the registration center; the second row shows the T1c, FLAIR, and T1 images after atlas_registration.

Environment

Ubuntu 20.04.6 LTS and Microsoft Windows [version 10.0.19045.6466]

NVIDIA drivers and GPUs

Linux

+-----------------------------------------------------------------------------------------+
| NVIDIA-SMI 550.107.02             Driver Version: 550.107.02     CUDA Version: 12.4     |
|-----------------------------------------+------------------------+----------------------+
| GPU  Name                 Persistence-M | Bus-Id          Disp.A | Volatile Uncorr. ECC |
| Fan  Temp   Perf          Pwr:Usage/Cap |           Memory-Usage | GPU-Util  Compute M. |
|                                         |                        |               MIG M. |
|=========================================+========================+======================|
|   0  NVIDIA GeForce RTX 3090        Off |   00000000:18:00.0 Off |                  N/A |
| 30%   49C    P2            106W /  350W |    6106MiB /  24576MiB |      1%      Default |
|                                         |                        |                  N/A |
+-----------------------------------------+------------------------+----------------------+
|   1  NVIDIA GeForce RTX 3090        Off |   00000000:3B:00.0 Off |                  N/A |
| 30%   26C    P8             27W /  350W |    5114MiB /  24576MiB |      0%      Default |
|                                         |                        |                  N/A |
+-----------------------------------------+------------------------+----------------------+
|   2  NVIDIA GeForce RTX 3090        Off |   00000000:86:00.0 Off |                  N/A |
| 30%   39C    P2            130W /  350W |    7560MiB /  24576MiB |      0%      Default |
|                                         |                        |                  N/A |
+-----------------------------------------+------------------------+----------------------+
|   3  NVIDIA GeForce RTX 3090        Off |   00000000:AF:00.0 Off |                  N/A |
| 30%   37C    P2            113W /  350W |    8724MiB /  24576MiB |      0%      Default |
|                                         |                        |                  N/A |
+-----------------------------------------+------------------------+----------------------+

Windows

+-----------------------------------------------------------------------------+
| NVIDIA-SMI 528.92       Driver Version: 528.92       CUDA Version: 12.0     |
|-------------------------------+----------------------+----------------------+
| GPU  Name            TCC/WDDM | Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  NVIDIA GeForce ... WDDM  | 00000000:01:00.0  On |                  N/A |
| N/A   39C    P8     7W /  94W |   2100MiB /  4096MiB |      1%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+

Python environment and version?

Miniconda3 environment with Python 3.10.19

version of brainles_preprocessing ?

brainles_preprocessing==0.6.8

[neuronflow]

Interesting, I recommend trying a different registration backend for your data.
Elastix support is a bit experimental. I recommend trying ANTs or niftyreg.

[neuronflow]

Hmm, I am not sure whether this qualifies as a bug, for me it looks like a resampling artifact. If there is a bug it is likely within the elastix codebase:
https://github.com/SuperElastix/elastix

What's your pov @MarcelRosier ?

[huashengmao]

I've identified the root cause—it's an issue with elastix. My original DICOM data consists partly of integer-type images and partly of float-type images. After converting them to NIfTI format using dcm2niix, they are saved as int16 and float32, respectively.

Since MRI images from different modalities of the same patient share the same data type, elastix produces correct results during the coregistration step. However, during atlas registration—when the float32-type center modalities are registered to the int-type brats_sri24 atlas—elastix generates an erroneous displacement field. This ultimately leads to "snowflake"-like artifacts in both the center and moving modalities.

Initially, I tried using the default ANTs registration tool, but it threw an error and failed to produce any registration result.