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.

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()

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()

I’m not sure why we need spline interpolation to do this. Using the default GR backend, seems like surface naturally interpolates via Delaunay triangulation. The following code

using Plots
x = randn(1000)
y = randn(1000)
z = @. sin(x)*sin(y)
surface(x,y,z,camera=(0,90))