Paso 1: Cargue imágenes

Este ejemplo utiliza dos imágenes 3-D de la cabeza del mismo paciente. En los problemas de registro, consideramos que una imagen es la imagen fija y que la otra imagen es la imagen en movimiento. El objetivo del registro es alinear la imagen en movimiento con la imagen fija. En este ejemplo, la imagen fija es una imagen de resonancia magnética ponderada T1. La imagen en movimiento que queremos registrar es una imagen CT. Los datos se almacenan en el formato de archivo utilizado por el proyecto de evaluación de registro de imágenes retrospectivas (RIRE). Se usa para leer los archivos binarios que contienen datos de imagen.multibandread Utilice la función para obtener los metadatos asociados a cada imagen.helperReadHeaderRIRE Puede utilizar el siguiente enlace para encontrar más información sobre el formato de archivo RIRE:Formato de datos RIRE

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

La función es una función auxiliar que se proporciona para ayudar a juzgar la calidad de los resultados del registro en 3-D.helperVolumeRegistration Puede rotar la vista de forma interactiva y ambos ejes permanecerán sincronizados.

helperVolumeRegistration(fixedVolume,movingVolume);

También puede utilizar para ver planos individuales de los volúmenes fijos y en movimiento para obtener una idea de la alineación general de los volúmenes.imshowpair En la imagen superpuesta de, las áreas grises corresponden a áreas que tienen intensidades similares, mientras que las áreas magenta y verde muestran lugares donde una imagen es más brillante que la otra.imshowpair Se utiliza para observar el registro erróneo de los volúmenes de imagen a lo largo de una rodaja axial tomada a través del centro de cada volumen.imshowpair

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

Paso 2: configurar el registro inicial

La función hace que sea fácil elegir el optimizador y la configuración de métricas correctas para usar con.imregconfigimregister Estas dos imágenes provienen de dos modalidades diferentes, MRI y CT, por lo que la opción ' multimodal ' es apropiada.

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

El algoritmo utilizado por convergerá para obtener mejores resultados más rápidamente cuando se especifica información de referencia espacial sobre la resolución y/o la ubicación de las imágenes de entrada.imregister En este caso, la resolución de los datasets CT y MRI se define en los metadatos de la imagen. Utilice estos metadatos para construir objetos de referencia espacial que se pasarán como argumentos de entrada para el registro.imref3d

Rfixed  = imref3d(size(fixedVolume),fixedHeader.PixelSize(2),fixedHeader.PixelSize(1),fixedHeader.SliceThickness); Rmoving = imref3d(size(movingVolume),movingHeader.PixelSize(2),movingHeader.PixelSize(1),movingHeader.SliceThickness);

Las propiedades de los objetos de referencia espacial definen dónde están los volúmenes de imagen asociados en el sistema de coordenadas del mundo y cuál es la extensión de píxel en cada dimensión. La propiedad XWorldLimits de rmoving define la posición del volumen móvil en la cota X. La propiedad PixelExtentInWorld define el tamaño de cada píxel en las unidades del mundo en la dimensión X (a lo largo de las columnas). El volumen móvil se extiende de 0,3269 a 334,97 en el sistema de coordenadas X del mundo y cada píxel tiene una extensión de 0.6536 mm. Las unidades están en milímetros porque la información del encabezado utilizada para construir la referencia espacial estaba en milímetros. La propiedad ImageExtentInWorldX determina la extensión completa de las unidades del mundo del volumen de imágenes en movimiento en unidades del mundo.

Rmoving.XWorldLimits

ans = 1×2

    0.3268  334.9674

Rmoving.PixelExtentInWorldX

ans = 0.6536 

Rmoving.ImageExtentInWorldX

ans = 334.6406 

La desalineación entre los dos volúmenes incluye la traducción, el escalado y la rotación. Utilice una transformación de similitud para registrar las imágenes.

Empiece por usar para obtener un volumen de imagen de salida registrado que puede ver y observar directamente para acceder a la calidad de los resultados del registro.imregister

