k-Wave

1. maktjik
   

   Member
   Joined: Mar '13
   Posts: 18

   Dear Dr. Treeby and Dr. Cox,

   I have some questions regarding my simulation. In my case, I want to simulate plane wave ultrasound (128 elements array) propagation through heterogeneous medium. My medium is steel cylinder in water (3D) that I generate by integrating 2D images (see the 2nd section of my code and the link for the image source).

   From several example I got confuse if I should use c in steel, c in water or c average of the medium to define:
   1. time array
   2. input_signal
   3. transducer sound speed

   Thank you.

   Here is the code:
   

   % =========================================================================
   % ULTRASOUND SIMULATION
   % Linear array transducer (128 element)
   % Transmit: Planewave, 128 active elements
   % Receive : 128 active elements
   % Medium  : steel cylinder in water (256 x 512 x 256)
   % =========================================================================
   
   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 = 1498;                    % c in water [m/s]
   rho0 = 1000;                  % density of water [kg/m^3]
   c_needle = 5940;              % c in steel [m/s]
   rho_needle = 7800;            % density of steel [kg/m^3]
   c_avg = 3719;                 % average c  of the medium (needle and water) [m/s]
   rho_avg = 4400;               % average density  of the medium (needle and water) [kg/m^3]
   
   % create the time array
   t_end = (Nx*dx)*2.2/c0;                % [s]
   kgrid.t_array = makeTime(kgrid, c_needle, [], 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 [SOURCE AND SENSOR]
   % =========================================================================
   
   % 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]);
   
   transducer.sound_speed = c0;               % sound speed [m/s]
   transducer.elevation_focus_distance = 16e-3;  % focus distance in the elevation plane [m]
   
   % 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 MEDIUM PROPERTIES
   % =========================================================================
   
   % define a random distribution of scatterers for the medium
   background_map_mean = 1;
   background_map_std = 0.008;
   
   % % define the properties of the propagation medium
   load c_needle_water_171103;
   medium.sound_speed = c_needle_water_171103;      % [m/s]
   
   load rho_needle_water_171103;
   medium.density = rho_needle_water_171103;        % [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
   
      RF_data_PW_needle_171104 = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
   
      % save the scan lines to disk
      save RF_data_PW_needle_171104 RF_data_PW_needle_171104;
   else
   %    % load the scan lines from disk
      load RF_data_PW_needle_171104
   end

   =======================================================================
   and this is the code for creating the medium:

   clear all
   close all
   
   [data]=imread('needle_2D.bmp');
   data_bw=double(data(:,:,1));
   
   [xx,yy]=size(data_bw);
   clear data
   
   for i=1:xx
       for j=1:yy
           if data_bw(i,j)==1
               c(i,j)=5940;       % average c  of the medium (phantom and water) [m/s]
               rho(i,j)=7800;     % average density  of the medium (phantom and water) [kg/m^3]
           elseif data_bw(i,j)==0
               c(i,j)=1498;       % c in water [m/s]
               rho(i,j)=1000;     % density of water [kg/m^3]
           end;
       end;
   end;
   
   cit_3d=zeros(xx,yy,492);
   
   for i=1:492
       c_3d(:,:,i)=c;
       rho_3d(:,:,i)=rho;
   end;
   
   n=zeros(236,492,236);
   m=zeros(236,492,236);
   
   [x,y,z]=size(n);
   
   for i=1:x
       for j=1:y
           for k=1:z
               c_needle_water_171103(i,j,k)=c_3d(i,k,j);
               rho_needle_water_171103(i,j,k)=rho_3d(i,k,j);
           end
       end
   end
   
   save c_needle_water_171103  c_needle_water_171103
   save rho_needle_water_171103  rho_needle_water_171103

   =============================================================
   and here is the image:

   https://drive.google.com/file/d/1nLgtWFzVAlHxV_gxr_dJjuELNnMxi38h/view?usp=sharing

   Posted 7 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi maktijk,

   1. When you use makeTime, you can directly pass the heterogeneous sound speed. This will use the highest sound speed to calculate the time step, and the lowest sound speed to calculate the number of time steps (although in your example you set this manually).

   Something to note, by default, k-Wave will use the highest sound speed in the k-space correction term to make the simulation unconditionally stable (at least in the case of no absorption). However, this means dispersion will accumulate in other regions of the domain where the sound speed is lower. In your case, I would suggest setting medium.sound_speed_ref to the background sound speed (water), so the propagation is dispersion free in these regions. In this case, the simulation will no longer be unconditionally stable. You can query the largest stable time step using the checkStability function.

   2. For the input signal, use the sound speed where the transducer is (water in your case).

   3. The transducer sound speed is only used to set the focus delays, so I would use the sound speed in water.

   Hope that helps,

   Brad.

   Posted 7 years ago #

RSS feed for this topic

Reply

You must log in to post.