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 NRGF (Normalizing Radial Graded Filter), tel qu'implémenté dans
#   le package Python `sunkit-image` (PyPI, version 0.7.0, module sunkit_image.radial, fonction
#   `nrgf()`, et fonctions utilitaires de sunkit_image.utils.utils : `find_pixel_radii`,
#   `get_radial_intensity_summary`, `equally_spaced_bins`, `bin_edge_summary`).
#   https://docs.sunpy.org/projects/sunkit-image/en/stable/api/sunkit_image.radial.nrgf.html
#
#   Référence de l'algorithme (et non de Druckmüller, contrairement à FNRGF/RHEF dont certaines
#   variantes s'en inspirent - NRGF lui-même vient de Morgan et al.) :
#   Morgan, H., Habbal, S. R., & Woo, R. 2006, "The Depiction of Coronal Structure in White-Light
#   Images", Solar Physics, 236, 263. https://doi.org/10.1007/s11207-006-0113-6
#
# NOTE SUR LA FIDÉLITÉ DU PORTAGE :
#   Ce script est basé sur la lecture directe du code source réel de `sunkit-image` (package
#   téléchargé depuis PyPI). La logique (découpage en anneaux radiaux, calcul de la moyenne et de
#   l'écart-type par anneau, normalisation) est reproduite fidèlement.
#   DEUX DIFFÉRENCES IMPORTANTES avec le code Python original, volontaires :
#     1. sunkit-image s'appuie sur les coordonnées WCS d'un `sunpy.map.Map` (issues de données
#        spatiales type SDO/AIA) pour calculer la distance de chaque pixel au centre solaire, en
#        utilisant le rayon solaire apparent (`smap.rsun_obs`) comme échelle. Une photo d'éclipse
#        n'a pas ce système de coordonnées : ce script demande donc directement les coordonnées du
#        centre du Soleil/de la Lune et son rayon, EN PIXELS (voir `cx_sun_px`, `cy_sun_px`,
#        `r_sun_px` dans main() ci-dessous) - à ajuster à la main pour chaque image (repérage
#        possible avec un afficheur FITS comme SAOImage DS9, Siril, etc.).
#     2. Le code Python boucle sur chaque anneau et recalcule un masque sur l'image ENTIÈRE à
#        chaque itération (potentiellement des milliers d'anneaux) - beaucoup trop lent sur une
#        grande image (8000x5000 ou plus). Ce script utilise à la place un calcul en un seul
#        passage sur l'image (chaque pixel est affecté une fois pour toutes à son anneau, puis les
#        statistiques de chaque anneau sont accumulées en une seule boucle) : c'est EXACTEMENT le
#        même résultat mathématique, juste vectorisé pour la performance.
#
# ÉVOLUTIONS PAR RAPPORT AU NRGF "PUR" (ajoutées pour corriger des artefacts visuels propres aux
# photos d'éclipse, absents des données SDO/AIA d'origine, PAS pour changer la mesure elle-même -
# le but reste ici la visualisation de l'image, pas une mesure) :
#   - lissage de la queue du profil radial_intensity/radial_width par spline de Hermite cubique,
#     pour supprimer les cassures de pente dues à la tangence anneau/bord d'image (voir
#     smooth_tail_hermite!) ;
#   - filtre gaussien appliqué à l'image radial_luminosity (fond radial reconstruit en 2D) avant
#     utilisation, pour supprimer les effets de discrétisation dus à la largeur finie des anneaux
#     (bin_width) - sinon visibles comme de très légers cercles concentriques ;
#   - choix du "canal diviseur" (rgb/l/r/g/b) : par défaut chaque canal R,G,B est normalisé par son
#     propre profil radial (comme le NRGF d'origine), mais on peut aussi normaliser les 3 canaux
#     par le profil d'un seul (luminance combinée, ou un canal donné), ce qui préserve mieux les
#     rapports de couleur qu'une normalisation strictement indépendante par canal.
#
# Principe (algorithme d'origine) :
#   1. Pour chaque pixel, calcul de sa distance au centre du Soleil/de la Lune, exprimée en rayons
#      solaires (map_r = distance en pixels / r_sun_px).
#   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.
#   3. Pour chaque anneau, calcul de deux statistiques sur les pixels qu'il contient :
#        - la moyenne des valeurs (ignore les éventuels NaN, comme np.nanmean)
#        - l'écart-type des valeurs (n'ignore PAS les NaN, comme np.std - un seul pixel NaN dans
#          l'anneau rend l'écart-type de tout l'anneau égal à NaN, fidèle au comportement original)
#   4. Pour chaque pixel situé au-delà du rayon d'application (`application_radius`, 1 rayon
#      solaire par défaut - le NRGF n'a de sens qu'au-delà du limbe) :
#        nouvelle_valeur = (valeur - moyenne_de_l'anneau) / écart_type_de_l'anneau
#      (si l'écart-type de l'anneau est nul, la division est simplement sautée, comme dans le code
#      source). Les pixels en-deçà du rayon d'application restent à la valeur de remplissage
#      (`fill_value`, NaN par défaut, comme dans sunkit-image).
###################################################################################################################

