k-Wave User Forum » Topic: Photoacoustic Equation fwd solution for spherical excitation

k-Wave User Forum » Topic: Photoacoustic Equation fwd solution for spherical excitation http://www.k-wave.org/forum/topic/photoacoustic-equation-fwd-solution-for-spherical-excitation Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:51:10 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "Photoacoustic Equation fwd solution for spherical excitation" http://www.k-wave.org/forum/topic/photoacoustic-equation-fwd-solution-for-spherical-excitation#post-4252 Sat, 28 Dec 2013 12:43:27 +0000 bencox 4252@http://www.k-wave.org/forum/ 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 <code>kspaceFirstOrder3D</code>. I recommend that you first setup an example in 2D using <code>kspaceFirstOrder2D</code>, 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 <code>kspaceFirstOrder3DC</code> 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, <code>makeBall</code>. Also, what you coded up is also included as a k-Wave function <code>kspaceSecondOrder</code>. Good luck with your modelling! Ben prabodh_iitk on "Photoacoustic Equation fwd solution for spherical excitation" http://www.k-wave.org/forum/topic/photoacoustic-equation-fwd-solution-for-spherical-excitation#post-4251 Sat, 28 Dec 2013 08:30:43 +0000 prabodh_iitk 4251@http://www.k-wave.org/forum/ 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: <a href="http://imageshack.com/a/img203/8402/s61o.jpg" rel="nofollow">http://imageshack.com/a/img203/8402/s61o.jpg</a> Kindly let me know where am I doing it wrong. Thank you. Regards Prabodh