Создайте Фантом Головы

Тестовое изображение является фантомом головы Shepp-Logan, который может быть сгенерирован с помощью функции 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')