using Pkg

# -----------------------------------------------------------------------------------------------
# Installation des packages : mettre à true au premier lancement, puis à false une fois les
# packages installés (Pkg.add() vérifie le registre à chaque exécution, ce qui coûte du temps).
# -----------------------------------------------------------------------------------------------
const FIRST_RUN = false
#const FIRST_RUN = true

if FIRST_RUN
    for pkg in ["Statistics", "Printf", "NativeFileDialog", "TiffImages", "ColorTypes",
                "FixedPointNumbers", "FFTW", "Plots"]
        Pkg.add(pkg)
    end
end

using Statistics
using Printf

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

# lecture / écriture des TIFF 16 bits ; "import" pour éviter les conflits de noms avec Plots
import TiffImages
import ColorTypes
import FixedPointNumbers

# convolution par FFT (méthode :psf)
using FFTW      # multi-thread automatique si Julia est lancé avec plusieurs threads

# graphiques de contrôle (profil photométrique avant / après, carte de la lumière diffusée)
using Plots

##################################################################################################################
# ATTÉNUATION DE LA LUMIÈRE DIFFUSÉE PAR LE CROISSANT LUNAIRE (images de lumière cendrée)
#
# Version 1.3 (14/10/2026)
#   v1.3 : paramètre SAVE_PNG : enregistrement optionnel des graphiques d'analyse (PNG)
#   v1.3 : lecture des fichiers RAW Nikon NEF (D5300, Z 6II, Z 7II...), convertis en TIFF 16 bits
#          linéaire par dcraw_emu (LibRaw), comme dans register_star_earthshine v2.0 (voir "FICHIERS
#          RAW NIKON (NEF)" ci-dessous). TIFF et NEF peuvent être mélangés dans une même sélection.
#   v1.0 : méthodes :radial, :polar, :psf et :hybrid ; dosage de la correction PSF sur le disque
#   v1.1 : suppression du "double limbe" côté lumière cendrée. Le bord réel de la Lune (relief,
#          cratères) s'écarte du cercle ajusté de quelques pixels : la correction du disque
#          s'arrêtait sur le cercle et celle du ciel commençait sur le cercle, si bien que la
#          bande entre le cercle et le vrai bord recevait la mauvaise correction. Le vrai bord est
#          maintenant détecté sur tout le pourtour éclairé par la lumière cendrée (LIMB_EDGE_*),
#          la transition disque / ciel est placée sur ce vrai bord et adoucie sur LIMB_BLEND pixels.
#          Graphique de contrôle <image>_bord_<méthode>.png : écart vrai bord - cercle.
#   v1.2 : traitement de plusieurs fichiers, sélectionnés ensemble dans la fenêtre de dialogue
#          (pick_multi_file) ; graphiques affichés pour le premier fichier seulement (PNG
#          enregistrés pour tous) ; une erreur sur un fichier n'arrête pas les suivants ; bilan
#          final (statut et durée de chaque fichier). Résultats écrits à côté de chaque image.
#
# Le croissant, saturé, diffuse de la lumière dans l'optique : le fond de ciel est surélevé près du
# croissant et décroît en s'éloignant de la Lune jusqu'au niveau "normal" du ciel, mesuré du côté de
# la lumière cendrée (zone de référence). Le script estime cette composante diffusée et la soustrait,
# pour ramener le fond de ciel au niveau de référence. Le croissant n'est jamais modifié.
#
# Trois méthodes au choix (paramètre METHOD) :
#
#   :radial : modèle 1D de la coupe photométrique, Y = a + b ln(d), d = distance au bord de la Lune
#             (pixels), ajusté dans un secteur angulaire centré sur le croissant, puis soustrait dans
#             ce secteur (avec une transition douce sur les bords du secteur). C'est la
#             généralisation directe de la courbe Y = -897 ln(X) + 6672 relevée à la main.
#             Simple, mais suppose que le halo a la même forme dans tout le secteur.
#
#   :polar  : (recommandée) carte empirique du fond de ciel en coordonnées polaires autour du
#             centre de la Lune : médiane du ciel dans des cases (distance au bord x angle), moins
#             le ciel de référence, lissée puis interpolée sur chaque pixel. Suit l'asymétrie
#             réelle du halo (plus fort face au croissant, plus faible vers les cornes) sans
#             supposer de forme de profil. Les étoiles sont éliminées par la médiane et conservées
#             dans l'image corrigée. Ne corrige que le ciel (pas l'intérieur du disque lunaire).
#
#   :psf    : modèle physique : lumière diffusée = (image du croissant) ⊗ (ailes de la PSF),
#             ailes en loi de puissance K(ρ) = (ρ² + ρc²)^(-α/2). L'exposant α et l'amplitude sont
#             ajustés sur les pixels de ciel (calcul sur une image réduite, convolution par FFT).
#             Comme le croissant est saturé, son intensité réelle est inconnue : l'amplitude
#             ajustée absorbe ce facteur. Avantage : le modèle prédit aussi la lumière diffusée
#             SUR le disque (lumière cendrée), ce qui est utile en photométrie de la lumière
#             cendrée (paramètre PSF_CORRECT_DISK). Plus lent, et dépend de la validité du modèle.
#             Sur le disque, la correction est dosée selon la distance au terminateur
#             (PSF_DISK_SCALE, PSF_DISK_W0, PSF_DISK_RAMP) : le croissant étant saturé, sa partie
#             proche du bord, en réalité beaucoup plus brillante, est sous-estimée par rapport à la
#             partie proche du terminateur ; le modèle surestime alors la diffusion près du
#             terminateur (facteur ~2 au contact, sur une image simulée).
#
#   :hybrid : :polar pour le ciel + :psf dosée pour la lumière cendrée (le meilleur des deux).
#             Graphique et tableau supplémentaires : lumière cendrée en fonction de la distance
#             au croissant, avant / après / correction soustraite, pour régler le dosage.
#
# Étapes communes :
#   1) lecture du TIFF (16 bits, RVB ou niveaux de gris) ou du NEF, un ou plusieurs fichiers choisis via IHM
#      (pick_folder / pick_multi_file)
#   2) détection automatique du centre et du rayon de la Lune : ajustement robuste d'un cercle
#      sur le bord extérieur du croissant (paramètres MOON_CENTER / MOON_RADIUS pour forcer)
#   3) niveau du ciel de référence (côté lumière cendrée, bande au bord de l'image, ou REF_BOX)
#   4) modèle de la lumière diffusée (méthode choisie), soustrait canal par canal (R, V, B)
#   5) enregistrement : <image>_corr_<méthode>.tif (même format que l'original),
#      <image>_diffusion_<méthode>.tif (lumière diffusée soustraite, pour contrôle),
#      <image>_diffusion_<méthode>.png (coupe photométrique avant / après et carte du modèle)
#
# Convention de coordonnées : x vers la droite, y vers le BAS, pixel en haut à gauche = (1, 1),
#   comme dans GIMP, Photoshop... (les TIFF sont stockés de haut en bas). Les angles sont comptés
#   dans le sens trigonométrique à l'écran : 0° = vers la droite, 90° = vers le haut.
#
# FICHIERS RAW NIKON (NEF) :
#   Julia ne sait pas décoder les NEF (compression sans perte Nikon, matrice de Bayer) : chaque NEF est
#   converti par le programme dcraw_emu de LibRaw en un TIFF 16 bits temporaire (effacé aussitôt), puis
#   lu normalement.
#   La conversion est LINÉAIRE (option -4 : gamma 1, pas d'éclaircissement automatique), le noir du
#   capteur est soustrait, l'image n'est pas tournée (-t 0 : orientation du capteur, identique pour
#   toutes les images même si l'appareil a enregistré une orientation "portrait").
#   Installation de dcraw_emu :
#     - Windows : télécharger l'archive LibRaw "Win64" sur https://www.libraw.org/download, et copier
#       bin\dcraw_emu.exe ET les DLL du dossier bin (libraw.dll...) dans le dossier de ce script
#       (ou indiquer le chemin complet dans RAW_CONVERTER)
#     - Linux : paquet libraw-bin (Debian / Ubuntu) ; macOS : brew install libraw
#   Version requise : LibRaw 0.21 ou plus récent pour les Z 6II / Z 7II (le D5300 est reconnu par des
#   versions plus anciennes). L'ancien programme dcraw (9.28, 2018) est accepté en secours mais ne
#   reconnaît PAS les Z 6II / Z 7II.
#   Les corrections de noirs (darks) et de plages de lumière uniforme (flats) ne sont pas appliquées.
#   Pour des fichiers RAW d'autres marques (CR2, ARW...), ajouter leur extension à RAW_EXTENSIONS.
#   Résultats d'un NEF : TIFF 16 bits (<image>_corr_<méthode>.tif, <image>_diffusion_<méthode>.tif) ;
#   les valeurs restent LINÉAIRES (proportionnelles à la lumière reçue), ce qui convient au modèle de
#   lumière diffusée. CRESCENT_SIDE se rapporte à l'image dans l'orientation du capteur (après
#   RAW_FLIP_VERTICAL / RAW_FLIP_HORIZONTAL), pas à l'orientation affichée par l'appareil.
#
# Limites :
#   - le niveau "normal" du ciel est supposé uniforme (pas de gradient de pollution lumineuse) ;
#     un gradient de ciel indépendant de la Lune serait en partie absorbé dans la correction
#   - le côté lumière cendrée reçoit lui aussi un peu de lumière diffusée (ailes lointaines) :
#     la correction ramène le ciel au niveau de la zone de référence, pas au "vrai" ciel
#   - les métadonnées du TIFF (profil ICC, EXIF) ne sont pas recopiées dans le fichier corrigé
#   - mémoire : prévoir ~10 fois la taille de l'image en Float32 (≈ 1,5 Go pour 45 Mpixels RVB)
##################################################################################################################

# -----------------------------------------------------------------------------------------------
# PARAMETRES GÉNÉRAUX
# -----------------------------------------------------------------------------------------------
const METHOD = :hybrid            # :radial, :polar, :psf ou :hybrid (voir ci-dessus)
const CRESCENT_SIDE = :right     # côté du croissant dans l'image : :right, :left, :top ou :bottom
                                 # (sert à détecter le bord de la Lune et à placer la zone de référence)
const MOON_CENTER = nothing      # nothing = détection automatique ; sinon (x, y) en pixels, ex. (3012.5, 2210.0)
const MOON_RADIUS = nothing      # nothing = détection automatique ; sinon rayon en pixels, ex. 1650.0
const LIMB_THRESHOLD = 0.5       # seuil de détection du bord du croissant (fraction du maximum de l'image)
const SAT_FRACTION = 0.98        # pixels >= SAT_FRACTION * maximum : croissant saturé, jamais modifiés
const LIMB_MARGIN = 20           # (pixels) ciel ignoré à moins de cette distance du bord (flou du bord,
                                 # erreur sur le rayon) ; la correction y est prolongée à valeur constante
const REF_BOX = nothing          # zone de ciel de référence (x1, x2, y1, y2), ex. (1, 400, 1, 4000) ;
                                 # nothing = bande au bord de l'image du côté opposé au croissant
