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", "GLMakie"]
        Pkg.add(pkg)
    end
end

using Statistics
using Printf
using Random                    # bibliothèque standard (tirages de la recherche du cercle du limbe)
import LinearAlgebra            # bibliothèque standard (systèmes linéaires des ajustements)

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

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

# fenêtre interactive de pointage de l'étoile (déplacement, zoom, seuil, réticule)
import GLMakie
const GM = GLMakie

##################################################################################################################
# RECENTRAGE D'IMAGES DE LUMIÈRE CENDRÉE SUR UNE ÉTOILE DE RÉFÉRENCE
#
# Version 2.0 (14/10/2026)
#   - option STACK_MOON : somme des images recentrées sur la LUNE => <1re image>_stack_reg_moon.tiff
#     (en plus de la somme recentrée sur l'étoile), mesures dans <1re image>_reg_moon.csv,
#     images recentrées <image>_reg_moon.tiff si SAVE_REG = true
#   - la position grossière de la Lune (centre du croissant, v1.2) sert de point de départ ; elle est
#     affinée à mieux que le pixel par l'une des 3 méthodes suivantes (MOON_METHOD) :
#       :correlation  corrélation par moindres carrés (Lucas-Kanade) dans une fenêtre MOON_WIN x MOON_WIN
#                     placée par l'utilisateur (clic DROIT dans la fenêtre de pointage, carré cyan), de
#                     préférence à cheval sur le limbe côté lumière cendrée (limbe + détails faibles).
#                     Les deux fenêtres sont filtrées (passe-haut MOON_CORR_HIGHPASS : le fond lent, dont le
#                     halo de lumière diffusée, est retiré), puis : image(x + d) = g x référence(x) + fond plan
#                     => insensible aux variations de transparence (g) et du halo. Pixels saturés exclus.
#                     Recherche entière (corrélation normalisée) puis affinage sous-pixel (Gauss-Newton,
#                     interpolation bicubique), rejet des pixels aberrants (étoiles, pixels chauds).
#       :limb         cercle ajusté sur le limbe : bord sous-pixel (chute d'intensité la plus forte) le long
#                     de MOON_LIMB_RAYS rayons, côté lumière cendrée (défaut) et/ou côté croissant
#                     (MOON_LIMB_SIDE) ; le terminateur (montée d'intensité vers l'extérieur) est ignoré.
#                     Le bord du croissant saturé se déplace avec la transparence et la turbulence : il
#                     n'est utilisé qu'avec MOON_LIMB_SIDE = :both ou :bright.
#     Essais sur images simulées (croissant saturé, lumière cendrée faible avec détails, halo variable,
#     transparence de 70 à 110 %, turbulence de 1,2 à 2,2 px) : erreur typique 0,02 à 0,05 px pour
#     :correlation et :limb (côté lumière cendrée), jusqu'à 0,12 px si la turbulence change beaucoup ;
#     0,2 à 0,3 px pour :crescent et pour :limb avec le côté croissant.
#       :crescent     barycentre de l'isophote MOON_CRESCENT_LEVEL x plateau saturé (rapide, mais sensible
#                     aux variations de transparence et de turbulence : sert surtout de comparaison)
#     MOON_COMPARE = true : les 3 méthodes sont calculées sur toutes les images et comparées (terminal et
#     CSV) ; seule MOON_METHOD est utilisée pour la somme. Indicateur de précision : le déplacement de la
#     Lune par rapport à l'étoile (mouvement propre de la Lune, ~0,5"/s) doit varier régulièrement d'une
#     image à l'autre ; l'écart à une parabole ajustée en fonction du numéro d'image est affiché pour
#     chaque méthode (valable pour une cadence de prise de vue régulière).
#   - dispersion atmosphérique (DISPERSION = true) : écarts R - V et B - V de la Lune mesurés dans chaque
#     image par la méthode choisie, moyennés sur toutes les images, comme pour l'étoile
#   - le cercle du limbe trouvé sur la 1re image est tracé en magenta dans la fenêtre de pointage ;
#     la fenêtre de corrélation est placée par défaut sur le limbe côté lumière cendrée
#
# Version 1.2 (14/10/2026)
#   - pré-recentrage grossier optionnel sur la Lune (MOON_PRE_ALIGN) : dans chaque image, le centre du
#     croissant saturé (barycentre de la plus grande zone brillante, sur une image réduite MOON_BIN x MOON_BIN)
#     est mesuré ; la fenêtre de recherche de l'étoile est centrée sur la position PRÉDITE
#        étoile prédite = étoile (image d'ancrage) + Lune (image courante) - Lune (image d'ancrage)
#     l'image d'ancrage étant l'image précédente (TRACK_PREVIOUS = true) ou la 1re image (false).
#     La Lune se déplace par rapport aux étoiles (~0,5"/s) : avec TRACK_PREVIOUS = true, seul le mouvement
#     de la Lune entre deux images consécutives entre dans l'erreur de prédiction.
#     Si la Lune n'est pas trouvée dans une image, la fenêtre est centrée comme en version 1.1.
#     Le pré-recentrage sert uniquement à placer la fenêtre : la translation appliquée reste celle mesurée
#     sur l'étoile. Positions de la Lune affichées et ajoutées au fichier CSV.
#   - fichiers RAW : retournement vertical (RAW_FLIP_VERTICAL) et/ou horizontal (RAW_FLIP_HORIZONTAL) de
#     l'image convertie par dcraw_emu (l'image décodée était à l'envers haut / bas)
#
# Version 1.1 (14/10/2026)
#   - lecture des fichiers RAW Nikon NEF (D5300, Z 6II, Z 7II...) : conversion en TIFF 16 bits LINÉAIRE
#     par le programme externe dcraw_emu (LibRaw), voir "FICHIERS RAW NIKON (NEF)" ci-dessous.
#     TIFF et NEF peuvent être mélangés dans une même sélection (mêmes dimensions requises).
#   - le fichier somme pondérée s'appelle désormais <1re image>_stack_reg_star.tiff
#     (au lieu de <1re image>_reg_star.tiff) ; le fichier CSV garde le nom <1re image>_reg_star.csv
#
# Version 1.0 (13/10/2026)
#
#   1) sélection de n images TIFF 16 bits RVB ou NEF (pick_folder / pick_multi_file, comme
#      diffusion_lumiere_cendree v1.2) ; la 1re image (ordre alphabétique) est l'image de référence
#   2) pointage de l'étoile de référence sur la 1re image dans une fenêtre GLMakie (comme HFD_star v9) :
#      curseurs X / Y, zoom (curseur ou molette), seuil haut de visualisation, réticule au clic gauche,
#      bouton "Mesurer" (touche M), bouton "Valider" (touche V)
#   3) mesure de la position de l'étoile avec une précision sous-pixel :
#        - fond de ciel local ajusté par un PLAN dans un anneau autour de l'étoile (le fond n'est pas
#          uniforme près de la Lune : halo de lumière diffusée)
#        - délimitation de l'étoile (composante connexe au-dessus de DETECTION_SIGMA), barycentre
#        - affinage par barycentre à fenêtre gaussienne itératif (méthode "XWIN" de SExtractor),
#          non biaisé pour une étoile symétrique, précision typique 0,02 à 0,1 pixel
#        - incertitude estimée (bruit du fond seul) affichée : avertissement si > MAX_POS_ERROR
#   4) mesure automatique de l'étoile dans les autres images, dans une fenêtre de recherche de
#      SEARCH_BOX x SEARCH_BOX pixels centrée sur la position trouvée dans l'image précédente
#      (TRACK_PREVIOUS) ; parmi les étoiles de la fenêtre, on retient celle dont le flux est proche
#      de celui de l'étoile de référence (FLUX_TOLERANCE) et la plus proche de la position attendue
#   5) affichage dans le terminal des translations delta X, delta Y de chaque image
#   6) translation sous-pixel de chaque image (SHIFT_METHOD) => <image>_reg.tiff
#   7) somme pondérée des images recentrées (WEIGHT_MODE), normalisée en 16 bits (NORMALIZE)
#      => <1re image>_stack_reg_star.tiff ; tableau des mesures => <1re image>_reg_star.csv
#
# Dispersion atmosphérique (paramètre DISPERSION) :
#   DISPERSION = false : la position de l'étoile est mesurée sur la luminance
#                        L = 0,2126 R + 0,7152 V + 0,0722 B ; les 3 couches reçoivent la même translation
#   DISPERSION = true  : la position est mesurée dans les couches R, V et B. Les écarts R - V et B - V
#                        (constants d'une image à l'autre : ils ne dépendent que de la dispersion) sont
#                        moyennés sur toutes les images (moyenne avec rejet à DISP_CLIP sigma) et affichés.
#                        Translation de la couche V : delta mesuré sur V.
#                        Translation des couches R et B : delta V + correction de dispersion, la
#                        correction étant l'écart moyen V - R (resp. V - B), c'est-à-dire l'opposé de
#                        l'écart R - V (resp. B - V) : l'étoile est ramenée, dans chaque couche, sur la
#                        position de l'étoile dans la couche V de la 1re image.
#
# 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, 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.
#
# Convention de coordonnées : x vers la droite, y vers le BAS, pixel en haut à gauche = (1, 1),
#   comme dans GIMP, Photoshop... et dans diffusion_lumiere_cendree (les TIFF sont stockés de haut en bas).
#   delta X > 0 : l'image est déplacée vers la droite ; delta Y > 0 : vers le bas.
#
# Remarques :
#   - les zones qui sortent du champ après translation sont mises à 0 dans les fichiers _reg.tiff ; dans
#     la somme, chaque pixel est la moyenne pondérée des seules images qui le couvrent
#   - les métadonnées du TIFF (profil ICC, EXIF) ne sont pas recopiées
#   - les images sont lues deux fois (mesure, puis translation) pour ne garder en mémoire qu'une image à
#     la fois en plus de la somme : prévoir ~8 fois la taille d'une image en Float32. Les NEF sont donc
#     convertis deux fois, sauf si RAW_CACHE = true (TIFF temporaires conservés jusqu'à la fin)
##################################################################################################################

# -----------------------------------------------------------------------------------------------
# PARAMETRES PRINCIPAUX
# -----------------------------------------------------------------------------------------------
const DISPERSION = true         # true : mesure dans R, V, B et correction de la dispersion atmosphérique
                                 # false : mesure sur la luminance L, même translation pour R, V, B
const SEARCH_BOX = 200           # (pixels) côté de la fenêtre de recherche de l'étoile dans les images 2 à n
const TRACK_PREVIOUS = true      # true : fenêtre centrée sur la position trouvée dans l'image précédente
                                 # (suit une dérive lente) ; false : centrée sur la position de la 1re image
const SHIFT_METHOD = :bicubic    # translation sous-pixel : :bicubic (Keys, a = -0,5 : conserve la netteté),
                                 # :bilinear (lisse un peu le bruit) ou :integer (pixel entier, aucune interpolation)
const WEIGHT_MODE = :equal       # poids de chaque image dans la somme :
                                 # :equal = poids égaux (moyenne simple)
                                 # :noise = 1 / bruit² (bruit pixel à pixel de l'image, mesuré sur L ou V)
                                 # :flux  = flux de l'étoile (favorise les images les plus transparentes)
                                 # :snr   = flux / bruit² (optimal si le signal suit la transparence)
const NORMALIZE = :none           # normalisation de la somme en 16 bits :
                                 # :max  = moyenne pondérée multipliée pour que son maximum (ou le centile
                                 #         NORM_PERCENTILE) vaille 65535, MÊME facteur pour R, V, B (couleurs conservées)
                                 # :none = moyenne pondérée (même échelle que les images d'origine)
const NORM_PERCENTILE = 100.0    # centile utilisé par NORMALIZE = :max (100 = maximum ; ex. 99.99 pour
                                 # ignorer quelques pixels chauds)
const SAVE_REG = false            # true : enregistre chaque image recentrée (<image>_reg.tiff)
const SAVE_CSV = true            # true : enregistre le tableau des mesures (<1re image>_reg_star.csv)
const CSV_SEP = ";"              # séparateur de colonnes du fichier CSV
const CSV_DECIMAL_COMMA = true   # true : virgule décimale (Excel / LibreOffice en français)
const SKIP_OUTPUT_FILES = true   # true : ignore les fichiers produits par ce script (*_reg.tiff,
                                 # *_reg_star.tiff, *_stack_reg_star.tiff) s'ils ont été sélectionnés par erreur

# ---- fichiers RAW Nikon (NEF), convertis par dcraw_emu (LibRaw)
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
const RAW_CACHE = false          # true : les TIFF convertis sont gardés (dossier temporaire) entre la mesure et
                                 # la translation : 2 fois plus rapide, mais ~3 x 2 octets par pixel par image
                                 # sur le disque (~270 Mo par image de Z 7II), effacés à la fin

# ---- pointage de l'étoile sur la 1re image
const USE_GUI = true             # true : pointage dans une fenêtre ; false : position (X_STAR, Y_STAR) ci-dessous
const X_STAR = 1000.0            # position approchée de l'étoile (pixels, origine en haut à gauche)
const Y_STAR = 1000.0            # (USE_GUI = true : position initiale du réticule seulement)
const GUI_SEARCH_RADIUS = 15     # (pixels) l'étoile est cherchée à moins de cette distance du réticule
const VIEW_WIDTH = 1000          # largeur de la zone d'affichage de l'image (pixels écran)
const VIEW_HEIGHT = 700          # hauteur de la zone d'affichage de l'image (pixels écran)
const MAX_ZOOM = 16              # zoom maximal (pixels écran par pixel image)

# ---- mesure de la position de l'étoile
const STAR_BOX_HALF = 50         # demi-taille de la zone de mesure autour de l'étoile (étoile + anneau de fond)
const MAX_STAR_RADIUS = 20       # (pixels) rayon maximal de l'étoile ; au-delà, objet étendu rejeté (bord de la Lune...)
const DETECTION_SIGMA = 3.0      # seuil de délimitation de l'étoile, en sigma du bruit (image lissée 3x3)
const MIN_PEAK_SNR = 8.0         # rapport signal / bruit minimal du pic (image lissée 3x3)
const BG_PLANE = true            # true : fond de ciel local = plan ajusté dans l'anneau (gradient du halo)
                                 # false : fond constant (médiane de l'anneau)
const WINDOWED_CENTROID = true   # true : affinage par barycentre à fenêtre gaussienne (XWIN) ;
                                 # false : barycentre simple des pixels de l'étoile
const MAX_POS_ERROR = 0.5        # (pixels) avertissement si l'incertitude de position dépasse cette valeur
const FLUX_TOLERANCE = 5.0       # images 2 à n : flux de l'étoile entre flux_réf / 5 et 5 x flux_réf
const MAX_CANDIDATES = 30        # nombre maximal d'étoiles candidates examinées dans la fenêtre de recherche
const SAT_LEVEL = 65000.0        # (ADU) étoile signalée saturée si un pixel atteint ce niveau
const LUMA_WEIGHTS = (0.2126, 0.7152, 0.0722)   # luminance R, V, B (Rec. 709)

# ---- dispersion atmosphérique (DISPERSION = true)
const DISP_SEARCH_RADIUS = 10    # (pixels) recherche de l'étoile dans R et B autour de sa position dans V
const DISP_CLIP = 3.0            # moyenne des écarts R - V et B - V : rejet des valeurs à plus de N sigma
                                 # (0 = moyenne simple)

# ---- pré-recentrage grossier sur la Lune (croissant saturé)
const MOON_PRE_ALIGN = true      # true : fenêtre de recherche de l'étoile placée grâce au déplacement de la Lune
const MOON_BIN = 4               # réduction de l'image (blocs MOON_BIN x MOON_BIN moyennés) pour la détection
const MOON_LEVEL = 0.5           # seuil de la zone lunaire : fond + MOON_LEVEL x (maximum - fond)
                                 # (le croissant saturé dépasse largement ce seuil)
const MOON_MIN_AREA = 5000       # (pixels de l'image d'origine) surface minimale de la zone lunaire

