Évaluation de la qualité d’interpolation dans le pré-traitement des RAW

En soi : une image dark « idéale » devrait être noire… Si ce n’est pas le cas, c’est dû :

• Au bruit (inhérent à la capture), qui est lié au capteur (APN, CCD) et aux conditions de la capture.
• Au traitement des « plans » de l’image, notamment la méthode pour « interpoler » des valeurs RAW (brute) afin de reconstruire des valeurs (RGB) conventionnelles.

Rappelons les phases de la capture sur APN (ou capteur similaire) https://fr.wikipedia.org/wiki/Dématriçage

Le fichier RAW est (normalement) issu de la modification minimale des informations issues du capteur. En soi, il devrait être porteur de beaucoup d’informations.
Dans le cas d’une CCD, c’est généralement le cas. La caméra est simplifiée pour « sortir » les données au plus vite et ne passe pas trop de temps à les traiter.
Dans les cas des APN, prévus pour la photo « non astronomique » à la base, le fichier RAW fourni est souvent déjà le résultat d’une première transformation.

Une fois stocké extérieurement au boitier, le fichier RAW est ensuite retraité, via un algorithme qui va réduire la qualité des informations disponibles (en réduisant les plans images) = dématriçage (« debayerisation »)

Il est à remarquer que cet algorithme peut se trouver implémenté à deux endroits :

• Dans le boitier même (le setup format « RAW+JPEG »)
• Dans des programmes extérieurs (fourni par le fournisseur, open-source : DCRAW, LibRaw, Rawtherapee, DarkTable, etc…)

Mais rien ne dit que ce sont les mêmes qui sont disponibles, sauf (peut-être) dans les logiciels fournis par les constructeurs eux-mêmes… Il y a donc une double problématique avec l’astrophoto :

• Chaque image RAW doit être « réduite » (Dark, offset, flat) avant empilement et calcul final. Si chaque dématriçage engendre une erreur, celles-ci vont s’accumuler aux signaux utiles.
• Le format RAW est loin d’être « transparent », il subit de plus en plus de traitement (analyse, suppression d’artefacts de capture, correction, compression,…) avant d’être sauvé sur fichier. Cela va de la simple suppression de « débordement de pixel » (blooming) à la suppression pure et simple d’informations jugées « non significatives ». Donc, le choix de la méthode de dématriçage peut égaler ajouter des erreurs d’interprétations.

Comme montré dans le topic sur la lecture/écriture RAW,

• on dispose des « modes » de dématriçage : LINEAR, VNG, PPG, AHD, DCB, MODIFIED_AHD, AFD, VCD, VCD_MODIFIED_AHD, LMMSE, AMAZE, DHT, AAHD
• les APN actuels (ex: SONY) ont rajouté (dès 2015) une couche logicielle en compressant ou pas, les RAW pour optimiser la taille (en mode supposé LossLess).

Donc.. En considérant que désormais, on ne donne même plus le choix à l’utilisateur (le RAW est de-facto compressé…) et que l’on rajoute à cela le choix possible de méthode de dématriçage, quel est le meilleur setup pour l’astrophoto ?

En gros, ces mécanismes de dématriçage sont de deux sortes : adaptatif ou non-adaptatif.
Les non-adaptatifs sont les plus simples et plus rapides, car ils considèrent chaque « plan » image indépendamment et ne réalisent la fusion qu’à l’étape finale (via une « moyenne ») sur les pixels existants.
Les adaptatifs vont examiner tous les plans images simultanément pour en détecter des caractéristiques (seuil, limites, etc…) et « adapter » chaque plan à l’étape de fusion finale. Ainsi, certains éléments (contours d’un visage, par exemple) seront renforcés, d’autres supprimés (trame, défauts, etc…).
Ces algorithmes sont adaptés aux images et courleurs « humainement visibles » (Human Vision System).

Mais dans l’astrophotographie, la grille de Bayer fait déjà perdre une grande quantité d’information utile (entre deux pixels : il y un « trou » d’informations). C’est pour cela que prendre des images monochromes est plus performant et précis que d’utiliser des captures couleur native.

Mais le prix et la sensibilité des capteur CMOS (et Full Frame) actuels renforce l’usage des boitiers APN pour l’astrophoto, donc : on va chercher à minimiser les défauts…

Le sujet fait polémique depuis aussi longtemps que les capteurs numériques couleurs existent… Aujourd’hui, les méthodes suivantes sont les plus utilisées :

• AHD : Adaptive Homogeneity-Directed , l’interpolation par défaut de DCRAW (Adaptative, la plus utilisée)
  Cette méthode détermine une « direction » pour l’interpolation et fera ses estimations pour les pixels manquant sur cette base.
