Рисование эпициклами

"""
Вспомогательные функции учебного скрипта fourier_epicycles.ngscript.
Порог, обход границ, соединение участков, равномерная выборка и визуализация.
ДПФ и построение цепочки явно показаны в математических ячейках скрипта.
"""
module EpicycleLesson
using Images, FileIO, FFTW, Plots
using Statistics, Printf
const Point = Tuple{Int,Int}
const Edge = Tuple{Point,Point}

function otsu_threshold(gray)
    lo, hi = extrema(gray)
    hi > lo || error("Изображение однотонное: контур не найден.")
    hist = zeros(Int, 256)
    for x in gray
        hist[clamp(floor(Int, x * 255) + 1, 1, 256)] += 1
    end
    total = length(gray)
    weighted = sum((i - 1) * hist[i] for i in 1:256)
    wb = 0; sb = 0.0; score = -Inf; threshold = 127
    for i in 1:255
        wb += hist[i]; sb += (i - 1) * hist[i]
        wf = total - wb
        (wb == 0 || wf == 0) && continue
        candidate = wb * wf * (sb / wb - (weighted - sb) / wf)^2
        if candidate > score
            score = candidate; threshold = i - 1
        end
    end
    return threshold / 255
end

function remove_small_components!(mask, min_pixels)
    h, w = size(mask)
    seen = falses(h, w)
    for c in 1:w, r in 1:h
        (!mask[r,c] || seen[r,c]) && continue
        component = Point[(r,c)]; seen[r,c] = true; head = 1
        while head <= length(component)
            rr, cc = component[head]; head += 1
            for dr in -1:1, dc in -1:1
                r2, c2 = rr + dr, cc + dc
                if 1 <= r2 <= h && 1 <= c2 <= w && mask[r2,c2] && !seen[r2,c2]
                    seen[r2,c2] = true; push!(component, (r2,c2))
                end
            end
        end
        if length(component) < min_pixels
            for p in component; mask[p...] = false; end
        end
    end
    return mask
end

function load_mask(input_path; foreground=:auto, threshold=nothing,
                   max_dimension=800, min_component_pixels=10)
    lowercase(splitext(input_path)[2]) in (".png", ".jpg", ".jpeg", ".bmp") ||
        error("Поддерживаются PNG, JPEG и BMP.")
    filesize(input_path) <= 10 * 1024^2 || error("Файл больше 10 МиБ. Уменьшите изображение.")
    img = FileIO.load(input_path)
    ndims(img) == 2 && eltype(img) <: Colorant || error("Нужно одно растровое изображение.")
    length(img) <= 12_000_000 || error("Изображение больше 12 миллионов пикселей.")
    max_dimension = min(max_dimension, 800)
    original_size = size(img)
    h, w = original_size
    if max(h,w) > max_dimension
        ratio = max_dimension / max(h,w)
        rows = round.(Int, range(1, h; length=max(2,round(Int,h*ratio))))
        cols = round.(Int, range(1, w; length=max(2,round(Int,w*ratio))))
        img = img[rows,cols]
    end
    gray = Float64.(Gray.(img))
    opacity = Float64.(alpha.(img))
    has_transparency = minimum(opacity) < 0.05 && maximum(opacity) > 0.5
    mode = foreground
    if mode == :alpha || (mode == :auto && has_transparency)
        mode = :alpha
        used_threshold = threshold === nothing ? 0.5 : Float64(threshold)
        mask = BitMatrix(opacity .>= used_threshold)
    else
        # Полупрозрачные тёмные линии компонуются на белом фоне.
        gray = gray .* opacity .+ (1 .- opacity)
        if mode == :auto
            border = vcat(gray[1,:],gray[end,:],gray[:,1],gray[:,end])
            mode = median(border) >= 0.5 ? :dark : :light
        end
        used_threshold = threshold === nothing ? otsu_threshold(gray) : Float64(threshold)
        mask = mode == :dark ? BitMatrix(gray .<= used_threshold) : BitMatrix(gray .> used_threshold)
    end
    remove_small_components!(mask,min_component_pixels)
    any(mask) || error("После порога/удаления шума контур пуст. Проверьте foreground, threshold и min_component_pixels.")
    all(mask) && error("Всё изображение выделено как объект. Проверьте foreground и threshold.")
    return mask, (original_size=original_size, processed_size=size(mask),
                  foreground=mode, threshold=used_threshold)
end

# Границы пиксельных ячеек. Все внешние контуры и отверстия сохраняются.
# В неоднозначной вершине диагонального касания выбирается правый поворот.
function boundary_paths(mask)
    h, w = size(mask)
    edges = Edge[]
    for r in 1:h, c in 1:w
        mask[r,c] || continue
        x, y = c - 1, r - 1
        (r == 1 || !mask[r-1,c]) && push!(edges,((x,y),(x+1,y)))
        (c == w || !mask[r,c+1]) && push!(edges,((x+1,y),(x+1,y+1)))
        (r == h || !mask[r+1,c]) && push!(edges,((x+1,y+1),(x,y+1)))
        (c == 1 || !mask[r,c-1]) && push!(edges,((x,y+1),(x,y)))
    end
    length(edges) <= 6000 || error("Граница слишком сложная: больше 6000 рёбер. Уменьшите размер обработки.")
    outgoing = Dict{Point,Vector{Point}}()
    for (a,b) in edges; push!(get!(outgoing,a,Point[]),b); end
    used = Set{Edge}(); paths = Vector{ComplexF64}[]
    for first_edge in edges
        first_edge in used && continue
        start, nxt = first_edge
        vertices = Point[start]; previous = start; current = nxt
        push!(used,first_edge)
        while current != start
            push!(vertices,current)
            options = [q for q in get(outgoing,current,Point[]) if !((current,q) in used)]
            isempty(options) && error("Незамкнутая граница: повреждён порядок рёбер.")
            dx, dy = current[1]-previous[1], current[2]-previous[2]
            dirs = ((-dy,dx),(dx,dy),(dy,-dx),(-dx,-dy))
            chosen = options[1]
            for d in dirs
                q = (current[1]+d[1],current[2]+d[2])
                if q in options; chosen = q; break; end
            end
            push!(used,(current,chosen)); previous,current = current,chosen
        end
        push!(paths,[ComplexF64(x,-y) for (x,y) in vertices])
    end
    return paths, trues(length(paths))
end

# Zhang–Suen: середина толстого штриха, полезно для линейных рисунков.
function thin_mask(mask)
    h,w = size(mask); a = falses(h+2,w+2); a[2:h+1,2:w+1] .= mask
    changed = true
    while changed
        changed = false
        for phase in 1:2
            deleted = Point[]
            for r in 2:h+1, c in 2:w+1
                a[r,c] || continue
                p = (a[r-1,c],a[r-1,c+1],a[r,c+1],a[r+1,c+1],
                     a[r+1,c],a[r+1,c-1],a[r,c-1],a[r-1,c-1])
                2 <= sum(p) <= 6 || continue
                sum(!p[i] && p[mod1(i+1,8)] for i in 1:8) == 1 || continue
                if phase == 1
                    (p[1] && p[3] && p[5]) && continue
                    (p[3] && p[5] && p[7]) && continue
                else
                    (p[1] && p[3] && p[7]) && continue
                    (p[1] && p[5] && p[7]) && continue
                end
                push!(deleted,(r,c))
            end
            for p in deleted; a[p...] = false; end
            changed |= !isempty(deleted)
        end
    end
    return BitMatrix(a[2:h+1,2:w+1])
end

function centerline_paths(mask)
    a = thin_mask(mask); h,w = size(a)
    pixels = [(r,c) for r in 1:h for c in 1:w if a[r,c]]
    length(pixels) <= 6000 || error("Скелет слишком сложный: больше 6000 точек. Уменьшите размер обработки.")
    index = Dict(p=>i for (i,p) in enumerate(pixels))
    neighbors = [Int[] for p in pixels]
    for (i,(r,c)) in enumerate(pixels), dr in -1:1, dc in -1:1
        (dr == 0 && dc == 0) && continue
        j = get(index,(r+dr,c+dc),0); j == 0 && continue
        # Диагональ не дублирует уже существующий прямоугольный путь.
        if dr != 0 && dc != 0 && (get(index,(r+dr,c),0) != 0 || get(index,(r,c+dc),0) != 0)
            continue
        end
        push!(neighbors[i],j)
    end
    edgekey(i,j) = minmax(i,j)
    used = Set{Tuple{Int,Int}}(); paths = Vector{ComplexF64}[]; closed = Bool[]
    starts = vcat(findall(x->length(x)!=2,neighbors),findall(x->length(x)==2,neighbors))
    for i in starts, j in neighbors[i]
        edgekey(i,j) in used && continue
        vertices = Int[i]; previous,current = i,j; push!(used,edgekey(i,j))
        while true
            push!(vertices,current)
            (current == i || length(neighbors[current]) != 2) && break
            k = neighbors[current][1] == previous ? neighbors[current][2] : neighbors[current][1]
            edgekey(current,k) in used && break
            push!(used,edgekey(current,k)); previous,current = current,k
        end
        isclosed = vertices[end] == vertices[1]
        isclosed && pop!(vertices)
        push!(paths,[ComplexF64(pixels[k][2]-0.5,-pixels[k][1]+0.5) for k in vertices])
        push!(closed,isclosed)
    end
    return paths,closed
end

path_length(p,closed=true) = sum(abs.(diff(p))) + (closed ? abs(p[end]-p[1]) : 0.0)
close_path(p,closed) = closed ? vcat(p,p[1]) : vcat(p,reverse(p[1:end-1]))

# Вставка каждого следующего контура в ближайшую точку уже построенного
# маршрута. Мост проходится туда и обратно; границы фигур сохраняются.
function join_paths(paths,closed)
    length(paths) <= 40 && sum(length, paths) <= 6000 ||
        error("Слишком много участков: допустимо до 40 контуров и 6000 вершин. Уменьшите размер обработки или удалите шум.")
    seed = argmax([path_length(p,c) for (p,c) in zip(paths,closed)])
    route = close_path(paths[seed],closed[seed])
    remaining = [i for i in eachindex(paths) if i != seed]
    bridges = Tuple{ComplexF64,ComplexF64}[]
    while !isempty(remaining)
        best = Inf; route_idx = 1; path_idx = remaining[1]; point_idx = 1
        for j in remaining, b in eachindex(paths[j]), a in 1:length(route)-1
            d = abs2(route[a]-paths[j][b])
            if d < best
                best=d; route_idx=a; path_idx=j; point_idx=b
            end
        end
        p = paths[path_idx]
        if closed[path_idx]
            rotated = vcat(p[point_idx:end],p[1:point_idx-1])
            excursion = close_path(rotated,true)
        else
            # Открытый штрих обходится из точки входа к обоим концам и назад.
            excursion = vcat(p[point_idx:end],reverse(p[1:end-1]),p[2:point_idx])
        end
        anchor = route[route_idx]
        push!(bridges,(anchor,excursion[1]))
        route = vcat(route[1:route_idx],excursion,anchor,route[route_idx+1:end])
        filter!(i->i!=path_idx,remaining)
    end
    return route,bridges
end

function resample_path(path,n)
    n <= 4096 || error("Для учебной анимации допускается до 4096 отсчётов.")
    # Удаляем нулевые рёбра; затем равномерно семплируем по длине пути.
    p = ComplexF64[path[1]]
    for z in path[2:end]; abs(z-p[end]) > 1e-12 && push!(p,z); end
    abs(p[end]-p[1]) > 1e-12 && push!(p,p[1])
    cumulative = vcat(0.0,cumsum(abs.(diff(p))))
    total = cumulative[end]; total > 0 || error("Траектория имеет нулевую длину")
    result = Vector{ComplexF64}(undef,n); edge = 1
    for k in 0:n-1
        s = k * total / n
        while edge < length(p)-1 && cumulative[edge+1] <= s; edge += 1; end
        f = (s-cumulative[edge])/(cumulative[edge+1]-cumulative[edge])
        result[k+1] = p[edge] + f*(p[edge+1]-p[edge])
    end
    return result
end

"""Минимальное K по энергии отброшенных гармоник; предел защищает GIF от перегрузки."""
function select_harmonics(C; accuracy_percent=99.0, n_segments=nothing, max_segments=300)
    order = sortperm(abs2.(C[2:end]); rev=true) .+ 1
    energy = sum(abs2, C[2:end])
    if n_segments === nothing
        budget = energy * (1 - accuracy_percent / 100)^2
        tail = vcat(reverse(cumsum(reverse(abs2.(C[order])))), 0.0)
        needed = something(findfirst(k -> tail[k+1] <= budget, 1:length(order)), length(order))
    else
        needed = n_segments
    end
    k = min(needed, max_segments, length(order))
    return (indices=order[1:k], needed=needed, capped=k < needed)
end

chain(model,t) = vcat(model.offset, model.offset .+
    cumsum(model.coefficients .* cis.(2π .* model.frequencies .* t)))

