k-Wave

1. Tim Wang
   

   Member
   Joined: Feb '12
   Posts: 2

   Hello Dr. Cox,

   The alpha_coeff in the code now is set to 0.75, it is the propagation medium absorption, right? Where can I find the target source absorption coefficient right now?

   Thanks,
   Tim

   Posted 12 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Hi Tim,

   Instead of being a single number, alpha_coeff can be a matrix the same size as your image, so alpha_coeff takes a different value for each pixel (in 2D) or voxel (in 3D). If you want just two values - one or the background and one for an object - you have to define values for all pixels. Here's an example for a disc with higher absorption:
   
object = makeDisc(kgrid.Nx, kgrid.Ny, kgrid.Nx/2, kgrid.Ny/2, 5);
medium.alpha_coeff = ones(kgrid.Nx,kgrid.Ny)*0.75;
medium.alpha_coeff(object==1) = 2;


   Is something like this what you're trying to do?

   Ben

   Posted 12 years ago #

3. Tim Wang
   

   Member
   Joined: Feb '12
   Posts: 2

   Thanks Dr. Cox. This is what I was trying to do.

   So originally, if alpha_coeff =0.75, that means both disc source and the background meduim has the same absorption coeff.? (All pixels coeff = 0.75?)

   Tim

   Posted 12 years ago #

4. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   That's right, yes.

   Hope you get it working!

   Ben

   Posted 12 years ago #

5. DuuMushishi
   

   Member
   Joined: Aug '13
   Posts: 4

   Hello Dr.Cox
   I'm trying to run a simulation whit different absorption coeff. in the medium (i.e 0.75,0.5) and I'm using the B-mode image example file. However, when the alfa_coeff is change to a matrix, the simulation encounter the following error:

   Error in ==> kspaceFirstOrder_inputChecking at 83
   checkFieldNames(medium, {'sound_speed', 'sound_speed_ref', 'density', 'alpha_coeff', 'alpha_power',
   'alpha_mode', 'alpha_filter', 'alpha_sign', 'BonA'});

   Error in ==> kspaceFirstOrder3D at 495
   kspaceFirstOrder_inputChecking;

   Error in ==> simulation at 135
   sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});

   Is there any other value that I must set/correct in order to properly simulate a medium with two different absorption properties?

   Andre

   Posted 10 years ago #

6. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Andre,

   Is there an error message above where it says "Error in ..." ?

   Brad.

   Posted 10 years ago #

7. DuuMushishi
   

   Member
   Joined: Aug '13
   Posts: 4

   Oh, I'm sorry, I didn´t copy the entire error message

   ??? Error using ==> checkFieldNames at 51
   alfa_coeff is not a valid field for the structure medium

   Posted 10 years ago #

8. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   The problem is just a small typo: alfa_coeff must be spelt alpha_coeff.

   Posted 10 years ago #

9. DuuMushishi
   

   Member
   Joined: Aug '13
   Posts: 4

   OK, I've just corrected that typo and I got this error:

   ....medium.alpha_coeff = ones(kgrid.Nx,kgrid.Ny,kgrid.Nz)*0.2; % [dB/(MHz^y cm)]
   medium.alpha_power = 1.0;
   medium.alpha_mode = 'no_dispersion';.....

   ??? Attempt to grow array along ambiguous dimension.

   Error in ==> example_us_bmode_linear_transducer at 189
   medium.alpha_coeff(scattering_region1==1)=0.6;

   So I assumed that it was because the scattering_region1 had different size than kgrid (Ny_tot = Ny + number_scan_lines*transducer.element_width;) thus, I changed the dimension of the alpha_coeff matrix to:

   medium.alpha_coeff = ones(Nx_tot,Ny_tot,Nz_tot)*0.2;

   but then this error comes out:

   ??? Array dimensions must match for binary array op.

   Error in ==> kspaceFirstOrder_createAbsorptionVariables at 40
   absorb_tau = -2*medium.alpha_coeff.*c.^(medium.alpha_power -
   1);

   Error in ==> kspaceFirstOrder3D at 705
   kspaceFirstOrder_createAbsorptionVariables;

   Error in ==> example_us_bmode_linear_transducer at 247
   sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer,
   transducer, input_args{:});

   Am I missing something?
   Thanks,
   Andre

   Posted 10 years ago #

10. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,640

    Hi Andre,

    This particular simulation works by assigning larger heterogenous maps of sound and density, and then assigning a portion of these to medium.sound_speed and medium.density for the simulation of each scan line.

    If you want to follow the same format for the absorption, you should define an additional absorption coefficient variable that is Nx_tot by Ny_tot by Nz_tot, and then for the simulation of each scan line, assign the appropriate portion of this to medium.alpha_coeff. You can see how this is done for the sound speed and density at line 235/236.

    Hope that helps,

    Brad.

    Posted 10 years ago #

11. DuuMushishi
    

    Member
    Joined: Aug '13
    Posts: 4

    Hello Brad,

    I followed your instructions and it seems the simulation is now running without problems. Thank you for your quick reply. I really appreciate your help.

    Thanks again
    Andre

    Posted 10 years ago #

