using Pkg

# -----------------------------------------------------------------------------------------------
# Installation des packages : mettre à false une fois les packages installés une première fois.
# -----------------------------------------------------------------------------------------------
const FIRST_RUN = false

if FIRST_RUN
    for pkg in ["Statistics", "FITSIO", "NativeFileDialog"]
        Pkg.add(pkg)
    end
end

using Statistics
using Printf

# fits file  : https://juliaastro.org/FITSIO.jl/stable/
using FITSIO

# package pour sélection interactive du répertoire et du fichier (fenêtres de dialogue)
using NativeFileDialog

##################################################################################################################
# Fonctions du programme :
#   Portage Julia de l'algorithme FNRGF (Fourier Normalizing Radial Graded Filter), tel qu'implémenté
#   dans le package Python `sunkit-image` (PyPI, version 0.7.0, module sunkit_image.radial, fonction
#   `fnrgf()`, et fonctions utilitaires de sunkit_image.utils.utils : `find_radial_bin_edges`,
#   `_set_attenuation_coefficients`). Comme pour le portage de NRGF (voir TSE2026-4-eclipse-NRGF-v5.jl),
#   ce script est basé sur la lecture directe du code source réel de sunkit-image (package téléchargé
#   depuis PyPI), pas seulement de sa documentation.
#   https://docs.sunpy.org/projects/sunkit-image/en/stable/api/sunkit_image.radial.fnrgf.html
#
#   Référence de l'algorithme :
#   Morgan, H., Habbal, S. R., & Druckmüllerová, H. 2011, "Enhancing Coronal Structures with the
#   Fourier Normalizing-Radial-Graded Filter", Astrophysical Journal, 737, 88.
#   https://iopscience.iop.org/article/10.1088/0004-637X/737/2/88
#   L'implémentation de sunkit-image s'inspire aussi directement de la thèse de doctorat de
#   Druckmüllerová : "Application of adaptive filters in processing of solar corona images",
#   https://dspace.vutbr.cz/bitstream/handle/11012/34520/DoctoralThesis.pdf (chapitre 6).
#   Je n'ai trouvé aucun portage indépendant de FNRGF sur GitHub (Python, MATLAB ou IDL) en dehors de
#   sunkit-image - seulement l'article et la thèse ci-dessus.
#
# VERSION -v1 : premier jet, traitement RVB uniquement - PAS d'option de choix d'un "canal diviseur"
#   commun (L/R/G/B) comme dans NRGF v5. Chaque canal R, G, B est normalisé par son PROPRE profil de
#   Fourier (mean/std azimutal par anneau), exactement comme le fait NRGF v5 en mode :rgb.
#
# NOTE SUR LA FIDÉLITÉ DU PORTAGE (différences volontaires avec sunkit-image, mêmes raisons que pour
# NRGF v5) :
#   1. sunkit-image calcule l'angle et le rayon de chaque pixel à partir des coordonnées WCS d'un
#      `sunpy.map.Map` (Helioprojective, via `smap.pixel_to_world`). Une photo d'éclipse n'a pas ce
#      système de coordonnées : ce script calcule directement le rayon et l'angle de chaque pixel par
#      rapport au centre du Soleil/de la Lune EN PIXELS (voir `cx_sun_px`, `cy_sun_px`, `r_sun_px`
#      dans main(), à ajuster à la main pour chaque image - mêmes paramètres que dans NRGF v5).
#   2. Le code Python boucle sur chaque anneau et reconstruit un masque booléen sur l'image ENTIÈRE à
#      chaque itération pour retrouver les pixels de l'anneau - beaucoup trop lent sur une grande
#      image. Ce script pré-calcule une seule fois, pour chaque anneau, la LISTE des indices de pixels
#      qui lui appartiennent (tri par "seau", en une seule passe sur l'image - voir
#      build_ring_pixel_indices), puis réutilise directement cette liste dans la boucle sur les
#      anneaux : même résultat mathématique, juste vectorisé/indexé pour la performance.
#
# QUIRKS DE L'IMPLÉMENTATION D'ORIGINE, REPRODUITS FIDÈLEMENT ICI (volontairement, pour rester fidèle
# au comportement réel de sunkit-image 0.7.0, même là où ce comportement semble surprenant) :
#   a. Les coefficients d'atténuation du terme d'ordre 0 (a_0 pour la moyenne, c_0 pour l'écart-type)
#      utilisent, dans le code source, le MÊME coefficient d'atténuation que l'harmonique k=1
#      (`attenuation_coefficients[0, 1]`), et non un coefficient dédié à l'ordre 0
#      (`attenuation_coefficients[0, 0]`), alors qu'un tableau de coefficients dédiés existe bel et
#      bien pour cet indice 0. C'est très probablement une coquille dans sunkit-image, mais comme
#      pour les quirks déjà reproduits dans NRGF v5 (comportement de np.std vis-à-vis des NaN), on
#      reproduit ici le comportement RÉEL du code publié plutôt que ce qu'il "devrait" faire.
#   b. Un segment angulaire totalement vide (aucun pixel) reçoit une moyenne ET un écart-type de 0
#      (test explicite dans le code source, pour éviter un avertissement "mean of empty slice").
#      Un segment NON vide mais dont TOUS les pixels sont NaN donnera en revanche une moyenne NaN
#      (comportement de np.nanmean sur un tableau entièrement NaN) - cette distinction entre "vide"
#      et "non vide mais tout NaN" est intentionnellement reproduite telle quelle.
#   c. `ratio_mix` (par défaut `[15, 1]`) n'est PAS une paire de poids normalisés (qui sommeraient à
#      1) : le code source calcule littéralement `ratio_mix[0]*original + ratio_mix[1]*filtré`, SANS
#      diviser par leur somme. Avec les valeurs par défaut, le résultat est donc dominé par 15 fois
#      l'image originale (dans ses unités d'origine, ex. DN/pixel) plus 1 fois le champ normalisé
#      (d'ordre de grandeur ~1). Ce n'est pas un défaut du portage : c'est le comportement exact du
#      code source, pensé pour être suivi d'un étirement/une normalisation visuelle a posteriori
#      (affichage avec normalisation automatique, traitement dans Siril/PixInsight, etc.), pas pour
#      produire directement une image "prête à l'affichage".
#
# Principe (algorithme d'origine) :
#   1. Pour chaque pixel, calcul de sa distance au centre du Soleil/de la Lune (en rayons solaires) et
#      de son angle azimutal (dans [0, 2π)).
#   2. Découpage de l'image en anneaux radiaux (bins) régulièrement espacés entre 0 et le rayon
#      maximal atteint dans l'image, et en segments angulaires (par défaut 130 segments égaux sur le
#      cercle complet).
#   3. Pour chaque anneau ET chaque segment angulaire, calcul de deux statistiques sur les pixels
#      qu'il contient (au-delà du rayon d'application) : la moyenne (ignore les NaN, comme
#      np.nanmean) et l'écart-type (n'ignore PAS les NaN, comme np.std - voir quirk (b) ci-dessus).
#   4. Pour chaque anneau, ajustement d'une série de Fourier tronquée à l'ordre `order` (par défaut 3)
#      sur les moyennes/écarts-types des segments de cet anneau (fonction de l'angle du segment), avec
#      atténuation linéaire optionnelle des coefficients de haut ordre (voir quirk (a) ci-dessus).
#   5. Pour chaque pixel de l'anneau (au-delà du rayon d'application), la moyenne et l'écart-type
#      "locaux" sont approximés en évaluant cette série de Fourier À L'ANGLE PROPRE DU PIXEL (pas
#      celui du segment) :
#        champ_normalisé = (valeur - moyenne_approximée(angle)) / écart_type_approximé(angle)
#      (si l'écart-type approximé est nul, il est remplacé par 1, comme dans le code source).
#   6. Mélange linéaire (NON normalisé - voir quirk (c) ci-dessus) avec l'image d'origine :
#        nouvelle_valeur = ratio_mix[1]·valeur + ratio_mix[2]·champ_normalisé
#      Les pixels en-deçà du rayon d'application restent à la valeur de remplissage (`fill_value`).
###################################################################################################################

