The objective is to evaluate a method for Level Set to Signed Distance Functions using the Intepolations.jl package. This serves several purposes (non-exhaustive).
In [5]:
using Interpolations
using BenchmarkTools
In [6]:
circle(x,y,r) = sqrt(sum((x,y).^2)) - r;
In [10]:
n = 2000; x = LinRange(-1,1,n); y = x; @show length(x);
In [12]:
@benchmark ls = [circle(x[i], y[j], 0.5) for i = 1:n, j = 1:n]
Out[12]:
In [127]:
using Plots
contour(x,y,ls)
Out[127]:
In [128]:
itp_lin = interpolate((x,y), ls, Gridded(Linear()));
In [129]:
i1, t1 = itp_lin(0,0), circle(0,0,0.5);
@show i1 - t1;
i2, t2 = itp_lin(0.5,0), circle(0.5,0,0.5);
@show i2 - t2;
In [130]:
itp_con = interpolate((x,y), ls, Gridded(Constant()));
In [131]:
i1, t1 = itp_con(0,0), circle(0,0,0.5);
@show i1 - t1;
i2, t2 = itp_con(0.5,0), circle(0.5,0,0.5);
@show i2 - t2;
In [132]:
using StatsBase
In [133]:
n_samples = 100; sample_pts = [(rand(),rand()) for i in 1:n_samples];
In [134]:
sample_itp_con = map(x -> itp_con(x...), sample_pts);
sample_itp_lin = map(x -> itp_lin(x...), sample_pts);
sample_true = map(x -> circle(x..., 0.5), sample_pts);
error_itp_con = sample_true .- sample_itp_con;
error_itp_lin = sample_true .- sample_itp_lin;
In [135]:
hist_con1 = histogram(error_itp_con, title="Error (20x20, Constant)")
Out[135]:
In [136]:
describe(error_itp_con)
In [137]:
hist_lin1 = histogram(error_itp_lin, title="Error (20x20, Linear)")
Out[137]:
In [138]:
describe(error_itp_lin)
In [139]:
# Performance
using BenchmarkTools
@benchmark map(x -> itp_con(x...), $sample_pts)
Out[139]:
In [140]:
@benchmark map(x -> itp_lin(x...), $sample_pts)
Out[140]:
In [141]:
@benchmark map(x -> circle(x..., 0.5), $sample_pts)
Out[141]:
In [142]:
n = 500; x = collect(LinRange(-1,1,n)); y = x; @show length(x);
In [143]:
ls = [circle(x[i], y[j], 0.5) for i = 1:n, j = 1:n]
Out[143]:
In [144]:
itp_lin_500 = interpolate((x,y), ls, Gridded(Linear()));
In [145]:
# Same sample points
sample_itp_lin500 = map(x -> itp_lin_500(x...), sample_pts);
error_itp_lin500 = sample_true .- sample_itp_lin500;
In [146]:
hist_lin2 = histogram(error_itp_lin500, title="Error (500x500, Linear)")
plot(hist_lin2, hist_lin1)
Out[146]:
In [147]:
describe(error_itp_lin500)
In [148]:
describe(error_itp_lin)
In [149]:
@benchmark map(x -> itp_lin_500(x...), $sample_pts)
Out[149]:
In [150]:
@benchmark map(x -> itp_lin(x...), $sample_pts)
Out[150]:
Observations: The runtime scales minimally with the size of the level set
In [151]:
n = 200; x = collect(LinRange(-1,1,n)); y = x; z = x; @show length(x);
In [152]:
sphere(x,y,z,r) = sqrt(sum((x,y,z).^2)) - r;
In [153]:
@time ls_sphere = [sphere(x[i], y[j], z[k], 0.5) for i = 1:n, j = 1:n, k = 1:n];
In [154]:
size(ls_sphere)
Out[154]:
In [155]:
itp_sphere = interpolate((x,y,z), ls_sphere, Gridded(Linear())); # generate interpolation
In [156]:
n_samples = 100; sample_pts = [(rand(),rand(), rand()) for i in 1:n_samples];
In [157]:
sample_itp_sphere = map(x -> itp_sphere(x...), sample_pts);
sample_true = map(x -> sphere(x..., 0.5), sample_pts);
error_sphere = sample_itp_sphere .- sample_true;
describe(error_sphere)
In [158]:
@benchmark map(x -> itp_sphere(x...), $sample_pts)
Out[158]:
In [159]:
@benchmark map(x -> sphere(x..., 0.5), $sample_pts)
Out[159]:
Observation: with 250,000 points (2D) we had a median runtime of ~15 micro seconds and with 8,000,000 points (3D) we had a median runtime of ~20 micro seconds, so scaling is quite good. Similarly for the SDF of the circle the median runtime was ~250 nano seconds, and the sphere was ~590 nanoseconds.
The Interpolations.jl package provides a convenient way to convert a Level Set to an interpolation of a Signed Distance Function. This will be highly useful in DistMesh, and should provide an immediate path to supporting Level Set and Medical Data. There is a slight overhead compared to a direct sampling of the function, but for complex geometries this is likely to balance out. A good additional test case for this would be something involving trignometric functions.
In [ ]:
In [ ]:
In [ ]: