Шаг 1: Загрузка изображений

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

The 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: Настройка начальной регистрации

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

[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. Свойство PixelExtendInWorld определяет размер каждого пикселя в мировых единицах в размерность X (вдоль столбцов). Подвижный объем простирается от 0,3269 до 334,97 в мировой системе координат X, и каждый пиксель имеет протяженность 0,6536 мм. Модули указаны в миллиметрах, потому что информация о заголовке, используемая для построения пространственной ссылки, была в миллиметрах. Свойство ImageExtendInWorldX определяет полный объем в мировых единицах измерения объема движущегося изображения в мировых единицах измерения.

Rmoving.XWorldLimits

ans = 1×2

    0.3268  334.9674

Rmoving.PixelExtentInWorldX

ans = 0.6536

Rmoving.ImageExtentInWorldX

ans = 334.6406

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

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

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

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: Получите 3-D геометрическое преобразование, которое выравнивает движение с фиксированным.

The 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

The transformPointsForward функция может использоваться, чтобы определить, где точка [u, v, w] в движущемся изображении отображается в результате регистрации. Поскольку пространственно-ссылочные входы были заданы как imregtformгеометрическое преобразование сопоставляет точки в мировой системе координат с перемещением на фиксированные. The 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. Примените оценку геометрического преобразования к объему движущегося изображения.

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