# division régulière de l'intervalle [inner, outer] en nbins anneaux ; renvoie une matrice (2, nbins)
# avec les bords inférieur (ligne 1) et supérieur (ligne 2) de chaque anneau
function equally_spaced_bins(inner::Float64, outer::Float64, nbins::Int)
    edges = zeros(Float64, 2, nbins)
    for i in 0:(nbins-1)
        edges[1, i+1] = inner + i * (outer - inner) / nbins
        edges[2, i+1] = inner + (i + 1) * (outer - inner) / nbins
    end
    return edges
end

# Lisse la partie du profil au-delà de r_start par une spline de Hermite cubique reliant :
#   - (r_start, valeur du profil à r_start, pente estimée juste avant r_start - zone fiable liée à
#     la tangence Ymax)
#   - (r_end, valeur du profil à r_end, pente estimée autour de r_end - zone fiable liée à la
#     tangence Xmax)
# Au-delà de r_end (jusqu'à r_max), le profil est prolongé en ligne DROITE avec la pente calée en
# r_end (donc valeur ET pente restent continues à r_end - pas de nouveau coude), plutôt que de
# recaler sur la valeur/pente brute du dernier anneau : les coins de l'image (r proche de r_max)
# peuvent souffrir de défauts d'uniformité (flat imparfait) en plus d'être des zones à très faible
# nombre de pixels par anneau (donc statistiquement peu fiables). r_end (Xmax/2 - marge), pris avec
# la même marge de sécurité que r_start, est à la fois dans une zone fiable ET dans une zone où la
# couronne est plus intense (meilleur rapport signal/défaut).
#
# Contexte : un anneau centré sur le Soleil n'est un cercle COMPLET (couverture azimutale uniforme
# sur 360°) que tant que son rayon reste inférieur à la distance du centre au bord le plus proche
# de l'image. Au 1er anneau qui touche le bord haut/bas (rayon r_Y), puis au 1er anneau qui touche
# le bord gauche/droite (rayon r_X), l'anneau devient brutalement tronqué, ce qui casse la
# monotonie et la continuité de pente du profil radial_intensity - d'où l'artefact visible en
# cercle net dans l'image reconstruite. Le but ici est purement la visualisation de l'image (pas
# une mesure) : on accepte de modifier des valeurs de radial_intensity qui n'étaient pourtant pas
# "fausses" avant les seuils de tangence, au profit d'un profil visuellement sans cassure de pente.
function smooth_tail_hermite!(profile::Vector{Float64}, r_centers::Vector{Float64},
    r_start::Float64, r_end::Float64; n_slope::Int=50)

    nbins = length(profile)
    i_start = findfirst(r -> r >= r_start, r_centers)
    isnothing(i_start) && return   # r_start au-delà du rayon max atteint dans l'image : rien à faire
    i_end = findfirst(r -> r >= r_end, r_centers)
    isnothing(i_end) && (i_end = nbins)   # r_end au-delà du rayon max atteint : on va jusqu'au bout
    i_end <= i_start && return

    r_lo, r_hi = r_centers[i_start], r_centers[i_end]
    y_lo, y_hi = profile[i_start], profile[i_end]

    # pente estimée par régression linéaire sur n_slope points (plus robuste qu'une différence
    # finie sur 2 points seuls, sensible au bruit)
    function local_slope(idx_range)
        rr = r_centers[idx_range]
        yy = profile[idx_range]
        X = hcat(ones(length(rr)), rr)
        coeffs = X \ yy
        return coeffs[2]
    end
    m_lo = local_slope(max(1, i_start - n_slope + 1):i_start)                     # pente avant r_start
    half = n_slope ÷ 2
    m_hi = local_slope(max(1, i_end - half):min(nbins, i_end + half))             # pente autour de r_end

    h = r_hi - r_lo
    @inbounds for i in (i_start+1):(i_end-1)
        t = (r_centers[i] - r_lo) / h
        # base de Hermite cubique standard (t in [0,1])
        h00 = 2t^3 - 3t^2 + 1
        h10 = t^3 - 2t^2 + t
        h01 = -2t^3 + 3t^2
        h11 = t^3 - t^2
        profile[i] = h00 * y_lo + h10 * h * m_lo + h01 * y_hi + h11 * h * m_hi
    end

    # au-delà de r_end : prolongation en ligne droite avec la pente m_hi calée en r_end (valeur ET
    # pente continues à r_end, donc pas de nouveau coude à cette jonction)
    @inbounds for i in (i_end+1):nbins
        profile[i] = y_hi + m_hi * (r_centers[i] - r_hi)
    end
