using Pkg

# -----------------------------------------------------------------------------------------------
# Installation des packages : mettre à false une fois les packages installés une première fois.
# Pkg.add() vérifie le registre à chaque exécution, ce qui coûte plusieurs secondes/minutes
# à chaque lancement du script pour rien si les packages sont déjà présents.
# -----------------------------------------------------------------------------------------------
const FIRST_RUN = false
#const FIRST_RUN = true

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

using Statistics
using Printf

# calcul algébrique, pour calcul angle entre deux vecteurs
# (retiré : LinearAlgebra n'est pas utilisé dans ce script)

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

# package pour filtrage des images (imfilter, Kernel.gaussian)
using Images

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

# package pour l'affichage du nuage de points ADU rouge vs vert (s'affiche dans l'onglet
# "Plots" de VSCode, qui joue le rôle de fenêtre graphique dédiée)
using Plots

# Threads.@threads : fait partie de Base, pas besoin de Pkg.add / using séparé.
# Important : le nombre de threads utilisés est fixé au lancement de Julia, pas dans le script.
# Lancer Julia avec, par exemple : julia --threads=4 nom_du_script.jl
# (ou variable d'environnement JULIA_NUM_THREADS=4 avant de lancer Julia / VSCode)

##################################################################################################################
# Fonctions du programme:
#   read the HDR image (32-bit fit file, 3 color)
#   read the HDR_starless image (32-bit fit file, 3 color)
#   mask = HDR_starless convoluted by gaussian
#   result = HDR / (mask + offset)
#   write result in 32-fit   
#   => bon réglage pour traitement avec quatre valeurs de sigma = 30,20, 10, 5 et division avec constante
#
# Entrées/paramètres :
#      directory et nom de l'image HDR
#      image HDR starless (postfixe _starless, son nom est construit automatiquement)
#      valeur des différents sigma des gaussiennes et seuils ADU associés (en dur dans chaque fonction de pondération)
#      valeur de l'offset

#      maj pour TSE 2026
#           seuils ADU
#           ajout convolution partielle (le mask est calculé dans PS : noir sur disque lunaire, blanc hors disque lunaire)
#           claude : optimisation temps de calcul (main(), @views)
#           choix directory et fichier entrée par IHM : attention, penser à créer le fichier starless associé
#           boolean en début de script pour charger ou pas les packages
#           parallélisation (Threads.@threads) de la boucle sur les 4 sigmas + chronométrage
#      v9 :
#           suppression de l'ancienne méthode de pondération (4 fonctions FW_TPZ_sigma_1..4, 14 seuils
#             ADU) au profit d'une unique fonction générique à 6 points d'ancrage (adu_breakpoints_optionA)
#           ajout de la correction colorimétrique (canaux R et B recalés sur G via régression linéaire,
#             fichier de sortie "..._color.fit")
#
# TSE2026-3-multi-gaussian-mask-v1 (ce script) :
#      ajout du choix read_moon_mask = OUI/NON (au lieu de toujours lire "mask_small.fit") :
#        - OUI (true)  : comportement inchangé, lecture du fichier fit contenant le masque
#        - NON (false) : masque calculé directement dans le script, même principe que le masque
#          lunaire de radial-unsharp-v2 (centre X_centre_moon_mask/Y_centre_moon_mask, demi-axes r_mask_h
#          (horizontal) et r_mask_v (vertical) d'une ellipse, en dur) : ADU = 0 à l'intérieur de
#          l'ellipse, ADU = 1 à l'extérieur (le masque brut n'est pas utilisé pour une division
#          directe, seule sa version floutée l'est, déjà protégée contre les 0)
#      ajout d'une carte de contrôle "..._R_map.fit" (ellipse (r_mask_h, r_mask_v) surimprimée sur l'image
#        de sortie finale Image_div_f, épaisseur circle_thickness_px réglable) - générée seulement
#        si le masque est calculé (read_moon_mask=false), la géométrie n'étant pas connue sinon
###################################################################################################################

