Making a surface from scattered data

I would to transform scattered data into a surface.

Google told me that is I could do it in Matlab[1,2] or Python [3] using Delaunay triangulation. But I cannot find a similar example in Julia with VoronoiDelaunay.jl and its related pages.

Does anybody has an example ?

Thank you :smiley:

[1] https://www.mathworks.com/matlabcentral/fileexchange/5105-making-surface-plots-from-scatter-data
[2] https://www.mathworks.com/help/matlab/ref/delaunay.html
[3] https://fabrizioguerrieri.com/blog/2017/9/7/surface-graphs-with-irregular-dataset

You can use a Thin Plate Spline / PolyharmonicSpline to do it.

Here is an old gist that has an Julia v0.5ish implementation and example and here is an updated version that works with Julia v1.0

1 Like

just for anyone that needs an example with Julia v1.x:

ENV["MPLBACKEND"]="tkagg" # if you need
using PyPlot
using LinearAlgebra

include("polyharmonic_spline.jl") 
# contain all code from https://github.com/lstagner/OrbitTomography.jl/blob/master/src/polyharmonic.jl
# and `interpolate()`definition from https://gist.github.com/lstagner/04a05b120e0be7de9915


x,y = randn(500),randn(500)
z = exp.(-(x.^2 .+ y.^2))
S2 = PolyharmonicSpline(2,[x y],z)

n=20
xgrid = ones(n)*range(-3,stop=3,length=n)'
ygrid = range(-3,stop=3,length=n)*ones(n)'

xx = reshape(xgrid,n*n)
yy = reshape(ygrid,n*n)

zz = interpolate(S2,xx,yy)
zgrid = reshape(zz,n,n);

plot_surface(xgrid,ygrid,zgrid,alpha=0.5)
scatter3D(x,y,z,color="r")
show()

Figure_1

1 Like

With the updated version you don’t need the interpolate method. The constructor for the spline returns a interpolating “function”. You can simplify the code to be

using PyPlot
using LinearAlgebra

include("polyharmonic_spline.jl") 
# contain all code from https://github.com/lstagner/OrbitTomography.jl/blob/master/src/polyharmonic.jl


x,y = randn(500),randn(500)
z = exp.(-(x.^2 .+ y.^2))
S2 = PolyharmonicSpline(2,[x y],z)

n=20
xgrid = ones(n)*range(-3,stop=3,length=n)'
ygrid = range(-3,stop=3,length=n)*ones(n)'

zgrid = S2.(xgrid,ygrid) #Take advantage of dot broadcasting

plot_surface(xgrid,ygrid,zgrid,alpha=0.5)
scatter3D(x,y,z,color="r")
show()
2 Likes