end

# Filtre gaussien séparable (flou 2D), avec gestion de bord par réplication (clamp des indices) et
# remplacement préalable des NaN éventuels par 0.0 (zones hors traitement - fill_value) pour éviter
# leur propagation à travers le noyau de convolution. Implémenté directement (sans dépendance
# supplémentaire type ImageFiltering.jl) : le noyau reste petit (rayon ~3*sigma_px), donc rapide.
function gaussian_blur(data::AbstractMatrix{Float64}, sigma_px::Float64)
    sigma_px <= 0 && return copy(data)

    radius = max(1, ceil(Int, 3 * sigma_px))
    kernel = [exp(-0.5 * (k / sigma_px)^2) for k in -radius:radius]
    kernel ./= sum(kernel)

    dim1, dim2 = size(data)
    src = replace(data, NaN => 0.0)

    # convolution séparable : passage horizontal puis vertical
    tmp = similar(src)
    @inbounds for j in 1:dim2, i in 1:dim1
        s = 0.0
        for k in -radius:radius
            ii = clamp(i + k, 1, dim1)
            s += kernel[k+radius+1] * src[ii, j]
        end
        tmp[i, j] = s
    end

    out = similar(src)
    @inbounds for j in 1:dim2, i in 1:dim1
        s = 0.0
        for k in -radius:radius
            jj = clamp(j + k, 1, dim2)
            s += kernel[k+radius+1] * tmp[i, jj]
        end
        out[i, j] = s
    end

    return out
end

# Calcule radial_intensity/radial_width (moyenne/écart-type par anneau, cf. sunkit-image) pour une
# image 2D (un seul canal, ou une image de luminance combinée), avec lissage de la queue du profil
# (voir smooth_tail_hermite!) pour éviter les cassures de pente dues à la tangence anneau/bord.
function compute_radial_profile(data::AbstractMatrix{T}, bin_idx::Matrix{Int}, nbins::Int,
    r_centers::Vector{Float64};
    r_start_smooth::Union{Float64,Nothing}=nothing,
    r_end_smooth::Union{Float64,Nothing}=nothing,
    n_slope::Int=50) where {T<:AbstractFloat}

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

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

    radial_intensity = fill(NaN, nbins)
    radial_width = fill(NaN, nbins)
    for i in 1:nbins
        if count_valid[i] > 0
            radial_intensity[i] = sum_valid[i] / count_valid[i]
        end
        if count_all[i] > 0
            m = sum_all[i] / count_all[i]
            # écart-type "population" (ddof=0), comme np.std par défaut ; max(...,0) protège contre
            # un léger dépassement négatif dû aux arrondis flottants quand l'écart-type est ~0
            radial_width[i] = sqrt(max(sumsq_all[i] / count_all[i] - m^2, 0.0))
        end
    end

    if r_start_smooth !== nothing && r_end_smooth !== nothing
        smooth_tail_hermite!(radial_intensity, r_centers, r_start_smooth, r_end_smooth; n_slope=n_slope)
        smooth_tail_hermite!(radial_width, r_centers, r_start_smooth, r_end_smooth; n_slope=n_slope)
    end

    return radial_intensity, radial_width
