http://dx.doi.org/10.1109/TCE.2010.5506014 and everything before that

koreandebander.jl

# False contour detection is described in http://dx.doi.org/10.1109/ICCE.2006.1598468
# Smoothing is described in http://dx.doi.org/10.1109/TCE.2010.5506014
# Direction estimation is described in http://dx.doi.org/10.1109/TCE.2006.1706513
#
# These papers seem to all be a part of the same project.
# Their end result uses the original false contour detection and uses a NN along with the directional smoothing used here to choose which pixels to smooth and how.
# Keep in mind it works in your given image's bit depth.
#
# Usage in REPL:
# include("/path/to/script.jl")
# using ImageView
# img = load("/path/to/image.png")
# img = channelview(img)[1, :, :] # red channel - you can ycbcr to get luma instead
# cmap = candidatemap(img);
# dmap = directionmap(img, cmap);
# img = float.(img)
# smoothed = smooth(img, cmap, dmap);
# imshow(smoothed)


using Images, Statistics, ImageFiltering

"""
Re-quantize a value `v` by 2^`bits`.
"""
function requantize(v::Number, bits::Int)
    b = 2 ^ bits
    round(v / b) * b
end

function gradientkernel()
    (
     [0.0 0.0 0.0;
      -0.5 0.0 0.5;
      0.0 0.0 0.0],
     [0.0 -0.5 0.0;
      0.0 0.0 0.0;
      0.0 0.5 0.0]
    )
end

# this is a bad idea
"""
Find candidate contours and edges by comparing the input image to its re-quantized version and finding the edges of the result.

For some reason, the papers don't specify which algorithm they use to find edges here, so this just uses a simple [-1, 0, 1] gradient.
"""
function candidatemap(img::Matrix{}, bits::Int = 2)
    rq = requantize.(img, bits)
    threshold = (2 ^ bits) / 256
    canmap = @. abs(img - rq)
    canmap = imfilter(canmap, gradientkernel())
    canmap .< threshold
end

"""
Get direction index map. The actual window used is -window:window. The threshold is the maximum difference between highest and lowest directional contrast feature.

Returns a matrix of direction indices with the assumption of the following direction set:

1. horizontal
2. vertical
3. diagonal
4. anti-diagonal
5. unbiased
"""
function directionmap(img::Matrix{}, candidates::BitMatrix, window::Int = 5, threshold::Number = 16)
    H, W = size(img)
    if window % 2 == 0
        throw(DomainError("window must be uneven"))
    end
    padded = padarray(img, Pad(:replicate, window + 1, window + 1))

    # this is just a matrix of indices
    dmap = zeros(Int, H, W)
    
    # directional contrast features
    dcf = Vector{Number}(undef, 4)

    for h ∈ 1:H
        for w ∈ 1:W
            if ~(candidates[h, w])
                continue
            end

            hwin = h - window:h + window
            wwin = w - window:w + window

            # horizontal
            dcf[1] = sum([abs(padded[i, j] - padded[i + 1, j]) for i ∈ hwin for j ∈ wwin])
            # vertical
            dcf[2] = sum([abs(padded[i, j] - padded[i, j + 1]) for i ∈ hwin for j ∈ wwin])
            # diagonal
            dcf[3] = sum([abs(padded[i, j] - padded[i + 1, j - 1]) for i ∈ hwin for j ∈ wwin])
            # anti-diagonal
            dcf[4] = sum([abs(padded[i, j] - padded[i - 1, j + 1]) for i ∈ hwin for j ∈ wwin])

            hi = argmax(dcf)
            lo = argmin(dcf)

            if abs(dcf[hi] - dcf[lo]) < threshold
                # unbiased
                dmap[h, w] = 5
            else
                dmap[h, w] = hi
            end
        end
    end

    dmap
end

"""
The direction set.
"""
function directionset(i::Int)
    k = zeros(Int, 2)
    if i == 1
        k[1] = 1
    elseif i == 2
        k[2] = 1
    elseif i == 3
        k[1] = -1
        k[2] = 1
    elseif i == 4
        k[1] = 1
        k[2] = -1
    else
        throw(DomainError("Only four directions"))
    end
    k, k .* -1
end

"""
Traverse from position in direction until either maximum number of travels is hit or the hitmap returns false.
"""
function traversedirection(img, hitmap, position::Vector{Int}, d::Vector{Int}, max_travel::Int = 10, threshold::Number = 0.1)
    i = 1
    pos = position
    next = position .+ d
    while i ≤ max_travel && all(next .< size(img)) && all(next .> [1, 1]) && hitmap[next...] && abs(img[position...] .- img[next...]) < threshold
        pos = copy(next)
        next .+= d
        i += 1
    end
    img[pos...], i
end

"""
Directionally smooth contours in image according to direction. If direction is unbiased (5), use a Gaussian blur with σ.
"""
function smooth(img::Matrix{}, contours::BitMatrix, directionmap::Matrix{Int}, σ::Number = 3.0, max_travel::Int = 10, max_diff::Number = 0.1)
    H, W = size(img)
    out = copy(img)
    pad = padarray(img, Pad(:reflect, 1, 1))
    cpad = padarray(contours, Pad(:reflect, 1, 1))
    gauss = copy(img)
    gauss = imfilter(gauss, Kernel.gaussian(σ))
    v = Vector{typeof(img[1])}(undef, 2)
    dist = Vector{Int}(undef, 2)

    for h ∈ 1:H
        for w ∈ 1:W
            if contours[h, w]
                di = directionmap[h, w]
                if di == 5
                    out[h, w] = gauss[h, w]
                    continue
                end
                ds = directionset(di)
                for i in 1:2
                    v[i], dist[i] = traversedirection(pad, cpad, [h, w], ds[i], max_travel, max_diff)
                end
                t = sum(dist)
                out[h, w] = sum([v[i] * (1 - dist[i] / t) for i ∈ 1:2])
            end
        end
    end

    out
end