using Pkg # ----------------------------------------------------------------------------------------------- # Installation des packages : mettre à false une fois les packages installés une première fois. # Pkg.add() vérifie le registre à chaque exécution, ce qui coûte plusieurs secondes/minutes # à chaque lancement du script pour rien si les packages sont déjà présents. # ----------------------------------------------------------------------------------------------- const FIRST_RUN = false #const FIRST_RUN = true if FIRST_RUN for pkg in ["Statistics", "FITSIO", "Images", "NativeFileDialog"] Pkg.add(pkg) end end using Statistics using Printf # fits file : https://juliaastro.org/FITSIO.jl/stable/ using FITSIO # package pour filtrage des images (imfilter, Kernel.gaussian) using Images # package pour sélection interactive du répertoire et du fichier (fenêtres de dialogue) using NativeFileDialog ################################################################################################################## # Fonctions du programme: # read the input image (fit file, 3 couleurs) # pour chaque pixel, calcul de 3 masques flous complets (sigma_i, k_i, n_sigma_threshold_i, # i=1,2,3), combinés par interpolation linéaire en fonction de la distance au centre lunaire # (exprimée en unités de rayon solaire R_sol) - voir détail plus bas # write result in fit, postfixe "_u" + k1_k2_k3 sur 2 digits chacun (+ "_partial" si convolution partielle activée) # (+ "_L" pour le fichier correspondant à la couche de luminance, si process_L=true) # # Entrées/paramètres : # directory et nom de l'image d'entrée : choisis via IHM (pick_folder / pick_file) # # 3 triplets (sigma_i, k_i, n_sigma_threshold_i), i=1,2,3, en dur dans le script : chacun # définit un masque flou complet, comme en V1 : # masque flou_i = Input_image + k_i * diff_i * mask_seuil_i # diff_i = Input_image - Gauss(Input_image, sigma_i) [/ G_masque_lune_i si partial_conv] # seuil ("niveau") : exprimé en multiple (n_sigma_threshold_i) de l'écart-type du bruit de # diff_i, estimé par un estimateur robuste MAD (Median Absolute Deviation), insensible aux # vrais contours à fort gradient contrairement à un simple écart-type. Un k_i = 0 désactive # simplement le renforcement pour le masque i (pas de souci numérique). # # combinaison radiale des 3 masques flous, en fonction de la distance r au centre lunaire # (x_center_moon, y_center_moon), exprimée en unités de rayon solaire R_sol (en pixels) : # - r <= r1 : 100% masque 1 # - r1 < r < r2 : interpolation linéaire masque 1 -> masque 2 # - r = r2 : 100% masque 2 # - r2 < r < r3 : interpolation linéaire masque 2 -> masque 3 # - r >= r3 : 100% masque 3 # avec r1 < r2 < r3 (croissant, du centre vers l'extérieur). # # convolution partielle (partial_conv true/false) : masque_lune généré directement dans le # code (disque lunaire = 0, extérieur = 1) à partir de son centre (x_center_moon, y_center_moon, # partagé avec le calcul de r) et de son rayon (r_moon), tous modifiables en dur ; même # principe que dans le script v9, mais masque_lune généré ici par calcul plutôt que lu # depuis un fichier "mask_small.fit". Une gaussienne du masque lunaire est calculée par # sigma_i (3 au total), puisque chaque masque flou a son propre sigma. # # mode de traitement (process_RGB / process_L, chacun true/false, indépendants) : le masque # flou combiné peut être calculé sur les 3 canaux RGB, et/ou sur la couche de luminance L # extraite de R,G,B via des coefficients de pondération classiques (ITU-R BT.601), pour # améliorer le rapport signal/bruit. Le temps de calcul étant très court, les deux modes # peuvent être activés simultanément (deux fichiers de sortie distincts sont alors écrits). # L'extraction de la luminance se fait sur l'image d'entrée, avant tout calcul de masque flou. # # historique (TSE2026-5-radial-unsharp-V1.jl) : # v1 : première version, un seul masque flou (k, sigma, n_sigma_threshold) # v2 : ajout de la convolution partielle (masque_lune généré, G_masque_lune, partial_conv) # v3 : ajout du choix du mode de traitement RGB et/ou L (process_RGB, process_L) # # V2 (ce script, TSE2026-5-radial-unsharp-V2.jl) : # combinaison de 3 masques flous avec interpolation radiale (R_sol, r1/r2/r3) ################################################################################################################### ############################################################################################# # calcule le masque flou (avec seuillage MAD et convolution partielle optionnelle) sur UN # canal/couche d'image (2D). Réutilisée pour chacun des 3 canaux RGB et pour la couche L, afin # d'éviter de dupliquer la logique de calcul. # # Entrées : # img_channel : canal/couche d'entrée (2D) # gauss_kernel : noyau gaussien (dépend de sigma, construit une seule fois dans main()) # k, n_sigma_threshold : paramètres du masque flou (voir en-tête du script) # partial_conv : booléen, active la convolution partielle # G_masque_lune : gaussienne du masque lunaire (utilisée seulement si partial_conv=true, # sinon ignorée : peut valoir `nothing`) # const_MAD_to_sigma : facteur de conversion MAD -> écart-type équivalent gaussien # label : texte affiché en console (ex: "canal 1", "couche L") # # Sortie : # canal/couche résultat du masque flou (2D) ############################################################################################# function process_channel(img_channel, gauss_kernel, k, n_sigma_threshold, partial_conv, G_masque_lune, const_MAD_to_sigma, label::String) # calcul de la version floutée (gaussienne) du canal/couche courant img_gauss = imfilter(img_channel, gauss_kernel) # diff = Input_image - Gauss(Input_image, sigma) diff = img_channel .- img_gauss # convolution partielle : diff = (Input_image - Gauss(Input_image, sigma)) / G_masque_lune if partial_conv diff = diff ./ G_masque_lune end # estimation robuste de l'écart-type du bruit de diff via la MAD (Median Absolute Deviation) : # insensible aux vrais contours à fort gradient, contrairement à un simple std() median_diff = median(diff) mad_diff = median(abs.(diff .- median_diff)) robust_sigma = const_MAD_to_sigma * mad_diff # seuil ("niveau") pour le canal/couche courant threshold = n_sigma_threshold * robust_sigma println(" ", label, " : robust_sigma = ", @sprintf("%.3f", robust_sigma), " ; threshold = ", @sprintf("%.3f", threshold)) # masque booléen (0/1) : renforcement actif uniquement là où |diff| >= threshold mask_seuil = abs.(diff) .>= threshold # masque flou = Input_image + k * diff * mask_seuil return img_channel .+ k .* diff .* mask_seuil end ############################################################################################# # calcule, pour chaque pixel de r_map (distance au centre lunaire, en unités de rayon solaire), # les 3 poids d'interpolation radiale entre les 3 masques flous (w1+w2+w3 = 1 partout) : # r <= r1 : w1=1, w2=0, w3=0 (100% masque 1) # r1 < r < r2 : interpolation linéaire masque 1 -> 2 (w1 décroît, w2 croît) # r2 <= r < r3 : interpolation linéaire masque 2 -> 3 (w2 décroît, w3 croît) # r >= r3 : w1=0, w2=0, w3=1 (100% masque 3) # nécessite r1 < r2 < r3 (croissant, du centre vers l'extérieur). ############################################################################################# function radial_weights(r_map, r1, r2, r3) w1 = similar(r_map) w2 = similar(r_map) w3 = similar(r_map) @inbounds for idx in eachindex(r_map) r = r_map[idx] if r <= r1 w1[idx] = 1.0; w2[idx] = 0.0; w3[idx] = 0.0 elseif r < r2 t = (r - r1) / (r2 - r1) w1[idx] = 1.0 - t; w2[idx] = t; w3[idx] = 0.0 elseif r < r3 t = (r - r2) / (r3 - r2) w1[idx] = 0.0; w2[idx] = 1.0 - t; w3[idx] = t else w1[idx] = 0.0; w2[idx] = 0.0; w3[idx] = 1.0 end end return w1, w2, w3 end ############################################################################################# # calcule le masque flou combiné (3 masques, interpolation radiale) sur UN canal/couche # d'image (2D). Chaque masque i est calculé complètement (Gauss, convolution partielle # éventuelle, seuillage MAD, application de k_i) via process_channel, puis les 3 résultats # sont combinés pixel par pixel selon les poids radiaux (w1, w2, w3). # # Entrées : # kernels_3 : tuple/vecteur des 3 noyaux gaussiens (Kernel.gaussian(sigma_i)) # ks_3 : tuple/vecteur des 3 valeurs de k # thresholds_3 : tuple/vecteur des 3 valeurs de n_sigma_threshold # G_masque_lune_3 : tuple/vecteur des 3 G_masque_lune (un par sigma_i), ou `nothing` # si partial_conv=false # weights_3 : tuple (w1, w2, w3) des cartes de poids radiaux (2D) ############################################################################################# function combine_radial_masks(img_channel, kernels_3, ks_3, thresholds_3, partial_conv, G_masque_lune_3, weights_3, const_MAD_to_sigma, label) w1, w2, w3 = weights_3 result1 = process_channel(img_channel, kernels_3[1], ks_3[1], thresholds_3[1], partial_conv, partial_conv ? G_masque_lune_3[1] : nothing, const_MAD_to_sigma, "$(label) - masque 1") result2 = process_channel(img_channel, kernels_3[2], ks_3[2], thresholds_3[2], partial_conv, partial_conv ? G_masque_lune_3[2] : nothing, const_MAD_to_sigma, "$(label) - masque 2") result3 = process_channel(img_channel, kernels_3[3], ks_3[3], thresholds_3[3], partial_conv, partial_conv ? G_masque_lune_3[3] : nothing, const_MAD_to_sigma, "$(label) - masque 3") return w1 .* result1 .+ w2 .* result2 .+ w3 .* result3 end # tout le programme principal est placé dans une fonction : en Julia, le code # écrit directement au niveau global est beaucoup plus lent (variables non typées). # Une fonction est compilée avec des types connus => gain de vitesse important, sans # changer la logique du programme. function main() # chronométrage du temps total de calcul t_start = time() println("Nombre de threads Julia disponibles : ", Threads.nthreads()) if Threads.nthreads() == 1 println("(pour utiliser plusieurs threads, relancer avec : julia --threads=4 nom_du_script.jl)") end # ----------------------------------------------------------------------------------------------- # paramètres des 3 masques flous à combiner (en dur pour cette V2) # chaque triplet (sigma_i, k_i, n_sigma_threshold_i) définit un masque flou complet, comme en V1 # ----------------------------------------------------------------------------------------------- version_couleur = false binning = true if version_couleur==true # reglage pour version couleur offset = 250 sigma1 = 6.0 k1 = 2.0 n_sigma_threshold1 = 0.0 sigma2 = 12.0 k2 = 1.0 n_sigma_threshold2 = 0.75 sigma3 = 30.0 k3 = 0.0 n_sigma_threshold3 = 0.0 else if binning == false # reglage pour version NB offset = 0 et binning 1x1 sigma1 = 6.0 k1 = 2.0 n_sigma_threshold1 = 0.0 sigma2 = 12.0 k2 = 1.5 n_sigma_threshold2 = 0.50 sigma3 = 30.0 k3 = 1.0 n_sigma_threshold3 = 0.50 else # reglage pour version NB offset = 0 et binning 2x2 #sigma1 = 3.0 sigma1 = 10.0 k1 = 2.0 n_sigma_threshold1 = 0.0 #sigma2 = 6.0 sigma2 = 10.0 k2 = 1.5 n_sigma_threshold2 = 0.50 sigma3 = 15.0 k3 = 1.0 n_sigma_threshold3 = 0.50 end end # ----------------------------------------------------------------------------------------------- # interpolation radiale entre les 3 masques flous, en fonction de la distance au centre lunaire # exprimée en unités de rayon solaire (voir en-tête du script) - nécessite r1 < r2 < r3 # ----------------------------------------------------------------------------------------------- R_sol = 656.0 # rayon solaire, en pixels if binning == true R_sol = R_sol/2.0 end r1 = 1.4 # en dessous de r1 (en unités de R_sol) : 100% masque 1 r2 = 2.0 # à r2 : 100% masque 2 (interpolation 1->2 entre r1 et r2) r3 = 3.0 # au-dessus de r3 : 100% masque 3 (interpolation 2->3 entre r2 et r3) @assert r1 < r2 < r3 "r1 < r2 < r3 requis (ordre croissant, du centre vers l'extérieur)" # ----------------------------------------------------------------------------------------------- # convolution partielle : réduit les effets de rebond au limbe lunaire (voir en-tête du script) # ----------------------------------------------------------------------------------------------- partial_conv = false # géométrie du disque lunaire (masque_lune), en dur, modifiable #x_center_moon = 4055.0 #y_center_moon = 2698.0 # version initiale x_center_moon = 4055.0 y_center_moon = 2700.0 r_moon = 635.0 if binning == true x_center_moon = x_center_moon/2.0 y_center_moon = y_center_moon/2.0 r_moon = r_moon/2.0 end # ----------------------------------------------------------------------------------------------- # mode de traitement : sur les canaux RGB, et/ou sur la couche de luminance L (voir en-tête du # script). Les deux peuvent être activés simultanément. # ----------------------------------------------------------------------------------------------- process_RGB = false process_L = true # ----------------------------------------------------------------------------------------------- # carte de contrôle (debug) : écrit un fichier supplémentaire "..._L_Rmap.fit", copie de l'image # u3_L avec les 3 cercles r1, r2, r3 (rayons pivots, en pixels = r_i * R_sol) tracés en # surimpression (valeur = ADU max de l'image), pour repérer visuellement les zones de transition. # Nécessite process_L=true (le canal L suffit pour ce repérage). # ----------------------------------------------------------------------------------------------- export_Rmap = true circle_thickness_px = 4 # épaisseur des cercles tracés, en pixels # sélection interactive du répertoire de travail (remplace le cd() en dur) work_dir = pick_folder(pwd()) isempty(work_dir) && error("Aucun répertoire sélectionné - programme arrêté") cd(work_dir) println(pwd()) # nom du fichier image en entrée : sélection interactive (remplace le nom de fichier en dur) selected_input = pick_file(work_dir; filterlist="fit,fits") isempty(selected_input) && error("Aucun fichier sélectionné - programme arrêté") file_name_input = basename(selected_input) # lecture de l'image d'entrée (fit) f_in = FITS(file_name_input) # lecture complète du header header = read_header(f_in[1]) # lecture de l'image d'entrée (tableau avec trois couches RGB : X,Y,couleurs), résultat en Float64 image_input = Float64.(read(f_in[1])) close(f_in) dim_x = header["NAXIS1"] dim_y = header["NAXIS2"] println("dim_x = ", dim_x) println("dim_y = ", dim_y) println("masque 1 : sigma1 = ", sigma1, " ; k1 = ", k1, " ; n_sigma_threshold1 = ", n_sigma_threshold1) println("masque 2 : sigma2 = ", sigma2, " ; k2 = ", k2, " ; n_sigma_threshold2 = ", n_sigma_threshold2) println("masque 3 : sigma3 = ", sigma3, " ; k3 = ", k3, " ; n_sigma_threshold3 = ", n_sigma_threshold3) println("R_sol = ", R_sol, " px ; r1 = ", r1, " ; r2 = ", r2, " ; r3 = ", r3, " (unités R_sol)") println("partial_conv = ", partial_conv) println("process_RGB = ", process_RGB) println("process_L = ", process_L) println("export_Rmap = ", export_Rmap, " (épaisseur cercles = ", circle_thickness_px, " px)") if export_Rmap && !process_L println("ATTENTION : export_Rmap=true nécessite process_L=true - carte de contrôle non générée") end # stat sur l'image d'entrée println("max R = ", maximum(image_input[:,:,1])) println("max V = ", maximum(image_input[:,:,2])) println("max B = ", maximum(image_input[:,:,3])) println("mean R = ", mean(image_input[:,:,1])) println("mean G = ", mean(image_input[:,:,2])) println("mean B = ", mean(image_input[:,:,3])) # ----------------------------------------------------------------------------------------------- # extraction de la couche de luminance L, si nécessaire (avant tout calcul de masque flou), # via les coefficients de pondération classiques R,G,B -> L (ITU-R BT.601) # ----------------------------------------------------------------------------------------------- if process_L w_R = 0.299 w_G = 0.587 w_B = 0.114 println("extraction de la couche L (poids : R=", w_R, ", G=", w_G, ", B=", w_B, ")") image_L = w_R .* image_input[:,:,1] .+ w_G .* image_input[:,:,2] .+ w_B .* image_input[:,:,3] end # les 3 noyaux gaussiens, un par masque flou (chacun dépend de son propre sigma_i) kernels_3 = (Kernel.gaussian(sigma1), Kernel.gaussian(sigma2), Kernel.gaussian(sigma3)) ks_3 = (k1, k2, k3) thresholds_3 = (n_sigma_threshold1, n_sigma_threshold2, n_sigma_threshold3) # ----------------------------------------------------------------------------------------------- # convolution partielle : construction du masque lunaire (identique quel que soit le masque # flou / canal) et de sa version floutée G_masque_lune_i, une par sigma_i (3 gaussiennes # différentes du même masque_lune) # ----------------------------------------------------------------------------------------------- if partial_conv println("calcul du masque lunaire (centre = (", x_center_moon, ",", y_center_moon, "), rayon = ", r_moon, " px)") # masque_lune : 0 sur le disque lunaire, 1 en dehors masque_lune = ones(dim_x, dim_y) r_moon2 = r_moon^2 @inbounds for j in 1:dim_y, i in 1:dim_x if (i - x_center_moon)^2 + (j - y_center_moon)^2 <= r_moon2 masque_lune[i,j] = 0.0 end end epsilon = 1e-6 println("calcul des gaussiennes du masque lunaire (sigma1, sigma2, sigma3)") # une gaussienne du masque lunaire par sigma_i, remplacement des valeurs proches de 0 # (< epsilon) par 1 pour éviter division par 0 G_masque_lune_list = Vector{Matrix{Float64}}(undef, 3) for i in 1:3 G_i = imfilter(masque_lune, kernels_3[i]) G_masque_lune_list[i] = ifelse.(G_i .< epsilon, 1.0, G_i) end G_masque_lune_3 = Tuple(G_masque_lune_list) else G_masque_lune_3 = nothing end # ----------------------------------------------------------------------------------------------- # carte de la distance au centre lunaire, exprimée en unités de rayon solaire (R_sol), et poids # d'interpolation radiale entre les 3 masques flous (voir en-tête du script) # ----------------------------------------------------------------------------------------------- println("calcul de la carte de distance radiale (centre = (", x_center_moon, ",", y_center_moon, "), R_sol = ", R_sol, " px)") r_map = zeros(dim_x, dim_y) @inbounds for j in 1:dim_y, i in 1:dim_x r_map[i,j] = sqrt((i - x_center_moon)^2 + (j - y_center_moon)^2) / R_sol end weights_3 = radial_weights(r_map, r1, r2, r3) println("calcul du masque flou combiné (3 masques, interpolation radiale)") # facteur de conversion MAD -> écart-type équivalent gaussien const_MAD_to_sigma = 1.4826 # nom de base commun aux fichiers de sortie : postfixe "_u" + valeurs de k1, k2, k3 # (chacune codée sur 2 digits, valeur*10 arrondie, ex: k1=2.0, k2=0.7, k3=0.0 -> "_u20_07_00") # (+ "_partial" si convolution partielle activée) ; le mode L ajoute en plus le postfixe "_L" k_str = @sprintf("u%02d_%02d_%02d", round(Int, k1*10), round(Int, k2*10), round(Int, k3*10)) file_name_base = chop(file_name_input, tail=4) * "_" * k_str if partial_conv file_name_base = file_name_base * "_partial" end # G_masque_lune_3 n'est défini que si partial_conv=true ; combine_radial_masks n'y accède # que dans ce cas (transmis tel quel, `nothing` sinon) if process_RGB println("traitement du masque flou combiné sur les canaux RGB :") image_unsharp = zeros(dim_x, dim_y, 3) # @views évite de copier chaque canal couleur avant de le filtrer @views for c in [1,2,3] image_unsharp[:,:,c] = combine_radial_masks(image_input[:,:,c], kernels_3, ks_3, thresholds_3, partial_conv, G_masque_lune_3, weights_3, const_MAD_to_sigma, "canal $c") end file_name_output_RGB = file_name_base * ".fit" # écriture du fichier FIT avec le masque flou appliqué (3 couleurs, header d'origine conservé) println("écriture du fichier de sortie RGB : ", file_name_output_RGB) f_out = FITS(file_name_output_RGB, "w") write(f_out, image_unsharp[:,:,:]; header=header) close(f_out) end if process_L println("traitement du masque flou combiné sur la couche L :") image_L_unsharp = combine_radial_masks(image_L, kernels_3, ks_3, thresholds_3, partial_conv, G_masque_lune_3, weights_3, const_MAD_to_sigma, "couche L") file_name_output_L = file_name_base * "_L.fit" # écriture du fichier FIT avec le masque flou appliqué (image 2D, une seule couche) : pas de # recopie du header d'origine (prévu pour 3 couleurs, NAXIS3=3) afin d'éviter toute # incohérence avec des données 2D ; un header minimal est généré automatiquement println("écriture du fichier de sortie L : ", file_name_output_L) f_out = FITS(file_name_output_L, "w") write(f_out, image_L_unsharp) close(f_out) # ------------------------------------------------------------------------------------------- # carte de contrôle (debug) : copie de image_L_unsharp avec les 3 cercles r1, r2, r3 # (rayons en pixels = r_i * R_sol) tracés en surimpression, valeur = ADU max de l'image # ------------------------------------------------------------------------------------------- if export_Rmap println("génération de la carte de contrôle (cercles r1, r2, r3)") image_L_Rmap = copy(image_L_unsharp) max_val = maximum(image_L_unsharp) half_thick = circle_thickness_px / 2 # r_map est déjà calculé plus haut (en unités de R_sol) ; distance en pixels = r_map * R_sol dist_px = r_map .* R_sol for r_circle_px in (r1 * R_sol, r2 * R_sol, r3 * R_sol) @inbounds for idx in eachindex(dist_px) if abs(dist_px[idx] - r_circle_px) <= half_thick image_L_Rmap[idx] = max_val end end end file_name_output_Rmap = file_name_base * "_L_Rmap.fit" println("écriture du fichier de contrôle : ", file_name_output_Rmap) f_out = FITS(file_name_output_Rmap, "w") write(f_out, image_L_Rmap) close(f_out) end end # affichage du temps de calcul total elapsed = time() - t_start println(@sprintf("Temps de calcul total : %.1f s (%.1f min)", elapsed, elapsed/60)) end # fin de la fonction main() # lancement du programme main()