12. zengqw2021
    

    Member
    Joined: Jun '21
    Posts: 26

    Dear treeby,
    I’m a new user of k-Wave toolbox. It is great. Now, I'm trying to run a simulation with different absorption coeff. The alfa_coeff is change to a matrix, the simulation encounter the following error:
    Matrix dimensions must agree.
    Error in kspaceFirstOrder2D (line 864)
    (rhox + rhoy) ...
    Error in PAT_Heterogeneous (line 93)
    p0_recon = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:});

    %%The code is as follows:

    clearvars;clear;close all;clc;

    % =========================================================================
    % SIMULATION
    % =========================================================================

    % assign the grid size and create the computational grid
    PML_size = 20; % size of the PML in grid points
    Nx = 256 - 2 * PML_size; % number of grid points in the x direction
    Ny = 256 - 2 * PML_size; % number of grid points in the y direction
    x = 20e-3; % total grid size [m]
    y = 20e-3; % total grid size [m]
    dx = x / Nx; % grid point spacing in the x direction [m]
    dy = y / Ny; % grid point spacing in the y direction [m]
    kgrid = kWaveGrid(Nx, dx, Ny, dy);

    % create initial pressure distribution
    medium.sound_speed = 330; % [m/s]
    medium.alpha_coeff = ones(kgrid.Nx,kgrid.Ny)*0.75; % [dB/(MHz^y cm)] %ones(kgrid.Nx,kgrid.Ny)*
    medium.alpha_power = 1.5;
    medium.alpha_mode = 'no_dispersion';

    Cp_heat_cap = 1.002; % the specific heat capacity at constant pressure
    Heat_beta = 1/273.15; % the coefficient of the thermal expansion

    disc_H0 = 30; % [kJ/m3]
    disc_x_pos = floor(Nx/2); % [grid points]
    disc_y_pos = floor(Ny/2); % [grid points]
    disc_radius = 4/2/sqrt(log(2)); % [grid points]
    disc_1 = Heat_beta*medium.sound_speed^2/Cp_heat_cap*disc_H0*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); % Heat_beta*medium.sound_speed^2/Cp_heat_cap*

    p0 = disc_1 ;
    % imagesc(p0);
    source.p0 = disc_1 ;
    init_p0 = source.p0; %for output

    % smooth the initial pressure distribution and restore the magnitude
    source.p0 = smooth(source.p0, true);

    % define a centered Cartesian circular sensor
    sensor_radius = 2e-3; % [m]
    sensor_angle = 2*pi ; % [rad]
    sensor_pos = [0, 0]; % [m]
    num_sensor_points = 30;
    cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle);

    % assign to sensor structure
    sensor.mask = cart_sensor_mask;

    % create the time array
    t1=kgrid.makeTime(medium.sound_speed);

    % set the input arguements: force the PML to be outside the computational
    % grid; switch off p0 smoothing within kspaceFirstOrder2D
    input_args1 = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};

    % %% save the simulation animations as a movie
    % input_args2 = {'RecordMovie', true, 'MovieName', 'example_movie'};
    % input_args = [input_args1 input_args2];
    input_args = input_args1;

    % run the simulation
    sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

    % add noise to the recorded sensor data
    signal_to_noise_ratio = 40; % [dB]
    sensor_data = addNoise(sensor_data, signal_to_noise_ratio, 'peak');

    % create a second computation grid for the reconstruction to avoid the
    % inverse crime
    Nx = 300; % number of grid points in the x direction
    Ny = 300; % number of grid points in the y direction
    dx = x/Nx; % grid point spacing in the x direction [m]
    dy = y/Ny; % grid point spacing in the y direction [m]
    kgrid_recon = kWaveGrid(Nx, dx, Ny, dy);

    % use the same time array for the reconstruction
    kgrid_recon.setTime(kgrid.Nt, kgrid.dt);

    % reset the initial pressure
    source.p0 = 0;

    % assign the time reversal data
    sensor.time_reversal_boundary_data = sensor_data;

    % run the time-reversal reconstruction
    p0_recon = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:});

    Could you give me some advice?
    I appreciate any help in this regard.
    Zeng

    Posted 2 years ago #

13. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,640

    Hi Zeng,

    The second simulation is run on a different sized grid (to avoid the inverse crime). So you will either need to use the same grid, or redefine medium.alpha_coeff to be Nx by Ny.

    Note though, you probably want to turn absorption off anyway for the reconstruction (or invert it as described in the attenuation compensation examples).

    Brad.

    Posted 2 years ago #

14. zengqw2021
    

    Member
    Joined: Jun '21
    Posts: 26

    Hi Brad,
    Ok, I see. Is there a way for me to turn absorption off during reconstruction?

    Zeng

    Posted 2 years ago #

15. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,640

    If you don't specify the absorption, the simulation will be lossless. If you have already specified the medium properties, you can remove it by:

    medium = rmfield(medium, 'alpha_coeff');
    medium = rmfield(medium, 'alpha_power');

    Posted 2 years ago #

RSS feed for this topic

Reply

You must log in to post.