Histogramme

Pour étudier les caractéristiques d’une image, rien de mieux qu’une vue des données la composant sous forme d’un histogramme.

Je propose ici une fonction simplifiée (et qui évoluera certainement) au fur et à mesure du temps…

Première version (v1)

BuildHistogram_V1

est une fonction qui va :

  • Permettre d’invoquer la génération d’un histogramme des valeurs (ADU) détectées dans une image (peut importe sa position dans celle-ci), par canal couleur (R,G,B) et de la somme de celles-ci.
  • Sur base d’une image fournie : RAW, TIFF, PNG, JPEG (en se basant sur l’extension)
  • Au maximum de résolution possible ou ADU (Analog Digital Unit) de 16 bits (puisque RAW s’arrête à 14bits) ou de 0 à 65535.
  • En travaillant sur les valeurs RGB interpolées (et PAS les plans d’images non traités, dans cette version)
  • Pour les images RAW (ARW, CR2, CRW à ce stade), on peut préciser la méthode d’interpolation choisie (voir la documentation de demosaic)
  • Va sauver ces valeurs dans
    • Des images : 4 vues ( rouge, vert, bleu et générale)
    • Un fichier txt descriptif
    • L’affichage, pour un résumé rapide.

Appel de la V1 BuildHistogram(fn,path,header,min_value=0,max_value=65535,mode=None,val=16)

Avec à l’entrée
fn : path complet du fichier image à analyser
path : l’endroit où seront stockés les résultats
header : une valeur (texte) qui sera affichée sur toutes les courbes, afin de clairement identifier un cas.
min_value : la valeur minimale à partir de laquelle on veut afficher les courbes. Par défaut 0, mais si on ne vise que les “valeurs significatives”, on peut indiquer 1 par exemple, car sinon les bits “noirs” satureront les mesures affichées (cas d’un offset, par exemple).
max_value : la valeur maximale à partir de laquelle on arrête d’afficher les courbes. Par défaut 65535, mais réduire la valeur permet d’analyser une distribution plus réduite (bits non sur-exposé, par exemple).
mode : méthode de “démosaicage” choisie pour interpréter les plans CFA d’une image RAW, donc transformer RGGB en RGB. Par défaut : AHD (= None).
val : nombre de bits de chaque pixel, par défaut 16.

Et à la sortie :
Chaque appel génèrera (dans la subdir path indiquée) 4 fichiers

3 courbes :
<type> <header><mode> <type histogramme><min_value><max_value>.jpg

avec <type histogramme> :
– Scatter (les trois valeurs en mode “scatter” sur la courbe)
– Red, -Green, – Blue : auto-explicite

1 fichier texte
<type> <header><mode> <type histogramme><min_value><max_value>”-Global.txt”

Ce fichier contient pour toutes les valeurs analysées, le nombre de fois qu’elle est présente au total et par type de couleur… Ex :

Value;Glob;Red;Green;Blue
0;19707989;2348628;9798481;7560880
1;0;0;0;0
2;0;0;0;0
3;0;0;0;0
4;0;0;0;0
5;0;0;0;0
...
65535;0;21;0;0

Ce qui se décode par :
La valeur “0 ADU” se trouve présente dans 197077989 pixels (sur ) et respectivement (RGB) dans 23468628, 9798481,7560880 pixels, soit principalement dans le vert.

De même, la valeur “max” se trouve 21 fois présente dans le vert.

# -*- coding: utf-8 -*-
"""
Created on Sun Aug  5 06:38:22 2018

Build Histogram on RGB values (interpoled or directly delivered in image)

@author: ttf
"""
#☻sample 3

import rawpy
import cv2 as cv
import imageio
import numpy as np
from matplotlib import pyplot as plt

