Шаг 1: Загрузить изображения

В этом примере используются два 3-D изображения головы одного и того же пациента. При проблемах регистрации мы рассматриваем одно изображение как фиксированное изображение, а другое изображение как движущееся изображение. Целью регистрации является выравнивание движущегося изображения с фиксированным изображением. В этом примере фиксированное изображение представляет собой T1 взвешенное МРТ-изображение. Движущееся изображение, которое мы хотим зарегистрировать, является CT-изображением. Данные хранятся в формате файла, используемом в проекте ретроспективной оценки регистрации изображений (RIRE). Использовать multibandread для чтения двоичных файлов, содержащих данные изображения. Используйте helperReadHeaderRIRE получение метаданных, связанных с каждым изображением. Для получения дополнительной информации о формате файла 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' );

helperVolumeRegistration функция - это вспомогательная функция, которая предназначена для оценки качества результатов регистрации 3-D. Вид можно повернуть в интерактивном режиме, и обе оси останутся синхронизированными.

helperVolumeRegistration(fixedVolume,movingVolume);

Также можно использовать imshowpair для просмотра отдельных плоскостей из фиксированных и движущихся объемов, чтобы получить представление об общем выравнивании объемов. На перекрывающемся изображении из imshowpairсерые области соответствуют областям с одинаковой интенсивностью, в то время как пурпурные и зеленые области показывают места, где одно изображение ярче другого. Использовать imshowpair для наблюдения за неправильной регистрацией объемов изображения вдоль осевого среза, проходящего через центр каждого объема.

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

Шаг 2: Настройка начальной регистрации

imregconfig функция упрощает выбор правильной конфигурации оптимизатора и метрики для использования с imregister. Эти два изображения из двух различных методов, МРТ и КТ, поэтому «мультимодальный» вариант является подходящим.

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

Алгоритм, используемый imregister будет сходиться к лучшим результатам более быстро, когда пространственно-ссылочная информация о разрешении и/или местоположении входных изображений определена. В этом случае разрешение наборов данных КТ и МРТ определяется в метаданных изображения. Использовать эти метаданные для построения 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);

Свойства пространственных ссылочных объектов определяют, где находятся связанные объемы изображения в мировой системе координат и какова протяженность пикселей в каждом измерении. Свойство XWorldLimits команды Rmoving определяет положение движущегося объема в измерении X. Свойство PixelInWorld определяет размер каждого пикселя в единицах мира в измерении X (вдоль столбцов). Движущийся объем простирается от 0,3269 до 334,97 в мировой системе координат X, и каждый пиксель имеет протяженность 0,6536 мм. Единицы измерения выражены в миллиметрах, поскольку информация заголовка, используемая для построения пространственной ссылки, была в миллиметрах. Собственность ImageExtentInWorldX определяет полный объем в мировых единицах объема движущегося изображения в мировых единицах.

Rmoving.XWorldLimits

ans = 1×2

    0.3268  334.9674

Rmoving.PixelExtentInWorldX

ans = 0.6536

Rmoving.ImageExtentInWorldX

ans = 334.6406

Несоответствие между двумя томами включает перемещение, масштабирование и поворот. Используйте преобразование подобия для регистрации изображений.

Начать с использования imregister для получения зарегистрированного объема выходного изображения, который можно просматривать и наблюдать непосредственно, чтобы получить доступ к качеству результатов регистрации.

Для достижения лучшей сходимости в результатах регистрации задайте значение свойства InitityRadius оптимизатора, отличное от значения по умолчанию.

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

Использовать imshowpair снова и повторить процесс исследования выравнивания осевого среза, взятого через центр зарегистрированных объемов, чтобы получить представление о том, насколько успешна регистрация.

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

Судя по приведенному выше осевому срезу, регистрация прошла успешно. Использовать helperVolumeRegistration снова просмотреть зарегистрированный том, чтобы продолжить оценку успешности регистрации.

helperVolumeRegistration(fixedVolume,movingRegisteredVolume);

Шаг 3: получите 3D геометрическое преобразование, которое выравнивает перемещение с фиксированным.

imregtform функция может использоваться, когда вы заинтересованы в оценке геометрического преобразования, которая используется imregister для формирования зарегистрированного выходного изображения. imregtform использует тот же алгоритм, что и imregister и принимает те же входные аргументы, что и imregister. С момента визуального осмотра полученного объема от imregister указал, что регистрация прошла успешно, вы можете позвонить imregtform с теми же входными аргументами, чтобы получить геометрическое преобразование, связанное с этим результатом регистрации.

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

geomtform = 
  affine3d with properties:

                 T: [4x4 double]
    Dimensionality: 3

Результат imregtform является объектом геометрического преобразования. Этот объект включает свойство T, которое определяет 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

transformPointsForward функция может использоваться для определения, где точка [u, v, w] в картах движущегося изображения в результате регистрации. Поскольку пространственно привязанные входные данные были указаны дляimregtform, геометрическое преобразование отображает точки в мировой системе координат от перемещения к фиксированному. transformPointsForward используется ниже для определения преобразованного местоположения центра движущегося изображения в мировой системе координат.

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

Вы можете использовать worldToSubscript для определения элемента фиксированного объема, который выравнивается по центру движущегося объема.

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

r = 116

c = 132

p = 13

Шаг 4. Применить оценку геометрического преобразования к объему движущегося изображения.

imwarp функция может использоваться для применения оценки геометрического преобразования из imregtform к 3-D объему. Аргумент имя-значение «OutputView» используется для определения аргумента пространственной ссылки, который определяет мировые пределы и разрешение выходного повторно дискретизированного изображения. Можно получить те же результаты, что и imregister путем использования объекта пространственной ссылки, связанного с фиксированным изображением. Это создает выходной объем, в котором мировой предел и разрешение фиксированного и движущегося изображения одинаковы. Как только мировые пределы и разрешение обоих объемов одинаковы, между каждым образцом движущегося и фиксированного объемов имеется соответствие пиксел-пиксель.

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

Использовать imshowpair снова для просмотра осевого среза через центр зарегистрированного объема, производимого imwarp.

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