const REF_STRIP_FRAC = 0.08      # largeur de cette bande automatique (fraction de la largeur de l'image)
const STAT_STEP = 2              # statistiques sur 1 pixel sur STAT_STEP en x et en y (vitesse, mémoire)
const N_DIST_BINS = 60           # nombre d'intervalles de distance au bord (espacés logarithmiquement)
const MIN_PIXELS = 25            # nombre minimal de pixels pour calculer une médiane dans une case
const CLAMP_POSITIVE = true      # true : on ne fait que retirer de la lumière (jamais d'ajout)
const LUMA_WEIGHTS = (0.2126, 0.7152, 0.0722)   # luminance R, V, B (graphiques uniquement)
const PROFILE_HALF_ANGLE = 5.0   # demi-largeur (degrés) du secteur de la coupe photométrique tracée
const SAVE_MODEL = false          # true : enregistre aussi la lumière diffusée soustraite (TIFF)
const SAVE_PNG = false            # true : enregistre les graphiques d'analyse (<image>_diffusion_<méthode>.png,
                                 # <image>_bord_<méthode>.png) ; false : graphiques seulement affichés (voir
                                 # SHOW_PLOTS_FIRST_ONLY), ni calculés ni enregistrés pour les autres images

# -----------------------------------------------------------------------------------------------
# TRAITEMENT DE PLUSIEURS FICHIERS (v1.2)
# Tous les paramètres ci-dessus s'appliquent à chaque image ; la géométrie de la Lune (centre,
# rayon, vrai bord) et le ciel de référence sont recalculés pour chaque image, sauf si
# MOON_CENTER / MOON_RADIUS / REF_BOX sont imposés (ils valent alors pour toutes les images).
# -----------------------------------------------------------------------------------------------
const SHOW_PLOTS_FIRST_ONLY = true   # true : graphiques affichés pour le premier fichier seulement
                                     # (les PNG sont enregistrés pour tous les fichiers)
const STOP_ON_ERROR = false      # false : une erreur sur un fichier n'arrête pas le traitement des suivants
const SKIP_OUTPUT_FILES = true   # true : ignore les fichiers produits par ce script (*_corr_*.tif,
                                 # *_diffusion_*.tif) s'ils ont été sélectionnés par erreur

# -----------------------------------------------------------------------------------------------
# FICHIERS RAW NIKON (NEF), convertis par dcraw_emu (LibRaw) - v1.3
# -----------------------------------------------------------------------------------------------
const RAW_EXTENSIONS = ["nef", "nrw"]   # extensions (en minuscules) traitées comme fichiers RAW
const RAW_CONVERTER = "C:\\Users\\viladrich\\Documents\\Astro\\documentation\\logiciels\\Julia\\autres\\bin\\dcraw_emu.exe"         # chemin complet de dcraw_emu (ex. raw"C:\LibRaw\bin\dcraw_emu.exe") ;
                                 # "" : recherche dans le dossier du script, puis dans le PATH
                                 # (dcraw_emu en priorité, dcraw en secours)
const RAW_WB = :camera           # balance des blancs : :camera   = celle de l'appareil (option -w)
                                 #                      :daylight = "lumière du jour" de l'appareil (défaut de LibRaw)
                                 #                      :none     = aucune (-r 1 1 1 1 : signal brut de chaque couche)
const RAW_COLOR = 0              # espace couleur de sortie : 0 = couleurs brutes du capteur, 1 = sRGB
                                 # (primaires sRGB, mais toujours LINÉAIRE), 4 = ProPhoto, 5 = XYZ...
const RAW_DEMOSAIC = 3           # dématriçage : 0 = bilinéaire, 1 = VNG, 2 = PPG, 3 = AHD
const RAW_HALF_SIZE = false      # true : demi-résolution (option -h : chaque bloc 2x2 de la matrice de Bayer
                                 # donne un pixel RVB, sans interpolation)
const RAW_EXTRA_ARGS = String[]  # options supplémentaires de dcraw_emu, ex. ["-H", "1"]
const RAW_FLIP_VERTICAL = false   # true : retournement haut / bas de l'image convertie (image décodée à l'envers)
const RAW_FLIP_HORIZONTAL = false  # true : retournement gauche / droite de l'image convertie

# -----------------------------------------------------------------------------------------------
# VRAI BORD DE LA LUNE (v1.1) : transition entre correction du disque et correction du ciel
# -----------------------------------------------------------------------------------------------
const LIMB_EDGE_DETECT = true    # true : détection du vrai bord (relief) ; false : cercle ajusté
const LIMB_EDGE_SEARCH = 25      # (pixels) recherche du vrai bord à ± cette distance du cercle
const LIMB_EDGE_ARC = 3.0        # (pixels) résolution le long du bord (largeur d'un intervalle d'angle)
const LIMB_EDGE_LEVEL = 0.5      # bord = niveau ciel + LIMB_EDGE_LEVEL x (lumière cendrée - ciel)
const LIMB_EDGE_MIN_SNR = 3.0    # contraste lumière cendrée / ciel minimal (en sigma du bruit du ciel)
const LIMB_EDGE_SMOOTH = 2       # lissage du bord détecté : médiane glissante sur ±N intervalles
                                 # (2 => 15 pixels avec LIMB_EDGE_ARC = 3 : élimine étoiles et bruit,
                                 # conserve les cratères ; 0 = aucun lissage)
const LIMB_BLEND = 3.0           # (pixels) largeur du fondu disque / ciel centré sur le bord
                                 # (0 = transition franche). Avec LIMB_EDGE_DETECT = false, fondu
                                 # centré sur le cercle ajusté.
const PLOT_WINDOW = :display     # :display = volet Plots de VS Code ; :gui = fenêtre séparée

# -----------------------------------------------------------------------------------------------
# PARAMETRES DE LA MÉTHODE :radial
# -----------------------------------------------------------------------------------------------
const RADIAL_HALF_ANGLE = 45.0   # demi-ouverture (degrés) du secteur centré sur le croissant
const RADIAL_TAPER = 45.0        # transition (degrés) au-delà du secteur : correction -> 0 (cos²)
const RADIAL_FIT_DMAX = nothing  # distance maximale au bord (pixels) pour l'ajustement ; nothing =
                                 # automatique (tant que l'excès dépasse 2 % de l'excès près du bord)

# -----------------------------------------------------------------------------------------------
# PARAMETRES DE LA MÉTHODE :polar
# -----------------------------------------------------------------------------------------------
const POLAR_N_SECTORS = 120       # nombre de secteurs angulaires (72 => 5°)
const POLAR_SMOOTH_THETA = 1     # lissage (médiane glissante) sur ±N secteurs
const POLAR_SMOOTH_R = 1         # lissage (médiane glissante) sur ±N intervalles de distance

# -----------------------------------------------------------------------------------------------
# PARAMETRES DE LA MÉTHODE :psf
# -----------------------------------------------------------------------------------------------
const PSF_DOWNSAMPLE = 8         # facteur de réduction de l'image pour l'ajustement (8 => blocs 8x8)
const PSF_CORE = 3.0             # (pixels) rayon de cœur du noyau, évite la singularité en ρ = 0
const PSF_ALPHA_RANGE = (1.0, 4.0)   # plage de recherche de l'exposant α (pas 0,1 puis affinage 0,02)
const PSF_SOURCE_FRACTION = 0.10 # source = pixels du disque au-dessus de cette fraction du maximum
const PSF_CLIP = 3.0             # ajustement robuste : rejet des blocs à plus de N sigma (étoiles...)
const PSF_CORRECT_DISK = true   # true : corrige aussi la lumière cendrée (disque hors croissant)
const PSF_ZERO_LEVEL = :reference # :reference = ciel ramené au niveau de la zone de référence (comme
                                 # :radial et :polar) ; :model = retire TOUTE la diffusion modélisée,
                                 # y compris celle qui atteint le côté lumière cendrée (ciel ramené au
                                 # fond ajusté B, légèrement plus bas que la référence)

# dosage de la correction PSF SUR LE DISQUE (lumière cendrée), en fonction de la distance t au
# croissant (au terminateur) : poids(t) = PSF_DISK_SCALE * (W0 + (1 - W0) * sin²(π/2 · t / RAMP)),
# poids = PSF_DISK_SCALE au-delà de RAMP. Avec W0 = 1 (ou RAMP = 0) : dosage uniforme.
# Exemple : sur-correction sur ~800 pixels => RAMP = 800, puis ajuster W0 (0,3 à 0,7) pour que la
# bande près du terminateur ne soit ni trop sombre ni trop claire (voir la coupe "lumière cendrée").
const PSF_DISK_SCALE = 0.5       # facteur global sur la correction du disque (1 = modèle brut, < 1 sous-corrigé)
const PSF_DISK_W0 = 0.5          # poids au contact du terminateur (0 = aucune correction, 1 = complète)
const PSF_DISK_RAMP = 800.0      # (pixels) largeur de la zone de transition depuis le terminateur

# -----------------------------------------------------------------------------------------------
# MÉTHODE :hybrid = :polar pour le ciel + :psf (dosé ci-dessus) pour la lumière cendrée
# -----------------------------------------------------------------------------------------------
const DISK_PROFILE_STEP = 50.0   # (pixels) pas de la coupe de contrôle sur la lumière cendrée

##################################################################################################################
# LECTURE / ÉCRITURE DES TIFF
##################################################################################################################

# valeur d'un pixel en ADU (0..65535 en 16 bits, 0..255 en 8 bits)
to_adu(v::FixedPointNumbers.Normed) = Float32(reinterpret(v))
to_adu(v::AbstractFloat) = Float32(v) * 65535f0     # TIFF flottant 0..1 => échelle 16 bits
to_adu(v::Integer) = Float32(v)

#############################################################################################
# matrice de couleurs (TiffImages) -> vecteur de plans Float32 [y, x], y = 1 en haut
# (fonction séparée : le type des pixels y est connu, la boucle est compilée efficacement)
#############################################################################################
function extract_channels(raw::AbstractMatrix)
    ny, nx = size(raw)
    c0 = ColorTypes.color(raw[1, 1])                # couleur sans le canal alpha éventuel
    if c0 isa ColorTypes.AbstractGray
        g = Matrix{Float32}(undef, ny, nx)
        @inbounds for j in 1:nx, i in 1:ny
            g[i, j] = to_adu(ColorTypes.gray(ColorTypes.color(raw[i, j])))
        end
        return [g]
    elseif c0 isa ColorTypes.AbstractRGB
        r = Matrix{Float32}(undef, ny, nx)
        v = similar(r)
        b = similar(r)
        @inbounds for j in 1:nx, i in 1:ny
            c = ColorTypes.color(raw[i, j])
            r[i, j] = to_adu(ColorTypes.red(c))
            v[i, j] = to_adu(ColorTypes.green(c))
            b[i, j] = to_adu(ColorTypes.blue(c))
        end
        return [r, v, b]
    else
        error("type de pixel TIFF non géré : ", typeof(raw[1, 1]))
    end
end

