quick and dirty bilateral and noise aware bilateral (supposedly you should gen a LUT): Jang, Sung-Joon & Hwang, Youngbae. (2020). Noise-Aware and Light-Weight VLSI Design of Bilateral Filter for Robust and Fast Image Denoising in Mobile Systems. Sensors. 20. 4722. 10.3390/s20174722.

NABF.jl

# https://www.researchgate.net/publication/343806997_Noise-Aware_and_Light-Weight_VLSI_Design_of_Bilateral_Filter_for_Robust_and_Fast_Image_Denoising_in_Mobile_Systems
using LinearAlgebra, Images, SpecialFunctions

function fs(xi, x, σ_s)
    exp((norm(xi - x) ^ 2) / (-2 * σ_s ^ 2))
end

function fr(xi, x, σ_r)
    exp((norm(xi - x) ^ 2) / (-2 * σ_r ^ 2))
end

function bf(img, σ_s, σ_r, Ω)
    new = copy(img)
    height = size(img, 1)
    width = size(img, 2)
    Threads.@threads for h in 1:height
        for w in 1:width
            inner = [0.0, 0.0, 0.0]
            W = 0.0
            for i in -Ω:Ω
                ch = min(max(h + i, 1), height)
                for j in -Ω:Ω
                    cw = min(max(w + j, 1), width)
                    temp = fs([ch, cw], [h, w], σ_s) * fr(img[ch, cw, :], img[h, w, :], σ_r)
                    inner .+= img[ch, cw, :] * temp
                    W += temp
                end
            end
            new[h, w, :] = inner / W
        end
    end
    new
end

function μ(c0, c1, I::Vector{})
    I = norm(I)
    c0 * I + c1
end

function μ(c0, c1, I::Number)
    c0 * I + c1
end

function frNABF(xi, x, μ)
    # https://specialfunctions.juliamath.org/latest/functions_list/#SpecialFunctions.besseli
    exp(-2 * μ) * besseli(norm(xi - x),  2 * μ)
end

function nabf(img, σ_s, c::Tuple{Float64, Float64}, Ω)
    new = copy(img)
    height = size(img, 1)
    width = size(img, 2)
    Threads.@threads for h in 1:height
        for w in 1:width
            inner = [0.0, 0.0, 0.0]
            W = 0.0
            for i in -Ω:Ω
                ch = min(max(h + i, 1), height)
                for j in -Ω:Ω
                    cw = min(max(w + j, 1), width)
                    temp = fs([ch, cw], [h, w], σ_s) * frNABF(img[ch, cw, :], img[h, w, :], μ(c[1], c[2], img[ch, cw, :]))
                    inner .+= img[ch, cw, :] * temp
                    W += temp
                end
            end
            new[h, w, :] = inner / W
        end
    end
    new
end

function main()
    σ_s = 10.0
    σ_r = 2
    Ω = 5
    c = (0.4, 20.0)
    img = load("test.png")
    proper = permutedims(channelview(img), (2, 3, 1))
    bfed = bf(proper, σ_s, σ_r, Ω)
    bfed = permutedims(bfed, (3, 1, 2))
    save("test_bf.png", colorview(RGB, bfed))
    println("Bilateral done.")
    nabfed = nabf(proper, σ_s, c, Ω)
    nabfed = permutedims(nabfed, (3, 1, 2))
    save("test_nabf.png", colorview(RGB, nabfed))
    println("NABF done.")
end

main()