Tehnici de reconstructie a imaginilor

in tomografia computerizata

Reconstructia imaginilor folosind algoritmii

Kaczmarz-Tanabe

Cimmino extins

Anul 2012

Etapa I: Algoritmul Kaczmarz – Tanabe

Reconstructia imaginii pornind de la imaginea exacta

Imaginea reconstruita Imaginea reconstruita cu constrangeri

Imaginea cu xls Imaginea exacta

Erori

Erori fara constrangere

Erori cu constrangere

Etapa II: Algoritmul Cimmino extins

Erori

Kaczmarz – Tanabe

clear all;clc;% N=144;alpha=12;M=36;fm = fopen('a_36x144.mat');%M=144;fm = fopen('a_144x144.mat');%M=576;fm = fopen('a_576x144.mat');A = zeros(M,N);normAi2 = zeros(M,1);for i=1:M for j=1:N A(i,j)=fscanf(fm,'%f',1); end normAi2(i) = norm(A(i,:))^2;endfclose(fm);%%xex=zeros(N,1);%% Litera E% xex(29)=0.7;xex(30)=0.7;xex(31)=0.7;xex(41)=0.7;% xex(53)=0.7;xex(65)=0.7;xex(66)=0.7;xex(77)=0.7;% xex(89)=0.7;xex(101)=0.7;xex(102)=0.7;xex(103)=0.7; %Litera S Imaginea exactaxex(5)=0.7;xex(6)=0.7;xex(7)=0.7;xex(8)=0.7;xex(17)=0.7;xex(29)=0.7;xex(41)=0.7;xex(53)=0.7;xex(54)=0.7;xex(55)=0.7;xex(56)=0.7;xex(68)=0.7;xex(80)=0.7;xex(92)=0.7;xex(104)=0.7;xex(103)=0.7;xex(102)=0.7;xex(101)=0.7; %%b=A*xex;%%xls=pinv(A)*b;%%cont = 0;NITER = 1000;xk = zeros(N,1);xk_c = zeros(N,1);%%dist = zeros(NITER,1); dist_c = zeros(NITER,1);rel_err = zeros(NITER,1); rel_err_c = zeros(NITER,1);st_dev = zeros(NITER,1); st_dev_c = zeros(NITER,1);rez_err = zeros(NITER,1); rez_err_c = zeros(NITER,1);