• VNG : Variable Number of Gradients (Adaptative, assez ancienne).
  La méthode calcule des gradients autour de pixels considérés intéressants et utilise le résultat pour estimer la couleur finale.
• VNG 4C : VNG « four color », utilisée pour supprimer des artefacts de couleur (Adaptative, plus récente)
• PPG : Patterned Pixel Grouping, aussi meilleur que les autres, mais plus rapide (Adaptative)
  Travaille sur la « détection » de scènes naturelles pour estimer les valeurs manquantes. Elle produit moins d’artefacts de couleur que la méthode VNG.
• Bilinear : la plus simple, la plus ancienne, pour ordinateur lent…
  La valeur d’une couleur manquante est construite sur la base de la moyenne des valeurs adjacentes (Non adaptative)
• DHT : Distritubed Hash Table, une méthode de « hashing » générale qui peut se trouver implémentée via plusieurs algorithmes (Chord, Pastry, etc…) et donc varier par librairie. Comme on utilise libraw ici, cela sera donc identique pour beaucoup de logiciels.

Pour l’astronomie… Un pixel est un pixel… Sa valeur peut être qualitative (une planète, une étoile ?) et toute information est la bienvenue…
De plus, aucune des scènes traitées ne peut être considérée comme « naturelle », ce qui peut induire des interprétations fautives (un nuage dans le ciel n’est pas un nuage de nébuleuses).

Comment évaluer cela ?

De (très) nombreux articles existent, un exemple : http://www.ruevski.com/rawhistogram/40D_Demosaicing/40D_DemosaicingArtifacts.html
ou
https://www.libraw.org/node/2306

Et certains nouveaux algorithmes ont perturbé la détection d’objets…
https://petapixel.com/2018/06/08/sonys-star-eater-problem-has-been-defeated-in-the-a7-iii/

Donc : inutiles de tout faire en détail si c’est pour avoir des artefacts inutiles à la fin… Mais on ne peut pas passer son temps à vérifier chaque image. J’ai donc voulu tester plusieurs programmes pour avoir une « idée » de ce que l’on pouvait obtenir…

Via un histogramme « adapté » ?

Très facile d’observer les différences entre modèles d’interpolation…

A l'entrée : fichiers issus d'un Sony A7S
Zero_A7S_2848x4256.tif = une image (0,0,0) de taille équivalente au RAW

Dark_1600iso_1s.ARW
Dark_1600iso_2s.ARW
Dark_1600iso_4s.ARW
Dark_1600iso_8s.ARW
Dark_1600iso_10s.ARW
Dark_1600iso_20s.ARW

Dark30s1600_3879.ARW
Dark30s1600_3880.ARW
Dark30s1600_3881.ARW
Dark30s1600_3882.ARW
Dark30s1600_3883.ARW
= des captures de "dark", à ISO constant, format RAW en mode "standard" 

Flat_12800iso.ARW
Flat_6400iso.ARW
Flat_3200iso.ARW
Flat_1600iso.ARW
= des captures de "float", à ISO variable, format RAW en mode "standard"

A l'entrée : fichiers issus d'un Canon 550D
Dark30s1600_0055.CR2
Dark30s1600_0056.CR2
Dark30s1600_0057.CR2
Dark30s1600_0058.CR2
Dark30s1600_0059.CR2
= des captures de "dark", à ISO constant, format RAW en mode "standard" 

Prenons les algorithmes d’interpolation AHD (par défaut dans de nombreux logiciels) et DHT (mis en avant par certaines librairies) appliqués sur la même image et affichons les statistiques des deux types de RAW (Canon et Sony, même exposition et ISO).

Une petite fonction d’histogramme « personnalisé » devrait aider pour la suite… On compte les valeurs obtenues pour chaque pixel, on en fait une table de ces valeurs et on affiche.

import rawpy
import imageio
from matplotlib import pyplot as plt