# Calcule les coefficients d'atténuation utilisés pour pondérer les coefficients de Fourier de haut
# ordre (portage direct de _set_attenuation_coefficients de sunkit-image). Renvoie une matrice
# (2, order+1) : ligne 1 = coefficients pour la moyenne, ligne 2 = coefficients pour l'écart-type.
# Chaque ligne varie linéairement entre son "range" [haut, bas] sur les order+1 coefficients
# (indices 0 à order en convention Python, 1 à order+1 en convention Julia 1-based).
function set_attenuation_coefficients(order::Int;
    mean_attenuation_range::Vector{Float64}=[1.0, 0.0],
    std_attenuation_range::Vector{Float64}=[1.0, 0.0],
    cutoff::Int=0)

    cutoff > (order + 1) && error("cutoff ne peut pas dépasser order + 1")

    coeffs = zeros(Float64, 2, order + 1)
    coeffs[1, :] = collect(range(mean_attenuation_range[1], mean_attenuation_range[2], length=order + 1))
    coeffs[2, :] = collect(range(std_attenuation_range[1], std_attenuation_range[2], length=order + 1))

    if cutoff != 0
        coeffs[:, (end-cutoff+1):end] .= 0.0
    end

    return coeffs
end

# Construit les matrices cos_matrix/sin_matrix (nseg x order) utilisées pour projeter les profils
# moyenne/écart-type par segment sur la base de Fourier - portage direct des cos_matrix/sin_matrix du
# code source. cos_matrix[s,k] = cos(2π·k·(s-0.5)/nseg), même convention que le "(i+0.5)" du code
# Python (angle au CENTRE de chaque segment). Ne dépend que de nseg et order : calculé une seule fois,
# indépendamment du canal ou de l'anneau.
function build_fourier_basis(nseg::Int, order::Int)
    cos_matrix = Matrix{Float64}(undef, nseg, order)
    sin_matrix = Matrix{Float64}(undef, nseg, order)
    @inbounds for k in 1:order, s in 1:nseg
        arg = 2π * k * (s - 0.5) / nseg
        cos_matrix[s, k] = cos(arg)
        sin_matrix[s, k] = sin(arg)
    end
    return cos_matrix, sin_matrix