%cursor = 0;while (cont < NITER) if (mod(cont, NITER/100) == 0) cursor = cursor + 1; disp(cursor); end cont = cont + 1; % algoritm K-T for i=1: M if (normAi2(i) == 0) break end Ai = A(i,:); xk = (xk' - ((dot(xk,Ai) - b(i))/normAi2(i))*Ai)'; xk_c = (xk_c' - ((dot(xk_c,Ai) - b(i))/normAi2(i))*Ai)'; end % end algoritm K-T % % constangere for j=1:N if (xk_c(j) <= 0) xk_c(j) = 0; elseif (xk_c(j) >= 1) xk_c(j) = 1; end end % end contrangere % % calcul erori n=144; x_ex = sum(xex)/n; x_k = sum(xk)/n; x_k_c = sum(xk_c); % disanta dist(cont) = sqrt(sum((xex-xk).^2)/sum((xex-x_ex).^2)); dist_c(cont) = sqrt(sum((xex-xk_c).^2)/sum((xex-x_ex).^2)); % eroare relativa rel_err(cont) = sum(abs(xex-xk))/sum(xex); rel_err_c(cont) = sum(abs(xex-xk_c))/sum(xex); % deviatia standard st_dev(cont) = (1/sqrt(n))*sqrt(sum((xk-x_k).^2)); st_dev_c(cont) = (1/sqrt(n))*sqrt(sum((xk_c-x_k_c).^2)); % eroare reziduala rez_err(cont) = norm(A*xk - b); rez_err_c(cont) = norm(A*xk_c - b); %end figure('Name','Reconstructie - Kaczmarz-Tanabe','NumberTitle','off'); tmp = vtm(xk,alpha,alpha);dispmatrix(tmp,sprintf('KT/144/K-T_xk.png',M));%dispmatrix(tmp,'KT/K-T_xk.png');title('Imaginea reconstruita'); tmp = vtm(xk_c,alpha,alpha);

dispmatrix(tmp,sprintf('KT/144/K-T_xk_c.png',M));%dispmatrix(tmp,'KT/K-T_xk_c.png');title('Imaginea reconstruita cu Constrangeri'); tmp = vtm(xls,alpha,alpha);dispmatrix(tmp,sprintf('KT/144/K-T_xls.png',M));%dispmatrix(tmp,'KT/K-T_xls.png');title('Imaginea solutiei de norma minima'); tmp = vtm(xex,alpha,alpha);dispmatrix(tmp,sprintf('KT/144/K-T_xex.png',M));%dispmatrix(tmp,'KT/K-T_xex.png');title('Imaginea exacta'); % %figure('Name','Erori - Kaczmarz-Tanabe','NumberTitle','off')subplot(2,2,1); plot(rez_err); title('Reziduu')subplot(2,2,2); plot(dist); title('Distanta')subplot(2,2,3); plot(rel_err); title('Eroare Relativa')subplot(2,2,4); plot(st_dev); title('Deviatie standard')saveas(gcf,sprintf('KT/144/Erori_-_Kaczmarz-Tanabe.png',M));%saveas(gcf,'KT/Erori_-_Kaczmarz-Tanabe.png'); figure('Name', 'Erori Constrangere - Kaczmarz-Tanabe','NumberTitle','off')subplot(2,2,1); plot(rez_err_c); title('Reziduu')subplot(2,2,2); plot(dist_c); title('Distanta')subplot(2,2,3); plot(rel_err_c); title('Eroare Relativa')subplot(2,2,4); plot(st_dev_c); title('Deviatie standard')saveas(gcf,sprintf('KT/144/Erori_Constrangere_-_Kaczmarz-Tanabe.png',M));%saveas(gcf,'KT/Erori_Constrangere_-_Kaczmarz-Tanabe.png');

Cimmino Extins

clear all;clc;% M=144;N=144;alpha=12;fm = fopen('a_144x144.mat');A = zeros(M,N);normAi2 = zeros(N,1);normAj2 = zeros(1,M);for i=1:M for j=1:N A(i,j)=fscanf(fm,'%f',1); end normAi2(i) = norm(A(i,:))^2;endfclose(fm);%for j=1:N normAj2(j) = norm(A(:,j))^2;end

%xex=zeros(N,1); % Litera E% xex(29)=0.7;xex(30)=0.7;xex(31)=0.7;xex(41)=0.7;% xex(53)=0.7;xex(65)=0.7;xex(66)=0.7;xex(77)=0.7;% xex(89)=0.7;xex(101)=0.7;xex(102)=0.7;xex(103)=0.7; %Litera S Imaginea exactaxex(5)=0.7;xex(6)=0.7;xex(7)=0.7;xex(8)=0.7;xex(17)=0.7;xex(29)=0.7;xex(41)=0.7;xex(53)=0.7;xex(54)=0.7;xex(55)=0.7;xex(56)=0.7;xex(68)=0.7;xex(80)=0.7;xex(92)=0.7;xex(104)=0.7;xex(103)=0.7;xex(102)=0.7;xex(101)=0.7; %b=A*xex;%%xls=pinv(A)*b;%%cont = 0;NITER = 1000;xk = zeros(N,1);xk_c= zeros (N,1);rand('state',0);p=rand(M,1);p=p/norm(p);bb = b+p*0.10;y=bb;omega = ones(1,M); Omega = sum(omega);gamma = ones(1,M); Gamma = sum(gamma); %%dist = zeros(NITER,1);rel_err = zeros(NITER,1);st_dev = zeros(NITER,1);rez_err = zeros(NITER,1);% while (cont < NITER) cont = cont + 1; disp(cont) %algoritm Cimino extins sum_y = zeros(1,M); for j=1:N Aj = A(:,j); sum_y = sum_y + (gamma(j)/Gamma)*(y - 2*(dot(y,Aj)/normAj2(j))*Aj)'; end y = sum_y'; c = bb - y; sum_x = zeros(N,1); for i=1:M Ai = A(i,:);

sum_x = sum_x + (omega(i)/Omega)*(xk' - 2*((dot(xk,Ai) - c(i))/normAi2(j))*Ai)'; end xk = sum_x; %end algoritm Cimino extins % % constangere for j=1:N if (xk(j) < 0) xk(j) = 0; elseif (xk(j) > 1) xk(j) = 1; end end % end contrangere % calcul erori n=144; x_ex = sum(xex)/n; x_k = sum(xk)/n; x_k_c = sum(xk_c); % distance dist(cont) = sqrt(sum((xex-xk).^2)/sum((xex-x_ex).^2)); dist_c(cont) = sqrt(sum((xex-xk_c).^2)/sum((xex-x_ex).^2)); % relative error rel_err(cont) = sum(abs(xex-xk))/sum(xex); rel_err_c(cont) = sum(abs(xex-xk_c))/sum(xex); % standart deviation st_dev(cont) = (1/sqrt(n))*sqrt(sum((xk-x_k).^2)); st_dev_c(cont) = (1/sqrt(n))*sqrt(sum((xk_c-x_k).^2)); % residual error rez_err(cont) = norm(A*xk - b); rez_err_c(cont) = norm(A*xk_c - b); end tmp = vtm(xk,alpha,alpha);dispmatrix(tmp,'CIM/144/CIM_xk.png');title('Cimmino - Imaginea reconstruita'); tmp = vtm(xk_c,alpha,alpha);dispmatrix(tmp,'CIM/144/CIM_xk_c.png');title('Cimmino - Imaginea reconstruita cu constrangeri'); tmp = vtm(xls,alpha,alpha);dispmatrix(tmp,'CIM/144/CIM_xls.png');title('Cimmino - Imaginea solutiei de norma minima'); tmp = vtm(xex,alpha,alpha);dispmatrix(tmp,'CIM/144/CIM_xex.png');title('Cimmino - Imaginea exacta'); figure('Name','Erori - Cimmino','NumberTitle','off')subplot(2,2,1); plot(rez_err); title('Reziduu')subplot(2,2,2); plot(dist); title('Distanta')subplot(2,2,3); plot(rel_err); title('Eroare Relativa')

subplot(2,2,4); plot(st_dev); title('Deviatie standard')saveas(gcf,'CIM/Erori_-_Cimmino.png'); figure('Name', 'Erori Constrangere - Cimmino','NumberTitle','off')subplot(2,2,1); plot(rez_err_c); title('Reziduu')subplot(2,2,2); plot(dist_c); title('Distanta')subplot(2,2,3); plot(rel_err_c); title('Eroare Relativa')subplot(2,2,4); plot(st_dev_c); title('Deviatie standard')saveas(gcf,'CIM/Erori_Constrangere_-_Cimmino.png');