k-Wave

1. guns2roses
   

   Member
   Joined: Jul '18
   Posts: 3

   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.

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   https://imageshack.com/i/pl0JRRZ1j

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   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);

   Posted 5 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   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?

   Posted 5 years ago #

3. Yuxin
   

   Member
   Joined: Apr '21
   Posts: 2

   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)]'

   Posted 2 years ago #

RSS feed for this topic

Reply

You must log in to post.