Making a surface from scattered data

[Noel_Araujo]

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

[1] Making Surface Plots From Scatter Data - File Exchange - MATLAB Central
[2] Delaunay triangulation - MATLAB delaunay
[3] Use Python to plot Surface graphs of irregular Datasets

[lstagner]

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

[Noel_Araujo]

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

[lstagner]

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

[jlchan]

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

gives something like