ID: 6021

Название работы: Прямое и обратное преобразование Радона

Категория: Лабораторная работа

Предметная область: Коммуникация, связь, радиоэлектроника и цифровые приборы

Описание: Цель работы: Ознакомление с прямым и обратным преобразованием Радона изображений. Реализация прямого и обратного преобразований Радона с помощью функций в среде MatLAB (ImageProcessingToolbox). Короткие теоретические сведения ...

Язык: Русский

Дата добавления: 2012-12-27

Размер файла: 1.24 MB

Работу скачали: 99 чел.

Цель работы: Ознакомление с прямым и обратным преобразованием Радона изображений. Реализация прямого и обратного преобразований Радона  с помощью функций в среде MatLAB (Image Processing Toolbox).

Короткие теоретические сведения

 Преобразование Радона.

С помощью преобразования Радона изображение представляется в виде набора проекций вдоль различных направлений. В результате получается совокупность теней, т.е. трёхмерная структура объекта сводится к набору двумерных изображений. При этом, проекция функции двух переменных f(x,y) представляет собой интеграл в определённом направлении.  Например, интеграл от f(x,y) в вертикальном направлении является проекцией  f(x,y) на ось x; интеграл в горизонтальном направлении является проекцией на ось  y.

Рис.1

Проекции могут быть вычислены вдоль любого угла  θ. Так, проекция  функции двух переменных  f(x,y) на ось  задаётся  интегралом

где оси   и   задаются поворотом против часовой стрелки на угол θ с использованием следующего выражения:

Геометрическое представление  преобразования Радона приведено на рисунке 2.

Рис.2

Преобразование Радона для большого количества углов чаще всего отображается в виде изображения. Например, преобразование Радона для прямоугольника при изменении θ от  0 до 180° с шагом 1° имеет вид (рис.3):

Рис.3

Обратное преобразование Радона реконструирует изображение по его матрице проекций. В компьютерной томографии (рис.4) осуществляется восстановление изображения сечения человеческого тела с использованием облучения рентгеновскими лучами под различными углами. Задача восстановления f(x,y) сводится к решению конечного числа уравнений  при различных значениях угла θ.

Рис.4

В большинстве случаев, не существует исходного изображения, от которого получают проекции. Например, при томографии, проекции формируются путём измерения интенсивности излучения, проходящего через физический объект под различными углами. Значения проекций накапливаются в специальном оборудовании, а затем с помощью обратной функции Радона выполняется неинвазивное (без вторжения во внутрь) восстановление внутренней структуры объекта (человека).

 

Порядок работы

1.  Исследования преобразований Радона на модели

С помощью функции phantom сымитировать срез головы человека. Выполнить прямое преобразование Радона (radon), выбрав различные углы. Восстановить изображение с помощью обратного преобразования Радона. Исследовать влияние параметров команды iradon на сходство восстановленного изображения с оригиналом, а также на время выполнения обратного преобразования. Оценить качество восстановленного изображения.

close all; clear all; clc;

P=phantom(256); %Создание искусственного изображения среза головы

figure(1)

imagesc(P);

title('Head Phantom');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

%Прямые преобразования Радона с различными углами theta

theta1 = 0:10:170;

[R1,xp1] = radon(P,theta1);

theta2 = 0:5:175;

[R2,xp2] = radon(P,theta2);

theta3 = 0:2:178;

[R3,xp3] = radon(P,theta3);

theta4 = 0:0.1:179;

[R4,xp4] = radon(P,theta4);

figure(2)

subplot(2,2,1)

imagesc(theta1,xp1,R1)

title('Num angles = 18');

xlabel('Rotation Angle - \theta (degrees)');

ylabel('Sensor Position - x\prime (pixels)');

colormap(hot)

subplot(2,2,2)

imagesc(theta2,xp2,R2)

title('Num angles = 36');

xlabel('Rotation Angle - \theta (degrees)');

ylabel('Sensor Position - x\prime (pixels)');

colormap(hot)

subplot(2,2,3)

imagesc(theta3,xp3,R3)

title('Num angles = 90');

xlabel('Rotation Angle - \theta (degrees)');

ylabel('Sensor Position - x\prime (pixels)');

colormap(hot)

subplot(2,2,4)

imagesc(theta4,xp4,R4)

title('Num angles = 1791');

xlabel('Rotation Angle - \theta (degrees)');

ylabel('Sensor Position - x\prime (pixels)');

