Skip to content

Instantly share code, notes, and snippets.

@c42f
Last active August 12, 2020 14:33
Show Gist options
  • Save c42f/454c1cbaed66354b77912a1cf41a0d96 to your computer and use it in GitHub Desktop.
Save c42f/454c1cbaed66354b77912a1cf41a0d96 to your computer and use it in GitHub Desktop.
Makie.jl based visualization of 3D 1/f noise
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)
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment