Smoothing surface in GLMakie

[Alicanter]

I have data from a simulation (test.csv.gz) and I am trying to plot it as a surface.
The equations have a singularity at a certain threshold, so there is a cutoff in the values of the surface.
The problem I am facing is that as the values are irregular near the threshold, there is a “sawtooth” at the top instead of being smooth.

Is there a way to change the algorithm used to compute the surface so it is smoother? Or some way to interpolate it before?

I could always compute more points, but I had similar problems in other occasions so I’d like to find a solution if possible.

Here is the script to reproduce it:

using CSV
using DataFrames
using GLMakie

df = DataFrame(CSV.File("test.csv.gz"))

fig = Figure()


ax = Axis3(fig[1, 1]; aspect = :equal, xlabel = "x", ylabel = "y", zlabel = "z")
cmap = categorical_colors([:blue], 1)
surface!(ax, df[!, :x], df[!, :y], df[!, :z], colormap = cmap)
scatter!(ax, df[!, :x], df[!, :y], df[!, :z], color = :black, markersize = 3)
save("test.png", fig, px_per_unit = 5)

Thanks!

[sdanisch]

Maybe use filters from ImageFiltering.jl · ImageFiltering?
They should work with generic arrays…

[Alicanter]

I would prefer not smoothing by changing the values and instead computing a surface that is like a convex hull of the points, or adding points so the surface looks better.

[Nathan_Boyer]

Regardless of how you smooth the data, you don’t have to change the original data. Create a separate df_smooth object and plot that beside df until you are happy with the look of df_smooth.

[rocco_sprmnt21]

I don’t understand what your end goal is. But perhaps, as a naive idea, the following scheme can give you some ideas.

julia> m,M= extrema(df.z)
(0.0, 0.1087877834426117)
# 80%M
df008=findall( >(0.08),df.z)
df008m1=df008.-1
df8m=(dfsmooth.z[df008].+dfsmooth.z[df008m1])/2

dfsmooth=copy(df)
dfsmooth.z[df008].=df8m
dfsmooth.z[df008n].=df8m

fig07 = Figure()

ax = Axis3(fig07[1, 1]; aspect = :equal, xlabel = "x", ylabel = "y", zlabel = "z")
cmap = categorical_colors([:blue], 1)
surface!(ax, dfsmooth[!, :x], dfsmooth[!, :y], dfsmooth[!, :z], colormap = cmap)
scatter!(ax, dfsmooth[!, :x], dfsmooth[!, :y], dfsmooth[!, :z], color = :black, markersize = 3)
save("test.png", fig, px_per_unit = 5)

[rocco_sprmnt21]

A different naive idea could be to smooth the saw tips between the different minima of the projections on the X plane (for a part) or on the Y plane (for a part) or perhaps it would be better (for this figure) to use a projection on a cylinder.
placed the part of the X plane as an illustrative example of the idea.

insertcols!(df,:nr=>1:nrow(df))
gdf=groupby(df,:x)
cdf=combine(gdf, [:y,:z,:nr]=>((y,z,n)->[(y[argmax(z)],z[argmax(z)],n[argmax(z)])])=>identity)



tcdf=transform(cdf,:z=>(z->[0;diff(z)])=>:zd)


sinks=findall(<(-0.01),tcdf.zd)

julia> for (i,s) in enumerate(sinks)
           #cdf[tcdf.nr[s:s[i+1]]].=tcdf.z[s[i+1]]
           tcdf.z[s:sinks[i+1]].=tcdf.z[sinks[i+1]]
       end
ERROR: BoundsError: attempt to access 14-element Vector{Int64} at index [15]
Stacktrace:
# error but makes the wanted  transformation 
# It's not worth fixing this code just to exemplify the idea

dfsmooth=copy(df)
dfsmooth.z[tcdf.nr].=tcdf.z

figsm = Figure()

