ĐÀM ĐỨc cưỜng ứng dụng kỹ thuật kết hợp tần số nhằm nâng cao chất lưỢng ảnh siêU Âm cắt lớP



tải về 362.16 Kb.
trang6/6
Chuyển đổi dữ liệu29.10.2017
Kích362.16 Kb.
1   2   3   4   5   6

PHỤ LỤC 2: CODE MATLAB DBIM ĐỀ XUẤT


Ở phương pháp đề xuất này vẫn sử dụng hàm chính như ở Phụ Lục 1, nhưng ta thực hiện DBIM 2 lần với tần số f1 và f2 vì thế ta có chương trình cho phương pháp đề xuất như sau:

Phương pháp đề xuất

clear,clc,close all;

f0=2e6;% la tan so f2

f=1e6;% tuong ung la f1

N=22

Niter=4 %la so lan lap x



Niter1=4%so lan lap cua f2
co=1540; % m/s

c1=co*(1-0.02); % Khoa hoc tre

landa=co/f; % m/s

lamda=co/f0;

ko=2*pi/landa; % increase when f is increased

D=4*lamda

h=D/(N-1);

x0=[(N+1)/2;(N+1)/2];

fprintf('Ratio landa/h = %f \n',landa/h);
%tao muc ham muc tieu ly tuong

SC=zeros(N,N);

SC1=SC;

xo=(N+1)/2;yo=(N+1)/2;



for m=1:N

for n=1:N

dis=sqrt((xo-m)^2+(yo-n)^2);

%if dis

if dis>3.5355*N/10

SC(m,n)=0;

else

SC(m,n)=(2*pi*f0)^2*(1/(c1^2)-1/(co^2));



%SC(m,n)=dis;

end;


end;

end;


figure;

[X,Y] = meshgrid(linspace(0,D/landa,length(SC)));

mesh(X,Y,SC);

xlabel('\lambda')

ylabel('\lambda')
RRE=2^(-4);

%RRE=2^(-6);

L=N*N; % number of detector

phi=linspace(-pi,pi,L);

No=10*N; % at first N=11, No can be changed when changing N

% however it is affected to distance from tranceivers to object

% because distance=No*h= constant

K2=cos(phi)*(No+.5)+x0(2);

K1=sin(phi)*(1.5*No-.5) + x0(1);

KK2=K2;KK1=K1;

%noise_flag: option 0: noise init, 2: no noise

noise_flag=0;


transmiter=1:N:L;%may phat

detector=1:2*N:L;%may thu


plot(KK2(detector),KK1(detector),'s')

hold on


plot(K2(transmiter),K1(transmiter),'r*')

hold on;


mesh(abs(SC))

legend('detector','transmiter','scatter area')

calculate_PINC_matrix_cavichi;

figure(100)

subplot(211)

imagesc(abs(SC))

subplot(212)

mesh(abs(SC))

Niter;

err=zeros(1,Niter);



MSE=zeros(1,Niter);

PSNR=zeros(1,Niter);

gama=zeros(1,Niter);

for iter=1:Niter


iter

Mt=[];


delta_sc_t=[];

X=[];


p_sc_exact_t=[];

p_sc_t=[];

calculate_B_matrix_DBIM;

calculate_C_matrix_DBIM;

% update C matrix for each iteration

uu=-1; % transmit index

for dec1=transmiter % for each transmiter and detector

uu=uu+1;


u=0; % detector index, reset for each detector

for dec2=detector %

% create data for a single transmision and receving

% otain the p_sc_exact : scatter field measured at detector; (1 value)

% p_sc : scatter field predicted at detector; (1 value)

% p: presure calculated (predict at each pixel);matrix (NxN)

B_SC=BB_SC(1+u*N:N+u*N,:);

B=BB(1+u*N:N+u*N,:);

p_inc=PINC(1+uu*N:N+uu*N,:);

u3=0; % index cho C, reset for each new pixel

for n1=1:N % di vao tung pixel, xac dinh boi 2D: n1,n2

for n2=1:N

C_SC=CC_SC(1+u3*N:N+u3*N,:);

C_SC1=CC_SC1(1+u3*N:N+u3*N,:);

pp(n1,n2)=sum(sum(C_SC.*SC.*p_inc)); % C va p_inc trong moi se khac

pp1(n1,n2)=sum(sum(C_SC1.*SC1.*p_inc));

u3=u3+1;

end;


end;

p =p_inc+pp; % pp<>0; wave equation; % IDEAL

p1=p_inc+pp1; % at the first step p_inc=p1 because SC=zeros, pp1=0; PRDEICT

p_sc_exact=sum(sum(B.*SC.*p)); % using p matrix, IDEAL,1 point, IDEAL

p_sc=sum(sum(B.*SC1.*p)); % using p matrix, Predict,1 point

u=u+1;


% B change depends on detector

delta_sc=p_sc_exact-p_sc; % 1 valuse

M=reshape(B.*p1,1,N*N);

Mt=[Mt;M]; % add more detector

p_sc_exact_t=[p_sc_exact_t ; p_sc_exact]; % NCG

p_sc_t=[p_sc_t p_sc];

delta_sc_t=[delta_sc_t;delta_sc]; % add more scatter field in detector

end; % end one cycle of transmit abd detect

end;

if noise_flag==0;



n=0.05*sqrt(var(delta_sc_t))*randn(size(delta_sc_t));

noise_flag=1;

end;

if noise_flag==2;



n=zeros(size(delta_sc_t));