end

# *** Optimisation (voir note en en-tête) *** : pré-calcule, pour chaque anneau radial, la liste des
# indices linéaires des pixels qui lui appartiennent ET qui sont au-delà du rayon d'application
# (valid_mask). Tri par seau ("counting sort") en deux passes sur l'image entière : une pour compter
# la taille de chaque anneau, une pour remplir les listes - complexité linéaire en nombre de pixels,
# sans jamais reconstruire de masque booléen sur l'image entière à répétition.
function build_ring_pixel_indices(bin_idx::Matrix{Int}, valid_mask::BitMatrix, nbins::Int)
    counts = zeros(Int, nbins)
    @inbounds for idx in eachindex(bin_idx)
        if valid_mask[idx]
            counts[bin_idx[idx]] += 1
        end
    end

    ring_pixel_indices = Vector{Vector{Int}}(undef, nbins)
    for b in 1:nbins
        ring_pixel_indices[b] = Vector{Int}(undef, counts[b])
    end

    fill_pos = zeros(Int, nbins)
    @inbounds for idx in eachindex(bin_idx)
        if valid_mask[idx]
            b = bin_idx[idx]
            fill_pos[b] += 1
            ring_pixel_indices[b][fill_pos[b]] = idx
        end
    end

    return ring_pixel_indices
end