# ---- somme des images recentrées sur la Lune (version 2.0)
const STACK_MOON = true          # true : crée aussi <1re image>_stack_reg_moon.tiff (images recentrées sur la Lune)
const MOON_METHOD = :limb # recentrage fin sur la Lune : :correlation, :limb ou :crescent (voir en tête)
const MOON_COMPARE = false        # true : les 3 méthodes sont calculées et comparées ; la somme utilise MOON_METHOD
const MOON_WEIGHT_MODE = :equal  # poids dans la somme Lune : :equal (poids égaux) ou :noise (1 / bruit²)
const MOON_SAT_FRAC = 0.95       # pixel saturé : au-dessus de MOON_SAT_FRAC x niveau du plateau du croissant
                                 # (plateau = maximum de l'image réduite MOON_BIN x MOON_BIN, mesuré par image et par couche)
# :correlation
const MOON_WIN = 500             # (pixels) côté de la fenêtre de corrélation
const MOON_WIN_X = 0.0           # centre de la fenêtre si USE_GUI = false (avec USE_GUI = true : position initiale
const MOON_WIN_Y = 0.0           # du carré cyan) ; 0 : automatique, sur le limbe côté lumière cendrée
const MOON_CORR_SEARCH = 6       # (pixels) recherche entière à ± cette distance de la position grossière
const MOON_CORR_SMOOTH = 1.0     # (pixels) lissage gaussien des fenêtres avant l'affinage (réduit le bruit ; 0 = aucun)
const MOON_CORR_HIGHPASS = 20.0  # (pixels) passe-haut : soustraction du fond lissé à cette échelle (pixels saturés
                                 # exclus) ; rend la corrélation insensible aux variations du halo de lumière diffusée
                                 # (0 = aucun ; essais sur images simulées : 20 réduit fortement le biais dû au halo)
const MOON_CORR_MIN_FRAC = 0.2   # fraction minimale de pixels non saturés dans la fenêtre
# :limb
const MOON_LIMB_SIDE = :earthshine  # points du limbe : :earthshine (côté lumière cendrée, recommandé), :both (les deux
                                 # côtés), :bright (croissant). Le bord du croissant saturé se déplace avec la transparence
                                 # et la turbulence (essais sur images simulées : jusqu'à 0,3 px avec :both)
const MOON_LIMB_RAYS = 1440      # nombre de rayons (pas de 0,25°)
const MOON_LIMB_RANGE = 15       # (pixels) bord cherché à ± cette distance du cercle prédit
const MOON_LIMB_MIN_SNR = 5.0    # saut d'intensité minimal au bord, en sigma
const MOON_MIN_RADIUS = 50.0     # (pixels) rayon minimal de la Lune (recherche initiale du cercle)
# :crescent
const MOON_CRESCENT_LEVEL = 0.5  # isophote utilisée (fraction du plateau ; rampe de ± 0,1 pour un résultat continu)

##################################################################################################################
# LECTURE / ÉCRITURE DES TIFF (reprises de diffusion_lumiere_cendree v1.2)
##################################################################################################################

# 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
#############################################################################################
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).
#############################################################################################
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
    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é)")

    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))
        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)
    ny = get1(257, 0)
    bps = get1(258, 1)
    comp = get1(259, 1)
    spp = get1(277, 1)
    rps = min(get1(278, ny), ny)
    planar = get1(284, 1)
    sfmt = get1(339, 1)
    @printf("TIFF : %d x %d, %d canal(aux), %d bits, compression %d, %s\n",
            nx, ny, spp, bps, comp, planar == 1 ? "entrelacé" : "planaire")
    (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

    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
            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)
        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)
const RAW_CACHE_DIRS = Dict{String, String}()         # RAW_CACHE : fichier RAW => dossier du TIFF converti

#############################################################################################
# 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)
    if RAW_CACHE && haskey(RAW_CACHE_DIRS, path)
        d = RAW_CACHE_DIRS[path]
        tifs = filter(f -> occursin(r"\.tiff?$"i, f), readdir(d))
        if !isempty(tifs)
            println("RAW : TIFF converti relu depuis le cache")
            chans, _ = read_tiff(joinpath(d, first(tifs)))
            return raw_orient(chans), 16
        end
    end
    conv = raw_converter()
    tmp = mktempdir(; prefix="raw_reg_")
    ok = false
    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)))
        ok = true
        return raw_orient(chans), 16
    finally
        if ok && RAW_CACHE
            RAW_CACHE_DIRS[path] = tmp
        else
            try
                rm(tmp; recursive=true, force=true)     # Windows : peut échouer si un fichier est encore ouvert
            catch
            end
        end
    end
end

# effacement des TIFF temporaires gardés par RAW_CACHE
function clear_raw_cache()
    for d in values(RAW_CACHE_DIRS)
        try
            rm(d; recursive=true, force=true)
        catch
        end
    end
    empty!(RAW_CACHE_DIRS)
    return nothing
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

#############################################################################################
# plan de fond de ciel v = a + gx * dx + gy * dy, ajusté par moindres carrés avec rejet itératif
# à 3 sigma (étoiles voisines, pixels chauds). Sortie : (a, gx, gy, sigma)
#############################################################################################
function fit_plane(dx::Vector{Float64}, dy::Vector{Float64}, v::Vector{Float64})
    a, s = robust_background(v)
    gx = 0.0
    gy = 0.0
    keep = abs.(v .- a) .< 3.0 * s            # départ : fond constant, étoiles voisines retirées
    for _ in 1:8
        n = count(keep)
        n >= 30 || break
        M = hcat(ones(n), dx[keep], dy[keep])
        a, gx, gy = M \ v[keep]
        r = v .- (a .+ gx .* dx .+ gy .* dy)
        rk = r[keep]
        s = 1.4826 * median(abs.(rk .- median(rk)))
        s > 0 || break
        newkeep = abs.(r) .< 3.0 * s
        newkeep == keep && break
        keep = newkeep
    end
    return a, gx, gy, max(s, 1e-6)
end

#############################################################################################
# bruit pixel à pixel d'une image : MAD des différences entre pixels voisins (divisée par √2).
# Insensible aux gradients (halo, Lune) : sert à pondérer les images dans la somme.
#############################################################################################
function noise_sigma(img::Matrix{Float32})
    ny, nx = size(img)
    stp = max(1, ny ÷ 400)
    d = Float64[]
    sizehint!(d, (ny ÷ stp + 1) * (nx ÷ 2 + 1))
    @inbounds for y in 1:stp:ny, x in 1:2:nx-1
        a = img[y, x]
        b = img[y, x+1]
        (a >= SAT_LEVEL || b >= SAT_LEVEL) && continue      # croissant saturé
        (a <= 0 && b <= 0) && continue                      # bords noirs éventuels
        push!(d, b - a)
    end
    length(d) >= 100 || return NaN
    m = median(d)
    return 1.4826 * median(abs.(d .- m)) / sqrt(2.0)
end

# luminance L = wr R + wv V + wb B (image couleur) ; image N&B renvoyée telle quelle
function luminance(chans::Vector{Matrix{Float32}})
    length(chans) == 1 && return chans[1]
    wr, wv, wb = Float32.(LUMA_WEIGHTS)
    return wr .* chans[1] .+ wv .* chans[2] .+ wb .* chans[3]
end

##################################################################################################################
# MESURE DE LA POSITION DE L'ÉTOILE
##################################################################################################################

#############################################################################################
# moyenne glissante 3x3 (sur les bords : moyenne des voisins disponibles)
#############################################################################################
function box_mean3(a::AbstractMatrix{Float64})
    h, w = size(a)
    out = similar(a)
    @inbounds for j in 1:w, i in 1:h
        s = 0.0
        n = 0
        for dj in -1:1, di in -1:1
            ii = i + di
            jj = j + dj
            if 1 <= ii <= h && 1 <= jj <= w
                s += a[ii, jj]
                n += 1
            end
        end
        out[i, j] = s / n
    end
    return out
end

#############################################################################################
# composante connexe (8-connexité) du masque contenant le pixel (i0, j0)
#############################################################################################
function connected_component(mask::AbstractMatrix{Bool}, i0::Int, j0::Int)
    h, w = size(mask)
    comp = falses(h, w)
    mask[i0, j0] || return comp
    comp[i0, j0] = true
    stack = [(i0, j0)]
    while !isempty(stack)
        i, j = pop!(stack)
        for dj in -1:1, di in -1:1
            ii = i + di
            jj = j + dj
            if 1 <= ii <= h && 1 <= jj <= w && mask[ii, jj] && !comp[ii, jj]
                comp[ii, jj] = true
                push!(stack, (ii, jj))
            end
        end
    end
    return comp
end

#############################################################################################
# maxima locaux de l'image lissée 3x3 dans un carré de demi-côté half centré sur (x0, y0),
# au-dessus de fond + MIN_PEAK_SNR sigma. Sortie : [(valeur, x, y)] du plus brillant au plus faible
# (maxima distants d'au moins 3 pixels : un plateau saturé ne donne qu'un candidat)
#############################################################################################
function find_peaks(img::Matrix{Float32}, x0::Real, y0::Real, half::Int)
    ny, nx = size(img)
    xa = max(1, round(Int, x0) - half)
    xb = min(nx, round(Int, x0) + half)
    ya = max(1, round(Int, y0) - half)
    yb = min(ny, round(Int, y0) + half)
    (xb - xa >= 4 && yb - ya >= 4) || return Tuple{Float64, Int, Int}[]
    cut = Float64.(img[ya:yb, xa:xb])
    sm = box_mean3(cut)
    bg, sd = robust_background(vec(cut))
    thr = bg + MIN_PEAK_SNR * sd / 3.0           # le lissage 3x3 divise le bruit par ~3
    h, w = size(sm)
    raw = Tuple{Float64, Int, Int}[]
    @inbounds for j in 2:w-1, i in 2:h-1
        v = sm[i, j]
        v > thr || continue
        ismax = true
        for dj in -1:1, di in -1:1
            sm[i+di, j+dj] > v && (ismax = false)
        end
        ismax && push!(raw, (v, xa + j - 1, ya + i - 1))
    end
    sort!(raw; by=first, rev=true)
    peaks = Tuple{Float64, Int, Int}[]
    for p in raw
        any(q -> hypot(q[2] - p[2], q[3] - p[3]) < 3.0, peaks) && continue
        push!(peaks, p)
        length(peaks) >= MAX_CANDIDATES && break
    end
    return peaks
end

#############################################################################################
# mesure précise de l'étoile dont le pic (image lissée) est au pixel (xp, yp).
#   - fond local : plan ajusté dans un anneau autour de l'étoile (itéré avec la délimitation)
#   - étoile : composante connexe des pixels au-dessus de DETECTION_SIGMA (aucune hypothèse de forme)
#   - position : barycentre, puis barycentre à fenêtre gaussienne itératif (XWIN)
#   - incertitude : propagation du bruit du fond (le bruit de photons de l'étoile n'est pas compté)
# Sortie : NamedTuple, ou nothing si l'objet n'est pas une étoile mesurable (trop faible, étendu...)
#############################################################################################
function refine_star(img::Matrix{Float32}, xp::Int, yp::Int)
    ny, nx = size(img)
    half = STAR_BOX_HALF
    xa = max(1, xp - half); xb = min(nx, xp + half)
    ya = max(1, yp - half); yb = min(ny, yp + half)
    cut = Float64.(img[ya:yb, xa:xb])
    h, w = size(cut)
    X(j) = Float64(xa + j - 1)                 # coordonnées image d'une colonne / ligne de la zone
    Y(i) = Float64(ya + i - 1)
    ip = yp - ya + 1
    jp = xp - xa + 1
    sm = box_mean3(cut)

    # fond initial : toute la zone (le sigma-clipping élimine l'étoile)
    a, sd = robust_background(vec(cut))
    gx = 0.0
    gy = 0.0
    bgf(i, j) = a + gx * (X(j) - xp) + gy * (Y(i) - yp)
    (sm[ip, jp] - bgf(ip, jp)) / (sd / 3.0) >= MIN_PEAK_SNR || return nothing

    xc = Float64(xp)
    yc = Float64(yp)
    comp = falses(h, w)
    r_ap = 0.0
    dmax = 0.0
    ring_short = false
    for iter in 1:4
        thr = DETECTION_SIGMA * sd / 3.0
        mask = falses(h, w)
        for j in 1:w, i in 1:h
            mask[i, j] = (sm[i, j] - bgf(i, j) > thr) && (hypot(j - jp, i - ip) <= MAX_STAR_RADIUS)
        end
        comp = connected_component(mask, ip, jp)
        any(comp) || return nothing

        sw = 0.0; sx = 0.0; sy = 0.0
        for j in 1:w, i in 1:h
            comp[i, j] || continue
            v = cut[i, j] - bgf(i, j)
            v > 0 || continue
            sw += v; sx += v * X(j); sy += v * Y(i)
        end
        sw > 0 || return nothing
        xc = sx / sw
        yc = sy / sw

        dmax = 0.0
        for j in 1:w, i in 1:h
            comp[i, j] && (dmax = max(dmax, hypot(X(j) - xc, Y(i) - yc)))
        end
        r_ap = 1.25 * dmax + 2.0
        r_in = r_ap + 3.0
        r_out = r_in + max(8.0, 0.5 * r_ap)
        ring_short = r_out > half

        # anneau de fond de ciel autour de l'étoile
        dxs = Float64[]; dys = Float64[]; vs = Float64[]
        for j in 1:w, i in 1:h
            r = hypot(X(j) - xc, Y(i) - yc)
            if r_in < r <= r_out
                push!(dxs, X(j) - xp); push!(dys, Y(i) - yp); push!(vs, cut[i, j])
            end
        end
        length(vs) >= 30 || break
        if BG_PLANE
            a, gx, gy, sd = fit_plane(dxs, dys, vs)
        else
            a, sd = robust_background(vs)
        end
    end
    # objet étendu (bord de la Lune, nébulosité...) : la composante atteint le rayon maximal
    dmax >= MAX_STAR_RADIUS - 1.5 && return nothing

    # forme : moments d'ordre 2 => largeur gaussienne équivalente (fenêtre du barycentre XWIN)
    sw = 0.0; mxx = 0.0; myy = 0.0; peak = -Inf; n_star = 0
    ex2 = 0.0; ey2 = 0.0
    for j in 1:w, i in 1:h
        comp[i, j] || continue
        n_star += 1
        peak = max(peak, cut[i, j])
        v = cut[i, j] - bgf(i, j)
        ex2 += (X(j) - xc)^2
        ey2 += (Y(i) - yc)^2
        v > 0 || continue
        sw += v; mxx += v * (X(j) - xc)^2; myy += v * (Y(i) - yc)^2
    end
    sw > 0 || return nothing
    s_gauss = sqrt(max((mxx + myy) / (2sw), 0.25))
    # incertitude du barycentre simple (bruit du fond)
    err_x = sd * sqrt(ex2) / sw
    err_y = sd * sqrt(ey2) / sw
    method = "barycentre"

    # affinage : barycentre à fenêtre gaussienne itératif (XWIN)
    if WINDOWED_CENTROID
        σw = s_gauss
        rmax2 = (4σw)^2
        xw = xc
        yw = yc
        conv = false
        for _ in 1:50
            swv = 0.0; sxv = 0.0; syv = 0.0
            for j in 1:w, i in 1:h
                ddx = X(j) - xw
                ddy = Y(i) - yw
                r2 = ddx^2 + ddy^2
                r2 <= rmax2 || continue
                wt = exp(-r2 / (2σw^2))
                v = cut[i, j] - bgf(i, j)
                swv += wt * v; sxv += wt * v * ddx; syv += wt * v * ddy
            end
            swv > 0 || break
            stx = 2.0 * sxv / swv
            sty = 2.0 * syv / swv
            xw += stx
            yw += sty
            hypot(xw - xc, yw - yc) > 3.0 && break           # divergence : on garde le barycentre simple
            if hypot(stx, sty) < 1e-4
                conv = true
                break
            end
        end
        if conv
            swv = 0.0; ex = 0.0; ey = 0.0
            for j in 1:w, i in 1:h
                ddx = X(j) - xw
                ddy = Y(i) - yw
                r2 = ddx^2 + ddy^2
                r2 <= rmax2 || continue
                wt = exp(-r2 / (2σw^2))
                swv += wt * (cut[i, j] - bgf(i, j))
                ex += wt^2 * ddx^2
                ey += wt^2 * ddy^2
            end
            if swv > 0
                xc = xw
                yc = yw
                err_x = 2.0 * sd * sqrt(ex) / swv
                err_y = 2.0 * sd * sqrt(ey) / swv
                method = "XWIN"
            end
        end
    end

    # flux dans l'ouverture (fond soustrait)
    flux = 0.0
    n_ap = 0
    for j in 1:w, i in 1:h
        hypot(X(j) - xc, Y(i) - yc) <= r_ap || continue
        flux += cut[i, j] - bgf(i, j)
        n_ap += 1
    end
    flux > 0 || return nothing

    return (x=xc, y=yc, ex=err_x, ey=err_y, err=hypot(err_x, err_y),
            flux=flux, snr=flux / (sd * sqrt(n_ap)), peak=peak, bg=a, sd=sd,
            fwhm=2.3548 * s_gauss, r_ap=r_ap, n_star=n_star,
            sat=peak >= SAT_LEVEL, ring_short=ring_short, method=method)
