Восстановление изображения из данных о проекции

В этом примере показано, как использовать radon, iradon, fanbeam, и ifanbeam сформировать проекции из демонстрационного изображения и затем восстановить изображение от проекций. В то время как radon и iradon используйте геометрию параллельного луча в проекциях, fanbeam и ifanbeam используйте геометрию луча вентилятора. Чтобы сравнить параллельный луч и конфигурации луча вентилятора, примеры ниже создают синтетические проекции для каждой геометрии и затем используют те синтетические проекции, чтобы восстановить оригинальное изображение.

Реальное приложение, которое требует реконструкции изображений, является томографией поглощения рентгена, где проекции формируются путем измерения затухания излучения, которое проходит через физический экземпляр под различными углами. Оригинальное изображение может считаться сечением через экземпляр, в котором значения интенсивности представляют плотность экземпляра. Проекции собраны специальными медицинскими устройствами обработки изображений, и затем внутреннее изображение экземпляра восстановлено с помощью iradon или ifanbeam.

Функциональный iradon восстанавливает изображение от проекций параллельного луча. В геометрии параллельного луча каждая проекция формируется путем объединения набора линейных интегралов через изображение под определенным углом. Функциональный ifanbeam восстанавливает изображение от проекций луча вентилятора, которые имеют один эмиттер и несколько датчиков.

См. Руководство пользователя Image Processing Toolbox™ для схем, которые иллюстрируют оба конфигураций.

Создайте главный фантом

Тестовое изображение является фантомом главы Shepp-Логана, который может быть сгенерирован с помощью функционального phantom. Фантомное изображение иллюстрирует много качеств, которые найдены в реальной томографической обработке изображений человеческих голов. Яркий эллиптический интерпретатор вдоль внешнего вида походит на череп, и много замещающих знаков внутри походят на мозговые функции или опухоли.

P = phantom(256);
imshow(P)

Параллельный луч - вычисляет синтетические проекции

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

theta1 = 0:10:170; 
[R1,~] = radon(P,theta1); 
num_angles_R1 = size(R1,2)
num_angles_R1 = 18
theta2 = 0:5:175;  
[R2,~] = radon(P,theta2);
num_angles_R2 = size(R2,2)
num_angles_R2 = 36
theta3 = 0:2:178;  
[R3,xp] = radon(P,theta3); 
num_angles_R3 = size(R3,2)
num_angles_R3 = 90

Обратите внимание на то, что для каждого угла, проекция вычисляется в N точек вдоль xp-оси, где N является константой, которая зависит от диагонального расстояния изображения, таким образом, что каждый пиксель будет спроектирован для всех возможных углов проекции.

N_R1 = size(R1,1)
N_R1 = 367
N_R2 = size(R2,1)
N_R2 = 367
N_R3 = size(R3,1)
N_R3 = 367

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

P_128 = phantom(128);
[R_128,xp_128] = radon(P_128,theta1);
N_128 = size(R_128,1)
N_128 = 185

Отобразите данные о проекции R3. Некоторые функции исходного фантомного изображения отображаются в изображении R3. Первый столбец R3 соответствует проекции в 0 градусах, которая объединяется в вертикальном направлении. Ближайший к центру столбец соответствует проекции в 90 градусах, которая объединяется в горизонтальных направлениях. Проекция в 90 градусах имеет более широкий профиль, чем проекция в 0 градусах из-за большой вертикальной полуоси наиболее удаленного эллипса фантома.

imagesc(theta3,xp,R3)
colormap(hot)
colorbar
xlabel('Parallel Rotation Angle - \theta (degrees)'); 
ylabel('Parallel Sensor Position - x\prime (pixels)');

Параллельный луч - восстанавливает главный фантом из данных о проекции

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

Следующие три реконструкции (I1, I2, и I3) покажите эффект варьирования количества углов, под которыми сделаны проекции. Для I1 и I2 некоторые функции, которые отобразились в исходном фантоме, не ясны. А именно, посмотрите на эти три замещающих знака в нижней части каждого изображения. Результат в I3 тесно напоминает оригинальное изображение, P.

Заметьте значительные артефакты, существующие в I1 и I2. Чтобы избежать этих артефактов, используйте большее число углов.

% Constrain the output size of each reconstruction to be the same as the
% size of the original image, |P|.
output_size = max(size(P));

dtheta1 = theta1(2) - theta1(1);
I1 = iradon(R1,dtheta1,output_size);

dtheta2 = theta2(2) - theta2(1);
I2 = iradon(R2,dtheta2,output_size);

dtheta3 = theta3(2) - theta3(1);
I3 = iradon(R3,dtheta3,output_size);

figure
montage({I1,I2,I3},'Size',[1 3])
title('Reconstruction from Parallel Beam Projection with 18, 24, and 90 Projection Angles')

Луч вентилятора - вычисляет синтетические проекции

Вычислите синтетические проекции с помощью геометрии луча вентилятора и варьируйтесь 'FanSensorSpacing'.

D = 250; 
dsensor1 = 2;
F1 = fanbeam(P,D,'FanSensorSpacing',dsensor1);

dsensor2 = 1;
F2 = fanbeam(P,D,'FanSensorSpacing',dsensor2);

dsensor3 = 0.25;
[F3, sensor_pos3, fan_rot_angles3] = fanbeam(P,D,...
                                             'FanSensorSpacing',dsensor3);

Отобразите данные о проекции F3. Заметьте, что диапазон углов поворота вентилятора от 0 до 360 градусов и тех же шаблонов происходит при смещении 180 градусов, потому что те же функции производятся с обеих сторон. Можно коррелировать функции в этом изображении проекций луча вентилятора с теми же функциями в изображении проекций параллельного луча, выше.

imagesc(fan_rot_angles3, sensor_pos3, F3)
colormap(hot)
colorbar
xlabel('Fan Rotation Angle (degrees)')
ylabel('Fan Sensor Position (degrees)')

Луч вентилятора - восстанавливает главный фантом из данных о проекции

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

Изменение значения 'FanSensorSpacing' эффективно изменяет количество датчиков, используемых при каждом угле поворота. Для каждой из этих реконструкций луча вентилятора используются те же углы поворота. Это в отличие от реконструкций параллельного луча, которые каждый использовал различные углы поворота.

Обратите внимание на то, что 'FanSensorSpacing' является только одним параметром нескольких, что можно управлять при использовании fanbeam и ifanbeam. Можно также преобразовать назад и вперед между параллелью - и данными о проекции луча вентилятора с помощью функций fan2para и para2fan.

Ifan1 = ifanbeam(F1,D,'FanSensorSpacing',dsensor1,'OutputSize',output_size);
Ifan2 = ifanbeam(F2,D,'FanSensorSpacing',dsensor2,'OutputSize',output_size);
Ifan3 = ifanbeam(F3,D,'FanSensorSpacing',dsensor3,'OutputSize',output_size);

figure
montage({Ifan1,Ifan2,Ifan3},'Size',[1 3])
title('Reconstruction from Fan Beam Projection with 18, 24, and 90 Projection Angles')