function FW_TPZ(x,ADU_1,ADU_2,ADU_3,ADU_4)

    # Fonctions :
    #   fonction de ponderation en forme de trapèze s'appliquant à HDR/masq gaussien (après normalisation)
    
    # Entrée :
    #   x = intensité du pixel
    #   ADU_1 à 4 : seuils ADU (du plus bas au plus haut) définissant les bornes du trapèze
    #              
    # Sortie :
    #   weight_coeff : de 0 à 1

    weight_coeff=0.
    
    if (ADU_1 < 100) && (x < ADU_3)
        #courbe horizontale jusqu'à ADU_3
            weight_coeff=1
    else
        if (x < ADU_1)
            weight_coeff=0
        elseif (x < ADU_2)
            #montée linéaire de ADU_1 à ADU_2
            weight_coeff=(x-ADU_1)/(ADU_2-ADU_1)
        elseif (x < ADU_3)
            weight_coeff=1
        elseif (ADU_4 < 100)
            weight_coeff=1
        elseif (x < ADU_4)
            #descente linéaire de ADU_3 à ADU_4
            weight_coeff=1-(x-ADU_3)/(ADU_4-ADU_3)
        else
            weight_coeff=0
        end
              
        # suppression des valeurs négatives ou nulles pour éviter division par 0
        weight_coeff=max(weight_coeff,0.0001)
    end

    return weight_coeff
end

#############################################################################################
# calcul du poids d'un masque gaussien (sigma[i_sigma]) en fonction de l'intensité locale du
# pixel, via UNE fonction générique paramétrée par 6 points d'ancrage (au lieu de 4 fonctions
# séparées à 14 seuils ADU au total, à maintenir manuellement cohérents entre elles - historique :
# voir versions v8c et antérieures). Le poids du masque i_sigma est évalué sur SA PROPRE valeur
# floutée Image_masq[:,:,k,i_sigma] (voir appel dans main()) : chaque masque "vote" sur la base de
# sa propre intensité locale, pas d'une intensité de référence commune à tous les masques.
#
# Les 6 points d'ancrage (adu_breakpoints, défini dans main()) correspondent aux limites de chaque
# plateau et de chaque transition, dans l'ordre croissant d'ADU :
#   breakpoints[1] = fin du plateau "sigma_1 (30px) pur"          / début transition sigma_1↔sigma_2
#   breakpoints[2] = fin de la transition sigma_1↔sigma_2         / début du plateau "sigma_2 pur"
#   breakpoints[3] = fin du plateau "sigma_2 (20px) pur"          / début transition sigma_2↔sigma_3
#   breakpoints[4] = fin de la transition sigma_2↔sigma_3         / début du plateau "sigma_3 pur"
#   breakpoints[5] = fin du plateau "sigma_3 (10px) pur"          / début transition sigma_3↔sigma_4
#   breakpoints[6] = fin de la transition sigma_3↔sigma_4         / début du plateau "sigma_4 pur"
#                    (plateau qui continue jusqu'à l'infini : sigma_4 est le masque le plus fin,
#                    il n'y a pas de masque encore plus fin après lui pour prendre le relais)
# Avec les valeurs par défaut : [3000, 10000, 20000, 50000, 90000, 150000].
#############################################################################################

# Reconstitue, pour le masque n° i_sigma (parmi nb_sigma), le quadruplet (ADU_1,ADU_2,ADU_3,ADU_4)
# attendu par FW_TPZ, à partir des 6 points d'ancrage communs - gère les cas particuliers du
# premier masque (plateau=1 depuis l'origine, pas de "montée") et du dernier masque (plateau=1
# jusqu'à l'infini, pas de "descente").
function sigma_thresholds_optionA(i_sigma::Int, nb_sigma::Int, breakpoints::Vector{Float64})
    if i_sigma == 1
        # 1er masque (sigma le plus grand) : ADU_1<100 déclenche le plateau=1 depuis l'origine
        # (voir FW_TPZ) ; ADU_2 n'intervient pas dans ce cas, seuls ADU_3/ADU_4 comptent (descente)
        return (1.0, 2.0, breakpoints[1], breakpoints[2])
    elseif i_sigma == nb_sigma
        # dernier masque (sigma le plus petit) : montée entre les 2 derniers points d'ancrage, puis
        # ADU_4<100 déclenche le plateau=1 au-delà de ADU_3 sans jamais redescendre (voir FW_TPZ)
        return (breakpoints[end-1], breakpoints[end], 1.0e12, 1.0)
    else
        # masques intermédiaires : 4 points d'ancrage consécutifs (montée, plateau, descente)
        j = 2 * i_sigma - 3
        return (breakpoints[j], breakpoints[j+1], breakpoints[j+2], breakpoints[j+3])
    end
end

function FW_TPZ_general(x, i_sigma::Int, nb_sigma::Int, breakpoints::Vector{Float64})
    ADU_1, ADU_2, ADU_3, ADU_4 = sigma_thresholds_optionA(i_sigma, nb_sigma, breakpoints)
    return FW_TPZ(x, ADU_1, ADU_2, ADU_3, ADU_4)
