Register Multimodal 3-D Medical Images

Load Images

This example uses a CT image and a T1 weighted MRI image collected from the same patient. The 3-D CT and MRI data sets are from The Retrospective Image Registration Evaluation (RIRE) Dataset, provided by Dr. Michael Fitzpatrick. For more information, see the RIRE Project homepage.

In this example, the MRI scan is the fixed image and the CT image is the moving image. The images are stored in the file formats used by the RIRE Project. Use the helperReadHeaderRIRE helper function to obtain the metadata associated with each image. The helper function is attached to this example as a supporting file.

fixedHeader = helperReadHeaderRIRE("rirePatient007MRT1.header");
movingHeader = helperReadHeaderRIRE("rirePatient007CT.header");

Use multibandread to read the binary files that contain the image data.

fixedVolume = multibandread("rirePatient007MRT1.bin", ...
    [fixedHeader.Rows fixedHeader.Columns fixedHeader.Slices], ...
    "int16=>single",0,"bsq","ieee-be");
                        
movingVolume = multibandread("rirePatient007CT.bin", ...
    [movingHeader.Rows movingHeader.Columns movingHeader.Slices], ...
    "int16=>single",0,"bsq","ieee-be");

Display Unregistered Images

Judge the alignment of the unregistered volumes by displaying the middle transverse slice of each volume using imshowpair. The MRI slice is magenta, and the CT slice is green. The scaling differences indicate that the images have different pixel spacing.

centerFixed = size(fixedVolume,3)/2;
centerMoving = size(movingVolume,3)/2;
figure
imshowpair(movingVolume(:,:,centerMoving),fixedVolume(:,:,centerFixed))
title("Unregistered Transverse Slice")

You can also display the volumes as 3-D objects. Create a viewer3d object, in which you can display multiple volumes.

viewerUnregistered = viewer3d(BackgroundColor="black",BackgroundGradient="off");

Display the medicalVolume objects as 3-D isosurfaces by using the volshow function. You can rotate the volumes by clicking and dragging in the viewer window. When you plot the 3-D isosurfaces in intrinsic coordinates, they appear vertically distorted because the slice thickness is different than the in-plane pixel spacing in world coordinates.

volshow(fixedVolume,Parent=viewerUnregistered,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1);
volshow(movingVolume,Parent=viewerUnregistered,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1);

Define Spatial Referencing

You can improve the display and registration results for images with different pixel spacing by using spatial referencing. For this data set, the header files define the pixel spacing of the CT and MRI images. Use this metadata, plus the size of each volume in voxels, to create imref3d spatial referencing objects. The imref3d object properties define the position of an image volume in its world coordinate system and the pixel spacing in each dimension. For example, the PixelExtentInWorldX properties of the fixed and moving imref3d objects indicate pixel spacing along the X-axes of the fixed and moving volumes of 1.25 mm and 0.6536 mm, respectively. Units are in millimeters because the header information used to construct the spatial referencing is in millimeters.

Rfixed = imref3d(size(fixedVolume),fixedHeader.PixelSize(2), ...
    fixedHeader.PixelSize(1),fixedHeader.SliceThickness)

Rfixed = 
  imref3d with properties:

           XWorldLimits: [0.6250 320.6250]
           YWorldLimits: [0.6250 320.6250]
           ZWorldLimits: [2 106]
              ImageSize: [256 256 26]
    PixelExtentInWorldX: 1.2500
    PixelExtentInWorldY: 1.2500
    PixelExtentInWorldZ: 4
    ImageExtentInWorldX: 320
    ImageExtentInWorldY: 320
    ImageExtentInWorldZ: 104
       XIntrinsicLimits: [0.5000 256.5000]
       YIntrinsicLimits: [0.5000 256.5000]
       ZIntrinsicLimits: [0.5000 26.5000]

Rmoving = imref3d(size(movingVolume),movingHeader.PixelSize(2), ...
    movingHeader.PixelSize(1),movingHeader.SliceThickness)

Rmoving = 
  imref3d with properties:

           XWorldLimits: [0.3268 334.9674]
           YWorldLimits: [0.3268 334.9674]
           ZWorldLimits: [2 114]
              ImageSize: [512 512 28]
    PixelExtentInWorldX: 0.6536
    PixelExtentInWorldY: 0.6536
    PixelExtentInWorldZ: 4
    ImageExtentInWorldX: 334.6406
    ImageExtentInWorldY: 334.6406
    ImageExtentInWorldZ: 112
       XIntrinsicLimits: [0.5000 512.5000]
       YIntrinsicLimits: [0.5000 512.5000]
       ZIntrinsicLimits: [0.5000 28.5000]

Approach 1: Register Images Using imregister

The imregister function performs registration and returns the registered moving image volume. Use this approach when you want the registered image data and do not need the geometric transformation used to perform the registration.