# Calcule, pour chaque paire (anneau, segment angulaire), la moyenne (ignore les NaN, comme
# np.nanmean) et l'écart-type "population" (ddof=0, n'ignore PAS les NaN, comme np.std) des valeurs
# des pixels qu'elle contient - portage direct de la boucle interne de fnrgf() (calcul de
# average_segments/std_dev), mais vectorisé en une seule passe sur l'image plutôt qu'en une boucle
# anneau x segment reconstruisant un masque à chaque itération (voir quirk (b) en en-tête pour le
# traitement des segments vides vs. segments non vides mais entièrement NaN).
function compute_ring_segment_stats(data::AbstractMatrix{T}, bin_idx::Matrix{Int}, seg_idx::Matrix{Int},
    valid_mask::BitMatrix, nbins::Int, nseg::Int) where {T<:AbstractFloat}

    sum_valid = zeros(Float64, nbins, nseg)
    count_valid = zeros(Int, nbins, nseg)
    sum_all = zeros(Float64, nbins, nseg)
    sumsq_all = zeros(Float64, nbins, nseg)
    count_all = zeros(Int, nbins, nseg)

    @inbounds for idx in eachindex(data)
        if valid_mask[idx]
            b = bin_idx[idx]
            s = seg_idx[idx]
            v = Float64(data[idx])
            count_all[b, s] += 1
            sum_all[b, s] += v
            sumsq_all[b, s] += v * v
            if !isnan(v)
                count_valid[b, s] += 1
                sum_valid[b, s] += v
            end
        end
    end

    average_segments = zeros(Float64, nbins, nseg)
    std_dev = zeros(Float64, nbins, nseg)
    @inbounds for s in 1:nseg, b in 1:nbins
        if count_all[b, s] == 0
            # segment vide : 0 et 0, comme le test explicite du code source (pas de NaN ici)
            average_segments[b, s] = 0.0
            std_dev[b, s] = 0.0
        else
            average_segments[b, s] = count_valid[b, s] == 0 ? NaN : sum_valid[b, s] / count_valid[b, s]
            m_all = sum_all[b, s] / count_all[b, s]   # NaN si un pixel du segment est NaN (non filtré)
            std_dev[b, s] = sqrt(max(sumsq_all[b, s] / count_all[b, s] - m_all^2, 0.0))
        end
    end

    return average_segments, std_dev
end