def BuildHistogram(fn,path,header,min_value=0,max_value=65535,mode=None,val=16):
    rgb_red   = []
    rgb_green = []  
    rgb_blue  = []
    rgb_redgreen =[]
    rgb_redgreenblue =[]
    RawSupportedType =["ARW","CR2","CRW"]
    ImgSupportedType =["JPG","JPEG","TIF","TIFF","PNG"]
        
    elem = fn.split("\\")
    imageName=elem[len(elem)-1]
    imageType=imageName.split(".")[len(imageName.split("."))-1].upper()
    
    if imageType in RawSupportedType:
        imageType = "RAW"
        header=header+" Mode "+str(mode)
        print("\nHistogram of "+imageName+" type : "+imageType)
    else:
        if imageType in ImgSupportedType:
            print("\nHistogram of "+imageName+" type : "+imageType)
            header=header+" Mode "+imageType
        else:
            print(imageType+" Not supported")
            imageType ="none"            
    
    print("Decoding image...")    
    
    if imageType == "RAW":
        #decode Raw image
        with rawpy.imread(fn) as raw:
            rgb = raw.postprocess(demosaic_algorithm=mode,output_bps=val)
            # init working arrays
            
    if imageType in ImgSupportedType:
        bgr = cv.imread(fn, -1)
        if bgr is None:
            print ("error reading "+fn+"\n abend")
            return
        rgb = np.zeros((bgr.shape[0],bgr.shape[1],3),dtype="uint16") 
        for i in range(0,bgr.shape[0]):
            for j in range(0,bgr.shape[1]):
                rgb[i,j,0] = int(bgr[i,j,2]*65535)  
                rgb[i,j,1] = int(bgr[i,j,1]*65535)  
                rgb[i,j,2] = int(bgr[i,j,0]*65535)  
    
    print("Image type                        : ",rgb.dtype)
    print("Image size                        : ",rgb.shape)
    print("Min,max pixel value in image      : ",rgb.min(),rgb.max())
    
    max_range=rgb.max()
    for x in range(0,max_range+1):
        rgb_red.append(0)
        rgb_green.append(0)  
        rgb_blue.append(0)
        rgb_redgreen.append(0)
        rgb_redgreenblue.append(0)
    
    
    print("Building histogram...")    
    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  
    
    max_range=rgb.max()
    
    if min_value < rgb.min(): 
        min_value = rgb.min()
    if max_value > max_range: 
        max_value = max_range
    
    print("Min,max pixel value in histogram  : ",min_value,max_value)
    
    #Build plot values 
    plotx_red=[]
    plotx_green=[]
    plotx_blue=[]
    plotx_redgreenblue=[]
    ploty_red=[]
    ploty_green=[]
    ploty_blue=[]
    ploty_redgreenblue=[]
    
    print("Plotting histogram...")    
    print(header+ " Mode " + str(mode))
        
    for v in range(min_value,max_value): 
        rgb_redgreenblue[v]  = rgb_red[v]+rgb_green[v]+rgb_blue[v]  
        if rgb_red[v] > 0:
            plotx_red.append(v)  
            ploty_red.append(rgb_red[v])  
            
        if rgb_green[v] > 0:
            plotx_green.append(v)  
            ploty_green.append(rgb_green[v])  
    
        if rgb_blue[v] > 0:
            plotx_blue.append(v)  
            ploty_blue.append(rgb_blue[v])  
    
        if rgb_redgreenblue[v] > 0:
            plotx_redgreenblue.append(v)  
            ploty_redgreenblue.append(rgb_redgreenblue[v])  
     
    max_y = max(ploty_red[int(min_value):int(max_value)])
    if max(ploty_green[int(min_value):int(max_value)]) > max_y :
        max_y = max(ploty_green[int(min_value):int(max_value)])
    if max(ploty_blue[int(min_value):int(max_value)]) > max_y :
        max_y = max(ploty_blue[int(min_value):int(max_value)])
    
    print("Min red value/pixel in image      : ",min(rgb_red))
    print("Min green value/pixel in image    : ",min(rgb_green))
    print("Min blue value/pixel in image     : ",min(rgb_blue))
    print("Max red value/pixel in image      : ",max(rgb_red))
    print("Max green value/pixel in image    : ",max(rgb_green))
    print("Max blue value/pixel in image     : ",max(rgb_blue))
    
    fig = plt.figure(figsize=(10,5))     
    plt.title(header+" "+imageType+" - Red "+"("+str(min_value)+","+str(max_value)+")")            
    plt.ylim(ymax = max_y, ymin = 0)
    plt.plot(plotx_red[int(min_value):int(max_value)],ploty_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+" "+imageType+" - Green "+"("+str(min_value)+","+str(max_value)+")")            
    plt.ylim(ymax = max_y, ymin = 0)
    plt.plot(plotx_green[int(min_value):int(max_value)],ploty_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+" - Blue "+"("+str(min_value)+","+str(max_value)+")")            
    plt.ylim(ymax = max_y, ymin = 0)
    plt.plot(plotx_blue[int(min_value):int(max_value)],ploty_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+" - Scatter "+"("+str(min_value)+","+str(max_value)+")")            
    plt.ylim(ymax = max(ploty_redgreenblue), ymin = 0)
    x=range(1,max_value)
    plt.scatter(plotx_red[1:max_value],ploty_red[1:max_value],c=['r']) 
    plt.scatter(plotx_green[1:max_value],ploty_green[1:max_value],c=['g']) 
    plt.scatter(plotx_blue[1:max_value],ploty_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+" "+imageType+" 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):
            if x == 0 :
                rgb_redgreenblue[x]=rgb_red[x]+rgb_green[x]+rgb_blue[x]
            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__":
    BuildHistogram("D:\\$python\\PSNR\\test dark 500D_2.CR2","D:\\$python\\PSNR\\","500D (2)",1,65535,rawpy.DemosaicAlgorithm.AAHD)
    

""" 
Some others calls...  
BuildHistogram("D:\\$python\\PSNR\\test dark 500D_3.CR2","D:\\$python\\PSNR\\","500D (3)",1,65535,rawpy.DemosaicAlgorithm.AAHD)
BuildHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.tiff","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,65535,rawpy.DemosaicAlgorithm.AAHD)
   BuildHistogram("D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\Dark30s1600_3879.ARW","D:\\$python\\PSNR\\A7S_Dark30s1600Iso\\","A7S",1,65535,rawpy.DemosaicAlgorithm.AAHD)
"""

Et les images générées…

Deuxième version (v2)

BuildRawHistogram_V1

Est une amélioration de la fonction précédente, qui va désormais traiter les plans d’image RAW pour présenter leur histogramme, avec de nouveaux paramètres de traitement.

Extraire les “bonnes” données RAW ?

Que faut-il extraire et traiter ?
Comment estimer une valeur extraite ?
Quels sont les valeurs fournies par le fabricant ?
Sont-elles justes ? Utiles ?
Quels sont les effets de la compression de RAW sur les images ?
Est-ce que cela peut être estimés ?


Beaucoup de questions, pour afficher une simple courbe… 🙂
On va donc y aller étape par étape… Il y a beaucoup à dire et à écrire…
Ce texte, quasi ce blog, va se publier au fur et à mesure…

Format RAW

Si on veut une solution “toute faite” pour examiner ses RAW’s, pas de problème, allez voir RawDigger qui utilise la librairie LIbRaw qui sera aussi utilisée ici.

Les discussions sur le forum sont à la hauteur de la complexité du sujet, surtout pour les photos “professionnelles”.

Pour les besoins des explications, je vais me baser sur deux types de RAW :
– un type “.CR2” de Canon (2ème version du format RAW)
– un type “.ARW” de Sony (format cRAW (c de compressed))

Pour les exemples, je me baserai sur l’un ou l’autre. Les concepts généraux restent évidemment valables pour l’un ou l’autre.

Une image RAW est dite “brute de capteur”… Il est évident que c’est faux !
Ce que l’on accède a déjà été traité, agencé, modifié par l’appareil et cela avec des algorithmes plus ou moins poussés selon le fabricant et les options activées.

Ex pratique : une image d’un Sony A7S accessible fait 2848 x 4256 pixels. Mais en réalité, l’image RAW sur le capteur fait 2848 x 4288 photosites, et les “rangées” se situant sur le bord du capteur sont soit ignorées, soit utilisées à d’autres fins (on y reviendra).

Si on se met du côté des données APN mises “accessibles” au niveau du fichier (telles que les logiciels standards les utilisent), la structure logique disponible est issue de l’utilisation d’une grille de Bayer et d’un point de vue des données “disponibles” dans le format “capteur”, on aura :

Image RAW 2848×4256

Remarque : Dans une CCD, la surface capteur “utile” = surface “image captée”. Si on annonce 8µm, c’est que chaque détail captable par cette surface via le chemin optique, dans une longueur d’onde donnée, est bien capté et analysé sur cette surface.

Dans un APN sur le modèle Bayer, la capture d’un même détail ponctuel est impossible dans les 3 longueurs d’onde simultanée. La “surface utile” est 4 x plus large (1 R au lieu de 4R ou 1 B au lieu de 4B ou 1G au lieu de 4G pour le même détail ponctuel du ciel). Cette “découpe / distribution” (RGGB) est prévue pour l’œil humain, pas pour les objets célestes, qui se fichent totalement ce cette règle… 🙂

La résolution finale à la sortie, si elle est bien celle annoncée par le fabricant (ici 12M Pix), elle ne sera que “virtuelle”, car issu d’une “estimation” plus ou moins juste. 1 pixel réel “capturé” = 4 pixels physiques (RGGB), donnant 4 pixels logiques “estimés” via fonct_interpolation(R,G,G,B) = (RGB, RGB, RGB, RGB)

Pour construire l’histogramme des valeurs réelles “brutes de capteur”, il faut donc

  • extraire la grille visible (avec toutes les valeurs)
  • extraire la grille complète (avec toutes les valeurs)
  • déduire les valeurs de corrections
  • appliquer la séparer par canal
  • appliquer les corrections nécessaires
  • faire des optimisations (ou pas, selon le cas)
  • et en finale : calculer l’histogramme.

Ensuite, il faudra déterminer l’usage de l’histogramme :

  • Soit pour la photographie “visuelle” (le but est d’avoir une image jolie)
  • Soit pour la science (photométrie, analyse, etc…)

Dans le premier cas, on va s’intéresser à un mélange d’informations : la mesure des bruits, des défauts, de la qualité de l’image, de la qualité de l’interpolation, etc…

Dans le deuxième cas, ce seront les informations les plus “pures” possibles qui seront extraites. Ex : en photométrie, on pourrait ne s’intéresser qu’aux pixels verts (proche du V Johnson).

Extraire les informations

Trois instructions rawpy disponibles

raw_image = raw.raw_image.copy()
=> l’image RAW complète, avec les marges. Dans le cas d’une image Bayer images c’est une ndarray 2D qui est accessible. Rem : si on est dans un mode “4 couleurs”, le 4ème pourrait être à zéro.

raw_visible= raw.raw_image_visible.copy()
=> l’image RAW visible, ou la “complète sans les marges”.

extract_thumb(self)
=> Extraction d’une version “thumbnail” (format JPEG généralement) comme un objet de type rawpy.Thumbnail. Les données pourront être sauvé “as is” dans un fichier.

with rawpy.imread('image.nef') as raw:
  try:
    thumb = raw.extract_thumb()
  except rawpy.LibRawNoThumbnailError:
    print('no thumbnail found')
  except rawpy.LibRawUnsupportedThumbnailError:
    print('unsupported thumbnail')
  else:
    if thumb.format == rawpy.ThumbFormat.JPEG:
      with open('thumb.jpg', 'wb') as f:
        f.write(thumb.data)
    elif thumb.format == rawpy.ThumbFormat.BITMAP:
      imageio.imsave('thumb.tiff', thumb.data)

Déduire les valeurs de correction

L’image RAW “utile” (en fonction du besoin) doit être corrigée en fonction du type de capture, donc, on va récupérer les valeurs utiles via Rawpy…

print ("Raw type                         :",raw.raw_type)
print ("Raw black level                  :",raw.black_level_per_channel)
print ("Raw white balance                :",raw.camera_whitebalance)
print ("Raw color desc                   :",raw.color_desc)
print ("Raw daylight balance             :",raw.daylight_whitebalance)
print ("Raw color number                 :",raw.num_colors)
print ("Raw sizes                        :",raw.sizes)
print ("Raw min, max values ",raw_visible.min(),raw_visible.max())

Black Level : un signal lié à la technologie utilisée (CCD, C-MOS, BSI-C-MOS, etc…). Cette valeur indique le niveau minimal à partir duquel une valeur (convertie en valeur numérique) devient utilisable.

Cette valeur représente un niveau de signal “minimal” émis par le capteur pour fonctionner. Tout capteur mis en activité en génère un, même sans convertir aucun photon.

Ce “black level” :
– est totalement lié à un matériel (APN, capteur, caméra, etc..), que la capture soit statique ou en vidéo
– est soit directement fourni par le fabricant, soit calculé depuis les “pixels masqués” (différence en total et visible) vu qu’ils n’ont pas été exposés aux photons (= signal de “fonctionnement” minimal)
– induit une limite à la “dynamique” disponible par le capteur.

Une dynamique est la “plage de valeurs numériques” que peut prendre un signal (photon, puis e-) convertit en une valeur “discrète”. Le nombre de valeurs “discrètes”, exprimée en ADU (Analog to Digital Unit) est exprimée selon un nombre de bits
– 8 bits : 0 – 255 valeurs (ADU)
– 9 bits : 0 -512 valeurs
– 10 bits : 0 – 1024 valeurs
– 11 bits : 0- 2048 valeurs
– 12 bits : 0- 4096 valeurs
– 13 bits : 0- 8192 valeurs
– 14 bits : 0-16384 valeurs
– 15 bits : 0-32767 valeurs
– 16 bits : 0-65535 valeurs

Et trois limites apparaissent :
– maximale, avec la saturation du photosite une fois passé un certain nombre de photons/e- générés (= val max permanente, en 12 bits = 4096)
– minimale, avec le black level, qui est la limite avant laquelle aucune valeur ne peut être détectée. Donc, qui doit être systématiquement soustrait aux valeurs lues
– l’étendue dynamique (Dynamic Range), qui est l’ensemble des valeurs réellement utilisables.

ex : si 12 bits et BL = 196, il ne reste que 4096-196 = 3900 valeurs pour exprimer réellement le signal numérisé. Et si ce capteur “sature” dès 50000 e- envoyés, avec un taux de conversion (GAIN) de 1 (1phot = 1e-), on a donc : 4096 ADU = 50000e- et 196 ADU = 0e-, 1 ADU = 50000/3900 = ~13 e- (dans lequel se trouvent mêlés tous les signaux parasites et bruits)
Et : un signal inférieur à 12e- n’a aucune chance d’être détecté…

Conversion vers la luminosité…

Outre les informations rendues par Rawpy, on peut extraire des informations par les métadonnées de l’image (EXIF : Exchangeable image file format). On va utiliser ceci dans l’ordre

  • Extraction du “Black level” (vie Rawpy ou EXIF 0xc61a :BlackLevel )
  • Estimation de la moyenne de cette valeur
  • L = imageRaw values – average(blacklevel)
  • Extraction de l’exposition via EXIF 0x829a : ExposureTime

La suite en écriture…