% Example of using an ultrasound transducer as part of a time reversal % simulation. % % author: Bradley Treeby % date: 3rd October 2014 % last update: 3rd October 2014 % % This function is part of the k-Wave Toolbox (https://bug.medphys.ucl.ac.uk/kwave) % Copyright (C) 2009-2014 Bradley Treeby and Ben Cox % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see . clear all; % simulation settings DATA_CAST = 'single'; % ========================================================================= % 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 = 64 - 2*PML_X_SIZE; % [grid points] Ny = 128 - 2*PML_Y_SIZE; % [grid points] Nz = 64 - 2*PML_Z_SIZE; % [grid points] % set desired grid size in the x-direction not including the PML x = 40e-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 % ========================================================================= % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] medium.density = 1000; % [kg/m^3] % create the time array t_end = 40e-6; % [s] kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end); % ========================================================================= % DEFINE THE SOURCE % ========================================================================= % create initial pressure distribution source_strength = 1e6; p0 = source_strength*(makeBall(Nx, Ny, Nz, round(25e-3/dx), Ny/2, Nz/2, 3)... + makeBall(Nx, Ny, Nz, round(8e-3/dx), Ny/4, Nz/2, 3)); % smooth p0 and assign p0 = smooth(kgrid, p0, true); source.p0 = p0; % ========================================================================= % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= % physical properties of the transducer transducer.number_elements = 72; % total number of transducer elements transducer.element_width = 1; % width of each element [grid points/voxels] transducer.element_length = 12; % length of each element [grid points/voxels] transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points/voxels] 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 = 1540; % sound speed [m/s] transducer.focus_distance = 25e-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 = 'Rectangular'; transducer.receive_apodization = 'Rectangular'; % define the transducer elements that are currently active transducer.active_elements = ones(transducer.number_elements, 1); % create the transducer using the defined settings transducer = makeTransducer(kgrid, transducer); % ========================================================================= % RUN THE SIMULATION % ========================================================================= % set the input settings input_args = {'DisplayMask', transducer.active_elements_mask, ... 'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE], ... 'DataCast', DATA_CAST, 'PlotScale', [-source_strength/4, source_strength/4], ... 'Smooth', false}; % run the simulation sensor_data = kspaceFirstOrder3D(kgrid, medium, source, transducer, input_args{:}); % ========================================================================= % TIME REVERSAL % ========================================================================= % clear the source and sensor structures clear source sensor; % create the source mask from the transducer mask source.p_mask = transducer.active_elements_mask; % get equivalent mask where the matrix values correspond to which physical % element the grid point belongs to indexed_mask = transducer.indexed_active_elements_mask; % reshape to a linear index mask, and remove the zeros indexed_mask = indexed_mask(:); indexed_mask(indexed_mask == 0) = []; % create an empty time varying source source.p = zeros(sum(source.p_mask(:)), size(sensor_data, 2)); % for each grid point in the source mask, assign the appropriate time % series from the sensor_data recorded using the transducer for index = 1:size(source.p, 1) source.p(index, :) = sensor_data(indexed_mask(index), :); end % reverse the source ready for time reversal source.p = flipdim(source.p, 2); % enfore as a dirichlet boundary condition source.p_mode = 'dirichlet'; % run time reversal simulation sensor.record = {'p_final'}; sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); % extract reconstruction p0_recon = sensor_data.p_final; % add first order compensation for only recording over a half plane p0_recon = p0_recon*2; % ========================================================================= % VISUALISATION % ========================================================================= figure; subplot(2, 1, 1); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, squeeze(p0(:, :, Nz/2)), [-source_strength/4, source_strength/4]); axis image; ylabel('x position [mm]') xlabel('y position [mm]') title('Initial Pressure'); subplot(2, 1, 2); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, squeeze(p0_recon(:, :, Nz/2)), [-source_strength/4, source_strength/4]); colormap(getColorMap); axis image; ylabel('x position [mm]') xlabel('y position [mm]') title('Reconstruction');