k-Wave User Forum » Topic: Simulating medium with two acoustic properties

k-Wave User Forum » Topic: Simulating medium with two acoustic properties http://www.k-wave.org/forum/topic/simulating-medium-with-two-acoustic-properties Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 03:09:21 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulating medium with two acoustic properties" http://www.k-wave.org/forum/topic/simulating-medium-with-two-acoustic-properties#post-6041 Tue, 04 Jul 2017 09:01:31 +0000 Bradley Treeby 6041@http://www.k-wave.org/forum/ Hi Mohalmadi, See my comment above. Brad. Mohalmadi on "Simulating medium with two acoustic properties" http://www.k-wave.org/forum/topic/simulating-medium-with-two-acoustic-properties#post-6039 Mon, 03 Jul 2017 16:56:39 +0000 Mohalmadi 6039@http://www.k-wave.org/forum/ Hi Bradley .. Thank you for your response. What I did is creating a medium which has the following properties: c0, rho0, medium.alpha_coeff. Then I wanted to create another medium within it which has the following properties: scattering_c0, scattering_rho0, and scattering_alpha_coeff. By doing that I was able to change the properties of the cavity independently from the medium. Even changing the value of the smaller medium's (scatterer) attenuation doesn't seem to affect the final result. So basically when a sound signal travels from a medium with high attenuation (Brain) and passes through a medium with a lower attenuation (water) we should see a hyperechoic signal coming from the base of the water medium. From the simulation, I wasn't able to see that. So my question is, how can I change the value of the attenuation? Thank you Bradley Treeby on "Simulating medium with two acoustic properties" http://www.k-wave.org/forum/topic/simulating-medium-with-two-acoustic-properties#post-6022 Thu, 29 Jun 2017 15:03:53 +0000 Bradley Treeby 6022@http://www.k-wave.org/forum/ Hi Moe, The field name for the attenuation is called <code>medium.alpha_coeff</code> not <code>medium.attenuation</code>, so it looks like you're not assigning the heterogeneous absorption value. You might also want to try a stronger attenuation contrast. Brad. Mohalmadi on "Simulating medium with two acoustic properties" http://www.k-wave.org/forum/topic/simulating-medium-with-two-acoustic-properties#post-6017 Wed, 28 Jun 2017 04:04:52 +0000 Mohalmadi 6017@http://www.k-wave.org/forum/ Hello I'm trying to mimic some of the artifact associated with the intraoperative ultrasound system. basically, I created a medium that has homogeneous acoustic properties. and within the medium, I created a cavity that has different acoustic properties. I started by changing the code of one of the examples (b-mode linear transducer), however, I couldn't change the attenuation coefficient of the cavity (Ball). So what I did is changing the matrix of the medium. I was able to change the attenuation value, but unfortunately, the B-mode image doesn't change when I change the value. It seems that it only depends on c0 (sound speed) and rho0 (density) but not attenuation. Can anyone explain to me how to do that? Thank you for your help Moe My code is here if you want to check it: close all clear all clc % 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 %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % DEFINE THE K-WAVE GRID %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % set the size of the perfectly matched layer (PML) pml_x_size = 5; % Depth [grid points] pml_y_size = 3; % Width [grid points] pml_z_size = 3; % Length[grid points] % set total number of grid points not including the PML Nx = 64 - 2*pml_x_size; % [grid points] Ny = 32 - 2*pml_y_size; % [grid points] Nz = 32 - 2*pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML x = 20e-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); % Calculations for plotting depth N_whole_x = Nx + (2*pml_x_size); Whole_x_cartesian = N_whole_x * dx; Start_medium_depth = (Whole_x_cartesian-x)/2; %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % DEFINE THE MEDIUM PARAMETERS %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % define the properties of the propagation medium c0 = 1560; % Speed of sound in medium [m/s] rho0 = 1046; % Denisty of medium [kg/m^3] medium.alpha_coeff = 0.60; % Attenuation coefficent of medium [dB/(MHz^y cm)] medium.alpha_power = 1.35; % Power law absorption exponent medium.BonA = 6; % Parameter of nonlinearity medium.alpha_mode = 'no_dispersion'; % create the time array t_end = (Nx*dx)*2.2/c0; % Total simulation time default number is 2.2 however, % when need more time to capture signal you % nedd to increase it [s] kgrid.t_array = makeTime(kgrid, c0, [], t_end); %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % DEFINE THE INPUT SIGNAL %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % The maximum frequency that can be represented on the spatial grid is the % Nyquist limit of c_min/(2*dx) where dx is the grid spacing. % When using higher frequencies in the simulation, decrease the grid spacing. source_strength = 1e6; % Intensity of source [Pa] tone_burst_freq = 0.4e6; % Signal frequesncy [Hz] tone_burst_cycles = 4; % Number of cycles [#] % create the input signal using toneBurst % <becomes our ultrasound pulse> input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Envelope','Rectangular'); % 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 = 8; % Total number of transducer elements transducer.element_width = 2; % Width of each element [grid points] transducer.element_length = 6; % 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; % Speed of sound [m/s] transducer.focus_distance = Inf; % Focus distance [m] transducer.elevation_focus_distance = Inf; % Focus distance in the elevation plane [m] transducer.steering_angle = 0; % Steering angle [degrees] % apodization transducer.transmit_apodization = 'Rectangular'; % Apodization used on transmit transducer.receive_apodization = 'Rectangular'; % Apodization used on receive % define the transducer elements that are currently active number_active_elements = transducer.number_elements; transducer.active_elements = ones(transducer.number_elements, 1); % 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; %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % DEFINE THE MEDIUM PROPERTIES %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % define a large image size to move across number_scan_lines = 30; Nx_tot = Nx; Ny_tot = Ny + number_scan_lines * transducer.element_width; Nz_tot = Nz; % define a random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.00; background_map = background_map_mean + background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]); % define a random distribution of scatterers for the highly scattering region scattering_map_bin = zeros(Nx_tot, Ny_tot, Nz_tot); % The position of the highly scattering object on the x-axis x_min = round(Nx_tot/4); % Starting point x_center = round(Nx_tot/2); % Middle point x_max = round(Nx_tot/2);% End point % The position of the highly scattering object on the y-axis y_min = round(Ny_tot/4); % Starting point y_center = round(Ny_tot/2); % Middle point y_max = round(Ny_tot/2);% End point % The position of the highly scattering object on the z-axis z_min = round(Nz_tot/8); % Starting point z_center = round(Nz_tot/4); % Middle point z_max = round(Nz_tot/4);% End point for k = z_max:Nz_tot for i = x_center:x_max for j = y_center:y_max scattering_map_bin(i,j,k) = 1; end end for i = x_min:x_center for j = y_min:y_center scattering_map_bin(i,j,k) = 1; end end for i = x_center:x_max for j = y_min:y_center scattering_map_bin(i,j,k) = 1; end end for i = x_min:x_center for j = y_center:y_max scattering_map_bin(i,j,k) = 1; end end end transducer_bin = zeros(Nx_tot, Ny_tot, Nz_tot); for i = transducer.position(1); for j = transducer.position(2):transducer.position(2)+transducer_width for k = transducer.position(3):transducer.position(3)+ transducer.element_length transducer_bin(i,j,k) = 1; end end end % Define the properties of the highly scatterer object (Water) scattering_c0 = 1450 *scattering_map_bin; % Speed of sound in medium [m/s] scattering_rho0 = 993 *scattering_map_bin; % Denisty of medium [kg/m^3] scattering_alpha_coeff = 1 *scattering_map_bin; % Attenuation coefficent of medium [dB/(MHz^y cm)] % 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; attenuation_map= medium.alpha_coeff *ones(Nx_tot, Ny_tot, Nz_tot).*background_map; % Assign region sound_speed_map(scattering_map_bin == 1) = scattering_c0(scattering_map_bin == 1); density_map(scattering_map_bin == 1) = scattering_rho0(scattering_map_bin == 1); attenuation_map(scattering_map_bin == 1) = scattering_alpha_coeff(scattering_map_bin == 1); %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % RUN THE SIMULATION %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % preallocate the storage scan_lines = zeros(number_scan_lines, kgrid.Nt); % 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}; % run the simulation if set to true, otherwise, load previous results from % disk if RUN_SIMULATION % set medium position medium_position = 1; for scan_line_index = 1:number_scan_lines % update the command line status disp(''); disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]); % load the current section of the medium medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :); medium.density = density_map(:, medium_position:medium_position + Ny - 1, :); medium.attenuation = attenuation_map(:, medium_position:medium_position + Ny - 1, :); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data); % update medium position medium_position = medium_position + transducer.element_width; end % save the scan lines to disk save example_us_bmode_scan_lines scan_lines; else % load the scan lines from disk load example_us_bmode_scan_lines end %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % PROCESS THE RESULTS %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % 1) 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(end/2, :); % 2) 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.4; % [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, :); % 3) 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); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(3, :) = scan_lines_fund(end/2, :); % 4) Envelope Detection % envelope detection scan_lines_fund = envelopeDetection(scan_lines_fund); scan_lines_harm = envelopeDetection(scan_lines_harm); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(4, :) = scan_lines_fund(end/2, :); % 5) Log Compression % normalised log compression compression_ratio = 3; scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true); scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(5, :) = scan_lines_fund(end/2, :); % 6) Scan Conversion % upsample the image using linear interpolation scale_factor = 2; scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % VISUALISATION %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ % 1) Plotting the transducer: transducer.plot title('Schematic Of The Transducer'); xlabel('x-dimension (voxel)'); ylabel('y-dimension (voxel)'); zlabel('z-dimension (voxel)'); % Plotting the input signal figure; plot (input_signal) title('Input Signal (MHz)'); xlabel('Time'); ylabel('Amplitude'); % Plottung the highly scattering object voxelPlot(single(scattering_map_bin)) title('Schematic Of The Highly Scattering Object'); xlabel('x-dimension (voxel)'); ylabel('y-dimension (voxel)'); zlabel('z-dimension (voxel)'); % Plotting transducer and scattering object voxelPlot(single(scattering_map_bin | transducer_bin)); title('Schematic Of The Highly Scattering Object And The Transducer'); xlabel('x-dimension (voxel)'); ylabel('y-dimension (voxel)'); zlabel('z-dimension (voxel)'); colormap(gray); % Extracting a single scan line from the sensor data using the current beamforming settings scan_line = transducer.scan_line(sensor_data); figure; plot (scan_line) title('Received Signal'); xlabel('Time'); ylabel('Amplitude'); % Plotting the sound speed map figure; imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, sound_speed_map(:, 1 + Ny/2:end - Ny/2, Nz/2)); axis image; colormap(gray); set(gca, 'YLim', [0, 8]); title('The Sound Speed Map'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); colorbar; % Plotting the density map figure; imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, density_map(:, 1 + Ny/2:end - Ny/2, Nz/2)); axis image; colormap(gray); set(gca, 'YLim', [0, 8]); title('The Density Map'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); colorbar; % Plotting the attenuation map figure; imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, attenuation_map(:, 1 + Ny/2:end - Ny/2, Nz/2)); axis image; colormap(gray); set(gca, 'YLim', [0, 8]); title('The Attenuation Map'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); colorbar; % Plotting the processed b-mode ultrasound image figure; horz_axis = (0:length(scan_lines_fund(:, 1))-1)*transducer.element_width*dy/scale_factor*1e3; imagesc(horz_axis, r*1e3, scan_lines_fund.'); axis image; colormap(gray); set(gca, 'YLim', [0, 8]); title('The B-mode Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); colorbar; % plot the processed harmonic ultrasound image figure; imagesc(horz_axis, r*1e3, scan_lines_harm.'); axis image; colormap(gray); set(gca, 'YLim', [0, 8]); title('The Harmonic Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); colorbar;