ax = Axis3(figsm[1, 1]; aspect = :equal, xlabel = "x", ylabel = "y", zlabel = "z")
cmap = categorical_colors([:blue], 1)
surface!(ax, dfsmooth[!, :x], dfsmooth[!, :y], dfsmooth[!, :z], colormap = cmap)
scatter!(ax, dfsmooth[!, :x], dfsmooth[!, :y], dfsmooth[!, :z], color = :black, markersize = 3)
save("byX_smooth.png", figsm, px_per_unit = 5)

to “cut” the tips projected onto a cylinder in a more graceful way you could use some properly done interpolation technique
Could the projection on a cylinder be simulated via a groupby done on the x/y ratio?

PS
Obviously, all this while waiting for someone who knows how to seriously do what the OP requests.

[rocco_sprmnt21]

Another smoothing exercise

smoother

using CSV,DataFrames,Plots
df = DataFrame(CSV.File("test.csv.gz"))
insertcols!(df,:nr=>1:nrow(df))
gdf=groupby(df,:x) # proiezione sul piano X(x,0,0)
cdf=combine(gdf, [:y,:z,:nr]=>((y,z,n)->[(y[argmax(z)],z[argmax(z)],n[argmax(z)])])=>identity) # prende i valori di cresta

function findmin(cdf)
    insertcols!(cdf,1,:r=>1:nrow(cdf))
    insertcols!(cdf,4, :dz=> [0;diff(cdf.z)]) # serve a vedere dove "sale" e dove "scende"
    cdf.min.=true
    st=true
    for i in eachindex(cdf.dz) # cerca i punti di minimo da interpolare per tagliare via le creste
        if st 
            cdf.min[i]= true
            if cdf.dz[i] < 0 
                st=false
            end
        else 
            cdf.min[i]= cdf.dz[i] < 0 ? true :  false
        end
    end
end
function smoother(cdf)
    cdf.csm=cumsum(cdf.min)
    gcdf=groupby(cdf,:csm)

    function interpol(gr) #calcola i valori lineari tra due mimimi consecutivi e li sostituisce ai valori di "cresta"
        n=nrow(gr)+1
        zf=cdf.z[gr.r[n-1]+1]
        zi=cdf.z[gr.r[1]]
        for i in 2:n-1
            gr.z[i]=zi+(zf-zi)*i/n
        end
        gr
    end
    dfsm=combine(gr->nrow(gr)>1 ? interpol(gr) : gr, gcdf)
    Plots.plot!(dfsm.x, dfsm.z,xlims=(0, 10), ylims=(0,0.12))
end

using Plots
Plots.plot(cdf.x, cdf.z,xlims=(0, 10), ylims=(0,0.12))
ccdf=copy(cdf)
findmin(ccdf)
smoother(ccdf)

PS
Perhaps, in some cases, more passes should be made.
Like when you shave

[Alicanter]

I am not looking for a way to change the points that are computed so the surface looks better. I am looking for a way for the surface to be more connected near the top.

A maybe better and simpler example is this code:

using GLMakie

# Gaussian with a maximum at (3, 3).
f(x, y) = 10*exp(-(x - 3)^2/2 - (y - 3)^2/2)

X = ones(7)' .* (0:6)
Y = ones(7) .* (0:6)'
Z = f.(X, Y)

# Convert to array and remove middle point to simulate singularity.
# An alternative would be to set Z[4, 4] = Inf but messes up with the colors.
X = reduce(vcat, X)
Y = reduce(vcat, Y)
Z = reduce(vcat, Z)

deleteat!(X, 25)
deleteat!(Y, 25)
deleteat!(Z, 25)

fig = Figure()
ax = Axis3(fig[1, 1]; aspect = :equal)
surface!(ax, X, Y, Z)
fig

This is a 2D Gaussian with a “singularity” at the middle (maximum).
The resulting surface has 4 triangles at the top, but I’d like it to also have egdes between the top points so they are connected and the surface looks like a Gaussian clipped at some threshold instead

In this case, Plots.jl works as I expect:

using Plots
Plots.surface(X, Y, Z)