Colorisation de l'Empire Russe

Sergei Mikhailovich Prokudin-Gorskii (1863-1944)

Alignement et restauration d'images du siècle dernier.

L'Empire Russe

Au début du siècle dernier, le chimiste et photographe Sergei Prokudin-Gorskii fut mandaté par le Tsar Nicholas II pour prendre des photos couleur de l'Empire Russe. Le procédé employé par Prokudin-Gorskii capturait trois images différentes, soit une par teinte de couleur primaire.

Projet

Le but du projet est de recréer les images couleurs à partir des trois captures de Prokudin-Gorskii. Ces dernières ont été numérisées par la Librairie du Congrès aux États-Unis en 2004 et sont disponibles publiquement.

Une fois les trois images représentant chacune des couleurs primaires extraites, il suffit de les aligner pour en composer l'image couleur. Diverses améliorations seront ensuite discutées par la suite.

Code

Le code est disponible dans l'archive au même niveau que cette page. Le code est disponible en format Python (.py), IPython Notebook (.ipynb) et HTML (q1 et q2, reste)


Question 1

Afin d'aligner les images représentant chacun des canaux, une métrique de comparaison des images est nécessaire. Pour ce travail, la somme des différences aux carré (SDC) a été employée:

$$\tag{1} \epsilon_{\hat{Y},Y} = \sum_{i=1}^{n}(\hat{Y}_i - Y_i)^2 \quad ,$$

où $ n $ est la quantité de pixels dans l'image, $ Y $ et $ \hat{Y} $ sont les deux images comparées et $ \epsilon_{\hat{Y}Y} $ la somme des différences au carré.

Afin d'améliorer la performance de l'algorithme, la comparaison des images se fait sur leurs arêtes, correspondant aux hautes fréquences de l'images. Le filtre appliqué à l'image $ f $ pour obtenir ses hautes fréquences est le gradient $ \nabla f $, obtenu en calculant: $$\tag{2} \nabla f = \frac{\partial f}{\partial x} \hat{x} + \frac{\partial f}{\partial y} \hat{y} \quad .$$

Pour trouver l'alignement entre deux images, celles-ci sont superposées puis comparées à l'aide de la métrique SDC. Le processus est répété tout en décalant la seconde image d'une plage de ±15 pixels selon les deux axes. Le décalage ayant obtenu la plus petite erreur $ \epsilon $ est considéré comme le déplacement inverse à effectuer pour aligner les deux images. L'opération revient à effectuer l'optimisation suivante : $$\label{eq:3}\tag{3} \min \left(\epsilon_{\hat{Y},Y}\right) \quad .$$

La référence de l'alignement est le canal contenant l'information sur la couleur rouge (R). L'alignement de l'image couleur (RGB) est obtenu en effectuant les alignements des paires de canaux rouge et vert (G), puis rouge et bleu (B).

Cette technique est appelée «force brute» dans les résultats d'alignement.

Question 2

Un des problèmes de la technique employée à la question 1 est la performance : un très grand nombre de calculs est nécessaire afin d'effectuer la minimisation d'erreur $\eqref{eq:3}$. Une façon d'accélérer cette optimisation est de l'effectuer sur une version simplifiée du problème, soit une image plus petite. On peut réduire la taille des deux images à comparer avant d'effectuer l'algorithme d'alignement décrit à la question 1. Une fois l'alignement effectué sur les images réduites, il suffit de recommencer avec des versions plus grandes des images, jusqu'à retourner à la taille originale après quelques itérations. Cette technique permet d'aligner des images dont le décalage est beaucoup plus grand que ce que permettait l'algorithme de la question 1 pour un même temps de calcul.

Cette technique est appelée «pyramide» dans les résultats d'alignement.

Question 3

Quatre images «maison» furent prises pour le travail: deux de magasines sur des étagères, une de notre dépanneur et une d'un tableau blanc avec mon nom écrit en trois couleurs. Lors du traitement de ces images, deux problèmes firent surface:

Pour corriger ces problèmes, il serait nécessaire de prendre les photos dans un laps de temps très court et d'employer une technique d'alignement à 8 degrés de liberté pouvant réaligner toute transformation de perspective.

Bells and Whistles: Alignment par corrélation de phase

La superposition des images suivi de leur comparaison est une technique appelée la corrélation croisée en traitement de signal. Cette opération est l'opération duale de la convolution. Une façon d'accélérer cette opération est de l'effectuer dans le domaine fréquentiel. En effet, une convolution dans le domaine spatial se traduit en une multiplication de Hadamard (élément par élément) dans le domaine fréquentiel. Puisque la translation en domaine spatial se traduit en déphasage dans le domaine spectral $\left( \mathcal{F}\{f\left(x-a\right)\} = e^{-2\pi i a \xi} \cdot \hat{f}\left(\xi\right) \right)$, il suffit de trouver la différence de phase entre les deux images dans le domaine spectral pour obtenir le décalage nécessaire pour effectuer leur alignement. Pour se faire, on peut multiplier la première image spectrale $f_a$ avec le conjugé complexe de la deuxième image spectrale $f_b^*$ (afin d'obtenir une corrélation et non une convolution). Cette densité spectrale croisée est habituellement normalisée dans la littérature, comme suit: $$\label{eq:4}\tag{4} \Gamma_{a,b} = \frac{ f_a \circ f_b^*}{|f_a \circ f_b^*|} \quad ,$$ où $ \circ $ dénote la multiplication de Hadamard et $ \Gamma_{a,b} $ la densité spectrale croisée. Cette dernière représente la transformée de Fourier de la corrélation croisée. Une fois la transformée de Fourier inverse effectuée sur $ \Gamma_{a,b} $, l'index du pixel le plus lumineux donne la translation à effectuer pour aligner les deux images. Cette technique est expliquée dans [1].

