k-Wave

1. prabodh_iitk
   

   Member
   Joined: Nov '13
   Posts: 6

   Dear Sir,

   I have modified the matlab code given in the following chapter for modeling photoacoustic pressure in 3D for spherical excitation.

   B. T. Cox and P. C. Beard, "Modeling Photoacoustic Propagation in Tissue Using k-Space Techniques," in Photoacoustic Imaging and Spectroscopy, L. V. Wang, Ed. London: CRC Press, 2009, pp. 25-34. ISBN-10:1420059912, ISBN-13:978-1420059915.

   Here is my code:
   '
   tic
   clear all
   clear
   clc

   %MATLAB code for k-space model in 3-dimensions acoustically homogenious media

   N=100; n=(-N/2:N/2-1)'; % 3- dim grid
   v=1500; %mm/ms ----- sound velocity in medium
   dx=4e-2; dy=4e-2; dz= 4e-2; %mm ---- size of voxel
   [x,y,z]=meshgrid(n*dx,n*dy,n*dz); % grid construction
   [kx,ky,kz] = meshgrid(n*2*pi/(N*dx),n*2*pi/(N*dy),n*2*pi/(N*dz)); % k-space grid construction
   k = sqrt(kx.^2 + ky.^2+kz.^2); %

   x0= 1; % radius of spherical absorber in mm
   dd= 30; % detector to center of absorber distance
   rd= N/2; cd= N/2; zd= N/2+dd; % position of detector (voxel)

   p1= zeros(N,N,N);
   for i= 1:N %**************************************** SPHERICAL OBJECT
   for j= 1:N
   for l=1:N
   if ((i-N/2)^2+(j-N/2)^2+(l-N/2)^2)<(x0/dx)^2
   p1(i,j,l)=1; % initital pressure in the object
   end
   end
   end
   end

   p0fft=(fftshift(fftn(p1)));% 3-d FFT of initial pressure
   t_samp= 100; % # of time samples
   t_s= 0;% starting time for data acquisiation
   t_e=1.6e-3; % ending time
   dt= (t_e-t_s)/t_samp; % sampling period
   t1= (t_s:dt:t_e-dt);
   for t_t=1:length(t1)
   p = (ifftn(fftshift(p0fft.*cos(v*k*t1(t_t)))));
   P_d (t_t)=p (rd,cd,zd);
   end

   %% rectangular windowing
   % X= P_d';
   % Y= fft (X);
   % Y1= (fftshift(Y));
   % Y2= zeros(t_samp,1);
   % for i=1:t_samp
   % if i<t_samp/2-0.1*t_samp
   % Y2(i,1)= 0;
   % elseif i>t_samp/2+0.1*t_samp+2
   % Y2(i,1)=0;
   % else
   % Y2(i,1)= Y1(i,1);
   % end
   % end
   %
   % Z= ifftshift((Y2));
   % Zr= ifft(Z);

   %% analytical pressure at detector
   P_an= zeros(1,100);
   P_an= (heaviside(x0-abs(dd*dz-v*t1))).*(dd*dz-v*t1)/(2*dz*dd);

   %% Comparision of results
   figure,plot (t1,P_d,'r') % simulated result
   hold on
   % plot(t1,Zr,'k') % simulated result after windowing
   % hold on
   plot(t1,P_an) % analytical result
   hold off

   toc
   '

   But the solution is not matching with the analytical solution.
   Please have a look at the figure:
   http://imageshack.com/a/img203/8402/s61o.jpg

   Kindly let me know where am I doing it wrong.

   Thank you.

   Regards
   Prabodh

   Posted 10 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Prabodh,

   I think the problem is caused by the periodic boundary condition which is implicit when the FFT is used. It causes the waves leaving one side of the domain to appear immediately on the opposite side. One way to avoid it is to use a very large domain, but this quickly becomes impractical in 3D. An alternative is to put an absorbing boundary condition around the edge, which is what we've done in the function kspaceFirstOrder3D. I recommend that you first setup an example in 2D using kspaceFirstOrder2D, as it is easier to visualise and you can check it more easily, and then progress to a 3D example. It can be helpful to use one of the examples that comes with the toolbox as a template to get you started quicker.

   For the 3D case you may want to try the version which uses C++ code in the background as it runs several times faster than the Matlab version. The function is called kspaceFirstOrder3DC and the details can be found in the manual under the documentation tab. Note that the binaries (or source code) need to be downloaded separately from the Toolbox.

   By the way, there is a k-Wave function for creating balls in 3D, makeBall. Also, what you coded up is also included as a k-Wave function kspaceSecondOrder.

   Good luck with your modelling!

   Ben

   Posted 10 years ago #

RSS feed for this topic

Reply

You must log in to post.