colormap(hot)

%Вычисление обратного преобразования Радона с различными углами theta

IR1=iradon(R1,theta1,'nearest','Ram-Lak',1,256);

IR2=iradon(R2,theta2,'nearest','Ram-Lak',1,256);

IR3=iradon(R3,theta3,'nearest','Ram-Lak',1,256);

IR4=iradon(R4,theta4,'nearest','Ram-Lak',1,256);

figure(3)

subplot(2,2,1)

imagesc(theta1,xp1,IR1)

title('Num angles = 18');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,2,2)

imagesc(theta2,xp2,IR2)

title('Num angles = 36');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,2,3)

imagesc(theta3,xp3,IR3)

title('Num angles = 90');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,2,4)

imagesc(theta4,xp4,IR4)

title('Num angles = 1791');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

%Среднеквадратическая относительная погрешность (18 углов)

close all; clear all; clc;

P=phantom(256);

theta1 = 0:10:170;

[R1,xp1] = radon(P,theta1);

IR1=iradon(R1,theta1,'nearest','Ram-Lak',1,256);

d1_18=((abs(sum(sum((P-IR1).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность (18 углов)

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_18=((abs(sum(sum((P-IR1).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность (18 углов)

d3_18=max(max(abs(P-IR1)))

%Максимальная относительная погрешность (18 углов)

d4_18=((max(max(abs(P-IR1))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность (18 углов)

d5_18=((sum(sum(abs(P-IR1))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность (18 углов)

d6_18=(sum(sum(abs(P-IR1))))

%Среднеквадратическая абсолютная погрешность (18 углов)

d7_18=sqrt(sum(sum((P-IR1).^2)))

Погрешность	Количество углов
		18	36	90	1791
d1 %	96.0877	60.1466	37.0012	32.1464
d2 %	110.8376	69.3794	42.6810	37.0810
d3	1.4261	1.1539	1.0088	0.9241
d4 %	142.6084	115.3943	100.8821	92.4145
d5 %	131.8101	81.1806	38.1935	17.2767
d6	1.0603e+004	6.5302e+003	3.0723e+003	1.3897e+003
d7	60.5740	37.9166	23.3256	20.2652

%Вычисление обратного преобразования Радона с различными типами интерполяции

close all; clear all; clc;

P=phantom(256);

theta = 0:0.5:179;

[R,xp] = radon(P,theta);

IR1=iradon(R,theta,'nearest','Ram-Lak',1,256);

%Среднеквадратическая относительная погрешность

d1_nearest=((abs(sum(sum((P-IR1).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_nearest=((abs(sum(sum((P-IR1).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_nearest=max(max(abs(P-IR1)))

%Максимальная относительная погрешность

d4_nearest=((max(max(abs(P-IR1))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_nearest=((sum(sum(abs(P-IR1))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_nearest=(sum(sum(abs(P-IR1))))

%Среднеквадратическая абсолютная погрешность

d7_nearest=sqrt(sum(sum((P-IR1).^2)))

IR2=iradon(R,theta,'linear','Ram-Lak',1,256);

%Среднеквадратическая относительная погрешность

d1_linear=((abs(sum(sum((P-IR2).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_linear=((abs(sum(sum((P-IR2).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_linear=max(max(abs(P-IR2)))

%Максимальная относительная погрешность

d4_linear=((max(max(abs(P-IR2))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_linear=((sum(sum(abs(P-IR2))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_linear=(sum(sum(abs(P-IR2))))

%Среднеквадратическая абсолютная погрешность

d7_linear=sqrt(sum(sum((P-IR2).^2)))

IR3=iradon(R,theta,'spline','Ram-Lak',1,256);

%Среднеквадратическая относительная погрешность

d1_spline=((abs(sum(sum((P-IR3).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_spline=((abs(sum(sum((P-IR3).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_spline=max(max(abs(P-IR3)))

%Максимальная относительная погрешность

d4_spline=((max(max(abs(P-IR3))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_spline=((sum(sum(abs(P-IR3))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_spline=(sum(sum(abs(P-IR3))))

%Среднеквадратическая абсолютная погрешность

d7_spline=sqrt(sum(sum((P-IR3).^2)))

figure(1)

subplot(1,3,1)

imagesc(theta,xp,IR1)

title('Интерполяция по ближайшей окрестности');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(1,3,2)

imagesc(theta,xp,IR2)

title('Линейная интерполяция');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(1,3,3)

imagesc(theta,xp,IR3)

title('Сплайновая интерполяция');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

Погрешность	Тип интерполяции
		nearest	linear	spline
d1 %	32.4182	31.2885	31.6798
d2 %	37.3945	36.0914	36.5428
d3	0.9206	0.8663	0.9012
d4 %	92.0633	86.6276	90.1225
d5 %	20.0611	15.4540	15.1226
d6	1.6137e+003	1.2431e+003	1.2165e+003
d7	20.4365	19.7244	19.9711

%Вычисление обратного преобразования Радона с различными типами фильтра

close all; clear all; clc;

P=phantom(256);

theta = 0:0.5:179;

[R,xp] = radon(P,theta);

IR1=iradon(R,theta,'linear','Ram-Lak',1,256);

%Среднеквадратическая относительная погрешность

d1_Ram=((abs(sum(sum((P-IR1).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_Ram=((abs(sum(sum((P-IR1).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_Ram=max(max(abs(P-IR1)))

%Максимальная относительная погрешность

d4_Ram=((max(max(abs(P-IR1))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_Ram=((sum(sum(abs(P-IR1))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_Ram=(sum(sum(abs(P-IR1))))

%Среднеквадратическая абсолютная погрешность

d7_Ram=sqrt(sum(sum((P-IR1).^2)))

IR2=iradon(R,theta,'linear','Shepp-Logan',1,256);

%Среднеквадратическая относительная погрешность

d1_Shepp=((abs(sum(sum((P-IR2).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_Shepp=((abs(sum(sum((P-IR2).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_Shepp=max(max(abs(P-IR2)))

%Максимальная относительная погрешность

d4_Shepp=((max(max(abs(P-IR2))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_Shepp=((sum(sum(abs(P-IR2))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_Shepp=(sum(sum(abs(P-IR2))))

%Среднеквадратическая абсолютная погрешность

d7_Shepp=sqrt(sum(sum((P-IR2).^2)))

IR3=iradon(R,theta,'linear','Cosine',1,256);

%Среднеквадратическая относительная погрешность

d1_Cosine=((abs(sum(sum((P-IR3).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_Cosine=((abs(sum(sum((P-IR3).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_Cosine=max(max(abs(P-IR3)))

%Максимальная относительная погрешность

d4_Cosine=((max(max(abs(P-IR3))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_Cosine=((sum(sum(abs(P-IR3))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_Cosine=(sum(sum(abs(P-IR3))))

%Среднеквадратическая абсолютная погрешность

d7_Cosine=sqrt(sum(sum((P-IR3).^2)))

IR4=iradon(R,theta,'linear','Hamming',1,256);

%Среднеквадратическая относительная погрешность

d1_Hamming=((abs(sum(sum((P-IR4).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_Hamming=((abs(sum(sum((P-IR4).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_Hamming=max(max(abs(P-IR4)))

%Максимальная относительная погрешность

d4_Hamming=((max(max(abs(P-IR4))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_Hamming=((sum(sum(abs(P-IR4))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_Hamming=(sum(sum(abs(P-IR4))))

%Среднеквадратическая абсолютная погрешность

d7_Hamming=sqrt(sum(sum((P-IR4).^2)))

IR5=iradon(R,theta,'linear','Hann',1,256);

%Среднеквадратическая относительная погрешность

d1_Hann=((abs(sum(sum((P-IR5).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_Hann=((abs(sum(sum((P-IR5).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_Hann=max(max(abs(P-IR5)))

%Максимальная относительная погрешность

d4_Hann=((max(max(abs(P-IR5))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_Hann=((sum(sum(abs(P-IR5))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_Hann=(sum(sum(abs(P-IR5))))

%Среднеквадратическая абсолютная погрешность

d7_Hann=sqrt(sum(sum((P-IR5).^2)))

figure(1)

subplot(2,3,1)

imagesc(theta,xp,IR1)

title('Усеченный фильтр Рама-Лака');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,3,2)

imagesc(theta,xp,IR2)

title('Фильтр Шеппа-Логана');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,3,3)

imagesc(theta,xp,IR3)

title('Косинусный фильтр');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,3,4)

imagesc(theta,xp,IR4)

title('Фильтр Хэмминга');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(2,3,5)

imagesc(theta,xp,IR5)

title('Фильтр Ханна');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)'); 

Погрешность	Тип фильтра
		Рама-Лака	Шеппа-Логана	Косинусный	Хэмминга	Ханна
d1 %	31.2885	31.1055	31.0972	31.5070	31.6213
d2 %	36.0914	35.8803	35.8707	36.3435	36.4753
d3	0.8663	0.8476	0.8124	0.7917	0.7852
d4 %	86.6276	84.7566	81.2392	79.1673	78.5186
d5 %	15.4540	15.3622	15.8115	16.9154	17.1456
d6	1.2431e+003	1.2357e+003	1.2719e+003	1.3607e+003	1.3792e+003
d7	19.7244	19.6090	19.6037	19.8621	19.9341

%Вычисление обратного преобразования Радона с различными параметрами d (сдвиг по частотной области)

close all; clear all; clc;

P=phantom(256);

theta = 0:0.5:179;

[R,xp] = radon(P,theta);

IR1=iradon(R,theta,'linear','Cosine',0.1,256);

%Среднеквадратическая относительная погрешность

d1_01=((abs(sum(sum((P-IR1).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_01=((abs(sum(sum((P-IR1).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_01=max(max(abs(P-IR1)))

%Максимальная относительная погрешность

d4_01=((max(max(abs(P-IR1))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_01=((sum(sum(abs(P-IR1))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_01=(sum(sum(abs(P-IR1))))

%Среднеквадратическая абсолютная погрешность

d7_01=sqrt(sum(sum((P-IR1).^2)))

IR2=iradon(R,theta,'linear','Cosine',0.5,256);

%Среднеквадратическая относительная погрешность

d1_05=((abs(sum(sum((P-IR2).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_05=((abs(sum(sum((P-IR2).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_05=max(max(abs(P-IR2)))

%Максимальная относительная погрешность

d4_05=((max(max(abs(P-IR2))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_05=((sum(sum(abs(P-IR2))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_05=(sum(sum(abs(P-IR2))))

%Среднеквадратическая абсолютная погрешность

d7_05=sqrt(sum(sum((P-IR2).^2)))

IR3=iradon(R,theta,'linear','Cosine',0.9,256);

%Среднеквадратическая относительная погрешность

d1_09=((abs(sum(sum((P-IR3).^2))/sum(sum(P.^2))))^0.5)*100

%Нормализованная среднеквадратическая погрешность

P_sr=(sum(sum(P)))/(size(P,1)*size(P,2));

d2_09=((abs(sum(sum((P-IR3).^2))/sum(sum((P-P_sr).^2))))^0.5)*100

%Максимальная погрешность

d3_09=max(max(abs(P-IR3)))

%Максимальная относительная погрешность

d4_09=((max(max(abs(P-IR3))))/(max(max(P))))*100

%Нормализированная абсолютная средняя погрешность

d5_09=((sum(sum(abs(P-IR3))))/(sum(sum(abs(P)))))*100

%Абсолютная средняя погрешность

d6_09=(sum(sum(abs(P-IR3))))

%Среднеквадратическая абсолютная погрешность

d7_09=sqrt(sum(sum((P-IR3).^2)))

figure(1)

subplot(1,3,1)

imagesc(theta,xp,IR1)

title('D=0.1');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(1,3,2)

imagesc(theta,xp,IR2)

title('D=0.5');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

subplot(1,3,3)

imagesc(theta,xp,IR3)

title('D=0.9');

xlabel('x\prime (pixels)');

ylabel('y\prime (pixels)');

Погрешность	Сдвиг по частотной области (параметр D)
		D=0.1	D=0.5	D=0.9
d1 %	61.8194	33.7143	31.2002
d2 %	71.3089	38.8895	35.9896
d3	0.7558	0.7386	0.8032
d4 %	75.5842	73.8645	80.3183
d5 %	53.7335	20.3204	16.1436
d6	4.3223e+003	1.6346e+003	1.2986e+003
d7	38.9711	21.2536	19.6687

1.  Преобразование реального биомедицинского изображения

Выполнить прямое и обратное преобразование Радона, используя реальное биомедицинское изображение.

close all; clear all; clc;

Ima=rgb2gray(imread('C:\24','jpg'));

Im=double(Ima);

figure(1)

imshow(Ima)

theta = 0:0.5:179;

[R_Im, xp]=radon(Im,theta);

figure(5)

imagesc(theta,xp,R_Im)

IR_Im=iradon(R_Im,theta,'linear','Cosine',1);

figure(6)

imagesc(theta,xp,IR_Im)

PAGE  2