end

#############################################################################################
# recherche et mesure de l'étoile dans le carré de demi-côté half centré sur (x0, y0).
#   ref_flux = nothing : étoile la plus brillante (mesurable) du carré
#   ref_flux = flux    : parmi les étoiles de flux compatible (FLUX_TOLERANCE), la plus proche de
#                        (x0, y0), position attendue
#############################################################################################
function locate_star(img::Matrix{Float32}, x0::Real, y0::Real, half::Int, ref_flux)
    best = nothing
    dbest = Inf
    for (_, xp, yp) in find_peaks(img, x0, y0, half)
        r = refine_star(img, xp, yp)
        r === nothing && continue
        ref_flux === nothing && return r
        ratio = r.flux / ref_flux
        (1 / FLUX_TOLERANCE <= ratio <= FLUX_TOLERANCE) || continue
        d = hypot(r.x - x0, r.y - y0)
        if d < dbest
            best = r
            dbest = d
        end
    end
    return best
end

##################################################################################################################
# PRÉ-RECENTRAGE GROSSIER SUR LA LUNE
##################################################################################################################

# image réduite par blocs b x b (moyenne), y = 1 en haut (les dernières lignes / colonnes incomplètes sont ignorées)
function bin_image(img::Matrix{Float32}, b::Int)
    ny, nx = size(img)
    my_ = ny ÷ b
    mx_ = nx ÷ b
    sm = Matrix{Float32}(undef, my_, mx_)
    winv = 1f0 / Float32(b * b)
    Threads.@threads for J in 1:mx_
        @inbounds for I in 1:my_
            s = 0f0
            for jj in (J-1)*b+1:J*b, ii in (I-1)*b+1:I*b
                s += img[ii, jj]
            end
            sm[I, J] = s * winv
        end
    end
    return sm
end

# niveau du plateau saturé du croissant : maximum de l'image réduite MOON_BIN x MOON_BIN (pixels chauds atténués)
plateau_level(img::Matrix{Float32}) = Float64(maximum(bin_image(img, max(1, MOON_BIN))))

#############################################################################################
# centre du croissant lunaire : image réduite par blocs MOON_BIN x MOON_BIN (moyenne), seuil
# fond + MOON_LEVEL x (maximum - fond), plus grande composante connexe (8-connexité, les étoiles et
# pixels chauds isolés sont écartés), barycentre géométrique de cette zone (le croissant étant saturé,
# un barycentre non pondéré est le plus stable d'une image à l'autre).
# Sortie : (x, y, aire en pixels d'origine, edge = zone touchant le bord du champ, box = (xa, xb, ya, yb)
#          rectangle englobant la zone (avec une marge), top = niveau du plateau), ou nothing
#############################################################################################
function moon_center(img::Matrix{Float32})
    ny, nx = size(img)
    b = max(1, MOON_BIN)
    (ny ÷ b >= 8 && nx ÷ b >= 8) || return nothing
    sm = bin_image(img, b)
    my_, mx_ = size(sm)
    stp = max(1, round(Int, sqrt(length(sm) / 1.0e6)))
    bg = Float64(median(vec(sm[1:stp:end, 1:stp:end])))
    top = Float64(maximum(sm))
    top > bg || return nothing
    thr = Float32(bg + MOON_LEVEL * (top - bg))

    # plus grande composante connexe des pixels au-dessus du seuil
    lab = zeros(Int32, my_, mx_)
    best = Int32(0)
    bestn = 0
    nl = Int32(0)
    stack = Tuple{Int, Int}[]
    @inbounds for J in 1:mx_, I in 1:my_
        (sm[I, J] > thr && lab[I, J] == 0) || continue
        nl += Int32(1)
        lab[I, J] = nl
        push!(stack, (I, J))
        cnt = 0
        while !isempty(stack)
            i, j = pop!(stack)
            cnt += 1
            for dj in -1:1, di in -1:1
                ii = i + di
                jj = j + dj
                if 1 <= ii <= my_ && 1 <= jj <= mx_ && lab[ii, jj] == 0 && sm[ii, jj] > thr
                    lab[ii, jj] = nl
                    push!(stack, (ii, jj))
                end
            end
        end
        if cnt > bestn
            best = nl
            bestn = cnt
        end
    end
    bestn * b * b >= MOON_MIN_AREA || return nothing

    sx = 0.0
    sy = 0.0
    edge = false
    Imin, Imax, Jmin, Jmax = my_, 1, mx_, 1
    @inbounds for J in 1:mx_, I in 1:my_
        lab[I, J] == best || continue
        sx += J
        sy += I
        Imin = min(Imin, I); Imax = max(Imax, I)
        Jmin = min(Jmin, J); Jmax = max(Jmax, J)
        (I == 1 || I == my_ || J == 1 || J == mx_) && (edge = true)
    end
    # le bloc J couvre les colonnes (J-1) b + 1 ... J b : centre (J-1) b + (b+1)/2
    xm = (sx / bestn - 1) * b + (b + 1) / 2
    ym = (sy / bestn - 1) * b + (b + 1) / 2
    mg = 2b + 10
    box = (max(1, (Jmin - 1) * b + 1 - mg), min(nx, Jmax * b + mg),
           max(1, (Imin - 1) * b + 1 - mg), min(ny, Imax * b + mg))
    return (x=xm, y=ym, area=bestn * b * b, edge=edge, box=box, top=top)
end

# mesure de la Lune dans une image (couche V si DISPERSION, luminance sinon) ; résultat dans rec
function measure_moon!(rec, chans::Vector{Matrix{Float32}})
    det = DISPERSION ? chans[2] : luminance(chans)
    m = moon_center(det)
    if m === nothing
        rec.okm = false
        return rec
    end
    rec.mx = m.x
    rec.my = m.y
    rec.okm = true
    rec.medge = m.edge
    rec.mbox = m.box
    rec.mtop = m.top
    return rec
end

##################################################################################################################
# RECENTRAGE FIN SUR LA LUNE (version 2.0)
##################################################################################################################

const MOON_METHODS = (:correlation, :limb, :crescent)
method_name(m::Symbol) = m == :correlation ? "corrélation" : m == :limb ? "limbe" : "croissant"

# données de référence (1re image) du recentrage fin sur la Lune
Base.@kwdef mutable struct MoonRef
    methods::Vector{Symbol} = Symbol[]           # méthodes calculées
    mx::Float64 = NaN                            # centre grossier du croissant
    my::Float64 = NaN
    xr::UnitRange{Int} = 1:0                     # fenêtre de corrélation (colonnes, lignes)
    yr::UnitRange{Int} = 1:0
    T::Matrix{Float64} = zeros(0, 0)             # fenêtre de référence lissée
    tbad::BitMatrix = falses(0, 0)               # pixels exclus (saturés)
    cx::Float64 = NaN                            # cercle du limbe
    cy::Float64 = NaN
    R::Float64 = NaN
end

#############################################################################################
# lissage gaussien séparable (bords : poids renormalisés)
#############################################################################################
function gauss_smooth(A::Matrix{Float64}, σ::Float64)
    σ > 0 || return A
    r = max(1, ceil(Int, 3σ))
    k = [exp(-t^2 / (2σ^2)) for t in -r:r]
    h, w = size(A)
    Tm = similar(A)
    @inbounds for j in 1:w, i in 1:h
        s = 0.0
        sw = 0.0
        for t in -r:r
            jj = j + t
            (1 <= jj <= w) || continue
            s += k[t+r+1] * A[i, jj]
            sw += k[t+r+1]
        end
        Tm[i, j] = s / sw
    end
    B = similar(A)
    @inbounds for j in 1:w, i in 1:h
        s = 0.0
        sw = 0.0
        for t in -r:r
            ii = i + t
            (1 <= ii <= h) || continue
            s += k[t+r+1] * Tm[ii, j]
            sw += k[t+r+1]
        end
        B[i, j] = s / sw
    end
    return B
end

# dilatation d'un masque par un carré de demi-côté r
function dilate_mask(M::BitMatrix, r::Int)
    r <= 0 && return copy(M)
    h, w = size(M)
    Tm = falses(h, w)
    @inbounds for j in 1:w, i in 1:h
        M[i, j] || continue
        for jj in max(1, j - r):min(w, j + r)
            Tm[i, jj] = true
        end
    end
    D = falses(h, w)
    @inbounds for j in 1:w, i in 1:h
        Tm[i, j] || continue
        for ii in max(1, i - r):min(h, i + r)
            D[ii, j] = true
        end
    end
    return D
end

#############################################################################################
# passe-bas "normalisé" par 3 flous de boîte successifs (≈ gaussienne d'écart-type σ), pixels exclus (bad)
# ignorés : fond = flou(A x valide) / flou(valide) ; les bords de la fenêtre sont traités de la même façon
#############################################################################################
function box_pass(A::Matrix{Float64}, r::Int, dim::Int)
    h, w = size(A)
    B = similar(A)
    if dim == 1
        @inbounds for j in 1:w
            s = 0.0
            for i in 1:min(h, 1 + r)
                s += A[i, j]
            end
            for i in 1:h
                B[i, j] = s
                i + r + 1 <= h && (s += A[i+r+1, j])
                i - r >= 1 && (s -= A[i-r, j])
            end
        end
    else
        @inbounds for i in 1:h
            s = 0.0
            for j in 1:min(w, 1 + r)
                s += A[i, j]
            end
            for j in 1:w
                B[i, j] = s
                j + r + 1 <= w && (s += A[i, j+r+1])
                j - r >= 1 && (s -= A[i, j-r])
            end
        end
    end
    return B
end

function masked_lowpass(A::Matrix{Float64}, bad::BitMatrix, σ::Float64)
    r = max(1, round(Int, (sqrt(1 + 4σ^2) - 1) / 2))      # 3 boîtes de demi-largeur r : σ² = r (r + 1)
    num = similar(A)
    den = similar(A)
    @inbounds for i in eachindex(A)
        num[i] = bad[i] ? 0.0 : A[i]
        den[i] = bad[i] ? 0.0 : 1.0
    end
    for _ in 1:3, d in 1:2
        num = box_pass(num, r, d)
        den = box_pass(den, r, d)
    end
    @inbounds for i in eachindex(num)
        num[i] = den[i] > 1e-9 ? num[i] / den[i] : 0.0
    end
    return num
end

# passe-haut MOON_CORR_HIGHPASS (fenêtre moins son fond lissé)
highpass(A::Matrix{Float64}, bad::BitMatrix) =
    MOON_CORR_HIGHPASS > 0 ? A .- masked_lowpass(A, bad, MOON_CORR_HIGHPASS) : A

# demi-largeur de la dilatation des pixels saturés (le lissage étale leur valeur)
sat_dilation() = MOON_CORR_SMOOTH > 0 ? ceil(Int, 3 * MOON_CORR_SMOOTH) + 1 : 1

#############################################################################################
# valeurs de img aux points (x + sx, y + sy), x ∈ xr, y ∈ yr : interpolation bicubique séparable (même
# noyau que la translation des images). Sortie : (fenêtre Float64, masque des pixels dont le voisinage
# d'interpolation contient un pixel saturé (>= sat)), ou nothing si la fenêtre décalée sort de l'image
#############################################################################################
function sample_window(img::Matrix{Float32}, xr::UnitRange{Int}, yr::UnitRange{Int},
                       sx::Float64, sy::Float64, sat::Float64)
    ny, nx = size(img)
    ox, wx = shift_taps(-sx, :bicubic)        # sortie(x) = Σ w img(x + o) ≈ img(x + sx)
    oy, wy = shift_taps(-sy, :bicubic)
    (first(xr) + minimum(ox) >= 1 && last(xr) + maximum(ox) <= nx &&
     first(yr) + minimum(oy) >= 1 && last(yr) + maximum(oy) <= ny) || return nothing
    r0 = first(yr) + minimum(oy)
    r1 = last(yr) + maximum(oy)
    h = length(yr)
    w = length(xr)
    nrow = r1 - r0 + 1
    tmp = Matrix{Float64}(undef, nrow, w)
    stmp = falses(nrow, w)
    satf = Float32(sat)
    @inbounds for j in 1:w
        x = xr[j]
        for row in r0:r1
            s = 0.0
            sat_any = false
            for k in eachindex(ox)
                v = img[row, x + ox[k]]
                s += wx[k] * v
                v >= satf && (sat_any = true)
            end
            tmp[row - r0 + 1, j] = s
            stmp[row - r0 + 1, j] = sat_any
        end
    end
    out = Matrix{Float64}(undef, h, w)
    m = falses(h, w)
    @inbounds for j in 1:w, i in 1:h
        base = yr[i] - r0 + 1
        s = 0.0
        sat_any = false
        for k in eachindex(oy)
            s += wy[k] * tmp[base + oy[k], j]
            stmp[base + oy[k], j] && (sat_any = true)
        end
        out[i, j] = s
        m[i, j] = sat_any
    end
    return out, m
end

#############################################################################################
# fenêtre de corrélation de référence (1re image), MOON_WIN x MOON_WIN centrée sur (xw, yw), ramenée
# dans l'image avec une marge. Sortie : (xr, yr, T lissée, tbad, frac = fraction de pixels utilisables)
#############################################################################################
function moon_template(det::Matrix{Float32}, xw::Float64, yw::Float64, sat::Float64)
    ny, nx = size(det)
    margin = MOON_CORR_SEARCH + 8
    side = min(MOON_WIN, nx - 2margin, ny - 2margin)
    side >= 50 || return nothing
    xa = clamp(round(Int, xw) - side ÷ 2, 1 + margin, nx - margin - side + 1)
    ya = clamp(round(Int, yw) - side ÷ 2, 1 + margin, ny - margin - side + 1)
    xr = xa:xa+side-1
    yr = ya:ya+side-1
    T0 = Float64.(det[yr, xr])
    tbad = dilate_mask(BitMatrix(T0 .>= sat), sat_dilation())
    frac = 1 - count(tbad) / length(tbad)
    return (xr=xr, yr=yr, T=highpass(gauss_smooth(T0, MOON_CORR_SMOOTH), tbad), tbad=tbad, frac=frac)
end