L'algorithme est donc le suivant:

  1. Appliquer une fenêtre aux images afin de limiter le phénomène de Gibbs;
  2. Calculer les transformées de Fourier 2D des deux images;
  3. Calculer $ \eqref{eq:4} $;
  4. Calculer $ \gamma_{a,b} = \mathcal{F}^{-1} \{ \Gamma_{a,b} \} $;
  5. $ \left( \Delta x, \Delta y \right) = \arg \max \left(\gamma_{a,b}\right) $

Cette technique d'alignement comporte plusieurs avantages. Par exemple, elle permet le calcul de l'alignement à une précision sous-pixel (il suffit d'ajouter des zéros à la fin de $ \Gamma_{a,b} $) et elle possède une complexité algorithmique en $ n\log n $, ce qui est mieux que la complexité quadratique de la question 1. On note que cette technique calcule la corrélation sur toute l'image, tandis que l'algorithme de la question 1 ne calcule qu'une infime partie de celle-ci.

[1] Reddy, B. Srinivasa, and Biswanath N. Chatterji. "An FFT-based technique for translation, rotation, and scale-invariant image registration." IEEE transactions on image processing 5.8 (1996): 1266-1271.

Résultats d'alignement

On remarque que les images sont toutes alignées correctement. Les valeurs de translation sont différentes entre l'approche par pyramide et la corrélation de phase parce que le découpage initial a été effectué différemment. Uniquement une partie des résultats sont présentés ici, tous les résultats sont disponibles dans le répertoire /resultats/.

Les temps d'exécution

Exemples d'images fournies

Aucun Alignement Pyramide / Force brute Corrélation de phase
Pyramide: R->G: (2, -19), R->B: (21, -39) R->G: (1, -19), R->B: (121, -33)
Pyramide: R->G: (40, -13), R->B: (62, -38) R->G: (40, -13) R->B: (165, -37)
Pyramide: R->G: (-13, -8), R->B: (-14, -34)
Pyramide: R->G: (5, -11), R->B: (49, -20)

Tableau de tous les déplacements:

Image ΔX ΔY
00757v.jpg (2, -2) (5, -5)
00888v.jpg (1, 0) (2, -1)
00889v.jpg (1, -1) (4, -3)
00907v.jpg (1, 1) (2, 1)
00911v.jpg (-6, 0) (-1, 1)
01031v.jpg (1, 0) (4, -2)
01657v.jpg (1, -1) (2, -1)
01880v.jpg (1, -2) (4, -4)
00029u.tif ( 2, -19) ( 21, -39)
00087u.tif ( 6, -17) ( 21, -68)
00128u.tif ( 40, -13) ( 62, -38)
00458u.tif ( 18, -26) ( 37, -32)
00737u.tif ( 1, -9) ( 22, -13)
00822u.tif (-13, -8) (-14, -34)
00892u.tif (11, -1) (32, -6)
01043u.tif ( 5, -11) ( 49, -20)
01047u.tif (-22, -13) (-19, -34)

Images maison

Voici quatre images alignées en utilisant la corrélation de phase:

Bells and Whistles: Retrait du contour

La première étape du retrait de contour est de s'assurer que chaque image représentant un canal est entière. Pour se faire, il est nécessaire de prendre des zones de l'image d'origine se chevauchant et d'avoir un algorithme d'alignement robuste aux images de différentes tailles.

Pour retirer le contour, la méthode employée calcule les arêtes de l'image en utilisant le gradient tel que précédemment mentionné, puis la différence des trois images est effectuée. Ensuite, une opération de fermeture est appliquée sur tous les pixels qui sont au-dessus d'un seuil manuellement déterminé à $ 0.2 $. Le résultat est ensuite fusionné en un seul masque où chacun des pixels booléens représente lorsque uniquement deux images sources sont d'accord (pour confirmer le fait qu'une seule image impacte toujours sur deux différences). Puis, une somme cumulative est appliquée dans les deux axes afin de trouver la densité d'arêtes dans chacun des axes. Après seuillage, les coordonnées (x, y) des 4 coins intérieurs de l'image, correspondant à l'image sans contour noir.

Résultats finaux

En utilisant la corrélation de phase suivi du retrait de la bordure, voici les résultats (plus disponible dans /resultats/):