k-Wave

1. TristanWang
   

   Member
   Joined: Dec '21
   Posts: 3

   I always get many NaN values in sensor_data in the linear array B mode simulation when I put two scatterers closed to each other as shown in the below figure. My simulation medium is very simple: a homogeneous medium with sound of speed as 1540 m/s and density as 1000 kg/m3 with two scatterers simulating air bubbles with sound of speed as 340 m/s and density as 1 kg/m3. When there is only one scatterer, everything is fine. However, many NaN values in the sensor data appearred when another scatterer was added to the adjacent position to the initial scatterer. If the distance between two scatterers is not zero, everything is fine, too...
   Could anyone please help to solve the problem?

   Below is my code:

   clear
   
   % 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 = 10;                % [grid points]
   pml_y_size = 5;                % [grid points]
   pml_z_size = 5;                % [grid points]
   
   % set total number of grid points not including the PML
   Nx = 128 - 2 * pml_x_size;      % [grid points]
   Ny = 64 - 2 * pml_y_size;      % [grid points]
   Nz = 64 - 2 * pml_z_size;      % [grid points]
   
   % set desired grid size in the x-direction not including the PML
   x = 20e-3;                      % [m]
   
   % calculate the spacing between the grid points
   dx = x / Nx;                    % [m]
   dy = dx;                        % [m]
   dz = dx;                        % [m]
   
   % create the k-space grid
   kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
   
   % =========================================================================
   % DEFINE THE MEDIUM PARAMETERS
   % =========================================================================
   
   % define the properties of the propagation medium
   c0 = 1540;                  % [m/s]
   rho0 = 1000;                % [kg/m^3]
   medium.alpha_coeff = 0; 	% [dB/(MHz^y cm)]
   medium.alpha_power = 0;
   medium.BonA = 0;
   
   % % create the time array
   t_end = (Nx*dx)*2.5/c0;     % [s]
   kgrid.makeTime(c0, [], t_end);
   
   % =========================================================================
   % DEFINE THE INPUT SIGNAL
   % =========================================================================
   
   % define properties of the input signal
   source_strength = 1e6;          % [Pa]
   tone_burst_freq = 1.5e6;        % [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 = 16;  	% total number of transducer elements
   transducer.element_width = 1;       % width of each element [grid points]
   transducer.element_length = 1;  	% 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 = 10e-3;          % focus distance [m]
   transducer.elevation_focus_distance = 10e-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
   number_active_elements = 16;
   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
   number_scan_lines = 16;
   Nx_tot = Nx;
   Ny_tot = Ny + number_scan_lines * transducer.element_width;
   Nz_tot = Nz;
   
   % define a random distribution of scatterers for the medium
   background_map_mean = 1;
   background_map_std = 0;
   background_map = background_map_mean + background_map_std * randn([Nx_tot, Ny_tot, Nz_tot]);
   
   % define properties
   sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
   density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
   
   simstart =  0.01;
   phantom_positions = [0 0 simstart];
   N_scatter = size(phantom_positions,1);
   
   i = 1;
   current_pos = phantom_positions(i,3);
   x_pos = round(current_pos/dx);
   
   sound_speed_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) =  sound_speed_map(x_pos, round(Ny_tot/2), round(Nz_tot/2))*340/1540;
   sound_speed_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) =  sound_speed_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2))*340/1540;
   
   density_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) =  density_map(x_pos, round(Ny_tot/2), round(Nz_tot/2)) *0.001;
   density_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) =  density_map(x_pos+1, round(Ny_tot/2), round(Nz_tot/2)) *0.001;
   
   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================
   
   % preallocate the storage
   scan_lines = zeros(number_scan_lines, kgrid.Nt);
   
   % set the input settings
   input_args = {...
       'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
       'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
   
   % run the simulation if set to true, otherwise, load previous results from
   % disk
   if RUN_SIMULATION
   
       % set medium position
       medium_position = 1;
   
       % loop through the scan lines
   %     for scan_line_index = 1:number_scan_lines
       for scan_line_index = 1
   
           % update the command line status
           disp('');
           disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
   
           % load the current section of the medium
           medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :);
           medium.density = density_map(:, medium_position:medium_position + Ny - 1, :);
   
           % run the simulation
           sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
   
           % extract the scan line from the sensor data
           scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
   
           % update medium position
           medium_position = medium_position + transducer.element_width;
   
       end
   
   else
   
       % load the scan lines from disk
       load example_us_bmode_scan_lines
   
   end
   
   % =========================================================================
   % 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(kgrid.Nt * 2, 'Tukey', 'Param', 0.05).';
   scan_line_win = [zeros(1, length(input_signal) * 2), scan_line_win(1:end/2 - length(input_signal) * 2)];
   
   % apply the window to each of the scan lines
   scan_lines = bsxfun(@times, scan_line_win, scan_lines);
   
   % store a copy of the middle scan line to illustrate the effects of each
   % processing step
   % scan_line_example(1, :) = scan_lines(end/2, :);
   
   % -----------------------------
   % Time Gain Compensation
   % -----------------------------
   
   % create radius variable assuming that t0 corresponds to the middle of the
   % input signal
   t0 = length(input_signal) * kgrid.dt / 2;
   r = c0 * ( (1:length(kgrid.t_array)) * kgrid.dt - 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 = bsxfun(@times, tgc, scan_lines);
   
   % -----------------------------
   % Envelope Detection
   % -----------------------------
   
   % envelope detection
   scan_lines_fund = envelopeDetection(scan_lines);
   
   % -----------------------------
   % Log Compression
   % -----------------------------
   
   % normalised log compression
   compression_ratio = 3;
   scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
   
   % -----------------------------
   % Scan Conversion
   % -----------------------------
   
   % upsample the image using linear interpolation
   scale_factor = 1;
   
   % =========================================================================
   % VISUALISATION
   % =========================================================================
   
   % plot the medium, truncated to the field of view
   figure;
   imagesc((0:number_scan_lines * transducer.element_width - 1) * dy * 1e3, (0:Nx_tot-1) * dx * 1e3, sound_speed_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
   axis image;
   colormap(gray);
   colorbar;
   set(gca, 'YLim', [0, 20]);
   title('Scattering Phantom');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   
   % plot the processed b-mode ultrasound image
   figure;
   horz_axis = (0:length(scan_lines_fund(:, 1)) - 1) * transducer.element_width * dy / scale_factor * 1e3;
   imagesc(horz_axis, r * 1e3, scan_lines_fund.');
   axis image;
   colormap(gray);
   set(gca, 'YLim', [0, 20]);
   title('B-mode Image');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');

   Posted 2 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi TristanWang,

   Is the simulation stable and converged? Does the simulation work if you make the density and sound speed in the bubble unrealistically large, ie. multiply by a factor of 10 or 100? (You could potentially do this as a trick to make the simulations more tractable without affecting the reflections very much.)

   Best wishes
   Ben

   Posted 2 years ago #

RSS feed for this topic

Reply

You must log in to post.