Résumé

Les caméras modernes sont incapables de capturer la plage dynamique du monde réel. Dans certaines scènes, même la meilleure photo sera partiellement sous ou sur-exposée. Les chercheurs et les photographes contournent généralement cette limitation en combinant l'information provenant de plusieurs images de la même scène à différentes expositions. Pour ce travail pratique, vous devez écrire un logiciel qui combinera automatiquement plusieurs expositions en une seule image de radiance à haute de plage dynamique, puis convertir cette image de radiance en une image appropriée à l'affichage grâce au processus de reproduction tonale ("tone mapping").

Overview

Modern cameras are unable to capture the full dynamic range of commonly encountered real-world scenes. In some scenes, even the best possible photograph will be partially under or over-exposed. Researchers and photographers commonly get around this limitation by combining information from multiple exposures of the same scene. For this homework, you will write software to automatically combine multiple exposures into a single high dynamic range radiance map, and then convert this radiance map to an image suitable for display through tone mapping.

Détails

Il y a deux composantes majeures à ce projet:

1. Récupérer une image de radiance d'une série d'images;
2. Convertir cette image radiance à une image qui peut être affichée.

Pour la première partie, vous devez récupérer l'inverse de la fonction de l'exposition en fonction de la valeur des pixels, g. Pour la deuxième partie, vous allez commencer par implémenter un opérateur global de reproduction tonale, pour ensuite l'améliorer avec une procédure locale. Chacune des étapes de l'algorithme sont décrites plus bas; deux articles pertinents pour réaliser ces étapes sont Debevec et Malik 1997 et Durand et Dorsey 2002.

Details

There are two major components to this project:

1. Recovering a radiance map from a collection of images;
2. Converting this radiance map into a display image.

For the first part, you will have to recover the inverse of the function mapping exposure to pixel value, g. For the second part, you will start by implementing a global tone mapping operator, and you then have to improve it to a local tone mapping procedure. All of these steps are described further down; two relevant papers for these steps are Debevec and Malik 1997 and Durand 2002.

Pour cette étape, notre but est de construire une image de radiance à partir de plusieurs exposition à faible plage dynamique. Il est grandement recommandé de lire les sections 2.1 et 2.2 de Debevec and Malik 1997 pour bien comprendre ce processus. Vous trouverez ci-dessous un résumé de ces sections.

La valeur de pixel observée Z_ij pour le pixel i d'une image j est une fonction de la radiance d'une scène (inconnue) et d'une durée d'exposition (connue):

E_i est la radiance de la scène au pixel i, et l'intégrale dans le temps de la radiance de la scène E_i Δ t_j est l'exposition a un pixel donné. En général, f, la courbe de réponse des pixels, peut être assez compliquée. C'est pour cette raison que nous ne résolvons pas pour f, mais pour g=ln(f^-1) qui place les valeurs des pixels (de 0 à 255) sur le log des valeurs d'exposition:

(equation 2 dans l'article de Debevec). Résoudre g peut sembler impossible (et, en effet, nous pouvons uniquement retrouver g jusqu'à un facteur d'échelle) puisque nous ne connaissons ni g, ni E_i. La solution provient du fait que la scène est statique, et bien que nous ne pourrions pas savoir la valeur absolue de E_i à chaque pixel i, nous savons que la valeur relative reste constante sur toutes les images d'une séquence.

Résoudre pour g est la partie compliquée du problème. En effet, une fois que l'on a obtenu g, nous pouvons associer les valeurs des pixels observés et le temps d'exposition à la radiance grâce à l'équation suivante:

Ceci est l'équation 5 dans Debevec.

Pour rendre les résultats robustes, nous voulons considérer deux détails supplémentaires:

• Nous allons faire l'hypothèse que g est lisse. Debevec ajoute une contrainte à notre système linéaire pénalisant g par rapport à l'ampleur de sa dérivée seconde. Puisque g est une fonction discrète (definie uniquement à des valeurs entières de g(0) à g(255)), nous pouvons approximer la dérivée seconde avec des soustractions sur des valeurs entières, e.g. g''(x) = (g(x-1) - g(x)) - (g(x) - g(x+1)) = g(x-1) + g(x+1) - 2*g(x). Nous aurons une équation de ce genre pour chaque nombre entier dans le domaine de g, mis-à-part pour g(0) et g(255) où la deuxième dérivée est indéfinie.
• Chaque temps d'exposition nous donne uniquement de l'information fiable sur certains pixels (i.e. les pixels correctement exposés pour cette image). Pour les pixels sombres, la contribution relative du bruit est haute et, pour les pixels très clairs, le capteur aurait pu être saturé. Pour que nos estimations de E_i soient plus précises, nous devons mesurer la contribution de chaque pixel selon l'équation 6 dans Debevec. Un exemple d'une fonction mesurant la contribution de chaque pixel w est une fonction triangulaire avec une amplitude maximale à Z=127.5 et dont la valeur est zéro à Z=0 et Z=255. Cette pondération devrait être utilisée pour résoudre g et lorsque vous utilisez g pour créer la carte de radiance HDR pour tous les pixels.

We want to build an HDR radiance map from several LDR exposures. It is highly recommended that you read Sections 2.1 and 2.2 in Debevec and Malik 1997 to help understand this process. Below is a summary.

The observed pixel value Z_ij for pixel i in image j is a function of unknown scene radiance and known exposure duration:

E_i is the unknown scene radiance at pixel i, and scene radiance integrated over some time E_i Δ t_j is the exposure at a given pixel. In general, f might be a somewhat complicated pixel response curve. We will not solve for f, but for g=ln(f^-1) which maps from pixel values (from 0 to 255) to the log of exposure values:

(equation 2 in Debevec). Solving for g might seem impossible (and indeed, we only recover g up to a scale factor) because we know neither g nor E_i. The key observation is that the scene is static, and while we might not know the absolute value of E_i at each pixel i, we do know that the value remains constant across the image sequence.

Solving for g is the tricky part. After we have g, is it straightforward to map from the observed pixels values and shutter times to radiance by the following equation (rearranged from above):

This is Equation 5 in Debevec.

To make the results robust, we consider two additional things:

• We expect g to be smooth. Debevec adds a constraint to our linear system which penalizes g according to the magnitude of its second derivative. Since g is discrete (defined only at integer values from g(0) to g(255) we can approximate the second derivative with finite differences, e.g. g''(x) = (g(x-1) - g(x)) - (g(x) - g(x+1)) = g(x-1) + g(x+1) - 2*g(x). We will have one such equation for each integer in the domain of g, except for g(0) and g(255) where the second derivative would be undefined.
• Each exposure only gives us trustworthy information about certain pixels (i.e. the well exposed pixels for that image). For dark pixels the relative contribution of noise is high and for bright pixels the sensor may have been saturated. To make our estimates of E_i more accurate we need to weight the contribution of each pixel according to Equation 6 in Debevec. An example of a weighting function w is a triangle function that peaks at Z=127.5, and is zero at Z=0 and Z=255. This weighting should be used both when solving for g and when using g to build the HDR radiance map for all pixels.

Trouver l'image de radiance est la moitié du problème. Vous voulez maintenant afficher votre image sur un écran à 8 bits. Il existe quelques opérateurs de reproduction tonale avec lesquels vous pouvez expérimenter, comme log(L), sqrt(L), et L/(1+L). Peu importe quelle transformation vous utilisez (n'oubliez pas de l'indiquer dans le rapport), vous devez étendre les valeurs d'intensité dans l'image résultante pour remplir l'écart de [0 255] pour maximiser votre contraste.

Une fois que votre méthode de tone mapping globale fonctionne bien, vous devez aussi implémenter la méthode de reproduction tonale locale suivante. Il s'agit d'une version simplifiée de la méthode proposée par Durand et Dorsey en 2002. Les étapes vont comme suit:

1. Vos entrées sont des valeurs RGB linéaires de radiance.
2. Calculez l'intensité (I).
3. Calculez la chrominance: (R/I, G/I, B/I)
4. Calculez le log de l'intensité: L = log2(I)
5. Filtrez cette valeur avec un filtre bilatéral (vous devez l'implémenter vous-même): B = bf(L)
6. Calculez les "détails" (detail layer): D = L - B
7. Normalisez la base avec une mise à l'échelle et une constante: B' = (B - o) * s
   1. La constante est telle que le maximum d'intensité de la base est 1. Puisque les valeurs sont dans le domaine logarithmique, o = max(B).
   2. Le facteur d'échelle est choisi de sorte que la base ait une plage dynamique (dR) intéressante, i.e., s = dR / (max(B) - min(B)). Essayez des valeurs entre 2 et 8 pour dR, cela devrait couvrir une plage intéressante. Une valeur autour de 4 ou 5 devrait être correcte.
8. Reconstruisez l'intensité logarithmique: O = 2^(B' + D)
9. Remettez les couleurs: R',G',B' = O * (R/I, G/I, B/I)
10. Appliquez une compression gamma. Sans la compression gamma le résultat sera trop sombre. Une valeur autour de 0.5 devrait être correcte (e.g. resultat.^0.5).

Votre image HDR va contenir des 0 ce qui va apporter des problèmes avec log. Vous pouvez régler ce problème en remplaçant tous les zéros par la plus petite valeur non-nulle dans l'image, multipliée par un petit facteur d'échelle (e.g. 0.001).

Getting the radiance map is only half the battle. You want to be able to show off your image clearly. There are a few global tone-mapping operators to play with, such as log(L), sqrt(L), and L / (1+L). Regardless of which transform you use, you'll want to stretch the intensity values in the resulting image to fill the [0 255] range for maximum contrast.

Once you see something reasonable using a global tone mapping operation, you'll need to implement a local method for full credit.

You'll be implementing a simplified version of Durand and Dorsey 2002. The steps are roughly as follows:

1. Your input is linear RGB values of radiance.
2. Compute the intensity (I).
3. Compute the chrominance channels: (R/I, G/I, B/I)
4. Compute the log intensity: L = log2(I)
5. Filter that with a bilateral filter (you have to implement it yourself): B = bf(L)
6. Compute the detail layer: D = L - B
7. Apply an offset and a scale to the base: B' = (B - o) * s
   1. The offset is such that the maximum intensity of the base is 1. Since the values are in the log domain, o = max(B).
   2. The scale is set so that the output base has dR stops of dynamic range, i.e., s = dR / (max(B) - min(B)). Try values between 2 and 8 for dR, that should cover an interesting range. Values around 4 or 5 should look fine.
8. Reconstruct the log intensity: O = 2^(B' + D)
9. Put back the colors: R',G',B' = O * (R/I, G/I, B/I)
10. Apply gamma compression. Without gamma compression the result will look too dark. Values around 0.5 should look fine (e.g. result.^0.5).

Your HDR images have zero values which will cause problems with log. You can fix this problem by replacing all zeroes by some factor times the smallest non-zero values.

Conseils d'implémentation

• Pour créer des résultats clairement meilleurs que des images à une seule exposition, assurez-vous que votre scène contient réellement une plage dynamique élevée! Les images que l'on vous a fournies contiennent des combinaisons d'éléments très lumineux et très sombres. Donc, si vous prenez vos propres photos, assurez-vous qu'elles contiennent combinaisons d'éléments intérieurs et extérieurs, ou des scènes qui sont fortement rétro-éclairées, ou bien des images où les sources lumineuses (ampoules) sont directement visibles.
• Vous pouvez obtenir le temps d'exposition pour chaque image à partir de l'entête EXIF dans les fichiers JPG directement dans Matlab. Pour ce faire, utilisez la fonction imfinfo. Lorsque cette information n'est pas disponible, nous vous l'avons fournie dans un fichier texte.
• Si vous prenez les images sans trépied, alignez géométriquement vos images au préalables pour corriger le mouvement de la caméra avec une homographie. Cela risque d'améliorer vos résultats.
• L'implémentation naïve d'un filtre bilatéral est un peu lente (quelques minutes d'exécution). Vous pouvez tenter d'optimiser votre implémentation (voir crédits supplémentaires).
• Le code Matlab pour retrouver g est disponible dans Debevec et Malik 1997. Bien entendu, vous pouvez le regarder et même le tester comme référence. Par contre, le code que vous retournez devrait être le vôtre.
• Le filtre bilatéral est sensible aux valeurs de sigma pour les gaussiennes (pour l'intensité et pour les coordonnées).
• La composante spatiale du filtre devrait, idéalement, avoir une grande étendue spatiale, mais cela rendrait une implémentation naïve très lente.

Implementation Tips

• To create results which are clearly better than any single exposure, make sure your scene actually has a high dynamic range! The data we give you have combinations of very bright (sunlight or bright lights) and very dark (shadowed) elements. If you take your own photos, take photos that have combinations of indoor and outdoor elements, or scenes that are heavily back-lit, possibly with the light source directly visible.
• You can obtain the exposure time for each photo from the EXIF header in the JPG files from within Matlab. Look at the imfinfo function. When the information is not available, we wrote it in a text file.
• For your own photos, and for any photos with long exposures, geometric alignment to correct for camera motion will probably improve your results.
• The brute force implementation of a bilateral filter is a little slow (couple of minutes for the test images). You can try to optimize it (see bells and whistle).
• Matlab code for recovering g is available in Debevec and Malik 1997. You can of course look at it and even test it as a reference. However, the code you turn in must be your own.
• The bilateral filter is sensitive to the sigma values for the intensity Gaussian and spatial Gaussian.
• The spatial component of the bilateral filter should, ideally, have a large spatial extent but that will make a naive implementation quite slow.

Crédits supplémentaires

Essayez une de ces idées pour approfondir vos connaissances (et augmenter votre note):

• (jusqu'à 5 points) Utilisez vos propres images! Dans la plupart des caméras, vous pouvez facilement obtenir plusieurs expositions grâce à la fonction "Exposure Bracketing". Si vous avez besoin d'aide, n'hésitez pas à demander.
• (jusqu'à 10 points) Implémentez une version optimisée du filtrage bilatéral dans lequel les convolutions gaussiennes sont effectuées dans un espace 3d sous-échantillonnée. Discutez du gain en performance atteint par rapport à la méthode standard.
• (jusqu'à 20 points) Implémentez une autre méthode de reproduction tonale. Par exemple, une version simplifiée de cet algorithme.
• Certaines des images de test ne sont pas parfaitement alignées. Essayez de corriger automatiquement l'alignement géométrique des images avant d'appliquer votre algorithme de HDR. Si vous voulez tester sur des images non-alignées, vous pouvez essayer celles-ci, ou bien prendre vos propres photos.
  • (jusqu'à 2.5 points): Pour une implémentation en pyramides similaire à celle du TP1.
  • (jusqu'à 5 points): Pour une méthode d'alignement local similaire à celle du TP4.

Bells and Whistles

Try one of these ideas to increase your understanding on this subject (and increase your score):

• (up to 5 points) Use your own images! In most cameras, this can be done easily with the "Exposure Bracketing" feature. If you need help, ask us.
• (up to 10 points) Implement an optimized bilateral filter, in which Gaussian convolutions are performed in a subsampled 3d space. Discuss the speedup achieved over the baseline method.
• (up to 20 points) Implement any other local tone mapping algorithm. For example a simplified version of this algorithm.
• Some of the test image stacks are not well aligned. Try to automatically correct the geometric alignment of the images before running your HDR pipeline. If you want to test on images that are not well aligned, you can try on those ones, or take your own photos.
  • (up to 2.5 points): For a pyramid implementations similar to Homework 1.
  • (up to 5 points): For local alignment methods (local features).

Livrables

Comme lors des travaux précédents, celui-ci sera remis dans un format page Web. Rappel: le site Web n'a pas besoin d'être esthétiquement agréable; ne faite que décrire ce que vous avez fait.

Plus précisément, la page devrait contenir:

• Explication de l'algorithme que vous avez implémentez. Illustrez toutes les étapes de l'algorithme. Entre autres, pour chaque image HDR, montrez l'image de radiance et la fonction g estimée.
• Explication des résultats produits sur toutes les séquences d'images dans le fichier suivant: images.zip. Discutez des bons et des mauvais résultats.
• Expliquez les crédits supplémentaires que vous avez implémentés, et illustrez ces explications. Aussi, si cela s'applique, montrez correctement le "avant" et le "après".

Pour la remise de votre travail, créez un fichier tp5.zip qui contient

• Votre rapport en format HTML dans un dossier tp5/web. Vos images doivent être dans un dossier tp5/web/images.
• Votre code matlab doit être dans un dossier tp5/code. De préférence, ne pas inclure les images que vous avez utilisées pour produire vos résultats dans ce dossier dans le but de ne pas alourdir le fichier.

Finalement, veuillez téléverser votre fichier tp5.zip sur pixel (http://pixel.fsg.ulaval.ca) avant la date limite. Bien entendu, la politique des retards mentionnée dans le plan de cours sera appliquée. Pour toutes questions concernant la procédure de remise ou le travail en tant que tel, envoyez vos questions à l'adresse courriel du cours.

Deliverables

As in the previous homework, this one will be handed in a webpage format. Remember: the aesthetics of the website will not be evaluated, but it is important that the information be presented clearly.

More precisely, the webpage should:

• Explain the algorithm you have implemented. Illustrate every step of the algorithm. Among other things, for each HDR image, show the radiance map and the recovered response function g.
• Explanation of the generated results obtained on all the images in the following file: images.zip. Discuss good and bad results.
• Explain all the bells and whistles you have implemented. Try to illustrate these explanations. Also, if applicable, show correctly the "before" and "after".

For this homework, you must create a tp5.zip file. In this file you'll put:

• Your report in the HTML format inside a folder named tp5/web. Your images for this web page should be inside a folder named tp5/web/images.
• Your matlab code should be put inside the folder tp5/code. Please do not include the images you have used to produce your results inside this folder.

Finally, you should upload the file tp5.zip on Pixel (http://pixel.fsg.ulaval.ca) before the deadline. Naturally, the late submission policy in the course plan will be applied. For any question regarding the submission process or the homework as such, send your question to the course's email address.

Évaluation

Ce travail est évalué sur 100 points. La répartition des points va comme suit:

• (40 pts, 30 pour les étudiants gradués) pour l'implémentation de l'algorithme pour reconstruire l'image de radiance à haute plage dynamique et de la fonction de réponse inverse à partir de plusieurs expositions (explications, code, images). Appliquez votre implémentation sur toutes les séquences d'images fournies.
• (20 pts, 20 pour les étudiants gradués) pour l'implémentation de la reproduction tonale globale (explications, code, images). Appliquez votre implémentation sur toutes les séquences d'images fournies.
• (40 pts, 30 pour les étudiants gradués) pour l'implémentation de la reproduction tonale locale (explications, code, images). Appliquez votre implémentation sur toutes les séquences d'images fournies.
• (N pts) pour les crédits supplémentaire. Rappel: les étudiants gradués doivent livrer au moins 20 pts de crédits supplémentaires.

Evaluation

This assignment is evaluated on 100 points, as follows:

• (40 pts, 30 pts for graduate students) for the implementation of the HDR radiance map and inverse response function estimation from multiple exposures (code, results, explanations). Use this implementation on all the given images sequences.
• (20 pts, 20 pts for graduate students) for the implementation of the global tone mapping (code, results, explanations). Use this implementation on all the given images sequences.
• (40 pts, 30 pts for graduate students) for the implementation of the local tone mapping (code, results, explanations). Use this implementation on all the given images sequences.
• (N pts) for any bells and/or whistles you've added in. Reminder: graduate students must add to their assignment at least 20 pts worth of bells and whistles.

Liens rapides

• Ébauche de page web;
• "Recovering High Dynamic Range Radiance Maps from Photographs" par Debevec et Malik;
• "Fast bilateral filtering for the display of high-dynamic-range images" par Durand et coll.;
• Séquences d'images sur lesquelles vous devez tester vos algorithmes.
• Séquences d'images non-alignées pour les crédits supplémentaires.

Quick links

• Webpage template;
• "Recovering High Dynamic Range Radiance Maps from Photographs" by Debevec and Malik;
• "Fast bilateral filtering for the display of high-dynamic-range images" by Durand and al.;
• Images sequences on which to test your algorithms.
• Non-aligned image sequences for bells & whistles

Remerciements

Merci à James Hays d'avoir créé le TP original qui a servi d'inspiration pour celui-ci!

Thanks

Many thanks to James Hays for creating the assignment which inspired this one!