k-Wave User Forum » User Favorites: mj

k-Wave User Forum » User Favorites: mj_262 mj_262 Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 09:41:06 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6673 Wed, 05 Dec 2018 01:10:22 +0000 bencox 6673@http://www.k-wave.org/forum/ Hi Mohsen, Try turning the 'smooth' option off, by adding the optional argument <code>...,'Smooth', false)</code> in kspaceFirstOrder3D. Best wishes, Ben bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6672 Tue, 04 Dec 2018 13:34:56 +0000 bencox 6672@http://www.k-wave.org/forum/ Hi maryam.parto, Good question. I put that 4 in there to make the output equal 1, but now I look as it again it's not clear to me what that factor should be in the general case. I will investigate further.... Thanks for the question! Best wishes, Ben mohsen on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6661 Mon, 26 Nov 2018 01:51:04 +0000 mohsen 6661@http://www.k-wave.org/forum/ Hello People to find a matrix model to build a linear system is in my interest. as you said I put the source.p0 = 1 which means an impulse in my pixel. for each pixel I put it and the solution which is impulse response will be assigned in the columns of the interested matrix. but when I multiply this matrix in the initial value pressure field although the behaviour of the plot is the same as the solution of k wave the magnitude is different. I will be thankful if somebody guides me the Best Mohsen maryam.parto on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6573 Tue, 28 Aug 2018 10:17:16 +0000 maryam.parto 6573@http://www.k-wave.org/forum/ Hello dear bencox, My question is how do you scale the delta function for setting the initial pressure gradient? Actually, I can't understand the line "source.dp0dt=delta./4dt". why did you divide it by 4? Generally, I would be grateful if you explain me how can I obtain the initial pressure gradient from a known initial pressure. thanks a lot Bradley Treeby on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5265 Fri, 09 Oct 2015 12:28:11 +0000 Bradley Treeby 5265@http://www.k-wave.org/forum/ Hi Manyou, (1) Because your input signal is both time varying and spatially distributed, the pressure wave radiated from the source will be the summation of the field from each of the source grid points in the simulation. In turn, the impulse response from one source grid point will be the convolution of the band-limited interpolant (an inherent part of the numerical model), and the Green's function. Finally, the signal reflected back to the transducer surface from the single point scatterer will be integrated over the spatial distribution of the source. Hence, the recorded signal will not be the same as the input. Conceptually, the steps that happen inside k-Wave can be considered in the same way as Field II. (2) Because your transducer is focused, additional zeros are added to the input signal such that the same signal can be used for all the source points. This number is printed to the screen when the simulation runs (see Fig. 3.2 in the manual). This will introduce an additional offset in your simulation. You should also consider using grid dimensions for your simulations that have small prime factors - it will make them much faster. There is some more explanation in Section 3.2 of the k-Wave Manual. Hope that helps, Brad. manyouma on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5251 Sun, 20 Sep 2015 07:54:58 +0000 manyouma 5251@http://www.k-wave.org/forum/ % and here is my code <pre><code>clear all; % simulation settings DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations RUN_SIMULATION = true; % set to false to reload previous results instead of running simulation pitch = 0.208/1000; % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= N_element = 64; % set the size of the perfectly matched layer (PML) pml_x_size = 20; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 10; % [grid points] % calculate the spacing between the grid points grid_per_ele = 1; dx = pitch/grid_per_ele; % [m] dy = dx; % [m] dz = dx; % [m] % set desired grid size in the x-direction not including the PML x = 15/1000; % [m] y = pitch*(N_element+1); z = 0.0135; Nx = ceil(x/dx); Ny = ceil(y/dy); Nz = ceil(z/dz); % create the k-space grid kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium c0 = 1540; % [m/s] rho0 = 1000; % [kg/m^3] medium.alpha_coeff = 0; % [dB/(MHz^y cm)] medium.alpha_power = 0; medium.BonA = 0; % % create the time array % t_end = (Nx*dx)*2.2/c0; % [s] % kgrid.t_array = makeTime(kgrid, c0, [], t_end); fs = 40e6; t_end = (Nx*dx)*2.2/c0; % [s] kgrid.t_array = 0:1/fs:t_end; % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 1e6; % [Hz] tone_burst_cycles = 3; % create the input signal using toneBurst input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength./(c0*rho0)).*input_signal; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer transducer.number_elements = N_element; % total number of transducer elements transducer.element_width = grid_per_ele; % width of each element [grid points] transducer.element_length = 24; % length of each element [grid points] transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points] transducer.radius = inf; % radius of curvature of the transducer [m] % calculate the width of the transducer in grid points transducer_width = transducer.number_elements*transducer.element_width ... + (transducer.number_elements - 1)*transducer.element_spacing; % use this to position the transducer in the middle of the computational grid transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]); % properties used to derive the beamforming delays transducer.sound_speed = c0; % sound speed [m/s] transducer.focus_distance = 20e-3; % focus distance [m] transducer.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] % apodization transducer.transmit_apodization = 'Hanning'; % define the transducer elements that are currently active active_array = ones(transducer.number_elements, 1); transducer.active_elements = active_array; % append input signal used to drive the transducer transducer.input_signal = input_signal; % create the transducer using the defined settings transducer = makeTransducer(kgrid, transducer); % print out transducer properties transducer.properties; %============================================================================================== % physical properties of the transducer receiver.number_elements = 64; % total number of transducer elements receiver.element_width = grid_per_ele; % width of each element [grid points] receiver.element_length = 24; % length of each element [grid points] receiver.element_spacing = 0; % spacing (kerf width) between the elements [grid points] receiver.radius = inf; % radius of curvature of the transducer [m] % calculate the width of the transducer in grid points receiver_width = receiver.number_elements*receiver.element_width ... + (receiver.number_elements - 1)*receiver.element_spacing; % use this to position the transducer in the middle of the computational grid receiver.position = round([1, Ny/2 - receiver_width/2, Nz/2 - receiver.element_length/2]); % properties used to derive the beamforming delays receiver.sound_speed = c0; % sound speed [m/s] receiver.focus_distance = 20e-3; % focus distance [m] receiver.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m] receiver.steering_angle = 0; % steering angle [degrees] % define the transducer elements that are currently active number_active_elements = receiver.number_elements; active_array = zeros(receiver.number_elements, 1); active_array(1:number_active_elements) = 1; receiver.active_elements = active_array; % create the transducer using the defined settings receiver = makeTransducer(kgrid, receiver); % print out transducer properties receiver.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across Nx_tot = Nx; Ny_tot = Ny; Nz_tot = Nz; % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.0; background_map = background_map_mean + background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]); % define properties sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map; density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot).*background_map; simstart = 0.005; phantom_positions = [0 0 simstart]; % 0 0 simstart+0.01]; N_scatter = size(phantom_positions,1); for i = 1:N_scatter current_pos = phantom_positions(i,3); x_pos = round(current_pos/dx); density_map(x_pos, round(Ny/2), round(Nz/2)) = density_map(x_pos, round(Ny/2), round(Nz/2)) *10; end medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings input_args = {... 'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ... 'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false}; sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, receiver, input_args{:}); scan_line_win = getWin(kgrid.Nt*2, 'Tukey', 'Param', 0.05).'; scan_line_win = [zeros(1, length(input_signal)*2), scan_line_win(1:end/2 - length(input_signal)*2)]; scan_lines_k = bsxfun(@times, scan_line_win, sensor_data); scan_lines_fund = gaussianFilter(scan_lines_k, 1/kgrid.dt, tone_burst_freq, 100, true); if mod(N_element,2)==0 extreme = (N_element/2-1/2)*pitch; else extreme = (N_element-1)/2*pitch; end aperture_center_array = linspace(-extreme, extreme, N_element)'; t0 = length(input_signal)*kgrid.dt/2; r = c0*( (1:length(kgrid.t_array))*kgrid.dt/2 - t0); close all; imagesc(aperture_center_array,r*1000,scan_lines_fund'); xlabel('Lateral position (m)'); ylabel('Axial position(mm)');</code></pre> manyouma on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5250 Sun, 20 Sep 2015 04:55:16 +0000 manyouma 5250@http://www.k-wave.org/forum/ Hi Brad, Thanks for your response. I tried your method by changing the grid size. I chose a input signal of 2MHz, and a grid size of 0.104mm. I set the absorption and non-linear BonA factor to 0. This is only a point source (in-continuity in density ) at the depth 5mm. <a href="https://www.dropbox.com/s/5p47w6045y6ip7y/k_wave.pdf?dl=0" rel="nofollow">https://www.dropbox.com/s/5p47w6045y6ip7y/k_wave.pdf?dl=0</a> The result converges since using a grid size of 0.208 renders the same result. see Figure (a) and (b). The input signal I used was in Figure (c). One of the RF line in the aperture center is shown in Figure (d). The Field II output using the same excitation and aperture settings is shown in Figure (e) and (f). I have two questions on the results and hope you can point me to where to find the answers: (1) Why is the impulse response in figure (d) have a different shape compared to the input signal? (I know that if this was in Field II, it would depend on the impulse response of the transducer. But what is the equivalent setting in K-wave) (2) I gave an axial offset to the output signal by shifting the data by half the length of the input signal (as shown in the diagnostic US example), however the peak of the output signal is still not at the exact position of the scatter. I was wondering if there are other things I need to consider to align the origin of the phantom and t_array. Thank you very much! Thanks, Manyou Bradley Treeby on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5231 Fri, 04 Sep 2015 17:00:53 +0000 Bradley Treeby 5231@http://www.k-wave.org/forum/ Hi Manyou, Unfortunately I'm not precisely sure off the top of my head. You could try running a simple simulation with a single-point scatterer in the same geometric position relative your source, and gradually increase the grid resolution to check when the results converge (that is, check when the scattered signal that you receive at the source no longer changes as you increase the grid size). Brad. manyouma on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5220 Wed, 26 Aug 2015 19:36:50 +0000 manyouma 5220@http://www.k-wave.org/forum/ Hi Brad, thank you for your answer! I have a further question on this topic. In theory, should the PSF generated by K-wave be similar to other spatial-impulse-response based software (such as Field II)? If so, How small does the scatterer needs to be to generate similar results? Thanks, Manyou Bradley Treeby on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5173 Tue, 28 Jul 2015 12:27:04 +0000 Bradley Treeby 5173@http://www.k-wave.org/forum/ Hi Majid, To compute the point spread function, I would suggest using a homogeneous medium with a single point scatterer in the desired position. The scatterer should be much smaller than the acoustic wavelength (you may need to adjust the grid spacing to achieve this). Brad. mj_262 on "phantom spread function" http://www.k-wave.org/forum/topic/phantom-spread-function#post-5169 Sat, 25 Jul 2015 16:57:45 +0000 mj_262 5169@http://www.k-wave.org/forum/ hi there my thesis is about "improve resolution in us image" for this work first i need to have psf from b-mode image. for this destination i used "Simulating B-mode Ultrasound Images Example" into examples... how can i make psf with this file??? <a href="http://www.k-wave.org/documentation/example_us_bmode_linear_transducer.php" rel="nofollow">http://www.k-wave.org/documentation/example_us_bmode_linear_transducer.php</a> best regards majid bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-342 Thu, 01 Mar 2012 00:56:37 +0000 bencox 342@http://www.k-wave.org/forum/ Hi Chao, Yes, very similar. Both will remove the high frequency components. Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-341 Thu, 01 Mar 2012 00:45:08 +0000 huangchao 341@http://www.k-wave.org/forum/ Hi, Dr. Cox, I have read your paper, and it's clearer to me now about the oscillation. In your paper, the oscillations in the recorded time series were reduced by smoothing the initial pressure distribution (Kronecker delta function) by using windows in spatial frequency domain (k-space), and I am wondering if this is equivalent to smooth the recorded time series by use of windows in temporal frequency domain. Thanks, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-330 Fri, 24 Feb 2012 10:21:37 +0000 bencox 330@http://www.k-wave.org/forum/ Hi Chao, The oscillations come from the fact that what you are seeing is the propagation of the underlying bandlimited function, which (as the code in my previous post showed) is an oscillating function. So you are getting what you should expect to get. It is not that it is inaccurate, it is that what you thought was originally a delta function wasn't (as explained above) it was just that at t=0 you were sampling it mostly at zeros so it looked like a delta function. At other times you might be sampling at different points on the wave form and so you will now see oscillations. Have a look at <a href="http://www.k-wave.org/papers/2011-Treeby-JASA.pdf">this paper</a> for more explanation. As I said before, you can only get a bandlimited version of the Green's function when calculating it discretely. The only way to get closer to a Dirac delta function is to user a finer mesh. If you define your discrete delta properly (so its height grows as the pixel size reduces) then in the limit of an infinitesimally small pixel size it will converge towards the Dirac delta function, and the output will converge towards the Green's function (without any bandlimiting). Hope that's clearer, Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-329 Thu, 23 Feb 2012 22:06:53 +0000 huangchao 329@http://www.k-wave.org/forum/ Hi Dr. Cox, Thank you very much for the detailed explanation. It's very helpful. I have one more question, that is when I set source.p0 to be 1 at a single grid point to calculate the impulse response ('sensor_data' in the code), or the time derivative of the Green's function, whatever you call it, it has a lot of oscillations. Of course, this can be avoid by smoothing the source.p0, but it will compromise the accuracy of the computed impulse response, because the source is not Kronecker δ after smoothing. Any ideas to how to avoid the oscillation and also maintain the impulse response accuracy at the same time? Thanks, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-328 Thu, 23 Feb 2012 12:45:54 +0000 bencox 328@http://www.k-wave.org/forum/ Hi Chao, On point (1): what you think is a perfect delta function, because you only see a single 1 and many zeros, can be thought of as a discrete representation of a continuous signal, i.e. you are only looking at a sampled version of the underlying continuous signal. The underlying continuous signal depends on the basis functions you choose (i.e. how you choose to interpolate between the samples) and for a Fourier basis, the underlying function that fits these samples is a bandlimited, oscillatory function. Run this code and you will see what I mean. <code> % create the computational grid Nx = 64; % number of grid points in the x (row) direction dx = 1; % grid point spacing in the x direction [m] x = (-Nx/2:Nx/2-1)*dx; % define a single source pixel (numerical delta function) delta = zeros(Nx,1); delta(Nx/2+1) = 1; % calculate underlying bandlimited delta function using a much finer mesh dx_fine = 0.001; x_fine = min(x):dx_fine:max(x); underlying_func = (1/Nx)*cot(x_fine*pi/(dx*Nx)).*sin(x_fine*pi/dx); underlying_func(x_fine==0) = 1; % plots figure subplot(3,1,1) stem(x,delta,'o') axis([min(x) max(x) -0.22 1]) title('''discrete'' delta function') subplot(3,1,2) plot(x_fine,underlying_func,'k') axis([min(x) max(x) -0.22 1]) title('underlying bandlimited function') subplot(3,1,3) stem(x,delta,'o') hold on plot(x_fine,underlying_func,'k') axis([min(x) max(x) -0.22 1]) title('''discrete'' delta is sampled version of bandlimited function') </code> (2) Yes. You could think of this as a linear system (as long as you stick to the linear wave equation and not the nonlinear version). You can imagine a matrix equation with a vector of p(x,t) stacked up on the left, a matrix made up by discretising the wave operator, and an initial condition (or a source, depending on how you have written the operator matrix) on the right. You haven't avoided solving the wave equation: this just is a method for solving the wave equation (and probably not a very efficient one). You need to choose some way of discretizing the wave operator to form the system matrix, and that might be using finite differences, finite elements, spectral approaches, or whatever you like! (3) These are equivalent in the sense that p_0 = \Gamma H is true for photoacoustics (when the optical pulse is short enough), and dp/dt = 0 when u = 0 (which can be seen by combining the mass conservation equation and pressure-density relation). Hope that helps, Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-327 Thu, 23 Feb 2012 07:02:39 +0000 huangchao 327@http://www.k-wave.org/forum/ Dear Dr. Cox, Thanks for your reply, and I have some additional questions: 1. I don't understand why Kronecker delta is a band-limited version of the Dirac delta function. In my opinion, the spectrum of Kronecker delta is still infinitely wide. For example, if you plot the amplitude spectrum of the Kronecker delta using MATALB commands 'plot(abs(fft([1, zeros(1,N)])))', the amplitude of all frequency components is 1. 2. I agree that setting source.p0 is 1 at a single grid point will give the time derivative of the Green's function. This leads to another relevant question: Is solving the linear wave equation (no matter the 2nd order equation, or the 3 coupled 1st order equations) is equivalent to compute the output of a linear system where source is the input? If so, for the photoacoustic wave equation: <img src="http://i.imgur.com/C3jyF.jpg" /> we can treat the source H as the input, p as the output, and other parts of the equation (including the time derivative on the right side) as the linear system operator. In this case, because H=A(r)δ(t) is used as the source term for photoacoustic wave equation, just setting source.p0 is 1 at a single grid point (A(r) becomes Kronecker δ) is equivalent to setting the input as an impulse (if Kronecker δ is indeed a good approximation of Dirac δ). Then using k-Wave to solve the equation will give us the impulse response of the linear system, which can be used to compute the output (p) just like any other linear systems if somehow we don't want to solve the wave equation to compute p. 3. I noticed that when solving the wave equation, 2nd order equation or 1st order equations, there are 2 kinds of initial conditions: <img src="http://i.imgur.com/ZVpQ9.jpg" /> <img src="http://i.imgur.com/fvZT5.jpg" /> Could you tell me what's the relation between these 2 initial conditions? Are they equivalent? Thanks again for your help! Best, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-325 Wed, 22 Feb 2012 14:02:37 +0000 bencox 325@http://www.k-wave.org/forum/ Dear Chao, The frequency spectrum of a Dirac delta function is infinitely wide, whereas any numerical model will have some upper limit on the frequencies it can deal with which is set by the spacing of the grid, so the best you can hope for is a band-limited version of the Green's function. (When you input a Kronecker delta you are effectively using a band-limited version of the Dirac delta function). In practice, you just need to ensure the bandwidth is wide enough for your application. Setting the initial condition that <code>source.p0</code> is 1 at a single pixel will not give you the Green's function but the time derivative of the Green's function, as that initial condition is equivalent to using the time derivative of the delta function as a source term. The Green's function is defined as the solution that arises from using a delta function as the source term which you can get by setting the initial pressure gradient. Something like this example below (which easily extends to 2D and 3D) gives the gist of it. <code> % create the computational grid Nx = 512; % number of grid points in the x (row) direction dx = 0.1e-3; % grid point spacing in the x direction [m] kgrid = makeGrid(Nx, dx); k = kgrid.k; % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] % define a time array [kgrid.t_array dt] = makeTime(kgrid,medium.sound_speed); kgrid.t_array = (0:length(kgrid.t_array)/2)*dt; % define a single source pixel delta = zeros(Nx,1); delta(Nx/2+1) = 1; %scale the delta function to have unit area (in time) and assign to the %initial pressure gradient source.dp0dt = delta./(4*dt); % define a single sensor point sensor.mask = zeros(Nx,1); sensor.mask(Nx/4) = 1; % run the simulation output = kspaceSecondOrder(kgrid, medium, source, sensor,'PlotScale',[0 1.5]); %plot the output, which is the bandlimited Green's function. figure plot(kgrid.t_array, output) </code> On your other point, k-Wave does the spatial calculations using Fourier basis functions (it is a type of Fourier pseudospectral method), but there won't be any additional higher frequency information between the pixels in the grid so you could just interpolate linearly and it will be ok. Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-323 Wed, 22 Feb 2012 07:56:48 +0000 huangchao 323@http://www.k-wave.org/forum/ Hi, I want to compute the wave filed from a point source, which is called Green's function or point spread function, by using k-Wave. Theoretically, the point source is described by a Dirac delta function, so I don't know if it's appropriate to use Kronecker delta function to approximate Green's function in k-Wave to compute Green's function. Using Kronecker delta function, I mean to set 1 pixel value of the source.p0 to be 1, and other pixel values to be zeros. Another relevant question is what the basis or expansion function is used in k-Wave. I mean if I have the initial pressure distribution source.p0 as a matrix, but I want to know the pressure between grid points, how to compute that from source.p0? Thanks in advanced! Best, Chao