# *** Cœur de l'algorithme *** : applique le FNRGF à un canal, à partir des statistiques par
# (anneau, segment) déjà calculées (average_segments, std_dev) et des listes de pixels par anneau déjà
# pré-calculées (ring_pixel_indices). Pour chaque anneau : projection des profils sur la base de
# Fourier tronquée à l'ordre `order` (avec atténuation - voir quirk (a) en en-tête), puis évaluation
# de cette approximation à l'angle propre de CHAQUE pixel de l'anneau (pas celui d'un segment) pour
# obtenir moyenne_approximée(angle)/écart_type_approximé(angle), et enfin normalisation + mélange
# (non normalisé - voir quirk (c) en en-tête) avec l'image d'origine.
# Renvoie deux tableaux : le champ final mélangé (image_mixed) et le champ normalisé seul, avant
# mélange (image_residual, utile pour juger l'effet du filtre isolément).
function fnrgf_apply_channel(data::AbstractMatrix{T}, angles::Matrix{Float64},
    average_segments::Matrix{Float64}, std_dev::Matrix{Float64}, ring_pixel_indices::Vector{Vector{Int}},
    cos_matrix::Matrix{Float64}, sin_matrix::Matrix{Float64}, attenuation_coefficients::Matrix{Float64},
    order::Int, nseg::Int, ratio_mix::Vector{Float64}; fill_value::Float64=0.0) where {T<:AbstractFloat}

    dims = size(data)
    image_mixed = fill(fill_value, dims)
    image_residual = fill(fill_value, dims)

    a_k = Vector{Float64}(undef, order)
    b_k = Vector{Float64}(undef, order)
    c_k = Vector{Float64}(undef, order)
    d_k = Vector{Float64}(undef, order)

    @inbounds for i in eachindex(ring_pixel_indices)
        idxs = ring_pixel_indices[i]
        isempty(idxs) && continue

        avg_row = @view average_segments[i, :]
        std_row = @view std_dev[i, :]

        # coefficients d'ordre 0 - voir quirk (a) : utilisent attenuation_coefficients[.,2] (= indice
        # python 1, celui de l'harmonique k=1), pas attenuation_coefficients[.,1] (indice python 0)
        a0 = sum(avg_row) * (2 / nseg) * attenuation_coefficients[1, 2]
        c0 = sum(std_row) * (2 / nseg) * attenuation_coefficients[2, 2]

        for k in 1:order
            ck_col = @view cos_matrix[:, k]
            sk_col = @view sin_matrix[:, k]
            atten_mean_k = attenuation_coefficients[1, k+1]
            atten_std_k = attenuation_coefficients[2, k+1]
            a_k[k] = (avg_row' * ck_col) * (2 / nseg) * atten_mean_k
            b_k[k] = (avg_row' * sk_col) * (2 / nseg) * atten_mean_k
            c_k[k] = (std_row' * ck_col) * (2 / nseg) * atten_std_k
            d_k[k] = (std_row' * sk_col) * (2 / nseg) * atten_std_k
        end

        θ = angles[idxs]
        mean_approx = fill(a0 / 2, length(idxs))
        std_approx = fill(c0 / 2, length(idxs))
        for k in 1:order
            kθ = k .* θ
            ck = cos.(kθ)
            sk = sin.(kθ)
            mean_approx .+= a_k[k] .* ck .+ b_k[k] .* sk
            std_approx .+= c_k[k] .* ck .+ d_k[k] .* sk
        end
        std_approx[std_approx.==0] .= 1.0   # comme dans le code source : évite la division par 0

        vals = Float64.(data[idxs])
        resid = (vals .- mean_approx) ./ std_approx
        image_residual[idxs] .= resid
        image_mixed[idxs] .= ratio_mix[1] .* vals .+ ratio_mix[2] .* resid
    end

    return image_mixed, image_residual
end

###############################################################################################
#
## début du programme principal
#
###############################################################################################

function main()

    # chronométrage du temps total de calcul
    t_start = time()

    println("Nombre de threads Julia disponibles : ", Threads.nthreads())
    if Threads.nthreads() == 1
        println("ATTENTION : un seul thread disponible - le calcul des 3 canaux R,G,B restera")
        println("            séquentiel. Relancez avec 'julia -t 3 ...' (ou -t auto) pour paralléliser.")
    end

    # sélection interactive du répertoire de travail
    work_dir = pick_folder(pwd())
    isempty(work_dir) && error("Aucun répertoire sélectionné - programme arrêté")
    cd(work_dir)

    println(pwd())

    # sélection interactive du fichier fit à traiter
    selected_file = pick_file(work_dir; filterlist="fit,fits")
    isempty(selected_file) && error("Aucun fichier sélectionné - programme arrêté")
    file_name_in = basename(selected_file)

    # lecture de l'image d'entrée (fit, 3 couches RGB)
    f_in = FITS(file_name_in)
    header = read_header(f_in[1])
    image_in = read(f_in[1])
    close(f_in)

    dim_x = header["NAXIS1"]
    dim_y = header["NAXIS2"]
    println("dim_x = ", dim_x)
    println("dim_y = ", dim_y)

    # conversion en Float32 (au cas où l'image d'entrée serait dans un autre type)
    image_in = Float32.(image_in)

    # -----------------------------------------------------------------------------------------------
    # géométrie du disque solaire/lunaire dans l'image : ces 3 valeurs DOIVENT être ajustées à la
    # main pour CHAQUE image (mêmes paramètres, mêmes conventions que dans NRGF v5). Convention
    # 1-based : x = colonne (1 à dim_x), y = ligne (1 à dim_y).
    # -----------------------------------------------------------------------------------------------
    cx_sun_px = 4050   # <-- À AJUSTER : coordonnée x du centre du Soleil/de la Lune, en pixels
    cy_sun_px = 2700   # <-- À AJUSTER : coordonnée y du centre du Soleil/de la Lune, en pixels
    r_sun_px = 600.0       # <-- À AJUSTER : rayon apparent du Soleil (ou de la Lune), en pixels

    # -----------------------------------------------------------------------------------------------
    # paramètres FNRGF (noms et valeurs par défaut alignés sur ceux de sunkit-image 0.7.0)
    # -----------------------------------------------------------------------------------------------
    nbins_fnrgf = dim_x ÷ 2                 # nombre d'anneaux radiaux (= shape[0]//2, comme le défaut de sunkit-image)
    number_angular_segments_fnrgf = 130     # nombre de segments angulaires par anneau (défaut sunkit-image)
    order_fnrgf = 4                         # ordre de la série de Fourier (défaut sunkit-image ; minimum 1)
    application_radius_fnrgf = 1.0          # en rayons solaires ; FNRGF appliqué seulement au-delà
    mean_attenuation_range_fnrgf = [1.0, 0.5]  # <-- À AJUSTER si besoin (défaut sunkit-image)
    std_attenuation_range_fnrgf = [1.0, 0.5]   # <-- À AJUSTER si besoin (défaut sunkit-image)
    cutoff_fnrgf = 0                        # <-- À AJUSTER si besoin (défaut sunkit-image : pas de cutoff)
    ratio_mix_fnrgf = [15.0, 1.0]            # <-- À AJUSTER : [K1, K2] NON normalisés (voir quirk (c) en en-tête)
    fill_value_fnrgf = 0.0                  # valeur hors de la zone d'application (0.0 plutôt que NaN,
    # comme dans NRGF v5, pour faciliter les traitements en aval)

    @assert order_fnrgf >= 1 "order_fnrgf doit être >= 1"

    # -----------------------------------------------------------------------------------------------
    # calcul du rayon (en rayons solaires) et de l'angle azimutal (dans [0, 2π)) de chaque pixel,
    # par rapport au centre du Soleil/de la Lune - mêmes conventions géométriques que NRGF v5.
    # -----------------------------------------------------------------------------------------------
    println("Calcul des distances et des angles au centre solaire pour chaque pixel...")
    flush(stdout)
    map_r = Matrix{Float64}(undef, dim_x, dim_y)
    angles = Matrix{Float64}(undef, dim_x, dim_y)
    @inbounds for y in 1:dim_y, x in 1:dim_x
        dx = Float64(x) - cx_sun_px
        dy = Float64(y) - cy_sun_px
        map_r[x, y] = sqrt(dx^2 + dy^2) / r_sun_px
        θ = atan(dy, dx)
        angles[x, y] = θ < 0 ? θ + 2π : θ
    end
    println(@sprintf("Rayon maximal atteint dans l'image : %.2f rayons solaires (vérifier que ceci est cohérent avec cx_sun_px/cy_sun_px/r_sun_px)", maximum(map_r)))
    flush(stdout)

    rmax = maximum(map_r)
    bin_width = rmax / nbins_fnrgf
    bin_idx = clamp.(floor.(Int, map_r ./ bin_width) .+ 1, 1, nbins_fnrgf)

    segment_angle = 2π / number_angular_segments_fnrgf
    seg_idx = clamp.(floor.(Int, angles ./ segment_angle) .+ 1, 1, number_angular_segments_fnrgf)

    valid_mask = map_r .> application_radius_fnrgf

    # -----------------------------------------------------------------------------------------------
    # pré-calculs communs aux 3 canaux (indépendants des données elles-mêmes) : liste des pixels par
    # anneau, base de Fourier, coefficients d'atténuation
    # -----------------------------------------------------------------------------------------------
    println("Pré-calcul des listes de pixels par anneau...")
    flush(stdout)
    ring_pixel_indices = build_ring_pixel_indices(bin_idx, valid_mask, nbins_fnrgf)

    cos_matrix, sin_matrix = build_fourier_basis(number_angular_segments_fnrgf, order_fnrgf)
    attenuation_coefficients = set_attenuation_coefficients(order_fnrgf;
        mean_attenuation_range=mean_attenuation_range_fnrgf,
        std_attenuation_range=std_attenuation_range_fnrgf,
        cutoff=cutoff_fnrgf)

    image_fnrgf = similar(image_in)
    image_fnrgf_residual = similar(image_in)

    println("Démarrage du calcul FNRGF...")
    flush(stdout)

    # -----------------------------------------------------------------------------------------------
    # traitement indépendant de chaque canal couleur (R,G,B), en parallèle sur des threads Julia
    # distincts (voir 'julia -t auto ...' pour activer plusieurs threads) - traitement RVB uniquement
    # (voir en-tête) : chaque canal est normalisé par son PROPRE profil de Fourier.
    # -----------------------------------------------------------------------------------------------
    Threads.@threads for i_ch in [1, 2, 3]
        println("calcul FNRGF, canal ", i_ch, " (thread ", Threads.threadid(), ")")
        flush(stdout)

        channel_data = image_in[:, :, i_ch]
        average_segments, std_dev = compute_ring_segment_stats(channel_data, bin_idx, seg_idx, valid_mask,
            nbins_fnrgf, number_angular_segments_fnrgf)

        chan_mixed, chan_residual = Base.invokelatest(fnrgf_apply_channel, channel_data, angles,
            average_segments, std_dev, ring_pixel_indices, cos_matrix, sin_matrix,
            attenuation_coefficients, order_fnrgf, number_angular_segments_fnrgf, ratio_mix_fnrgf;
            fill_value=fill_value_fnrgf)

        image_fnrgf[:, :, i_ch] = Float32.(chan_mixed)
        image_fnrgf_residual[:, :, i_ch] = Float32.(chan_residual)

        println("-> canal ", i_ch, " terminé")
        flush(stdout)
    end

    # nom des fichiers de sortie : ajout de "_FNRGF"/"_FNRGF-residual" juste avant l'extension ".fit",
    # du nombre d'anneaux (nb), du nombre de segments angulaires (ns), de l'ordre de Fourier (ord), du
    # rayon d'application x10 sur 2 chiffres (ar), des coefficients ratio_mix (rm), et de l'intervalle
    # mean_attenuation_range (mar, x10 sur 2 chiffres chacun - ex. [1.0,0.0] -> "mar_10_00") pour
    # distinguer facilement les essais faits en jouant sur ce paramètre
    params_str = @sprintf("nb%d_ns%d_ord%d_ar%02d_rm%dx%d_mar_%02d_%02d_", nbins_fnrgf, number_angular_segments_fnrgf,
        order_fnrgf, round(Int, application_radius_fnrgf * 10),
        round(Int, ratio_mix_fnrgf[1]), round(Int, ratio_mix_fnrgf[2]),
        round(Int, mean_attenuation_range_fnrgf[1] * 10), round(Int, mean_attenuation_range_fnrgf[2] * 10))

    file_name_out = chop(file_name_in, tail=4) * "_" * params_str * "FNRGF_rgb.fit"
    f_out = FITS(file_name_out, "w")
    write(f_out, image_fnrgf[:, :, :]; header=header)
    close(f_out)
    println("fichier écrit : ", file_name_out)

    # fichier FNRGF-residual : champ normalisé seul (avant mélange avec l'image d'origine), utile
    # pour juger l'effet du filtre isolément indépendamment du choix de ratio_mix
    file_name_out_residual = chop(file_name_in, tail=4) * "_" * params_str * "FNRGF-residual_rgb.fit"
    f_out = FITS(file_name_out_residual, "w")
    write(f_out, image_fnrgf_residual[:, :, :]; header=header)
    close(f_out)
    println("fichier écrit : ", file_name_out_residual)

    # affichage du temps de calcul total
    elapsed = time() - t_start
    println(@sprintf("Temps de calcul total : %.1f s (%.1f min)", elapsed, elapsed / 60))

end # fin de la fonction main()

# invokelatest en défense contre les avertissements "world age" de Julia 1.12 (peuvent apparaître
# selon la façon dont le script est exécuté, ex. VS Code plutôt qu'un terminal classique)
Base.invokelatest(main)