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

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

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. Эти два изображения от двух различной модальности, MRI и CT, таким образом, 'многомодальная' опция является соответствующей.

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

Алгоритм используется imregister будет сходиться к лучшим результатам более быстро, когда пространственная информация ссылки о разрешении и/или местоположении входного формирования изображений будет указана. В этом случае разрешение наборов данных CT и MRI задано в метаданных изображения. Используйте эти метаданные, чтобы создать 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 размерностях. Свойство PixelExtentInWorld задает размер каждого пикселя в мировых единицах измерения в 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 получить зарегистрированный объем выходного изображения, который можно просмотреть и наблюдать непосредственно, чтобы получить доступ к качеству регистрационных результатов.

Задайте установку не по умолчанию для свойства 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 геометрическое преобразование, которое выравнивает перемещение с фиксированным.

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

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

geomtform = 
  affine3d with properties:

    Dimensionality: 3
                 T: [4x4 double]

Результатом 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 метод Rfixed, чтобы определить элемент фиксированного объема, который выравнивается с центром движущегося объема.

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