A MATLAB framework for simulating MRI sequences in the presence of metal implants, with a focus on metal artifact reduction techniques
"Open-Source Simulator of Imaging Near Metal at Arbitrary Magnetic Field Strengths" by Kübra Keskin, Ana R. Sanson, Brian A. Hargreaves, Krishna S. Nayak.
- Clone or download this repository
- Download
data.zipfrom - Create a folder
./dataand place the data inside - Install required MATLAB toolboxes:
- Statistics and Machine Learning Toolbox
- Signal Processing Toolbox
- Image Processing Toolbox
Default simulation parameters are defined in config/getSimulationConfig.m. These parameters can also be overwritten in the main scripts.
The simulator allows configuration of key MRI parameters including:
- B0 field strength
- Readout bandwidth
- RF bandwidth
- TR, TE
- Number of averages
- Resolution
- Slice thickness
- MRI sequence type
- Number of spectral encodings
- Noise parameters
The simulator uses phantom data that includes tissue property maps (T1, T2, proton density) and B0 field maps. The shared phantom data is available in different configurations:
-
Implant Materials:
CoCr: THA with CoCr cup, Polyethylene liner, CoCr head, Titanium stemTiCer: THA with Titanium cup, Polyethylene liner, Ceramic head, Titanium stemPlastic: THA with plastic materials (control case, only 0.5 mm version available)
-
Phantom Resolutions:
02mm: High resolution (0.2 mm isotropic)05mm: Standard resolution (0.5 mm isotropic)
-
Phantom Type:
combined: phantom with realistic human model combined with implantimplant: phantom with only implant model (only 0.5 mm version available)
Phantom files follow the naming pattern: hip_{type}_phantom_{material}_{resolution}.mat
- Examples:
hip_combined_phantom_CoCr_05mm.mat,hip_implant_phantom_TiCer_05mm.mat
Each phantom file contains:
phantom.t1map: T1 relaxation time map (ms)phantom.t2map: T2 relaxation time map (ms)phantom.pdmap: Proton density mapphantom.dfmap: B0 field offset map (Hz/T)phantom.susmap: Susceptibility map (ppm)phantom.mask: Tissue segmentation mask with tissue IDsphantom.resolution: Phantom resolution [x, y, z] (mm)phantom_config.materials: Material specifications
Note: If users want to use their own data, phantom data should follow the above file structure to run the simulation without any error.
Tissue parameters (T1, T2, PD) can be customized from simulation/helpers/prepareParameterMaps.m. The function supports:
- B0-dependent T1 modeling: T1 exponential model based on field strength
- Custom tissue values: Override default parameters for specific tissues by editing the customization section
-
Change the
phantom_fileparameter insidemain_simulation.mbased on the desired phantom data filename. The filename can be changed based on implant material (e.g., "TiCer" or "Plastic" instead of "CoCr") and phantom resolution (e.g., "02mm" instead of "05mm"):phantom_file = dir('./data/hip_combined_phantom*CoCr*05mm*.mat');
-
Customize simulation parameters by editing the default config values or editing overrides in
main_simulation.m:sim_config.system.B0 = 3.0; % Field strength (T) sim_config.imaging.readBWpix = 600; % Readout bandwidth (Hz/px) sim_config.imaging.rfBW = 1.5; % RF bandwidth (kHz)
-
Run simulation:
main_simulation.m
-
Change the
phantom_fileparameter insidemain_parameter_sweep.mbased on the desired phantom data filename. The filename can be changed based on implant material (e.g., "TiCer" or "Plastic" instead of "CoCr") and phantom resolution (e.g., "02mm" instead of "05mm"):phantom_file = dir('./data/hip_combined_phantom*CoCr*05mm*.mat');
-
Customize simulation parameters by editing the default config values or editing overrides in
run_parameter_sweep.m:sim_config.system.B0 = 3.0; % Field strength (T) sim_config.imaging.readBWpix = 600; % Readout bandwidth (Hz/px) sim_config.imaging.rfBW = 1.5; % RF bandwidth (kHz)
-
Define sweep parameters in
main_parameter_sweep.m:sweep_parameter = 'B0'; sweep_values = [0.55, 1.5, 3.0, 7.0];
-
Run the parameter sweep:
main_parameter_sweep.m
The parameter sweep script supports sweeping any of these parameters:
B0: Field strengthreadBWpix: Readout bandwidthrfBW: RF bandwidthNbins: Spectral binsseqname: SequenceTR: Repetition timeTE: Echo time
See main_parameter_sweep.m for examples of different parameter sweep configurations.
The simulator supports the following sequences:
- TSE (Turbo Spin Echo)
- VAT (View Angle Tilting)
- SEMAC (Slice Encoding for Metal Artifact Correction)
- MAVRIC (Multi-Acquisition Variable-Resonance Image Combination)
- MAVRIC-SL (MAVRIC-SeLective)
Note: Simulations can take several minutes to hours depending on the number of slices, spectral bins, chosen sequence, and phantom resolution.
For quick testing, use the 0.5 mm phantom and simulate TSE or VAT for only a few slices. The parameters set in main_simulation.m and main_parameter_sweep.m are specifically chosen for rapid testing.
Each simulation produces sample images visualized via the visualizeSimulationResults function called inside the main scripts.
Each simulation also generates various outputs saved inside a .mat file:
Filename Format: {sequence}_{materials}_{B0}T_Nbins{#bins}_readBW{readout_bw}_rfBW{rf_bw}_res{phantom_resolution}_{timestamp}.mat
Example: SEMAC_CoCr_0.55T_Nbins6_readBW400_rfBW1.0_res0.5_20250716_104931.mat
final_im: Final reconstructed images (struct)im: Raw simulated images (before final reconstruction) (matrix)params: Complete simulation parameters used (struct)kspace: K-space data without noise (matrix)noise_params: Noise characteristics (struct)phantom_config: Phantom configuration details (struct)sim_config: Original simulation configuration (struct)phantom_file: Source phantom file information (struct)
Results are organized as follows:
./results/: Single simulation results./results/parameter_sweep/: Parameter sweep results
The script to analyze signal percentages inside all spectral bins for a specific slice is available in analysis/main_spectral_bins_analysis.m. This analysis is recommended to be used with only implant data (without digital human model) and for simulations performed with high number of spectral bins.
To run with already available simulated results:
-
Place .mat files inside
./results/parameter_sweep -
Inside the script choose implant material, field strength, slice to be analyzed, and specific distance from the boundary of implant to be analyzed, and run.
-
If
save_resultsflag is set to true, results of the analysis will be saved in./results/analysisdirectory.
mri_metal_simulation/
├── config/ # Default simulation configuration
├── simulation/
│ ├── sequences/ # MRI sequence simulation functions
│ ├── reconstruction/ # Image reconstruction functions
│ ├── helpers/ # Helper functions for simulation
│ └── utils/ # General utility functions
├── analysis/ # Data analysis scripts
├── data/ # Implant phantom data
├── results/ # Simulation results save folder
│ ├── parameter_sweep/ # Parameter sweep simulation results save folder
│ └── analysis/ # Data analysis results
├── main_simulation.m # Main simulation script
├── main_parameter_sweep.m # Main parameter sweep script
└── run_parameter_sweep.m # Helper for parameter sweep script
-
Lu W, Pauly KB, Gold GE, Pauly JM, Hargreaves BA. SEMAC: Slice Encoding for Metal Artifact Correction in MRI. Magn Reson Med. 2009;62(1):66-76. doi:10.1002/mrm.21967
-
Koch KM, Lorbiecki JE, Hinks RS, King KF. A multispectral three-dimensional acquisition technique for imaging near metal implants. Magn Reson Med. 2009;61(2):381-390. doi:10.1002/mrm.21856
-
Koch KM, Brau AC, Chen W, et al. Imaging near metal with a MAVRIC-SEMAC hybrid. Magn Reson Med. 2011;65(1):71-82. doi:10.1002/mrm.22523
-
Koch KM, Hargreaves BA, Pauly KB, Chen W, Gold GE, King KF. Magnetic resonance imaging near metal implants. J Magn Reson Imaging. 2010;32(4):773-787. doi:10.1002/jmri.22313
-
Hargreaves BA, Worters PW, Pauly KB, Pauly JM, Koch KM, Gold GE. Metal-induced artifacts in MRI. AJR Am J Roentgenol. 2011;197(3):547-555. doi:10.2214/AJR.11.7364
-
Shi X, Yoon D, Koch KM, Hargreaves BA. Metallic implant geometry and susceptibility estimation using multispectral B0 field maps. Magn Reson Med. 2017;77(6):2402-2413. doi:10.1002/mrm.26313
This code is tested in MATLAB 2021a and 2025a.
For questions or issues please contact (keskin AT usc DOT edu).
Copyright (c) 2025 Kübra Keskin, Magnetic Resonance Engineering Laboratory. MIT license.