k-Wave

1. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Hi all,

   I would like to simulate US beam by 'shifting' the medium. I have 32 elements transducer (with 32 active elements). For each new scanline I updated the value of medium position by define local_medium.soundspeed and local_medium.density. My medium is 4 layers phantom (with different soundspeed and density in each layer) above water. I try to copy the idea from the 'Simulating B-mode Ultrasound Images Example'.

   In this example, scan_line_index iterates from 1:96. My question: it aims to model 'pseudo' transducer with 128 elements? Because it shifts 1 element_width in each iteration. Does the number of scanlines affect the output 'sensor_data' from the function 'kspaceFirstOrder3D'? or It only affects this beamforming step?:

   % extract the scan line from the sensor data
   scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);

   In my case, if I want to simulate 'pseudo' transducer with 128 elements, is it valid to make iteration from 0:(97-1)? Because it makes me easier to adjust the zero padding which I did later on. I just need the output from the function 'kspaceFirstOrder3D' (not the output after beamformed by using 'transducer.scan_line(sensor_data)').

   I generate the final output in 3D matrix with dimensions 128x1078x97 = (desired_num_channels, total_time_index, number_scan_lines). So I need to do zero padding in every iteration.

   This is the copy of the script:

   % =========================================================================
   % ULTRASOUND SIMULATION (MEDIUM 'SHIFTING')
   % Linear array transducer (32 element)
   % Transmit: Sound beam
   % Medium : Layer phantom (4 layers) in water
   % =========================================================================

   clear all;

   %simulation setting
   data_cast = 'single';

   % run the simulation
   run_simulation = true;

   % =========================================================================
   % 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 = 256 - 2*PML_X_SIZE; % [grid points]
   Ny = 256 - 2*PML_Y_SIZE; % [grid points]
   Nz = 128 - 2*PML_Z_SIZE; % [grid points]

   % set desired grid size in the x-direction not including the PML
   x = 40e-3; % [m]
   y = 50e-3; % [m]
   z = 60e-3; % [m]

   % calculate the spacing between the grid points
   dx = x/Nx; % [m]
   dy = y/Ny; % [m]
   dz = z/Nz; % [m]

   % create the k-space grid
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

   % =========================================================================
   % DEFINE THE MEDIUM PARAMETERS
   % =========================================================================

   c0 = 1497; % c in water [m/s]
   rho0 = 1000; % density of water [kg/m^3]

   % create the time array
   t_end = 40e-6; % [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 = 5e6; % [Hz]
   tone_burst_cycles = 10;

   % 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 = 32; % total number of transducer elements
   transducer.element_width = 1; % width of each element [grid points]
   transducer.element_length = 8; % 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 = 11e-3; % focus distance [m]
   transducer.elevation_focus_distance = 11e-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 = 32;
   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
   % =========================================================================

   % defining a large image size to move across
   number_scan_lines = 97;

   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.008;

   % define the thicknesses of layer phantom
   d_layer_1 = 5.4e-3; % [m]
   d_layer_2 = 4.9e-3; % [m]
   d_layer_3 = 4.4e-3; % [m]
   d_layer_4 = 6.3e-3; % [m]

   % define the properties of the propagation medium
   % water
   global_medium.sound_speed = c0*ones(Nx, Ny_tot, Nz); % [m/s]
   global_medium.density = rho0*ones(Nx, Ny_tot, Nz); % [kg/m^3]

   % layer_1
   background_map_layer_1 = background_map_mean + background_map_std*randn([round(Nx*(d_layer_1/x)), Ny_tot, Nz]);
   global_medium.sound_speed(1:round(Nx*(d_layer_1/x)), :, :) = 1520.*background_map_layer_1; % [m/s]
   global_medium.density(1:round(Nx*(d_layer_1/x)), :, :) = 1281.*background_map_layer_1; % [kg/m^3]

   % layer_2
   background_map_layer_2 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1218.*background_map_layer_2; % [m/s]
   global_medium.density(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1510.*background_map_layer_2; % [kg/m^3]

   % layer_3
   background_map_layer_3 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1325.*background_map_layer_3; % [m/s]
   global_medium.density(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1643.*background_map_layer_3; % [kg/m^3]

   % layer_4
   background_map_layer_4 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1429.*background_map_layer_4; % [m/s]
   global_medium.density(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1684.*background_map_layer_4; % [kg/m^3]

   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================

   % set the input settings
   input_args = {...
   'PlotLayout', true, 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
   'DataCast', data_cast};

   % run the simulation if set to true
   if run_simulation
   desired_num_channels = 128;
   total_time_index = 1078;

   % set medium position
   medium_position = 1;

   sensor_data_output = zeros(desired_num_channels, total_time_index, number_scan_lines);

   for scan_line_index = 0:number_scan_lines-1

   % update the command line status
   disp('');
   disp(['Computing scan line ' num2str(scan_line_index+1) ' of ' num2str(number_scan_lines)]);

   % load the current section of the medium
   local_medium.sound_speed = global_medium.sound_speed(:, medium_position:medium_position + Ny - 1, :);
   local_medium.density = global_medium.density(:, medium_position:medium_position + Ny - 1, :);

   % run the simulation
   sensor_data = kspaceFirstOrder3D(kgrid, local_medium, transducer, transducer, input_args{:});

   % Zero padding to fill 3D matrix output
   padding_before = zeros(scan_line_index, total_time_index);
   padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index);

   % Populate 3D matrix output with the output of each iteration
   sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after];

   % update medium position
   medium_position = medium_position + transducer.element_width;

   end

   % save the scan lines to disk
   save sensor_data_output sensor_data_output;

   else
   % load the scan lines from disk
   load sensor_data_output
   end
   % =========================================================================

   Thanks

   Makcik

   Posted 11 years ago #

2. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Here I post again the code with better arrangement:

   % =========================================================================
   % ULTRASOUND SIMULATION (MEDIUM 'SHIFTING')
   % Linear array transducer (32 element)
   % Transmit: Sound beam
   % Medium  : Layer phantom (4 layers) in water
   % =========================================================================
   
   clear all;
   
   %simulation setting
   data_cast = 'single';
   
   % run the simulation
   run_simulation = true;
   
   % =========================================================================
   % 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 = 256 - 2*PML_X_SIZE;   % [grid points]
   Ny = 256 - 2*PML_Y_SIZE;   % [grid points]
   Nz = 128 - 2*PML_Z_SIZE;   % [grid points]
   
   % set desired grid size in the x-direction not including the PML
   x = 40e-3;                  % [m]
   y = 50e-3;                  % [m]
   z = 60e-3;                  % [m]
   
   % calculate the spacing between the grid points
   dx = x/Nx;                  % [m]
   dy = y/Ny;                  % [m]
   dz = z/Nz;                  % [m]
   
   % create the k-space grid
   kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
   
   % =========================================================================
   % DEFINE THE MEDIUM PARAMETERS
   % =========================================================================
   
   c0 = 1497;                    % c in water [m/s]
   rho0 = 1000;                  % density of water [kg/m^3]
   
   % create the time array
   t_end = 40e-6;                  % [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 = 5e6;     % [Hz]
   tone_burst_cycles = 10;
   
   % 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 = 32;      % total number of transducer elements
   transducer.element_width = 1;       % width of each element [grid points]
   transducer.element_length = 8;  	% 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 = 11e-3;           % focus distance [m]
   transducer.elevation_focus_distance = 11e-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 = 32;
   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
   % =========================================================================
   
   % defining a large image size to move across
   number_scan_lines = 97;
   
   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.008;
   
   % define the thicknesses of layer phantom
   d_layer_1 = 5.4e-3;          % [m]
   d_layer_2 = 4.9e-3;          % [m]
   d_layer_3 = 4.4e-3;          % [m]
   d_layer_4 = 6.3e-3;          % [m]
   
   % define the properties of the propagation medium
   % water
   global_medium.sound_speed = c0*ones(Nx, Ny_tot, Nz);      % [m/s]
   global_medium.density = rho0*ones(Nx, Ny_tot, Nz);        % [kg/m^3]
   
   % layer_1
   background_map_layer_1 = background_map_mean + background_map_std*randn([round(Nx*(d_layer_1/x)), Ny_tot, Nz]);
   global_medium.sound_speed(1:round(Nx*(d_layer_1/x)), :, :) = 1520.*background_map_layer_1;        % [m/s]
   global_medium.density(1:round(Nx*(d_layer_1/x)), :, :) = 1281.*background_map_layer_1;            % [kg/m^3]
   
   % layer_2
   background_map_layer_2 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1218.*background_map_layer_2;    % [m/s]
   global_medium.density(round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x)), :, :) = 1510.*background_map_layer_2;    % [kg/m^3]
   
   % layer_3
   background_map_layer_3 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1325.*background_map_layer_3;       % [m/s]
   global_medium.density(round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x)), :, :) = 1643.*background_map_layer_3;       % [kg/m^3]
   
   % layer_4
   background_map_layer_4 = background_map_mean + background_map_std*randn([round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), Ny_tot, Nz]);
   global_medium.sound_speed(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1429.*background_map_layer_4; % [m/s]
   global_medium.density(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x)), :, :) = 1684.*background_map_layer_4; % [kg/m^3]
   
   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================
   
   % set the input settings
   input_args = {...
       'PlotLayout', true, 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
       'DataCast', data_cast};
   
   % run the simulation if set to true
   if run_simulation
       desired_num_channels = 128;
       total_time_index = 1078;
   
       % set medium position
       medium_position = 1;
   
       sensor_data_output = zeros(desired_num_channels, total_time_index, number_scan_lines);
   
       for scan_line_index = 0:number_scan_lines-1
   
           % update the command line status
           disp('');
           disp(['Computing scan line ' num2str(scan_line_index+1) ' of ' num2str(number_scan_lines)]);
   
           % load the current section of the medium
           local_medium.sound_speed = global_medium.sound_speed(:, medium_position:medium_position + Ny - 1, :);
           local_medium.density = global_medium.density(:, medium_position:medium_position + Ny - 1, :);
   
           % run the simulation
           sensor_data = kspaceFirstOrder3D(kgrid, local_medium, transducer, transducer, input_args{:});
   
           % Zero padding to fill 3D matrix output
           padding_before = zeros(scan_line_index, total_time_index);
           padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index);
   
           % Populate 3D matrix output with the output of each iteration
           sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after];
   
           % update medium position
           medium_position = medium_position + transducer.element_width;
   
       end
   
      % save the scan lines to disk
      save sensor_data_output sensor_data_output;
   
   else
      % load the scan lines from disk
      load sensor_data_output
   end

   Posted 11 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Makcik,

   For each simulation, it is only the total number of active elements that matters, as this defines where the acoustic waves will be added to the medium. Similarly for the beamforming step - it is only the number of active transducer elements that matters. In the example you reference, 96 separate simulations are performed, mimicking the collection of scan lines from a transducer with a large number of elements. These are only combined in the final scan conversion step.

   I noticed in your simulation you have a different grid spacing in X, Y, and Z, giving rise to different maximum supported frequencies (3.2MHz by 2.8MHz by 1.1MHz). Is there any reason for this? Also, your source frequency always be less than the maximum supported frequency.

   Hope that helps,

   Brad.

   Posted 11 years ago #

4. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Hi Brad,

   Sorry for my late response. I am quite new to this area and this toolbox, so I didn't realize that the spacing grid affects the supported frequencies. How does it work (how to calculate those 3.2MHz by 2.8MHz by 1.1MHz frequencies)?

   And for the source frequency, is it :
   tone_burst_freq = 5e6; % [Hz] which you mean? Is it suppose to be less or more than supported frequencies?

   Sorry, if I repeat these questions:
   In the example I reference: 96 separate simulations are performed, mimicking the collection of scan lines from a transducer with a large number of elements. What is exactly the number of elements of the transducer that you mimic? Is it 128?

   In my case: if I want to simulate 'pseudo' transducer with 128 elements, is it valid to make iteration from 0:(97-1)? And do this step to rearrange the final output:

   % Zero padding to fill 3D matrix output
           padding_before = zeros(scan_line_index, total_time_index);
           padding_after = zeros(desired_num_channels - number_active_elements - scan_line_index, total_time_index);
   
           % Populate 3D matrix output with the output of each iteration
           sensor_data_output(:, :, scan_line_index + 1) = [padding_before; sensor_data; padding_after];

   so my final otput will be 3D matrix with dimensions 128x1078x97 = (desired_num_channels x total_time_index x number_scan_lines).

   Actually my goal is to simulate ultrasound beam with 32 elements transmit aperture and it shifts (overlap 1 element in each iteration), so If I shift it 96 times I get 97 beam until the last element of my transducer (real or 'pseudo' 128 elements transducer). Therefore, in another post:
   http://www.k-wave.org/forum/topic/simulating-us-beam-by-shifting-the-transducer
   I try to do it in different way, Is it also possible? Or do you have suggestion how to do it properly?

   I will do beamforming later on outside this toolbox, so in this case I just use this toolbox to generate (simulate) raw RF ultrasound data.

   Thank you very much.

   Makcik

   Posted 11 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Makcik,

   The maximum supported temporal frequency in each direction can be calculated based on the grid spacing and the sound speed via the relation f_max = min(c_0)/(2*dx). These frequencies are also reported on the command line when the simulation runs. There is some background information that you may find useful in the k-Wave Manual.

   The entire frequency content of your source (including nonlinear harmonics) should be less than the maximum frequency to avoid aliasing errors. In this particular example, the frequency is set by tone_burst_freq as you say, but keep in mind that a tone-burst will contain a range of frequencies (you can see this by looking at the frequency spectrum of the input signal).

   If you are moving the medium, the number of mimicked transducer elements will be N_elements + N_simulations - 1, so in the example you mention, this will be 127. If you want to simulate the entire medium, the alternative way would be to define the transducer with 128 elements, then within each iteration, update which elements are active.

   Hope that helps,

   Brad.

   Posted 11 years ago #

6. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Hi Brad,

   I try to implement your idea to simulate the whole medium by update which elements are active within each iteration. My total number of elements is 128.
   I also define new computational grid and its properties, and the position of my medium.

   Could you take a look at this? whether my assumptions is already valid or not. Here is the copy of my new script.

   Thank you very much

   Makcik

   % =========================================================================
   % ULTRASOUND SIMULATION
   % Linear array transducer (128 element)
   % Transmit: Sound beam (Focusing)
   % Medium  : Layer phantom (4 layers) in water
   % =========================================================================
   
   clear all;
   
   %simulation setting
   data_cast = 'single';
   
   % run the simulation
   % true : run the simulation
   % false   : load previous results from disk
   run_simulation = true;
   
   % =========================================================================
   % DEFINE THE K-WAVE GRID
   % =========================================================================
   
   % set the size of the perfectly matched layer (PML)
   PML_X_SIZE = 10;            % [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 = 256 - 2*PML_X_SIZE;   % [grid points]
   Ny = 256 - 2*PML_Y_SIZE;   % [grid points]
   Nz = 256 - 2*PML_Z_SIZE;   % [grid points]
   
   % set desired grid size in the x-direction not including the PML
   x = 60e-3;                  % [m]
   y = 60e-3;                  % [m]
   z = 60e-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
   % =========================================================================
   
   c0 = 1497;                    % c in water [m/s]
   rho0 = 1000;                  % density of water [kg/m^3]
   medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
   medium.alpha_power = 1.5;
   medium.BonA = 6;            % represent nonlinearity properties of the medium
   
   % create the time array
   t_end = 40e-6;                  % [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 = 10;
   
   % 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 = 1;       % width of each element [grid points]
   transducer.element_length = 8;  	% 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 = 1530;                % sound speed [m/s]
   transducer.focus_distance = 10.5e-3;          % focus distance [m]
   transducer.elevation_focus_distance = 19e-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
   transducer.active_elements = zeros(transducer.number_elements, 1);
   transducer.active_elements(1:32)=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 random distribution of scatterers for the medium
   background_map_mean = 1;
   background_map_std = 0.008;
   
   % define the thicknesses of layer phantom
   d_layer_1 = 5.4e-3;          % [m]
   d_layer_2 = 4.9e-3;          % [m]
   d_layer_3 = 4.4e-3;          % [m]
   d_layer_4 = 6.3e-3;          % [m]
   
   % define the position of each layer in x direction (4 layers: 1 until 4)
   position_x1 = 1:round(Nx*(d_layer_1/x));
   position_x2 = round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x));
   position_x3 = round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x));
   position_x4 = round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x));
   
   % define the width of phantom in y and z direction (in this case it's equal to transducer width)
   width_y = 32.5424e-3; % [m]
   width_z = 32.5424e-3; % [m]
   
   % define position of each layer in y and z direction
   position_y = round(Ny*(((y-width_y)/2)/y))+1:round(Ny*(y-(y-width_y)/2)/y);
   position_z = round(Nz*(((z-width_z)/2)/z))+1:round(Nz*(z-(z-width_z)/2)/z);
   
   % define the dimension of each layer (and its background map)
   width_backgroundmap_layer_1 = [round(Nx*(d_layer_1/x)), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
   width_backgroundmap_layer_2 = [round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
   width_backgroundmap_layer_3 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
   width_backgroundmap_layer_4 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
   
   % % define the properties of the propagation medium
   % % water
   medium.sound_speed = c0*ones(Nx, Ny, Nz);      % [m/s]
   medium.density = rho0*ones(Nx, Ny, Nz);        % [kg/m^3]
   
   % layer_1
   background_map_layer_1 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_1);
   medium.sound_speed(position_x1, position_y, position_z) = 1520.*background_map_layer_1;        % [m/s]
   medium.density(position_x1, position_y, position_z) = 1281.*background_map_layer_1;            % [kg/m^3]
   
   % layer_2
   background_map_layer_2 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_2);
   medium.sound_speed(position_x2, position_y, position_z) = 1218.*background_map_layer_2;    % [m/s]
   medium.density(position_x2, position_y, position_z) = 1510.*background_map_layer_2;    % [kg/m^3]
   
   % layer_3
   background_map_layer_3 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_3);
   medium.sound_speed(position_x3, position_y, position_z) = 1325.*background_map_layer_3;       % [m/s]
   medium.density(position_x3, position_y, position_z) = 1643.*background_map_layer_3;       % [kg/m^3]
   
   % layer_4
   background_map_layer_4 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_4);
   medium.sound_speed(position_x4, position_y, position_z) = 1429.*background_map_layer_4; % [m/s]
   medium.density(position_x4, position_y, position_z) = 1684.*background_map_layer_4; % [kg/m^3]
   
   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================
   
   % set the input settings
   % plot medium map
   input_args = {...
       'PlotLayout', true, 'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
       'DataCast', data_cast};
   
   % run the simulation if set to true, otherwise, load previous results from disk
   if run_simulation
      num_channels = transducer.number_elements;
      total_time_index = 1078;
      total_scan_lines = 97;
   
       RF_data_SB_transd_shift = zeros(num_channels, total_time_index, total_scan_lines);
   
       for scan_line = 1 : total_scan_lines
           fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines);
   
           transducer_active_elements_position = transducer.active_elements(1:32);
   
           RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
   
           % update the position of transducer active elements
           transducer_active_elements_position = transducer_active_elements_position + transducer.element_width;
       end
   
       save RF_data_SB_transd_shift RF_data_SB_transd_shift;
   else
      % load the scan lines from disk
      load RF_data_SB_transd_shift;
   end

   Posted 11 years ago #

7. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Makcik,

   I haven't run your code, but it looks like you haven't updated the active transducer elements with the for loop. You could use something like:

   % define the transducer elements that are currently active
   transducer.active_elements = zeros(transducer.number_elements, 1);
   transducer.active_elements(scan_line:scan_line+31)=1;

   Here transducer is an object of the kWaveTransducer class which contains the properties of the transducer. There are more details in the k-Wave manual if you get stuck.

   Brad.

   Posted 11 years ago #

8. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Hi Brad,

   I implemented your suggestion in my loop:
   

   % =========================================================================
   % RUN THE SIMULATION
   % =========================================================================
   
   % set the input settings
   % plot medium map
   input_args = {...
       'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
       'DataCast', data_cast};
   
   % run the simulation if set to true, otherwise, load previous results from disk
   if run_simulation
      num_channels = transducer.number_elements;
      total_time_index = 1078;
      total_scan_lines = 97;
   
       RF_data_SB_transd_shift = zeros(num_channels, total_time_index); 
   
       for scan_line = 1 : total_scan_lines
           fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines);
   
           % define the transducer elements that are currently active
             transducer.active_elements = zeros(transducer.number_elements, 1);
             transducer.active_elements(scan_line:scan_line+31)=1;
   
           RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
       end
   
       save RF_data_SB_transd_shift RF_data_SB_transd_shift;

   My questions:

   1. In this case, I have only 2D data, which previously I expected to have the information about total number of elements and time index, but, I get the result which only contains 32 number of elements instead of 128. Why? And also where I can find the information about scanline in my data (in this case 97 beams or scanlines)?

   Actually, in this simulation, I try to mimic my laboratory experiment, which produce 3D DAQ raw RF data, which contains the information about total number of elements, total number of sample, and number of scanlines (beams).

   2. In KSpaceFirstOrder3D function, Is it in each iteration, the data which are filled in each transducer element being overwrited or being summed?

   3. Why does it crash every time I call transducer.plot? I run this toolbox in MATLAB 2012a.

   4. Why is the maximum supported frequency is not match with the formula: f_max = min(c_0)/(2*dx). In my case: It's supposed to be 1218/(2*254.2373um) = 2.3954 MHz, but on the command line when the simulation runs it was 2.3086MHz.

   Thank you very much.

   Makcik

   Posted 11 years ago #

9. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Hi Brad,

   Regarding my 1st question above, I also have tried to run the example from the toolbox (example_us_bmode_linear_transducer.m).

   * I changed the total number of elements to 128, and number of active elements was still the same (32).
   * I run the loop, only until scanline number 2.
   * And I get the result sensor_data with dimension 128 x total time step. All of 128 element is not empty, eventhough the scanline was only 2.
   * But why in my case (in my latest code), I get the data with dimension only 32 x total time index, eventhough I have 128 total number of elements and any number of scanlines?

   Thank you very much,

   Makcik

   Posted 11 years ago #

10. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    Hi Makcik,

    To answer your questions:

    1. If you use an object of the kWaveTransducer class (created using makeTransducer) to replace the sensor input, the number of time signals returned will be equal to the number of active elements. With the current code, if you wish to transmit using 32 elements, and receive using a different number of elements, you will need to create a second copy of your transducer, set the desired number of active elements, and use this in place of the sensor input.

    2. If your physical transducer elements are made up from multiple grid points, the signal received by the physical element is calculated at each time point by averaging over the grid points. There is more information about how this is done within the k-Wave manual.

    3. What is the error your receive, is it a memory error? Unfortunately the function voxelPlot does not deal very well with very large 3D matrices.

    4. The number displayed should be correct. Perhaps double check your grid spacing (kgrid.dx, kgrid.dt, kgrid.dz) and your sound speeds (min(medium.sound_speed(:)).

    Brad.

    Posted 11 years ago #

11. maktjik
    

    Member
    Joined: Mar '13
    Posts: 18

    Hi Brad,

    Thanks for your suggestion, and currently my simulation is still running.

    For question number 2:
    I still don't really understand how it works. I read in the manual:

    "The way in which the signals across each element are calculated depends
    on the setting for transducer.elevation_focus_distance. If this is set to Inf, the signals across the grid points within each sensor element are averaged at each time step,and only the average is stored. If a finite elevation focus distance is defined, a buffer the length of the longest beamforming delay is filled using a FIFO queue (first-in, first-out). The elevation beamforming is then computed on the fly once the buffer is filled."

    In my case: I define the elevation focus distance to be finite and each element has 2 grid points. In the first iteration I transmit the signal from the first 32 active elements and receive in 128 active elements, and so on until the last 32 active elements transmit the 97th beam. So from the first iteration, I have first 128 data filled in each transducer element. My question: for the final sensor_data which I name RF_data_SB_transd_shift, is it the average or summing of the 97 set data?

    For question number 3: My MATLAB was simply crash and didn't work at all every time I call transducer.plot, I think it's the issue of memory because I have relatively large 3D computational grid.

    Another question:
    I have another plan to simulate another type of US propagation (plane wave), which:
    - in transmit: each element (total: 128 elements), propagate the signal in the same time, focus_distance is infinite, and elevation_focus_distance is finite.
    - in receive: total 128 elements also receive the echo in the same time.

    Can I simulate it by simply define 128 total elements transducer (with also 128 active elements), define focus_distance to be inf, and set them as KSpaceFirstOrder3D input arguments in the places of source and sensor? Another parameters, such as computational grid, medium parameters, medium properties, input signal and time array are still the same as before.

    Thank you very much,

    Makcik

    Posted 11 years ago #

12. maktjik
    

    Member
    Joined: Mar '13
    Posts: 18

    PS:

    I just get the result from my simulation (transducer shifting, 32 active elements transmit and 128 active elements receive), but I just try it for 1 scanline and 2 scanlines.

    Why the first 2 columns and the last 2 columns of my data are always filled by zero?

    Here is my latest code:
    

    % =========================================================================
    % ULTRASOUND SIMULATION
    % Linear array transducer (128 element)
    % Transmit: Sound beam (Focusing), 32 active elements
    % Receive : 128 active elements
    % Medium  : Layer phantom (4 layers) in water
    % =========================================================================
    
    clear all;
    
    %simulation setting
    data_cast = 'single';
    
    % run the simulation
    % true : run the simulation
    % false   : load previous results from disk
    run_simulation = true;
    
    % =========================================================================
    % DEFINE THE K-WAVE GRID
    % =========================================================================
    
    % set the size of the perfectly matched layer (PML)
    PML_X_SIZE = 10;            % [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 = 256 - 2*PML_X_SIZE;   % [grid points]
    Ny = 512 - 2*PML_Y_SIZE;   % [grid points]
    Nz = 256 - 2*PML_Z_SIZE;   % [grid points]
    
    % set desired grid size in the x-direction not including the PML
    x = 38e-3;                  % [m]
    y = 50e-3;                  % [m]
    z = 38e-3;                  % [m]
    
    % calculate the spacing between the grid points
    dx = z/Nz;                  % [m]
    dy = dx;                  % [m]
    dz = dx;                  % [m]
    
    % create the k-space grid
    kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
    
    % =========================================================================
    % DEFINE THE MEDIUM PARAMETERS
    % =========================================================================
    
    c0 = 1497;                    % c in water [m/s]
    rho0 = 1000;                  % density of water [kg/m^3]
    medium.alpha_coeff = 0.75;    % [dB/(MHz^y cm)]
    medium.alpha_power = 1.5;
    medium.BonA = 6;              % represent nonlinearity properties of the medium
    
    % create the time array
    t_end = 40e-6;                  % [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 = 3.5e6;     % [Hz]
    tone_burst_cycles = 10;
    
    % 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 [SOURCE]
    % =========================================================================
    
    % 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 = 25;  	% 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 = 1530;                % sound speed [m/s]
    transducer.focus_distance = 10.5e-3;          % focus distance [m]
    transducer.elevation_focus_distance = 16e-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 = 128;
    % 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 ULTRASOUND TRANSDUCER [SENSOR]
    % =========================================================================
    
    % physical properties of the transducer
    transducer2.number_elements = 128;  	% total number of transducer elements
    transducer2.element_width = 2;       % width of each element [grid points]
    transducer2.element_length = 25;  	% length of each element [grid points]
    transducer2.element_spacing = 0;  	% spacing (kerf  width) between the elements [grid points]
    transducer2.radius = inf;            % radius of curvature of the transducer [m]
    
    % calculate the width of the transducer in grid points
    transducer_width2 = transducer2.number_elements*transducer2.element_width ...
        + (transducer2.number_elements - 1)*transducer2.element_spacing;
    
    % use this to position the transducer in the middle of the computational grid
    transducer2.position = round([1, Ny/2 - transducer_width2/2, Nz/2 - transducer2.element_length/2]);
    
    % properties used to derive the beamforming delays
    transducer2.sound_speed = 1530;                % sound speed [m/s]
    transducer2.focus_distance = 10.5e-3;          % focus distance [m]
    transducer2.elevation_focus_distance = 16e-3;  % focus distance in the elevation plane [m]
    transducer2.steering_angle = 0;                % steering angle [degrees]
    
    % apodization
    transducer2.transmit_apodization = 'Hanning';
    transducer2.receive_apodization = 'Rectangular';
    
    % define the transducer elements that are currently active
    number_active_elements2 = 128;
    transducer2.active_elements = ones(transducer2.number_elements, 1);
    
    % create the transducer using the defined settings
    transducer2 = makeTransducer(kgrid, transducer2);
    
    % print out transducer properties
    transducer2.properties;
    
    % =========================================================================
    % DEFINE THE MEDIUM PROPERTIES
    % =========================================================================
    
    % define a random distribution of scatterers for the medium
    background_map_mean = 1;
    background_map_std = 0.008;
    
    % define the thicknesses of layer phantom in x direction
    d_layer_1 = 5.4e-3;          % [m]
    d_layer_2 = 4.9e-3;          % [m]
    d_layer_3 = 4.4e-3;          % [m]
    d_layer_4 = 6.3e-3;          % [m]
    
    % define the thicknesses of layer phantom in x direction
    position_x1 = 1:round(Nx*(d_layer_1/x));
    position_x2 = round(Nx*(d_layer_1/x))+1:round(Nx*((d_layer_1+d_layer_2)/x));
    position_x3 = round(Nx*((d_layer_1+d_layer_2)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x));
    position_x4 = round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))+1:round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x));
    
    % define the width of phantom in y and z direction
    width_y = 29e-3; % width_phantom_in y direction [grid points]*dy = [m]
    width_z = 35e-3; % width_phantom_in z direction [grid points]*dz = [m]
    
    % define position of each layer in y and z direction
    position_y = round(Ny*(((y-width_y)/2)/y))+1:round(Ny*(y-(y-width_y)/2)/y);
    position_z = round(Nz*(((z-width_z)/2)/z))+1:round(Nz*(z-(z-width_z)/2)/z);
    
    % define the dimension of each layer (equal to its background map)
    width_backgroundmap_layer_1 = [round(Nx*(d_layer_1/x)), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
    width_backgroundmap_layer_2 = [round(Nx*((d_layer_1+d_layer_2)/x))-(round(Nx*(d_layer_1/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
    width_backgroundmap_layer_3 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))-(round(Nx*((d_layer_1+d_layer_2)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
    width_backgroundmap_layer_4 = [round(Nx*((d_layer_1+d_layer_2+d_layer_3+d_layer_4)/x))-(round(Nx*((d_layer_1+d_layer_2+d_layer_3)/x))), round(Ny*(y-(y-width_y)/2)/y)-round(Ny*(((y-width_y)/2)/y)), round(Nz*(z-(z-width_z)/2)/z)-round(Nz*(((z-width_z)/2)/z))];
    
    % % define the properties of the propagation medium
    % % water
    medium.sound_speed = c0*ones(Nx, Ny, Nz);      % [m/s]
    medium.density = rho0*ones(Nx, Ny, Nz);        % [kg/m^3]
    
    % layer_1
    background_map_layer_1 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_1);
    medium.sound_speed(position_x1, position_y, position_z) = 1520.*background_map_layer_1;        % [m/s]
    medium.density(position_x1, position_y, position_z) = 1281.*background_map_layer_1;            % [kg/m^3]
    
    % layer_2
    background_map_layer_2 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_2);
    medium.sound_speed(position_x2, position_y, position_z) = 1218.*background_map_layer_2;    % [m/s]
    medium.density(position_x2, position_y, position_z) = 1510.*background_map_layer_2;    % [kg/m^3]
    
    % layer_3
    background_map_layer_3 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_3);
    medium.sound_speed(position_x3, position_y, position_z) = 1325.*background_map_layer_3;       % [m/s]
    medium.density(position_x3, position_y, position_z) = 1643.*background_map_layer_3;       % [kg/m^3]
    
    % layer_4
    background_map_layer_4 = background_map_mean + background_map_std*randn(width_backgroundmap_layer_4);
    medium.sound_speed(position_x4, position_y, position_z) = 1429.*background_map_layer_4; % [m/s]
    medium.density(position_x4, position_y, position_z) = 1684.*background_map_layer_4; % [kg/m^3]
    
    % =========================================================================
    % RUN THE SIMULATION
    % =========================================================================
    
    % set the input settings
    % plot medium map
    input_args = {...
        'PMLInside', false, 'PlotSim', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ...
        'DataCast', data_cast};
    
    % run the simulation if set to true, otherwise, load previous results from disk
    if run_simulation
        total_scan_lines = 97;
        total_scan_lines_test = 2;
    
        for scan_line = 1 : total_scan_lines_test
            fprintf('Calculating scan line %d of %d\n', scan_line, total_scan_lines_test);
    
            % define the transducer elements that are currently active
              transducer.active_elements = zeros(transducer.number_elements, 1);
              transducer.active_elements(scan_line:scan_line+31)=1;
    
            RF_data_SB_transd_shift = kspaceFirstOrder3D(kgrid, medium, transducer, transducer2, input_args{:});
    
        end
    
    end

    Thanks

    Makcik

    Posted 11 years ago #

13. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    Hi Makcik,

    For each individual simulation you run, you will obtain 128 time signals. These are equivalent to the raw RF data from an ultrasound machine. If you run 97 simulations changing the transmit elements each time, in the end you will finally obtain 97*128 time series. You could store these using something like:

    RF_data(sim_index, :, :) = kspaceFirstOrder3D(...).

    Glancing at your code, it looks like you set 128 elements on both transmit and receive. Also, see my previous content about the frequency of the source and the frequency supported by the grid.

    Brad.

    Posted 11 years ago #

14. maktjik
    

    Member
    Joined: Mar '13
    Posts: 18

    Hi Brad,

    Thanks for your suggestion, now I can simply save the final output in 3D matrix.
    I also reduce tone_burst center frequency, because I just check the spectrum exceed the maximum supported frequency.

    And what about the plane wave case I asked before?
    - in transmit: each element (total: 128 elements), propagate the signal in the same time, focus_distance is infinite, and elevation_focus_distance is finite.
    - in receive: total 128 elements also receive the echo in the same time.

    Can I simulate it by simply define 128 total elements transducer (with also 128 active elements), define focus_distance to be inf, and set them as KSpaceFirstOrder3D input arguments in the places of source and sensor? Another parameters, such as computational grid, medium parameters, medium properties, input signal and time array are still the same as before.

    Thanks,

    Makcik

    Posted 11 years ago #

15. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    Hi Makcik,

    Your approach sounds right for simulating plane wave imaging using a conventional diagnostic ultrasound transducer. In this case, you could use the same transducer object for both source and sensor as the active elements will be 128 in both cases.

    Brad.

    Posted 11 years ago #

RSS feed for this topic

Reply

You must log in to post.