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.

The 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 medicalVolume and imregmoment to register multimodal MRI and CT image of the head.

Load Images

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.

Use medicalVolume to read the MRI image. By looking at the Voxels and 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 medicalVolume objects.

fixedMRIVoxels = fixedMRIVolume.Voxels;
movingCTVoxels = movingCTVolume.Voxels;

Display Unregistered Images

Use the 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.

helperVolumeRegistration(fixedMRIVoxels,movingCTVoxels);

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")

Register Images

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 VolumeSize and VoxelSpacing properties, respectively, of the medicalVolume objects.

Rfixed  = imref3d(fixedVolumeSize,fixedMRIVolume.VoxelSpacing(2),fixedMRIVolume.VoxelSpacing(1),fixedMRIVolume.VoxelSpacing(3));
Rmoving = imref3d(movingVolumeSize,movingCTVolume.VoxelSpacing(2),movingCTVolume.VoxelSpacing(1),movingCTVolume.VoxelSpacing(3));

Use 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);

The 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.

geomtform.T

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

The 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")

Use helperVolumeRegistration to view the results of the 3-D registration. The axes limits are the same, because of scaling the CT volume during registration.

helperVolumeRegistration(fixedMRIVoxels,movingCTRegisteredVoxels);

Create 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 fixedMRIVolume.

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 imshowpair.

figure 
imshowpair(movingRegisteredVolume.Voxels(:,:,centerFixed(3)), fixedMRIVoxels(:,:,centerFixed(3)))
title("Axial Slice of Registered Volume")