I'm working on ultrasound reflection tomography. K-wave is used in order to generate reference data for a circular ring of sensors. There is an issue arising while trying to simulate an elliptic boundary with regular grid.

A reflective region of elliptic shape in the centre of 512x512 grid is surrounded by the ring of 1500 sensors. Each sensor emits the pulse and records the reflection from the boundary of the ellipse. Thus we get an array of 1500 signals which we combine into 2d image:

https://www.dropbox.com/s/k7o8052xogaoexi/ellipse_reflection_data_full.png?dl=0

The problem is that the boundary visible on the image is not smooth and consists of disconnected curve segments:

https://www.dropbox.com/s/j8fbvbf6zn3k3pw/ellipse_reflection_data.png?dl=0

The issue remains even if the grid of size 1024x1024 is used (~50 hours of computation on CPU for 1500 sensors). How can I get a nice continuous boundary? Are there any options for smoothing while keeping the simulation physically reasonable? The code I use is cited below.

clc; clearvars all; close all;

nSource = 1500;

N = 512-1;

speed = zeros(N,N);

density = zeros(N,N);

alpha_power = zeros(N,N);

alpha_coeff = zeros(N,N);

BonA = zeros(N,N);

S = zeros(N,N);

r_h = round(N/4);

r_w = round(N/3);

c_h = floor(N/2)+1;

c_w = floor(N/2)+1;

[IH,IW] = ndgrid(1:N,1:N);

d = ((IH-c_h)./r_h).^2 + ((IW-c_w)./r_w).^2;

ellipse = d < 1.0;

figure, imagesc(ellipse);

water = ~ellipse;

medium = [];

sensor = [];

source = [];

for i=1:nSource

speed(water) = 1509;

speed(ellipse) = 1800;

density(water) = 995;

density(ellipse) = 1700;

alpha_coeff(water) = 0.002;

alpha_coeff(ellipse) = 5.0;

alpha_power(water) = 1.2;

alpha_power(ellipse) = 1.5;

BonA(water) = 6.0;

BonA(ellipse) = 8.0;

medium.sound_speed = reshape(speed,[N N]);

medium.sound_speed_ref = 1460;

medium.density = reshape(density,[N N]);

medium.alpha_coeff = reshape(alpha_coeff,[N N]);

medium.alpha_power = 1.5;

medium.BonA = reshape(BonA,[N N]);

%medium.alpha_mode = 'no_dispersion';

dx = 0.1e-3; % grid point spacing in the x direction [m]

dy = 0.1e-3; % grid point spacing in the y direction [m]

kgrid = makeGrid(N, dx, N, dy);

sensors = makeCartCircle((N/2-22)*dx,nSource);

sensor.mask = sensors;

source.p0 = cart2grid(kgrid, sensors(:,i));

[t_array, dt] = makeTime(kgrid, 2000);

% set the input arguments

kgrid.t_array = t_array;

input_args = {'PlotLayout',false, 'PMLInside', true, 'PlotScale', [0 0.05],...

'RecordMovie', true, 'MovieArgs', {'compression','None'}, 'MovieType', 'image', 'MovieName', ['~/tmp/example_movie' num2str(i)]};

% run the simulation

sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

if (i==1)

num_t_steps = numel(t_array);

data = zeros(num_t_steps,nSource);

end

data(:,i) = sensor_data(i,:);

end

figure, imagesc(abs(data)), caxis([0 0.025]);