#############################################################################################
# lecture d'un TIFF. Sortie : (plans Float32 [y, x], nombre de bits)
#############################################################################################
function read_tiff(path::AbstractString)
    raw = try
        TiffImages.load(path)
    catch err
        msg = sprint(showerror, err)
        println("TiffImages n'a pas pu lire le fichier : ", first(msg, 150))
        println("=> lecture avec le lecteur TIFF de secours (TIFF non compressé)")
        return read_tiff_basic(path)
    end
    if ndims(raw) == 3
        println("TIFF à plusieurs pages : seule la première est utilisée")
        raw = raw[:, :, 1]
    end
    ndims(raw) == 2 || error("format TIFF non géré (", ndims(raw), " dimensions)")
    T = eltype(ColorTypes.color(raw[1, 1]))
    bits = T <: FixedPointNumbers.Normed ? 8 * sizeof(T) : 16
    chans = extract_channels(raw)
    println("TIFF ", length(chans) == 3 ? "couleur RVB" : "niveaux de gris", ", ", bits, " bits",
            T <: AbstractFloat ? " (flottant, mis à l'échelle 0..65535)" : "")
    bits < 16 && println("ATTENTION : image ", bits, " bits, dynamique faible")
    return chans, bits
end

#############################################################################################
# lecteur TIFF minimal (secours) : TIFF NON COMPRESSÉ, 8 ou 16 bits entiers ou 32 bits flottants,
# niveaux de gris ou RVB (canal alpha ignoré), entrelacé ou planaire, en bandes (strips).
# Tolère l'absence du mot-clé StripByteCounts (tailles recalculées), absence qui fait échouer
# TiffImages ("BoundsError: attempt to access 0-element Vector{TiffImages.Tag}").
# Sortie : (plans Float32 [y, x], y = 1 en haut ; nombre de bits)
#############################################################################################
function read_tiff_basic(path::AbstractString)
    data = read(path)
    length(data) >= 8 || error("fichier TIFF trop court")
    if data[1] == 0x49 && data[2] == 0x49          # "II" : petit-boutiste (Intel)
        le = true
    elseif data[1] == 0x4d && data[2] == 0x4d      # "MM" : gros-boutiste (Motorola)
        le = false
    else
        error("ce fichier n'est pas un TIFF")
    end
    # lecture d'entiers à l'octet p (compté à partir de 0, comme dans la norme TIFF)
    rd16(p) = le ? (UInt16(data[p+1]) | (UInt16(data[p+2]) << 8)) :
                   ((UInt16(data[p+1]) << 8) | UInt16(data[p+2]))
    rd32(p) = le ? (UInt32(data[p+1]) | (UInt32(data[p+2]) << 8) | (UInt32(data[p+3]) << 16) | (UInt32(data[p+4]) << 24)) :
                   ((UInt32(data[p+1]) << 24) | (UInt32(data[p+2]) << 16) | (UInt32(data[p+3]) << 8) | UInt32(data[p+4]))
    rd16(2) == 42 || error("TIFF non standard (BigTIFF non géré)")

    # première IFD (répertoire des mots-clés) : seuls les types SHORT (3) et LONG (4) sont utiles ici
    ifd = Int(rd32(4))
    ntags = Int(rd16(ifd))
    tags = Dict{Int, Vector{Int}}()
    for e in 0:ntags-1
        p = ifd + 2 + 12 * e
        tag = Int(rd16(p))
        typ = Int(rd16(p + 2))
        cnt = Int(rd32(p + 4))
        typ in (3, 4) || continue
        sz = typ == 3 ? 2 : 4
        q = cnt * sz <= 4 ? p + 8 : Int(rd32(p + 8))     # valeur dans l'entrée ou à une adresse
        tags[tag] = [Int(typ == 3 ? rd16(q + sz * i) : rd32(q + sz * i)) for i in 0:cnt-1]
    end
    get1(t, default) = haskey(tags, t) ? tags[t][1] : default

    nx = get1(256, 0)              # ImageWidth
    ny = get1(257, 0)              # ImageLength
    bps = get1(258, 1)             # BitsPerSample
    comp = get1(259, 1)            # Compression
    spp = get1(277, 1)             # SamplesPerPixel
    rps = min(get1(278, ny), ny)   # RowsPerStrip
    planar = get1(284, 1)          # PlanarConfiguration : 1 = entrelacé (RVBRVB...), 2 = plans séparés
    sfmt = get1(339, 1)            # SampleFormat : 1 = entier non signé, 3 = flottant
    @printf("TIFF : %d x %d, %d canal(aux), %d bits, compression %d, %s, %d bande(s)%s\n",
            nx, ny, spp, bps, comp, planar == 1 ? "entrelacé" : "planaire",
            haskey(tags, 273) ? length(tags[273]) : 0,
            haskey(tags, 279) ? "" : " (mot-clé StripByteCounts absent)")
    (nx > 0 && ny > 0) || error("dimensions de l'image absentes")
    haskey(tags, 273) || error("TIFF en tuiles ou sans StripOffsets : non géré par le lecteur de secours")
    comp == 1 || error("TIFF compressé (code ", comp, ") : non géré par le lecteur de secours. ",
                       "Réenregistrer l'image en TIFF 16 bits non compressé.")
    ((bps == 8 || bps == 16) && sfmt == 1) || (bps == 32 && sfmt == 3) ||
        error("format de pixel non géré : ", bps, " bits, SampleFormat ", sfmt)
    offs = tags[273]
    bsz = bps ÷ 8
    nc = spp >= 3 ? 3 : 1          # canaux utiles (alpha ignoré)

    sample(p) = bps == 8 ? Float32(data[p+1]) :
                bps == 16 ? Float32(rd16(p)) : reinterpret(Float32, rd32(p)) * 65535f0

    chans = [Matrix{Float32}(undef, ny, nx) for _ in 1:nc]
    if planar == 1
        rowbytes = nx * spp * bsz
        for y in 1:ny
            k = (y - 1) ÷ rps                                  # bande contenant la ligne (0..)
            base = offs[k+1] + ((y - 1) - k * rps) * rowbytes
            base + rowbytes <= length(data) || error("fichier TIFF tronqué")
            for x in 1:nx, c in 1:nc
                chans[c][y, x] = sample(base + ((x - 1) * spp + (c - 1)) * bsz)
            end
        end
    else
        rowbytes = nx * bsz
        sp = cld(ny, rps)                                      # bandes par plan
        length(offs) >= nc * sp || error("TIFF planaire : nombre de bandes incohérent")
        for c in 1:nc, y in 1:ny
            k = (y - 1) ÷ rps
            base = offs[(c-1)*sp + k + 1] + ((y - 1) - k * rps) * rowbytes
            base + rowbytes <= length(data) || error("fichier TIFF tronqué")
            for x in 1:nx
                chans[c][y, x] = sample(base + (x - 1) * bsz)
            end
        end
    end
    println("TIFF ", nc == 3 ? "couleur RVB" : "niveaux de gris", " lu par le lecteur de secours")
    return chans, (bps == 32 ? 16 : bps)
end

#############################################################################################
# écriture d'un TIFF (8 ou 16 bits entiers, RVB ou niveaux de gris) à partir de plans Float32
#############################################################################################
function save_tiff(path::AbstractString, chans::Vector{Matrix{Float32}}, bits::Int)
    if bits <= 8
        save_tiff_typed(path, chans, FixedPointNumbers.N0f8, UInt8)
    else
        save_tiff_typed(path, chans, FixedPointNumbers.N0f16, UInt16)
    end
end

function save_tiff_typed(path, chans, ::Type{T}, ::Type{U}) where {T, U}
    ny, nx = size(chans[1])
    maxcode = Int(typemax(U))
    q(v) = reinterpret(T, U(clamp(round(Int, v), 0, maxcode)))
    if length(chans) == 3
        out = Matrix{ColorTypes.RGB{T}}(undef, ny, nx)
        @inbounds for j in 1:nx, i in 1:ny
            out[i, j] = ColorTypes.RGB{T}(q(chans[1][i, j]), q(chans[2][i, j]), q(chans[3][i, j]))
        end
    else
        out = Matrix{ColorTypes.Gray{T}}(undef, ny, nx)
        @inbounds for j in 1:nx, i in 1:ny
            out[i, j] = ColorTypes.Gray{T}(q(chans[1][i, j]))
        end
    end
    TiffImages.save(path, out)
    println("enregistré : ", path)
    return nothing
end

##################################################################################################################
# LECTURE DES FICHIERS RAW NIKON (NEF) : conversion par dcraw_emu (LibRaw) en TIFF 16 bits linéaire
##################################################################################################################

is_raw_file(path::AbstractString) = lowercase(lstrip(splitext(path)[2], '.')) in RAW_EXTENSIONS

const RAW_CONV_PATH = Ref{String}("")                 # convertisseur trouvé (cherché une seule fois)

#############################################################################################
# recherche du convertisseur : RAW_CONVERTER, puis dossier du script, puis PATH
#############################################################################################
function raw_converter()
    isempty(RAW_CONV_PATH[]) || return RAW_CONV_PATH[]
    found = ""
    if !isempty(RAW_CONVERTER)
        isfile(RAW_CONVERTER) || error("RAW_CONVERTER : programme introuvable : ", RAW_CONVERTER)
        found = RAW_CONVERTER
    else
        for name in ("dcraw_emu.exe", "dcraw_emu", "dcraw.exe", "dcraw")
            p = joinpath(@__DIR__, name)
            if isfile(p)
                found = p
                break
            end
        end
        if isempty(found)
            for name in ("dcraw_emu", "dcraw")
                p = Sys.which(name)
                if p !== nothing
                    found = p
                    break
                end
            end
        end
    end
    isempty(found) && error("lecture des fichiers RAW : programme dcraw_emu (LibRaw) introuvable.\n",
        "  Windows : copier dcraw_emu.exe et les DLL du dossier bin de l'archive LibRaw Win64\n",
        "            (https://www.libraw.org/download) dans le dossier du script : ", @__DIR__, "\n",
        "  Linux : sudo apt install libraw-bin ; macOS : brew install libraw\n",
        "  ou indiquer le chemin complet dans le paramètre RAW_CONVERTER")
    occursin("dcraw_emu", lowercase(basename(found))) ||
        println("ATTENTION : convertisseur dcraw (ancien) : les Nikon Z 6II / Z 7II ne sont pas reconnus, ",
                "installer dcraw_emu (LibRaw 0.21 ou plus récent)")
    RAW_CONV_PATH[] = found
    return found
end

# options de conversion : TIFF (-T), 16 bits linéaire sans éclaircissement automatique (-4),
# pas de rotation (-t 0), dématriçage, espace couleur, balance des blancs
function raw_args()
    a = ["-T", "-4", "-t", "0", "-q", string(RAW_DEMOSAIC), "-o", string(RAW_COLOR)]
    if RAW_WB == :camera
        push!(a, "-w")
    elseif RAW_WB == :none
        append!(a, ["-r", "1", "1", "1", "1"])
    elseif RAW_WB != :daylight
        error("RAW_WB inconnu : ", RAW_WB, " (attendu :camera, :daylight ou :none)")
    end
    RAW_HALF_SIZE && push!(a, "-h")
    append!(a, RAW_EXTRA_ARGS)
    return a
end

# retournements de l'image RAW convertie (RAW_FLIP_VERTICAL, RAW_FLIP_HORIZONTAL)
function raw_orient(chans::Vector{Matrix{Float32}})
    if RAW_FLIP_VERTICAL && RAW_FLIP_HORIZONTAL
        return [c[end:-1:1, end:-1:1] for c in chans]
    elseif RAW_FLIP_VERTICAL
        return [c[end:-1:1, :] for c in chans]
    elseif RAW_FLIP_HORIZONTAL
        return [c[:, end:-1:1] for c in chans]
    end
    return chans