def BuildRawHistogram(fn,path,header,min_value=0,max_value=255,mode=None):
    rgb_red   = []
    rgb_green = []  
    rgb_blue  = []
    rgb_glob  = []
    with rawpy.imread(fn) as raw:
        rgb = raw.postprocess(demosaic_algorithm=mode)
        # init working arrays
               
        for x in range(0,rgb.max()+1):
            rgb_red.append(0)
            rgb_green.append(0)  
            rgb_blue.append(0)
            rgb_glob.append(0)
        
        if min_value < rgb.min(): 
            min_value = rgb.min()
        if max_value > rgb.max(): 
            max_value = rgb.max()
       
        print("Min pixel value in image : ",rgb.min())
        print("Max pixel value in image : ",rgb.max())
        print("Min pixel value in histogram : ",min_value)
        print("Max pixel value in histogram : ",max_value)
        
        print(header+" Mode "+str(mode))
        for i in range(0,rgb.shape[0]):
            for j in range(0,rgb.shape[1]):
                rgb_red[rgb[i,j,0]] +=1  
                rgb_green[rgb[i,j,1]] +=1  
                rgb_blue[rgb[i,j,2]] +=1  
                rgb_glob[rgb[i,j,0]] +=1  
                rgb_glob[rgb[i,j,1]] +=1  
                rgb_glob[rgb[i,j,2]] +=1  

        print("Min pixel count by value in image : ",min(rgb_glob))
        print("Max pixel count by value in image : ",max(rgb_glob))
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Red "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_red[int(min_value):int(max_value)])            
        fig.savefig(path+header+" Mode "+str(mode)+"- Red ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()

        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Green "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_green[int(min_value):int(max_value)])            
        fig.savefig(path+header+" Mode "+str(mode)+"- Green ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Blue "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_blue[int(min_value):int(max_value)])            
        fig.savefig(path+header+" Mode "+str(mode)+"- Blue ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Global "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_glob[int(min_value):int(max_value)])            
        fig.savefig(path+header+" Mode "+str(mode)+"- Global ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
    return     

if __name__ == "__main__":
    BuildRawHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,255,rawpy.DemosaicAlgorithm.AAHD)
    BuildRawHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,255,rawpy.DemosaicAlgorithm.DHT)
    BuildRawHistogram("D:\\$python\\PSNR\\550D_Dark30sec1600iso\\Dark30s1600_0055.CR2","D:\\$python\\PSNR\\550D_Dark30sec1600iso\\","550D",1,255,rawpy.DemosaicAlgorithm.AAHD)
    BuildRawHistogram("D:\\$python\\PSNR\\550D_Dark30sec1600iso\\Dark30s1600_0055.CR2","D:\\$python\\PSNR\\550D_Dark30sec1600iso\\","550D",1,255,rawpy.DemosaicAlgorithm.DHT)

Et on obtient déjà des résultats intéressants…

Min pixel value in image :  0
Max pixel value in image :  255
Min pixel value in histogram :  1
Max pixel value in histogram :  255
A7S Mode DemosaicAlgorithm.AAHD
Min pixel count by value in image :  0
Max pixel count by value in image :  15133161

Min pixel value in image :  0
Max pixel value in image :  255
Min pixel value in histogram :  1
Max pixel value in histogram :  255
A7S Mode DemosaicAlgorithm.DHT
Min pixel count by value in image :  0
Max pixel count by value in image :  18104827

Min pixel value in image :  0
Max pixel value in image :  255
Min pixel value in histogram :  1
Max pixel value in histogram :  255
550D Mode DemosaicAlgorithm.AAHD
Min pixel count by value in image :  0
Max pixel count by value in image :  21416408

Min pixel value in image :  0
Max pixel value in image :  255
Min pixel value in histogram :  1
Max pixel value in histogram :  255
550D Mode DemosaicAlgorithm.DHT
Min pixel count by value in image :  2648
Max pixel count by value in image :  26837392

Déjà à ce stade, le niveau et les valeurs varient clairement selon la méthode choisie. Avec les graphiques générés, on voit clairement l’impact au niveau global (addition des trois canaux par valeur) et les influences des boitiers !

Dire que le choix de l’algorithme d’interpolation n’a pas d’impact serait largement incorrect… Mais ici, cela ne renseigne que sur le « bruit » présent (que l’on devrait normalement supprimer lors de la réduction).

Apparemment,

• le DHT augmente le nombre de pixels avec des valeurs « basses » et minimize le nombre de valeurs élevées.
• le AAHD distribue plus les valeurs sur la dynamique complète, ce qui augmente le niveau global.
• Si on examine par canal, la différence est flagrante dans le cas du rouge :

Par contre, si on distingue une « cohérence » (courbe assez semblable) dans le cas des RAW Canon avec des modes différents, cependant, cette tendance ne se retrouve pas du tout dans les RAW de Sony.

De plus, une bizarrerie qui doit trouver une explication logique : les « trous » de valeurs dans l’image. Si on fait un « zoom » sur une zone (par exemple de 10 à 50), on trouve :

Pourquoi ces « trous » dans les valeurs dans le cas A7S ? Alors que c’est largement moins fréquent dans le cas 550D…

Maintenant, reste à déterminer quel est la manière optimale de pré-traiter une image dans toute la chaîne.

Dans la documentation de Rawtherapee, il y a quelques exemples de traitements selon les modes supportés par le produit :

Pour la photo courante, ils conseillent d’utiliser AMaZE.