end

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

# tout le programme principal est placé dans une fonction : en Julia, le code
# écrit directement au niveau global est beaucoup plus lent (variables non typées).
# Une fonction est compilée avec des types connus => gain de vitesse important, sans
# changer la logique du programme.
function main()

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

# =================================================================================================
# PARAMÈTRES UTILISATEUR — à ajuster selon la run (regroupés ici, cf. discussion sur les bonnes
# pratiques : ces valeurs ne dépendent d'aucune donnée chargée plus loin dans le script)
# =================================================================================================

# convolution partielle évite les rebonds au limbe lunaire
partial_conv = true

# écriture optionnelle des fichiers masque gaussien (file_name_masq_gauss[i_sigma]) : ces fichiers
# ne servent que pour des tests/vérifications visuelles, pas au calcul lui-même => écriture désactivable
write_gaussian_masque = false

# écriture optionnelle des fichiers HDR/masque gaussien (file_name_div_gauss[i_sigma]) : ces fichiers
# ne servent que pour des tests/vérifications visuelles, pas au calcul lui-même => écriture désactivable
write_div_gauss = false

# -----------------------------------------------------------------------------------------------
# choix read_moon_mask = OUI(true)/NON(false) : lire le masque depuis un fichier fit
# (comportement d'origine), ou le calculer directement dans le script (voir en-tête du script,
# même principe que le masque lunaire de radial-unsharp-v2).
# il s'agit du masque définissant le disque lunaire, utilisé dans le calcul de la convolution partielle
# -----------------------------------------------------------------------------------------------
read_moon_mask = false

# géométrie du masque lunaire (utilisée seulement si read_moon_mask = false), en dur, modifiable :
# ellipse définie par son centre et ses 2 demi-axes (horizontal, vertical)
X_centre_moon_mask = 4055.0
Y_centre_moon_mask = 2698.0
r_mask_h = 655.0   # demi-grand axe horizontal (ellipse)
r_mask_v = 654.0   # demi-grand axe vertical (ellipse)

# ADU à l'intérieur de l'ellipse masquée, si calculé : 0 (le masque brut n'est pas utilisé pour une
# division directe - seule sa version floutée Image_mask_conv l'est, déjà protégée contre les
# valeurs proches de 0, voir plus bas)
valeur_disque_mask = 0.0

# -----------------------------------------------------------------------------------------------
# carte de contrôle (debug) : écrit un fichier supplémentaire "..._R_map.fit", copie de l'image
# de sortie finale (Image_div_f, avant correction colorimétrique) avec l'ellipse (r_mask_h,
# r_mask_v) tracée en surimpression (valeur = ADU max de l'image), pour vérifier visuellement
# l'alignement du masque calculé avec la Lune réelle. Nécessite read_moon_mask=false
# (géométrie du masque non connue si lu depuis un fichier). Générée plus bas, une fois Image_div_f
# calculée.
# -----------------------------------------------------------------------------------------------
export_Rmap = false
circle_thickness_px = 2   # épaisseur du cercle tracé, en pixels

# nombre de masques/sigmas calculés, et valeur (en pixels) de chaque sigma gaussien
nb_sigma = 4
sigma = zeros(nb_sigma)
sigma[1] = 30
sigma[2] = 20
sigma[3] = 10
sigma[4] = 5

# -----------------------------------------------------------------------------------------------
# points d'ancrage définissant le poids de chaque masque gaussien en fonction de l'intensité locale
# (voir commentaire détaillé au-dessus de FW_TPZ_general/sigma_thresholds_optionA)
# -----------------------------------------------------------------------------------------------
adu_breakpoints_optionA = [3000.0, 10000.0, 20000.0, 50000.0, 90000.0, 150000.0]
@assert length(adu_breakpoints_optionA) == 2 * (nb_sigma - 1) "adu_breakpoints_optionA doit contenir 2*(nb_sigma-1) valeurs"

# réglage la force du masque par ajout d'un offset au masque : 0 en noir&blanc, 250 en couleur
couleur = false
if couleur
    c_offset = 250
else
    c_offset = 0
end

println("Nombre de threads Julia disponibles : ", Threads.nthreads())
if Threads.nthreads() == 1
    println("(pour utiliser plusieurs threads, relancer avec : julia --threads=4 nom_du_script.jl)")
end

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

println(pwd())