end

#############################################################################################
# conversion d'un fichier RAW dans un dossier temporaire (le fichier d'origine n'est pas modifié ;
# dcraw_emu écrit <nom>.NEF.tiff, dcraw écrit <nom>.tiff). Sortie : (plans Float32 [y, x], 16)
#############################################################################################
function read_raw(path::AbstractString)
    conv = raw_converter()
    tmp = mktempdir(; prefix="raw_dlc_")
    try
        src = joinpath(tmp, basename(path))
        cp(path, src; force=true)
        errf = joinpath(tmp, "conversion_stderr.txt")
        t0 = time()
        cmd = ignorestatus(Cmd(`$conv $(raw_args()) $src`; dir=tmp))
        proc = run(pipeline(cmd; stdout=devnull, stderr=errf))
        rm(src; force=true)
        tifs = filter(f -> occursin(r"\.tiff?$"i, f), readdir(tmp))
        if proc.exitcode != 0 || isempty(tifs)
            msg = isfile(errf) ? strip(read(errf, String)) : ""
            error("échec de la conversion RAW par ", basename(conv), " (code ", proc.exitcode, ")",
                  isempty(msg) ? "" : " : " * first(replace(msg, '\n' => ' '), 200),
                  " - appareil non reconnu ? (LibRaw 0.21 ou plus récent requis pour les Z 6II / Z 7II)")
        end
        @printf("RAW converti par %s en %.1f s (linéaire 16 bits, balance des blancs %s, espace couleur %d)\n",
                basename(conv), time() - t0, RAW_WB, RAW_COLOR)
        chans, _ = read_tiff(joinpath(tmp, first(tifs)))
        return raw_orient(chans), 16
    finally
        try
            rm(tmp; recursive=true, force=true)     # Windows : peut échouer si un fichier est encore ouvert
        catch
        end
    end
end

# lecture d'une image : RAW (NEF...) ou TIFF selon l'extension
read_image(path::AbstractString) = is_raw_file(path) ? read_raw(path) : read_tiff(path)

##################################################################################################################
# OUTILS STATISTIQUES
##################################################################################################################

#############################################################################################
# fond de ciel et bruit robustes : médiane et MAD (convertie en sigma), sigma-clipping à 3 sigma
# (repris de HFD_star) - Sortie : (fond, sigma)
#############################################################################################
function robust_background(values::AbstractVector{Float64})
    w = filter(isfinite, collect(values))
    isempty(w) && error("aucun pixel valide pour estimer le fond de ciel")
    m = median(w)
    s = 1.4826 * median(abs.(w .- m))
    for _ in 1:5
        s > 0 || break
        keep = w[abs.(w .- m) .< 3.0 * s]
        (length(keep) < 10 || length(keep) == length(w)) && break
        w = keep
        m = median(w)
        s = 1.4826 * median(abs.(w .- m))
    end
    if s <= 0
        s = length(w) > 1 ? std(w) : 0.0
    end
    return m, max(s, 1e-6)
end

# luminance d'un pixel (graphiques)
@inline function lumval(chans::Vector{Matrix{Float32}}, y::Int, x::Int)
    length(chans) == 3 || return Float64(chans[1][y, x])
    wr, wv, wb = LUMA_WEIGHTS
    return wr * chans[1][y, x] + wv * chans[2][y, x] + wb * chans[3][y, x]
end

##################################################################################################################
# GÉOMÉTRIE DE LA LUNE
##################################################################################################################

struct MoonGeom
    cx::Float64      # centre (pixels, origine en haut à gauche)
    cy::Float64
    R::Float64       # rayon (pixels)
    θc::Float64      # direction du croissant vue du centre (degrés, 0 = droite, 90 = haut)
end

# distance au BORD de la Lune (négative à l'intérieur du disque) et angle (degrés, -180..180)
@inline function dist_angle(g::MoonGeom, x::Real, y::Real)
    dx = x - g.cx
    dy = g.cy - y                      # y vers le bas dans l'image => signe inversé
    return sqrt(dx * dx + dy * dy) - g.R, atand(dy, dx)
end

@inline wrap180(a) = mod(a + 180.0, 360.0) - 180.0

#############################################################################################
# position (sous-pixel) où une ligne de pixels franchit le seuil, en partant d'une extrémité.
# Exige 3 pixels consécutifs au-dessus du seuil (élimine pixels chauds et petites étoiles).
#############################################################################################
function edge_crossing(line::AbstractVector, thr, from_end::Bool)
    n = length(line)
    if from_end
        k = findlast(>=(thr), line)
        (k === nothing || k >= n || k < 3) && return nothing
        (line[k-1] >= thr && line[k-2] >= thr) || return nothing
        return k + (line[k] - thr) / (line[k] - line[k+1])
    else
        k = findfirst(>=(thr), line)
        (k === nothing || k <= 1 || k > n - 2) && return nothing
        (line[k+1] >= thr && line[k+2] >= thr) || return nothing
        return k - (line[k] - thr) / (line[k] - line[k-1])
    end
end

# points du bord extérieur (éclairé) du croissant : pixel brillant le plus extrême de chaque ligne
function limb_points(lum::Matrix{Float32}, thr, side::Symbol)
    ny, nx = size(lum)
    xs = Float64[]
    ys = Float64[]
    if side in (:right, :left)
        for y in 1:ny
            p = edge_crossing(view(lum, y, :), thr, side == :right)
            p === nothing || (push!(xs, p); push!(ys, y))
        end
    elseif side in (:top, :bottom)
        for x in 1:nx
            p = edge_crossing(view(lum, :, x), thr, side == :bottom)
            p === nothing || (push!(xs, x); push!(ys, p))
        end
    else
        error("CRESCENT_SIDE doit valoir :right, :left, :top ou :bottom")
    end
    return xs, ys
end

# ajustement algébrique d'un cercle (méthode de Kåsa) : x² + y² + D x + E y + F = 0
function fit_circle(xs::Vector{Float64}, ys::Vector{Float64})
    mx = mean(xs)
    my = mean(ys)
    u = xs .- mx
    v = ys .- my
    M = hcat(u, v, ones(length(u)))
    D, E, F = M \ (-(u .^ 2 .+ v .^ 2))
    cu = -D / 2
    cv = -E / 2
    return cu + mx, cv + my, sqrt(cu^2 + cv^2 - F)
end

#############################################################################################
# centre et rayon de la Lune : cercle ajusté sur le bord du croissant, avec rejet itératif des
# points aberrants (étoiles brillantes, irrégularités du bord)
#############################################################################################
function detect_moon(lum::Matrix{Float32}, thr, side::Symbol)
    xs, ys = limb_points(lum, thr, side)
    length(xs) >= 20 || error("bord du croissant non trouvé (", length(xs), " points) : vérifier ",
                              "CRESCENT_SIDE et LIMB_THRESHOLD, ou donner MOON_CENTER et MOON_RADIUS")
    keep = trues(length(xs))
    cx, cy, R = 0.0, 0.0, 0.0
    for _ in 1:10
        cx, cy, R = fit_circle(xs[keep], ys[keep])
        res = hypot.(xs .- cx, ys .- cy) .- R
        σ = 1.4826 * median(abs.(res[keep]))
        newkeep = abs.(res) .< max(3σ, 1.0)
        (newkeep == keep || count(newkeep) < 20) && break
        keep = newkeep
    end
    rms = sqrt(mean(abs2, hypot.(xs[keep] .- cx, ys[keep] .- cy) .- R))
    @printf("bord du croissant : %d points (%d retenus), écart quadratique au cercle %.2f pixels\n",
            length(xs), count(keep), rms)
    return cx, cy, R
end

# direction du croissant vue du centre : moyenne des directions des pixels saturés
function crescent_direction(lum::Matrix{Float32}, sat_thr, cx, cy)
    sx = 0.0
    sy = 0.0
    ny, nx = size(lum)
    for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
        lum[y, x] >= sat_thr || continue
        dx = x - cx
        dy = cy - y
        r = hypot(dx, dy)
        r > 0 || continue
        sx += dx / r
        sy += dy / r
    end
    if sx == 0 && sy == 0
        return Dict(:right => 0.0, :top => 90.0, :left => 180.0, :bottom => -90.0)[CRESCENT_SIDE]
    end
    return atand(sy, sx)
end

#############################################################################################
# intervalles de distance au bord, espacés logarithmiquement de LIMB_MARGIN au coin le plus
# éloigné (fins près du bord où le halo varie vite, larges loin de la Lune)
#############################################################################################
function dist_edges(nx::Int, ny::Int, g::MoonGeom)
    dmin = Float64(max(LIMB_MARGIN, 1))
    dmax = maximum(hypot(x - g.cx, y - g.cy) for (x, y) in ((1, 1), (nx, 1), (1, ny), (nx, ny))) - g.R
    dmax > 1.5 * dmin || error("la Lune remplit l'image : pas de ciel à analyser")
    n = N_DIST_BINS
    return [dmin * (dmax / dmin)^(k / n) for k in 0:n]
end

# numéro de l'intervalle contenant d (0 si hors limites)
@inline function dist_bin(d::Float64, edges::Vector{Float64})
    d < edges[1] && return 0
    n = length(edges) - 1
    k = floor(Int, log(d / edges[1]) / log(edges[end] / edges[1]) * n) + 1
    return k <= n ? k : 0
end

# centres (moyenne géométrique) des intervalles
bin_centers(edges) = [sqrt(edges[k] * edges[k+1]) for k in 1:length(edges)-1]

#############################################################################################
# zone de ciel de référence : fonction (x, y) -> Bool
#############################################################################################
function make_ref_test(nx::Int, ny::Int, g::MoonGeom)
    outside(x, y) = dist_angle(g, x, y)[1] >= LIMB_MARGIN
    if REF_BOX !== nothing
        x1, x2, y1, y2 = REF_BOX
        return (x, y) -> (x1 <= x <= x2) && (y1 <= y <= y2) && outside(x, y)
    end
    w = REF_STRIP_FRAC * nx
    h = REF_STRIP_FRAC * ny
    if CRESCENT_SIDE == :right
        return (x, y) -> x <= w && outside(x, y)
    elseif CRESCENT_SIDE == :left
        return (x, y) -> x > nx - w && outside(x, y)
    elseif CRESCENT_SIDE == :top
        return (x, y) -> y > ny - h && outside(x, y)
    else
        return (x, y) -> y <= h && outside(x, y)
    end
end

##################################################################################################################
# VRAI BORD DE LA LUNE (v1.1)
##################################################################################################################

# écart (pixels) entre le vrai bord et le cercle ajusté, par intervalle d'angle (vide = cercle)
const LIMB_EDGE = Ref(Float64[])

# écart vrai bord - cercle à l'angle θ (degrés), interpolé entre centres des intervalles
@inline function edge_offset(θ::Float64)
    off = LIMB_EDGE[]
    n = length(off)
    n == 0 && return 0.0
    t = (θ + 180.0) / (360.0 / n) + 0.5
    i0 = floor(Int, t)
    w = t - i0
    return (1 - w) * off[mod1(i0, n)] + w * off[mod1(i0 + 1, n)]