"""Плотная интерполяция ряда: геометрия следа не зависит от FPS."""
function reconstruct(model; count=8192)
    count = max(count, length(model.samples))
    spectrum = zeros(ComplexF64, count)
    spectrum[1] = model.offset
    for (c,f) in zip(model.coefficients, model.frequencies)
        spectrum[mod(f,count)+1] += c
    end
    return vcat(ifft(spectrum) * count, chain(model,1.0)[end])
end

function separated_xy(paths)
    x=Float64[]; y=Float64[]
    for p in paths
        append!(x,real.(p)); append!(y,imag.(p)); push!(x,NaN); push!(y,NaN)
    end
    return x,y
end

function geometry_plot(mask, paths, closed, bridges)
    gr(format=:png); default(fontfamily="DejaVu Sans")
    p1=heatmap(1:size(mask,2), 1:size(mask,1), Float64.(mask);
        color=[:white,:black], yflip=true, aspect_ratio=:equal,
        colorbar=false, axis=false, title="Бинарная маска", titlefontsize=12, legend=false)
    p2=plot(;aspect_ratio=:equal, axis=false, legend=false,
        title="Контуры: $(length(paths)) | мосты: $(length(bridges))",titlefontsize=12)
    x,y=separated_xy([c ? vcat(p,p[1]) : p for (p,c) in zip(paths,closed)])
    plot!(p2,x,y;color="#0284c7",linewidth=2)
    bx,by=separated_xy([[a,b] for (a,b) in bridges])
    plot!(p2,bx,by;color="#f59e0b",linewidth=2,linestyle=:dash)
    return plot(p1,p2;layout=(1,2),size=(1000,460))
end

function spectrum_plot(C, frequencies)
    gr(format=:png); default(fontfamily="DejaVu Sans")
    order=sortperm(abs.(C[2:end]);rev=true).+1
    p1=scatter(frequencies[2:end],abs.(C[2:end]);markersize=2,
        xlabel="Частота, обороты за период",ylabel="Амплитуда, пиксели",legend=false,
        title="Спектр комплексной траектории",titlefontsize=11,color="#0284c7")
    p2=bar(1:min(30,length(order)),abs.(C[order[1:min(30,length(order))]]);
        xlabel="Порядок по амплитуде",ylabel="Амплитуда, пиксели",legend=false,
        title="30 крупнейших гармоник",titlefontsize=11,color="#fb7185")
    return plot(p1,p2;layout=(1,2),size=(1040,480),margin=6Plots.mm)
end

function comparison_plot(route, reconstructed, model)
    gr(format=:png); default(fontfamily="DejaVu Sans")
    p=plot(real.(route),imag.(route);color="#94a3b8",linewidth=2,
        label="Целевой маршрут с мостами",aspect_ratio=:equal,
        xlabel="x, пиксели",ylabel="y, пиксели",size=(880,560),legend=:outerbottom,
        margin=6Plots.mm,titlefontsize=11,title=@sprintf("K = %d | точность = %.3f%% | RMS = %.3f пикселя",
            length(model.coefficients),model.accuracy,model.rms_error))
    plot!(p,real.(reconstructed),imag.(reconstructed);
        color="#e11d48",linewidth=1.5,label="Восстановление рядом Фурье")
    return p
end

"""Габариты цепочки и окружностей с запасом по производной между проверяемыми фазами."""
function canvas_limits(model, route, points)
    z=vcat(route,points)
    xmin,xmax=extrema(real.(z)); ymin,ymax=extrema(imag.(z))
    count=max(512,min(8192,4maximum(abs,model.frequencies)))
    times=range(0,1;length=count+1)
    positions=fill(model.offset,count+1); derivative_bound=0.0
    for (c,f) in zip(model.coefficients,model.frequencies)
        radius=abs(c); gap=derivative_bound/(2count)
        xmin=min(xmin,minimum(real,positions)-radius-gap)
        xmax=max(xmax,maximum(real,positions)+radius+gap)
        ymin=min(ymin,minimum(imag,positions)-radius-gap)
        ymax=max(ymax,maximum(imag,positions)+radius+gap)
        positions .+= c .* cis.(2π .* f .* times)
        derivative_bound += 2π * abs(f) * radius
    end
    gap=derivative_bound/(2count)
    xmin=min(xmin,minimum(real,positions)-gap); xmax=max(xmax,maximum(real,positions)+gap)
    ymin=min(ymin,minimum(imag,positions)-gap); ymax=max(ymax,maximum(imag,positions)+gap)
    margin=0.06max(xmax-xmin,ymax-ymin,1)
    return (xmin-margin,xmax+margin),(ymin-margin,ymax+margin)