Oui, mais les RAW, c’est 16 bits !

Artefact de compression vers 8bits, me dira-t-on…. Alors que le « vrai » RAW est en 10,12 ou 14 bits ! Soit, il suffit d’adapter quelques lignes dans le code, et on va améliorer les graphiques fournis + garder les données pour disposer d’un « profil » par image….

import rawpy
import imageio
from matplotlib import pyplot as plt

def BuildRawHistogram(fn,path,header,min_value=0,max_value=65535,mode=None,val=16):
    rgb_red   = []
    rgb_green = []  
    rgb_blue  = []
    rgb_redgreen =[]
    rgb_redgreenblue =[]
     
    with rawpy.imread(fn) as raw:
        rgb = raw.postprocess(demosaic_algorithm=mode,output_bps=val)
        # init working arrays
           
        for x in range(0,rgb.max()+1):
            rgb_red.append(0)
            rgb_green.append(0)  
            rgb_blue.append(0)
            rgb_redgreen.append(0)
            rgb_redgreenblue.append(0)
        
        if min_value < rgb.min(): 
            min_value = rgb.min()
        if max_value > rgb.max(): 
            max_value = rgb.max()
           
        print("Min pixel value in image : ",rgb.min())
        print("Max pixel value in image : ",rgb.max())
        print("Min pixel value in histogram : ",min_value)
        print("Max pixel value in histogram : ",max_value)
        
        print(header+" Mode "+str(mode))
        for i in range(0,rgb.shape[0]):
            for j in range(0,rgb.shape[1]):
                rgb_red[rgb[i,j,0]]   +=1  
                rgb_green[rgb[i,j,1]] +=1  
                rgb_blue[rgb[i,j,2]]  +=1  
    
        for v in range(0,rgb.max()+1): 
                rgb_redgreen[v]      = rgb_red[v]+rgb_green[v]  
                rgb_redgreenblue[v]  = rgb_red[v]+rgb_green[v]+rgb_blue[v]  
        
        print("Min pixel count by value in image : ",min(rgb_redgreenblue))
        print("Max pixel count by value in image : ",max(rgb_redgreenblue))
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Red "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_red[int(min_value):int(max_value)],color='r')            
        fig.savefig(path+header+" Mode "+str(mode)+"- Red ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Green "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_green[int(min_value):int(max_value)],color='g')            
        fig.savefig(path+header+" Mode "+str(mode)+"- Green ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Blue "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_blue[int(min_value):int(max_value)],color='b')            
        fig.savefig(path+header+" Mode "+str(mode)+"- Blue ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
    
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Global "+"("+str(min_value)+","+str(max_value)+")")            
        plt.plot(rgb_redgreenblue[int(min_value):int(max_value)],color='k')            
        fig.savefig(path+header+" Mode "+str(mode)+"- Global ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
        
        fig = plt.figure(figsize=(10,5))     
        plt.title(header+" Mode "+str(mode)+"- Scatter "+"("+str(min_value)+","+str(max_value)+")")            
        x=range(1,max_value)
        plt.scatter(x,rgb_red[1:max_value],c=['r']) 
        plt.scatter(x,rgb_green[1:max_value],c=['g']) 
        plt.scatter(x,rgb_blue[1:max_value],c=['b']) 
        fig.savefig(path+header+" Mode "+str(mode)+"- Scatter ("+str(min_value)+","+str(max_value)+").jpg", bbox_inches='tight', dpi=150)
        plt.show()
               
        
        with open(path+header+" Mode "+str(mode)+"- Global ("+str(min_value)+","+str(max_value)+").txt", 'w') as f:
            f.write("Value;Glob;Red;Green;Blue\n")
            for x in range(0,rgb.max()+1):
                line=str(x) +";"+ str(rgb_redgreenblue[x]) +";"+ str(rgb_red[x]) +";"+ str(rgb_green[x]) +";"+ str(rgb_blue[x]) + "\n"                
                f.writelines(line)
        f.close()   
    return rgb    

if __name__ == "__main__":
    BuildRawHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,65535,rawpy.DemosaicAlgorithm.AAHD)
    BuildRawHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,65535,rawpy.DemosaicAlgorithm.DHT)

Et le résultat est tout aussi éloquent…

A regarder ceci globalement…
Dans le cas d’un « dark », le « bruit » résultat en DHT est plus faible qu’en AHD…

Peut-on estimer la qualité via PSNR ?

Par curiosité, j’ai eu l’idée d’appliquer le PSNR (chapitre précédant) sur plusieurs comparaisons, en considérant (peut-être faussement) que :

