anisotropic_diffusion.jl

using Plots, Statistics, ImageFiltering

##

N = 200
img = zeros(N,N) 
img[80:120,80:120] .= 1
img = img .+ rand(N,N)*0.5

heatmap(img)

grad = ImageFiltering.imgradients(img, KernelFactors.bickley, "replicate")
G = @. sqrt(grad[1]^2 + grad[2]^2)
gradient_map = exp.(-10*G)

plot(
    heatmap(img, clim=(0,1)),
    heatmap(G, clim=(0,1)),
    heatmap(gradient_map, clim=(0,1)),
)

##


# based on : https://www.youtube.com/watch?v=QoN6I_5yOFs

diffuse!(img, speed; Δt = 0.1) = begin

    N = size(img,1)

    for i in 2:N-1, j in 2:N-1
        gradient_up    = img[i+1,j] - img[i,j]
        gradient_down  = img[i-1,j] - img[i,j]
        gradient_right = img[i,j+1] - img[i,j]
        gradient_left  = img[i,j-1] - img[i,j]

        speed_up    = speed[i+1,j]
        speed_down  = speed[i-1,j]
        speed_right = speed[i,j+1]
        speed_left  = speed[i,j-1]

        img[i,j] += Δt*(
            gradient_up * speed_up +
            gradient_down * speed_down +
            gradient_right * speed_right +
            gradient_left * speed_left 
        )
    end

end

out = copy(img)
speed = gradient_map

for i in 1:100
    diffuse!(out, speed)
end

iso = copy(img)
speed = ones(size(img))/2

for i in 1:100
    diffuse!(iso, speed)
end

plot(
    heatmap(img, legend = false, clim=(0,1), title="Original"),
    heatmap(iso, legend = false, clim=(0,1), title="Isotropic"),
    heatmap(out, legend = false, clim=(0,1), title="Anisotropic"),
    size = (1400,400),
    layout = @layout [a b c]
)

## try to reverse blurred one

out = copy(iso)
speed = ones(size(img))/2

for i in 1:100
    diffuse!(out, speed; Δt = -0.1)
end

plot(
    heatmap(img, legend = false, clim=(0,1), title="Original"),
    heatmap(iso, legend = false, clim=(0,1), title="Isotropic"),
    heatmap(out, legend = false, clim=(0,1), title="Reversed"),
    size = (1400,400),
    layout = @layout [a b c]
)

##