The imregister function sets the value of fill pixels added to the registered volume to 0. To improve the display of registration results, scale the CT intensities to the range [0, 1], so that the fill value is equal to the minimum of the image data range.

rescaledMovingVolume = rescale(movingVolume);

The imregister function requires you to specify an optimizer and metric configuration for the registration calculation. Define the configuration by using the imregconfig function. Specify the modality as "multimodal" because the CT and MRI images are from different modalities.

[optimizer,metric] = imregconfig("multimodal");

Change the value of the InitialRadius property of the optimizer to achieve better convergence during registration.

optimizer.InitialRadius = 0.004;

Align the images using imregister. Specify the spatial referencing information so the function converges to better results more quickly. The misalignment between the two spatially referenced volumes includes translation and rotation, so use a rigid transformation to register the images.

movingRegisteredVolume = imregister(rescaledMovingVolume,Rmoving,fixedVolume,Rfixed, ...
    "rigid",optimizer,metric);

To assess the registration results, display the middle transverse slices of the fixed volume and registered moving volume. The images appear more aligned. The imregister function also adjusts the spatial referencing of the moving image to match the fixed image, so the images have the same number of voxels in each dimension without scaling differences.

figure
imshowpair(movingRegisteredVolume(:,:,centerMoving),fixedVolume(:,:,centerFixed))
title("Registered Transverse Slice - imregister")

Display the registered volumes as 3-D isosurfaces using volshow. You can zoom and rotate the display to assess the registration results.

viewerRegistered1 = viewer3d(BackgroundColor="black",BackgroundGradient="off");
volFixed1 = volshow(fixedVolume,Parent=viewerRegistered1,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1);
volRegistered1 = volshow(movingRegisteredVolume,Parent=viewerRegistered1,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1);

Optionally, you can view the registered volumes in world coordinates with correct voxel spacing by setting the Transformation property of the Volume objects created by volshow. Specify the Transformation value to scale the volumes based on the voxel spacing from the fixed volume header file. Use the spacing from the fixed volume file for both volumes because the registered moving volume has been resampled in the voxel grid of the fixed volume.

sx = fixedHeader.PixelSize(1);
sy= fixedHeader.PixelSize(2);
sz = fixedHeader.SliceThickness;

A = [sx 0 0 0; 0 sy 0 0; 0 0 sz 0; 0 0 0 1];
tformVol = affinetform3d(A);

volFixed1.Transformation = tformVol;
volRegistered1.Transformation = tformVol;

Approach 2: Estimate and Apply 3-D Geometric Transformation

The imregister function registers images, but does not return information about the geometric transformation applied to the moving image. When you are interested in the estimated geometric transformation, you can use the imregtform function to get a geometric transformation object that stores information about the transformation. The imregtform function has the same input arguments and uses the same algorithm as imregister.

tform = imregtform(rescaledMovingVolume,Rmoving,fixedVolume,Rfixed, ...
    "rigid",optimizer,metric)

tform = 
  rigidtform3d with properties:

    Dimensionality: 3
       Translation: [-15.8648 -17.5692 29.1830]
                 R: [3×3 double]
                 A: [4×4 double]

The A property of tform specifies the 3-D affine transformation matrix used to align the moving image to the fixed image. Because the inputs to imregtform are spatially referenced, the geometric transformation maps points in the world coordinate system from moving to fixed.

tform.A

ans = 4×4

    0.9704    0.0228    0.2404  -15.8648
   -0.0143    0.9992   -0.0369  -17.5692
   -0.2410    0.0324    0.9700   29.1830
         0         0         0    1.0000

You can apply the geometric transformation from imregtform to the moving image volume by using the imwarp function. Specify the moving volume, spatial referencing for the moving volume, and transformation from imregtform. The OutputView name-value argument specifies the spatial referencing for the transformed output image volume. To produce the same results as imregister, specify the OutputView as the imref3d object for the fixed image. This creates a registered image volume with the same spatial referencing as the fixed volume.

movingRegisteredVolume = imwarp(rescaledMovingVolume,Rmoving,tform, ...
    "bicubic",OutputView=Rfixed);

To assess the registration results, display the middle transverse slices of the fixed volume and transformed moving volume. The results are the same as the results from imregister.

figure 
imshowpair(movingRegisteredVolume(:,:,centerFixed), ...
    fixedVolume(:,:,centerFixed))
title("Registered Transverse Slice - imregtform")

Display the registered volumes as 3-D isosurfaces using volshow. View the volume in world coordinates by setting the Transformation property as a name-value argument.

viewerRegistered2 = viewer3d(BackgroundColor="black",BackgroundGradient="off");
volshow(fixedVolume,Parent=viewerRegistered2,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1,Transformation=tformVol);
volshow(movingRegisteredVolume,Parent=viewerRegistered2,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1,Transformation=tformVol);