Especifique una configuración no predeterminada para la propiedad InitialRadius del optimizador para lograr una mejor convergencia en los resultados del registro.

optimizer.InitialRadius = 0.004; movingRegisteredVolume = imregister(movingVolume,Rmoving, fixedVolume,Rfixed, 'rigid', optimizer, metric);

Use de nuevo y repita el proceso de examinar la alineación de una rebanada axial tomada a través del centro de los volúmenes registrados para tener una idea de lo exitoso que es el registro.imshowpair

figure imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3))); title('Axial Slice of Registered Volume')

Desde el sector axial anterior, parece que el registro se realizó correctamente. Utilice de nuevo para ver el volumen registrado para seguir juzgando el éxito del registro.helperVolumeRegistration

helperVolumeRegistration(fixedVolume,movingRegisteredVolume);

Paso 3: obtener transformación geométrica 3-D que alinea el movimiento con fijo.

La función se puede utilizar cuando esté interesado en la estimación de transformación geométrica que se utiliza para formar la imagen de salida registrada. utiliza el mismo algoritmo que y toma los mismos argumentos de entrada que.imregtformimregisterimregtformimregisterimregister Dado que la inspección visual del volumen resultante indica que el registro se realizó correctamente, puede llamar con los mismos argumentos de entrada para obtener la transformación geométrica asociada a este resultado de registro.imregisterimregtform

geomtform = imregtform(movingVolume,Rmoving, fixedVolume,Rfixed, 'rigid', optimizer, metric)

geomtform =    affine3d with properties:      Dimensionality: 3                  T: [4x4 double]  

El resultado de imregtform es un objeto de transformación geométrica. Este objeto incluye una propiedad, T, que define la matriz de transformación de afín 3-D.

geomtform.T

ans = 4×4

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

El método transformPointsForward de la transformación geométrica se puede utilizar para determinar dónde se asigna un punto [u, v, w] en la imagen en movimiento como resultado del registro. Dado que las entradas referenciadas espacialmente se especificaron para imregtform, la transformación geométrica asigna puntos en el sistema de coordenadas del mundo de mover a fijo. El método transformPointsForward se utiliza a continuación para determinar la ubicación transformada del centro de la imagen en movimiento en el sistema de coordenadas del mundo.

centerXWorld = mean(Rmoving.XWorldLimits); centerYWorld = mean(Rmoving.YWorldLimits); centerZWorld = mean(Rmoving.ZWorldLimits); [xWorld,yWorld,zWorld] = transformPointsForward(geomtform,centerXWorld,centerYWorld,centerZWorld);

Puede utilizar el método worldToSubscript de Rfixed para determinar el elemento del volumen fijo que se alinea con el centro del volumen en movimiento.

[r,c,p] = worldToSubscript(Rfixed,xWorld,yWorld,zWorld)

r = 116 

c = 132 

p = 13 

Paso 4: aplicar estimación de transformación geométrica a volumen de imagen en movimiento.

La función se puede utilizar para aplicar la estimación de transformación geométrica de imregtform a un volumen 3-D.imwarp El nombre/valor ' OutputView ' se utiliza para definir un argumento de referencia espacial que determina los límites mundiales y la resolución de la imagen remuestreada de salida. Puede producir los mismos resultados dados mediante el objeto de referencia espacial asociado a la imagen fija.imregister Esto crea un volumen de salida en el que los límites del mundo y la resolución de la imagen fija y en movimiento son los mismos. Una vez que los límites del mundo y la resolución de ambos volúmenes son los mismos, hay correspondencia de píxel a píxel entre cada muestra de los volúmenes móviles y fijos.

movingRegisteredVolume = imwarp(movingVolume,Rmoving,geomtform,'bicubic','OutputView',Rfixed);

Utilice imshowpair de nuevo para ver un corte axial a través del centro del volumen registrado producido por.imwarp

figure  imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3))); title('Axial Slice of Registered Volume')