k-Wave User Forum » Topic: B-mode Ultrasound Imaging Test

k-Wave User Forum » Topic: B-mode Ultrasound Imaging Test http://www.k-wave.org/forum/topic/b-mode-ultrasound-imaging-test Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:39:46 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Yuxin on "B-mode Ultrasound Imaging Test" http://www.k-wave.org/forum/topic/b-mode-ultrasound-imaging-test#post-8519 Tue, 26 Apr 2022 19:12:18 +0000 Yuxin 8519@http://www.k-wave.org/forum/ I think it's because of: 1. the focus depth 'transducer.focus_distance = x*0.5; % focus distance [m]' 2. the attenuation and the time gain compensation 'medium.alpha_power=1.05;' 'tgc_alpha = 0.1; % [dB/(MHz cm)]' Bradley Treeby on "B-mode Ultrasound Imaging Test" http://www.k-wave.org/forum/topic/b-mode-ultrasound-imaging-test#post-6916 Tue, 25 Jun 2019 08:51:37 +0000 Bradley Treeby 6916@http://www.k-wave.org/forum/ I haven't run your code, but a couple of things to check: Do the raw a-lines look as would expect before applying any further image processing? For the dark band at the bottom, have you tried increasing the simulation time? guns2roses on "B-mode Ultrasound Imaging Test" http://www.k-wave.org/forum/topic/b-mode-ultrasound-imaging-test#post-6847 Fri, 10 May 2019 13:25:53 +0000 guns2roses 6847@http://www.k-wave.org/forum/ Hi Bradley, I recently ran a test simulation for B-mode ultrasound imaging. The medium is just a cube with randomly distributed scatters. However, the resulting images look a bit weird. I can see a very bright strap zone above 40 mm depth in both fundamental and harmonic images (especially in the harmonic image, this strap zone is very distinctive). Besides, the signal seems lost at the bottom (20 - 80 mm horizontal; around 90 mm depth) of the two images. So I wonder if there is something wrong with my code. Could you please help me with it? Thanks very much. The images and the Matlab code are pasted below. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <a href="https://imageshack.com/i/pl0JRRZ1j" rel="nofollow">https://imageshack.com/i/pl0JRRZ1j</a> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearvars; % 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=10; % [grid points] % set total number of grid points not including the PML Nx = 1000 - 2*pml_x_size; % [grid points] Ny = 1000 - 2*pml_y_size; % [grid points] Nz = 100 - 2*pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML x = Nx*0.1e-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 = 1480; % [m/s] rho0 = 1000; % [kg/m^3] medium.alpha_power=1.05; % create the time array t_end = (Nx*dx)*2.2/c0; % [s] kgrid.t_array = makeTime(kgrid, c0, [], t_end); % ========================================================================= % DEFINE THE INPUT SIGNAL % ========================================================================= % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 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; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer transducer.number_elements = 128; % total number of transducer elements transducer.element_width = 2; % width of each element [grid points] transducer.element_length = 78; % 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 = x*0.5; % focus distance [m] transducer.elevation_focus_distance = x*0.5;% focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] transducer.steering_angle_max = 35; % 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 = makeTransducer(kgrid, transducer); % print out transducer properties transducer.properties; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across steering_angles = -32:1:32; number_scan_lines = length(steering_angles); np2db=0.2*log10(exp(1)); ppts=[1561 1081 3.917*np2db 1.2598 6.756]; % heart muscle medium.alpha_coeff=ppts(3); medium.BonA=ppts(5); reconsf=zeros(Nx, Ny, Nz, 2); reconsf(:,:,:,1)=ppts(1); reconsf(:,:,:,2)=ppts(2); muc=ppts(1); sigmac=15.4673; scattering_c0=normrnd(muc,sigmac,[Nx, Ny, Nz, 5]); scattering_c0(scattering_c0 > 1572) = 1572; scattering_c0(scattering_c0 < 1529) = 1529; mup=ppts(2); sigmap=36.09; scattering_rho0=normrnd(mup,sigmap,[Nx, Ny, Nz, 5]); scattering_rho0(scattering_rho0 > 1143) = 1143; scattering_rho0(scattering_rho0 < 1059) = 1059; reconsf(reconsf == ppts(1)) = scattering_c0(reconsf == ppts(1)); reconsf(reconsf == ppts(2)) = scattering_rho0(reconsf == ppts(2)); % load the current section of the medium medium.sound_speed=reconsf(:, :, :,1); medium.density=reconsf(:, :, :,2); clearvars scattering_rho0 scattering_c0 reconsf % ========================================================================= % 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}; %sensor.record = {'p_max_all'}; %sensor.mask=[]; % run the simulation if set to true, otherwise, load previous results from % disk if RUN_SIMULATION for angle_index = 1:number_scan_lines %: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 = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer, input_args{:}); sensor_data = transducer.combine_sensor_data(sensor_data); % 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 us_ph_bmode_scan_mc_scatter_only scan_lines; else % load the scan lines from disk load us_ph_bmode_scan_mc_scatter_only 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] % create time gain compensation function based on attenuation value, % transmit frequency, and round trip distance tgc_alpha = 0.1; % [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); % ----------------------------- % 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 = 8; scan_lines_fund_cp = logCompression(scan_lines_fund, compression_ratio, true); scan_lines_harm_cp = 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_cp, steering_angles, image_size, c0, kgrid.dt); b_mode_harm = scanConversion(scan_lines_harm_cp, 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] 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);