end

#############################################################################################
# fraction "Lune" d'un pixel (1 = disque, 0 = ciel), d = distance au cercle (négative dedans) :
# fondu linéaire de largeur LIMB_BLEND centré sur le vrai bord. Sert à répartir la correction
# entre celle du disque (lumière cendrée) et celle du ciel.
#############################################################################################
@inline function moon_weight(d::Float64, θ::Float64)
    band = LIMB_EDGE_SEARCH + LIMB_BLEND
    d < -band && return 1.0
    d > band && return 0.0
    dE = edge_offset(θ)
    LIMB_BLEND <= 0 && return d < dE ? 1.0 : 0.0
    return clamp(0.5 - (d - dE) / LIMB_BLEND, 0.0, 1.0)
end

#############################################################################################
# détection du vrai bord : pour chaque intervalle d'angle (LIMB_EDGE_ARC pixels le long du bord),
# profil radial (pas de 1 pixel, image lissée 3x3) à ± LIMB_EDGE_SEARCH du cercle ; le bord est
# le premier passage, de l'intérieur vers l'extérieur, sous le niveau
#     ciel + LIMB_EDGE_LEVEL x (lumière cendrée - ciel)
# (niveaux = médianes des profils sur la moitié intérieure et la moitié extérieure de la zone).
# Côté croissant (intérieur brillant) ou contraste trop faible : pas de détection, écart = 0
# (cercle), sauf trous courts (<= 5 intervalles) interpolés. Puis médiane glissante.
#############################################################################################
function detect_true_limb(chans, g::MoonGeom, src_thr, σlum)
    ny, nx = size(chans[1])
    L = LIMB_EDGE_SEARCH
    nθ = max(360, round(Int, 2π * g.R / LIMB_EDGE_ARC))
    wθ = 360.0 / nθ
    nd = 2L + 1
    sums = zeros(nd, nθ)
    cnt = zeros(Int, nd, nθ)
    x1 = max(2, floor(Int, g.cx - g.R - L - 1))
    x2 = min(nx - 1, ceil(Int, g.cx + g.R + L + 1))
    y1 = max(2, floor(Int, g.cy - g.R - L - 1))
    y2 = min(ny - 1, ceil(Int, g.cy + g.R + L + 1))
    @inbounds for x in x1:x2, y in y1:y2
        d, θ = dist_angle(g, x, y)
        abs(d) <= L + 0.49 || continue
        s = 0.0
        for j in -1:1, i in -1:1
            s += lumval(chans, y + i, x + j)
        end
        id = clamp(round(Int, d) + L + 1, 1, nd)          # d = id - L - 1
        it = clamp(floor(Int, (θ + 180.0) / wθ) + 1, 1, nθ)
        sums[id, it] += s / 9
        cnt[id, it] += 1
    end

    h = max(3, L ÷ 2)
    off = fill(NaN, nθ)
    ncres = 0
    nlow = 0
    for it in 1:nθ
        p = [cnt[k, it] > 0 ? sums[k, it] / cnt[k, it] : NaN for k in 1:nd]
        inner = filter(isfinite, p[1:h])
        outer = filter(isfinite, p[nd-h+1:nd])
        (length(inner) >= 3 && length(outer) >= 3) || continue
        Iin = median(inner)
        Iout = median(outer)
        if Iin >= src_thr
            ncres += 1                          # croissant : on garde le cercle
            continue
        end
        C = Iin - Iout
        if C < LIMB_EDGE_MIN_SNR * σlum
            nlow += 1
            continue
        end
        T = Iout + LIMB_EDGE_LEVEL * C
        for k in 2:nd                           # de l'intérieur vers l'extérieur
            if isfinite(p[k-1]) && isfinite(p[k]) && p[k-1] >= T > p[k]
                off[it] = (k - 2 - L) + (p[k-1] - T) / (p[k-1] - p[k])
                break
            end
        end
    end
    ndet = count(isfinite, off)

    # trous courts interpolés, longs => cercle (écart 0)
    gap = 5
    filled = similar(off)
    for it in 1:nθ
        if isfinite(off[it])
            filled[it] = off[it]
            continue
        end
        a = 0
        b = 0
        for k in 1:gap
            (a == 0 && isfinite(off[mod1(it - k, nθ)])) && (a = k)
            (b == 0 && isfinite(off[mod1(it + k, nθ)])) && (b = k)
        end
        filled[it] = (a > 0 && b > 0) ?
                     (b * off[mod1(it - a, nθ)] + a * off[mod1(it + b, nθ)]) / (a + b) : 0.0
    end
    hs = LIMB_EDGE_SMOOTH
    sm = [median([filled[mod1(it + j, nθ)] for j in -hs:hs]) for it in 1:nθ]
    sm .= clamp.(sm, -(L - 2), L - 2)

    det = filter(isfinite, off)
    @printf("vrai bord : détecté sur %.0f %% du pourtour (%d intervalles de %.1f px ; croissant %.0f %%, contraste insuffisant %.0f %%)\n",
            100 * ndet / nθ, nθ, LIMB_EDGE_ARC, 100 * ncres / nθ, 100 * nlow / nθ)
    if !isempty(det)
        @printf("   écart vrai bord - cercle : médiane %+.1f px, de %+.1f à %+.1f px (centiles 1 %% et 99 %%)\n",
                median(det), quantile(det, 0.01), quantile(det, 0.99))
    end
    return sm, off
end

function plot_limb_edge(path, sm, raw, g::MoonGeom, title_str)
    n = length(sm)
    θs = [-180.0 + (i - 0.5) * 360.0 / n for i in 1:n]
    ok = findall(isfinite, raw)
    p = scatter(θs[ok], raw[ok]; ms=1.5, msw=0, color=:gray, label="détection brute",
                xlabel="angle (°, 0 = droite, 90 = haut)", ylabel="vrai bord - cercle (pixels)",
                title=title_str, xlims=(-180, 180), size=(1600, 550), legend=:topleft,
                left_margin=5Plots.mm, bottom_margin=5Plots.mm)
    plot!(p, θs, sm; color=:red, lw=1.5, label="bord retenu (lissé)")
    hline!(p, [0.0]; color=:black, lw=0.5, label="cercle ajusté")
    vline!(p, [g.θc]; color=:orange, ls=:dash, label="direction du croissant")
    save_plot(p, path)
    show_plot(p)
    return nothing
end

##################################################################################################################
# MÉTHODE :radial - Y = a + b ln(d) dans un secteur centré sur le croissant
##################################################################################################################
function model_radial(chans, lum, sat_thr, g::MoonGeom, edges, refs)
    ny, nx = size(chans[1])
    nch = length(chans)
    nr = length(edges) - 1
    dc = bin_centers(edges)
    cells = [Float32[] for _ in 1:nr, _ in 1:nch]
    for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
        lum[y, x] >= sat_thr && continue
        d, θ = dist_angle(g, x, y)
        abs(wrap180(θ - g.θc)) <= RADIAL_HALF_ANGLE || continue
        k = dist_bin(d, edges)
        k == 0 && continue
        for c in 1:nch
            push!(cells[k, c], chans[c][y, x])
        end
    end

    E = [zeros(Float32, ny, nx) for _ in 1:nch]
    for c in 1:nch
        prof = [length(cells[k, c]) >= MIN_PIXELS ? Float64(median(cells[k, c])) : NaN for k in 1:nr]
        ok = findall(isfinite, prof)
        length(ok) >= 3 || error("profil radial : pas assez de ciel dans le secteur du croissant")
        exc = prof .- refs[c]
        if RADIAL_FIT_DMAX === nothing
            e0 = maximum(exc[ok])
            sel = [k for k in ok if exc[k] > 0.02 * e0]
        else
            sel = [k for k in ok if dc[k] <= RADIAL_FIT_DMAX]
        end
        length(sel) >= 3 || (println("canal $c : pas assez de points pour l'ajustement, pas de correction"); continue)
        X = hcat(ones(length(sel)), log.(dc[sel]))
        a, b = X \ prof[sel]
        @printf("canal %d : Y = %.1f ln(d) + %.1f  (d = distance au bord en pixels, de %.0f à %.0f)\n",
                c, b, a, dc[sel[1]], dc[sel[end]])
        if b >= 0
            println("   pente positive : pas de halo décroissant détecté, pas de correction")
            continue
        end
        d0 = exp((refs[c] - a) / b)
        @printf("   le modèle rejoint le ciel de référence (%.1f ADU) à d = %.0f pixels\n", refs[c], d0)

        Ec = E[c]
        H = RADIAL_HALF_ANGLE
        T = RADIAL_TAPER
        @inbounds for x in 1:nx, y in 1:ny
            lum[y, x] >= sat_thr && continue
            d, θ = dist_angle(g, x, y)
            wm = moon_weight(d, θ)                     # part "Lune" : pas de correction du ciel
            wm >= 1 && continue
            Δ = abs(wrap180(θ - g.θc))
            w = Δ <= H ? 1.0 : (T > 0 && Δ <= H + T ? cos(π / 2 * (Δ - H) / T)^2 : 0.0)
            w == 0 && continue
            e = a + b * log(max(d, edges[1])) - refs[c]
            e > 0 && (Ec[y, x] = (1 - wm) * w * e)
        end
    end
    return E
end

##################################################################################################################
# MÉTHODE :polar - carte empirique du halo (distance au bord x angle)
##################################################################################################################

# remplace les NaN d'un vecteur par interpolation linéaire (valeur constante aux extrémités)
function fill_along!(v::AbstractVector{Float64})
    idx = findall(!isnan, v)
    isempty(idx) && return false
    for k in eachindex(v)
        isnan(v[k]) || continue
        j = searchsortedfirst(idx, k)
        if j > length(idx)
            v[k] = v[idx[end]]
        elseif j == 1
            v[k] = v[idx[1]]
        else
            a = idx[j-1]
            b = idx[j]
            t = (k - a) / (b - a)
            v[k] = (1 - t) * v[a] + t * v[b]
        end
    end
    return true
end

# cases vides : interpolation le long de la distance, puis secteur valide le plus proche
function fill_grid!(G::Matrix{Float64})
    nr, ns = size(G)
    has = [fill_along!(view(G, :, s)) for s in 1:ns]
    any(has) || error("aucune mesure du fond de ciel autour de la Lune")
    for s in 1:ns
        has[s] && continue
        best = 0
        bd = typemax(Int)
        for t in 1:ns
            has[t] || continue
            dd = min(abs(t - s), ns - abs(t - s))
            dd < bd && (bd = dd; best = t)
        end
        G[:, s] .= G[:, best]
    end
    return G
end

# médianes glissantes : sur l'angle (circulaire) puis sur la distance
function smooth_grid(G::Matrix{Float64}, hr::Int, hθ::Int)
    nr, ns = size(G)
    H = similar(G)
    for k in 1:nr, s in 1:ns
        H[k, s] = median([G[k, mod1(s + j, ns)] for j in -hθ:hθ])
    end
    K = similar(G)
    for s in 1:ns, k in 1:nr
        K[k, s] = median(H[max(1, k - hr):min(nr, k + hr), s])
    end
    return K
end

