k-Wave User Forum » Topic: Ultrasound B-Mode Phased Array: Image not detected

k-Wave User Forum » Topic: Ultrasound B-Mode Phased Array: Image not detected http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 03:11:31 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php kpalac on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8225 Wed, 23 Jun 2021 17:18:47 +0000 kpalac 8225@http://www.k-wave.org/forum/ I'm back with the code. Still cannot see the objects. It's quite interesting since i used three different impedances and densities. Can you help me find and fix the issue? % Simulating B-mode Ultrasound Images Example % This example illustrates how k-Wave can be used for the simulation of % B-mode ultrasound images using a phased-array or sector transducer. It % builds on the Simulating B-mode Ultrasound Images Example. % % To allow the simulated scan line data to be processed multiple times with % different settings, the simulated RF data is saved to disk. This can be % reloaded by setting RUN_SIMULATION = false within the example m-file. The % data can also be downloaded from % <a href="http://www.k-wave.org/datasets/example_us_phased_array_scan_lines.mat" rel="nofollow">http://www.k-wave.org/datasets/example_us_phased_array_scan_lines.mat</a> % % author: Bradley Treeby % date: 7th September 2012 % last update: 22nd January 2020 % % This function is part of the k-Wave Toolbox (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % Copyright (C) 2012-2020 Bradley Treeby % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see <http://www.gnu.org/licenses/>. %#ok<*UNRCH clearvars; addpath('/home/kpalac/k-Wave') % simulation settings DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations RUN_SIMULATION = false; % 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 = 15; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 10; % [grid points] % set total number of grid points not including the PML sc = 1; Nx = 256/sc - 2*pml_x_size; % [grid points] Ny = 256/sc - 2*pml_y_size; % [grid points] Nz = 256/sc - 2*pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML x = 50e-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 = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium c0 = 1496.727708; % [m/s] rho0 = 997.0751; % [kg/m^3] medium.alpha_coeff = 0.5; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % create the time array t_end = (Nx * dx) * 2.2 / c0; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 2e6; % [Pa] tone_burst_freq = 2e6 / sc; % [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; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % define the physical properties of the phased array transducer transducer.number_elements = 32 / sc; % total number of transducer elements transducer.element_width = 2; % width of each element [grid points] transducer.element_length = 85 / sc; % length of each element [grid points] transducer.element_spacing = 2; % spacing (kerf width) between the elements [grid points] % 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 = 25e-3; % focus distance [m] transducer.elevation_focus_distance = 25e-3; % focus distance in the elevation plane [m] transducer.steering_angle = 30; % steering angle [degrees] transducer.steering_angle_max = 30; % maximum 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 random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; background_map = background_map_mean + background_map_std * randn([Nx, Ny, Nz]); % define a random distribution of scatterers for the highly scattering % region scattering_map = randn([Nx, Ny, Nz]); scattering_c0 = c0 + 25 + 50 * scattering_map; % scattering_c0(scattering_c0 > 1900) = 1900; scattering_c0(scattering_c0 < 1400) = 1400; scattering_rho0 = scattering_c0 / 1.5; scattering_c1 = 1550 + 10 + 60 * scattering_map; %bone scattering_rho1 = scattering_c1/1; scattering_c2 = 1570 + 10 + 60 * scattering_map; scattering_rho2 = scattering_c2/2; %scattering_c3 = 1573 + 10 + 100 * scattering_map; %scattering_rho3 = scattering_c3/1.6; %scattering_c4 = 1570 + 10 + 50 * scattering_map; %scattering_rho4 = scattering_c4/1.35; % define properties sound_speed_map = c0 * ones(Nx, Ny, Nz) .* background_map; density_map = rho0 * ones(Nx, Ny, Nz) .* background_map; % define a sphere for a highly scattering region radius = 25e-3; x_pos = 25e-3; y_pos = dy * Ny/2; scattering_region1 = makeBall(Nx, Ny, Nz, round(x_pos / dx), round(y_pos / dx), Nz/2, round(radius / dx)); scattering_region2 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(25e-3 / dx), Nz/2, round(7e-3 / dx)); scattering_region3 = makeBall(Nx, Ny, Nz, round(20e-3 / dx), round(15e-3 / dx), Nz/2, round(4e-3 / dx)); %scattering_region4 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(30 / dx), Nz/2, round(4e-3 / dx)); %scattering_region5 = makeBall(Nx, Ny, Nz, round(45e-3 / dx), round(40 / dx), Nz/2, round(2e-3 / dx)); % assign region sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1); density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1); sound_speed_map(scattering_region2 == 1) = scattering_c1(scattering_region2 == 1); density_map(scattering_region2 == 1) = scattering_rho1(scattering_region2 == 1); sound_speed_map(scattering_region3 == 1) = scattering_c2(scattering_region3 == 1); density_map(scattering_region3 == 1) = scattering_rho2(scattering_region3 == 1); %sound_speed_map(scattering_region4 == 1) = scattering_c3(scattering_region4 == 1); %density_map(scattering_region4 == 1) = scattering_rho3(scattering_region4 == 1); %sound_speed_map(scattering_region5 == 1) = scattering_c4(scattering_region5 == 1); %density_map(scattering_region5 == 1) = scattering_rho4(scattering_region5 == 1); % assign to the medium inputs medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % range of steering angles to test steering_angles = -30:10:30; % preallocate the storage number_scan_lines = length(steering_angles); 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 if RUN_SIMULATION % loop through the range of angles to test for angle_index = 1:number_scan_lines % update the command line status disp(''); disp(['Computing scan line ' num2str(angle_index) ' of ' num2str(number_scan_lines)]); % update the current steering angle transducer.steering_angle = steering_angles(angle_index); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(angle_index, :) = transducer.scan_line(sensor_data); end % save the scan lines to disk save example_us_phased_array_scan_lines scan_lines; else % load the scan lines from disk load example_us_phased_array_scan_lines end % trim the delay offset from the scan line data t0_offset = round(length(input_signal) / 2) + (transducer.appended_zeros - transducer.beamforming_delays_offset); scan_lines = scan_lines(:, t0_offset:end); % get the new length of the scan lines Nt = length(scan_lines(1, :)); % ========================================================================= % 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(Nt * 2, 'Tukey', 'Param', 0.05).'; scan_line_win = [zeros(1, t0_offset * 2), scan_line_win(1:end/2 - t0_offset * 2)]; % apply the window to each of the scan lines scan_lines = bsxfun(@times, scan_line_win, scan_lines); % ----------------------------- % Time Gain Compensation % ----------------------------- % create radius variable r = c0 * (1:Nt) * kgrid.dt / 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 = bsxfun(@times, tgc, scan_lines); % ----------------------------- % 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); % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund = envelopeDetection(scan_lines_fund); scan_lines_harm = envelopeDetection(scan_lines_harm); % ----------------------------- % 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); % ----------------------------- % Scan Conversion % ----------------------------- % set the desired size of the image image_size = [Nx * dx, Ny * dy]; % convert the data from polar coordinates to Cartesian coordinates for % display b_mode_fund = scanConversion(scan_lines_fund, steering_angles, image_size, c0, kgrid.dt); b_mode_harm = scanConversion(scan_lines_harm, steering_angles, image_size, c0, kgrid.dt); % ========================================================================= % VISUALISATION % ========================================================================= % create the axis variables x_axis = [0, Nx * dx * 1e3]; % [mm] y_axis = [0, Ny * dy * 1e3]; % [mm] % plot the data before and after scan conversion figure; subplot(1, 3, 1); imagesc(steering_angles, x_axis, scan_lines.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Raw Scan-Line Data'); subplot(1, 3, 2); imagesc(steering_angles, x_axis, scan_lines_fund.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Processed Scan-Line Data'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_fund); axis square; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); colormap(gray); scaleFig(2, 1); % plot the medium and the B-mode images figure; subplot(1, 3, 1); imagesc(y_axis, x_axis, medium.sound_speed(:, :, end/2)); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Scattering Phantom'); subplot(1, 3, 2); imagesc(y_axis, x_axis, b_mode_fund); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_harm); colormap(gray); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Harmonic Image'); scaleFig(2, 1); Bradley Treeby on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8094 Thu, 25 Mar 2021 16:48:31 +0000 Bradley Treeby 8094@http://www.k-wave.org/forum/ I didn't run your simulation to completion, but the following is given on the command line: <code>WARNING: time step may be too large for a stable simulation.</code> If the sensor data is filled with nans, then it sounds like classic instability. You will need to make the time step smaller. kpalac on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8056 Mon, 08 Feb 2021 19:37:03 +0000 kpalac 8056@http://www.k-wave.org/forum/ Sure. I will add code again. "$" will be put at the beginning of each line that was modified: clearvars; % 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 = 15; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 10; % [grid points] % set total number of grid points not including the PML sc = 1; $ Nx = 300/sc - 2*pml_x_size; % [grid points] $ Ny = 300/sc - 2*pml_y_size; % [grid points] $ Nz = 300/sc - 2*pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML $ x = 50e-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 = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium $ c0 = 1496.727708; % [m/s] $ rho0 = 997.0751; % [kg/m^3] $ medium.alpha_coeff = 0.7642; % [dB/(MHz^y cm)] $ medium.alpha_power = 1.5; medium.BonA = 6; % create the time array t_end = (Nx * dx) * 2.2 / c0; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal $ source_strength = 2e6; % [Pa] $ tone_burst_freq = 2e6 / sc; % [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; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % define the physical properties of the phased array transducer $ transducer.number_elements = 32 / sc; % total number of transducer elements $ transducer.element_width = 5; % width of each element [grid points] $ transducer.element_length = 170 / sc; % length of each element [grid points] $ transducer.element_spacing = 3; % spacing (kerf width) between the elements [grid points] % 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 = 25e-3; % focus distance [m] $transducer.elevation_focus_distance = 25e-3; % focus distance in the elevation plane [m] $transducer.steering_angle = 20; % steering angle [degrees] $transducer.steering_angle_max = 20; % maximum 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 random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; background_map = background_map_mean + background_map_std * randn([Nx, Ny, Nz]); % define a random distribution of scatterers for the highly scattering % region scattering_map = randn([Nx, Ny, Nz]); $scattering_c0 = c0 + 25 + 150 * scattering_map; % $scattering_c0(scattering_c0 > 1600) = 1600; $scattering_c0(scattering_c0 < 1400) = 1400; $scattering_rho0 = scattering_c0 / 1.5; $scattering_c1 = 4080 + 10 + 100 * scattering_map; %bone scattering_rho1 = scattering_c1/1.5; %scattering_c2 = 1530 + 10 + 75 * scattering_map; %scattering_rho2 = scattering_c2/1.5; %scattering_c3 = 1573 + 10 + 100 * scattering_map; %scattering_rho3 = scattering_c3/1.6; %scattering_c4 = 1570 + 10 + 50 * scattering_map; %scattering_rho4 = scattering_c4/1.35; % define properties sound_speed_map = c0 * ones(Nx, Ny, Nz) .* background_map; density_map = rho0 * ones(Nx, Ny, Nz) .* background_map; % define a sphere for a highly scattering region $radius = 25e-3; $x_pos = 25e-3; y_pos = dy * Ny/2; scattering_region1 = makeBall(Nx, Ny, Nz, round(x_pos / dx), round(y_pos / dx), Nz/2, round(radius / dx)); $scattering_region2 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(25e-3 / dx), Nz/2, round(2e-3 / dx)); %scattering_region3 = makeBall(Nx, Ny, Nz, round(20e-3 / dx), round(15 / dx), Nz/2, round(3e-3 / dx)); %scattering_region4 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(30 / dx), Nz/2, round(4e-3 / dx)); %scattering_region5 = makeBall(Nx, Ny, Nz, round(45e-3 / dx), round(40 / dx), Nz/2, round(2e-3 / dx)); % assign region sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1); density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1); $sound_speed_map(scattering_region2 == 1) = scattering_c1(scattering_region2 == 1); $density_map(scattering_region2 == 1) = scattering_rho1(scattering_region2 == 1); %sound_speed_map(scattering_region3 == 1) = scattering_c2(scattering_region3 == 1); %density_map(scattering_region3 == 1) = scattering_rho2(scattering_region3 == 1); %sound_speed_map(scattering_region4 == 1) = scattering_c3(scattering_region4 == 1); %density_map(scattering_region4 == 1) = scattering_rho3(scattering_region4 == 1); %sound_speed_map(scattering_region5 == 1) = scattering_c4(scattering_region5 == 1); %density_map(scattering_region5 == 1) = scattering_rho4(scattering_region5 == 1); % assign to the medium inputs medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % range of steering angles to test $steering_angles = -20:10:20; % preallocate the storage number_scan_lines = length(steering_angles); 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 if RUN_SIMULATION % loop through the range of angles to test for angle_index = 1:number_scan_lines % update the command line status disp(''); disp(['Computing scan line ' num2str(angle_index) ' of ' num2str(number_scan_lines)]); % update the current steering angle transducer.steering_angle = steering_angles(angle_index); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(angle_index, :) = transducer.scan_line(sensor_data); end % save the scan lines to disk save example_us_phased_array_scan_lines scan_lines; else % load the scan lines from disk load example_us_phased_array_scan_lines end % trim the delay offset from the scan line data t0_offset = round(length(input_signal) / 2) + (transducer.appended_zeros - transducer.beamforming_delays_offset); scan_lines = scan_lines(:, t0_offset:end); % get the new length of the scan lines Nt = length(scan_lines(1, :)); % ========================================================================= % 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(Nt * 2, 'Tukey', 'Param', 0.05).'; scan_line_win = [zeros(1, t0_offset * 2), scan_line_win(1:end/2 - t0_offset * 2)]; % apply the window to each of the scan lines scan_lines = bsxfun(@times, scan_line_win, scan_lines); % ----------------------------- % Time Gain Compensation % ----------------------------- % create radius variable r = c0 * (1:Nt) * kgrid.dt / 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 = bsxfun(@times, tgc, scan_lines); % ----------------------------- % 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); % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund = envelopeDetection(scan_lines_fund); scan_lines_harm = envelopeDetection(scan_lines_harm); % ----------------------------- % 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); % ----------------------------- % Scan Conversion % ----------------------------- % set the desired size of the image image_size = [Nx * dx, Ny * dy]; % convert the data from polar coordinates to Cartesian coordinates for % display b_mode_fund = scanConversion(scan_lines_fund, steering_angles, image_size, c0, kgrid.dt); b_mode_harm = scanConversion(scan_lines_harm, steering_angles, image_size, c0, kgrid.dt); % ========================================================================= % VISUALISATION % ========================================================================= % create the axis variables x_axis = [0, Nx * dx * 1e3]; % [mm] y_axis = [0, Ny * dy * 1e3]; % [mm] % plot the data before and after scan conversion figure; subplot(1, 3, 1); imagesc(steering_angles, x_axis, scan_lines.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Raw Scan-Line Data'); subplot(1, 3, 2); imagesc(steering_angles, x_axis, scan_lines_fund.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Processed Scan-Line Data'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_fund); axis square; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); colormap(gray); scaleFig(2, 1); % plot the medium and the B-mode images figure; subplot(1, 3, 1); imagesc(y_axis, x_axis, medium.sound_speed(:, :, end/2)); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Scattering Phantom'); subplot(1, 3, 2); imagesc(y_axis, x_axis, b_mode_fund); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_harm); colormap(gray); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Harmonic Image'); scaleFig(2, 1); Bradley Treeby on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8042 Sat, 06 Feb 2021 12:19:13 +0000 Bradley Treeby 8042@http://www.k-wave.org/forum/ Hi kpalac, Can you highlight which lines of code you changed from the example? Brad. kpalac on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8041 Sat, 06 Feb 2021 10:32:16 +0000 kpalac 8041@http://www.k-wave.org/forum/ kpalac on "Ultrasound B-Mode Phased Array: Image not detected" http://www.k-wave.org/forum/topic/ultrasound-b-mode-phased-array-image-not-detected#post-8040 Wed, 03 Feb 2021 20:50:01 +0000 kpalac 8040@http://www.k-wave.org/forum/ Greetings users, I am experiencing some issues with my code. I'll be writing master thesis about b-mode phased ultrasound imaging using k-Wave. Unfortunately, I cannot get a proper image on the output. The output image is almost completely black and it is not the fault of contrast. I ask You for Your help. Here is my code (the array size is locked, I am obliged to keep kerf, width and length of transducer array) Code: clearvars; % 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 = 15; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 10; % [grid points] % set total number of grid points not including the PML sc = 1; Nx = 300/sc - 2*pml_x_size; % [grid points] Ny = 300/sc - 2*pml_y_size; % [grid points] Nz = 300/sc - 2*pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML x = 50e-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 = kWaveGrid(Nx, dx, Ny, dy, Nz, dz); % ========================================================================= % DEFINE THE MEDIUM PARAMETERS % ========================================================================= % define the properties of the propagation medium c0 = 1496.727708; % [m/s] rho0 = 997.0751; % [kg/m^3] medium.alpha_coeff = 0.7642; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % create the time array t_end = (Nx * dx) * 2.2 / c0; % [s] kgrid.makeTime(c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 2e6; % [Pa] tone_burst_freq = 2e6 / sc; % [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; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % define the physical properties of the phased array transducer transducer.number_elements = 32 / sc; % total number of transducer elements transducer.element_width = 5; % width of each element [grid points] transducer.element_length = 170 / sc; % length of each element [grid points] transducer.element_spacing = 3; % spacing (kerf width) between the elements [grid points] % 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 = 25e-3; % focus distance [m] transducer.elevation_focus_distance = 25e-3; % focus distance in the elevation plane [m] transducer.steering_angle = 15; % steering angle [degrees] transducer.steering_angle_max = 25; % maximum 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 random distribution of scatterers for the medium background_map_mean = 1; background_map_std = 0.008; background_map = background_map_mean + background_map_std * randn([Nx, Ny, Nz]); % define a random distribution of scatterers for the highly scattering % region scattering_map = randn([Nx, Ny, Nz]); scattering_c0 = c0 + 25 + 150 * scattering_map; % %scattering_c0(scattering_c0 > 1600) = 1600; %scattering_c0(scattering_c0 < 1400) = 1400; scattering_rho0 = scattering_c0 / 1.5; scattering_c1 = 4080 + 10 + 100 * scattering_map; %bone scattering_rho1 = scattering_c1/1.5; %scattering_c2 = 1530 + 10 + 75 * scattering_map; %scattering_rho2 = scattering_c2/1.5; %scattering_c3 = 1573 + 10 + 100 * scattering_map; %scattering_rho3 = scattering_c3/1.6; %scattering_c4 = 1570 + 10 + 50 * scattering_map; %scattering_rho4 = scattering_c4/1.35; % define properties sound_speed_map = c0 * ones(Nx, Ny, Nz) .* background_map; density_map = rho0 * ones(Nx, Ny, Nz) .* background_map; % define a sphere for a highly scattering region radius = 25e-3; x_pos = 25e-3; y_pos = dy * Ny/2; scattering_region1 = makeBall(Nx, Ny, Nz, round(x_pos / dx), round(y_pos / dx), Nz/2, round(radius / dx)); scattering_region2 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(25e-3 / dx), Nz/2, round(2e-3 / dx)); %scattering_region3 = makeBall(Nx, Ny, Nz, round(20e-3 / dx), round(15 / dx), Nz/2, round(3e-3 / dx)); %scattering_region4 = makeBall(Nx, Ny, Nz, round(25e-3 / dx), round(30 / dx), Nz/2, round(4e-3 / dx)); %scattering_region5 = makeBall(Nx, Ny, Nz, round(45e-3 / dx), round(40 / dx), Nz/2, round(2e-3 / dx)); % assign region sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1); density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1); sound_speed_map(scattering_region2 == 1) = scattering_c1(scattering_region2 == 1); density_map(scattering_region2 == 1) = scattering_rho1(scattering_region2 == 1); %sound_speed_map(scattering_region3 == 1) = scattering_c2(scattering_region3 == 1); %density_map(scattering_region3 == 1) = scattering_rho2(scattering_region3 == 1); %sound_speed_map(scattering_region4 == 1) = scattering_c3(scattering_region4 == 1); %density_map(scattering_region4 == 1) = scattering_rho3(scattering_region4 == 1); %sound_speed_map(scattering_region5 == 1) = scattering_c4(scattering_region5 == 1); %density_map(scattering_region5 == 1) = scattering_rho4(scattering_region5 == 1); % assign to the medium inputs medium.sound_speed = sound_speed_map; medium.density = density_map; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % range of steering angles to test steering_angles = -25:5:25; % preallocate the storage number_scan_lines = length(steering_angles); 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 if RUN_SIMULATION % loop through the range of angles to test for angle_index = 1:number_scan_lines % update the command line status disp(''); disp(['Computing scan line ' num2str(angle_index) ' of ' num2str(number_scan_lines)]); % update the current steering angle transducer.steering_angle = steering_angles(angle_index); % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:}); % extract the scan line from the sensor data scan_lines(angle_index, :) = transducer.scan_line(sensor_data); end % save the scan lines to disk save example_us_phased_array_scan_lines scan_lines; else % load the scan lines from disk load example_us_phased_array_scan_lines end % trim the delay offset from the scan line data t0_offset = round(length(input_signal) / 2) + (transducer.appended_zeros - transducer.beamforming_delays_offset); scan_lines = scan_lines(:, t0_offset:end); % get the new length of the scan lines Nt = length(scan_lines(1, :)); % ========================================================================= % 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(Nt * 2, 'Tukey', 'Param', 0.05).'; scan_line_win = [zeros(1, t0_offset * 2), scan_line_win(1:end/2 - t0_offset * 2)]; % apply the window to each of the scan lines scan_lines = bsxfun(@times, scan_line_win, scan_lines); % ----------------------------- % Time Gain Compensation % ----------------------------- % create radius variable r = c0 * (1:Nt) * kgrid.dt / 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 = bsxfun(@times, tgc, scan_lines); % ----------------------------- % 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); % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund = envelopeDetection(scan_lines_fund); scan_lines_harm = envelopeDetection(scan_lines_harm); % ----------------------------- % 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); % ----------------------------- % Scan Conversion % ----------------------------- % set the desired size of the image image_size = [Nx * dx, Ny * dy]; % convert the data from polar coordinates to Cartesian coordinates for % display b_mode_fund = scanConversion(scan_lines_fund, steering_angles, image_size, c0, kgrid.dt); b_mode_harm = scanConversion(scan_lines_harm, steering_angles, image_size, c0, kgrid.dt); % ========================================================================= % VISUALISATION % ========================================================================= % create the axis variables x_axis = [0, Nx * dx * 1e3]; % [mm] y_axis = [0, Ny * dy * 1e3]; % [mm] % plot the data before and after scan conversion figure; subplot(1, 3, 1); imagesc(steering_angles, x_axis, scan_lines.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Raw Scan-Line Data'); subplot(1, 3, 2); imagesc(steering_angles, x_axis, scan_lines_fund.'); axis square; xlabel('Steering angle [deg]'); ylabel('Depth [mm]'); title('Processed Scan-Line Data'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_fund); axis square; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); colormap(gray); scaleFig(2, 1); % plot the medium and the B-mode images figure; subplot(1, 3, 1); imagesc(y_axis, x_axis, medium.sound_speed(:, :, end/2)); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Scattering Phantom'); subplot(1, 3, 2); imagesc(y_axis, x_axis, b_mode_fund); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('B-Mode Image'); subplot(1, 3, 3); imagesc(y_axis, x_axis, b_mode_harm); colormap(gray); axis image; xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); title('Harmonic Image'); scaleFig(2, 1);