KẾT LUẬN
Luận văn này đã thành công trong việc nâng cao chất lượng ảnh chụp siêu âm cắt lớp bằng cách sử dụng kết hợp 2 tần số f1 và f2. Ảnh khôi phục theo phương pháp đề xuất cho chất lượng tốt hơn ảnh theo phương pháp truyền thống.
Tác giả cũng đã tìm được số bước lặp tối ưu với f1 sao cho việc kết hợp f1 và f2 cho chất lượng tốt nhất.
Đánh giá được tham số chất lượng Q được trình bày ở phần 2.3. Từ đó kết luận được ảnh tái tạo bởi việc sử dụng kết hợp 2 tần số f1 và f2, cho kết quả đánh giá về mặt sai số toán học thông dụng hay có xét cả đến vấn đề VHS (visual human system) đều tốt hơn so với chỉ sử dụng một tần số đơn.
Như vậy việc sử dụng kết hợp 2 tần số trong việc cải thiện chất lượng ảnh đã thành công, tạo điều kiện áp dụng trong lĩnh vực Y – Sinh. Bước tiếp theo của đề xuất này là việc thử nghiệm đề xuất trong tạo ảnh với những dữ liệu thực tế để có thể áp dụng trong ngành chuẩn đoán y khoa.
TÀI LIỆU THAM KHẢO
[1] C. F. Schueler, H. Lee, and G. Wade, “Fundamentals of digital ultrasonic processing,” IEEE Transactions on Sonics and Ultrasonics, vol. 31, no. 4, pp. 195–217, July 1984.
[2] N. Duric, P. Littrup, A. Babkin, D. Chambers, S. Azevedo, A. Kalinin, R.Pevzner, M. Tokarev, E. Holsapple, O. Rama, and R. Duncan, “Development of ultrasound tomography for breast imaging: Technical assessment,” Medical Physics, vol. 32, no. 5, pp. 1375–1386, May 2005.
[3] J.-W. Jeong, T.-S. Kim, D. C. Shin, S. Do, M. Singh, and V. Z. Marmarelis, “Soft tissue differentiation using multiband signatures of high resolution ul-trasonic transmission tomography,” IEEE Transactions on Medical Imaging, vol. 24, no. 3, pp. 399–408, March 2005.
[4] S. A. Johnson, T. Abbott, R. Bell, M. Berggren, D. Borup, D. Robinson, J. Wiskin, S. Olsen, and B. Hanover, “Noninvasive breast tissue charac-terization using ultrasound speed and attenuation,” in Acoustical Imaging, vol. 28, 2007, pp. 147–154.
[5] C. Li, N. Duric, and L. Huang, “Breast imaging using transmission ultra-sound: Reconstructing tissue parameters of sound speed and attenuation,” in International Conference on BioMedical Engineering and Informatics, vol. 2, 2008, pp. 708–712.
[6] R. J. Lavarello and M. L. Oelze: Tomographic Reconstruction of Three-Dimensional Volumes Using the Distorted Born Iterative Method. IEEE Transactions on Medical Imaging, 28, 2009, pp. 1643-1653.
[7] Lavarello Robert: New Developments on Quantitative Imaging Using Ultrasonic Waves. University of Illinois at Urbana-Champaign, 2009.
[8] http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method
[9] M. T. Heath, Scientific Computing: An Introductory Survey. New York, NY: McGraw-Hill, 2002.
[10] Martin, R., Noise power spectral density estimation based on optimal smoothing and minimum statistics, IEEE Transactions on Speech and Audio Processing, Vol. 9, 2001, pp. 504 - 512.
[11] http://www-stat.stanford.edu/~susan/courses/s60/split/node60.html
[12] Tran Duc Tan, N. Linh-Trung, M. L. Oelze, M. N. Do, Application of L1 regularization for high-quality reconstruction of ultrasound tomography, International Federation for Medical and Biological Engineering (IFMBE), NXB SPRINGER, ISSN: 1680-0737, Volume 40, 2013, pp. 309-312.
[13] Tran Duc Tan, Nguyen Linh-Trung, Minh N. Do, Modified Distorted Born Iterative Method for Ultrasound Tomography by Random Sampling, The 12th International Symposium on Communications and Information Technologies (ISCIT 2012), Australia, 2012, pp. 1065-1068.
[14] Tran Duc Tan, Automated Regularization Parameter Selection in Born Iterative Method for Ultrasound Tomography, Vietnam Conference on Control and Automation (VCCA-2011), ISBN 978-604-911-020-7, 2011, pp.786-791.
[15] Tran Duc Tan, Gian Quoc Anh, Improvement of Distorted Born Iterative Method for Reconstructing of Sound Speed, Vietnam Conference on Control and Automation (VCCA-2011), ISBN 978-604-911-020-7, 2011, pp.798-803.
[16] Zhou Wang, Student Member : A Universal Image Quality Index, IEEE Signal Processing Letters, Vol. 9, No. 3, March 2002
Để tính toán được sai số của phép khôi phục, cũng như để so sánh ảnh tạo được với vật thể ta cần có một giá trị tham chiếu của thí nghiệp gọi là hàm mục tiêu lý tưởng. Hàm này có thể được tạo bởi phương trình (2.1) như vậy ta có thể tạo hàm mục tiêu lý tưởng cho vật thể hình trụ với tần số đo f0 như sau:
Hàm mục tiêu lý tưởng
SC=zeros(N,N);
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>3.5355*N/10
SC(m,n)=0;
else
SC(m,n)=(2*pi*f0)^2*(1/(c1^2)-1/(co^2));%da doi f thanh f0
end;
end;
end;
Sau khi đã có hàm mục tiêu lý tưởng ta tạo cấu hình hệ đo với việc bố trí các máy phát và máy thu xung quanh vật thể:
L=N*N;
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 co the thay doi duoc
detector=1:2*N:L;%may thu co the thay doi duoc
plot(KK2(detector),KK1(detector),'s')
hold on
plot(K2(transmiter),K1(transmiter),'r*')
hold on;
mesh(abs(SC))
legend('detector','transmiter','scatter area')
Như vậy ta đã có hàm mục tiêu lý tưởng cần khôi phục và một hệ đo, áp dụng thuật toán 1: Lặp vi phân Born, để viết chương trình trên Matlab ta có các hàm con phải tính sau:
Đầu tiên
Tính , , ,và tương ứng sử dụng (2.6) và (2.7)
Tính tín hiệu của sóng tới viết trên Matlab ta có
pix=[];
k=1:N;
for i=1:N
pix=[pix;k];
end;
PINC=[];
for l=transmiter
pinc=besselj(0,ko*h*sqrt((K1(l)-pix').^2+(K2(l)-pix).^2));
%pinc=besselj(0,4.1e3*h*sqrt((K1(l)-pix').^2+(K2(l)-pix).^2));
PINC=[PINC ; pinc];
end;
save PINC_2D_matrix PINC
Tín hiệu trong thực tế có thể đo được bằng cách lấy hiệu số của tín hiệu tại máy thu khi có đối tượng và khi không có đối tượng. Còn trong mô phỏng thì lại có thể tính bằng phương trình (2.7) sử dụng hàm mục tiêu lý tưởng. Như vậy theo phương trình ta còn phải tính hai ma trận B và C, ma trận B và C tính là ma trận hệ số của hàm Green từ các pixel tới máy thu và hệ số Green giữa các pixel:
calculate_B_matrix_DBIM:
BB=[];
ko_SC=sqrt(ko*ko+abs(SC)); % matrix
ko_SC1=sqrt(ko*ko+abs(SC1)); % matrix
BB_SC=[];BB_SC1=[];
for l=detector
B=-.25*j*h*h*besselj(0,ko*h*sqrt((KK1(l)-pix').^2+(KK2(l)-pix).^2));% no update
B_SC=-.25*j*h*h*besselj(0,ko_SC.*h*sqrt((KK1(l)-pix').^2+(KK2(l)-pix).^2));
B_SC1=-.25*j*h*h*besselj(0,ko_SC1.*h*sqrt((KK1(l)-pix').^2+(KK2(l)-pix).^2));
BB=[BB ; B];
BB_SC=[BB_SC ; B_SC];BB_SC1=[BB_SC1 ; B_SC1];
end;
calculate_C_matrix_DBIM:
CC=[];
CC_SC=[];CC_SC1=[];
for l1=1:N
for l2=1:N
C=-.25*j*h*h*besselj(0,ko*h*sqrt((l1-pix').^2+(l2-pix).^2)); %no update
C_SC =-.25*j*h*h*besselj(0,ko_SC.*h*sqrt((l1-pix').^2+(l2-pix).^2));
C_SC1=-.25*j*h*h*besselj(0,ko_SC1.*h*sqrt((l1-pix').^2+(l2-pix).^2));
CC=[CC ; C];
CC_SC=[CC_SC ; C_SC];CC_SC1=[CC_SC1 ; C_SC1];
end;
end;
Sau khi đã tính những tham số và hàm con trên ta tính được từ giá trị đo được và giá trị tiên đoán và tính RRE tương ứng sử dụng công thức (2.11), Tính giá trị mới sử dụng (2.9). Việc tính sử dụng công thức (2.11) ta phải áp dụng phương pháp NCG như sau:
Áp Dụng NCG Để Tính
function[delta_sound]=test_NCG(Mt,delta_sc_t,ni,RRE,gama)
[n1,n2]=size(Mt);
b=Mt'*delta_sc_t;
x=b;
r=b;
delta_sound=zeros(n2,1);
%delta_sound=zz;
for i=1:ni
q=Mt*x;
%alpha=transpose(r)*r/(transpose(q)*q+gama*transpose(x)*x);
alpha=r'*r/(q'*q+gama*x'*x);
%s=transpose(Mt)*q;
s=Mt'*q;
r_update=r-alpha*(s+gama*x);
%beta=(transpose(r_update)*r_update)/(transpose(r)*r);
beta=(r_update'*r_update)/(r'*r);
delta_sound=delta_sound+alpha*x;
x=r_update+beta*x;
r=r_update;
%e=sum(abs(delta_sc_t))/sum(abs(p_sc_exact_t));
temp=delta_sc_t-Mt*delta_sound;
e=temp'*temp/(delta_sc_t'*delta_sc_t);
%tol=delta_sc_t'*delta_sc_t/(p_sc_exact_t'*p_sc_exact_t)
if e
e
fprintf('convergent at step %d \n',i)
break
end
end;
Như vậy ta có chương chính chính như sau với kết quả đầu ra là hàm khôi phục và các giá trị sai số:
Chương Trình Chính
clear,clc,close all;
f0=2e6;%
f=1e6;% tan so don co the thay doi duoc
N=22
Niter=8 %so lan lap
co=1540; % m/s
c1=co*(1-0.02); % Khoa hoc tre
contrast
landa=co/f; % m/s
lamda=co/f0;
ko=2*pi/landa; % increase when f is increased
D=4*lamda
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));%da doi f thanh f0
%SC(m,n)=dis;
end;
% if dis<3.5355*N/20
% SC(m,n)=2*(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;
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 co the thay doi duoc
detector=1:2*N:L;%may thu co the thay doi duoc
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;
err1=zeros(1,Niter);
MSE1=zeros(1,Niter);
PSNR1=zeros(1,Niter);
gama1=zeros(1,Niter);
Quality=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));
%n=0.05*sqrt(var(p_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(iter)=sum(sum(abs(SC-abs(SC1))))/sum(sum(SC));
MSE1(iter)=(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(iter)=Q;
%het
end; % interation for convergen of sound contrast
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:Niter,err1,'-ro');
xlabel('interation');ylabel('error');
figure
plot(1:Niter,MSE1,'-ro');
xlabel('interation');ylabel('MSE');
Поделитесь с Вашими друзьями: |