Seam Carving

begin
    using PlutoUI, WasmMakie
end

begin
    # ── helpers (row-major: index (i-1)*w + j, row 1 = top) ─────────────────────
    # clamped pixel access, so gradients at the border are well-defined
    function pat(img::Vector{Float64}, nr::Int, w::Int, i::Int, j::Int)
        ii = i; jj = j
        ii < 1 && (ii = 1); ii > nr && (ii = nr)
        jj < 1 && (jj = 1); jj > w && (jj = w)
        return img[(ii - 1) * w + jj]
    end

    # render a row-major image as a WasmMakie grayscale figure (convert to the
    # column-major, row-1-on-top layout that image! wants)
    function show_gray(img::Vector{Float64}, nr::Int, w::Int; px::Int = 300)
        disp = Vector{NTuple{4,Float64}}(undef, nr * w)
        for i in 1:nr
            for j in 1:w
                v = img[(i - 1) * w + j]
                disp[j + (nr - i) * w] = (v, v, v, 1.0)
            end
        end
        fig = Figure(size = (px, max(40, round(Int, px * nr / w))))
        ax = Axis(fig[1, 1])
        hidedecorations!(ax); hidespines!(ax)
        image!(ax, (0.0, Float64(w)), (0.0, Float64(nr)), disp,
               Int64(w), Int64(nr); interpolate = false)
        fig
    end
end

# A synthetic scene: a soft left-to-right background gradient (LOW energy) with a bright
# disk (HIGH-energy edges) a little left of centre. Seam carving should eat the bland
# background and leave the disk intact.
function scene(nr::Int, nc::Int)
    img = Vector{Float64}(undef, nr * nc)
    ci = nr * 0.5
    cj = nc * 0.42
    rad = nr * 0.30
    for i in 1:nr
        for j in 1:nc
            v = 0.25 + 0.40 * (j / nc)
            di = i - ci
            dj = j - cj
            if di * di + dj * dj < rad * rad
                v = 0.95
            end
            img[(i - 1) * nc + j] = v
        end
    end
    return img
end

# Energy = magnitude of the image gradient (how fast brightness changes). Flat areas are
# cheap to cut; edges are expensive, so seams avoid them.
function energy(img::Vector{Float64}, nr::Int, w::Int)
    e = Vector{Float64}(undef, nr * w)
    for i in 1:nr
        for j in 1:w
            gx = pat(img, nr, w, i, j + 1) - pat(img, nr, w, i, j - 1)
            gy = pat(img, nr, w, i + 1, j) - pat(img, nr, w, i - 1, j)
            e[(i - 1) * w + j] = sqrt(gx * gx + gy * gy)
        end
    end
    return e
end

# Find the minimum-energy vertical seam with dynamic programming: dp[i,j] = e[i,j] +
# cheapest of the three cells above. Then trace the cheapest bottom cell upward.
# Returns the column of the seam at each row.
function min_seam(e::Vector{Float64}, nr::Int, w::Int)
    dp = copy(e)
    for i in 2:nr
        for j in 1:w
            best = dp[(i - 2) * w + j]
            if j > 1 && dp[(i - 2) * w + j - 1] < best
                best = dp[(i - 2) * w + j - 1]
            end
            if j < w && dp[(i - 2) * w + j + 1] < best
                best = dp[(i - 2) * w + j + 1]
            end
            dp[(i - 1) * w + j] = e[(i - 1) * w + j] + best
        end
    end
    seam = Vector{Int}(undef, nr)
    sj = 1
    for j in 2:w
        dp[(nr - 1) * w + j] < dp[(nr - 1) * w + sj] && (sj = j)
    end
    seam[nr] = sj
    for i in nr:-1:2
        j = seam[i]
        nj = j
        best = dp[(i - 2) * w + j]
        if j > 1 && dp[(i - 2) * w + j - 1] < best
            best = dp[(i - 2) * w + j - 1]
            nj = j - 1
        end
        if j < w && dp[(i - 2) * w + j + 1] < best
            nj = j + 1
        end
        seam[i - 1] = nj
    end
    return seam
end

# Remove `k` seams one after another. Each removal recomputes the energy and the best
# seam on the now-narrower image — that re-use is what makes it dynamic programming.
function carve(img0::Vector{Float64}, nr::Int, nc::Int, k::Int)
    cur = copy(img0)
    w = nc
    steps = k
    steps > (nc - 4) && (steps = nc - 4)
    for _s in 1:steps
        e = energy(cur, nr, w)
        seam = min_seam(e, nr, w)
        nw = w - 1
        nxt = Vector{Float64}(undef, nr * nw)
        for i in 1:nr
            sj = seam[i]
            col = 0
            for j in 1:w
                if j != sj
                    col += 1
                    nxt[(i - 1) * nw + col] = cur[(i - 1) * w + j]
                end
            end
        end
        cur = nxt
        w = nw
    end
    return cur, w
end

let
    nr, nc = 48, 64
    img = scene(nr, nc)
    show_gray(img, nr, nc; px = 320)
end

let
    nr, nc = 48, 64
    img = scene(nr, nc)
    e = energy(img, nr, nc)
    # normalise energy to 0..1 for display
    hi = e[1]
    for t in 2:(nr * nc)
        e[t] > hi && (hi = e[t])
    end
    hi <= 0.0 && (hi = 1.0)
    for t in 1:(nr * nc)
        e[t] = e[t] / hi
    end
    show_gray(e, nr, nc; px = 320)
end

let
    nr, nc = 48, 64
    img = scene(nr, nc)
    carved, w = carve(img, nr, nc, nseams)
    show_gray(carved, nr, w; px = 320)
end