Makie.jl based visualization of 3D 1/f noise

fourier_isosurface.jl

using Makie
using Meshing
using GeometryTypes
using ColorSchemes
using FFTW
using LinearAlgebra
using Random
using Interpolations


# Fourier synthesis of an isosurface using 1/f^α noise
function one_on_f_noise(T, α, f)
    # 2α needed to convert powers to amplitudes here
    ϵ = sqrt(eps(real(T)))
    A = randn(T)
    return convert(T, A / (ϵ + norm(f))^2α)
end

N = 200
α = 1.2f0
Random.seed!(2)
# symmetric frequency space
fx = range(-1,1,length=N)
fy = reshape(fx,1,:)
fz = reshape(fx,1,1,:)
@time spectrum = one_on_f_noise.(ComplexF32, α, Vec3f0.(fx, fy, fz))
@time iso = real.(fft(fftshift(spectrum)))

# Compute an isosurface via Meshing.jl
#
# There's various choices of algorithm available.  One is to use marching cubes
# which seems to be quite fast for large volumes, but results in duplicate
# vertices which means that topology will be disconnected and any directly
# estimated surface normals will not be smooth.
#   mc_algo = MarchingCubes(iso=0, insidepositive=false)
# Another nice option is NaiveSurfaceNets which are still quite fast and have
# watertight topology.
mc_algo = NaiveSurfaceNets(iso=0, insidepositive=true)
@time m = HomogenousMesh{Point{3,Float32},Face{3,Int}}(iso, mc_algo)

# We can visualize the mesh `m` directly, but this will compute surface normals
# from the triangulation which can introduce discretization artifacts.
#
# To get more accurate normals, we can instead estimate the gradient of `iso`
# at each vertex using Interpolations.jl. (Note: `range` may not be entirely
# correct here)
@time itp = scale(interpolate(iso, BSpline(Quadratic(Periodic(OnGrid())))),
            range(-1,1,length=size(iso,1)),
            range(-1,1,length=size(iso,2)),
            range(-1,1,length=size(iso,3)))
@time normals = [normalize(Vec3f0(Interpolations.gradient(itp, Tuple(v)...))) for v in m.vertices]
# Construct a new mesh with these more accurate normals
m2 = HomogenousMesh(vertices=m.vertices, faces=m.faces, normals=normals)

# Plot with Makie
mesh(m2, color=last.(m.vertices), colormap=ColorSchemes.turbo.colors)

fourier_isosurface.png

iso_interactive.jl

using Makie
using Meshing
using GeometryTypes
using ColorSchemes
using FFTW
using LinearAlgebra
using Random
using Interpolations

# Interactive attempt with varying α
sl,α = textslider(0.5:0.1:2, "alpha", start=1)

m = lift(α) do α
    N = 50
    r = range(-1,1,length=N)
    f = [Vec3f0(x,y,z) for x in r, y in r, z in r]

    Random.seed!(2)
    spectrum = randn(ComplexF64, size(f)) ./ (1e-8 .+ norm.(f)).^2α
    iso = real.(fft(fftshift(spectrum)))
    # Isosurface via Meshing
    mc_algo = MarchingCubes(iso=0, insidepositive=true)
    m = HomogenousMesh{Point{3,Float32},Face{3,Int}}(iso, mc_algo)
    HomogenousMesh(vertices=m.vertices, faces=m.faces, color=last.(m.vertices))
end

# Plot using Makie.
# Color seems to be ignored??
pl = mesh(m, colormap=ColorSchemes.turbo.colors)

# Attempting the following explodes
# pl = mesh(m, color=lift(m->last.(m.vertices),m), colormap=ColorSchemes.turbo.colors)

hbox(sl, pl)