end

function frame_plot(model, route, points, t, style, limits; show_circles=true)
    nodes=chain(model,t)
    label=@sprintf("K = %d  |  точность = %.3f%%  |  нарисовано = %3.0f%%",
        length(model.coefficients),model.accuracy,100t)
    p=plot(;xlims=limits[1],ylims=limits[2],aspect_ratio=:equal,
        legend=false,axis=false,grid=false,background_color=style.background,
        foreground_color=style.segments,title=label,titlefontsize=11,
        size=(900,760),margin=5Plots.mm)
    plot!(p,real.(route),imag.(route);color=style.reference,linewidth=1)
    last_idx=min(length(points)-1,floor(Int,t*(length(points)-1))+1)
    trace=vcat(points[1:last_idx],nodes[end])
    plot!(p,real.(trace),imag.(trace);color=style.contour,linewidth=2.2)
    if show_circles
        angles=range(0,2π;length=49); cx=Float64[]; cy=Float64[]
        for j in eachindex(model.coefficients)
            circle=nodes[j] .+ abs(model.coefficients[j]) .* cis.(angles)
            append!(cx,real.(circle)); append!(cy,imag.(circle))
            push!(cx,NaN); push!(cy,NaN)
        end
        plot!(p,cx,cy;color=style.circles,linealpha=0.35,linewidth=0.7)
    end
    plot!(p,real.(nodes),imag.(nodes);color=style.segments,linewidth=1.3)
    scatter!(p,[real(nodes[end])],[imag(nodes[end])];color=style.contour,
        markersize=3,markerstrokewidth=0)
    return p
end

function save_results(model, mask, meta, output_dir)
    mkpath(output_dir)
    FileIO.save(joinpath(output_dir,"mask.png"), Gray.(1 .- Float64.(mask)))
    open(joinpath(output_dir,"coefficients.csv"),"w") do io
        println(io,"segment,frequency,real,imag,amplitude,phase")
        for (j,(c,f)) in enumerate(zip(model.coefficients,model.frequencies))
            println(io,join((j,f,real(c),imag(c),abs(c),angle(c)),","))
        end
    end
    open(joinpath(output_dir,"report.txt"),"w") do io
        println(io,"Входной размер: ",meta.original_size,"; размер обработки: ",meta.processed_size)
        println(io,"Режим: ",meta.foreground,"; порог: ",meta.threshold)
        println(io,"N: ",length(model.samples),"; K: ",length(model.coefficients))
        @printf(io,"Точность: %.6f%%; RMS: %.6f px; максимальная ошибка: %.6f px\n",
            model.accuracy,model.rms_error,model.max_error)
        println(io,"Метрика относится к отсчётам подготовленного маршрута, включая мосты.")
        println(io,"Сработал предел числа отрезков: ",model.capped)
    end
end

"""Один период рисования и пауза; возвращает родной объект Plots для показа GIF в Engee."""
function render_gif(model,route,output_path;style,fps=25,duration=12.0,
        hold_seconds=1.0,show_circles=true)
    frames=max(2,round(Int,duration*fps)); holds=round(Int,hold_seconds*fps)
    frames+holds <= 900 || error("Больше 900 кадров. Уменьшите FPS или длительность.")
    gr(format=:png); default(fontfamily="DejaVu Sans"); mkpath(dirname(output_path))
    points=reconstruct(model); limits=canvas_limits(model,route,points)
    result=mktempdir() do temp_dir
        animation=Plots.Animation(temp_dir)
        last_plot=nothing
        for i in 0:frames-1
            t=i/(frames-1)
            last_plot=frame_plot(model,route,points,t,style,limits;show_circles)
            Plots.frame(animation,last_plot)
        end
        for i in 1:holds; Plots.frame(animation,last_plot); end
        savefig(last_plot,splitext(output_path)[1]*"_final.png")
        Plots.gif(animation,output_path;fps,loop=0,show_msg=false)
    end
    return (animation=result,path=output_path,frames=frames+holds,
        requested_fps=fps,duration=(frames+holds)/fps)
end

end # module EpicycleLesson