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", "Printf", "Images", "FileIO", "TiffImages", "NativeFileDialog", "Plots"]
        Pkg.add(pkg)
    end
end

using Statistics
using Printf

# package pour manipulation d'images (colorview, N0f16, RGB, etc.)
using Images

# FileIO (load/save) + backend TIFF (lecture/écriture de TIF 16 bits)
using FileIO
using TiffImages

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

# package pour affichage des 3 courbes ADU_C = f(ADU_L)
using Plots

##################################################################################################################
# Fonctions du programme :
#   lecture de l'image N&B d'entrée (TIF 16 bits), par IHM (pick_folder / pick_file)
#   construction de 3 fonctions de colorisation ADU_C = f_C(ADU_L), C = R, V (Vert/G), B :
#     - pour ADU_L > 125            : ADU_C = ADU_L (pas de changement, au-dessus du point pivot)
#     - pour ADU_L <= 125           : Hermite cubique entre (L_dark, cible_C) et (125, 125)
#   application aux 3 canaux, écriture du fichier de sortie TIF 16 bits, suffixe "_color"
#   affichage (fenêtre) des 3 courbes ADU_C = f(ADU_L) pour vérification visuelle
#
# Entrées/paramètres :
#      directory et nom de l'image d'entrée : choisis via IHM (pick_folder / pick_file), TIF N&B
#
#      point pivot (en dur) : R = V = B = 125, tangente = 1 en ce point (continuité de tangente
#      avec la fonction identité utilisée au-dessus de 125)
#
#      point cible (en dur) : ADU_R = 36, ADU_V = 48, ADU_B = 70, atteint pour ADU_L = L_dark
#      (L_dark = médiane des ADU des 10% pixels les plus sombres de l'image d'entrée, calculée
#      dynamiquement pour chaque image)
#
# Choix de modélisation (v2, suite à comparaison avec la courbe "Courbes" de Photoshop) :
#      une parabole (v1) est entièrement déterminée par pivot + tangente-pivot + cible (3
#      conditions, 3 inconnues) : la tangente au point cible n'est donc pas libre, elle est
#      imposée par le calcul. Pour les canaux où le cible s'écarte beaucoup du pivot (typiquement
#      le Bleu), cette tangente imposée s'effondre, voire change de signe -> aplat non désiré
#      près du point cible, différent de l'allure obtenue avec l'outil "Courbes" de Photoshop.
#
#      v2 utilise donc un véritable Hermite cubique (2 points + 2 tangentes) :
#        - au pivot (125, 125)      : tangente fixée à 1 (continuité avec l'identité, imposée
#                                      par le cahier des charges)
#        - au point cible (L_dark, cible_C) : tangente RÉGLABLE (tangent_target_R/V/B ci-dessous),
#                                      un paramètre "en dur" de plus, à ajuster par canal pour
#                                      reproduire l'allure de la courbe Photoshop (forte pente
#                                      près du cible, qui se radoucit progressivement jusqu'au
#                                      pivot). C'est la 4e contrainte qui manquait à la parabole.
#
# historique :
#   v1 : parabole par canal (tangente au point cible imposée par le calcul)
#   v2 : Hermite cubique par canal, tangente au point cible réglable (voir ci-dessus)
#   v3 : point cible défini en coordonnées HSL, converties en RGB par hsl_to_rgb()
#   v4 (_HSB) : point cible défini en coordonnées HSB/HSV (Teinte/Saturation/Brillance), converties
#        en RGB par hsb_to_rgb() - pour comparaison directe avec la pipette + panneau Couleur de
#        Photoshop, qui affiche du HSB et non du HSL (voir échange avec l'utilisateur : les deux
#        modèles donnent des S et 3e composante différents pour un même RGB, d'où la confusion
#        initiale). Le graphe de vérification reste exprimé en coordonnées RVB.
###################################################################################################################

