k-Wave User Forum » Topic: Row of pins phantom

k-Wave User Forum » Topic: Row of pins phantom http://www.k-wave.org/forum/topic/row-of-pins-phantom Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:11:16 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Row of pins phantom" http://www.k-wave.org/forum/topic/row-of-pins-phantom#post-3653 Sat, 18 May 2013 09:06:55 +0000 Bradley Treeby 3653@http://www.k-wave.org/forum/ Hi Avinoam, I suspect the problem is related to spectral blocking. In your simulation, the maximum frequency supported by the grid is 3.85 MHz and your source frequency is set to 1.2 MHz. However, because you are using a tone-burst, the actual frequency content of your input signal varies over some range. You can check this by looking at the frequency spectrum of your input signal, e.g., <code>spect(input_signal, 1/kgrid.dt, 'Plot', [true, false]);</code> When you run a nonlinear simulation, higher frequency harmonics of your input signal will be generated due to the steepening waves, however, some of these frequencies are not supported by the grid, so are incorrectly aliased to lower frequencies (this behaviour is called spectral blocking). This is likely to be the cause of the high-frequency noise that you have observed. You could try using a lower source frequency, reducing the source amplitude, using a smaller grid spacing, etc. Brad. barz on "Row of pins phantom" http://www.k-wave.org/forum/topic/row-of-pins-phantom#post-3643 Thu, 16 May 2013 03:43:33 +0000 barz 3643@http://www.k-wave.org/forum/ Hi Brad, I've tried both your suggestions without success. When the non linearity is switched off the high frequency noise decreases but did not disappear. In addition I am trying to simulate harmonic imaging and therefore this is not a true solution. I am attaching my code hoping it could help. Thanks for your help! Avinoam clear all; % simulation settings data_cast = 'gpuArray-single'; run_simulation = true; % ========================================================================= % DEFINE THE K-WAVE GRID % ========================================================================= % 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] % set total number of grid points not including the PML Nx = 512 - 2*PML_X_SIZE; % [grid points] Ny = 128 - 2*PML_Y_SIZE; % [grid points] Nz = 64 - 2*PML_Z_SIZE; % [grid points] % set desired grid size in the x-direction not including the PML x = 94.4e-3; % [m] % calculate the spacing between the grid points dx = x/Nx; % [m] dy = dx; % [m] dz = dx; % [m] % 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] scattering_c0 = 1613; % [m/s] scattering_rho0 = 1120; % [kg/m^3] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across Nx_tot = Nx; Ny_tot = Ny; Nz_tot = Nz; % % create a binary sensor mask with seven detection positions pointMask = zeros(Nx_tot, Ny_tot, Nz_tot); pointMask(round([Nx/8, Nx/4, 3*Nx/8, Nx/2, 5*Nx/8, 6*Nx/8, 7*Nx/8]), Ny/2, Nz/2) = 1; % define background properties sound_speed_map = c0*ones(Nx_tot, Ny_tot, Nz_tot); density_map = rho0*ones(Nx_tot, Ny_tot, Nz_tot); % assign pins sound_speed_map(pointMask == 1) = scattering_c0; density_map(pointMask == 1) = scattering_rho0; % load the current section of the medium medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer transducer.number_elements = 32; % total number of transducer elements transducer.element_width = 2; % 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 = 70e-3; % focus distance [m] transducer.elevation_focus_distance = 70e-3;% focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] % apodization transducer.transmit_apodization = 'Hanning'; transducer.receive_apodization = 'Hanning'; % define the transducer elements that are currently active transducer.active_elements = ones(transducer.number_elements, 1); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % create the time array t_end = (Nx*dx)*2.2/c0; % [s] cfl = 0.2; kgrid.t_array = makeTime(kgrid, medium.sound_speed, cfl, t_end); % define properties of the input signal source_strength = 1.5e6; % [Pa] tone_burst_freq = 1.2e6; % [Hz] tone_burst_cycles = 4; % 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; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % preallocate the storage scan_lines = zeros(1, kgrid.Nt); % set the input settings input_args = {... 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ... 'DataCast', data_cast,'Smooth',[true,true,true]}; % run the simulation if set to true, otherwise, load previous results from % disk if run_simulation % 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; if strcmp(data_cast, 'gpuArray-single') % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(1, :) = gather(transducer.scan_line(sensor_data)); else % run the simulation sensor_data = kspaceFirstOrder3DC(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(1, :) = transducer.scan_line(sensor_data); end % save the scan lines to disk save example_us_bmode_scan_lines_pins scan_lines; else % load the scan lines from disk load example_us_bmode_scan_lines_pins end % ========================================================================= % PROCESS THE RESULTS % ========================================================================= % ----------------------------- % Remove Input Signal % ----------------------------- % create a window to set the first part of each scan line to zero to remove % interference from the input signal 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)]; % apply the window to each of the scan lines scan_lines = bsxfun(@times, scan_line_win, scan_lines); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(1, :) = scan_lines; % ----------------------------- % Time Gain Compensation % ----------------------------- % create radius variable assuming that t0 corresponds to the middle of the % input signal t0 = length(input_signal)*kgrid.dt/2; r = c0*( (1:length(kgrid.t_array))*kgrid.dt/2 - t0); % [m] % create time gain compensation function based on attenuation value, % transmit frequency, and round trip distance tgc_alpha = 0.3; % [dB/(MHz cm)] tgc = exp(2*tgc_alpha*tone_burst_freq/1e6*r*100); % apply the time gain compensation to each of the scan lines scan_lines = bsxfun(@times, tgc, scan_lines); % store a copy of the middle scan line to illustrate the effects of each % processing step % scan_line_example(2, :) = scan_lines(end/2, :); scan_line_example(2, :) = scan_lines(1, :); % ----------------------------- % Frequency Filtering % ----------------------------- % filter the scan lines using both the transmit frequency and the second % harmonic scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true); scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2*tone_burst_freq, 30, true); Bradley Treeby on "Row of pins phantom" http://www.k-wave.org/forum/topic/row-of-pins-phantom#post-3640 Wed, 15 May 2013 14:39:38 +0000 Bradley Treeby 3640@http://www.k-wave.org/forum/ Hi Avinoam, There are a couple of things to try. First, if you are using a strong sound-speed contrast, make sure that your time step is set appropriately. In the example you mention, you will need to move the line <code>kgrid.t_array = makeTime(kgrid, c0, [], t_end)</code> to after you have created your medium maps, and replace <code>c0</code> with <code>medium.sound_speed</code>. Second, you could try switching off the nonlinearity (by removing the line <code>medium.BonA = 6;</code>) to see if the high-frequency noise is related to spectral blocking (where the nonlinearity generates high-frequency harmonics not supported by the spatial grid). Brad. barz on "Row of pins phantom" http://www.k-wave.org/forum/topic/row-of-pins-phantom#post-3636 Wed, 15 May 2013 00:53:38 +0000 barz 3636@http://www.k-wave.org/forum/ Hi Brad, I am trying to modify the example "Simulating B-mode Ultrasound Images" to simulate a B-mode image of row of highly echogenic pins inside a tissue. The simulation results I got so far include a high frequency noise with a peak just below the maximal supported frequency. I have tried to use the smoothing option for both the pulse and the speed and density maps. This noise is not apparent when the speed and density maps of the original phantoms are used. Do you have any Idea about the source of this artifact and how I can remove it? Thanks, Avinoam