# nom du fichier HDR en entrée : sélection interactive (remplace le nom de fichier en dur)
selected_hdr = pick_file(work_dir; filterlist="fit,fits")
isempty(selected_hdr) && error("Aucun fichier sélectionné - programme arrêté")
file_name_HDR = basename(selected_hdr)

# lecture de l'image HDR (fit)
f_in_hdr=FITS(file_name_HDR)
# read complete header
header = read_header(f_in_hdr[1])
# lecture de l'image HDR (tableau avec trois couches RGB : X,Y,couleurs)
image_hdr=read(f_in_hdr[1])

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

# stat sur l'image HDR
println("max R = ",maximum(image_hdr[:,:,1]))
println("max V = ",maximum(image_hdr[:,:,2]))
println("max B = ",maximum(image_hdr[:,:,3]))
println("mean R = ",mean(image_hdr[:,:,1]))
println("mean G = ",mean(image_hdr[:,:,2]))
println("mean B = ",mean(image_hdr[:,:,3]))

println("median R = ",median(image_hdr[:,:,1]))

# analyse d'histogramme par canal : valeur ADU en-dessous de laquelle se trouve
# le X% de pixels les plus brillants (comptage du plus brillant vers le plus sombre)
# ex : "5% max ADU" = valeur ADU du 95e centile (95% des pixels sont plus sombres qu'elle)
# stockage dans des variables nommées explicites p5_/p10_/p50_/p90_ + canal (R,G,B)

pixels_R=vec(image_hdr[:,:,1])
pixels_G=vec(image_hdr[:,:,2])
pixels_B=vec(image_hdr[:,:,3])

p5_R=quantile(pixels_R,0.95)
p10_R=quantile(pixels_R,0.90)
p50_R=quantile(pixels_R,0.50)
p90_R=quantile(pixels_R,0.10)

p5_G=quantile(pixels_G,0.95)
p10_G=quantile(pixels_G,0.90)
p50_G=quantile(pixels_G,0.50)
p90_G=quantile(pixels_G,0.10)

p5_B=quantile(pixels_B,0.95)
p10_B=quantile(pixels_B,0.90)
p50_B=quantile(pixels_B,0.50)
p90_B=quantile(pixels_B,0.10)

println(@sprintf(" 5%% max ADU (R) = %.2f",p5_R))
println(@sprintf("10%% max ADU (R) = %.2f",p10_R))
println(@sprintf("50%% max ADU (R) = %.2f",p50_R))
println(@sprintf("90%% max ADU (R) = %.2f",p90_R))

println(@sprintf(" 5%% max ADU (G) = %.2f",p5_G))
println(@sprintf("10%% max ADU (G) = %.2f",p10_G))
println(@sprintf("50%% max ADU (G) = %.2f",p50_G))
println(@sprintf("90%% max ADU (G) = %.2f",p90_G))

println(@sprintf(" 5%% max ADU (B) = %.2f",p5_B))
println(@sprintf("10%% max ADU (B) = %.2f",p10_B))
println(@sprintf("50%% max ADU (B) = %.2f",p50_B))
println(@sprintf("90%% max ADU (B) = %.2f",p90_B))

# nom du fichier HDR starless en entrée
# chop = suppression des 4 derniers char du fichier HDR (.fit) + ajout extension
file_name_HDR_starless=chop(file_name_HDR,tail=4)*"_starless.fit"

# lecture de l'image HDR_starless (fit)
f_in_hdr_starless=FITS(file_name_HDR_starless)
# lecture de l'image HDR (tableau avec trois couches RGB : X,Y,couleurs)
image_hdr_starless=read(f_in_hdr_starless[1])

