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???
http://www.k-wave.org/documentation/example_us_bmode_linear_transducer.php
best regards
majid
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
phantom spread function
(7 posts) (3 voices)-
Posted 8 years ago #
-
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.
Posted 8 years ago # -
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,
ManyouPosted 8 years ago # -
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.
Posted 8 years ago # -
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.https://www.dropbox.com/s/5p47w6045y6ip7y/k_wave.pdf?dl=0
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,
ManyouPosted 8 years ago # -
% and here is my 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)');Posted 8 years ago # -
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.
Posted 8 years ago #
Reply
You must log in to post.