Register Multimodal Medical Image Volumes with Spatial Referencing
This example shows how to align two medical image volumes using moment-of-mass-based registration.
Multimodal image registration aligns images acquired using different imaging modalities, such as MRI and CT. Even when acquired for the same patient and region of interest, multimodal images can be misaligned due to differences in patient positioning (translation or rotation) and voxel size (scaling). Different imaging modalities often have different voxel sizes due to scanner variability and concerns about scan time and radiation dose.
imregmoment function enables fast registration of 3-D image volumes, including multimodal images. To register images with scaling differences using
imregmoment, you must define the spatial referencing for the images. Spatial referencing maps the row, column, and slice indices of the image array to the patient coordinate system. You can use the
medicalVolume object to automatically import image volumes and spatial referencing metadata from an image file.
In this example, you use
imregmoment to register multimodal MRI and CT image of the head.
The data used in this example is a modified version of the 3-D CT and MRI data sets from The Retrospective Image Registration Evaluation (RIRE) Dataset, provided by Dr. Michael Fitzpatrick. For more information, see the RIRE Project homepage. The modified data set contains one CT scan and one MRI scan stored in the NRRD file format. The size of the entire data set is approximately 35 MB. Download the data set from the MathWorks® website, then unzip the folder.
zipFile = matlab.internal.examples.downloadSupportFile("medical","MedicalRegistrationNRRDdata.zip"); filepath = fileparts(zipFile); unzip(zipFile,filepath)
In image registration, consider one image to be the fixed image and the other image to be the moving image. The goal of registration is to align the moving image with the fixed image. In this example, the fixed image is a T1 weighted MRI image. The moving image is a CT image from the same patient. The images are stored in the NRRD file format.
medicalVolume to read the MRI image. By looking at the
VoxelSpacing properties, you can determine that the MRI volume is a 256-by-256-by-26 voxel array, where each voxel is 1.25-by-1.25-by-4.0 mm.
filenameMRI = fullfile(filepath,"supportfilesNRRD","Patient007MRT1.nrrd"); fixedMRIVolume = medicalVolume(filenameMRI)
fixedMRIVolume = medicalVolume with properties: Voxels: [256×256×26 single] VolumeGeometry: [1×1 medicalref3d] SpatialUnits: "unknown" Orientation: "unknown" VoxelSpacing: [1.2500 1.2500 4] NormalVector: [0 0 1] NumCoronalSlices:  NumSagittalSlices:  NumTransverseSlices:  PlaneMapping: ["unknown" "unknown" "unknown"] Modality: "unknown" WindowCenters:  WindowWidths: 
Import the CT image. The size of each voxel in the CT image is 0.65-by-0.65-by-4 mm, which is smaller than in the MRI image. Therefore, the CT volume contains more voxels than the MRI scan for the same region of interest.
filenameCT = fullfile(filepath,"supportfilesNRRD","Patient007CT.nrrd"); movingCTVolume = medicalVolume(filenameCT)
movingCTVolume = medicalVolume with properties: Voxels: [512×512×28 single] VolumeGeometry: [1×1 medicalref3d] SpatialUnits: "unknown" Orientation: "unknown" VoxelSpacing: [0.6536 0.6536 4] NormalVector: [0 0 1] NumCoronalSlices:  NumSagittalSlices:  NumTransverseSlices:  PlaneMapping: ["unknown" "unknown" "unknown"] Modality: "unknown" WindowCenters:  WindowWidths: 
Extract the image voxel data from the
Voxels property of the
fixedMRIVoxels = fixedMRIVolume.Voxels; movingCTVoxels = movingCTVolume.Voxels;
Display Unregistered Images
helperVolumeRegistration helper function to judge image alignment. The output aces remain in syn if you interactively rotate the view.
The plot is in intrinsic coordinates, in units of voxels. The axes limits are different between the two images because of the differences in image array size.
You can also use
imshowpair to check alignment in single slices from the fixed and moving volumes. Use
imshowpair to view the middle transverse slice of each volume.
First, calculate the location of the middle slice of each volume. The
VolumeGeometry property of the medical volume object contains a
medicalref3d object, which specifies spatial details about the image. Use the
VolumeSize property of each
medicalref3d object to calculate the location of the middle slice of the corresponding volume, in voxels.
fixedVolumeSize = fixedMRIVolume.VolumeGeometry.VolumeSize; movingVolumeSize = movingCTVolume.VolumeGeometry.VolumeSize; centerFixed = fixedVolumeSize/2; centerMoving = movingVolumeSize/2;
Next, display the image slices using
imshowpair. Plot the slice in the third spatial dimension, which corresponds to the transverse anatomical plane. The MRI slice is magenta, and the CT slice is green. The images are misaligned due to translations, rotations, and scaling.
figure imshowpair(movingCTVoxels(:,:,centerMoving(3)),fixedMRIVoxels(:,:,centerFixed(3))) title("Unregistered Transverse Slice")
To obtain accurate scaling results using
imregmoment, you must specify spatial referencing information for each volume. Use the spatial referencing information stored in the
medicalVolume objects to define corresponding
imref3d spatial referencing objects. Specify the image array dimensions and voxel sizes using the
VoxelSpacing properties, respectively, of the
Rfixed = imref3d(fixedVolumeSize,fixedMRIVolume.VoxelSpacing(2),fixedMRIVolume.VoxelSpacing(1),fixedMRIVolume.VoxelSpacing(3)); Rmoving = imref3d(movingVolumeSize,movingCTVolume.VoxelSpacing(2),movingCTVolume.VoxelSpacing(1),movingCTVolume.VoxelSpacing(3));
imregmoment to register the moving volume to the fixed volume. Specify the
MedianThresholdBitmap name-value argument as
true, which is appropriate for multimodal images .
[geomtform,movingCTRegisteredVoxels] = imregmoment(movingCTVoxels,Rmoving,fixedMRIVoxels,Rfixed,MedianThresholdBitmap=true);
geomtform output is an
affinetform3d geometric transformation object. The
T property of
geomtform contains the 3-D affine transformation matrix that maps the moving CT volume to the fixed MRI volume.
ans = 4×4 0.5229 -0.0021 0.0007 0 0.0020 0.5229 0.0033 0 -0.0014 -0.0063 1.0000 0 -4.8094 -16.0063 -1.3481 1.0000
movingRegisteredVoxels output contains the registered CT image. The
imregmoment function resamples the registered image coordinates using the voxel grid of the fixed image. Therefore, the registered volume and fixed volume are the same size and have voxel-to-voxel correspondence.
whos movingCTRegisteredVoxels fixedMRIVoxels
Name Size Bytes Class Attributes fixedMRIVoxels 256x256x26 6815744 single movingCTRegisteredVoxels 256x256x26 6815744 single
Display Registered Volumes
To check the registration results, use
imshowpair to view the middle transverse slices of the fixed and registered volumes. The images are aligned and scaled to the same size.
figure imshowpair(movingCTRegisteredVoxels(:,:,centerFixed(3)),fixedMRIVoxels(:,:,centerFixed(3))) title("Registered Transverse Slice")
helperVolumeRegistration to view the results of the 3-D registration. The axes limits are the same, because of scaling the CT volume during registration.
medicalVolume Object for Registered Image
Create a new
medicalVolume object that contains the registered voxel data and its spatial referencing information. You can create a
medicalVolume object by specifying a voxel array and a
medicalref3d object containing spatial details. The registered CT volume has the same spatial referencing as the original MRI volume, so use the
medicalref3d object stored in the
VolumeGeometry property of
R = fixedMRIVolume.VolumeGeometry; movingRegisteredVolume = medicalVolume(movingCTRegisteredVoxels,R)
movingRegisteredVolume = medicalVolume with properties: Voxels: [256×256×26 single] VolumeGeometry: [1×1 medicalref3d] SpatialUnits: "unknown" Orientation: "unknown" VoxelSpacing: [1.2500 1.2500 4] NormalVector: [0 0 1] NumCoronalSlices:  NumSagittalSlices:  NumTransverseSlices:  PlaneMapping: ["unknown" "unknown" "unknown"] Modality: "unknown" WindowCenters:  WindowWidths: 
Verify that the new
medicalVolume object contains CT voxel data that is aligned with the MRI image by using
figure imshowpair(movingRegisteredVolume.Voxels(:,:,centerFixed(3)), fixedMRIVoxels(:,:,centerFixed(3))) title("Axial Slice of Registered Volume")