Hi,
I reconstructed initial pressure using 3D Time Reversal Reconstruction For A Planar Sensor in heterogeneous medium. But i can not get good reconstructed images.The codes I use are as following:
%% create the computational grid
PML_size = 10; % size of the PML in grid points
% Nx = 64 * scale - 2 * PML_size; % number of grid points in the x direction
% Ny = 64 * scale - 2 * PML_size; % number of grid points in the y direction
% Nz = 32 * scale - 2 * PML_size; % number of grid points in the z direction
Nx = 190 - 2 * PML_size; % number of grid points in the x direction
Ny = 270 - 2 * PML_size; % number of grid points in the y direction
Nz = 126 - 2 * PML_size; % number of grid points in the z direction
dx = 1.25e-3/2 ; % grid point spacing in the x direction [m]
dy = 1.25e-3/2; % grid point spacing in the y direction [m]
dz = 1.25e-3/2; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
%% load the speed of sound and the tissue density and the initial pressure distribution from an image and scale
soundspeed=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\a_segmentation\mat\soundspeed_160_240_96.mat');
tissdst=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\a_segmentation\mat\tissdst_160_240_96.mat');
p0=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\b_initial_pressure\mat\initial_dose_160_240_53.mat');
%% get the size of soundspeed
[r,c,h]=size(soundspeed);
%% define the speed of sound(m/s)
for k=1:h
for i=1:r
for j=1:c
if soundspeed(i,j,k)<=343
soundspeed(i,j,k)=1500;
end
end
end
end
%% define the tissue density(kg/m3)
for k=1:h
for i=1:r
for j=1:c
if tissdst(i,j,k)<=1.2
tissdst(i,j,k)=1000;
end
end
end
end
p0_binary=zeros(Nx,Ny,Nz);
p0_binary(6:165,6:245,6:58)=p0;
% p0_binary=zeros(Nx,Ny,Nz);
% p0_binary(1:32,1:32,5)=im;
% smooth the initial pressure distribution and restore the magnitude
% p0 = smooth(kgrid, p0_binary, true);
p0=p0_binary;
%% assign to the source structure
source.p0 = p0;
%% define the properties of the propagation medium
sound_speed=ones(Nx,Ny,Nz)*1500;
sound_speed(6:165,6:245,6:101)=soundspeed;
medium.sound_speed = sound_speed; % [m/s]
density=ones(Nx,Ny,Nz)*1000;
density(6:165,6:245,6:101)=tissdst;
medium.density=density; % [kg/m3]
% medium.alpha_coeff = 0.75; % power law absorption prefactor [dB/(MHz^y*cm)]
% medium.alpha_power = 1.5; % power law absorption exponent
%% define a binary planar sensor
sensor.mask = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);
sensor.mask(6:2:165,6:2:245,102) = 1;
% sensor.mask(1:32, 1:32,1) = 1;
%% define frequence and bandwidth of sensor
center_freq =1e6; %!!!!!!定义探测器的中心频率
bandwidth = 100; %带宽100%
sensor.frequency_response = [center_freq, bandwidth]; %定义传感器完整参数，包括中心频率和带宽
%% create the time array
kgrid.makeTime(medium.sound_speed);
%% set the input arguements
input_args = {'PMLSize', PML_size, 'PMLInside', false, ...
'PlotPML', false, 'Smooth', false, 'DataCast', 'single'};
%% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
%% 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 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
%% add first order compensation for only recording over a half plane
p0_recon = 2 * p0_recon;
%% apply a positivity condition
p0_recon(p0_recon < 0) = 0;
Thank you very much.
Sincerely,
Mengxiao