k-Wave User Forum » Topic: Increasing the number of active elements and grid points Ny

k-Wave User Forum » Topic: Increasing the number of active elements and grid points Ny http://www.k-wave.org/forum/topic/increasing-number-of-active-elements-and-grid-points-ny Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 03:12:26 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php marjanbn on "Increasing the number of active elements and grid points Ny" http://www.k-wave.org/forum/topic/increasing-number-of-active-elements-and-grid-points-ny#post-8226 Fri, 25 Jun 2021 17:21:57 +0000 marjanbn 8226@http://www.k-wave.org/forum/ Hello, I’m trying to simulate a B-mode image using a phantom with an inclusion. When I use 32 active transducer elements with Ny = 364, the simulation works well. Please find below the resultant B-mode image: <a href="https://photos.app.goo.gl/8LYRwA2CR41abxcK8" rel="nofollow">https://photos.app.goo.gl/8LYRwA2CR41abxcK8</a> However, when I change the number of active elements to 64 or 128 and also increase Ny to 652 (for 64 active elements), I get only a dark patch. I cannot find the reason, it seems the wave is not traveled properly. I’m only changing Ny and the number of active elements and not the Nx or travel time. Please find below the resultant B-mode image: <a href="https://photos.app.goo.gl/8LYRwA2CR41abxcK8" rel="nofollow">https://photos.app.goo.gl/8LYRwA2CR41abxcK8</a> I appreciate it if you can help me with this problem and clarify why this happens. Please see below the codes for both images. In the second code, as I mentioned before, only Ny and the number of active elements have been changed. Best, Marjan -------------------------------------------------------------------------------------------------------------------------- Code related to the first image: % Simulating B-mode Ultrasound Images clearvars; clc; % simulation settings DATA_CAST = 'gpuArray-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 = 20; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 1; % [grid points] c_background = 1530; %[m/s] c_target = 1600; %[m/s] c0 = [c_background c_target]; f_max = 10e6; %5e6; % Maximum frequency c0_min = min([c_background c_target]); % Mininmum apeed of sound c0_max = max([c_background c_target]); % Maximum apeed of sound c0_mean = mean([c_background c_target]); % % set desired grid size in the x- and y- direction not including the PML x = 30e-3; % [m] % y = 40e-3; % [m] points_per_wavelength = 5; %4; dx = c0_min/(points_per_wavelength*f_max); dy = dx; % [m] dz = dx; % [m] Nx = round(x/dx) % Ny = round(y/dy) Nx = 1024 - 2 * pml_x_size; % Check for Prime Factors! % Ny = 672 - 2 * pml_y_size; % 1350 - 2 * pml_y_size; % Check for Prime Factors! % set total number of grid points not including the PML % Nx = 512 - 2 * pml_x_size; % [grid points] Ny = 384 - 2 * pml_y_size; % [grid points] Nz = 2; % create the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium rho0= 1000; % [kg/m^3] medium.alpha_coeff = 0.5; %0.73; %0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.05; % medium.BonA = 9.21; %6; % medium.alpha_mode = 'no_dispersion'; % create the time array t_end = (Nx * dx) * 2.2 / c0_min; % [s] % t_end = (Nx * dx) * 4.2 / c0_min; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = f_max; %7.5e6; %1.5e6; % [Hz] tone_burst_cycles = 4; fs = 1/kgrid.dt; %sampling frequency [Hz] % figure() % create the input signal using toneBurst input_signal = toneBurst(fs, tone_burst_freq, tone_burst_cycles); % input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Plot',true); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength ./ (c_background * rho0)) .* input_signal; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer (considering L14-5/38 Linear % transducer) transducer.number_elements = 32; % total number of transducer elements transducer.element_width = 10; % width of each element [grid points] Pitch = 0.3 mm transducer.element_length = 1; % length of each element [grid points] Elevation height = 4 mm 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_mean; % sound speed [m/s] transducer.focus_distance = 10e-3; % focus distance [m] transducer.elevation_focus_distance = 16e-3; %19e-3; % focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] % apodization transducer.transmit_apodization = 'Hanning'; transducer.receive_apodization = 'Rectangular'; % define the transducer elements that are currently active 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 = kWaveTransducer(kgrid, transducer); % print out transducer properties transducer.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across y_FOV = 40e-3; % number_scan_lines = 96; number_scan_lines = round((dy+y_FOV)/(transducer.element_width*dy)); Nx_tot = Nx; Ny_tot = Ny + number_scan_lines * transducer.element_width; Nz_tot = Nz; % define a random distribution of scatterers for the medium and the % properties background_map_mean = 1; % background_map_std_background = 0.0026; background_map_std_background = 0.008; %0.04; background_map = background_map_mean + background_map_std_background * randn([Nx_tot, Ny_tot, Nz_tot]); sound_speed_map = c_background * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; % background_map_std_target = 0.0036; background_map_std_target = 0.001; %0; background_map_target = background_map_mean + background_map_std_target * randn([Nx_tot, Ny_tot, Nz_tot]); sound_speed_target = c_target * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target; density_target = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target; radius_x = 5e-3; radius_y = 5e-3; radius_z = 0.25e-3; r_x = round(radius_x/dx); r_y = round(radius_y/dy); r_z = round(radius_z/dz); x_pos = 15e-3; % [m] y_pos = 20e-3; % [m] c_x = round(x_pos/dx) ; c_y = round(y_pos/dy)+ Ny/2 + 1; c_z = Nz_tot/2; target_region = makeEllipsoid([Nx_tot, Ny_tot, Nz_tot], [c_x, c_y, c_z], [r_x, r_y, r_z]); sound_speed_map(target_region == 1) = sound_speed_target(target_region == 1); density_map(target_region == 1) = density_target(target_region == 1); y_FOV_mm = (Ny_tot-Ny-1)*dy*1000 % ========================================================================= % RUN THE SIMULATION % ========================================================================= % preallocate the storage scan_lines_pr= 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, ... 'Smooth', [true,true,true]}; % run the simulation if set to true, otherwise, load previous results from % disk if RUN_SIMULATION % set medium position medium_position = 1; %%% predeformation medium % loop through the scan lines 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) ' (Predeformation)']); % 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, :); % run the simulation sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines_pre(scan_line_index, :) = transducer.scan_line(sensor_data); save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_2.mat scan_lines_pre % update medium position medium_position = medium_position + transducer.element_width; end % save the scan lines to disk save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_2.mat scan_lines_pre else % load the scan lines from disk load Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus.mat; end % ========================================================================= % PROCESS THE RESULTS % ========================================================================= % ----------------------------- % Remove Input Signal % ----------------------------- % Precompression % 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_pre = bsxfun(@times, scan_line_win, scan_lines_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(1, :) = scan_lines_pre(end/2, :); else scan_line_example_pre(1, :) = scan_lines_pre((end+1)/2, :); end % ----------------------------- % 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_max * ( (1:length(kgrid.t_array)) * kgrid.dt - t0 ) / 2; % [m] % define absorption value and convert to correct units tgc_alpha_db_cm = medium.alpha_coeff * (tone_burst_freq * 1e-6)^medium.alpha_power; tgc_alpha_np_m = tgc_alpha_db_cm / 8.686 * 100; % create time gain compensation function based on attenuation value and % round trip distance tgc = exp(tgc_alpha_np_m * 2 * r); % apply the time gain compensation to each of the scan lines scan_lines_pre = bsxfun(@times, tgc, scan_lines_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(2, :) = scan_lines_pre(end/2, :); else scan_line_example_pre(2, :) = scan_lines_pre((end+1)/2, :); end % ----------------------------- % Frequency Filtering % ----------------------------- % filter the scan lines using both the transmit frequency and the second % harmonic % scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100, true); scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100); % set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]); % title('Predeformation'); % scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30, true); scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30); % set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]); % title('Predeformation'); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(3, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(3, :) = scan_lines_fund_pre((end+1)/2, :); end % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund_pre = envelopeDetection(scan_lines_fund_pre); scan_lines_harm_pre = envelopeDetection(scan_lines_harm_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(4, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(4, :) = scan_lines_fund_pre((end+1)/2, :); end % ----------------------------- % Log Compression % ----------------------------- % normalised log compression compression_ratio = 30; scan_lines_fund_pre = logCompression(scan_lines_fund_pre, compression_ratio, true); scan_lines_harm_pre = logCompression(scan_lines_harm_pre, compression_ratio, true); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(5, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(5, :) = scan_lines_fund_pre((end+1)/2, :); end % ----------------------------- % Scan Conversion % ----------------------------- % upsample the image using linear interpolation scale_factor = 30; scan_lines_fund_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); scan_lines_harm_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); % ========================================================================= % VISUALISATION % ========================================================================= 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)... .* density_map(:, 1 + Ny/2:end - Ny/2, Nz/2)); axis image; colormap(gray); title('Impedance Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); % plot the processing steps figure; stackedPlot(kgrid.t_array * 1e6, {'1. Beamformed Signal', '2. Time Gain Compensation', '3. Frequency Filtering', '4. Envelope Detection', '5. Log Compression'}, scan_line_example_pre); xlabel('Time [\mus]'); title('Predeformation'); set(gca, 'XLim', [5, t_end * 1e6]); % plot the processed b-mode ultrasound image figure; horz_axis = (0:length(scan_lines_fund_pre(:, 1)) - 1) * transducer.element_width * dy / scale_factor * 1e3; imagesc(horz_axis, r * 1e3, scan_lines_fund_pre.'); axis image; colormap(gray); set(gca, 'YLim', [2, 30]); title('B-mode Image (Predeformation)'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ Code related to the second image: % Simulating B-mode Ultrasound Images clearvars; clc; % simulation settings DATA_CAST = 'gpuArray-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 = 20; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 1; % [grid points] c_background = 1530; %[m/s] c_target = 1600; %[m/s] c0 = [c_background c_target]; f_max = 10e6; %5e6; % Maximum frequency c0_min = min([c_background c_target]); % Mininmum apeed of sound c0_max = max([c_background c_target]); % Maximum apeed of sound c0_mean = mean([c_background c_target]); % % set desired grid size in the x- and y- direction not including the PML x = 30e-3; % [m] % y = 40e-3; % [m] points_per_wavelength = 5; %4; dx = c0_min/(points_per_wavelength*f_max); dy = dx; % [m] dz = dx; % [m] Nx = round(x/dx) % Ny = round(y/dy) Nx = 1024 - 2 * pml_x_size; % Check for Prime Factors! Ny = 672 - 2 * pml_y_size; % 1350 - 2 * pml_y_size; % Check for Prime Factors! % set total number of grid points not including the PML % Nx = 512 - 2 * pml_x_size; % [grid points] % Ny = 384 - 2 * pml_y_size; % [grid points] Nz = 2; % create the k-space grid kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium rho0= 1000; % [kg/m^3] medium.alpha_coeff = 0.5; %0.73; %0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.05; % medium.BonA = 9.21; %6; % medium.alpha_mode = 'no_dispersion'; % create the time array t_end = (Nx * dx) * 2.2 / c0_min; % [s] % t_end = (Nx * dx) * 4.2 / c0_min; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = f_max; %7.5e6; %1.5e6; % [Hz] tone_burst_cycles = 4; fs = 1/kgrid.dt; %sampling frequency [Hz] % figure() % create the input signal using toneBurst input_signal = toneBurst(fs, tone_burst_freq, tone_burst_cycles); % input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Plot',true); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength ./ (c_background * rho0)) .* input_signal; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer (considering L14-5/38 Linear % transducer) transducer.number_elements = 64; % total number of transducer elements transducer.element_width = 10; % width of each element [grid points] Pitch = 0.3 mm transducer.element_length = 1; % length of each element [grid points] Elevation height = 4 mm 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_mean; % sound speed [m/s] transducer.focus_distance = 10e-3; % focus distance [m] transducer.elevation_focus_distance = 16e-3; %19e-3; % focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] % apodization transducer.transmit_apodization = 'Hanning'; transducer.receive_apodization = 'Rectangular'; % define the transducer elements that are currently active 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 = kWaveTransducer(kgrid, transducer); % print out transducer properties transducer.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across y_FOV = 40e-3; % number_scan_lines = 96; number_scan_lines = round((dy+y_FOV)/(transducer.element_width*dy)); Nx_tot = Nx; Ny_tot = Ny + number_scan_lines * transducer.element_width; Nz_tot = Nz; % define a random distribution of scatterers for the medium and the % properties background_map_mean = 1; % background_map_std_background = 0.0026; background_map_std_background = 0.008; %0.04; background_map = background_map_mean + background_map_std_background * randn([Nx_tot, Ny_tot, Nz_tot]); sound_speed_map = c_background * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map; % background_map_std_target = 0.0036; background_map_std_target = 0.001; %0; background_map_target = background_map_mean + background_map_std_target * randn([Nx_tot, Ny_tot, Nz_tot]); sound_speed_target = c_target * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target; density_target = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map_target; radius_x = 5e-3; radius_y = 5e-3; radius_z = 0.25e-3; r_x = round(radius_x/dx); r_y = round(radius_y/dy); r_z = round(radius_z/dz); x_pos = 15e-3; % [m] y_pos = 20e-3; % [m] c_x = round(x_pos/dx) ; c_y = round(y_pos/dy)+ Ny/2 + 1; c_z = Nz_tot/2; target_region = makeEllipsoid([Nx_tot, Ny_tot, Nz_tot], [c_x, c_y, c_z], [r_x, r_y, r_z]); sound_speed_map(target_region == 1) = sound_speed_target(target_region == 1); density_map(target_region == 1) = density_target(target_region == 1); y_FOV_mm = (Ny_tot-Ny-1)*dy*1000 % ========================================================================= % RUN THE SIMULATION % ========================================================================= % preallocate the storage scan_lines_pr= 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, ... 'Smooth', [true,true,true]}; % run the simulation if set to true, otherwise, load previous results from % disk if RUN_SIMULATION % set medium position medium_position = 1; % loop through the scan lines 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) ' (Predeformation)']); % 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, :); % run the simulation sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines_pre(scan_line_index, :) = transducer.scan_line(sensor_data); save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_3.mat scan_lines_pre % update medium position medium_position = medium_position + transducer.element_width; end % save the scan lines to disk save Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus_3.mat scan_lines_pre else % load the scan lines from disk load Bmode_scan_lines_CIRSphantom_1Inclusion_10MHz_10mm_focus.mat; end % ========================================================================= % PROCESS THE RESULTS % ========================================================================= % ----------------------------- % Remove Input Signal % ----------------------------- % Precompression % 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_pre = bsxfun(@times, scan_line_win, scan_lines_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(1, :) = scan_lines_pre(end/2, :); else scan_line_example_pre(1, :) = scan_lines_pre((end+1)/2, :); end % ----------------------------- % 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_max * ( (1:length(kgrid.t_array)) * kgrid.dt - t0 ) / 2; % [m] % define absorption value and convert to correct units tgc_alpha_db_cm = medium.alpha_coeff * (tone_burst_freq * 1e-6)^medium.alpha_power; tgc_alpha_np_m = tgc_alpha_db_cm / 8.686 * 100; % create time gain compensation function based on attenuation value and % round trip distance tgc = exp(tgc_alpha_np_m * 2 * r); % apply the time gain compensation to each of the scan lines scan_lines_pre = bsxfun(@times, tgc, scan_lines_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(2, :) = scan_lines_pre(end/2, :); else scan_line_example_pre(2, :) = scan_lines_pre((end+1)/2, :); end % ----------------------------- % Frequency Filtering % ----------------------------- % filter the scan lines using both the transmit frequency and the second % harmonic % scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100, true); scan_lines_fund_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, tone_burst_freq, 100); % set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]); % title('Predeformation'); % scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30, true); scan_lines_harm_pre = gaussianFilter(scan_lines_pre, 1/kgrid.dt, 2 * tone_burst_freq, 30); % set(gca, 'XLim', [0, 20]); %set(gca, 'XLim', [0, 6]); % title('Predeformation'); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(3, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(3, :) = scan_lines_fund_pre((end+1)/2, :); end % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund_pre = envelopeDetection(scan_lines_fund_pre); scan_lines_harm_pre = envelopeDetection(scan_lines_harm_pre); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(4, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(4, :) = scan_lines_fund_pre((end+1)/2, :); end scan_lines_ED = scan_lines_fund_pre; % scale_factor = 2; % scan_lines_ED = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_ED, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); % ----------------------------- % Log Compression % ----------------------------- % normalised log compression compression_ratio = 30; scan_lines_fund_pre = logCompression(scan_lines_fund_pre, compression_ratio, true); scan_lines_harm_pre = logCompression(scan_lines_harm_pre, compression_ratio, true); % store a copy of the middle scan line to illustrate the effects of each % processing step if mod(number_scan_lines,2)==0 scan_line_example_pre(5, :) = scan_lines_fund_pre(end/2, :); else scan_line_example_pre(5, :) = scan_lines_fund_pre((end+1)/2, :); end % ----------------------------- % Scan Conversion % ----------------------------- % upsample the image using linear interpolation scale_factor = 30; scan_lines_fund_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); scan_lines_harm_pre = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm_pre, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).'); % ========================================================================= % VISUALISATION % ========================================================================= 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)... .* density_map(:, 1 + Ny/2:end - Ny/2, Nz/2)); axis image; colormap(gray); title('Impedance Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); % plot the processing steps figure; stackedPlot(kgrid.t_array * 1e6, {'1. Beamformed Signal', '2. Time Gain Compensation', '3. Frequency Filtering', '4. Envelope Detection', '5. Log Compression'}, scan_line_example_pre); xlabel('Time [\mus]'); title('Predeformation'); set(gca, 'XLim', [5, t_end * 1e6]); % plot the processed b-mode ultrasound image figure; horz_axis = (0:length(scan_lines_fund_pre(:, 1)) - 1) * transducer.element_width * dy / scale_factor * 1e3; imagesc(horz_axis, r * 1e3, scan_lines_fund_pre.'); axis image; colormap(gray); set(gca, 'YLim', [2, 30]); title('B-mode Image (Predeformation)'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]');