k-Wave

1. Farnaz
   

   Member
   Joined: Oct '18
   Posts: 20

   Hi,

   First of all thanks for responding to the questions in the forum, it's really helpful.

   secondly, I am a bit confused about pressure and velocity sources. How these two are related? and what is the difference between value 'p', 'p_final' that sensor records.

   I simulate a 2D transducer for plane wave using sensor and source functions. I modified the input shape to a ricker wavelet but the result somehow looks strange.
   If I use a velocity source and don't set sensor.record, then load '/p' from .h5 file what am I exactly looking at?

   Posted 3 years ago #

2. Farnaz
   

   Member
   Joined: Oct '18
   Posts: 20

   Hi,

   Since I am struggling with setting a correct setup for my simulation for very long time, I would like to rephrase my question and post the code.

   I added a b-mode reconstructed image from my simulation here:
   

   There is medium with specific speed of sound and an inclusion with a lower speed of sound on the center which after reconstruction I can only see a reflection from above and below edges in the b-mode image while horizontal reflections are missing and there is an artifact which can be removed if I set the speed of sound in the inclusion to a higher value compared to background. However, in real setup you can see a similar inclusion even with single planwave shot.
   What am I doing wrong?

   I have transducer with 384 elements (only 192 elements in middle are activated). I use a plane wave. The pitch of the transducer is 200 micrometer. My medium is a medium of size 7.6*3.8 Cm.

   Here is the setup of my code:

   clearvars;
   tic
   % =========================================================================
   % Define Grid
   % =========================================================================
   
   % set the size of the perfectly matched layer (PML)
   pml_y_size = 64;                % [grid points]
   pml_x_size = 96;                % [grid points]
   
   % set total number of grid points not including the PML
   Ny = 3072;%
   Nx = 1536;%       % [grid points]
   
   dy = 2*3.84e-2/Ny;    	% grid point spacing in the y direction [m]
   dx = dy;            % grid point spacing in the x direction [m]
   
   kgrid = kWaveGrid(Nx, dx, Ny, dy);
   
   % =========================================================================
   % Define Medium
   % =========================================================================
   
   % define the properties of the propagation medium
   medium.sound_speed = 1540;  % [m/s]
   medium.alpha_coeff = 0.5;  % [dB/(MHz^y cm)]
   medium.alpha_power = 1.05;
   medium.BonA = 10;
   rho0 = 900;
   
   c0 = 1300 ;
   % create the time array
   t_end = (Nx * dx) * 2.2 / c0;   % [s]
   kgrid.makeTime(c0, [], t_end);
   
   % =========================================================================
   % Define Medium
   % =========================================================================
   % define a random distribution of scatterers for the medium
   background_map_mean = 1;
   SD = 1 ;
   num_run = 1 ;
   medium = repmat(medium,num_run,1) ; 
   
   rng('shuffle')
   background_map = (0.9 - 0.03) + (0.06).* sprand(Nx, Ny,0.01);
   density_map = rho0 * ones(Nx, Ny).* background_map;
   sound_speed_map = 1540 * ones(Nx,Ny) ;
   
   theta = 0;
   
   R_x = 150 ;
   R_y = 150 ; 
   
   C_x = Nx/2 ;
   C_y = Ny/2; 
   
   scattering_region = makeEllipsoid2D([Nx,Ny],[C_x , C_y ], [R_x, R_y],theta) ;
   % assign region
   sound_speed_map(scattering_region == 1) = 1300  ;
   
   medium.sound_speed = sound_speed_map(:, :);
   medium.density = density_map(:,:);
   
   % =========================================================================
   % Define Source middle
   % =========================================================================
   kerf = 1;
   groupspacing = 7;
   element_num = 192*2;
   source_shape = reshape((1:groupspacing)' + (0:element_num-1)*(kerf+groupspacing), 1, []);
   
   % % define source mask for a linear transducer with an odd number of elements
   x_offset = 1;          % [grid points]
   source_m.u_mask = zeros(Nx, Ny);
   source_m.u_mask(x_offset , source_shape ) = 1  ; 
   
   %activate middle elements
   active_elements = element_num/2  ;
   sensor_m_offset_mask = groupspacing*element_num/4+kerf*(element_num/4)-1 ; %groupspacing*element_num/7+kerf*(element_num/7)-1 ;
   source_m.u_mask(x_offset , 1:sensor_m_offset_mask) = 0 ;
   source_m.u_mask(x_offset , Ny-sensor_m_offset_mask:Ny)=0 ; 
   
   source_strength = 1e6 ;
   % define the properties of the tone burst used to drive the transducer
   sampling_freq = 1/kgrid.dt;     % [Hz]
   steering_angle_m = 0;            % [deg]
   tone_burst_freq = 5e6;          % [Hz]
   [rw,t] = ricker(tone_burst_freq,120, kgrid.dt);
   input_signal = rw ;
   input_signal = (source_strength ./ (c0 * rho0)) .* input_signal ;
   source_m.ux = input_signal ;
   % =========================================================================
   % Define Sensor
   % =========================================================================
   
   sensor_m.mask = source_m.u_mask ;
   %sensor_m.directivity_angle = ones(Nx,Ny)*pi/2 ;
   sensor_m.directivity_size = kgrid.dx;
   
   input_args = {...
       'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size], ...
        'DataRecast', true, 'PlotSim', true, 'PlotPML', true , 'PlotScale', 'auto','DataCast', 'gpuArray-single'} ; %,'RecordMovie', true , 'MovieName', 'Final_simulation_1phantome', 'MovieProfile', 'MPEG-4'};
   
   % run the simulation
   
   filename_m = sprintf('%s_%d.h5','test') ;
   kspaceFirstOrder2D(kgrid, medium, source_m,sensor_m, input_args{:},'SaveToDisk',filename_m ) ;

   I run the c++ code and save the output and again load it here do postprocessing and resample it and reconstruct the 'transducer_data_resample' with a standard delay and sum :

   %%
   sensor_data = h5read('test_out.h5','/p') ;
   size_data = size(sensor_data) ;
   transducer_data = zeros(active_elements , size_data(2)) ;
   element = 1 ;
   %
   % =========================================================================
   % Extract transducer data
   % =========================================================================
   
   sumnum = groupspacing-1 ; 
   
   for i = 1:groupspacing:size_data(1)- sumnum
       transducer_data(element , :) = (sensor_data(i,:)+ sensor_data(i+1,:)+ sensor_data(i+2,:)+ sensor_data(i+3,:)+sensor_data(i+4,:)+ sensor_data(i+5,:)+ sensor_data(i+6,:)) ;
       element = element+1 ;
   end 
   
   originalFs = 1/kgrid.dt ;
   desiredFs = 40e6;
   [p,q] = rat(desiredFs / originalFs); %originalFs * p / q ;
   
   for i = 1 : active_elements
   transducer_data_resample(i,:) = resample(transducer_data(i,:),p,q);
   end

   additional functions to run the code:
   =========================================================================
   

   function [rw,t] = ricker(f,n,dt,t0,t1)
   %RICKER creates an causal ricker wavelet signal
   %
   %   RICKER creates and plots a default causal ricker wavelet with:
   %
   %       peak frequency   = 20 Hz
   %       sampling time    = 0.001 seconds
   %       number of points = 100;
   %       peak location    = 1/F = 1/20Hz
   %
   %   RW = RICKER(...) returns the default wavelet in RW.
   %
   %   [RW,T] = RICKER(...) returns the time vector in T.
   %
   %   Specifying parameters of peak frequency (F, Hz), number of points (N),
   %   and sampling time (DT) are specified by the syntax:
   %
   %       [RW,T] = RICKER(F)
   %       [RW,T] = RICKER(F,N)
   %       [RW,T] = RICKER(F,N,DT)
   %
   %   [RW,T] = RICKER(F,N,DT,T0) creates a ricker wavelet with peak centered
   %   at T0.
   %
   %   [RW,T] = RICKER(F,N,DT,T0,T1) creates a 2 dimensional symmetric
   %   ricker wavelet with sift in 1st dimension of T0 and second dimension of
   %   T1.
   %
   %   Example 1:
   %       ricker % plots a 20 Hz Ricker Wavelet over 0.1 seconds
   %
   %   Example 2:
   %    % create a ricker wavelet with 40 Hz, 200 points, and 0.02 s between
   %    % samples
   %    [rw,t] = ricker(40,200,0.002);
   %    plot(t,rw), xlabel('Time'), ylabel('Amplitude')
   
   % Define inputs if needed
   switch nargin
       case 0
           f  = 20;
           n  = 100;
           dt = 0.001;
           t0 = 1/f;
           is2d = false;
       case 1
           n = 100;
           dt = 0.001;
           t0 = 1/f;
           is2d = false;
       case 2
           dt = 0.001;
           t0 = 1/f;
           is2d = false;
       case 3
           t0 = 1/f;
           is2d = false;
       case 4 % use all values
           is2d = false;
       case 5 % use all inputs
           is2d = true;
       otherwise
           warning('RICKER:tooManyInputs','Ignoring additional inputs')
   end
   
   % Create the wavelet and shift in time if needed
   T = dt*(n-1);
   t = 0:dt:T;
   tau = t-t0;
   if ~is2d
       s = (1-tau.*tau*f^2*pi^2).*exp(-tau.^2*pi^2*f^2);
   else
       [t1,t2] = meshgrid(tau,t-t1);
       s = (1-(t1.^2+t2.^2)*f^2*pi^2).*exp(-(t1.^2+t2.^2)*pi^2*f^2);
   end
   
   if nargout == 0
           plot(t,s)
           xlabel('Time (s)')
           ylabel('Amplitude')
           title('Ricker Wavelet')
   else
           rw = s;
   end

   =============

   function ellipsoid = makeEllipsoid2D(grid_size, position, radii,theta, plot_ellipsoid, binary)
   A = theta ;
   if nargin < 5 || isempty(plot_ellipsoid)
       plot_ellipsoid = false;
   end
   
   % check for binary input
   if nargin < 6 || isempty(binary)
       binary = false;
   end
   
   % force integer grid values
   Nx = round(grid_size(1));
   Ny = round(grid_size(2));
   
   % assign center and radius
   cx = position(1);
   cy = position(2);
   rx = radii(1);
   ry = radii(2);
   
   % check for zero values for center position, and set to middle of grid
   if cx == 0
       cx = floor(Nx / 2) + 1;
   end
   if cy == 0
       cy = floor(Ny / 2) + 1;
   end
   
   % create grid axes
   [x_index, y_index] = ndgrid(1:Nx, 1:Ny) ;
   
   % create ellipsoid based on code from https://stackoverflow.com/questions/36420023/generating-a-3d-binary-mask-of-geometric-shapes-in-matlab
   %ellipsoid = ((x_index - cx).^2 / (rx.^2)) + ((y_index - cy).^2 / (ry.^2)) <= 1;
   ellipsoid = (((x_index-cx)*cos(A) + (y_index-cy)*sin(A)).^2/(rx.^2))+(((x_index-cx)*sin(A) - (y_index-cy)*cos(A)).^2 /ry.^2)<= 1;
   % convert to double precision if not binary
   if ~binary
       ellipsoid = double(ellipsoid);
   end
   
   % plot results
   if plot_ellipsoid
       voxelPlot(double(ellipsoid));
   end

   Posted 3 years ago #

3. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Farnaz,

   There is a more detailed description of pressure and velocity sources in the k-Wave user manual, as well as simple examples in the first two Time Varying Source Problems examples (you can open these in MATLAB and run the code).

   I assume don't see the left and right edges of the inclusion, as the incident plane wave and the detector are placed above it, and the density of the scatterers is constant throughout your phantom (giving uniform speckle). This means there is no energy reflected from these edges back to the transducer. You could try changing the density (or magnitude) of the scatterers within the ellipse.

   Not sure about the artefact I'm afraid without studying your reconstruction more closely. Does it change if you run the simulation for longer?

   Brad

   Posted 3 years ago #

4. Farnaz
   

   Member
   Joined: Oct '18
   Posts: 20

   Dear Brad,

   Thank a lot for your response.
   For the artifact I was suspicious of PML size for a while, since in a post you mentioned bigger PML is not necessarily better...
   I realized there is a problem with speed of sound definition of the medium. I could remove the artifact with modifying the speed of sound, the variation of speed of sound in the medium affects the artifact. If I have an inclusion in which speed of sound drastically changes compared to background such artifact appears e.g. 1300 vs. 1700.
   About the inclusion shape, changing the distribution of speckles worked to get the whole shape.

   Best,
   Farnaz

   Posted 3 years ago #

RSS feed for this topic

Reply

You must log in to post.