function model_polar(chans, lum, sat_thr, g::MoonGeom, edges, refs)
    ny, nx = size(chans[1])
    nch = length(chans)
    nr = length(edges) - 1
    ns = POLAR_N_SECTORS
    wθ = 360.0 / ns
    nr >= 2 || error("N_DIST_BINS doit être >= 2")

    # 1) échantillons du ciel dans chaque case (distance, secteur) et chaque canal
    cells = [Float32[] for _ in 1:nr, _ in 1:ns, _ in 1:nch]
    for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
        lum[y, x] >= sat_thr && continue
        d, θ = dist_angle(g, x, y)
        k = dist_bin(d, edges)
        k == 0 && continue
        s = clamp(floor(Int, (θ + 180.0) / wθ) + 1, 1, ns)
        for c in 1:nch
            push!(cells[k, s, c], chans[c][y, x])
        end
    end

    ldmin = log(edges[1])
    lstep = (log(edges[end]) - ldmin) / nr
    E = Vector{Matrix{Float32}}(undef, nch)
    for c in 1:nch
        # 2) grille de l'excès de fond de ciel par rapport à la référence
        G = fill(NaN, nr, ns)
        for s in 1:ns, k in 1:nr
            v = cells[k, s, c]
            length(v) >= MIN_PIXELS && (G[k, s] = median(v) - refs[c])
        end
        @printf("canal %d : %d cases mesurées sur %d\n", c, count(!isnan, G), nr * ns)
        fill_grid!(G)
        G = smooth_grid(G, POLAR_SMOOTH_R, POLAR_SMOOTH_THETA)
        CLAMP_POSITIVE && (G .= max.(G, 0.0))

        # 3) interpolation bilinéaire (log(distance), angle) sur chaque pixel du ciel
        Ec = zeros(Float32, ny, nx)
        @inbounds for x in 1:nx, y in 1:ny
            lum[y, x] >= sat_thr && continue
            d, θ = dist_angle(g, x, y)
            wm = moon_weight(d, θ)                     # part "Lune" : pas de correction du ciel
            wm >= 1 && continue
            tr = (log(max(d, edges[1])) - ldmin) / lstep + 0.5     # = k au centre de l'intervalle k
            k0 = clamp(floor(Int, tr), 1, nr - 1)
            wr = clamp(tr - k0, 0.0, 1.0)
            ts = (θ + 180.0) / wθ + 0.5                             # = s au centre du secteur s
            s0 = floor(Int, ts)
            ws = ts - s0
            s0 = mod1(s0, ns)
            s1 = mod1(s0 + 1, ns)
            Ec[y, x] = (1 - wm) * ((1 - wr) * ((1 - ws) * G[k0, s0] + ws * G[k0, s1]) +
                                   wr * ((1 - ws) * G[k0+1, s0] + ws * G[k0+1, s1]))
        end
        E[c] = Ec
    end
    return E
end

##################################################################################################################
# MÉTHODE :psf - convolution du croissant par les ailes de la PSF
##################################################################################################################

# noyau K(ρ) = (ρ² + ρc²)^(-α/2) sur la grille réduite (ρ en pixels de l'image d'origine), centré en
# (1, 1) avec périodicité (convolution circulaire sur une grille 2x plus grande : pas de repliement)
function kernel_fft(α::Float64, Py::Int, Px::Int, f::Int)
    core2 = max(PSF_CORE, 0.5 * f)^2
    K = Matrix{Float64}(undef, Py, Px)
    @inbounds for j in 1:Px, i in 1:Py
        dy = min(i - 1, Py - i + 1) * f
        dx = min(j - 1, Px - j + 1) * f
        K[i, j] = (Float64(dx * dx + dy * dy) + core2)^(-α / 2)
    end
    return rfft(K)
end

conv_lowres(Sf, Kf, Py, nyl, nxl) = irfft(Sf .* Kf, Py)[1:nyl, 1:nxl]

# régression robuste y = B + A x (rejet itératif à PSF_CLIP sigma) ; score = moyenne des 90 %
# plus petits carrés des résidus (insensible aux étoiles restantes)
function fit_lin_robust(x::Vector{Float64}, y::Vector{Float64})
    keep = trues(length(x))
    A, B = 0.0, 0.0
    for _ in 1:6
        xs = x[keep]
        ys = y[keep]
        mx = mean(xs)
        my = mean(ys)
        sxx = sum(abs2, xs .- mx)
        sxx > 0 || break
        A = sum((xs .- mx) .* (ys .- my)) / sxx
        B = my - A * mx
        r = y .- (B .+ A .* x)
        rk = r[keep]
        σ = 1.4826 * median(abs.(rk .- median(rk)))
        σ > 0 || break
        newkeep = abs.(r) .< PSF_CLIP * σ
        (newkeep == keep || count(newkeep) < 10) && break
        keep = newkeep
    end
    r2 = sort((y .- (B .+ A .* x)) .^ 2)
    score = mean(view(r2, 1:max(1, floor(Int, 0.9 * length(r2)))))
    return A, B, score
end

#############################################################################################
# distance (en blocs) de chaque bloc au bloc "source" le plus proche : transformée de distance
# par chanfrein 3x3 (deux passes, erreur < 8 %), suffisante pour doser la correction
#############################################################################################
function chamfer_distance(mask::BitMatrix)
    ny, nx = size(mask)
    D = fill(Inf, ny, nx)
    D[mask] .= 0.0
    s2 = sqrt(2.0)
    for j in 1:nx, i in 1:ny                     # passe directe : voisins déjà traités
        v = D[i, j]
        i > 1 && (v = min(v, D[i-1, j] + 1))
        if j > 1
            v = min(v, D[i, j-1] + 1)
            i > 1 && (v = min(v, D[i-1, j-1] + s2))
            i < ny && (v = min(v, D[i+1, j-1] + s2))
        end
        D[i, j] = v
    end
    for j in nx:-1:1, i in ny:-1:1               # passe inverse
        v = D[i, j]
        i < ny && (v = min(v, D[i+1, j] + 1))
        if j < nx
            v = min(v, D[i, j+1] + 1)
            i < ny && (v = min(v, D[i+1, j+1] + s2))
            i > 1 && (v = min(v, D[i-1, j+1] + s2))
        end
        D[i, j] = v
    end
    return D
end

# poids de la correction sur le disque en fonction de la distance t au croissant (pixels)
function disk_weight(t::Float64)
    if PSF_DISK_RAMP <= 0 || t >= PSF_DISK_RAMP
        return PSF_DISK_SCALE
    end
    return PSF_DISK_SCALE * (PSF_DISK_W0 + (1 - PSF_DISK_W0) * sin(π / 2 * t / PSF_DISK_RAMP)^2)
end

# carte réduite de la distance au croissant (pixels), mémorisée pour la coupe de contrôle
const TERM_DIST = Ref(zeros(0, 0))

function model_psf(chans, lum, sat_thr, g::MoonGeom, ref_test, refs;
                   correct_disk::Bool=PSF_CORRECT_DISK, disk_only::Bool=false)
    ny, nx = size(chans[1])
    nch = length(chans)
    f = PSF_DOWNSAMPLE
    nyl = ny ÷ f
    nxl = nx ÷ f
    src_thr = PSF_SOURCE_FRACTION * maximum(lum)

    # 1) image réduite : source (croissant, moyenne des blocs) et ciel (médiane des blocs 100 % ciel)
    S = [zeros(nyl, nxl) for _ in 1:nch]
    Y = [fill(NaN, nyl, nxl) for _ in 1:nch]
    isref = falses(nyl, nxl)
    issrc = falses(nyl, nxl)
    buf = [Vector{Float32}(undef, f * f) for _ in 1:nch]
    for bj in 1:nxl, bi in 1:nyl
        allsky = true
        for j in 1:f, i in 1:f
            y = (bi - 1) * f + i
            x = (bj - 1) * f + j
            d, _ = dist_angle(g, x, y)
            lv = lum[y, x]
            if d < 0 && lv >= src_thr
                for c in 1:nch
                    S[c][bi, bj] += chans[c][y, x]
                end
                issrc[bi, bj] = true
                allsky = false
            elseif d < LIMB_MARGIN || lv >= sat_thr
                allsky = false
            end
            n = (j - 1) * f + i
            for c in 1:nch
                buf[c][n] = chans[c][y, x]
            end
        end
        for c in 1:nch
            S[c][bi, bj] /= f * f
        end
        if allsky
            for c in 1:nch
                Y[c][bi, bj] = median(buf[c])
            end
            isref[bi, bj] = ref_test((bj - 0.5) * f + 0.5, (bi - 0.5) * f + 0.5)
        end
    end
    idx = findall(!isnan, Y[1])
    length(idx) >= 100 || error("PSF : pas assez de blocs de ciel pour l'ajustement")
    @printf("image réduite %d x %d (blocs %dx%d), %d blocs de ciel pour l'ajustement\n",
            nxl, nyl, f, f, length(idx))

    # 2) FFT des sources sur une grille 2x plus grande
    Py = nextprod([2, 3, 5], 2 * nyl)
    Px = nextprod([2, 3, 5], 2 * nxl)
    Sf = map(S) do s
        P = zeros(Py, Px)
        P[1:nyl, 1:nxl] .= s
        rfft(P)
    end
    yv = [Y[c][idx] for c in 1:nch]

    # 3) recherche de α (grille au pas 0,1), amplitude A et fond B par régression robuste
    best = [(Inf, 0.0, 0.0, 0.0) for _ in 1:nch]       # (score, α, A, B)
    try_alpha!(c, α, Kf) = begin
        C = conv_lowres(Sf[c], Kf, Py, nyl, nxl)
        A, B, sc = fit_lin_robust(C[idx], yv[c])
        (A > 0 && sc < best[c][1]) && (best[c] = (sc, α, A, B))
    end
    αs = PSF_ALPHA_RANGE[1]:0.1:PSF_ALPHA_RANGE[2]
    print("ajustement de α : ")
    for α in αs
        Kf = kernel_fft(α, Py, Px, f)
        for c in 1:nch
            try_alpha!(c, α, Kf)
        end
        print(".")
    end
    println()
    # affinage au pas 0,02 autour du meilleur α de chaque canal
    for c in 1:nch
        α0 = best[c][2]
        α0 > 0 || error("PSF : ajustement impossible pour le canal $c (amplitude négative)")
        for α in (α0 - 0.08):0.02:(α0 + 0.08)
            try_alpha!(c, α, kernel_fft(α, Py, Px, f))
        end
    end

    # 4) modèle de la lumière diffusée, ramené à zéro dans la zone de référence, agrandi à la
    #    taille de l'image (interpolation bilinéaire entre centres des blocs)
    # distance au croissant (pixels) et poids de la correction sur le disque, sur la grille réduite
    Dt = chamfer_distance(issrc) .* f
    TERM_DIST[] = Dt
    Wd = disk_weight.(Dt)
    if correct_disk
        @printf("dosage sur le disque : ×%.2f au terminateur, ×%.2f au-delà de %.0f pixels\n",
                disk_weight(0.0), disk_weight(Inf), PSF_DISK_RAMP)
    end

    E = Vector{Matrix{Float32}}(undef, nch)
    for c in 1:nch
        sc, α, A, B = best[c]
        M = A .* conv_lowres(Sf[c], kernel_fft(α, Py, Px, f), Py, nyl, nxl)
        refblocks = M[isref]
        Mref_measured = isempty(refblocks) ? 0.0 : median(refblocks)
        isempty(refblocks) && println("ATTENTION : aucun bloc de ciel dans la zone de référence")
        Mref = PSF_ZERO_LEVEL == :model ? 0.0 : Mref_measured
        @printf("canal %d : α = %.2f, fond ajusté B = %.1f ADU (référence %.1f), diffusion dans la zone de référence %.1f ADU, résidu rms %.1f ADU\n",
                c, α, B, refs[c], Mref_measured, sqrt(sc))
        (α <= PSF_ALPHA_RANGE[1] + 0.05 || α >= PSF_ALPHA_RANGE[2] - 0.05) &&
            println("   ATTENTION : α en limite de PSF_ALPHA_RANGE, élargir la plage")
        El = M .- Mref
        CLAMP_POSITIVE && (El .= max.(El, 0.0))

        # estimation de la lumière diffusée sur la lumière cendrée (mi-rayon, côté opposé au croissant)
        xe = g.cx - 0.5 * g.R * cosd(g.θc)
        ye = g.cy + 0.5 * g.R * sind(g.θc)
        bi = clamp(round(Int, (ye - 0.5) / f + 0.5), 1, nyl)
        bj = clamp(round(Int, (xe - 0.5) / f + 0.5), 1, nxl)
        xi = clamp(round(Int, xe), 1, nx)
        yi = clamp(round(Int, ye), 1, ny)
        raw = chans[c][yi, xi] - refs[c]
        @printf("   lumière cendrée à mi-rayon (x=%d, y=%d) : signal %.1f ADU au-dessus du ciel, dont diffusion estimée %.1f ADU (%.0f %%), correction dosée %.1f ADU\n",
                xi, yi, raw, El[bi, bj], raw > 0 ? 100 * El[bi, bj] / raw : NaN, Wd[bi, bj] * El[bi, bj])

        Eld = El .* Wd                      # version dosée, utilisée sur le disque uniquement
        Ec = zeros(Float32, ny, nx)
        @inbounds for x in 1:nx
            tx = (x - 0.5) / f + 0.5
            j0 = clamp(floor(Int, tx), 1, nxl - 1)
            wx = clamp(tx - j0, 0.0, 1.0)
            for y in 1:ny
                lv = lum[y, x]
                lv >= sat_thr && continue
                d, θ = dist_angle(g, x, y)
                # répartition disque / ciel autour du vrai bord (fondu de LIMB_BLEND pixels)
                wm = moon_weight(d, θ)
                wdisk = (correct_disk && lv < src_thr) ? wm : 0.0
                wsky = disk_only ? 0.0 : 1.0 - wm
                (wdisk > 0 || wsky > 0) || continue
                ty = (y - 0.5) / f + 0.5
                i0 = clamp(floor(Int, ty), 1, nyl - 1)
                wy = clamp(ty - i0, 0.0, 1.0)
                v = 0.0
                if wsky > 0
                    v += wsky * ((1 - wy) * ((1 - wx) * El[i0, j0] + wx * El[i0, j0+1]) +
                                 wy * ((1 - wx) * El[i0+1, j0] + wx * El[i0+1, j0+1]))
                end
                if wdisk > 0
                    v += wdisk * ((1 - wy) * ((1 - wx) * Eld[i0, j0] + wx * Eld[i0, j0+1]) +
                                  wy * ((1 - wx) * Eld[i0+1, j0] + wx * Eld[i0+1, j0+1]))
                end
                Ec[y, x] = v
            end
        end
        E[c] = Ec
    end
    return E