#############################################################################################
# conversion HSB/HSV -> RGB (échelle 0-255). Mêmes principes que HSL, mais la 3e composante B
# (Brillance/Value) est définie comme B = Max(R,V,B) (et non (Max+Min)/2 comme en HSL) : c'est
# le modèle utilisé par le panneau Couleur / la pipette de Photoshop.
#
# Entrées :
#   H : teinte, en degrés (0-360)
#   S : saturation, fraction (0-1)
#   B : brillance (Value), fraction (0-1)
#
# Sortie :
#   (R, V, B), chacun en échelle 0-255 (Float64)
#############################################################################################
function hsb_to_rgb(H::Float64, S::Float64, B::Float64)
    C = B * S
    Hp = H / 60.0
    X = C * (1.0 - abs(mod(Hp, 2.0) - 1.0))
    m = B - C

    if 0.0 <= Hp < 1.0
        R1, G1, B1 = C, X, 0.0
    elseif 1.0 <= Hp < 2.0
        R1, G1, B1 = X, C, 0.0
    elseif 2.0 <= Hp < 3.0
        R1, G1, B1 = 0.0, C, X
    elseif 3.0 <= Hp < 4.0
        R1, G1, B1 = 0.0, X, C
    elseif 4.0 <= Hp < 5.0
        R1, G1, B1 = X, 0.0, C
    else
        R1, G1, B1 = C, 0.0, X
    end

    return (R1 + m) * 255.0, (G1 + m) * 255.0, (B1 + m) * 255.0
end

#############################################################################################
# interpolation de Hermite cubique entre deux points (x0,y0,m0) et (x1,y1,m1), avec continuité
# de valeur ET de tangente aux deux extrémités (f(x0)=y0, f'(x0)=m0, f(x1)=y1, f'(x1)=m1).
# Fonctionne aussi en extrapolation (x < x0 ou x > x1), la formule restant valable pour tout t.
#
# Entrées :
#   x            : ADU_L (échelle 0-255), scalaire
#   x0, y0, m0   : point + tangente au 1er point (ici : point cible, L_dark)
#   x1, y1, m1   : point + tangente au 2nd point  (ici : point pivot, 125)
#
# Sortie :
#   ADU_C (échelle 0-255), scalaire (non clampé, voir colorize_function)
#############################################################################################
function cubic_hermite(x::Float64, x0::Float64, y0::Float64, m0::Float64,
                        x1::Float64, y1::Float64, m1::Float64)
    h = x1 - x0
    t = (x - x0) / h
    t2 = t * t
    t3 = t2 * t

    h00 = 2t3 - 3t2 + 1
    h10 = t3 - 2t2 + t
    h01 = -2t3 + 3t2
    h11 = t3 - t2

    return h00 * y0 + h10 * h * m0 + h01 * y1 + h11 * h * m1
end

#############################################################################################
# fonction de colorisation pour UN canal (R, V ou B) : identité au-dessus du pivot, Hermite
# cubique en dessous (interpolation entre L_dark et pivot, extrapolation en dessous de L_dark).
# Un clamp final ramène le résultat dans [0, 255] par sécurité (extrapolation aux ADU très
# sombres, ou dépassement local dû à une tangente cible élevée).
#
# Entrées :
#   x             : ADU_L (échelle 0-255), scalaire ou tableau (fonction "broadcastée" via .)
#   L_dark        : abscisse du point cible
#   target        : ordonnée du point cible (ADU_C visé en L_dark)
#   tangent_target: tangente réglable au point cible (voir en-tête du script)
#   pivot         : abscisse/ordonnée du point pivot (125 par défaut)
#   tangent_pivot : tangente au pivot (1.0 par défaut, imposé par le cahier des charges)
#
# Sortie :
#   ADU_C (échelle 0-255), même forme que x
#############################################################################################
function colorize_function(x::Float64, L_dark::Float64, target::Float64, tangent_target::Float64,
                            pivot::Float64=125.0, tangent_pivot::Float64=1.0)
    if x > pivot
        return x
    else
        val = cubic_hermite(x, L_dark, target, tangent_target, pivot, pivot, tangent_pivot)
        return clamp(val, 0.0, 255.0)
    end
end

