k-Wave

1. Mohalmadi
   

   Member
   Joined: Sep '16
   Posts: 2

   Hello

   I'm trying to mimic some of the artifact associated with the intraoperative ultrasound system. basically, I created a medium that has homogeneous acoustic properties. and within the medium, I created a cavity that has different acoustic properties.
   I started by changing the code of one of the examples (b-mode linear transducer), however, I couldn't change the attenuation coefficient of the cavity (Ball). So what I did is changing the matrix of the medium. I was able to change the attenuation value, but unfortunately, the B-mode image doesn't change when I change the value. It seems that it only depends on c0 (sound speed) and rho0 (density) but not attenuation. Can anyone explain to me how to do that?

   Thank you for your help
   Moe

   My code is here if you want to check it:

   close all
   clear all
   clc

   % 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 = 5; % Depth [grid points]
   pml_y_size = 3; % Width [grid points]
   pml_z_size = 3; % Length[grid points]

   % set total number of grid points not including the PML
   Nx = 64 - 2*pml_x_size; % [grid points]
   Ny = 32 - 2*pml_y_size; % [grid points]
   Nz = 32 - 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 = makeGrid(Nx, dx, Ny, dy, Nz, dz);

   % Calculations for plotting depth
   N_whole_x = Nx + (2*pml_x_size);
   Whole_x_cartesian = N_whole_x * dx;
   Start_medium_depth = (Whole_x_cartesian-x)/2;

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % DEFINE THE MEDIUM PARAMETERS
   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   % define the properties of the propagation medium
   c0 = 1560; % Speed of sound in medium [m/s]
   rho0 = 1046; % Denisty of medium [kg/m^3]
   medium.alpha_coeff = 0.60; % Attenuation coefficent of medium [dB/(MHz^y cm)]
   medium.alpha_power = 1.35; % Power law absorption exponent
   medium.BonA = 6; % Parameter of nonlinearity
   medium.alpha_mode = 'no_dispersion';

   % create the time array
   t_end = (Nx*dx)*2.2/c0; % Total simulation time default number is 2.2 however,
   % when need more time to capture signal you
   % nedd to increase it [s]

   kgrid.t_array = makeTime(kgrid, c0, [], t_end);

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % DEFINE THE INPUT SIGNAL
   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   % The maximum frequency that can be represented on the spatial grid is the
   % Nyquist limit of c_min/(2*dx) where dx is the grid spacing.
   % When using higher frequencies in the simulation, decrease the grid spacing.

   source_strength = 1e6; % Intensity of source [Pa]
   tone_burst_freq = 0.4e6; % Signal frequesncy [Hz]
   tone_burst_cycles = 4; % Number of cycles [#]

   % create the input signal using toneBurst % <becomes our ultrasound pulse>
   input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'Envelope','Rectangular');

   % 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 = 8; % Total number of transducer elements
   transducer.element_width = 2; % Width of each element [grid points]
   transducer.element_length = 6; % 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; % Speed of sound [m/s]
   transducer.focus_distance = Inf; % Focus distance [m]
   transducer.elevation_focus_distance = Inf; % Focus distance in the elevation plane [m]
   transducer.steering_angle = 0; % Steering angle [degrees]

   % apodization
   transducer.transmit_apodization = 'Rectangular'; % Apodization used on transmit
   transducer.receive_apodization = 'Rectangular'; % Apodization used on receive

   % define the transducer elements that are currently active
   number_active_elements = transducer.number_elements;
   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
   number_scan_lines = 30;
   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.00;
   background_map = background_map_mean + background_map_std*randn([Nx_tot, Ny_tot, Nz_tot]);

   % define a random distribution of scatterers for the highly scattering region

   scattering_map_bin = zeros(Nx_tot, Ny_tot, Nz_tot);

   % The position of the highly scattering object on the x-axis
   x_min = round(Nx_tot/4); % Starting point
   x_center = round(Nx_tot/2); % Middle point
   x_max = round(Nx_tot/2);% End point

   % The position of the highly scattering object on the y-axis
   y_min = round(Ny_tot/4); % Starting point
   y_center = round(Ny_tot/2); % Middle point
   y_max = round(Ny_tot/2);% End point

   % The position of the highly scattering object on the z-axis
   z_min = round(Nz_tot/8); % Starting point
   z_center = round(Nz_tot/4); % Middle point
   z_max = round(Nz_tot/4);% End point

   for k = z_max:Nz_tot

   for i = x_center:x_max
   for j = y_center:y_max
   scattering_map_bin(i,j,k) = 1;
   end
   end

   for i = x_min:x_center
   for j = y_min:y_center
   scattering_map_bin(i,j,k) = 1;
   end
   end

   for i = x_center:x_max
   for j = y_min:y_center
   scattering_map_bin(i,j,k) = 1;
   end
   end

   for i = x_min:x_center
   for j = y_center:y_max
   scattering_map_bin(i,j,k) = 1;
   end
   end

   end

   transducer_bin = zeros(Nx_tot, Ny_tot, Nz_tot);
   for i = transducer.position(1);
   for j = transducer.position(2):transducer.position(2)+transducer_width
   for k = transducer.position(3):transducer.position(3)+ transducer.element_length
   transducer_bin(i,j,k) = 1;
   end
   end
   end

   % Define the properties of the highly scatterer object (Water)
   scattering_c0 = 1450 *scattering_map_bin; % Speed of sound in medium [m/s]
   scattering_rho0 = 993 *scattering_map_bin; % Denisty of medium [kg/m^3]
   scattering_alpha_coeff = 1 *scattering_map_bin; % Attenuation coefficent of medium [dB/(MHz^y cm)]

   % 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;
   attenuation_map= medium.alpha_coeff *ones(Nx_tot, Ny_tot, Nz_tot).*background_map;

   % Assign region
   sound_speed_map(scattering_map_bin == 1) = scattering_c0(scattering_map_bin == 1);
   density_map(scattering_map_bin == 1) = scattering_rho0(scattering_map_bin == 1);
   attenuation_map(scattering_map_bin == 1) = scattering_alpha_coeff(scattering_map_bin == 1);

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % 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;

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

   % 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, :);
   medium.attenuation = attenuation_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

   % save the scan lines to disk
   save example_us_bmode_scan_lines scan_lines;
   else
   % load the scan lines from disk
   load example_us_bmode_scan_lines
   end

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % PROCESS THE RESULTS
   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   % 1) 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, :);

   % 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/2 - t0); % [m]

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % create time gain compensation function based on attenuation value,
   % transmit frequency, and round trip distance
   tgc_alpha = 0.4; % [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);

   % store a copy of the middle scan line to illustrate the effects of each
   % processing step
   scan_line_example(2, :) = scan_lines(end/2, :);

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

   % store a copy of the middle scan line to illustrate the effects of each
   % processing step
   scan_line_example(3, :) = scan_lines_fund(end/2, :);

   % 4) Envelope Detection

   % envelope detection
   scan_lines_fund = envelopeDetection(scan_lines_fund);
   scan_lines_harm = envelopeDetection(scan_lines_harm);

   % store a copy of the middle scan line to illustrate the effects of each
   % processing step
   scan_line_example(4, :) = scan_lines_fund(end/2, :);

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

   % store a copy of the middle scan line to illustrate the effects of each
   % processing step
   scan_line_example(5, :) = scan_lines_fund(end/2, :);

   % 6) Scan Conversion

   % upsample the image using linear interpolation
   scale_factor = 2;
   scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
   scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');

   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   % VISUALISATION
   %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   % 1) Plotting the transducer:
   transducer.plot
   title('Schematic Of The Transducer');
   xlabel('x-dimension (voxel)');
   ylabel('y-dimension (voxel)');
   zlabel('z-dimension (voxel)');

   % Plotting the input signal
   figure;
   plot (input_signal)
   title('Input Signal (MHz)');
   xlabel('Time');
   ylabel('Amplitude');

   % Plottung the highly scattering object
   voxelPlot(single(scattering_map_bin))
   title('Schematic Of The Highly Scattering Object');
   xlabel('x-dimension (voxel)');
   ylabel('y-dimension (voxel)');
   zlabel('z-dimension (voxel)');

   % Plotting transducer and scattering object
   voxelPlot(single(scattering_map_bin | transducer_bin));
   title('Schematic Of The Highly Scattering Object And The Transducer');
   xlabel('x-dimension (voxel)');
   ylabel('y-dimension (voxel)');
   zlabel('z-dimension (voxel)');
   colormap(gray);

   % Extracting a single scan line from the sensor data using the current beamforming settings
   scan_line = transducer.scan_line(sensor_data);
   figure;
   plot (scan_line)
   title('Received Signal');
   xlabel('Time');
   ylabel('Amplitude');

   % Plotting the sound speed map
   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);
   set(gca, 'YLim', [0, 8]);
   title('The Sound Speed Map');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   colorbar;

   % Plotting the density map
   figure;
   imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, density_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
   axis image;
   colormap(gray);
   set(gca, 'YLim', [0, 8]);
   title('The Density Map');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   colorbar;

   % Plotting the attenuation map
   figure;
   imagesc((0:number_scan_lines*transducer.element_width-1)*dy*1e3, (0:Nx_tot-1)*dx*1e3, attenuation_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
   axis image;
   colormap(gray);
   set(gca, 'YLim', [0, 8]);
   title('The Attenuation Map');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   colorbar;

   % Plotting 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, 8]);
   title('The B-mode Image');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   colorbar;

   % plot the processed harmonic ultrasound image
   figure;
   imagesc(horz_axis, r*1e3, scan_lines_harm.');
   axis image;
   colormap(gray);
   set(gca, 'YLim', [0, 8]);
   title('The Harmonic Image');
   xlabel('Horizontal Position [mm]');
   ylabel('Depth [mm]');
   colorbar;

   Posted 6 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Moe,

   The field name for the attenuation is called medium.alpha_coeff not medium.attenuation, so it looks like you're not assigning the heterogeneous absorption value. You might also want to try a stronger attenuation contrast.

   Brad.

   Posted 6 years ago #

3. Mohalmadi
   

   Member
   Joined: Sep '16
   Posts: 2

   Hi Bradley .. Thank you for your response.

   What I did is creating a medium which has the following properties: c0, rho0, medium.alpha_coeff. Then I wanted to create another medium within it which has the following properties: scattering_c0, scattering_rho0, and scattering_alpha_coeff. By doing that I was able to change the properties of the cavity independently from the medium. Even changing the value of the smaller medium's (scatterer) attenuation doesn't seem to affect the final result.

   So basically when a sound signal travels from a medium with high attenuation (Brain) and passes through a medium with a lower attenuation (water) we should see a hyperechoic signal coming from the base of the water medium. From the simulation, I wasn't able to see that.

   So my question is, how can I change the value of the attenuation?

   Thank you

   Posted 6 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Mohalmadi,

   See my comment above.

   Brad.

   Posted 6 years ago #

RSS feed for this topic

Reply

You must log in to post.