#############################################################################################
# recalage fin d'une fenêtre par moindres carrés (Lucas-Kanade, Gauss-Newton) : recherche du
# déplacement s = (sx, sy) tel que
#        img(x + s) ≈ g T(x) + o + px X + py Y        (X, Y : coordonnées réduites -1..1 de la fenêtre)
#   T : fenêtre de référence lissée et filtrée (passe-haut MOON_CORR_HIGHPASS), tbad : pixels de T exclus
#   (saturés) ; g absorbe les variations de transparence ; le passe-haut et le plan (o, px, py) celles du
#   halo de lumière diffusée et du fond de ciel. La fenêtre de img subit les mêmes traitements.
#   1) si search > 0 : recherche entière (corrélation normalisée) à ± search pixels autour de s0
#   2) affinage sous-pixel (img interpolée en bicubique et lissée comme T)
#   3) rejet des pixels à plus de 5 sigma (étoiles, pixels chauds, rayons cosmiques), nouvel affinage
# Sortie : (sx, sy, ex, ey, g, frac, W = fenêtre img(x + s) lissée, bad = pixels exclus), ou nothing
# ex, ey : incertitudes estimées (bruit de l'image et de la référence, corrélation due au lissage)
#############################################################################################
function lk_register(T::Matrix{Float64}, tbad::BitMatrix, img::Matrix{Float32}, xr::UnitRange{Int},
                     yr::UnitRange{Int}, s0x::Float64, s0y::Float64, sat::Float64, search::Int)
    h, w = size(T)
    σs = MOON_CORR_SMOOTH
    rd = sat_dilation()

    # 1) recherche entière
    sx = s0x
    sy = s0y
    if search > 0
        ix0 = round(Int, s0x)
        iy0 = round(Int, s0y)
        best = -Inf
        bx, by = ix0, iy0
        for dy in -search:search, dx in -search:search
            smp = sample_window(img, xr, yr, Float64(ix0 + dx), Float64(iy0 + dy), sat)
            smp === nothing && continue
            Wi, wsat = smp
            Wi = highpass(Wi, wsat .| tbad)
            n = 0
            sa = 0.0
            sb = 0.0
            @inbounds for j in 1:w, i in 1:h
                (tbad[i, j] || wsat[i, j]) && continue
                n += 1
                sa += T[i, j]
                sb += Wi[i, j]
            end
            n >= 100 || continue
            ma = sa / n
            mb = sb / n
            sab = 0.0
            saa = 0.0
            sbb = 0.0
            @inbounds for j in 1:w, i in 1:h
                (tbad[i, j] || wsat[i, j]) && continue
                a = T[i, j] - ma
                b = Wi[i, j] - mb
                sab += a * b
                saa += a * a
                sbb += b * b
            end
            c = sab / sqrt(saa * sbb + 1e-30)
            if c > best
                best = c
                bx = ix0 + dx
                by = iy0 + dy
            end
        end
        isfinite(best) || return nothing
        (abs(bx - ix0) == search || abs(by - iy0) == search) &&
            println("      corrélation : maximum en bordure de la zone de recherche (augmenter MOON_CORR_SEARCH ?)")
        sx = Float64(bx)
        sy = Float64(by)
    end

    # 2) et 3) affinage sous-pixel
    Xc = [(j - (w + 1) / 2) / (w / 2) for j in 1:w]
    Yc = [(i - (h + 1) / 2) / (h / 2) for i in 1:h]
    p = [1.0, 0.0, 0.0, 0.0]                  # g, o, px, py
    rej = falses(h, w)
    rejected = false
    conv = false
    W = zeros(0, 0)
    bad = falses(h, w)
    A = zeros(6, 6)
    sig = NaN
    nuse = 0
    for it in 1:60
        smp = sample_window(img, xr, yr, sx, sy, sat)
        smp === nothing && return nothing
        bad = dilate_mask(smp[2], rd) .| tbad .| rej
        W = highpass(gauss_smooth(smp[1], σs), bad)
        fill!(A, 0.0)
        bv = zeros(6)
        ss = 0.0
        nuse = 0
        @inbounds for j in 3:w-2, i in 3:h-2
            bad[i, j] && continue
            gx = 0.5 * (W[i, j+1] - W[i, j-1])
            gy = 0.5 * (W[i+1, j] - W[i-1, j])
            r = W[i, j] - (p[1] * T[i, j] + p[2] + p[3] * Xc[j] + p[4] * Yc[i])
            J = (gx, gy, -T[i, j], -1.0, -Xc[j], -Yc[i])
            for a in 1:6
                bv[a] += J[a] * r
                for b in a:6
                    A[a, b] += J[a] * J[b]
                end
            end
            ss += r * r
            nuse += 1
        end
        nuse >= 100 || return nothing
        for a in 1:6, b in 1:a-1
            A[a, b] = A[b, a]
        end
        Δ = try
            -(A \ bv)
        catch
            return nothing
        end
        all(isfinite, Δ) || return nothing
        sig = sqrt(ss / max(nuse - 6, 1))
        sx += Δ[1]
        sy += Δ[2]
        p .+= Δ[3:6]
        hypot(sx - s0x, sy - s0y) > max(search, 2) + 5 && return nothing      # divergence
        if hypot(Δ[1], Δ[2]) < 1e-4
            if rejected
                conv = true
                break
            end
            # rejet des pixels aberrants, puis nouvel affinage
            @inbounds for j in 3:w-2, i in 3:h-2
                bad[i, j] && continue
                r = W[i, j] - (p[1] * T[i, j] + p[2] + p[3] * Xc[j] + p[4] * Yc[i])
                abs(r) > 5sig && (rej[i, j] = true)
            end
            rejected = true
        end
    end
    conv || return nothing
    C = try
        inv(A)
    catch
        return nothing
    end
    # bruit corrélé par le lissage (aire effective 4π σ²) et bruit de la référence (facteur √2)
    f = (σs > 0 ? sqrt(4π) * σs : 1.0) * sqrt(2.0)
    ex = f * sig * sqrt(max(C[1, 1], 0.0))
    ey = f * sig * sqrt(max(C[2, 2], 0.0))
    return (sx=sx, sy=sy, ex=ex, ey=ey, g=p[1], frac=nuse / ((h - 4) * (w - 4)), W=W, bad=bad)
end

# interpolation bilinéaire ; NaN hors de l'image
@inline function bilinear(img::Matrix{Float32}, x::Float64, y::Float64)
    ny, nx = size(img)
    (x >= 1 && y >= 1 && x < nx && y < ny) || return NaN
    j = floor(Int, x)
    i = floor(Int, y)
    fx = x - j
    fy = y - i
    @inbounds return (1 - fy) * ((1 - fx) * img[i, j] + fx * img[i, j+1]) +
                     fy * ((1 - fx) * img[i+1, j] + fx * img[i+1, j+1])
end

#############################################################################################
# points du limbe le long de MOON_LIMB_RAYS rayons issus de (xc, yc), autour du rayon R :
#   profil radial (moyenne sur 5 pixels perpendiculairement au rayon, pas de 0,25 px), dérivée lissée
#   (gaussienne, σ = 1 px) ; bord = chute d'intensité la plus forte vers l'extérieur à ± MOON_LIMB_RANGE
#   du cercle (le terminateur, montée d'intensité vers l'extérieur, n'est pas retenu) ; position
#   sous-pixel par parabole. Bord gardé si le saut d'intensité (moyennes de 3 à 8 px de part et d'autre)
#   dépasse MOON_LIMB_MIN_SNR sigma ; "brillant" si l'intérieur est saturé (>= sat).
# Sortie : (xs, ys, bright)
#############################################################################################
function limb_points(img::Matrix{Float32}, xc::Float64, yc::Float64, R::Float64, sat::Float64, sig::Float64)
    stp = 0.25
    rng = Float64(MOON_LIMB_RANGE)
    band = 8.0
    nr = round(Int, 2 * (rng + band + 2) / stp) + 1
    r0 = R - rng - band - 2
    rs = [r0 + (k - 1) * stp for k in 1:nr]
    gs = 4                                         # σ de la dérivée gaussienne : 1 px = 4 échantillons
    kt = collect(-3gs:3gs)
    dk = [t * exp(-t^2 / (2gs^2)) for t in kt]
    dk ./= sum(kt[q] * dk[q] for q in eachindex(kt)) * stp      # dérivée en ADU / pixel
    klo = findfirst(r -> r >= R - rng, rs)
    khi = findlast(r -> r <= R + rng, rs)
    prof = zeros(nr)
    der = zeros(nr)
    thr = MOON_LIMB_MIN_SNR * sig * sqrt(2 / 25)    # ~25 mesures indépendantes par bande
    xs = Float64[]
    ys = Float64[]
    br = Bool[]
    for a in 0:MOON_LIMB_RAYS-1
        θ = 2π * a / MOON_LIMB_RAYS
        c = cos(θ)
        s = sin(θ)
        ok = true
        for k in 1:nr
            acc = 0.0
            for t in -2:2
                v = bilinear(img, xc + rs[k] * c - t * s, yc + rs[k] * s + t * c)
                if isnan(v)
                    ok = false
                    break
                end
                acc += v
            end
            ok || break
            prof[k] = acc / 5
        end
        ok || continue
        best = Inf
        kb = 0
        for k in klo-1:khi+1
            d = 0.0
            for q in eachindex(kt)
                d += dk[q] * prof[k + kt[q]]
            end
            der[k] = d
            if klo <= k <= khi && d < best
                best = d
                kb = k
            end
        end
        (best < 0 && klo < kb < khi) || continue
        y0 = der[kb-1]
        y1 = der[kb]
        y2 = der[kb+1]
        den = y0 - 2y1 + y2
        δ = den > 0 ? clamp(0.5 * (y0 - y2) / den, -0.5, 0.5) : 0.0
        re = rs[kb] + δ * stp
        ins = 0.0; ni = 0; outs = 0.0; no = 0; pmax = -Inf
        for k in 1:nr
            d = rs[k] - re
            if -band <= d <= -3
                ins += prof[k]; ni += 1; pmax = max(pmax, prof[k])
            elseif 3 <= d <= band
                outs += prof[k]; no += 1
            end
        end
        (ni > 0 && no > 0) || continue
        ins / ni - outs / no > thr || continue
        bright = pmax >= sat
        (MOON_LIMB_SIDE == :bright && !bright) && continue
        (MOON_LIMB_SIDE == :earthshine && bright) && continue
        push!(xs, xc + re * c)
        push!(ys, yc + re * s)
        push!(br, bright)
    end
    return xs, ys, br
end

#############################################################################################
# cercle ajusté sur des points : ajustement algébrique (Kasa), puis géométrique (Gauss-Newton sur les
# distances au cercle) avec rejet itératif des points à plus de 3 sigma.
# Sortie : (xc, yc, R, ex, ey, rms, n, keep) ou nothing
#############################################################################################
function fit_circle(xs::Vector{Float64}, ys::Vector{Float64})
    n = length(xs)
    n >= 10 || return nothing
    mx = mean(xs)
    my = mean(ys)
    u = xs .- mx
    v = ys .- my
    sol = try
        hcat(u, v, ones(n)) \ (-(u .^ 2 .+ v .^ 2))
    catch
        return nothing
    end
    a = -sol[1] / 2
    b = -sol[2] / 2
    R2 = a^2 + b^2 - sol[3]
    R2 > 0 || return nothing
    R = sqrt(R2)
    keep = trues(n)
    JTJ = zeros(3, 3)
    for pass in 1:10
        for it in 1:30
            fill!(JTJ, 0.0)
            JTr = zeros(3)
            for k in 1:n
                keep[k] || continue
                dx = u[k] - a
                dy = v[k] - b
                d = hypot(dx, dy)
                d > 0 || continue
                r = d - R
                J = (-dx / d, -dy / d, -1.0)
                for p in 1:3
                    JTr[p] += J[p] * r
                    for q in 1:3
                        JTJ[p, q] += J[p] * J[q]
                    end
                end
            end
            Δ = try
                -(JTJ \ JTr)
            catch
                return nothing
            end
            a += Δ[1]
            b += Δ[2]
            R += Δ[3]
            maximum(abs, Δ) < 1e-7 && break
        end
        res = [hypot(u[k] - a, v[k] - b) - R for k in 1:n]
        rk = res[keep]
        s = max(1.4826 * median(abs.(rk .- median(rk))), 0.02)
        newkeep = abs.(res) .<= 3s
        count(newkeep) >= 10 || return nothing
        newkeep == keep && break
        keep = newkeep
    end
    nk = count(keep)
    s2 = sum(abs2, [hypot(u[k] - a, v[k] - b) - R for k in 1:n if keep[k]]) / max(nk - 3, 1)
    C = try
        inv(JTJ)
    catch
        return nothing
    end
    return (xc=a + mx, yc=b + my, R=R, ex=sqrt(s2 * max(C[1, 1], 0.0)), ey=sqrt(s2 * max(C[2, 2], 0.0)),
            rms=sqrt(s2), n=nk, keep=keep)
end

#############################################################################################
# cercle du limbe à partir d'un cercle approché (xc, yc, R) : points du limbe, ajustement, nouveaux points
# autour du cercle ajusté... (3 passes au plus). Sortie : fit_circle + (nb, ne) = nombres de points retenus
# côté croissant (brillant) et côté lumière cendrée, ou nothing
#############################################################################################
function limb_circle(img::Matrix{Float32}, xc::Float64, yc::Float64, R::Float64, sat::Float64, sig::Float64)
    fit = nothing
    nb = 0
    ne = 0
    for pass in 1:3
        xs, ys, br = limb_points(img, xc, yc, R, sat, sig)
        f = fit_circle(xs, ys)
        f === nothing && return nothing
        moved = hypot(f.xc - xc, f.yc - yc) + abs(f.R - R)
        xc, yc, R = f.xc, f.yc, f.R
        fit = f
        nb = count(br .& f.keep)
        ne = count(.!br .& f.keep)
        moved < 0.02 && break
    end
    return merge(fit, (nb=nb, ne=ne))
end

#############################################################################################
# cercle approché du limbe sur l'image réduite MOON_BIN x MOON_BIN, sans connaissance préalable :
# pixels de fort gradient (Sobel), cercle ayant le plus de points à moins de 1,5 pixel (tirages
# aléatoires de 3 points, RANSAC). Le limbe (côté croissant + côté lumière cendrée) l'emporte sur le
# terminateur, qui n'est pas un cercle. Sortie : (xc, yc, R) en pixels de l'image d'origine, ou nothing
#############################################################################################
function moon_circle_guess(img::Matrix{Float32})
    b = max(1, MOON_BIN)
    sm = bin_image(img, b)
    h, w = size(sm)
    (h >= 8 && w >= 8) || return nothing
    G = zeros(Float64, h, w)
    @inbounds for j in 2:w-1, i in 2:h-1
        gx = (sm[i-1, j+1] + 2sm[i, j+1] + sm[i+1, j+1]) - (sm[i-1, j-1] + 2sm[i, j-1] + sm[i+1, j-1])
        gy = (sm[i+1, j-1] + 2sm[i+1, j] + sm[i+1, j+1]) - (sm[i-1, j-1] + 2sm[i-1, j] + sm[i-1, j+1])
        G[i, j] = hypot(gx, gy) / 8
    end
    stp = max(1, round(Int, sqrt(h * w / 2.5e5)))
    gb, gsd = robust_background(vec(G[2:stp:h-1, 2:stp:w-1]))
    thr = gb + 8gsd
    px = Float64[]
    py = Float64[]
    for j in 2:w-1, i in 2:h-1
        if G[i, j] > thr
            push!(px, j)
            push!(py, i)
        end
    end
    np = length(px)
    np >= 20 || return nothing
    rng = MersenneTwister(12345)
    if np > 4000
        idx = randperm(rng, np)[1:4000]
        px = px[idx]
        py = py[idx]
        np = 4000
    end
    Rmin = MOON_MIN_RADIUS / b
    Rmax = hypot(h, w)
    bestn = 0
    bc = (0.0, 0.0, 0.0)
    for it in 1:3000
        i1 = rand(rng, 1:np); i2 = rand(rng, 1:np); i3 = rand(rng, 1:np)
        (i1 == i2 || i2 == i3 || i1 == i3) && continue
        x1, y1 = px[i1], py[i1]
        x2, y2 = px[i2], py[i2]
        x3, y3 = px[i3], py[i3]
        d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
        abs(d) < 1e-9 && continue
        q1 = x1^2 + y1^2
        q2 = x2^2 + y2^2
        q3 = x3^2 + y3^2
        ux = (q1 * (y2 - y3) + q2 * (y3 - y1) + q3 * (y1 - y2)) / d
        uy = (q1 * (x3 - x2) + q2 * (x1 - x3) + q3 * (x2 - x1)) / d
        R = hypot(x1 - ux, y1 - uy)
        (Rmin <= R <= Rmax) || continue
        cnt = 0
        @inbounds for k in 1:np
            abs(hypot(px[k] - ux, py[k] - uy) - R) < 1.5 && (cnt += 1)
        end
        if cnt > bestn
            bestn = cnt
            bc = (ux, uy, R)
        end
    end
    bestn >= 20 || return nothing
    inl = [k for k in 1:np if abs(hypot(px[k] - bc[1], py[k] - bc[2]) - bc[3]) < 1.5]
    f = fit_circle(px[inl], py[inl])
    xc, yc, R = f === nothing ? bc : (f.xc, f.yc, f.R)
    return (xc=(xc - 1) * b + (b + 1) / 2, yc=(yc - 1) * b + (b + 1) / 2, R=R * b)