end

# Reconstruit l'image 2D radial_luminosity (chaque pixel = valeur de son anneau), puis lui applique
# le filtre gaussien (voir gaussian_blur ci-dessus).
function reconstruct_and_smooth_luminosity(radial_intensity::Vector{Float64}, bin_idx::Matrix{Int},
    sigma_px::Float64)

    out = Array{Float64}(undef, size(bin_idx))
    @inbounds for idx in eachindex(bin_idx)
        out[idx] = radial_intensity[bin_idx[idx]]
    end
    return gaussian_blur(out, sigma_px)
end

# Applique la normalisation NRGF à un canal, à partir d'une image radial_luminosity déjà
# reconstruite (et lissée) et d'un profil radial_width (par anneau, via bin_idx). C'est la même
# opération que le NRGF d'origine (nouvelle_valeur = (valeur - moyenne_anneau) / écart_type_anneau,
# au-delà du rayon d'application), mais découplée du calcul du profil, pour pouvoir combiner
# librement "quel profil normalise quel canal" (voir division_mode dans main()).
function nrgf_apply(data::AbstractMatrix{T}, luminosity_map::AbstractMatrix{Float64},
    radial_width::Vector{Float64}, bin_idx::Matrix{Int}, map_r::AbstractMatrix{Float64};
    application_radius::Float64=1.0, fill_value::T=T(NaN)) where {T<:AbstractFloat}

    out = fill(fill_value, size(data))
    @inbounds for idx in eachindex(data)
        if map_r[idx] > application_radius
            v = T(data[idx]) - T(luminosity_map[idx])
            w = radial_width[bin_idx[idx]]
            if w != 0.0   # si w est NaN, cette condition est vraie (comme en Python) -> v/NaN = NaN
                v = v / T(w)
            end
            out[idx] = v
        end
    end
    return out
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 (repérage possible avec un afficheur FITS comme SAOImage DS9, Siril,
    # etc.). Convention 1-based (comme les indices du tableau Julia image_in[x, y, :]) : 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 NRGF
    # -----------------------------------------------------------------------------------------------
    #nbins_nrgf = dim_x ÷ 2            # nombre d'anneaux radiaux (valeur par défaut de sunkit-image)
    nbins_nrgf = dim_x ÷ 2
    application_radius_nrgf = 1.0     # en rayons solaires ; NRGF appliqué seulement au-delà
    fill_value_nrgf = Float32(0.0)    # valeur hors de la zone d'application (0.0 plutôt que NaN,
    # par défaut dans sunkit-image, pour faciliter les traitements en aval)

    # -----------------------------------------------------------------------------------------------
    # rayons (en rayons solaires) des deux points de calage de la spline de Hermite qui lisse la
    # queue du profil radial_intensity, pour éviter les cassures de pente dues à la tangence
    # anneau/bord (voir smooth_tail_hermite!) :
    #   - r_start, juste avant le 1er seuil de tangence (anneau touchant Ymax)
    #   - r_end, juste avant le 2e seuil de tangence (anneau touchant Xmax) - plutôt que r_max
    #     (coins de l'image), pour éviter les zones sujettes à des défauts d'uniformité (flat) et
    #     rester dans une zone où la couronne est plus intense
    # Au-delà de r_end, le profil est prolongé en ligne droite (pente conservée) jusqu'à r_max.
    # Même marge de sécurité utilisée dans les deux sens (X et Y) - valeur à ajuster par essai.
    # -----------------------------------------------------------------------------------------------
    marge_avant_tangence_px = 200.0   # <-- À AJUSTER par essai
    r_start_smooth_nrgf = (dim_y / 2 - marge_avant_tangence_px) / r_sun_px
    r_end_smooth_nrgf = (dim_x / 2 - marge_avant_tangence_px) / r_sun_px
    n_slope_smooth_nrgf = 50          # nb de points utilisés pour estimer les pentes de raccord
    println(@sprintf("Lissage de la queue du profil : r_start = %.3f R_sun, r_end = %.3f R_sun",
        r_start_smooth_nrgf, r_end_smooth_nrgf))

    # -----------------------------------------------------------------------------------------------
    # écart-type (en pixels) du filtre gaussien appliqué à l'image radial_luminosity (fond radial
    # reconstruit en 2D) avant utilisation dans le calcul du NRGF - supprime les effets de
    # discrétisation subtils dus à la largeur finie des anneaux (bin_width).
    # -----------------------------------------------------------------------------------------------
    sigma_gaussian_px = 2.5   # <-- À AJUSTER

    # -----------------------------------------------------------------------------------------------
    # choix du "canal diviseur" (image radial_luminosity utilisée pour normaliser chaque canal de
    # sortie) :
    #   1 = :rgb -> chaque canal R,G,B est normalisé par SON PROPRE profil radial (comportement
    #               d'origine du NRGF, indépendant par canal)
    #   2 = :l   -> les 3 canaux sont normalisés par le profil d'une image de LUMINANCE (Rec.709,
    #               combinaison de R,G,B), la même pour les trois
    #   3 = :r   -> les 3 canaux sont normalisés par le profil du seul canal ROUGE
    #   4 = :g   -> les 3 canaux sont normalisés par le profil du seul canal VERT
    #   5 = :b   -> les 3 canaux sont normalisés par le profil du seul canal BLEU
    # Utiliser un profil commun (l/r/g/b) plutôt qu'un profil indépendant par canal (rgb) préserve
    # mieux les rapports de couleur d'origine de la couronne.
    # -----------------------------------------------------------------------------------------------
    println()
    division_mode_nrgf = :rgb   # <-- À CHOISIR parmi :rgb, :l, :r, :g, :b (voir liste ci-dessus)
    @assert division_mode_nrgf in (:rgb, :l, :r, :g, :b) "division_mode_nrgf doit être :rgb, :l, :r, :g ou :b"
    println("Mode choisi : ", division_mode_nrgf)
    flush(stdout)

    println("Calcul des distances au centre solaire pour chaque pixel...")
    flush(stdout)
    map_r = [sqrt((Float64(x) - cx_sun_px)^2 + (Float64(y) - cy_sun_px)^2) / r_sun_px for x in 1:dim_x, y in 1:dim_y]
    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_nrgf
    bin_idx = clamp.(floor.(Int, map_r ./ bin_width) .+ 1, 1, nbins_nrgf)
    r_centers = [(i - 0.5) * bin_width for i in 1:nbins_nrgf]

    # -----------------------------------------------------------------------------------------------
    # calcul du/des profil(s) radial(aux) nécessaire(s) selon le mode choisi (un seul, sauf en mode
    # :rgb où il en faut un par canal)
    # -----------------------------------------------------------------------------------------------
    println("Calcul du/des profil(s) radial(aux)...")
    flush(stdout)

    profiles = Dict{Symbol,Tuple{Vector{Float64},Vector{Float64}}}()

    if division_mode_nrgf == :rgb
        for (i_ch, sym) in zip([1, 2, 3], [:r, :g, :b])
            profiles[sym] = compute_radial_profile(image_in[:, :, i_ch], bin_idx, nbins_nrgf, r_centers;
                r_start_smooth=r_start_smooth_nrgf, r_end_smooth=r_end_smooth_nrgf, n_slope=n_slope_smooth_nrgf)
        end
    elseif division_mode_nrgf == :l
        # luminance Rec.709 standard, calculée à partir des 3 canaux de l'image d'entrée
        image_L = 0.2126f0 .* image_in[:, :, 1] .+ 0.7152f0 .* image_in[:, :, 2] .+ 0.0722f0 .* image_in[:, :, 3]
        profiles[:l] = compute_radial_profile(image_L, bin_idx, nbins_nrgf, r_centers;
            r_start_smooth=r_start_smooth_nrgf, r_end_smooth=r_end_smooth_nrgf, n_slope=n_slope_smooth_nrgf)
    else # :r, :g, :b
        i_ch = division_mode_nrgf == :r ? 1 : (division_mode_nrgf == :g ? 2 : 3)
        profiles[division_mode_nrgf] = compute_radial_profile(image_in[:, :, i_ch], bin_idx, nbins_nrgf, r_centers;
            r_start_smooth=r_start_smooth_nrgf, r_end_smooth=r_end_smooth_nrgf, n_slope=n_slope_smooth_nrgf)
    end

    # -----------------------------------------------------------------------------------------------
    # reconstruction + lissage gaussien de l'image radial_luminosity, pour chaque profil calculé
    # (une seule fois par profil, même si réutilisé pour plusieurs canaux)
    # -----------------------------------------------------------------------------------------------
    println("Reconstruction et lissage gaussien (sigma = ", sigma_gaussian_px, " px) de radial_luminosity...")
    flush(stdout)
    luminosity_maps = Dict{Symbol,Matrix{Float64}}()
    for (sym, prof) in profiles
        radial_intensity, _ = prof
        luminosity_maps[sym] = reconstruct_and_smooth_luminosity(radial_intensity, bin_idx, sigma_gaussian_px)
    end

    # profil (clé du Dict) à utiliser pour chacun des 3 canaux de sortie, selon le mode choisi
    channel_profile_key = division_mode_nrgf == :rgb ? [:r, :g, :b] : fill(division_mode_nrgf, 3)

    image_nrgf = similar(image_in)
    image_radial_luminosity = similar(image_in)

    println("Démarrage du calcul NRGF...")
    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) - profiles/luminosity_maps
    # sont uniquement lus ici (calculés avant la boucle), donc sans risque de conflit entre threads
    # -----------------------------------------------------------------------------------------------
    Threads.@threads for i_ch in [1, 2, 3]
        sym = channel_profile_key[i_ch]
        println("calcul NRGF, canal ", i_ch, " (profil '", sym, "', thread ", Threads.threadid(), ")")
        flush(stdout)
        _, radial_width = profiles[sym]
        lum_map = luminosity_maps[sym]
        image_nrgf[:, :, i_ch] = Base.invokelatest(nrgf_apply, image_in[:, :, i_ch], lum_map, radial_width, bin_idx, map_r;
            application_radius=application_radius_nrgf, fill_value=fill_value_nrgf)
        image_radial_luminosity[:, :, i_ch] = Float32.(lum_map)
        println("-> canal ", i_ch, " terminé")
        flush(stdout)
    end

    # nom des fichiers de sortie : ajout de "_NRGF"/"_radial-luminosity" juste avant l'extension
    # ".fit", du nombre d'anneaux (nb), du rayon d'application x10 sur 2 chiffres (ar), et du
    # suffixe correspondant au mode de normalisation choisi (_rgb, _l, _r, _g ou _b)
    params_str = @sprintf("nb%d_ar%02d_", nbins_nrgf, round(Int, application_radius_nrgf * 10))
    suffix_mode = String(division_mode_nrgf)   # "rgb", "l", "r", "g" ou "b"

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

    # fichier radial-luminosity : fond radial (lissé) reconstruit en 2D, effectivement utilisé pour
    # normaliser chaque canal - utile pour vérifier visuellement l'absence de cassure de pente
    file_name_out_lum = chop(file_name_in, tail=4) * "_" * params_str * "radial-luminosity_" * suffix_mode * ".fit"
    f_out = FITS(file_name_out_lum, "w")
    write(f_out, image_radial_luminosity[:, :, :]; header=header)
    close(f_out)
    println("fichier écrit : ", file_name_out_lum)

    # 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)