register multimodal 3-凯发k8网页登录

load images

this example uses a ct image and a t1 weighted mr image collected from the same patient at different time. the 3-d ct and mri data sets were provided by as part of .

this example specifies the mri image as the fixed image and the ct image as the moving image. the data is stored in the file format used by the retrospective image registration evaluation (rire) project. use multibandread to read the binary files that contain image data. use the helperreadheaderrire function to obtain the metadata associated with each image.

fixedheader  = helperreadheaderrire("rirepatient007mrt1.header");
movingheader = helperreadheaderrire("rirepatient007ct.header");
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" );

the helpervolumeregistration function is a helper function that is provided to help judge the quality of 3-d registration results. you can interactively rotate the view and both axes will remain in sync.

figure
helpervolumeregistration(fixedvolume,movingvolume);

you can also use imshowpair to look at single planes from the fixed and moving volumes to get a sense of the overall alignment of the volumes. in the overlapping image from imshowpair, gray areas correspond to areas that have similar intensities, while magenta and green areas show places where one image is brighter than the other. use imshowpair to observe the misregistration of the image volumes along an axial slice taken through the center of each volume. it is clear that the images have different spatial referencing information, such as different world limits and pixel extents.

centerfixed = size(fixedvolume,3)/2;
centermoving = size(movingvolume,3)/2;
figure
imshowpair(movingvolume(:,:,centermoving),fixedvolume(:,:,centerfixed))
title("unregistered axial slice")

you can improve the display and registration results by incorporating spatial referencing information. for this data, the resolution of the ct and mri data sets is defined in the image metadata. use this metadata to create imref3d spatial referencing objects.

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]

the properties of the spatial referencing objects define where the associated image volumes are in the world coordinate system and what the pixel extent in each dimension is. the xworldlimits property of rmoving defines the position of the moving volume in the x dimension. the pixelextentinworld property defines the size of each pixel in world units in the x dimension (along columns). the moving volume extends from 0.3269 to 334.97 mm in the world x coordinate system and each pixel has an extent of 0.6536 mm. units are in millimeters because the header information used to construct the spatial referencing was in millimeters. the imageextentinworldx property determines the full extent of the moving image volume in world units.

approach 1: register images using imregister

the imregister function enables you to obtain a registered output image volume that you can view and observe directly to access the quality of registration results.

pick the correct optimizer and metric configuration to use with imregister by using the imregconfig function. these two images are from two different modalities, mri and ct, so the "multimodal" option is appropriate. change the value of the initialradius property of the optimizer to achieve better convergence in registration results.

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

the misalignment between the two volumes includes translation and rotation so use a rigid transformation to register the images. specify the spatial referencing information so that the algorithm used by imregister will converge to better results more quickly.

movingregisteredvolume = imregister(movingvolume,rmoving,fixedvolume,rfixed, ...
    "rigid",optimizer,metric);

display an axial slice taken through the center of the registered volumes to get a sense of how successful the registration is. in addition to aligning the images, the registration process adjusts the spatial referencing information of the moving image so that it is consistent with the fixed image. the images are now the same size and are successfully aligned.

figure
imshowpair(movingregisteredvolume(:,:,centerfixed),fixedvolume(:,:,centerfixed));
title("axial slice of registered volume")

use helpervolumeregistration again to view the registered volume to continue judging the success of registration.

helpervolumeregistration(fixedvolume,movingregisteredvolume);

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. imregtform uses the same algorithm as imregister and takes the same input arguments as imregister.

tform = imregtform(movingvolume,rmoving,fixedvolume,rfixed, ...
    "rigid",optimizer,metric)

tform = 
  rigidtform3d with properties:
    dimensionality: 3
       translation: [-15.8648 -17.5692 29.1830]
                 r: [3x3 double]
                 a: [4x4 double]

the property a defines the 3-d affine transformation matrix used to align the moving image to the fixed image.

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

the transformpointsforward function can be used to determine where a point [u,v,w] in the moving image maps as a result of the registration. because spatially referenced inputs were specified to imregtform, the geometric transformation maps points in the world coordinate system from moving to fixed. the transformpointsforward function is used below to determine the transformed location of the center of the moving image in the world coordinate system.

centerxworld = mean(rmoving.xworldlimits);
centeryworld = mean(rmoving.yworldlimits);
centerzworld = mean(rmoving.zworldlimits);
[xworld,yworld,zworld] = transformpointsforward(tform, ...
    centerxworld,centeryworld,centerzworld);

you can use the worldtosubscript function to determine the element of the fixed volume that aligns with the center of the moving volume.

[r,c,p] = worldtosubscript(rfixed,xworld,yworld,zworld)

r = 116

c = 132

p = 13

apply the geometric transformation estimate from imregtform to a 3-d volume using the imwarp function. the "outputview" name-value argument is used to define a spatial referencing argument that determines the world limits and resolution of the output resampled image. you can produce the same results given by imregister by using the spatial referencing object associated with the fixed image. this creates an output volume in which the world limits and resolution of the fixed and moving image are the same. once the world limits and resolution of both volumes are the same, there is pixel to pixel correspondence between each sample of the moving and fixed volumes.

movingregisteredvolume = imwarp(movingvolume,rmoving,tform, ...
    "bicubic",outputview=rfixed);

view an axial slice through the center of the registered volume produced by imwarp.

figure 
imshowpair(movingregisteredvolume(:,:,centerfixed), ...
    fixedvolume(:,:,centerfixed));
title("axial slice of registered volume")