end

#############################################################################################
# barycentre du croissant : poids 0 sous (MOON_CRESCENT_LEVEL - 0,1) x plateau, 1 au-dessus de
# (MOON_CRESCENT_LEVEL + 0,1) x plateau, rampe linéaire entre les deux, dans le rectangle box
#############################################################################################
function crescent_center(img::Matrix{Float32}, box::NTuple{4, Int}, top::Float64)
    xa, xb, ya, yb = box
    lo = (MOON_CRESCENT_LEVEL - 0.1) * top
    hi = (MOON_CRESCENT_LEVEL + 0.1) * top
    sw = 0.0
    sx = 0.0
    sy = 0.0
    @inbounds for x in xa:xb, y in ya:yb
        v = img[y, x]
        wgt = v <= lo ? 0.0 : v >= hi ? 1.0 : (v - lo) / (hi - lo)
        wgt > 0 || continue
        sw += wgt
        sx += wgt * x
        sy += wgt * y
    end
    sw > 0 || return nothing
    return (x=sx / sw, y=sy / sw, n=sw)
end

#############################################################################################
# position de la Lune dans un plan par une méthode. "Position" = (x, y, incertitude, qualité) :
#   :correlation  déplacement par rapport à la 1re image ; qualité = gain g (transparence relative)
#   :limb         centre du cercle du limbe ; qualité = rayon
#   :crescent     barycentre du croissant ; qualité = surface pondérée (pixels)
# d0x, d0y : déplacement grossier de la Lune par rapport à la 1re image. Sortie : (position, aux) ou nothing
#############################################################################################
function moon_measure(m::Symbol, img::Matrix{Float32}, rec, mref::MoonRef, d0x::Float64, d0y::Float64,
                      sat::Float64, top::Float64, sig::Float64)
    if m == :correlation
        lk = lk_register(mref.T, mref.tbad, img, mref.xr, mref.yr, d0x, d0y, sat, MOON_CORR_SEARCH)
        lk === nothing && return nothing
        return (lk.sx, lk.sy, hypot(lk.ex, lk.ey), lk.g), lk
    elseif m == :limb
        f = limb_circle(img, mref.cx + d0x, mref.cy + d0y, mref.R, sat, sig)
        f === nothing && return nothing
        return (f.xc, f.yc, hypot(f.ex, f.ey), f.R), f
    else
        c = crescent_center(img, rec.mbox, top)
        c === nothing && return nothing
        return (c.x, c.y, NaN, c.n), c
    end
end

#############################################################################################
# écart de position de la Lune entre la couche imgc et la couche de détection (même image), par la
# méthode m ; pos, aux : résultat de moon_measure sur la couche de détection. Sortie : (dx, dy) ou nothing
#############################################################################################
function channel_offset(m::Symbol, imgc::Matrix{Float32}, pos, aux, rec, mref::MoonRef)
    topc = plateau_level(imgc)
    satc = MOON_SAT_FRAC * topc
    if m == :correlation
        lk = lk_register(aux.W, aux.bad, imgc, mref.xr, mref.yr, aux.sx, aux.sy, satc, 0)
        lk === nothing && return nothing
        return (lk.sx - aux.sx, lk.sy - aux.sy)
    elseif m == :limb
        f = limb_circle(imgc, aux.xc, aux.yc, aux.R, satc, noise_sigma(imgc))
        f === nothing && return nothing
        return (f.xc - aux.xc, f.yc - aux.yc)
    else
        c = crescent_center(imgc, rec.mbox, topc)
        c === nothing && return nothing
        return (c.x - pos[1], c.y - pos[2])
    end
end

#############################################################################################
# mesure fine de la Lune dans une image : position par chacune des méthodes mref.methods, puis
# (DISPERSION) écarts R - V et B - V par la méthode MOON_METHOD.
# Prérequis : mesure grossière de la Lune (measure_moon!) et du bruit (rec.noise) sur cette image.
#############################################################################################
function measure_moon_fine!(rec, chans::Vector{Matrix{Float32}}, mref::MoonRef)
    if !rec.okm
        rec.mok = false
        rec.mstatus = "Lune non trouvée"
        return rec
    end
    det = DISPERSION ? chans[2] : luminance(chans)
    top = rec.mtop
    sat = MOON_SAT_FRAC * top
    sig = isfinite(rec.noise) ? rec.noise : noise_sigma(det)
    d0x = rec.mx - mref.mx
    d0y = rec.my - mref.my
    warn = String[]
    sel = nothing
    for m in mref.methods
        res = moon_measure(m, det, rec, mref, d0x, d0y, sat, top, sig)
        if res === nothing
            push!(warn, method_name(m) * " en échec")
            continue
        end
        rec.mpos[m] = res[1]
        m == MOON_METHOD && (sel = res)
    end
    if sel === nothing
        rec.mok = false
        rec.mstatus = "échec de la méthode " * method_name(MOON_METHOD)
        return rec
    end
    rec.mok = true
    if DISPERSION
        for c in (1, 3)
            off = channel_offset(MOON_METHOD, chans[c], sel[1], sel[2], rec, mref)
            if off === nothing
                push!(warn, (c == 1 ? "R" : "B") * " non mesurée")
            elseif c == 1
                rec.mrv = off
            else
                rec.mbv = off
            end
        end
    end
    rec.medge && push!(warn, "Lune au bord du champ")
    rec.mstatus = isempty(warn) ? "OK" : "OK (" * join(warn, ", ") * ")"
    return rec
end

##################################################################################################################
# TRANSLATION SOUS-PIXEL
##################################################################################################################

# noyau cubique de Keys (a = -0,5, "Catmull-Rom")
@inline function keys_cubic(t::Float64)
    a = -0.5
    t = abs(t)
    t <= 1 && return (a + 2) * t^3 - (a + 3) * t^2 + 1
    t < 2 && return a * t^3 - 5a * t^2 + 8a * t - 4a
    return 0.0
end

#############################################################################################
# décalage d (pixels) le long d'un axe : la valeur de sortie en x est celle de l'entrée en x - d,
# interpolée sur les pixels x + o[k] avec les poids w[k] (identiques pour tous les pixels : la
# translation est constante). Sortie : (o, w)
#############################################################################################
function shift_taps(d::Float64, method::Symbol)
    u = -d
    if method == :integer
        return [round(Int, u)], Float32[1]
    end
    iu = floor(Int, u)
    f = u - iu
    f < 1e-6 && return [iu], Float32[1]
    (1 - f) < 1e-6 && return [iu + 1], Float32[1]
    if method == :bilinear
        return [iu, iu + 1], Float32[1 - f, f]
    elseif method == :bicubic
        return [iu - 1, iu, iu + 1, iu + 2],
               Float32[keys_cubic(1 + f), keys_cubic(f), keys_cubic(1 - f), keys_cubic(2 - f)]
    else
        error("SHIFT_METHOD inconnue : ", method, " (attendu :bicubic, :bilinear ou :integer)")
    end
end

#############################################################################################
# translation d'un plan [y, x] de (dx, dy) pixels (dx > 0 : vers la droite, dy > 0 : vers le bas).
# Interpolation séparable (passe en x puis en y). Les pixels sans donnée sont mis à 0.
# Sortie : (plan translaté, zone valide (x1, x2, y1, y2))
#############################################################################################
function shift_plane(src::Matrix{Float32}, dx::Float64, dy::Float64, method::Symbol)
    ny, nx = size(src)
    ox, wx = shift_taps(dx, method)
    oy, wy = shift_taps(dy, method)
    x1 = max(1, 1 - minimum(ox)); x2 = min(nx, nx - maximum(ox))
    y1 = max(1, 1 - minimum(oy)); y2 = min(ny, ny - maximum(oy))
    (x1 <= x2 && y1 <= y2) || error("translation plus grande que l'image")

    tmp = zeros(Float32, ny, nx)
    Threads.@threads for x in x1:x2
        @inbounds for k in eachindex(ox)
            wk = wx[k]
            xs = x + ox[k]
            for y in 1:ny
                tmp[y, x] += wk * src[y, xs]
            end
        end
    end
    out = zeros(Float32, ny, nx)
    Threads.@threads for x in x1:x2
        @inbounds for y in y1:y2
            s = 0f0
            for k in eachindex(oy)
                s += wy[k] * tmp[y + oy[k], x]
            end
            out[y, x] = clamp(s, 0f0, 65535f0)      # dépassements de l'interpolation bicubique
        end
    end
    return out, (x1, x2, y1, y2)
end

##################################################################################################################
# FENÊTRE DE POINTAGE DE L'ÉTOILE (adaptée de HFD_star v9)
##################################################################################################################

# pyramide d'images réduites d'un facteur 2, 4, 8... par MAXIMUM sur des blocs 2x2
function reduce_max2(src::AbstractMatrix)
    s1, s2 = size(src)
    n1 = cld(s1, 2)
    n2 = cld(s2, 2)
    dst = Matrix{Float32}(undef, n1, n2)
    @inbounds for j in 1:n2, i in 1:n1
        m = -Inf32
        for jj in 2j-1:min(2j, s2), ii in 2i-1:min(2i, s1)
            v = Float32(src[ii, jj])
            (isfinite(v) && v > m) && (m = v)
        end
        dst[i, j] = m == -Inf32 ? NaN32 : m
    end
    return dst
end

function build_pyramid(img::AbstractMatrix, nlevels::Int)
    pyr = Matrix{Float32}[]
    nlevels >= 1 || return pyr
    push!(pyr, reduce_max2(img))
    for L in 2:nlevels
        push!(pyr, reduce_max2(pyr[L-1]))
    end
    return pyr
end

#############################################################################################
# tampon d'affichage buf[s, t] (s = colonne écran, t = ligne écran, t = 1 en haut) à partir d'une
# image source [y, x] (image d'origine si k = 1, niveau de la pyramide réduit d'un facteur k sinon)
#   xl, yt : bords gauche et haut de la vue (pixel i couvre [i-0.5, i+0.5]) ; z : zoom
#############################################################################################
function render_view!(buf::Matrix{Float32}, src::AbstractMatrix, k::Int,
                      xl::Float64, yt::Float64, z::Float64)
    vw, vh = size(buf)
    sy, sx = size(src)
    pidx = [floor(Int, (xl + (s - 0.5) / z - 0.5) / k) + 1 for s in 1:vw]
    qidx = [floor(Int, (yt + (t - 0.5) / z - 0.5) / k) + 1 for t in 1:vh]
    @inbounds for t in 1:vh
        q = qidx[t]
        if q < 1 || q > sy
            buf[:, t] .= NaN32
            continue
        end
        for s in 1:vw
            p = pidx[s]
            buf[s, t] = (1 <= p <= sx) ? Float32(src[q, p]) : NaN32
        end
    end
    return buf
end

circle_xy(x0, y0, r) = (x0 .+ r .* cos.(range(0, 2π, length=91)),
                        y0 .+ r .* sin.(range(0, 2π, length=91)))

fmt_adu(v::Real) = abs(v) >= 1000 ? @sprintf("%.0f", v) : @sprintf("%.4g", v)

function star_text(r)
    s = @sprintf("X = %.2f ; Y = %.2f\nincertitude : ± %.3f px (%s)\nFWHM ≈ %.1f px ; SNR ≈ %.0f\npic = %s ADU ; flux = %s",
                 r.x, r.y, r.err, r.method, r.fwhm, r.snr, fmt_adu(r.peak), fmt_adu(r.flux))
    r.sat && (s *= "\nATTENTION : étoile saturée")
    r.err > MAX_POS_ERROR && (s *= @sprintf("\nATTENTION : précision > %.1f px", MAX_POS_ERROR))
    return s
end

#############################################################################################
# fenêtre de pointage de l'étoile de référence sur la 1re image (image [y, x], y = 1 en haut) :
#   - curseurs X / Y : déplacement ; curseur ou molette : zoom ; curseur du seuil haut de visualisation
#   - clic gauche : réticule (cercle rouge : zone de recherche de l'étoile, carré jaune : fenêtre de
#     recherche SEARCH_BOX utilisée dans les images suivantes)
#   - "Mesurer" (touche M) : mesure de l'étoile la plus brillante près du réticule (croix verte)
#   - "Valider" (touche V) : fin du pointage, lancement du recentrage (mesure faite si nécessaire)
#   - "Quitter" ou fermeture de la fenêtre : abandon
#   - moon_win = (x, y) : clic droit = centre de la fenêtre de corrélation Lune (carré cyan MOON_WIN)
#   - moon_circle = (xc, yc, R) : cercle du limbe tracé en magenta (contrôle)
# Sortie : (star = résultat de refine_star, win = (x, y) ou nothing), ou nothing (abandon)
#############################################################################################
square_pts(p, hw) = [GM.Point2f(p[1] - hw, p[2] - hw), GM.Point2f(p[1] + hw, p[2] - hw),
                     GM.Point2f(p[1] + hw, p[2] + hw), GM.Point2f(p[1] - hw, p[2] + hw),
                     GM.Point2f(p[1] - hw, p[2] - hw)]

