CI5 : Retirer les yeux rouges
Retirer les yeux rouges d'une photo avec un algorithme de tri.
Retirer les yeux rouges (∼90mn, – moyen)
Dans cet exercice, nous allons implémenter un algorithme permettant de retirer les yeux rouges d'une photo. Pour ce faire, nous allons utiliser un score nous indiquant à quel point un pixel a des chances d'être un œil rouge. Ce scorage est déjà implémenté pour vous. Vous allez recevoir les scores et devoir les trier par ordre croissant pour savoir quels sont les pixels à modifier pour retirer les yeux rouges.
Chaque score est associé à une position. Quand vous triez les scores, il faudra faire suivre les positions.
Nous allons utiliser le tri Radix pour trier les scores. L'idée va être de construire un histogram pour chaque bit (0 ou 1). Ensuite, nous allons implémenter un compact afin de déterminer la position relative de chaque 0. Enfin, en utilisant l'histogramme et les position relative, nous allons effectuer les déplacements de l'algorithme Radix.
Commencez par faire un git pull pour être sur d'être sur la bonne version du code.
Rendez-vous ensuite dans problem_sets/problem_set4 où vous trouverez le fichier
student_func.cu. Exécutez cmake . pour générer le Makefile puis make pour
compiler le code.
Pour effectuer le tri Radix, nous allons devoir boucler sur tous les bits et exécuter les
kernels que nous allons implémenter à chaque tour de boucle. Écrivez cette boucle dans la fonction
your_sort.
for (int i = 0; i < sizeof(unsigned int) * 8; i++){
Écrivez maintenant un kernel calculant l'histogramme des 0 et de 1 situés à l'index i
(l'index de la boucle) de chaque score. N'oubliez pas de gérer la mémoire.
__global__ void histogram_kernel(unsigned int bit_index,
unsigned int * d_bins,
unsigned int* const d_input,
const int size) {
int pos = threadIdx.x + blockDim.x * blockIdx.x;
if(pos >= size) return;
int bin = ((d_input[pos] & (1<<bit_index)) >> bit_index);
atomicAdd(&d_bins[bin], 1);
}
Nous allons maintenant écrire un compact retournant la position relative de chaque élément ayant un
0 à l'index i (l'index de la boucle). Vous pouvez réutiliser et adapter le code du dernier TP.
Par contre, nous ne supposons plus que tout tient dans un seul block, ni que l'entrée a la taille une puissance
de deux. Nous rappelons que l'algorithme de Blelloch ne fonctionne que si nous travaillons sur des puissances de 2.
Une manière de contourner le problème consiste à définir la taille de la sortie comme étant la puissance de 2
supérieur à la taille de l'entrée.
Il va nous falloir procéder en trois étapes.
Il va nous falloir procéder en trois étapes.
- Effectuez un scan sur chaque block indépendamment. Ce scan retournera aussi un tableau de taille le nombre de blocks et où chaque case contient la somme des élements du block associé. On fera bien attention de définir la taille de ce tableau de sommes comme une puissance de 2.
- On effectue un scan exclusif sur ce tableau de sommes.
- On ajoute à tous les éléments d'un block la somme des totaux précédents.
int find_next_power_of_2(unsigned int x){
unsigned int power = 1;
while (power < x) power <<= 1;
return power;
}
__global__ void cumulative_sum(unsigned int* const d_in,
unsigned int * d_out,
unsigned int * d_block_sums,
int nElems,
int current_bit){
// Sum-sum reduce
int tidx = threadIdx.x;
int block_pos = blockDim.x * blockIdx.x;
int global_pos = tidx + block_pos;
if (global_pos >= nElems) {
d_out[global_pos] = 0; // On est dans la partie artificiellement mise à 0 pour atteindre une puissance de 2.
} else {
if (current_bit == -1){ // Le -1 permet de faire l'algorithme classique
d_out[global_pos] = d_in[global_pos];
} else { // Ici, on calcule le prédicat du compact en prenant le current_bit ième bit.
d_out[global_pos] = 1 - ((d_in[global_pos] & (1 << current_bit)) >> current_bit);
}
}
__syncthreads();
int last = 0; // On on fait un scan exclusif, on doit se rappeler du dernier élément pour le rajouter à la somme de tous les éléments
if (tidx == blockDim.x - 1){ // On travaille bien sur des blocks, et non sur tout le tableau
last = d_out[global_pos];
}
for (int i = 2; i < blockDim.x; i <<= 1){
if (tidx % i == i - 1) {
d_out[global_pos] += d_out[global_pos - (i >> 1)];
}
__syncthreads();
}
// Downsweep
if (tidx == blockDim.x - 1){
d_out[global_pos] = 0;
}
__syncthreads();
for (int i = blockDim.x; i > 0; i >>= 1){
int current = d_out[global_pos];
__syncthreads();
int next_shift = i >> 1;
if (tidx % i == i - 1) {
d_out[global_pos - next_shift] = current;
} else if ((tidx + next_shift) % i == i - 1) {
d_out[global_pos + next_shift] += current;
}
__syncthreads();
}
if (tidx == blockDim.x - 1){ // On met à jour le tableau de la somme sur tout le block
d_block_sums[blockIdx.x] = d_out[global_pos] + last;
}
}
__global__ void add_block_sum(unsigned int* d_scan,
unsigned int* d_block_sum,
int nElems) {
int tidx = threadIdx.x;
int block_pos = blockDim.x * blockIdx.x;
int global_pos = tidx + block_pos;
if (global_pos >= nElems) return;
d_scan[global_pos] += d_block_sum[blockIdx.x];
}
void compact(unsigned int* const d_inputVals,
unsigned int* const d_scan,
const size_t numElems,
int current_bit){
unsigned int threads_per_block = 1024;
int scan_size = find_next_power_of_2(numElems); // On fait bien attention d'avoir des puissances de 2.
unsigned int n_blocks_scan = scan_size / threads_per_block;
int block_sum_size = find_next_power_of_2(n_blocks_scan);
unsigned int* d_block_sum;
checkCudaErrors(cudaMalloc((void**) &d_block_sum, block_sum_size * sizeof(unsigned int)));
checkCudaErrors(cudaMemset((void**) d_block_sum, 0, block_sum_size * sizeof(unsigned int)));
unsigned int* d_block_sum_dummy; // Ne sert à rien, mais on va devoir le donner en entrée plus tard
checkCudaErrors(cudaMalloc((void**) &d_block_sum_dummy, sizeof(unsigned int)));
checkCudaErrors(cudaMemset((void**) d_block_sum_dummy, 0, sizeof(unsigned int)));
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
// On fait le scan par block
cumulative_sum<<<n_blocks_scan, threads_per_block>>>(d_inputVals, d_scan, d_block_sum, numElems, current_bit);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
// On fait un scan sur les totaux de chaque block
cumulative_sum<<<1, block_sum_size>>>(d_block_sum, d_block_sum, d_block_sum_dummy, block_sum_size, -1);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
// On remet les blocks au bon niveau
add_block_sum<<<n_blocks_scan, threads_per_block>>>(d_scan, d_block_sum, numElems);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
checkCudaErrors(cudaFree(d_block_sum));
checkCudaErrors(cudaFree(d_block_sum_dummy));
}
void your_sort(unsigned int* const d_inputVals,
unsigned int* const d_inputPos,
unsigned int* const d_outputVals,
unsigned int* const d_outputPos,
size_t numElems)
{
unsigned int threads_per_block = 1024;
unsigned int n_blocks = numElems / threads_per_block + 1;
printf("numElems: %d\nNum blocks: %d\n", (int) numElems, n_blocks);
unsigned int numBins = 2;
unsigned int *d_hist;
unsigned int h_hist[2];
checkCudaErrors(cudaMalloc((void**) &d_hist, numBins * sizeof(unsigned int)));
int scan_size = find_next_power_of_2(numElems);
unsigned int* d_scanned;
checkCudaErrors(cudaMalloc(&d_scanned, scan_size * sizeof(unsigned int)));
for (int i = 0; i < sizeof(unsigned int) * 8; i++){
checkCudaErrors(cudaMemset((void**) d_hist, 0, numBins * sizeof(unsigned int)));
histogram_kernel<<<n_blocks, threads_per_block>>>(i, d_hist, d_inputVals, numElems);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
compact(d_inputVals, d_scanned, numElems, i);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
checkCudaErrors(cudaMemcpy(h_hist, d_hist, 2 * sizeof(unsigned int), cudaMemcpyDeviceToHost));
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
} ...
Finalement, implémentez un kernel qui déplace les éléments par rapport au compact et à l'histogramme
calculé précédemment.
La nouvelle position d'un 0 est directement donnée par le compact. La position d'un 1 est donnée
par la première valeur de l'histogramme + la position actuelle - le compact à cette position.
__global__ void move(unsigned int* d_inputVals,
unsigned int* d_inputPos,
unsigned int* d_outputVals,
unsigned int* d_outputPos,
unsigned int* d_scanned,
int nElems,
int current_bit,
unsigned int shift){
int pos = threadIdx.x + blockDim.x * blockIdx.x;
if (pos >= nElems) return;
int new_pos = 0;
if ((d_inputVals[pos] & (1 << current_bit)) >> current_bit) {
new_pos = shift + pos - d_scanned[pos];
} else {
new_pos = d_scanned[pos];
}
d_outputVals[new_pos] = d_inputVals[pos];
d_outputPos[new_pos] = d_inputPos[pos];
}
void your_sort(unsigned int* const d_inputVals,
unsigned int* const d_inputPos,
unsigned int* const d_outputVals,
unsigned int* const d_outputPos,
size_t numElems)
{
...
checkCudaErrors(cudaMemcpy(h_hist, d_hist, 2 * sizeof(unsigned int), cudaMemcpyDeviceToHost));
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
move<<<n_blocks, threads_per_block>>>(d_inputVals, d_inputPos, d_outputVals, d_outputPos,
d_scanned, numElems, i, h_hist[0]);
cudaDeviceSynchronize(); checkCudaErrors(cudaGetLastError());
checkCudaErrors(cudaMemcpy(d_inputVals, d_outputVals, numElems * sizeof(unsigned int), cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemcpy(d_inputPos, d_outputPos, numElems * sizeof(unsigned int), cudaMemcpyDeviceToDevice));
}
checkCudaErrors(cudaFree(d_hist));
checkCudaErrors(cudaFree(d_scanned));
}
Testez votre code. Pour le test final, il faudra exécuter ./HW4 red_eye_effect_5.jpg red_eye_effect_template_5.jpg
(après avoir fait un make). Par contre, vous pouvez modifier l'entrée de your_sort pour y mettre un tableau polus simple
et tester vos algorithmes. Par exemple, vous pouvez mettre au début du your_sort:
Votre sortie est écrite dans le fichier HW4_output.png;
unsigned int threads_per_block = 1024;
numElems = threads_per_block * 3 + 3;
unsigned int h_in[numElems];
for (int i = 0; i < numElems; i++) h_in[i] = numElems - i - 1;
cudaMemcpy(d_inputVals, h_in, numElems * sizeof(unsigned int), cudaMemcpyHostToDevice);