end;


% add noise

if iter==1

delta_sc_t=delta_sc_t+n;

end;


%%%%%%%

[n1,n2]=size(Mt);

gama=1.8755e-021; % co dinh gama

%%%%%%%%


delta_sound = test_NCG(Mt,delta_sc_t,1e3,RRE,gama);

SC1=SC1+reshape(delta_sound,N,N); % update the sound contrast

figure;

subplot(211);imagesc(abs(SC1))



subplot(212);mesh(abs(SC1))

end; % interation for convergen of sound contrast

ko=2*pi/lamda; % increase when f is increased

calculate_PINC_matrix_cavichi;

Niter1;

err1=zeros(1,Niter1);



MSE1=zeros(1,Niter1);

PSNR1=zeros(1,Niter1);

gama1=zeros(1,Niter1);

Quality=zeros(1,Niter1);

for iter1=1:Niter1
iter

%iter1=iter1+1 %use for while

Mt=[];

delta_sc_t=[];



X=[];

p_sc_exact_t=[];

p_sc_t=[];

calculate_B_matrix_DBIM;

calculate_C_matrix_DBIM;

% update C matrix for each iteration

uu=-1; % transmit index

for dec1=transmiter % for each transmiter and detector

uu=uu+1;

u=0; % detector index, reset for each detector

for dec2=detector %

% create data for a single transmision and receving

% otain the p_sc_exact : scatter field measured at detector; (1 value)

% p_sc : scatter field predicted at detector; (1 value)

% p: presure calculated (predict at each pixel);matrix (NxN)

B_SC=BB_SC(1+u*N:N+u*N,:);

B=BB(1+u*N:N+u*N,:);

p_inc=PINC(1+uu*N:N+uu*N,:);

u3=0; % index cho C, reset for each new pixel

for n1=1:N % di vao tung pixel, xac dinh boi 2D: n1,n2

for n2=1:N

C_SC=CC_SC(1+u3*N:N+u3*N,:);

C_SC1=CC_SC1(1+u3*N:N+u3*N,:);

pp(n1,n2)=sum(sum(C_SC.*SC.*p_inc)); % C va p_inc trong moi se khac

pp1(n1,n2)=sum(sum(C_SC1.*SC1.*p_inc));

u3=u3+1;


end;

end;


p =p_inc+pp; % pp<>0; wave equation; % IDEAL

p1=p_inc+pp1; % at the first step p_inc=p1 because SC=zeros, pp1=0; PRDEICT

p_sc_exact=sum(sum(B.*SC.*p)); % using p matrix, IDEAL,1 point, IDEAL

p_sc=sum(sum(B.*SC1.*p)); % using p matrix, Predict,1 point

u=u+1;

% B change depends on detector



delta_sc=p_sc_exact-p_sc; % 1 valuse

M=reshape(B.*p1,1,N*N);

Mt=[Mt;M]; % add more detector

p_sc_exact_t=[p_sc_exact_t ; p_sc_exact]; % NCG

p_sc_t=[p_sc_t p_sc];

delta_sc_t=[delta_sc_t;delta_sc]; % add more scatter field in detector

end; % end one cycle of transmit abd detect

end;


if noise_flag==0;

n=0.05*sqrt(var(delta_sc_t))*randn(size(delta_sc_t))

noise_flag=1;

end;


if noise_flag==2;

n=zeros(size(delta_sc_t));

end;

% add noise



if iter==1

delta_sc_t=delta_sc_t+n;

end;

%%%%%%%


[n1,n2]=size(Mt);

gama=1.8755e-021; % co dinh gama

delta_sound = test_NCG(Mt,delta_sc_t,1e3,RRE,gama);
SC1=SC1+reshape(delta_sound,N,N); % update the sound contrast

figure;


subplot(211);imagesc(abs(SC1))

subplot(212);mesh(abs(SC1))


err1(iter1)=sum(sum(abs(SC-abs(SC1))))/sum(sum(SC));

MSE1(iter1)=(1/N^2)*sum(sum(abs(SC-abs(SC1))^2));


% tinh chi so Quality

[nn1,nn2]=size(SC1);

N=nn1

xx1=(1/N^2)*sum(sum(SC));



xx=xx1*ones(N,N);

yy1=(1/N^2)*sum(sum(abs(SC1)));

yy=xx1*ones(N,N);

detal_y=sum(sum((abs(SC1)-yy)^2/(N^2-1)));

detal_x=sum(sum((abs(SC)-xx)^2/(N^2-1)));

detal_xy=(1/(N^2-1))*sum(sum((abs(SC)-xx)*(abs(SC1)-yy)));

Q=4*detal_xy*xx1*yy1/((detal_x+detal_y)*(xx1^2+yy1^2));

Quality(iter1)=Q;

%het phan tinh chi so

%PSNR1(iter)=10*log10((2^N-1)^2/MSE1);

end;

newSC1=SC1(floor(length(SC1)/2),1:length(SC1));



newSC=SC(floor(length(SC)/2),1:length(SC));

figure()


plot(newSC);

hold on;


plot(abs(newSC1),'r');

figure


plot(1:Niter1,err1,'-ro');

xlabel('interation');ylabel('error');



figure

plot(1:Niter1,MSE1,'-ro');



xlabel('interation');ylabel('MSE');


1   2   3   4   5   6


Cơ sở dữ liệu được bảo vệ bởi bản quyền ©tieuluan.info 2017
được sử dụng cho việc quản lý

    Quê hương