• L’image « signal idéal » était une image créée artificiellement avec tous les pixels à 0… Le « noir » absolu.
• L’image « bruit réel » serait fournie par l’interpolation d’une capture au format RAW selon diverses méthodes.

Le résultat pourrait indiquer quoi ?
– on est certes en-dehors du cadre de PSNR, car aucune compression n’a été faite.
– mais le processus d’interpolation est, en soit, un traitement mathématique de l’image dans lequel on « déforme » la réalité
– on rajoute de plus en plus de traitement sur le « RAW » (compression, amélioration, suppression de bruit), qu’il finit par ne plus trop l’être…
Si on les traite pour « reconstruire » une réalité subjective, on n’est pas loin des techniques de compression/décompression.

Peut-on analyser les images de manière « automatisée » ?
Tout d’abord, créer deux images « signal idéal » avec RGB (0,0,0) équivalentes en taille aux images RAW respectives. Puis les appliquer au calcul du PSNR (signal = noir idéal, noise = image réelle) pour obtenir les valeurs caractéristiques.

import os
import cv2
import math
import rawpy
import numpy as np

def ComputePSNR(signal, noise):
    #Peak signal-to-noise ratio, PSNR, is the ratio between 
    #the maximum possible power of a signal and the power of corrupting noise 
    #that affects the fidelity of its representation. 
    #Because many signals have a very wide dynamic range, PSNR is usually expressed 
    #in terms of the logarithmic decibel scale.
    # Applied to 3 color planes (BGR Mode), compute globally and for each 
    D = np.array(noise - signal, dtype=np.int64)
    D[:,:,:] = D[:,:,:]**2
    #we assume 16bits stored image with a mex 14 bits RAW
    R= float(16384.**2)
    #compute for global, then for each plane RMSE and PSNR
    RMSE = D[:,:,:].sum()/signal.size
    psnr = 10*math.log10(R/RMSE)

    RMSE = D[:,:,2].sum()/signal.size
    psnrr = 10*math.log10(R/RMSE)

    RMSE = D[:,:,1].sum()/signal.size
    psnrg = 10*math.log10(R/RMSE)

    RMSE = D[:,:,0].sum()/signal.size
    psnrb = 10*math.log10(R/RMSE)

    print("PSNR :",psnr)
    print("PSNR R,G,B :",psnrr,psnrg,psnrb)
    return psnr,psnrr,psnrg,psnrb

def ReadRawImage(fn,demosaic=None,export=False):
    #Read a raw DSLR image, interpolate it and provide an RGB array
    #Possible demosaic are : 
    #None
    #rawpy.DemosaicAlgorithm.AHD
    #rawpy.DemosaicAlgorithm.AAHD
    #rawpy.DemosaicAlgorithm.DCB
    #rawpy.DemosaicAlgorithm.DHT
    #rawpy.DemosaicAlgorithm.LINEAR
    #rawpy.DemosaicAlgorithm.PPG
    #rawpy.DemosaicAlgorithm.VNG
        
    with rawpy.imread(fn) as raw: 
        raw_image = raw.raw_image.copy() 
        #interpolate CFA array to provide an RGB version
        #keep it as most linear as possible
        rgb = raw.postprocess(no_auto_bright=True,use_auto_wb =False,gamma=None,output_bps=16,demosaic_algorithm=demosaic)
        #extract filename and maximum value in image
        path, filename = os.path.split(fn)
        maxvalue = raw_image.max()
        if (export == True):
            imageio.imsave(fn+"_"+str(demosaic)+".TIFF", rgb)
            print(fn+"_"+str(demosaic)+".TIFF exported")
    print(fn+" size = "+str(rgb.shape))
    return (filename,maxvalue,rgb)

def ReadStdImage(fn):
    #Read a standard (jpg, png, tiff) image and provide an RGB array
    bgr = cv2.imread(fn)
    path, filename = os.path.split(fn)
    maxvalue = bgr.max()
    rgb = bgr 
    rgb[:,:,0] = bgr[:,:,2] #(bgr)R > (rgb)R
    rgb[:,:,1] = bgr[:,:,1] #(bgr)G > (rgb)G
    rgb[:,:,2] = bgr[:,:,0] #(bgr)B > (rgb)B
    print(fn+" size = "+str(rgb.shape))
    return (filename,maxvalue,rgb)

