Une image HDR (HDR) contient une variation d'intensité
et de couleur plus large que celle autorisée par le format RGB classique avec 1 octet par canal que nous
avons utilisé dans la tâche précédente.
Pour stocker ces informations supplémentaires, nous utilisons un flottant à simple précision
pour chaque canal au lieu d'un entier. Cela permet une plage de valeurs d'intensité extrêmement large.
L'image de cette tâche représente l'intérieur d'une église avec de la lumière entrant
à travers des vitraux. Les valeurs brutes données en entrée sont des floattants. Pour chaque
canal, elles vont de 0 à 275. Cependant, la moyenne est de 0,41 et 98 % des valeurs sont
inférieures à 3 ! Cela signifie que certaines zones (les fenêtres) sont
extrêmement lumineuses par rapport à partout ailleurs. Si nous mappons linéairement cette plage [0-275] dans la
plage [0-255] que nous avons utilisée, la plupart des valeurs seront mappées à zéro !
La seule chose que nous pourrons voir sont les zones les plus lumineuses - les
fenêtres - tout le reste apparaîtra noir.
Le problème est que, même si nous disposons de caméras capables d'enregistrer la large
gamme d'intensité qui existe dans le monde réel, nos écrans ne sont pas
capables de les afficher. Nos yeux sont également tout à fait capables d'observer une
gamme d'intensités bien plus large que celle que nos formats d'image/écran sont
capables d'afficher.
Le Tone Mapping est un processus qui transforme les intensités de l'image de sorte que
les valeurs les plus lumineuses ne soient pas si éloignées de la moyenne. De cette façon, lorsque
nous transformons les valeurs en [0-255], nous pouvons réellement voir l'image entière.
Il existe de nombreuses façons d'effectuer ce processus et c'est autant un art qu'une
science - il n'y a pas de "bonne" réponse unique. Dans ce devoir, nous
mettrons en œuvre une technique possible.
Chrominance-Luminance d'arrière-plan
L'espace RGB que nous avons utilisé pour représenter les images peut être considéré
comme un ensemble possible d'axes couvrant un espace tridimensionnel de couleur. Nous choisissons parfois d'autres axes pour représenter cet
espace car ils rendent certaines opérations plus pratiques.
Une autre manière possible de représenter une image couleur consiste à séparer les informations de couleur (chromaticité) des informations de luminosité. Il existe plusieurs méthodes différentes pour le faire - une méthode courante à l'époque de la télévision analogique était connue sous le nom de Chrominance-Luminance ou YUV.
Nous choisissons de représenter l'image de cette manière afin de pouvoir remapper uniquement le canal d'intensité, puis de recombiner les nouvelles valeurs d'intensité avec les informations de couleur pour former l'image finale.
Les anciens signaux de télévision étaient transmis de cette manière afin que les téléviseurs noir et blanc puissent afficher le canal de luminance tandis que les téléviseurs couleur affichaient les trois canaux.
Mappage de tons
Dans cette tâche, nous allons transformer le canal de luminance (en fait, le logarithme de la luminance, mais cela n'est pas important pour les parties de l'algorithme que vous allez implémenter) en compressant sa plage à [0, 1].
Pour ce faire, nous avons besoin de la distribution cumulative des valeurs de luminance.
Exemple
- input : [2 4 3 3 1 7 4 5 7 0 9 4 3 2]
- min / max / range: 0 / 9 / 9
- histo avec 3 bins : [4 7 3]
- cdf : [4 11 14]
Votre tâche consiste à calculer cette distribution cumulative en suivant ces
étapes.
Faites un git pull pour être sur d'être à jour.
Accédez au dossier du TP dans cs344/problem_sets/problem_set3. Le fichier à modifier se nomme student_func.cu. Dans ce TP, nous allons implémenter la fonction your_histogram_and_prefixsum qui appelera les différents kernels et gèrera la mémoire temporelle allouée sur le GPU.
Nous allons commencer par écrire l'étape permettant de calculer le minimum et le maximum dans la liste des flottants. De quel genre d'opérateur s'agit-il ?
En vous inspirant du code vu en cours (disponible dans
cs344/lesson_code_snippets/lesson3_code_snippets), écrivez un kernel permettant de calculer le minimum et le maximum donné dans
d_logLuminance en utilisant de la mémoire partagée. Afficher le minimum et le maximum à l'écran pour vérifier votre réponse (les valeurs attendues s'afficheront lors de l'exécution du programme). Tester votre code avec
./HW3 memorial_raw.png.
Maintenant, il va vous falloir gérer la mémoire et les appels aux kernels tout seul !
Au lieu de faire deux kernels pour le minimum et le maximum, rajoutez un argument à votre kernel pour traiter les deux cas dans une seule fonction.
Dans ce TP, contrairement à l'exemple du cours, nous ne supposons pas que la taille de l'image est un multiple de 1024. Par contre, vous supposons toujours que l'image est plus petite que 1024x1024.
Ajoutez cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError()); après les appels de kernel pour vérifiez les éventuelles erreurs.
__global__ void minmax_reduce(float * d_out, const float * d_in, int image_size, bool compute_min) {
extern __shared__ float sdata[];
int pos = blockDim.x * blockIdx.x + threadIdx.x;
int t_idx = threadIdx.x;
if (pos >= image_size) return;
sdata[t_idx] = d_in[pos];
__syncthreads();
for (int i = blockDim.x / 2; i > 0; i >>= 1){
int matching_idx = t_idx + i;
int matching_pos = pos + i;
if (t_idx < i && matching_pos < image_size) {
if (compute_min){
sdata[t_idx] = min(sdata[t_idx], sdata[matching_idx]);
} else {
sdata[t_idx] = max(sdata[t_idx], sdata[matching_idx]);
}
}
__syncthreads();
}
if (t_idx == 0) {
d_out[blockIdx.x] = sdata[0];
}
}
void reduce(float * d_max, float * d_min, float * d_intermediate, const float * d_in, int image_size){
const int max_threads = 1024;
const int n_blocks = image_size / max_threads;
minmax_reduce<<<n_blocks, max_threads, max_threads * sizeof(float)>>>(d_intermediate, d_in, image_size, true);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
minmax_reduce<<<1, 1024, n_blocks * sizeof(float)>>>(d_min, d_intermediate, n_blocks, true);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
minmax_reduce<<<n_blocks, max_threads, max_threads * sizeof(float)>>>(d_intermediate, d_in, image_size, false);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
minmax_reduce<<<1, 1024, n_blocks * sizeof(float)>>>(d_max, d_intermediate, n_blocks, false);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
}
void your_histogram_and_prefixsum(const float* const d_logLuminance,
unsigned int* const d_cdf,
float &min_logLum,
float &max_logLum,
const size_t numRows,
const size_t numCols,
const size_t numBins)
{
int image_size = numRows * numCols;
// Allocate memory
float * d_min;
float * d_max;
float * d_intermediate;
cudaMalloc((void **) &d_min, sizeof(float));
cudaMalloc((void **) &d_max, sizeof(float));
cudaMalloc((void **) &d_intermediate, image_size / 1024 * sizeof(float));
reduce(d_max, d_min, d_intermediate, d_logLuminance, image_size);
cudaMemcpy(&min_logLum, d_min, sizeof(float), cudaMemcpyDeviceToHost);
cudaMemcpy(&max_logLum, d_max, sizeof(float), cudaMemcpyDeviceToHost);
cudaFree(d_min);
cudaFree(d_max);
cudaFree(d_intermediate);
}
Maintenant, calculez l'histogramme des valeurs dans le tableau. Le bin d'une valeur sera donnée par la formule
bin = (lum[i] - lumMin) / lumRange * numBins, où numBins est passé en paramètre.
Pour vérifier vos calculs, affichez la valeur dans bin 0, 1014 et 1015.
Pour l'instant, vous pouvez utiliser une addition atomique.
__global__ void compute_histogram(int * d_bins, const float * d_in, int image_size, int min_logLum, float range, int numBins){
int pos = blockIdx.x * blockDim.x + threadIdx.x;
if (pos >= image_size) return;
int n_bin = ((d_in[pos] - min_logLum) / range) * numBins;
atomicAdd(&(d_bins[n_bin]), 1);
}
void your_histogram_and_prefixsum(const float* const d_logLuminance,
unsigned int* const d_cdf,
float &min_logLum,
float &max_logLum,
const size_t numRows,
const size_t numCols,
const size_t numBins)
//...
int * d_bins;
cudaMalloc((void **) &d_bins, numBins * sizeof(int));
int h_bins[numBins];
for (int i = 0; i < numBins; i++){
h_bins[i] = 0;
}
cudaMemcpy(d_bins, h_bins, numBins * sizeof(int), cudaMemcpyHostToDevice);
compute_histogram<<<image_size / 1024, 1024>>>(d_bins, d_logLuminance, image_size, min_logLum, range, numBins);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
cudaMemcpy(h_bins, d_bins, numBins * sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(d_bins);
}
Nous voulons maintenant calculer un scan exclusif de l'histogramme pour obtenir la distribution cumulative et assigner le resultat à d_cdf (donné en entrée). Quel algorithme permet de directement calculer un scan exclusif ?
L'algorithme de Blelloch.
Implémentez cet algorithme. Vérifiez vos résultats en affichant la valeur de la distribution en 0, 1 et numBins - 1. Vérifiez la sortie du programme.
On suppose dans cette question que l'histogramme fait moins de 1024 éléments (donc, nous n'avons qu'un seul block).
__global__ void cumulative_sum(unsigned int* const d_cdf, int * d_bins, int numBins){
// Sum-sum reduce
int pos = threadIdx.x;
d_cdf[pos] = d_bins[pos];
__syncthreads();
for (int i = 2; i < numBins; i <<= 1){
if (pos % i == i - 1) {
d_cdf[pos] += d_cdf[pos - (i >> 1)];
}
__syncthreads();
}
// Downsweep
if (pos == numBins - 1){
d_cdf[pos] = 0;
}
__syncthreads();
for (int i = numBins; i > 0; i >>= 1){
int current = d_cdf[pos];
__syncthreads();
int next_shift = i >> 1;
if (pos % i == i - 1) {
d_cdf[pos - next_shift] = current;
} else if ((pos + next_shift) % i == i - 1) {
d_cdf[pos + next_shift] += current;
}
__syncthreads();
}
}
void your_histogram_and_prefixsum(const float* const d_logLuminance,
unsigned int* const d_cdf,
float &min_logLum,
float &max_logLum,
const size_t numRows,
const size_t numCols,
const size_t numBins){
//...
cumulative_sum<<<1, numBins>>>(d_cdf, d_bins, numBins);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
}