function pick_star_gui(img::Matrix{Float32}, file_name::AbstractString; moon_win=nothing, moon_circle=nothing)
    ny, nx = size(img)
    VW = VIEW_WIDTH
    VH = VIEW_HEIGHT

    # ---- statistiques globales sur un sous-échantillon (niveau de noir, plage du seuil haut)
    stp = max(1, round(Int, sqrt(nx * ny / 1.0e6)))
    sample = Float64.(filter(isfinite, vec(img[1:stp:end, 1:stp:end])))
    bg, sd = robust_background(sample)
    vmax = Float64(maximum(img))
    lo = sd > 0 ? sd : 1e-6 * max(abs(vmax), 1e-12)
    hi = max(vmax - bg, 2.0 * lo)
    white_of(s) = bg + exp(log(lo) + s * (log(hi) - log(lo)))
    w0 = max(quantile(sample, 0.995), bg + 10.0 * sd)
    s0 = clamp((log(max(w0 - bg, lo)) - log(lo)) / (log(hi) - log(lo)), 0.0, 1.0)

    # ---- niveaux de zoom : image entière, puis puissances de 2 jusqu'à MAX_ZOOM
    z_fit = min(VW / nx, VH / ny)
    zooms = [z_fit; [2.0^p for p in -8:round(Int, log2(MAX_ZOOM)) if 2.0^p > 1.05 * z_fit]]
    nlevels = z_fit < 1 ? floor(Int, log2(1 / z_fit) + 1e-9) : 0
    level_of(z) = z >= 1 ? 0 : min(nlevels, floor(Int, log2(1 / z) + 1e-9))
    println("préparation de l'affichage (", nlevels, " niveaux de réduction)...")
    pyr = build_pyramid(img, nlevels)
    blacks = Float64[bg]
    for L in 1:nlevels
        P = pyr[L]
        sL = max(1, stp >> L)
        push!(blacks, median(filter(isfinite, vec(P[1:sL:end, 1:sL:end]))))
    end

    x0u = (1 <= X_STAR <= nx) ? X_STAR : Float64(cld(nx, 2))
    y0u = (1 <= Y_STAR <= ny) ? Y_STAR : Float64(cld(ny, 2))

    # ---- fenêtre
    fig = GM.Figure(size=(VW + 450, VH + 150))
    ax = GM.Axis(fig[1, 1], width=VW, height=VH, backgroundcolor=:black, yreversed=true,
                 title=file_name * "  —  clic gauche : étoile" *
                       (moon_win === nothing ? "" : " ; clic droit : fenêtre Lune") * " ; molette : zoom",
                 xlabel="X (px)", ylabel="Y (px, vers le bas)")
    for name in (:rectanglezoom, :dragpan, :scrollzoom, :limitreset)
        try
            GM.deregister_interaction!(ax, name)
        catch
        end
    end

    # curseur vertical : haut du curseur = haut de l'image (y = 1)
    sl_y = GM.Slider(fig[1, 2], range=ny:-1:1, startvalue=cld(ny, 2), horizontal=false, height=VH)
    sl_x = GM.Slider(fig[2, 1], range=1:nx, startvalue=cld(nx, 2), width=VW)

    panel = fig[1, 3] = GM.GridLayout(tellheight=false, valign=:top)
    GM.Label(panel[1, 1], "Zoom (ou molette)", halign=:left, font=:bold)
    sl_z = GM.Slider(panel[2, 1], range=1:length(zooms), startvalue=1, width=320)
    lab_z = GM.Label(panel[3, 1], "", halign=:left)
    GM.Label(panel[4, 1], "Seuil haut de visualisation", halign=:left, font=:bold)
    sl_w = GM.Slider(panel[5, 1], range=0.0:0.001:1.0, startvalue=s0, width=320)
    lab_w = GM.Label(panel[6, 1], "", halign=:left)
    GM.Label(panel[7, 1], "Réticule (clic gauche)", halign=:left, font=:bold)
    lab_c = GM.Label(panel[8, 1], "", halign=:left)
    GM.Label(panel[9, 1], "Sous le curseur", halign=:left, font=:bold)
    lab_h = GM.Label(panel[10, 1], "", halign=:left)
    btn_m = GM.Button(panel[11, 1], label="Mesurer l'étoile (touche M)", width=320)
    btn_v = GM.Button(panel[12, 1], label="Valider et recentrer les images (touche V)", width=320)
    btn_q = GM.Button(panel[13, 1], label="Quitter sans traiter", width=320)
    GM.Label(panel[14, 1], "Étoile de référence", halign=:left, font=:bold)
    lab_r = GM.Label(panel[15, 1], "-", halign=:left, justification=:left)
    lab_mw = nothing
    if moon_win !== nothing
        GM.Label(panel[16, 1], "Fenêtre Lune (clic droit)", halign=:left, font=:bold)
        lab_mw = GM.Label(panel[17, 1], @sprintf("X = %.0f ; Y = %.0f", moon_win[1], moon_win[2]), halign=:left)
    end

    # ---- image
    buf = fill(NaN32, VW, VH)
    buf_obs = GM.Observable(buf)
    xs_obs = GM.Observable(collect(range(0.0, 1.0, length=VW)))
    ys_obs = GM.Observable(collect(range(0.0, 1.0, length=VH)))
    cr_obs = GM.Observable((Float32(bg), Float32(white_of(s0))))
    GM.heatmap!(ax, xs_obs, ys_obs, buf_obs; colormap=:grays, colorrange=cr_obs,
                nan_color=:black, interpolate=false)

    # ---- réticule, zone de recherche autour du réticule, fenêtre SEARCH_BOX, étoile mesurée
    cross = GM.Observable(GM.Point2f(x0u, y0u))
    GM.vlines!(ax, GM.lift(p -> [p[1]], cross), color=:red, linewidth=1)
    GM.hlines!(ax, GM.lift(p -> [p[2]], cross), color=:red, linewidth=1)
    GM.lines!(ax, GM.lift(p -> GM.Point2f.(circle_xy(p[1], p[2], GUI_SEARCH_RADIUS)...), cross),
              color=:red, linewidth=1.5)
    hb = SEARCH_BOX / 2
    GM.lines!(ax, GM.lift(p -> [GM.Point2f(p[1] - hb, p[2] - hb), GM.Point2f(p[1] + hb, p[2] - hb),
                                GM.Point2f(p[1] + hb, p[2] + hb), GM.Point2f(p[1] - hb, p[2] + hb),
                                GM.Point2f(p[1] - hb, p[2] - hb)], cross),
              color=:yellow, linewidth=1, linestyle=:dash)
    # ---- cercle du limbe (1re image) et fenêtre de corrélation Lune
    if moon_circle !== nothing
        GM.lines!(ax, GM.Point2f.(circle_xy(moon_circle[1], moon_circle[2], moon_circle[3])...),
                  color=:magenta, linewidth=1)
    end
    mwin = GM.Observable(moon_win === nothing ? GM.Point2f(NaN, NaN) : GM.Point2f(moon_win[1], moon_win[2]))
    if moon_win !== nothing
        GM.lines!(ax, GM.lift(p -> square_pts(p, MOON_WIN / 2), mwin), color=:cyan, linewidth=1.5)
    end
    meas = GM.Observable([GM.Point2f(NaN, NaN)])
    GM.lines!(ax, meas, color=:lime, linewidth=2)
    meas_c = GM.Observable([GM.Point2f(NaN, NaN)])
    GM.scatter!(ax, meas_c, marker=:cross, color=:lime, markersize=16)

    set_cross_label!() = (lab_c.text[] = @sprintf("X = %.0f ; Y = %.0f", cross[][1], cross[][2]))
    set_cross_label!()

    function update_view!()
        z = zooms[sl_z.value[]]
        L = level_of(z)
        k = 2^L
        cx = Float64(sl_x.value[])
        cy = Float64(sl_y.value[])
        xl = cx - VW / (2z)
        yt = cy - VH / (2z)
        if L == 0
            render_view!(buf, img, 1, xl, yt, z)
        else
            render_view!(buf, pyr[L], k, xl, yt, z)
        end
        xs_obs[] = [xl + (s - 0.5) / z for s in 1:VW]
        ys_obs[] = [yt + (t - 0.5) / z for t in 1:VH]
        GM.notify(buf_obs)
        white = white_of(sl_w.value[])
        cr_obs[] = (Float32(blacks[L+1]), Float32(max(white, blacks[L+1] + 0.1 * lo)))
        GM.limits!(ax, xl, xl + VW / z, yt, yt + VH / z)
        lab_z.text[] = z >= 1 ? @sprintf("x %.0f  (1 pixel image = %.0f pixels écran)", z, z) :
                                @sprintf("x 1/%.1f%s", 1 / z, sl_z.value[] == 1 ? "  (image entière)" : "")
        lab_w.text[] = fmt_adu(white) * " ADU  (fond " * fmt_adu(bg) * ", max " * fmt_adu(vmax) * ")"
        return nothing
    end

    for sl in (sl_x, sl_y, sl_z, sl_w)
        GM.on(_ -> update_view!(), sl.value)
    end

    measure_req = Ref(false)
    validate_req = Ref(false)
    quit_req = Ref(false)
    moved = Ref(true)          # réticule déplacé depuis la dernière mesure

    GM.on(GM.events(fig).mousebutton) do ev
        if ev.button == GM.Mouse.left && ev.action == GM.Mouse.press && GM.is_mouseinside(ax.scene)
            p = GM.mouseposition(ax.scene)
            cross[] = GM.Point2f(round(p[1]), round(p[2]))
            moved[] = true
            set_cross_label!()
        elseif moon_win !== nothing && ev.button == GM.Mouse.right && ev.action == GM.Mouse.press &&
               GM.is_mouseinside(ax.scene)
            p = GM.mouseposition(ax.scene)
            mwin[] = GM.Point2f(round(p[1]), round(p[2]))
            lab_mw.text[] = @sprintf("X = %.0f ; Y = %.0f", mwin[][1], mwin[][2])
        end
        return nothing
    end

    GM.on(GM.events(fig).mouseposition) do _
        GM.is_mouseinside(ax.scene) || return nothing
        p = GM.mouseposition(ax.scene)
        i = round(Int, p[1])
        j = round(Int, p[2])
        lab_h.text[] = (1 <= i <= nx && 1 <= j <= ny) ?
                       @sprintf("X = %d ; Y = %d ; ", i, j) * fmt_adu(img[j, i]) * " ADU" :
                       "hors de l'image"
        return nothing
    end

    GM.on(GM.events(fig).scroll) do sc
        GM.is_mouseinside(ax.scene) || return nothing
        dy = sc[2]
        dy == 0 && return nothing
        p = GM.mouseposition(ax.scene)
        if dy > 0
            GM.set_close_to!(sl_x, p[1])
            GM.set_close_to!(sl_y, p[2])
        end
        GM.set_close_to!(sl_z, clamp(sl_z.value[] + (dy > 0 ? 1 : -1), 1, length(zooms)))
        return nothing
    end

    GM.on(GM.events(fig).keyboardbutton) do ev
        if ev.action == GM.Keyboard.press
            ev.key == GM.Keyboard.m && (measure_req[] = true)
            ev.key == GM.Keyboard.v && (validate_req[] = true)
        end
        return nothing
    end
    GM.on(_ -> (measure_req[] = true), btn_m.clicks)
    GM.on(_ -> (validate_req[] = true), btn_v.clicks)
    GM.on(_ -> (quit_req[] = true), btn_q.clicks)

    update_view!()
    screen = GM.Screen()
    display(screen, fig)
    println("fenêtre de pointage ouverte : clic gauche sur l'étoile, \"Mesurer\", ",
            moon_win === nothing ? "" : "clic droit pour placer la fenêtre Lune (carré cyan), ", "puis \"Valider\"")

    result = nothing
    accepted = false
    while isopen(screen) && !quit_req[]
        if measure_req[] || (validate_req[] && (result === nothing || moved[]))
            measure_req[] = false
            p = cross[]
            lab_r.text[] = "mesure en cours..."
            r = locate_star(img, Float64(p[1]), Float64(p[2]), GUI_SEARCH_RADIUS, nothing)
            moved[] = false
            if r === nothing
                result = nothing
                validate_req[] = false
                meas[] = [GM.Point2f(NaN, NaN)]
                meas_c[] = [GM.Point2f(NaN, NaN)]
                lab_r.text[] = @sprintf("aucune étoile mesurable à moins de\n%d px du réticule", GUI_SEARCH_RADIUS)
                println("pointage : aucune étoile mesurable près de X = ", p[1], ", Y = ", p[2])
            else
                result = r
                meas[] = GM.Point2f.(circle_xy(r.x, r.y, max(r.fwhm, 2.0))...)
                meas_c[] = [GM.Point2f(r.x, r.y)]
                lab_r.text[] = star_text(r)
                println("étoile mesurée : ", replace(star_text(r), '\n' => " | "))
            end
        end
        if validate_req[]
            validate_req[] = false
            if result !== nothing
                accepted = true
                break
            end
        end
        sleep(0.05)
    end
    isopen(screen) && close(screen)
    println("fenêtre de pointage fermée")
    accepted || return nothing
    win = moon_win === nothing ? nothing : (Float64(mwin[][1]), Float64(mwin[][2]))
    return (star=result, win=win)
end

##################################################################################################################
# MESURES SUR TOUTES LES IMAGES
##################################################################################################################

Base.@kwdef mutable struct ImgRec
    path::String
    # ---- étoile
    ok::Bool = false
    status::String = "non traitée"
    x::Float64 = NaN    # position sur L (ou V si DISPERSION)
    y::Float64 = NaN
    err::Float64 = NaN
    xr::Float64 = NaN    # position dans R (DISPERSION)
    yr::Float64 = NaN
    okr::Bool = false
    xb::Float64 = NaN    # position dans B (DISPERSION)
    yb::Float64 = NaN
    okb::Bool = false
    flux::Float64 = NaN
    fwhm::Float64 = NaN
    snr::Float64 = NaN
    sat::Bool = false
    noise::Float64 = NaN
    weight::Float64 = 0.0
    dx::Vector{Float64} = Float64[]                 # translation de chaque couche
    dy::Vector{Float64} = Float64[]
    xp::Float64 = NaN    # centre de la fenêtre de recherche (position prédite)
    yp::Float64 = NaN
    # ---- Lune, mesure grossière (croissant sur l'image réduite)
    mx::Float64 = NaN
    my::Float64 = NaN
    okm::Bool = false
    medge::Bool = false
    mbox::NTuple{4, Int} = (1, 1, 1, 1)             # rectangle englobant le croissant
    mtop::Float64 = NaN                             # niveau du plateau saturé
    # ---- Lune, mesure fine (STACK_MOON)
    mok::Bool = false
    mstatus::String = "non traitée"
    mpos::Dict{Symbol, NTuple{4, Float64}} = Dict{Symbol, NTuple{4, Float64}}()   # méthode => (x, y, ±, qualité)
    mrv::NTuple{2, Float64} = (NaN, NaN)            # écart R - V de la Lune (DISPERSION)
    mbv::NTuple{2, Float64} = (NaN, NaN)            # écart B - V
    mdx::Vector{Float64} = Float64[]                # translation de chaque couche (somme Lune)
    mdy::Vector{Float64} = Float64[]
    mweight::Float64 = 0.0
end

#############################################################################################
# mesure de l'étoile dans une image (chans = plans R, V, B), autour de (x0, y0), carré de demi-côté half
#############################################################################################
function measure_image!(rec::ImgRec, chans::Vector{Matrix{Float32}}, x0::Real, y0::Real, half::Int, ref_flux)
    det = DISPERSION ? chans[2] : luminance(chans)
    rec.noise = noise_sigma(det)
    r = locate_star(det, x0, y0, half, ref_flux)
    if r === nothing
        rec.ok = false
        rec.status = "étoile non trouvée"
        return rec
    end
    rec.x = r.x; rec.y = r.y; rec.err = r.err
    rec.flux = r.flux; rec.fwhm = r.fwhm; rec.snr = r.snr; rec.sat = r.sat
    rec.ok = true
    warn = String[]
    r.sat && push!(warn, "saturée")
    r.err > MAX_POS_ERROR && push!(warn, @sprintf("précision > %.1f px", MAX_POS_ERROR))
    if DISPERSION
        rr = locate_star(chans[1], r.x, r.y, DISP_SEARCH_RADIUS, nothing)
        rb = locate_star(chans[3], r.x, r.y, DISP_SEARCH_RADIUS, nothing)
        if rr !== nothing
            rec.xr = rr.x; rec.yr = rr.y; rec.okr = true
        else
            push!(warn, "R non mesurée")
        end
        if rb !== nothing
            rec.xb = rb.x; rec.yb = rb.y; rec.okb = true
        else
            push!(warn, "B non mesurée")
        end
    end
    rec.status = isempty(warn) ? "OK" : "OK (" * join(warn, ", ") * ")"
    return rec
end

#############################################################################################
# moyenne avec rejet itératif à DISP_CLIP sigma. Sortie : (moyenne, écart-type, nombre retenu)
#############################################################################################
function clipped_mean(v::Vector{Float64})
    isempty(v) && return NaN, NaN, 0
    w = copy(v)
    if DISP_CLIP > 0
        for _ in 1:5
            length(w) < 3 && break
            m = mean(w)
            s = std(w)
            s > 0 || break
            keep = filter(x -> abs(x - m) <= DISP_CLIP * s, w)
            length(keep) == length(w) && break
            w = keep
        end
    end
    return mean(w), (length(w) > 1 ? std(w) : NaN), length(w)
end

# ligne d'avancement de la mesure : position de l'étoile, écart à la position prédite (centre de la
# fenêtre de recherche), position de la Lune
function print_progress(k::Int, n::Int, rec::ImgRec, ref::ImgRec)
    s = @sprintf("   %3d / %d  X = %9.3f  Y = %9.3f  ± %.3f px", k, n, rec.x, rec.y, rec.err)
    if k > 1
        d = hypot(rec.x - rec.xp, rec.y - rec.yp)
        s *= isfinite(d) ? @sprintf("  (écart à la prédiction %6.1f px)", d) : "                                  "
    end
    MOON_PRE_ALIGN && (s *= rec.okm ? @sprintf("  Lune X = %7.1f Y = %7.1f", rec.mx, rec.my) :
                                      "  Lune non trouvée          ")
    println(s, "   ", rec.status)
    if STACK_MOON
        d = moon_disp(rec, ref, MOON_METHOD)
        println(rec.mok && isfinite(d[1]) ?
                @sprintf("               Lune (%s) : dX = %+9.3f  dY = %+9.3f  ± %.3f px   %s",
                         method_name(MOON_METHOD), d[1], d[2], d[3], rec.mstatus) :
                "               Lune : " * rec.mstatus)
    end
    return nothing
end

# déplacement de la Lune (méthode m) par rapport à la 1re image : (dx, dy, ±) ; NaN si non mesuré
function moon_disp(rec::ImgRec, ref::ImgRec, m::Symbol)
    (haskey(rec.mpos, m) && haskey(ref.mpos, m)) || return (NaN, NaN, NaN)
    p = rec.mpos[m]
    q = ref.mpos[m]
    return (p[1] - q[1], p[2] - q[2], p[3])