if __name__ == "__main__":
    # test avec algorithme rawpy.DemosaicAlgorithm.AHD 
    print("---- A7S ----")
    fn1,max1,rgb1 = ReadStdImage("D:\\$python\\PSNR\\Zero_A7S_2848x4256.tif")
    fn2,max2,rgb2 = ReadRawImage("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW",demosaic=rawpy.DemosaicAlgorithm.DHT)
    print("S:"+fn1+"\t N:"+fn2)
    print("S Max:"+str(max1)+"\t N Max :"+str(max2))
    psnr = ComputePSNR(rgb1, rgb2)

    print("---- 550D ----")
    fn1,max1,rgb1 = ReadStdImage("D:\\$python\\PSNR\\Zero_550D_5202x3465.tif")
    fn2,max2,rgb2 = ReadRawImage("D:\\$python\\PSNR\\550D_Dark30sec1600iso\\Dark30s1600_0055.CR2",demosaic=rawpy.DemosaicAlgorithm.DHT)
    print("S:"+fn1+"\t N:"+fn2)
    print("S Max:"+str(max1)+"\t N Max :"+str(max2))
    psnr = ComputePSNR(rgb1, rgb2)

Résultat: pour DHT
---- A7S ----
D:\$python\PSNR\Zero_A7S_2848x4256.tif size = (2848, 4256, 3)
D:\$python\PSNR\A7S_Dark30s1600Iso\Dark30s1600_3879.ARW size = (2848, 4256, 3)
S:Zero_A7S_2848x4256.tif         N:Dark30s1600_3879.ARW
S Max:0  N Max :5352
PSNR : 29.031744905541622
PSNR R,G,B : 36.61646701454189 39.32142371292062 30.386278011119288
---- 550D ----
D:\$python\PSNR\Zero_550D_5202x3465.tif size = (3465, 5202, 3)
D:\$python\PSNR\550D_Dark30sec1600iso\Dark30s1600_0055.CR2 size = (3465, 5202, 3)
S:Zero_550D_5202x3465.tif        N:Dark30s1600_0055.CR2
S Max:33         N Max :16382
PSNR : 17.539683692459473
PSNR R,G,B : 23.114462510715946 24.143550465319755 20.51200586227327

En exécutant cela pour diverses combinaisons, on obtient (en dB):

APN	AHD	AAHD	DCB	DHT	PPG	VNG	LINEAR
A7S	30.53	32.69	30.45	29.01	30.93	31.73	30.42
550D	20.04	20.41	19.99	17.53	20.06	20.26	19.98

Si on se baserai sur la valeur typique du PSNR, la meilleure reconstruction d’un signal « noir » idéal se trouve

• On regarde que les valeurs sont meilleures pour le A7S, ce qui pourrait indiquer que globalement le « bruit » est moindre (puisque « plus proche » du noir = signal recherché)
• Le mode AAHD fournit le meilleur PSNR, et DHT le moins bon
• La différence entre minimum et maximum PSNR est de 3.68 dB (A7S) et 2.88 dB (550D), ce qui est relativement important

Rappel : le décibel (dB) est une unité définie comme dix fois le logarithme décimal du rapport entre deux puissances.

Donc, le rapport de « qualité » entre deux images générées à partir de la même source varie entre 2 et 2.5x. C’est l’inverse de l’analyse intuitive sur base des graphes à l’étape précédente.

Que se passe-t-il à l’empilage ?

L’additions d’images peut se faire de plusieurs manières :

• Soit sur base des plans images (CFA) extraits de l’image RAW.
  Et en finale, une interpolation donnera l’image couleur finale.
• Soit sur base d’images couleurs issues de l’interpolation RAW et stockées dans des formats intermédiaires (TIFF, etc..)
• Soit avec un mélange des deux !

Évidemment, les effets cumulés de l’interpolation et de l’addition peut générer des résultats différents.

La forêt des paramètres

Une particularité du Raw est le nombre de paramètres disponibles pour obtenir la qualité d’image désirée…

Revoyons la liste :