end

##################################################################################################################
# CONTRÔLE : coupe photométrique (médiane dans un secteur étroit) et graphiques
##################################################################################################################
function wedge_profile(chans, lum, sat_thr, g::MoonGeom, θ0, halfw, edges)
    ny, nx = size(chans[1])
    nr = length(edges) - 1
    cells = [Float64[] for _ in 1:nr]
    for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
        lum[y, x] >= sat_thr && continue
        d, θ = dist_angle(g, x, y)
        abs(wrap180(θ - θ0)) <= halfw || continue
        k = dist_bin(d, edges)
        k == 0 && continue
        push!(cells[k], lumval(chans, y, x))
    end
    return [length(v) >= MIN_PIXELS ? median(v) : NaN for v in cells]
end

# affichage des graphiques autorisé (v1.2 : premier fichier seulement, voir SHOW_PLOTS_FIRST_ONLY)
const SHOW_PLOTS = Ref(true)

# enregistrement d'un graphique en PNG (si SAVE_PNG = true)
function save_plot(plt, path)
    SAVE_PNG || return nothing
    try
        savefig(plt, path)
        println("enregistré : ", path)
    catch err
        println("ERREUR d'enregistrement du graphique : ", sprint(showerror, err))
    end
    return nothing
end

function show_plot(plt)
    SHOW_PLOTS[] || return nothing
    try
        PLOT_WINDOW == :gui ? gui(plt) : display(plt)
    catch err
        println("ERREUR d'affichage du graphique : ", sprint(showerror, err))
    end
    return nothing
end

#############################################################################################
# coupe de contrôle SUR LA LUMIÈRE CENDRÉE : médiane (luminance) des pixels du disque hors
# croissant, en fonction de la distance t au croissant (pas DISK_PROFILE_STEP). Les pixels à
# moins de LIMB_MARGIN du bord du disque sont exclus.
#############################################################################################
function disk_profile(M::Vector{Matrix{Float32}}, lum, src_thr, g::MoonGeom, Dt::Matrix{Float64}, f::Int, nb::Int)
    ny, nx = size(M[1])
    nyl, nxl = size(Dt)
    cells = [Float64[] for _ in 1:nb]
    for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
        lum[y, x] >= src_thr && continue
        d, _ = dist_angle(g, x, y)
        d < -LIMB_MARGIN || continue
        t = Dt[clamp(cld(y, f), 1, nyl), clamp(cld(x, f), 1, nxl)]
        isfinite(t) || continue
        k = floor(Int, t / DISK_PROFILE_STEP) + 1
        k <= nb || continue
        push!(cells[k], lumval(M, y, x))
    end
    return [length(v) >= MIN_PIXELS ? median(v) : NaN for v in cells]
end

function plot_control(path, edges, before, after, opposite, ref_lum, E, g::MoonGeom, title_str;
                      disk=nothing)
    dc = bin_centers(edges)
    ok(v) = findall(isfinite, v)
    p1 = plot(dc[ok(before)], before[ok(before)]; xscale=:log10, lw=2, color=:red,
              label="côté croissant, avant", xlabel="distance au bord de la Lune (pixels)",
              ylabel="fond de ciel (ADU, luminance)", title="coupe photométrique", legend=:topright)
    plot!(p1, dc[ok(after)], after[ok(after)]; lw=2, color=:green, label="côté croissant, après")
    plot!(p1, dc[ok(opposite)], opposite[ok(opposite)]; lw=1.5, color=:blue, ls=:dot,
          label="côté lumière cendrée")
    hline!(p1, [ref_lum]; color=:black, ls=:dash, label="ciel de référence")

    # carte de la lumière soustraite (luminance, 1 pixel sur 8)
    ny, nx = size(E[1])
    st = 8
    Em = [lumval(E, y, x) for y in 1:st:ny, x in 1:st:nx]
    p2 = heatmap(1:st:nx, 1:st:ny, Em; yflip=true, aspect_ratio=:equal, color=:inferno,
                 title="lumière diffusée soustraite (ADU)", xlabel="x", ylabel="y",
                 xlims=(1, nx), ylims=(1, ny))
    plot!(p2, [g.cx .+ g.R .* cosd(t) for t in 0:2:360], [g.cy .- g.R .* sind(t) for t in 0:2:360];
          color=:cyan, lw=1, label="")
    if disk === nothing
        plt = plot(p1, p2; layout=(1, 2), size=(1600, 650), plot_title=title_str,
                   left_margin=5Plots.mm, bottom_margin=5Plots.mm)
    else
        # coupe sur la lumière cendrée en fonction de la distance au croissant
        okb = findall(isfinite, disk.before)
        p3 = plot(disk.t[okb], disk.before[okb]; lw=2, color=:red, label="avant",
                  xlabel="distance au croissant / terminateur (pixels)", ylabel="ADU (luminance)",
                  title="lumière cendrée", legend=:topright)
        oka = findall(isfinite, disk.after)
        plot!(p3, disk.t[oka], disk.after[oka]; lw=2, color=:green, label="après")
        okm = findall(isfinite, disk.model)
        plot!(p3, disk.t[okm], disk.model[okm]; lw=1.5, color=:orange, label="correction soustraite")
        PSF_DISK_RAMP > 0 && vline!(p3, [PSF_DISK_RAMP]; color=:gray, ls=:dash, label="PSF_DISK_RAMP")
        hline!(p3, [ref_lum]; color=:black, ls=:dash, label="ciel de référence")
        plt = plot(p1, p2, p3; layout=(1, 3), size=(2300, 650), plot_title=title_str,
                   left_margin=5Plots.mm, bottom_margin=5Plots.mm)
    end
    save_plot(plt, path)
    show_plot(plt)
    return nothing
end

##################################################################################################################
# PROGRAMME PRINCIPAL
##################################################################################################################
#############################################################################################
# traitement complet d'UN fichier : lecture, géométrie, modèle, correction, enregistrements.
#   show_plots : affichage des graphiques (premier fichier seulement par défaut) ; les fichiers
#                PNG des graphiques sont toujours enregistrés, pour chaque image
# Les résultats sont écrits dans le répertoire de l'image d'entrée.
#############################################################################################
function process_file(selected_input::AbstractString; show_plots::Bool=true)

SHOW_PLOTS[] = show_plots
out_dir = dirname(selected_input)
file_name_input = basename(selected_input)
stem = splitext(file_name_input)[1]
# états globaux réinitialisés pour chaque image
LIMB_EDGE[] = Float64[]
TERM_DIST[] = zeros(0, 0)

# -----------------------------------------------------------------------------------------------
# lecture de l'image
# -----------------------------------------------------------------------------------------------
println("lecture du fichier d'entrée : ", file_name_input)
chans, bits = read_image(selected_input)
ny, nx = size(chans[1])
nch = length(chans)
println("dim_x = ", nx, ", dim_y = ", ny)

# image de détection : maximum des canaux (un pixel est "saturé" dès qu'un canal l'est)
lum = nch == 3 ? max.(chans[1], chans[2], chans[3]) : copy(chans[1])
maxv = maximum(lum)
sat_thr = Float32(SAT_FRACTION * maxv)
@printf("maximum de l'image : %.0f ADU (%.1f %% de la pleine échelle)\n", maxv, 100 * maxv / (2.0^bits - 1))
maxv < 0.9 * (2.0^bits - 1) && println("ATTENTION : le croissant ne semble pas saturé")

# -----------------------------------------------------------------------------------------------
# géométrie de la Lune
# -----------------------------------------------------------------------------------------------
if MOON_CENTER === nothing || MOON_RADIUS === nothing
    cx, cy, R = detect_moon(lum, Float32(LIMB_THRESHOLD * maxv), CRESCENT_SIDE)