if (partial_conv==true) 

    println("read_moon_mask = ", read_moon_mask)
    println("export_Rmap = ", export_Rmap, " (épaisseur cercle = ", circle_thickness_px, " px)")
    if export_Rmap && read_moon_mask
        println("ATTENTION : export_Rmap=true nécessite read_moon_mask=false (géométrie du masque connue) - carte de contrôle non générée")
    end

    if read_moon_mask
        println("calcul des masques avec déconvolution partielle pour éviter les rebonds (masque lu depuis fichier)")
        # lecture du masque de deconvolution partielle
        # nom du masque déconvolution partielle
        file_name_mask="mask_small.fit"

        # lecture du masque de déconvolution partielle (fit)
        f_in_mask=FITS(file_name_mask)
        # lecture du masque de déconvolution partielle (tableau avec trois couches RGB : X,Y,couleurs)
        image_mask=Float64.(read(f_in_mask[1]))

        # normalisation à [0,1] du masque de déconvolution partielle (les trois canaux sont identiques)
        mask_norm=maximum(image_mask[:,:,3])
        image_mask[:,:,:]=image_mask[:,:,:]/mask_norm

        println("max mask R = ",maximum(image_mask[:,:,1]))

    else
        println("calcul des masques avec déconvolution partielle pour éviter les rebonds (masque calculé dans le script)")
        println("géométrie du masque : centre = (", X_centre_moon_mask, ",", Y_centre_moon_mask,
                "), r_mask_h = ", r_mask_h, " px, r_mask_v = ", r_mask_v, " px (ellipse)")

        # masque calculé : valeur_disque_mask (0) à l'intérieur de l'ellipse, 1 à l'extérieur (3
        # canaux identiques, comme le masque lu depuis fichier, déjà normalisé - pas besoin de
        # division par un max). Équation de l'ellipse : ((i-Xc)/r_mask_h)^2 + ((j-Yc)/r_mask_v)^2 <= 1
        image_mask = ones(dim_x, dim_y, 3)
        @inbounds for j in 1:dim_y, i in 1:dim_x
            if ((i - X_centre_moon_mask)/r_mask_h)^2 + ((j - Y_centre_moon_mask)/r_mask_v)^2 <= 1.0
                image_mask[i,j,1] = valeur_disque_mask
                image_mask[i,j,2] = valeur_disque_mask
                image_mask[i,j,3] = valeur_disque_mask
            end
        end

        println("max mask R = ",maximum(image_mask[:,:,1]))
    end

end

# tableau contenant le nom des fichiers masque et div gauss dépendant de sigma
file_name_masq_gauss = Array{String}(undef,nb_sigma)
file_name_div_gauss = Array{String}(undef,nb_sigma)

for i_sigma=1:nb_sigma
    # nom du fichier masque gaussien - sigma sur 2 caractères
    file_name_masq_gauss[i_sigma]=chop(file_name_HDR,tail=4)*"_masq_gauss_"*@sprintf("%02i",sigma[i_sigma])*"px.fit"
    # nom du fichier HDR/masq_gauss_1
    file_name_div_gauss[i_sigma]=chop(file_name_HDR,tail=4)*"_div_gauss_"*@sprintf("%02i",sigma[i_sigma])*"px.fit"

    if (partial_conv==true)
        # ajout du postfixe "_conv" juste avant l'extension ".fit"
        file_name_masq_gauss[i_sigma]=chop(file_name_masq_gauss[i_sigma],tail=4)*"_conv.fit"
        file_name_div_gauss[i_sigma]=chop(file_name_div_gauss[i_sigma],tail=4)*"_conv.fit"
    end
end

#création et initialisation à 0 des images masques gaussiens dépendant de sigma (Float64, 3 couleurs, nb_sigma)     
Image_masq=zeros(dim_x,dim_y,3,nb_sigma)
#création et initialisation à 0 des images HDR/masq dépendant de sigma (Float64, 3 couleurs)
Image_div=zeros(dim_x,dim_y,3,nb_sigma)

#création et initialisation à 0 des images mask gaussiens dépendant de sigma (Float64, 3 couleurs, nb_sigma) => utilisés pour déconvolution partielle (Lune noir, reste = 1)
Image_mask_conv=zeros(dim_x,dim_y,3,nb_sigma)

# chaîne "sigXX-XX-XX-XX" construite à partir des valeurs de sigma (2 digits chacune)
sig_str="sig"*join([@sprintf("%02i",sigma[i_sigma]) for i_sigma=1:nb_sigma],"-")

# nom du masque final obtenu par somme pondérée des masques et nom de l'image image div (gauss + offset) finale associée
file_name_masq_gauss_f=chop(file_name_HDR, tail=4)*"_masq_gauss_"*sig_str*"_final_offset"*@sprintf("%04i",c_offset)*".fit"
file_name_div_gauss_f=chop(file_name_HDR, tail=4)*"_div_gauss_"*sig_str*"_final_offset"*@sprintf("%04i",c_offset)*".fit"

if (partial_conv==true)
    # ajout du postfixe "_conv" devant l'extension ".fit"
    file_name_masq_gauss_f=chop(file_name_masq_gauss_f,tail=4)*"_conv.fit"
    file_name_div_gauss_f=chop(file_name_div_gauss_f,tail=4)*"_conv.fit"
end

#création et initialisation à 0 des images associées (Float64, 3 couleurs)  
Image_masq_f=zeros(dim_x,dim_y,3)
Image_div_f=zeros(dim_x,dim_y,3)