• demosaic_algorithm (rawpy.DemosaicAlgorithm) – default is AHD
• half_size (bool) – outputs image in half size by reducing each 2×2 block to one pixel instead of interpolating
• four_color_rgb (bool) – whether to use separate interpolations for two green channels
• dcb_iterations (int) – number of DCB correction passes, requires DCB demosaicing algorithm
• dcb_enhance (bool) – DCB interpolation with enhanced interpolated colors
• fbdd_noise_reduction (rawpy.FBDDNoiseReductionMode) – controls FBDD noise reduction before demosaicing
• noise_thr (float) – threshold for wavelet denoising (default disabled)
• median_filter_passes (int) – number of median filter passes after demosaicing to reduce color artifacts
• use_camera_wb (bool) – whether to use the as-shot white balance values
• use_auto_wb (bool) – whether to try automatically calculating the white balance
• user_wb (list) – list of length 4 with white balance multipliers for each color
• output_color (rawpy.ColorSpace) – output color space
• output_bps (int) – 8 or 16
• user_flip (int) – 0=none, 3=180, 5=90CCW, 6=90CW, default is to use image orientation from the RAW image if available
• user_black (int) – custom black level
• user_sat (int) – saturation adjustment (custom white level)
• no_auto_scale (bool) – Whether to disable pixel value scaling
• no_auto_bright (bool) – whether to disable automatic increase of brightness
• auto_bright_thr (float) – ratio of clipped pixels when automatic brighness increase is used (see no_auto_bright). Default is 0.01 (1%).
• adjust_maximum_thr (float) – see libraw docs
• bright (float) – brightness scaling
• highlight_mode (rawpy.HighlightMode | int) – highlight mode
• exp_shift (float) – exposure shift in linear scale. Usable range from 0.25 (2-stop darken) to 8.0 (3-stop lighter).
• exp_preserve_highlights (float) – preserve highlights when lightening the image with exp_shift. From 0.0 to 1.0 (full preservation).
• gamma (tuple) – pair (power,slope), default is (2.222, 4.5) for rec. BT.709
• chromatic_aberration (tuple) – pair (red_scale, blue_scale), default is (1,1), corrects chromatic aberration by scaling the red and blue channels
• bad_pixels_path (str) – path to dcraw bad pixels file. Each bad pixel will be corrected using the mean of the neighbor pixels. See the rawpy.enhance module for alternative repair algorithms, e.g. using the median.

Quels sont les effets exacts de chacun ?
Pour le déterminer plus précisément, car la documentation est plutôt limitée… Il suffit de comparer via programme !
Prenons une image prise par un APN et sauvée en format RAW et JPEG.

Le contenu JPEG

Le contenu RAW vers JPEG (8bits) via Darktable

Si on compare, le résultat est très proche… Et confirme simplement que le logiciel de l’appareil photo fait bien son boulot…

Donc, au minimum, le pre-processing vers RAW exécuté par le code doit se rapprocher de ceci… Voyons les effets liés par les options disponibles.

import cv2
import numpy as np
#default loop mode
loopstack=True

def writeimage(name,image):
    imagebgr = cv2.cvtColor(image,cv2.COLOR_RGB2BGR).astype(np.float32) / 65535
    output_image = (65535*imagebgr).astype(np.uint16)
    outputfile = "G:\$Python\livestack\output\\" + name + ".png"
    cv2.imwrite(outputfile,output_image)
    print(name,"saved")
    return    