#############################################################################################
# calcule L_dark : médiane des ADU (échelle 0-255) des 10% pixels les plus sombres de l'image
# d'entrée.
#
# Entrée :
#   image_L8 : image N&B (échelle 0-255), tableau 2D
#
# Sortie :
#   L_dark (Float64)
#############################################################################################
function compute_L_dark(image_L8)
    values_sorted = sort(vec(image_L8))
    n_dark = max(1, round(Int, 0.10 * length(values_sorted)))
    dark_values = values_sorted[1:n_dark]
    return median(dark_values)
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()

# -----------------------------------------------------------------------------------------------
# point pivot (fixe) et point cible (fixe), en dur (échelle 0-255, voir cahier des charges)
# -----------------------------------------------------------------------------------------------
pivot = 125.0
tangent_pivot = 1.0   # imposé par le cahier des charges (continuité avec l'identité au-dessus de 125)

# -----------------------------------------------------------------------------------------------
# point cible, défini directement en coordonnées HSB/HSV (Teinte/Saturation/Brillance) - le même
# modèle que le panneau Couleur / la pipette de Photoshop, pour comparaison directe avec l'image
# de sortie. Converti en RGB (échelle 0-255) par hsb_to_rgb() ci-dessous ; c'est ce triplet RGB
# qui sert ensuite aux fonctions de colorisation par canal (comme avant).
#
# valeurs de départ : équivalent HSB du triplet RGB (21.2, 42.4, 84.8) précédemment validé
# (H=220°, S=60%, L=20% en HSL <=> H=220°, S=75%, B=33% en HSB, voir échange avec l'utilisateur).
# -----------------------------------------------------------------------------------------------
H_target = 220.0   # teinte, en degrés
S_target = 0.50    # saturation, fraction (0-1)
B_target = 0.28    # brillance (Value), fraction (0-1)

target_R, target_V, target_B = hsb_to_rgb(H_target, S_target, B_target)
println("point cible (HSB) : H = ", H_target, "° ; S = ", @sprintf("%.0f%%", S_target * 100),
        " ; B = ", @sprintf("%.0f%%", B_target * 100))
println("point cible (RVB, échelle 0-255) : R = ", @sprintf("%.2f", target_R),
        " ; V = ", @sprintf("%.2f", target_V), " ; B = ", @sprintf("%.2f", target_B))

# -----------------------------------------------------------------------------------------------
# tangentes au point cible (L_dark), réglables par canal : c'est le paramètre qui manquait à la
# parabole (v1) pour reproduire l'allure de la courbe "Courbes" de Photoshop. Plus la valeur est
# grande, plus la courbe repart "raide" depuis le point cible avant de se radoucir vers le pivot.
# Ces tangentes sont exprimées dans le système RVB (comme les fonctions de colorisation
# elles-mêmes), pas en HSB. A ajuster en regardant le graphe affiché par le script.
# -----------------------------------------------------------------------------------------------
tangent_target_R = 0.35
tangent_target_V = 0.25
tangent_target_B = 0.2

# facteur de conversion ADU 16 bits (0..65535) <-> ADU 8 bits équivalent (0..255)
# 65535 = 255 * 257, conversion exacte dans les deux sens
const_16_to_8 = 257.0

# -----------------------------------------------------------------------------------------------
# sélection interactive du répertoire de travail et du fichier d'entrée (TIF N&B, 16 bits)
# -----------------------------------------------------------------------------------------------
work_dir = pick_folder(pwd())
isempty(work_dir) && error("Aucun répertoire sélectionné - programme arrêté")
cd(work_dir)

println(pwd())

selected_input = pick_file(work_dir; filterlist="tif,tiff")
isempty(selected_input) && error("Aucun fichier sélectionné - programme arrêté")
file_name_input = basename(selected_input)

# -----------------------------------------------------------------------------------------------
# lecture de l'image d'entrée (TIF 16 bits, N&B) et passage en Float64, échelle 0-65535 (ADU réels)
# -----------------------------------------------------------------------------------------------
println("lecture du fichier d'entrée : ", file_name_input)
img_in = load(file_name_input)

# image_L16 : ADU réels (0..65535), tableau 2D Float64, quel que soit le type de pixel d'origine
# (Gray{N0f16}, N0f16 brut, etc.)
image_L16 = Float64.(Gray.(img_in)) .* 65535.0