end

function csvnum(v::Real, fmt::String="%.3f")
    s = isfinite(v) ? Printf.format(Printf.Format(fmt), v) : ""
    return CSV_DECIMAL_COMMA ? replace(s, '.' => ',') : s
end

#############################################################################################
# translation de chaque image (une translation par couche : tx[k][c], ty[k][c]), fichiers recentrés
# <image><reg_suffix>.tiff si SAVE_REG, somme pondérée (poids wts) normalisée => out_path
#############################################################################################
function build_stack(paths::Vector{String}, tx::Vector{Vector{Float64}}, ty::Vector{Vector{Float64}},
                     wts::Vector{Float64}, nch::Int, ny::Int, nx::Int, out_path::String, reg_suffix::String)
    acc = [zeros(Float32, ny, nx) for _ in 1:nch]
    wmap = zeros(Float32, ny, nx)
    t_all = time()
    for k in eachindex(paths)
        println()
        @printf("translation %d / %d : %s\n", k, length(paths), basename(paths[k]))
        chans, bits = read_image(paths[k])
        out = Vector{Matrix{Float32}}(undef, nch)
        x1, x2, y1, y2 = 1, nx, 1, ny
        for c in 1:nch
            out[c], (a1, a2, b1, b2) = shift_plane(chans[c], tx[k][c], ty[k][c], SHIFT_METHOD)
            x1 = max(x1, a1); x2 = min(x2, a2); y1 = max(y1, b1); y2 = min(y2, b2)
        end
        chans = nothing
        SAVE_REG && save_tiff(joinpath(dirname(paths[k]), splitext(basename(paths[k]))[1] * reg_suffix * ".tiff"),
                              out, bits)
        w = Float32(wts[k])
        for c in 1:nch
            A = acc[c]
            O = out[c]
            Threads.@threads for x in x1:x2
                @inbounds for y in y1:y2
                    A[y, x] += w * O[y, x]
                end
            end
        end
        Threads.@threads for x in x1:x2
            @inbounds for y in y1:y2
                wmap[y, x] += w
            end
        end
        out = nothing
        GC.gc()
    end

    # moyenne pondérée (chaque pixel : images qui le couvrent), puis normalisation 16 bits
    for c in 1:nch
        A = acc[c]
        @inbounds for i in eachindex(A)
            A[i] = wmap[i] > 0 ? A[i] / wmap[i] : 0f0
        end
    end
    ncov = count(>(0.999f0), wmap)
    @printf("\nzone couverte par toutes les images : %.1f %% du champ\n", 100 * ncov / length(wmap))
    factor = 1.0
    if NORMALIZE == :max
        stp = max(1, round(Int, sqrt(nx * ny / 4.0e6)))
        vals = Float64[]
        for c in 1:nch
            append!(vals, Float64.(vec(acc[c][1:stp:end, 1:stp:end])))
        end
        top = NORM_PERCENTILE >= 100 ? maximum(maximum(A) for A in acc) : quantile(vals, NORM_PERCENTILE / 100)
        top > 0 && (factor = 65535.0 / top)
        @printf("normalisation 16 bits : niveau %.1f ADU ramené à 65535 (facteur %.4f, identique pour toutes les couches)\n",
                top, factor)
    elseif NORMALIZE != :none
        error("NORMALIZE inconnu : ", NORMALIZE, " (attendu :max ou :none)")
    end
    if factor != 1.0
        f32 = Float32(factor)
        for A in acc
            A .*= f32
        end
    end
    save_tiff(out_path, acc, 16)
    @printf("%d image(s) recentrée(s) et additionnée(s) en %.0f s\n", length(paths), time() - t_all)
    return nothing
end

# poids (non normalisé) d'une image dans la somme étoile
function star_weight(r)
    WEIGHT_MODE == :equal && return 1.0
    WEIGHT_MODE == :noise && return 1.0 / r.noise^2
    WEIGHT_MODE == :flux && return r.flux
    return r.flux / r.noise^2
end

# poids normalisés (somme = 1) ; erreur si tous nuls
function normalize_weights!(w::Vector{Float64}, what::String)
    for i in eachindex(w)
        (isfinite(w[i]) && w[i] > 0) || (w[i] = 0.0)
    end
    s = sum(w; init=0.0)
    s > 0 || error("somme des poids nulle (", what, ") : essayer des poids égaux (:equal)")
    w ./= s
    return w
end

# ajustement d'une parabole v(k) = a + b k + c k² (k = numéro d'image) ; écart-type des résidus
function parabola_rms(k::Vector{Float64}, v::Vector{Float64})
    n = length(v)
    n >= 5 || return NaN
    kc = k .- mean(k)
    M = hcat(ones(n), kc, kc .^ 2)
    coef = M \ v
    r = v .- M * coef
    return sqrt(sum(abs2, r) / (n - 3))
end

##################################################################################################################
# PROGRAMME PRINCIPAL
##################################################################################################################
function main()

println("register_star_earthshine - version 2.0")
println("pré-recentrage sur la Lune : ", MOON_PRE_ALIGN ? "oui" : "non")
println("somme recentrée sur la Lune : ", STACK_MOON ? "oui (méthode " * method_name(MOON_METHOD) *
        (MOON_COMPARE ? ", comparaison des 3 méthodes)" : ")") : "non")
println("dispersion atmosphérique : ", DISPERSION ? "oui (mesure dans R, V, B)" : "non (mesure sur la luminance L)")
println("threads Julia : ", Threads.nthreads())
MOON_METHOD in MOON_METHODS || error("MOON_METHOD inconnue : ", MOON_METHOD, " (attendu :correlation, :limb ou :crescent)")
MOON_LIMB_SIDE in (:both, :bright, :earthshine) || error("MOON_LIMB_SIDE inconnu : ", MOON_LIMB_SIDE)
MOON_WEIGHT_MODE in (:equal, :noise) || error("MOON_WEIGHT_MODE inconnu : ", MOON_WEIGHT_MODE, " (attendu :equal ou :noise)")

# -----------------------------------------------------------------------------------------------
# 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
    # *_reg.tiff, *_reg_star.tiff, *_reg_moon.tiff, *_stack_reg_star.tiff, *_stack_reg_moon.tiff
    is_output = f -> occursin(r"_reg(_star|_moon)?\.tiff?$"i, basename(f))
    for f in filter(is_output, files)
        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) sélectionné(s) ; image de référence : ", basename(files[1]))
n == 1 && println("ATTENTION : une seule image, rien à recentrer")
for (k, f) in enumerate(files)
    @printf("   %3d  %s\n", k, basename(f))
end

# 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(), " "), RAW_CACHE ? "  (cache des TIFF convertis activé)" : "")
    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

try
# -----------------------------------------------------------------------------------------------
# 1re image : Lune (cercle du limbe, fenêtre de corrélation par défaut), pointage de l'étoile
# -----------------------------------------------------------------------------------------------
println()
println("lecture de l'image de référence : ", basename(files[1]))
chans1, _ = read_image(files[1])
ny, nx = size(chans1[1])
nch = length(chans1)
println("dim_x = ", nx, ", dim_y = ", ny)
DISPERSION && nch != 3 && error("DISPERSION = true demande des images RVB (image en niveaux de gris)")
det1 = DISPERSION ? chans1[2] : luminance(chans1)

mref = MoonRef()
win0 = nothing
circ = nothing
if STACK_MOON
    moon1 = moon_center(det1)
    moon1 === nothing && error("STACK_MOON : Lune (croissant saturé) non détectée sur la 1re image")
    mref.mx, mref.my = moon1.x, moon1.y
    mref.methods = MOON_COMPARE ? collect(MOON_METHODS) : [MOON_METHOD]
    @printf("Lune (1re image) : centre du croissant X = %.1f ; Y = %.1f ; plateau %.0f ADU\n",
            moon1.x, moon1.y, moon1.top)
    # cercle du limbe (méthode :limb, contrôle visuel, position par défaut de la fenêtre de corrélation)
    print("recherche du cercle du limbe sur la 1re image... ")
    g = moon_circle_guess(det1)
    lc = g === nothing ? nothing : limb_circle(det1, g.xc, g.yc, g.R, MOON_SAT_FRAC * moon1.top, noise_sigma(det1))
    if lc === nothing
        println("échec")
        MOON_METHOD == :limb && error("méthode :limb : limbe non trouvé sur la 1re image")
        filter!(!=(:limb), mref.methods)
        :limb in MOON_METHODS && MOON_COMPARE && println("   méthode limbe retirée de la comparaison")
    else
        mref.cx, mref.cy, mref.R = lc.xc, lc.yc, lc.R
        circ = (lc.xc, lc.yc, lc.R)
        @printf("X = %.2f ; Y = %.2f ; R = %.2f px (%d points : %d côté croissant, %d côté lumière cendrée ; rms %.2f px)\n",
                lc.xc, lc.yc, lc.R, lc.n, lc.nb, lc.ne, lc.rms)
    end
    if :correlation in mref.methods
        if 1 <= MOON_WIN_X <= nx && 1 <= MOON_WIN_Y <= ny
            win0 = (MOON_WIN_X, MOON_WIN_Y)
        elseif circ !== nothing && hypot(moon1.x - circ[1], moon1.y - circ[2]) > 1
            # par défaut : sur le limbe, à l'opposé du croissant (côté lumière cendrée)
            ux = (circ[1] - moon1.x) / hypot(circ[1] - moon1.x, circ[2] - moon1.y)
            uy = (circ[2] - moon1.y) / hypot(circ[1] - moon1.x, circ[2] - moon1.y)
            rr = max(circ[3] - MOON_WIN / 4, 0.0)
            win0 = (clamp(circ[1] + rr * ux, 1.0, Float64(nx)), clamp(circ[2] + rr * uy, 1.0, Float64(ny)))
        else
            win0 = (moon1.x, moon1.y)
        end
    end
end

if USE_GUI
    sel = pick_star_gui(det1, basename(files[1]); moon_win=win0, moon_circle=circ)
    sel === nothing && (println("pointage abandonné - programme arrêté"); return nothing)
    x0, y0 = sel.star.x, sel.star.y
    win = sel.win
else
    x0, y0 = X_STAR, Y_STAR
    win = win0
end

# fenêtre de corrélation de référence
if STACK_MOON && win !== nothing
    tpl = moon_template(det1, win[1], win[2], MOON_SAT_FRAC * plateau_level(det1))
    if tpl === nothing || tpl.frac < MOON_CORR_MIN_FRAC
        msg = tpl === nothing ? "image trop petite" :
              @sprintf("%.0f %% seulement de pixels non saturés", 100 * tpl.frac)
        MOON_METHOD == :correlation && error("fenêtre de corrélation inutilisable (", msg, ") : la placer côté lumière cendrée")
        println("ATTENTION : fenêtre de corrélation inutilisable (", msg, ") : méthode corrélation retirée de la comparaison")
        filter!(!=(:correlation), mref.methods)
    else
        mref.xr, mref.yr, mref.T, mref.tbad = tpl.xr, tpl.yr, tpl.T, tpl.tbad
        @printf("fenêtre de corrélation Lune : X = %d..%d ; Y = %d..%d (%.0f %% de pixels non saturés)\n",
                first(tpl.xr), last(tpl.xr), first(tpl.yr), last(tpl.yr), 100 * tpl.frac)
    end
end
det1 = nothing

# -----------------------------------------------------------------------------------------------
# mesure de l'étoile (et de la Lune) dans toutes les images
# -----------------------------------------------------------------------------------------------
recs = [ImgRec(path=f) for f in files]
println()
println("mesure de l'étoile de référence", STACK_MOON ? " et de la Lune" : "", " dans les ", n,
        " images (couche ", DISPERSION ? "V" : "L", ")")
measure_image!(recs[1], chans1, x0, y0, GUI_SEARCH_RADIUS, nothing)
recs[1].ok || error("étoile de référence non mesurable sur la 1re image près de X = $x0, Y = $y0")
recs[1].xp, recs[1].yp = x0, y0
ref_flux = recs[1].flux
xlast, ylast = recs[1].x, recs[1].y

# pré-recentrage sur la Lune : image d'ancrage = 1re image (puis image précédente si TRACK_PREVIOUS)
anchor = nothing                     # (X étoile, Y étoile, X Lune, Y Lune) de l'image d'ancrage
(MOON_PRE_ALIGN || STACK_MOON) && measure_moon!(recs[1], chans1)
if MOON_PRE_ALIGN
    if recs[1].okm
        anchor = (recs[1].x, recs[1].y, recs[1].mx, recs[1].my)
        @printf("Lune (1re image) : X = %.1f ; Y = %.1f%s\n", recs[1].mx, recs[1].my,
                recs[1].medge ? "  (ATTENTION : croissant touchant le bord du champ)" : "")
    else
        println("ATTENTION : Lune non détectée sur la 1re image : pas de pré-recentrage sur la Lune")
    end
end
if STACK_MOON
    measure_moon_fine!(recs[1], chans1, mref)
    recs[1].mok || error("Lune non mesurable sur la 1re image (référence de la somme Lune) : ", recs[1].mstatus)
end
chans1 = nothing
print_progress(1, n, recs[1], recs[1])
for k in 2:n
    rec = recs[k]
    try
        chans, _ = read_image(files[k])
        if size(chans[1]) != (ny, nx) || length(chans) != nch
            rec.status = "dimensions ou nombre de couches différents de la 1re image"
            rec.mstatus = rec.status
        else
            xc, yc = TRACK_PREVIOUS ? (xlast, ylast) : (recs[1].x, recs[1].y)
            moon_note = ""
            (MOON_PRE_ALIGN || STACK_MOON) && measure_moon!(rec, chans)
            if MOON_PRE_ALIGN
                if rec.okm && anchor !== nothing
                    xc = anchor[1] + (rec.mx - anchor[3])
                    yc = anchor[2] + (rec.my - anchor[4])
                    rec.medge && (moon_note = "Lune au bord du champ")
                elseif !rec.okm
                    moon_note = "Lune non trouvée"
                end
            end
            rec.xp, rec.yp = xc, yc
            measure_image!(rec, chans, xc, yc, SEARCH_BOX ÷ 2, ref_flux)
            isempty(moon_note) || (rec.status *= " ; " * moon_note)
            if rec.ok
                xlast, ylast = rec.x, rec.y
                TRACK_PREVIOUS && rec.okm && (anchor = (rec.x, rec.y, rec.mx, rec.my))
            end
            if STACK_MOON
                try
                    measure_moon_fine!(rec, chans, mref)
                catch err
                    err isa InterruptException && rethrow()
                    rec.mok = false
                    rec.mstatus = "ERREUR : " * first(replace(sprint(showerror, err), '\n' => ' '), 120)
                end
            end
        end
    catch err
        err isa InterruptException && rethrow()
        rec.ok = false
        rec.status = "ERREUR : " * first(replace(sprint(showerror, err), '\n' => ' '), 120)
    end
    print_progress(k, n, rec, recs[1])
    GC.gc()
end
okrecs = filter(r -> r.ok, recs)
stem1 = splitext(basename(files[1]))[1]
out_dir = dirname(files[1])

##################################################################################################################
# SOMME RECENTRÉE SUR L'ÉTOILE
##################################################################################################################

# -----------------------------------------------------------------------------------------------
# translations : delta = position de référence (1re image) - position mesurée
# -----------------------------------------------------------------------------------------------
xref, yref = recs[1].x, recs[1].y
corr_r = (0.0, 0.0)
corr_b = (0.0, 0.0)
if DISPERSION
    rgx = [r.xr - r.x for r in okrecs if r.okr]; rgy = [r.yr - r.y for r in okrecs if r.okr]
    bgx = [r.xb - r.x for r in okrecs if r.okb]; bgy = [r.yb - r.y for r in okrecs if r.okb]
    mrx, srx, nrx = clipped_mean(rgx); mry, sry, nry = clipped_mean(rgy)
    mbx, sbx, nbx = clipped_mean(bgx); mby, sby, nby = clipped_mean(bgy)
    (nrx > 0 && nbx > 0) || error("dispersion : étoile non mesurée dans R ou dans B")
    println()
    println("-"^100)
    println("DISPERSION ATMOSPHÉRIQUE (ÉTOILE) : écart moyen de la position de l'étoile par rapport à la couche V")
    @printf("   R - V : dX = %+.3f px (σ = %.3f, %d images)   dY = %+.3f px (σ = %.3f, %d images)\n",
            mrx, srx, nrx, mry, sry, nry)
    @printf("   B - V : dX = %+.3f px (σ = %.3f, %d images)   dY = %+.3f px (σ = %.3f, %d images)\n",
            mbx, sbx, nbx, mby, sby, nby)
    @printf("   écart R - B : %.3f px\n", hypot(mrx - mbx, mry - mby))
    corr_r = (-mrx, -mry)
    corr_b = (-mbx, -mby)
    @printf("   correction ajoutée aux delta de V : couche R (%+.3f ; %+.3f) px, couche B (%+.3f ; %+.3f) px\n",
            corr_r[1], corr_r[2], corr_b[1], corr_b[2])
    println("-"^100)