# ===== MAIN =====
if __name__ == '__main__':
    fileraw    = "G:\$Python\livestack\input\_DSC0553.ARW"
    darkraw    = "G:\$Python\livestack\dark\_DSC0583.ARW"

    filejpeg   = "G:\$Python\livestack\jpeg\input\_DSC0553.JPG"
    darkjpeg    = "G:\$Python\livestack\jpeg\dark\_DSC0583.JPG"
    
    image = cv2.imread(filejpeg,1).astype(np.float32) / 255
    output_image = (255*image).astype(np.uint8)
    outputfile = "G:\$Python\livestack\output\\basejpeg.png"
    cv2.imwrite(outputfile,output_image)

    with rawpy.imread(fileraw) as raw:
        imagergb = raw.postprocess(output_bps=8)             
        imagebgr = cv2.cvtColor(imagergb,cv2.COLOR_RGB2BGR).astype(np.float32) / 255
        output_image = (255*imagebgr).astype(np.uint8)
        outputfile = "G:\$Python\livestack\output\\baseraw.png"
        cv2.imwrite(outputfile,output_image)

    with rawpy.imread(fileraw) as raw:

       imagergb = raw.postprocess(no_auto_bright=True, output_bps=16)             
       writeimage("bps16_no_auto_bright",imagergb)
       
       imagergb = raw.postprocess(auto_bright_thr=0.5, output_bps=16)             
       writeimage("bps16_auto_bright_thr0.5",imagergb)
       
       imagergb = raw.postprocess(use_camera_wb=True , output_bps=16)             
       writeimage("bps16_use_camera_wb ",imagergb)
        
       maxsat= round(np.max(imagergb)/2,0)

       imagergb = raw.postprocess(user_sat=maxsat,output_bps=16)             
       writeimage("bps16_maxsat",imagergb)

       imagergb = raw.postprocess(output_bps=16)             
       writeimage("bps16",imagergb)
 
       imagergb = raw.postprocess(half_size=True, output_bps=16)             
       writeimage("bps16_half_size",imagergb)
       
       imagergb = raw.postprocess(four_color_rgb=True, output_bps=16)             
       writeimage("bps16_four_color_rgb",imagergb)

       imagergb = raw.postprocess(dcb_iterations=5,demosaic_algorithm=rawpy.DemosaicAlgorithm.DCB, output_bps=16)             
       writeimage("bps16_dcb_iterations",imagergb)

       imagergb = raw.postprocess(dcb_enhance=True,demosaic_algorithm=rawpy.DemosaicAlgorithm.DCB, output_bps=16)             
       writeimage("bps16_dcb_enhance",imagergb)

       imagergb = raw.postprocess(fbdd_noise_reduction=rawpy.FBDDNoiseReductionMode.Full, output_bps=16)             
       writeimage("bps16_fbdd_noise_reduction",imagergb)

       imagergb = raw.postprocess(noise_thr=0.5, output_bps=16)             
       writeimage("bps16_noise_thr",imagergb)

       imagergb = raw.postprocess(median_filter_passes=5, output_bps=16)             
       writeimage("bps16_median_filter_passes",imagergb)
       
       
       imagergb = raw.postprocess(use_auto_wb=True, output_bps=16)             
       writeimage("bps16_use_auto_wb",imagergb)
       
       imagergb = raw.postprocess(user_wb=(1,1,1,1), output_bps=16)             
       writeimage("bps16_user_wb",imagergb)
       
       imagergb = raw.postprocess(output_color=rawpy.ColorSpace.raw, output_bps=16)             
       writeimage("bps16_output_colorraw",imagergb)
       imagergb = raw.postprocess(output_color=rawpy.ColorSpace.sRGB, output_bps=16)             
       writeimage("bps16_output_colorsrgb",imagergb)
       imagergb = raw.postprocess(output_color=rawpy.ColorSpace.Adobe, output_bps=16)             
       writeimage("bps16_output_coloradobe",imagergb)
       imagergb = raw.postprocess(output_color=rawpy.ColorSpace.Wide, output_bps=16)             
       writeimage("bps16_output_colorwide",imagergb)
       
       imagergb = raw.postprocess(user_black=0, output_bps=16)             
       writeimage("bps16_user_black",imagergb)
       
       imagergb = raw.postprocess(user_sat=32767, output_bps=16)             
       writeimage("bps16_user_sat32767",imagergb)
       
       imagergb = raw.postprocess(no_auto_scale=True, output_bps=16)             
       writeimage("bps16_no_auto_scale",imagergb)
       
       imagergb = raw.postprocess(adjust_maximum_thr=0.02, output_bps=16)             
       writeimage("bps16_adjust_maximum_thr0.02",imagergb)
       
       imagergb = raw.postprocess(bright=0.5, output_bps=16)             
       writeimage("bps16_bright0.5",imagergb)
       
       imagergb = raw.postprocess(highlight_mode=0, output_bps=16)             
       writeimage("bps16_highlight_mode0",imagergb)
       imagergb = raw.postprocess(highlight_mode=1, output_bps=16)             
       writeimage("bps16_highlight_mode1",imagergb)
       imagergb = raw.postprocess(highlight_mode=2, output_bps=16)             
       writeimage("bps16_highlight_mode2",imagergb)
       
       imagergb = raw.postprocess(exp_shift=0.25, output_bps=16)             
       writeimage("bps16_exp_shift0.25",imagergb)
       
       imagergb = raw.postprocess(exp_preserve_highlights=0.5, output_bps=16)             
       writeimage("bps16_exp_preserve_highlights0.5",imagergb)
       
       imagergb = raw.postprocess(gamma=(2.2,4.5),output_bps=16)             
       writeimage("bps16_gamm2.2_4.5",imagergb)
       imagergb = raw.postprocess(gamma=(1,1),output_bps=16)             
       writeimage("bps16_gamm1_1",imagergb)

       imagergb = raw.postprocess(chromatic_aberration=(0.2,0.5), output_bps=16)             
       writeimage("bps16_chromatic_aberration0.2_0.5",imagergb)
       
       imagergb = raw.postprocess(bad_pixels_path="G:\$Python\livestack\output\\badpixels.txt", output_bps=16)             
       writeimage("bps16_",imagergb)

       imagergb = raw.postprocess(no_auto_bright=True,use_camera_wb=True,user_sat=maxsat,output_bps=16)             
       writeimage("bps16_target ",imagergb) 

Ce programme n’a qu’un seul but : montrer les différences…
Et le moins que l’on puisse dire, c’est que c’est visible…

Pour être équivalent au pré-traitement automatique de Darktable, il faudra combiner plusieurs paramètres et valeurs.

Une première hypothèse, dans le cas de cette image, peut par exemple être de (en considérant que maxsat= round(np.max(imagergb)/2,0) ) :

Clairement, il faudra creuser le sujet pour obtenir un résultat optimal (donc analyser l’image). Le phénomène de « Hightlighting » qui est appliqué doit être contrôlé pour obtenir la valeur optimale.