Hi there,
First of all many congratulations on this very helpful toolbox.
I am working on travel-time ultrasound tomography. I am using frequency domain kernels (banana-doughnut kernels) to reconstruct sound speed profiles. I am building 2D simulations here to validate my experimental data and eventually to optimize them. My target is to create simulations data which could provide better reconstructions than the experimental ones (using the same system matrix - in a back-projection concept). At first, I used point-sensors but the results weren't that good. Now I am trying to give dimensions to my sensors (defining every sensor as a number of points(nodes of the grid)). The problem that I am facing is that I cannot completely understand in which order the sensor_data structure is created. As I am using all the massive data (coming from the point-sensors) to average them in order to create an averaged signal for every actual sensor.
For instance, in the attached example the amount of sensor_data is 368, while I have to link them to 32 actual sensors.
Can you help me with?
I put my code here:
%----%
clear all;
close all;
clc;
%% GRID - Centre Frequency
figure();
% computational grid
Nx = 256; % number of grid points in the x (row) direction
Ny = 256; % number of grid points in the y (column) direction
spacing=370/Nx;
dx = spacing*0.001; % grid point spacing in the x direction [m]
dy = spacing*0.001; % grid point spacing in the y direction [m]
pxl_mm=dx*1000;
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%% Define a centered circular sensor
sensor_radius = (329/2)*0.001; % [m]
num_sensor_points = 32;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);
% Sensors (grid)
a=(sensor.mask./dx)+(Nx/2);
%% Medium Porperties
% DEFINE MEDIUM
% WATER
medium.sound_speed = 1480 * ones(Nx, Ny); % [m/s]
medium.density = 9980 * ones(Nx, Ny); % [kg/m^3]
% INCLUSION
% Small_Inclusion
inc_centre_x=min(min(a))+(225/pxl_mm);
inc_centre_y=Nx-(min(min(a))+(275/pxl_mm));
radius=(14/pxl_mm);
%
plot(inc_centre_x,inc_centre_y,'*black');hold on;
[columnsInImage rowsInImage] = meshgrid(1:Nx, 1:Nx);
circlePixels = (rowsInImage - inc_centre_y).^2 ...
+ (columnsInImage - inc_centre_x).^2 <= radius.^2;
[posx,posy]=find((circlePixels)~=0);
%
% Air bottle
for i=1:length(posx)
medium.sound_speed(posx(i), posy(i)) =343; %343;
medium.density(posx(i), posy(i)) =18; %1.2;
end
%% Sources
mask=zeros(Nx);
for i=1:32
masks{i}=zeros(Nx);
end
mask3=zeros(Nx);
for i=1:32
source_x(i)=sensor.mask(2,i);
source_y(i)=sensor.mask(1,i);
%
disc_x_pos =(floor(source_x(i)/dx)+(Nx/2)); % [grid points]
disc_y_pos =(floor(source_y(i)/dy)+(Nx/2)); % [grid points]
%
masks{i}(disc_x_pos,disc_y_pos)=1;
%
arc_pos = [disc_x_pos, disc_y_pos]; % [grid points]
radius = 25; % [grid points]
diameter = 13; % [grid points]
focus_pos = [Nx/2, Nx/2]; % [grid points]
% source.p_mask = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos);
masks{i} = makeArc([Nx, Ny], arc_pos, radius, diameter, focus_pos);
%
sensor_node_num{i}=sum(masks{i}(:) == 1);
%
mask3=mask3+masks{i};
mask1=mask3;
mask1(mask1==1)=-1200;
end
%% Display
mask1=mask1+medium.sound_speed;
mask2=mask1;
mask2(mask2==280)=1480;
sensor1.mask=mask3;
medium.sound_speed=mask2;
%% Displays
imagesc(mask2);colorbar; view(0,-90); title('Domain');
figure;
imagesc(mask1); hold on; plot(a(1,:),a(2,:),'*r'); title('Geometry');
figure;
imagesc(mask3); hold on; plot(a(1,:),a(2,:),'*r'); title('Sensors');
%%
pause; figure();
% --- SIMULATION ---
% Running the Tomographic scenario
for i=1:length(sensor.mask)
% Sources
source_x(i)=sensor.mask(2,i);
source_y(i)=sensor.mask(1,i);
source.p_mask=masks{i};
% Source_Excitation_ToneBurst
tone_burst_freq=480e3;
sampling_freq=40 *tone_burst_freq;
tone_burst_cycles=3;
k = toneBurst(sampling_freq,tone_burst_freq,tone_burst_cycles,'Plot', true);
source.p=30*k;
figure();
% source.p_mode='dirichlet';
% run the simulation
kgrid.t_array = makeTime(kgrid, 1480, 0.3, 0.00042 ); %recording time 0.5 msec (0.5e-3) - auto t_end=0.005 (increasing t_end one decreases the recording time)
sensor_data{i} = kspaceFirstOrder2D(kgrid, medium, source, sensor1,'PMLInside',true);
%
% plot the simulated sensor data
figure;
imagesc(sensor_data{i}, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;
end
%----%