end
for r in okrecs
    dxg = xref - r.x
    dyg = yref - r.y
    if nch == 3
        r.dx = [dxg + corr_r[1], dxg, dxg + corr_b[1]]
        r.dy = [dyg + corr_r[2], dyg, dyg + corr_b[2]]
    else
        r.dx = [dxg]
        r.dy = [dyg]
    end
end

# -----------------------------------------------------------------------------------------------
# poids des images dans la somme
# -----------------------------------------------------------------------------------------------
WEIGHT_MODE in (:equal, :noise, :flux, :snr) ||
    error("WEIGHT_MODE inconnu : ", WEIGHT_MODE, " (attendu :equal, :noise, :flux ou :snr)")
wst = [star_weight(r) for r in okrecs]
normalize_weights!(wst, "somme étoile")
for (r, w) in zip(okrecs, wst)
    r.weight = w
end

# -----------------------------------------------------------------------------------------------
# tableau de contrôle dans le terminal
# -----------------------------------------------------------------------------------------------
println()
println("="^130)
@printf("ÉTOILE - TRANSLATIONS À APPLIQUER (position de référence : X = %.3f ; Y = %.3f, couche %s de la 1re image)\n",
        xref, yref, DISPERSION ? "V" : "L")
println("="^130)
if DISPERSION
    println("  n°   X étoile   Y étoile    ±(px)  |   R-V dX   R-V dY   B-V dX   B-V dY  |  dX R     dY R     dX V     dY V     dX B     dY B   |  poids  statut")
else
    println("  n°   X étoile   Y étoile    ±(px)  |  delta X   delta Y  |   FWHM     flux        bruit   |  poids   statut")
end
for (k, r) in enumerate(recs)
    if !r.ok
        @printf("  %3d  %s\n", k, r.status)
        continue
    end
    if DISPERSION
        @printf("  %3d  %9.3f  %9.3f  %7.3f  | %+8.3f %+8.3f %+8.3f %+8.3f | %+7.3f  %+7.3f  %+7.3f  %+7.3f  %+7.3f  %+7.3f | %6.3f  %s\n",
                k, r.x, r.y, r.err, r.xr - r.x, r.yr - r.y, r.xb - r.x, r.yb - r.y,
                r.dx[1], r.dy[1], r.dx[2], r.dy[2], r.dx[3], r.dy[3], r.weight, r.status)
    else
        @printf("  %3d  %9.3f  %9.3f  %7.3f  | %+8.3f  %+8.3f  | %6.2f  %10.0f  %8.2f   | %6.3f   %s\n",
                k, r.x, r.y, r.err, r.dx[1], r.dy[1], r.fwhm, r.flux, r.noise, r.weight, r.status)
    end
end
println("="^130)
println("delta X > 0 : image déplacée vers la droite ; delta Y > 0 : vers le bas (origine en haut à gauche)")
@printf("%d image(s) recentrée(s) sur %d ; translation %s ; poids %s\n", length(okrecs), n,
        SHIFT_METHOD, WEIGHT_MODE)
println()

# -----------------------------------------------------------------------------------------------
# fichier CSV des mesures
# -----------------------------------------------------------------------------------------------
if SAVE_CSV
    csv_path = joinpath(out_dir, stem1 * "_reg_star.csv")
    head = ["fichier", "statut", "X", "Y", "incertitude_px", "FWHM_px", "flux", "bruit", "poids"]
    nch == 3 ? append!(head, ["X_R", "Y_R", "X_B", "Y_B", "dX_R", "dY_R", "dX_V", "dY_V", "dX_B", "dY_B"]) :
               append!(head, ["dX", "dY"])
    MOON_PRE_ALIGN && append!(head, ["X_pred", "Y_pred", "X_Lune", "Y_Lune"])
    open(csv_path, "w") do io
        println(io, join(head, CSV_SEP))
        for r in recs
            row = [basename(r.path), r.status, csvnum(r.x), csvnum(r.y), csvnum(r.err), csvnum(r.fwhm, "%.2f"),
                   csvnum(r.flux, "%.0f"), csvnum(r.noise, "%.3f"), csvnum(r.weight, "%.5f")]
            if nch == 3
                append!(row, [csvnum(r.xr), csvnum(r.yr), csvnum(r.xb), csvnum(r.yb)])
                append!(row, r.ok ? [csvnum(v) for v in (r.dx[1], r.dy[1], r.dx[2], r.dy[2], r.dx[3], r.dy[3])] :
                                    fill("", 6))
            else
                append!(row, r.ok ? [csvnum(r.dx[1]), csvnum(r.dy[1])] : ["", ""])
            end
            MOON_PRE_ALIGN && append!(row, [csvnum(r.xp, "%.1f"), csvnum(r.yp, "%.1f"),
                                            csvnum(r.mx, "%.1f"), csvnum(r.my, "%.1f")])
            println(io, join(row, CSV_SEP))
        end
    end
    println("enregistré : ", csv_path)
end

##################################################################################################################
# MESURES FINES DE LA LUNE : DISPERSION, TRANSLATIONS, COMPARAISON DES MÉTHODES, CSV
##################################################################################################################
mrecs = ImgRec[]
if STACK_MOON
    mrecs = filter(r -> r.mok, recs)
    mm = MOON_METHOD
    ref1 = recs[1]
    mcorr_r = (0.0, 0.0)
    mcorr_b = (0.0, 0.0)
    if DISPERSION
        rvx = [r.mrv[1] for r in mrecs if isfinite(r.mrv[1])]; rvy = [r.mrv[2] for r in mrecs if isfinite(r.mrv[2])]
        bvx = [r.mbv[1] for r in mrecs if isfinite(r.mbv[1])]; bvy = [r.mbv[2] for r in mrecs if isfinite(r.mbv[2])]
        mrx, srx, nrx = clipped_mean(rvx); mry, sry, nry = clipped_mean(rvy)
        mbx, sbx, nbx = clipped_mean(bvx); mby, sby, nby = clipped_mean(bvy)
        println()
        println("-"^100)
        println("DISPERSION ATMOSPHÉRIQUE (LUNE, méthode ", method_name(mm), ") : écart moyen par rapport à la couche V")
        if nrx > 0 && nbx > 0
            @printf("   R - V : dX = %+.3f px (σ = %.3f, %d images)   dY = %+.3f px (σ = %.3f, %d images)\n",
                    mrx, srx, nrx, mry, sry, nry)
            @printf("   B - V : dX = %+.3f px (σ = %.3f, %d images)   dY = %+.3f px (σ = %.3f, %d images)\n",
                    mbx, sbx, nbx, mby, sby, nby)
            mcorr_r = (-mrx, -mry)
            mcorr_b = (-mbx, -mby)
            if DISPERSION && !isempty(okrecs)
                @printf("   (étoile : R - V = (%+.3f ; %+.3f), B - V = (%+.3f ; %+.3f))\n",
                        -corr_r[1], -corr_r[2], -corr_b[1], -corr_b[2])
            end
        else
            println("   Lune non mesurée dans R ou dans B : correction de la dispersion de l'ÉTOILE utilisée")
            mcorr_r = corr_r
            mcorr_b = corr_b
        end
        println("-"^100)
    end
    for r in mrecs
        d = moon_disp(r, ref1, mm)
        dxg = -d[1]
        dyg = -d[2]
        if nch == 3
            r.mdx = [dxg + mcorr_r[1], dxg, dxg + mcorr_b[1]]
            r.mdy = [dyg + mcorr_r[2], dyg, dyg + mcorr_b[2]]
        else
            r.mdx = [dxg]
            r.mdy = [dyg]
        end
    end
    wsm = [MOON_WEIGHT_MODE == :equal ? 1.0 : 1.0 / r.noise^2 for r in mrecs]
    isempty(mrecs) || normalize_weights!(wsm, "somme Lune")
    for (r, w) in zip(mrecs, wsm)
        r.mweight = w
    end

    # ---- tableau des translations Lune
    println()
    println("="^130)
    @printf("LUNE (méthode %s) - TRANSLATIONS À APPLIQUER (référence : 1re image, couche %s)\n",
            method_name(mm), DISPERSION ? "V" : "L")
    println("="^130)
    if DISPERSION
        println("  n°   dX Lune    dY Lune     ±(px)  |   R-V dX   R-V dY   B-V dX   B-V dY  |  dX R     dY R     dX V     dY V     dX B     dY B   |  poids  statut")
    else
        println("  n°   dX Lune    dY Lune     ±(px)  |  delta X   delta Y  |  poids   statut")
    end
    for (k, r) in enumerate(recs)
        if !r.mok
            @printf("  %3d  %s\n", k, r.mstatus)
            continue
        end
        d = moon_disp(r, ref1, mm)
        if DISPERSION
            @printf("  %3d  %+9.3f  %+9.3f  %7.3f  | %+8.3f %+8.3f %+8.3f %+8.3f | %+7.3f  %+7.3f  %+7.3f  %+7.3f  %+7.3f  %+7.3f | %6.3f  %s\n",
                    k, d[1], d[2], d[3], r.mrv[1], r.mrv[2], r.mbv[1], r.mbv[2],
                    r.mdx[1], r.mdy[1], r.mdx[2], r.mdy[2], r.mdx[3], r.mdy[3], r.mweight, r.mstatus)
        else
            @printf("  %3d  %+9.3f  %+9.3f  %7.3f  | %+8.3f  %+8.3f  | %6.3f   %s\n",
                    k, d[1], d[2], d[3], r.mdx[1], r.mdy[1], r.mweight, r.mstatus)
        end
    end
    println("="^130)
    @printf("%d image(s) recentrée(s) sur la Lune sur %d ; poids %s\n", length(mrecs), n, MOON_WEIGHT_MODE)

    # ---- comparaison des méthodes
    if length(mref.methods) > 1
        println()
        println("="^130)
        println("COMPARAISON DES MÉTHODES : déplacement de la Lune par rapport à la 1re image (pixels)")
        println("="^130)
        hdr = "  n°  "
        for m in mref.methods
            hdr *= @sprintf("| %-12s dX       dY       ±    ", method_name(m))
        end
        println(hdr)
        for (k, r) in enumerate(recs)
            line = @sprintf("  %3d ", k)
            for m in mref.methods
                d = moon_disp(r, ref1, m)
                line *= isfinite(d[1]) ? @sprintf("| %+12.3f %+8.3f %6.3f   ", d[1], d[2], d[3]) :
                                         "|            -                    "
            end
            println(line)
        end
        println("-"^130)
        println("écarts à la méthode ", method_name(mm), " (moyenne ; écart-type) :")
        for m in mref.methods
            m == mm && continue
            ex = Float64[]
            ey = Float64[]
            for r in recs
                a = moon_disp(r, ref1, m)
                b = moon_disp(r, ref1, mm)
                if isfinite(a[1]) && isfinite(b[1])
                    push!(ex, a[1] - b[1])
                    push!(ey, a[2] - b[2])
                end
            end
            length(ex) >= 2 || continue
            @printf("   %-12s : dX %+7.3f ; %6.3f    dY %+7.3f ; %6.3f   (%d images)\n",
                    method_name(m), mean(ex), std(ex), mean(ey), std(ey), length(ex))
        end
        if length(okrecs) >= 5
            println("indicateur de précision : déplacement Lune - déplacement étoile (mouvement propre de la Lune),")
            println("écart-type des résidus d'une parabole en fonction du numéro d'image (cadence régulière supposée ;")
            println("inclut le bruit de mesure de l'étoile) :")
            for m in mref.methods
                kk = Float64[]
                vx = Float64[]
                vy = Float64[]
                for (k, r) in enumerate(recs)
                    d = moon_disp(r, ref1, m)
                    (r.ok && isfinite(d[1])) || continue
                    push!(kk, k)
                    push!(vx, d[1] - (r.x - recs[1].x))
                    push!(vy, d[2] - (r.y - recs[1].y))
                end
                @printf("   %-12s : X %6.3f px   Y %6.3f px   (%d images)\n",
                        method_name(m), parabola_rms(kk, vx), parabola_rms(kk, vy), length(kk))
            end
        end
        println("="^130)
    end

    # ---- CSV Lune
    if SAVE_CSV
        csv_path = joinpath(out_dir, stem1 * "_reg_moon.csv")
        head = ["fichier", "statut_Lune", "X_croissant_grossier", "Y_croissant_grossier"]
        for m in mref.methods
            t = String(m)
            qn = m == :correlation ? "gain" : m == :limb ? "rayon" : "surface"
            append!(head, ["dX_" * t, "dY_" * t, "incertitude_" * t, qn * "_" * t])
        end
        nch == 3 ? append!(head, ["R-V_dX", "R-V_dY", "B-V_dX", "B-V_dY",
                                  "dX_R", "dY_R", "dX_V", "dY_V", "dX_B", "dY_B"]) : append!(head, ["dX", "dY"])
        push!(head, "poids")
        open(csv_path, "w") do io
            println(io, join(head, CSV_SEP))
            for r in recs
                row = [basename(r.path), r.mstatus, csvnum(r.mx, "%.1f"), csvnum(r.my, "%.1f")]
                for m in mref.methods
                    d = moon_disp(r, ref1, m)
                    q = haskey(r.mpos, m) ? r.mpos[m][4] : NaN
                    append!(row, [csvnum(d[1]), csvnum(d[2]), csvnum(d[3]), csvnum(q, "%.4f")])
                end
                if nch == 3
                    append!(row, [csvnum(r.mrv[1]), csvnum(r.mrv[2]), csvnum(r.mbv[1]), csvnum(r.mbv[2])])
                    append!(row, r.mok ? [csvnum(v) for v in (r.mdx[1], r.mdy[1], r.mdx[2], r.mdy[2], r.mdx[3], r.mdy[3])] :
                                         fill("", 6))
                else
                    append!(row, r.mok ? [csvnum(r.mdx[1]), csvnum(r.mdy[1])] : ["", ""])
                end
                push!(row, csvnum(r.mweight, "%.5f"))
                println(io, join(row, CSV_SEP))
            end
        end
        println("enregistré : ", csv_path)
    end
end

##################################################################################################################
# SOMMES : ÉTOILE, PUIS LUNE
##################################################################################################################
println()
println("#"^60, "  SOMME RECENTRÉE SUR L'ÉTOILE")
build_stack([r.path for r in okrecs], [r.dx for r in okrecs], [r.dy for r in okrecs], [r.weight for r in okrecs],
            nch, ny, nx, joinpath(out_dir, stem1 * "_stack_reg_star.tiff"), "_reg")
if STACK_MOON
    if isempty(mrecs)
        println("ATTENTION : aucune image recentrée sur la Lune : pas de somme Lune")
    else
        println()
        println("#"^60, "  SOMME RECENTRÉE SUR LA LUNE (", method_name(MOON_METHOD), ")")
        build_stack([r.path for r in mrecs], [r.mdx for r in mrecs], [r.mdy for r in mrecs],
                    [r.mweight for r in mrecs], nch, ny, nx, joinpath(out_dir, stem1 * "_stack_reg_moon.tiff"),
                    "_reg_moon")
    end
end
println("\nterminé")

finally
    clear_raw_cache()          # TIFF temporaires des fichiers RAW (RAW_CACHE = true)
end

if !isinteractive()
    println("appuyer sur Entrée pour terminer")
    readline()
end
return nothing

end # fin de la fonction main()

# lancement du programme
main()