#création et initialisation à 0 de du cumul des coefficients de pondération (Float64, 3 couleurs)
W_map_C_cumul=zeros(dim_x,dim_y,3)

#.......................................................................................................................

# calcul des différents masques gaussiens
# sur les kernel : https://juliaimages.org/ImageFiltering.jl/dev/kernels/
# https://juliaimages.org/ImageFiltering.jl/latest/function_reference.html#KernelFactors-1

# verrou pour sérialiser les écritures FITS : FITSIO/CFITSIO (bibliothèque C sous-jacente)
# n'est pas garanti thread-safe pour des écritures simultanées. Le calcul lourd (imfilter)
# reste, lui, parallélisé entre threads - seule l'écriture fichier est sérialisée (rapide).
io_lock = ReentrantLock()

Threads.@threads for i_sigma=1:nb_sigma

    println("[thread ", Threads.threadid(), "] calcul masque gaussien, sigma = ", sigma[i_sigma])

    # on utilise ici l'image HDR_starless pour éviter les artefacts de type rebond noir autour des étoiles

    # le kernel gaussien ne dépend que de i_sigma : on le construit une seule fois
    # au lieu de le refaire pour chaque canal k=1,2,3 (identique à chaque fois)
    gauss_kernel = Kernel.gaussian(sigma[i_sigma])

    # @views évite de copier chaque canal couleur (image_hdr_starless[:,:,k], image_mask[:,:,k])
    # avant de le filtrer : imfilter travaille directement sur la tranche d'origine => moins d'allocations
    @views for k in [1,2,3]
        # calcul de la gaussienne pour la couche courante - résultat en Float64
        # par défaut, la taille du kernel est égale à deux fois sigma (?) => 25 valeurs pour sigma = 1
        Image_masq[:,:,k,i_sigma]=imfilter(image_hdr_starless[:,:,k], gauss_kernel)      
    
#       prise en compte éventuelle de la déconvolution partielle dans le calcul du masque
        if (partial_conv==true) 
          # calcul du mask [0,1] convolué par la gaussienne courante
             Image_mask_conv[:,:,k,i_sigma]=imfilter(image_mask[:,:,k], gauss_kernel)   
         
            # remplacement des valeurs proches de 0 (< epsilon) par la valeur 1, pour éviter division par 0
             epsilon=1e-6
             Image_mask_conv[:,:,k,i_sigma]=ifelse.(Image_mask_conv[:,:,k,i_sigma] .< epsilon, 1.0, Image_mask_conv[:,:,k,i_sigma])

             # calcul du masque modifié par prise en compte convolution partielle
             Image_masq[:,:,k,i_sigma]=Image_masq[:,:,k,i_sigma]./Image_mask_conv[:,:,k,i_sigma]
        end
    end

    # creation du fichier masque gaussienne FIT (écriture sérialisée via le verrou io_lock)
    # écriture optionnelle, pilotée par write_gaussian_masque (fichiers utilisés seulement pour des tests)
    if write_gaussian_masque
        lock(io_lock) do
            local f_out
            f_out = FITS(file_name_masq_gauss[i_sigma], "w")
            # le header initial complète le nouveau header, infos récupérées de l'image initiale = date  et directory
            write(f_out,Image_masq[:,:,:,i_sigma];header=header)
            close(f_out)
        end
    end
 
    println("[thread ", Threads.threadid(), "] division par masque gaussien, sigma = ", sigma[i_sigma])
    # calcul de l'image HDR divisée par le masque gaussien, et résultat en Float64
    # normalisation sur niv moyen de chaque couleur de l'image HDR

    for k in [1,2,3]
        local k_norm
        k_norm=mean(image_hdr[:,:,k])
        Image_div[:,:,k,i_sigma]=k_norm*image_hdr[:,:,k]./Image_masq[:,:,k,i_sigma]
    end

    # écriture de l'image HDR divisée par le masque gaussien (écriture sérialisée via le verrou io_lock)
    # écriture optionnelle, pilotée par write_div_gauss (fichiers utilisés seulement pour des tests)
    if write_div_gauss
        lock(io_lock) do
            f_out = FITS(file_name_div_gauss[i_sigma], "w")
            write(f_out,Image_div[:,:,:,i_sigma];header=header)
            close(f_out)
        end
    end

end

#................................................................................................................

println("calcul du masque gaussien final avec sigma variable")