else
    cx, cy, R = 0.0, 0.0, 0.0
end
MOON_CENTER === nothing || ((cx, cy) = Float64.(MOON_CENTER))
MOON_RADIUS === nothing || (R = Float64(MOON_RADIUS))
θc = crescent_direction(lum, sat_thr, cx, cy)
g = MoonGeom(cx, cy, R, θc)
@printf("Lune : centre x = %.1f, y = %.1f (origine en haut à gauche), rayon %.1f pixels (diamètre %.0f)\n",
        cx, cy, R, 2R)
@printf("direction du croissant : %.1f° (0° = droite, 90° = haut)\n", θc)

edges = dist_edges(nx, ny, g)

# -----------------------------------------------------------------------------------------------
# ciel de référence (côté lumière cendrée)
# -----------------------------------------------------------------------------------------------
ref_test = make_ref_test(nx, ny, g)
refvals = [Float64[] for _ in 1:nch]
for x in 1:STAT_STEP:nx, y in 1:STAT_STEP:ny
    lum[y, x] >= sat_thr && continue
    ref_test(x, y) || continue
    for c in 1:nch
        push!(refvals[c], chans[c][y, x])
    end
end
length(refvals[1]) >= 500 || error("zone de référence trop petite (", length(refvals[1]),
                                   " pixels) : régler REF_BOX ou REF_STRIP_FRAC")
refs = Float64[]
sigs = Float64[]
for c in 1:nch
    m, s = robust_background(refvals[c])
    push!(refs, m)
    push!(sigs, s)
    @printf("ciel de référence, canal %d : %.1f ADU (bruit %.1f ADU, %d pixels)\n", c, m, s, length(refvals[c]))
end
ref_lum = nch == 3 ? sum(LUMA_WEIGHTS .* Tuple(refs)) : refs[1]
σ_lum = nch == 3 ? sqrt(sum((LUMA_WEIGHTS .* Tuple(sigs)) .^ 2)) : sigs[1]
src_thr = PSF_SOURCE_FRACTION * maxv

# -----------------------------------------------------------------------------------------------
# vrai bord de la Lune (v1.1) : position de la transition correction du disque / du ciel
# -----------------------------------------------------------------------------------------------
LIMB_EDGE[] = Float64[]                    # (réinitialisé à chaque lancement dans le REPL)
if LIMB_EDGE_DETECT
    edge_sm, edge_raw = detect_true_limb(chans, g, src_thr, σ_lum)
    LIMB_EDGE[] = edge_sm
else
    edge_sm, edge_raw = Float64[], Float64[]
    println("vrai bord non détecté (LIMB_EDGE_DETECT = false) : transition sur le cercle ajusté")
end

# coupe photométrique AVANT correction
before = wedge_profile(chans, lum, sat_thr, g, θc, PROFILE_HALF_ANGLE, edges)
opposite = wedge_profile(chans, lum, sat_thr, g, θc + 180.0, PROFILE_HALF_ANGLE, edges)

# -----------------------------------------------------------------------------------------------
# modèle de la lumière diffusée
# -----------------------------------------------------------------------------------------------
println("méthode : ", METHOD)
t0 = time()
E = if METHOD == :radial
    model_radial(chans, lum, sat_thr, g, edges, refs)
elseif METHOD == :polar
    model_polar(chans, lum, sat_thr, g, edges, refs)
elseif METHOD == :psf
    model_psf(chans, lum, sat_thr, g, ref_test, refs)
elseif METHOD == :hybrid
    println("ciel : méthode :polar")
    Ep = model_polar(chans, lum, sat_thr, g, edges, refs)
    println("lumière cendrée : méthode :psf (dosée)")
    Ed = model_psf(chans, lum, sat_thr, g, ref_test, refs; correct_disk=true, disk_only=true)
    [Ep[c] .+ Ed[c] for c in 1:nch]          # zones disjointes : ciel (polar) + disque (psf)
else
    error("METHOD inconnue : ", METHOD, " (attendu :radial, :polar, :psf ou :hybrid)")
end
@printf("modèle calculé en %.1f s\n", time() - t0)

# coupe de contrôle sur la lumière cendrée AVANT correction (méthodes :psf et :hybrid)
has_disk = METHOD in (:psf, :hybrid) && !isempty(TERM_DIST[])
nbd = ceil(Int, 2 * R / DISK_PROFILE_STEP)
disk_before = has_disk ? disk_profile(chans, lum, src_thr, g, TERM_DIST[], PSF_DOWNSAMPLE, nbd) : Float64[]
disk_model = has_disk ? disk_profile(E, lum, src_thr, g, TERM_DIST[], PSF_DOWNSAMPLE, nbd) : Float64[]

# -----------------------------------------------------------------------------------------------
# correction (le croissant saturé n'est jamais modifié : E = 0 sur ces pixels)
# -----------------------------------------------------------------------------------------------
for c in 1:nch
    ch = chans[c]
    Ec = E[c]
    @inbounds for i in eachindex(ch)
        ch[i] = max(ch[i] - Ec[i], 0f0)
    end
end
after = wedge_profile(chans, lum, sat_thr, g, θc, PROFILE_HALF_ANGLE, edges)

# tableau de contrôle
dc = bin_centers(edges)
println("\ncoupe côté croissant (luminance, ADU) :")
println("   distance    avant    après   référence")
for d in (25, 50, 100, 200, 500, 1000, 2000, 4000)
    k = dist_bin(Float64(d), edges)
    k == 0 && continue
    @printf("   %6.0f  %8.1f %8.1f %8.1f\n", dc[k], before[k], after[k], ref_lum)
end
println()

disk = nothing
if has_disk
    disk_after = disk_profile(chans, lum, src_thr, g, TERM_DIST[], PSF_DOWNSAMPLE, nbd)
    tcent = [(k - 0.5) * DISK_PROFILE_STEP for k in 1:nbd]
    disk = (t=tcent, before=disk_before, after=disk_after, model=disk_model)
    println("lumière cendrée en fonction de la distance au croissant (luminance, ADU) :")
    println("   distance    avant   correction    après")
    for k in 1:nbd
        (k <= 20 ? isodd(k) : k % 6 == 0) || continue
        isfinite(disk_before[k]) || continue
        @printf("   %6.0f  %8.1f %10.1f %8.1f\n", tcent[k], disk_before[k], disk_model[k], disk_after[k])
    end
    println()
end

# -----------------------------------------------------------------------------------------------
# enregistrements
# -----------------------------------------------------------------------------------------------
tag = string(METHOD)
save_tiff(joinpath(out_dir, stem * "_corr_" * tag * ".tif"), chans, bits)
SAVE_MODEL && save_tiff(joinpath(out_dir, stem * "_diffusion_" * tag * ".tif"), E, bits)
# graphiques : construits seulement s'ils sont enregistrés (SAVE_PNG) ou affichés
if SAVE_PNG || SHOW_PLOTS[]
    plot_control(joinpath(out_dir, stem * "_diffusion_" * tag * ".png"), edges, before, after, opposite,
                 ref_lum, E, g, file_name_input * " - méthode " * tag; disk=disk)
    LIMB_EDGE_DETECT && plot_limb_edge(joinpath(out_dir, stem * "_bord_" * tag * ".png"), edge_sm, edge_raw, g,
                                       file_name_input * " - vrai bord de la Lune")
end

return nothing

end # fin de la fonction process_file()

##################################################################################################################
# PROGRAMME PRINCIPAL : sélection de plusieurs fichiers, traitement l'un après l'autre
##################################################################################################################
function main()

println("diffusion_lumiere_cendree - version 1.3")

# -----------------------------------------------------------------------------------------------
# sélection interactive du répertoire de travail et des fichiers d'entrée (TIFF ou NEF)
# -----------------------------------------------------------------------------------------------
work_dir = pick_folder(pwd())
isempty(work_dir) && error("Aucun répertoire sélectionné - programme arrêté")
cd(work_dir)
println(pwd())

exts = vcat(["tif", "tiff"], RAW_EXTENSIONS)
filterlist = join(vcat(exts, uppercase.(exts)), ",")      # majuscules : sélecteurs sensibles à la casse (Linux)
files = pick_multi_file(work_dir; filterlist=filterlist)
isempty(files) && error("Aucun fichier sélectionné - programme arrêté")
files = sort(String.(files))
if SKIP_OUTPUT_FILES
    is_output = f -> occursin(r"_(corr|diffusion)_(radial|polar|psf|hybrid)\.tiff?$"i, basename(f))
    skipped = filter(is_output, files)
    for f in skipped
        println("ignoré (fichier produit par ce script) : ", basename(f))
    end
    files = filter(!is_output, files)
    isempty(files) && error("Aucun fichier à traiter - programme arrêté")
end
n = length(files)
println(n, " fichier(s) à traiter, méthode ", METHOD)

# fichiers RAW : recherche du convertisseur avant tout traitement (arrêt immédiat s'il manque)
nraw = count(is_raw_file, files)
if nraw > 0
    conv = raw_converter()
    println(nraw, " fichier(s) RAW, convertis par : ", conv)
    println("   options : ", join(raw_args(), " "))
    println("   retournement de l'image convertie : ",
            RAW_FLIP_VERTICAL && RAW_FLIP_HORIZONTAL ? "haut / bas et gauche / droite" :
            RAW_FLIP_VERTICAL ? "haut / bas" : RAW_FLIP_HORIZONTAL ? "gauche / droite" : "aucun")
end

# -----------------------------------------------------------------------------------------------
# traitement de chaque fichier (une erreur sur un fichier n'arrête pas les suivants,
# sauf si STOP_ON_ERROR = true)
# -----------------------------------------------------------------------------------------------
results = Tuple{String, String, Float64}[]
t_all = time()
for (i, f) in enumerate(files)
    println()
    println("="^100)
    @printf("fichier %d / %d : %s\n", i, n, basename(f))
    println("="^100)
    t_file = time()
    status = "OK"
    try
        process_file(f; show_plots=(!SHOW_PLOTS_FIRST_ONLY || i == 1))
    catch err
        err isa InterruptException && rethrow()
        status = "ERREUR : " * first(replace(sprint(showerror, err), '\n' => ' '), 150)
        println(status)
        STOP_ON_ERROR && rethrow()
    end
    push!(results, (basename(f), status, time() - t_file))
    GC.gc()                     # libère la mémoire de l'image précédente (images de grande taille)
end
SHOW_PLOTS[] = true

# -----------------------------------------------------------------------------------------------
# bilan
# -----------------------------------------------------------------------------------------------
println()
println("="^100)
@printf("bilan : %d fichier(s), %d traité(s) sans erreur, durée totale %.0f s\n",
        n, count(r -> r[2] == "OK", results), time() - t_all)
for (name, status, dt) in results
    @printf("   %-45s %6.0f s   %s\n", name, dt, status)
end
println("="^100)

# sinon la fenêtre Plots se ferme à la fin du programme (lancement en ligne de commande)
# (inutile et source de blocage dans le REPL de VS Code : on n'attend que hors session interactive)
if !isinteractive()
    println("appuyer sur Entrée pour terminer")
    readline()
end
println("terminé")
return nothing

end # fin de la fonction main()

# lancement du programme
main()