-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSimulation_data_blockbuster_sensors.m
131 lines (101 loc) · 4.31 KB
/
Simulation_data_blockbuster_sensors.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
clear
%%
% create the computational grid
scale = 1;
PML_size = 10; % size of the PML in grid points
Nx = 200 * scale - 2 * PML_size; % number of grid points in the x direction
Ny = 200 * scale - 2 * PML_size; % number of grid points in the y direction
Nz = 200 * scale - 2 * PML_size; % number of grid points in the z direction
dx = 1e-3 / scale; % grid point spacing in the x direction [m]
dy = 1e-3 / scale; % grid point spacing in the y direction [m]
dz = 1e-3 / scale; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium
medium.sound_speed = 1500; %[m/s]
% create initial pressure distribution using makeDisc
disc_radius = 3; % [grid points]
disc_magnitude = 10; %[Pa]
p0 = disc_magnitude * makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, disc_radius);
% smooth the initial pressure distribution and restore the magnitude
source.p0 = smooth(p0, true);
%% sensor mask creation
sensor.mask = zeros(kgrid.Nx,kgrid.Ny, kgrid.Nz);
sensor.mask (:,:,141:180) = 1;
%% visualisation
voxelPlot(double(p0 | sensor.mask));
%% run simulation
% create the time array
kgrid.makeTime(medium.sound_speed);
% input arguments **inclue PML
input_args = {'PMLSize', PML_size, 'PMLInside', false, ...
'PlotPML', true, 'Smooth', true, 'DataCast', 'single'};
% run simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
%% preparing for reconstruction
% creation of reshape sensor array for FFT
sensor_data_rs_pre = zeros(180,30,180);
sensor_data_rs_pre(:,:,180) = 1;
sensor_data_rs = reshape(a, Nx, Ny, kgrid.Nt);
%% reconstruction simulation in FFT
% reconstruct the initial pressure
p_xyz = kspacePlaneRecon(sensor_data_rs, kgrid.dx, kgrid.dy, kgrid.dt, ...
medium.sound_speed, 'DataOrder', 'yzt', 'PosCond', true, 'Plot', true);
%% pre-plots
% define a k-space grid using the dimensions of p_xyz
[Nx_recon, Ny_recon, Nz_recon] = size(p_xyz);
kgrid_recon = kWaveGrid(Nx_recon, kgrid.dt * medium.sound_speed, Ny_recon, kgrid.dy, Nz_recon, kgrid.dz);
% define a k-space grid with the same z-spacing as p0
[Nx_p0, Ny_p0, Nz_p0] = size(source.p0);
kgrid_interp = kWaveGrid(Nx_p0, kgrid.dx, Ny_p0, kgrid.dy, Nz_p0, kgrid.dz);
% resample the p_xyz to be the same size as p0; for a matrix indexed as
% [M, N, P], the axis variables passed to interp3 must be given in the
% order N, M, P
p_xyz_rs = interp3(kgrid_recon.y - min(kgrid_recon.y(:)), ...
kgrid_recon.x - min(kgrid_recon.x(:)), ...
kgrid_recon.z - min(kgrid_recon.z(:)), ...
p_xyz, ...
kgrid_interp.y - min(kgrid_interp.y(:)), ...
kgrid_interp.x - min(kgrid_interp.x(:)), ...
kgrid_interp.z - min(kgrid_interp.z(:)));
%% plots
% plot the initial pressure
voxelPlot(double(p0 | sensor.mask));
set(gca, 'Projection', 'perspective');
view([0, 99]);
% plot the initial pressure
figure;
plot_scale = [-10, 10];
subplot(2, 2, 1);
imagesc(kgrid_interp.y_vec * 1e3, kgrid_interp.x_vec * 1e3, squeeze(source.p0(:, :, Nz/2)), plot_scale);
title('x-y plane');
axis image;
subplot(2, 2, 2);
imagesc(kgrid_interp.z_vec * 1e3, kgrid_interp.x_vec * 1e3, squeeze(source.p0(:, Ny/2, :)), plot_scale);
title('x-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3);
imagesc(kgrid_interp.z_vec * 1e3, kgrid_interp.y_vec * 1e3, squeeze(source.p0(Nx/2, :, :)), plot_scale);
title('y-z plane');
axis image;
colormap(getColorMap);
sgtitle('Initial pressure source');
% plot the reconstructed initial pressure
figure;
subplot(2, 2, 1);
imagesc(kgrid_interp.y_vec * 1e3, kgrid_interp.x_vec * 1e3, squeeze(p_xyz_rs(:, :, Nz/2)), plot_scale);
title('x-y plane');
axis image;
subplot(2, 2, 2);
imagesc(kgrid_interp.z_vec * 1e3, kgrid_interp.x_vec * 1e3, squeeze(p_xyz_rs(:, Ny/2, :)), plot_scale);
title('x-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3);
imagesc(kgrid_interp.z_vec * 1e3, kgrid_interp.y_vec * 1e3, squeeze(p_xyz_rs(Nx/2, :, :)), plot_scale);
title('y-z plane');
axis image;
colormap(getColorMap);
sgtitle('Reconstructed Pressure Source');
%% view the reconstruction slice by slice
flyThrough(p_xyz_rs);