#création et initialisation à 0 du masque courant filtré par les coeff de pondération (Float64, 3 couleurs)
# alloué une seule fois en dehors de la boucle (au lieu de zeros(...) à chaque tour) => moins de pression sur le GC
image_W=zeros(dim_x,dim_y,3)
#création et initialisation à 0 de la carte courante des coefficients de pondération (Float64, 3 couleurs)
W_map_C=zeros(dim_x,dim_y,3)

# boucle sur le nombre de masques à pondérer
for i_sigma=1:nb_sigma

    println("indice du sigma courant = ",i_sigma)

    # calcul de chaque couche filtrée par les coefficients de pondération
    for k in [1,2,3]

        # carte des coeff de pondération pour le masque i_sigma courant (fonction générique à
        # 6 points d'ancrage - voir commentaire au-dessus de FW_TPZ_general)
        W_map_C[:,:,k] = FW_TPZ_general.(Image_masq[:,:,k,i_sigma], i_sigma, nb_sigma, Ref(adu_breakpoints_optionA))
        
        # masque courant filtré pour valeur courante sigma
        # le ".*" permet de faire une multiplication pixel par pixel
        image_W[:,:,k]=W_map_C[:,:,k].*Image_masq[:,:,k,i_sigma]

        # cumul (non normé) des masques filtrés
        Image_masq_f[:,:,k]=Image_masq_f[:,:,k] + image_W[:,:,k]

        # cumul de la carte des coeff de pondération 
        W_map_C_cumul[:,:,k]=W_map_C_cumul[:,:,k]+W_map_C[:,:,k]
    end

end

#..................................................................................................
# calcul du masque final pondéré et normé
for k in [1,2,3] 
    Image_masq_f[:,:,k]=Image_masq_f[:,:,k]./W_map_C_cumul[:,:,k]
end

# écriture du fichier FIT avec masque final pondéré et normé (fichier de test optionnel)
if write_gaussian_masque
    f_out = FITS(file_name_masq_gauss_f, "w")
    write(f_out,Image_masq_f[:,:,:];header=header)
    close(f_out)
end

println("division image HDR par masque gaussien final")
# calcul de l'image HDR divisée par le masque gaussien final, et résultat en Float64
# normalisation sur niv moyen de chaque couleur de l'image HDR

for k in [1,2,3]
    local k_norm
    k_norm=mean(image_hdr[:,:,k])

    # balance des couleurs réglées à la main
    #k_norm=26000
    #if k==3
    #    k_norm=27000
    #end

    # important : espace devant .+
    Image_div_f[:,:,k]=k_norm*image_hdr[:,:,k]./(Image_masq_f[:,:,k] .+c_offset)
end

# écriture de l'image HDR divisée par le masque gaussien final
f_out = FITS(file_name_div_gauss_f, "w")
write(f_out,Image_div_f[:,:,:];header=header)
close(f_out)

#..................................................................................................
# carte de contrôle (debug) : ellipse (r_mask_h, r_mask_v) surimprimée sur l'image de sortie finale
# (Image_div_f), pour vérifier visuellement l'alignement du masque calculé avec la Lune réelle
#..................................................................................................
if partial_conv && export_Rmap && !read_moon_mask
    println("génération de la carte de contrôle R_map (ellipse r_mask_h x r_mask_v)")

    image_Rmap = copy(Image_div_f)
    max_val = maximum(Image_div_f)

    # épaisseur du tracé, exprimée en unités "normalisées" de l'ellipse (rho=1 sur le pourtour) :
    # approximation à partir du rayon moyen, pour une épaisseur visuelle à peu près constante
    r_mask_avg = (r_mask_h + r_mask_v) / 2
    half_thick_rho = (circle_thickness_px / 2) / r_mask_avg

    @inbounds for j in 1:dim_y, i in 1:dim_x
        rho = sqrt(((i - X_centre_moon_mask)/r_mask_h)^2 + ((j - Y_centre_moon_mask)/r_mask_v)^2)
        if abs(rho - 1.0) <= half_thick_rho
            image_Rmap[i,j,1] = max_val
            image_Rmap[i,j,2] = max_val
            image_Rmap[i,j,3] = max_val
        end
    end

    file_name_Rmap = chop(file_name_HDR,tail=4)*"_R_map.fit"
    println("écriture du fichier de contrôle : ", file_name_Rmap)
    f_out = FITS(file_name_Rmap, "w")
    write(f_out, image_Rmap[:,:,:]; header=header)
    close(f_out)
end