dim_y, dim_x = size(image_L16)
println("dim_x = ", dim_x)
println("dim_y = ", dim_y)
println("max ADU (16 bits) = ", maximum(image_L16))
println("mean ADU (16 bits) = ", @sprintf("%.1f", mean(image_L16)))

# -----------------------------------------------------------------------------------------------
# passage à l'échelle 0-255 (par facilité, voir cahier des charges), calcul de L_dark
# -----------------------------------------------------------------------------------------------
image_L8 = image_L16 ./ const_16_to_8

L_dark = compute_L_dark(image_L8)
println("L_dark (médiane des 10% pixels les plus sombres, échelle 0-255) = ", @sprintf("%.2f", L_dark))
println("tangentes cible (réglables) : R = ", tangent_target_R,
        " ; V = ", tangent_target_V, " ; B = ", tangent_target_B)

# -----------------------------------------------------------------------------------------------
# affichage des 3 courbes ADU_C = f(ADU_L), pour vérification visuelle de l'allure des fonctions
# -----------------------------------------------------------------------------------------------
println("affichage des courbes de colorisation")
x_curve = 0.0:1.0:255.0
y_R = colorize_function.(x_curve, L_dark, target_R, tangent_target_R, pivot, tangent_pivot)
y_V = colorize_function.(x_curve, L_dark, target_V, tangent_target_V, pivot, tangent_pivot)
y_B = colorize_function.(x_curve, L_dark, target_B, tangent_target_B, pivot, tangent_pivot)

plt = plot(x_curve, y_R, label="ADU_R = f(ADU_L)", color=:red, linewidth=2,
           xlabel="ADU_L (échelle 0-255)", ylabel="ADU_C (échelle 0-255)",
           title="Fonctions de colorisation", legend=:topleft)
plot!(plt, x_curve, y_V, label="ADU_V = f(ADU_L)", color=:green, linewidth=2)
plot!(plt, x_curve, y_B, label="ADU_B = f(ADU_L)", color=:blue, linewidth=2)
plot!(plt, x_curve, x_curve, label="identité", color=:gray, linestyle=:dash)
scatter!(plt, [pivot], [pivot], label="pivot (125,125)", color=:black, markersize=5)
scatter!(plt, [L_dark, L_dark, L_dark], [target_R, target_V, target_B],
         label="cible (L_dark)", color=:orange, markersize=5)
display(plt)

# -----------------------------------------------------------------------------------------------
# application des 3 fonctions de colorisation à l'image (échelle 0-255), puis retour à
# l'échelle 16 bits (0-65535) pour l'écriture du fichier de sortie
# -----------------------------------------------------------------------------------------------
println("application des fonctions de colorisation aux canaux R, V, B")
image_R8 = colorize_function.(image_L8, L_dark, target_R, tangent_target_R, pivot, tangent_pivot)
image_V8 = colorize_function.(image_L8, L_dark, target_V, tangent_target_V, pivot, tangent_pivot)
image_B8 = colorize_function.(image_L8, L_dark, target_B, tangent_target_B, pivot, tangent_pivot)

image_R16 = clamp.(image_R8 .* const_16_to_8, 0.0, 65535.0)
image_V16 = clamp.(image_V8 .* const_16_to_8, 0.0, 65535.0)
image_B16 = clamp.(image_B8 .* const_16_to_8, 0.0, 65535.0)

# construction de l'image RGB 16 bits (N0f16 : type normalisé 0.0-1.0 <-> 0-65535)
image_color = RGB.(N0f16.(image_R16 ./ 65535.0),
                    N0f16.(image_V16 ./ 65535.0),
                    N0f16.(image_B16 ./ 65535.0))

# -----------------------------------------------------------------------------------------------
# nom du fichier de sortie : suffixe "_color" ajouté juste avant l'extension (.tif ou .tiff)
# -----------------------------------------------------------------------------------------------
base_name, ext = splitext(file_name_input)
file_name_output = base_name * "_color" * ext

println("écriture du fichier de sortie : ", file_name_output)
save(file_name_output, image_color)

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