Julia script to solve a Hamilton Jacobi Equation using the Rouy and Tourin Scheme.

hje.jl

using Pkg
Pkg.add("LinearAlgebra")
Pkg.add("GLMakie")
Pkg.add("Images")
Pkg.add("Colors")

# Requires these packages
using LinearAlgebra
using GLMakie
using Images, Colors

# Build Discretization
N = 200;
x = LinRange(-1,1,N+1)
y = LinRange(-1,1,N+1)
Δx = x[2] - x[1]
Δy = y[2] - y[1]

# Load image and convert to Grayscale
img = load("moz.png");
img = imresize(img, 201, 201)
imgg = Gray.(img);
Imatrix = imgg

# Compute functions for Perspective model
I(x,y) = Imatrix[x,y]
Q(x,y) = √(f^2/(f^2 + x^2 + y^2))
f = 2.32

# Rouy and Tourin/Godunov Scheme (Return Hamiltonian)
function H(i,j,U)
  #Iᵢⱼ = I(x[i], y[j])
  Iᵢⱼ = I(i,j)
  Iᵢⱼ = (Iᵢⱼ ≈ 0) ? 0.01 : Iᵢⱼ
  Qᵢⱼ = Q(x[i], y[j])

  D₊ˣ = (U[i+1,j] - U[i,j])/Δx
  D₋ˣ = (U[i,j] - U[i-1,j])/Δx
  D₊ʸ = (U[i,j+1] - U[i,j])/Δy
  D₋ʸ = (U[i,j] - U[i,j-1])/Δy
  Dˣ = max(0, -D₊ˣ, D₋ˣ)
  Dʸ = max(0, -D₊ʸ, D₋ʸ)

  ∇u = (Dˣ^2 + Dʸ^2)
  ∇ux = [Dˣ,Dʸ]⋅[x[i],y[j]]

  H = -exp(-2*U[i,j]) + (Iᵢⱼ*f^2/Qᵢⱼ)*√(f^2*∇u + (∇ux)^2 + Qᵢⱼ^2)
end

# Build Arrays
U = zeros(N+1,N+1)
U₁ = zeros(N+1,N+1)
U[1,:] .= 100
U[N+1,:] .= 100
U[:,1] .= 100
U[:,N+1] .= 100
U₁ .= U

# Choose Δt
Δt = 0.1*Δx
τ = 1e-5

# Solve. Should not take more than a minute ...
while 1 > 0
  for i = 2:N
    for j = 2:N
      U₁[i,j] = U[i,j] - Δt*H(i,j,U)
    end
  end
  δ = abs.(U₁[2:N,2:N] - U[2:N,2:N])
  ϵ = norm(δ,Inf)
  if(ϵ < τ)
    @show ϵ, norm(U, Inf)
    break
  end
  U .= U₁
end

# Plotting routine using Makie.jl
include("PlotSolution.jl")

# Return Figure handle
fig = PlotSolution(24);

moz.png

PlotSolution.jl

function PlotSolution(fontsize)
  fig = Figure(fontsize=fontsize, resolution=(1000,700))
  ax1 = GLMakie.Axis(fig[1, 1]);
  ax2 = GLMakie.Axis3(fig[1, 2]);
  ax3 = GLMakie.Axis3(fig[2, 1]);
  ax4 = GLMakie.Axis3(fig[2, 2]);

  image!(ax1, x, y, collect(imrotate(img,π/2)))
  GLMakie.surface!(ax2, x[2:end-1], y[2:end-1], -0.8*U[2:end-1,2:end-1],
                   colormap=(:gray,:gray),
                   lightposition = GLMakie.Vec3f0(0,0,20),
                   diffuse = GLMakie.Vec3f0(0.6,0.6,0.6),
                   ambient = GLMakie.Vec3f0(0.6,0.6,0.6),
                   specular = Vec3f0(0.1,0.1,0.1),
                   shininess=1024f0)
  GLMakie.surface!(ax3, x[2:end-1], y[2:end-1], -0.8*U[2:end-1,2:end-1],
                   colormap=(:gray,:gray),
                   lightposition = GLMakie.Vec3f0(0,0,20),
                   diffuse = GLMakie.Vec3f0(0.6,0.6,0.6),
                   ambient = GLMakie.Vec3f0(0.6,0.6,0.6),
                   specular = Vec3f0(0.1,0.1,0.1),
                   shininess=1024f0)
  GLMakie.surface!(ax4, x[2:end-1], y[2:end-1], -0.8*U[2:end-1,2:end-1],
                   colormap=(:gray,:gray),
                   lightposition = GLMakie.Vec3f0(0,0,20),
                   diffuse = GLMakie.Vec3f0(0.6,0.6,0.6),
                   ambient = GLMakie.Vec3f0(0.6,0.6,0.6),
                   specular = Vec3f0(0.1,0.1,0.1),
                   shininess=1024f0)
  ax1.title = "Intensity Image"
  ax2.title = "Reconstructed Surface"
  ax1.xlabel = ""
  ax1.ylabel = ""
  ax2.ylabel = ""
  ax2.xlabel = ""
  ax2.zlabel = ""
  ax3.ylabel = ""
  ax3.xlabel = ""
  ax3.zlabel = ""
  ax4.ylabel = ""
  ax4.xlabel = ""
  ax4.zlabel = ""
  fig
end

sol.jpeg