#..................................................................................................
# nuages de points ADU rouge (k=1) et ADU bleu (k=3) en fonction de ADU vert (k=2), pour visualiser
# une éventuelle relation (linéaire ou autre) entre canaux, sur l'image finale Image_div_f

println("affichage des nuages de points ADU rouge/bleu vs vert")

# sous-échantillonnage aléatoire : dim_x*dim_y pixels seraient beaucoup trop nombreux à afficher
# (nuage illisible + rendu très lent) ; n_sample points suffisent pour voir la tendance
n_pixels = dim_x*dim_y
n_sample = min(500_000, n_pixels)
idx_sample = rand(1:n_pixels, n_sample)

red_values   = vec(Image_div_f[:,:,1])[idx_sample]
green_values = vec(Image_div_f[:,:,2])[idx_sample]
blue_values  = vec(Image_div_f[:,:,3])[idx_sample]

# régression linéaire simple (moindres carrés) : y = a*x + b
function linreg(x, y)
    a = cov(x, y) / var(x)
    b = mean(y) - a*mean(x)
    return a, b
end

a_RG, b_RG = linreg(green_values, red_values)
a_BG, b_BG = linreg(green_values, blue_values)

println(@sprintf("Régression ADU rouge = %.4f * ADU_vert + %.4f", a_RG, b_RG))
println(@sprintf("Régression ADU bleu  = %.4f * ADU_vert + %.4f", a_BG, b_BG))

#..................................................................................................
# correction colorimétrique : on applique la pente ET l'ordonnée à l'origine de la régression
# pour que R_corr et B_corr coïncident avec G (pente 1, décalage 0). Seuls les canaux R et B
# sont modifiés, le canal G reste inchangé.
# R = a_RG*G + b_RG   =>   R_corr = (R - b_RG) / a_RG   =>   R_corr = G (pente 1, offset 0)
# B = a_BG*G + b_BG   =>   B_corr = (B - b_BG) / a_BG   =>   B_corr = G (pente 1, offset 0)

println("calcul de l'image avec correction colorimétrique (pente et décalage R/G et B/G ramenés à 1 et 0)")

Image_div_f_color = copy(Image_div_f)
Image_div_f_color[:,:,1] = (Image_div_f[:,:,1] .- b_RG) ./ a_RG
Image_div_f_color[:,:,3] = (Image_div_f[:,:,3] .- b_BG) ./ a_BG
# canal vert (k=2) inchangé

# nom du fichier de sortie : ajout de "_color" juste avant l'extension ".fit"
file_name_div_gauss_f_color = chop(file_name_div_gauss_f,tail=4)*"_color.fit"

f_out = FITS(file_name_div_gauss_f_color, "w")
write(f_out,Image_div_f_color[:,:,:];header=header)
close(f_out)

gr()   # backend GR : rendu rapide, compatible avec l'onglet "Plots" de VSCode

x_range = [minimum(green_values), maximum(green_values)]

# --- graphique rouge vs vert ---
plt_RG = scatter(green_values, red_values;
    xlabel = "ADU vert (k=2)",
    ylabel = "ADU rouge (k=1)",
    title  = @sprintf("ADU rouge vs vert  -  y = %.4f x + %.4f", a_RG, b_RG),
    markersize = 1.5,
    markerstrokewidth = 0,
    alpha = 0.25,
    label = "pixels (échantillon $(n_sample)/$(n_pixels))",
    legend = :topleft)

plot!(plt_RG, x_range, a_RG.*x_range .+ b_RG;
    linewidth = 2, color = :black, label = "régression linéaire")

display(plt_RG)

# --- graphique bleu vs vert ---
plt_BG = scatter(green_values, blue_values;
    xlabel = "ADU vert (k=2)",
    ylabel = "ADU bleu (k=3)",
    title  = @sprintf("ADU bleu vs vert  -  y = %.4f x + %.4f", a_BG, b_BG),
    markersize = 1.5,
    markerstrokewidth = 0,
    alpha = 0.25,
    label = "pixels (échantillon $(n_sample)/$(n_pixels))",
    legend = :topleft)

plot!(plt_BG, x_range, a_BG.*x_range .+ b_BG;
    linewidth = 2, color = :black, label = "régression linéaire")

display(plt_BG)

# pour enregistrer les graphiques en image, décommenter les lignes suivantes :
# savefig(plt_RG, chop(file_name_HDR,tail=4)*"_ADU_rouge_vs_vert.png")
# savefig(plt_BG, chop(file_name_HDR,tail=4)*"_ADU_bleu_vs_vert.png")

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

# lancement du programme
main()