Laplacian Filter in CUDA & OpenMp

Luca Sorace 1910722
Lorenzo Passaretti 1910723

Il laplacian filter è un filtro capace di fare tone mapping e di aumentare o diminuire i dettagli di un’immagine senza creare né aloni né distorsioni. È basato sulle piramidi laplaciane e gaussiane: crea una sottopiramide molto piccola dell’immagine originale ed applica la correzione solamente a quella piccola porzione.

L’algoritmo è diviso in tre grandi fasi:

Fase 1

Fase 2

Per ogni pixel G0, per ogni livello in inputGauss:

Fase 3

Strutture Dati

Image

Abbiamo due tipi di immagini: Image3 e Image4, le quali differiscono dal tipo di pixel che contengono (Pixel3 per le immagini RGB, o Pixel4 per le immagini RGBA). I due tipi di immagini sono definiti come segue:

typedef struct { uint32_t height; //altezza dell’immagine. uint32_t width; //larghezza dell’immagine. Pixel *pixels; //buffer dei pixel dell’immagine. } Image;

Vects

Abbiamo molti tipi di vettori che variano in base al tipo ed alla quantità di variabili che contengono. Due esempi sono:

typedef struct { float x; float y; float z; float w; } Vec4f; typedef struct { uint8_t x; uint8_t y; uint8_t z; uint8_t w; } Vec4u8;

In generale, la naming convention è la seguente:

Vec          -> è un vettore
[3|4]        -> numero di variabili del vettore
[u|i|f]      -> tipo di variabili del vettore (unsigned int, int, float)
[8|16|32|64] -> dimensione in bit delle variabili del vettore

WorkingBuffers

I working buffers sono delle strutture dati preallocate che contengono tutti i buffer utilizzati all’interno del rendering, in modo che non debbano essere riallocati ad ogni iterazione.

Ogni implementazione di LLF ha i suoi WorkingBuffer:

Single core (llf)

typedef struct { Kernel filter; //Kernel per il blur già precalcolato Pyramid gaussPyramid; //Piramide gaussiana che viene creata nella fase 1 dell'algoritmo Pyramid outputLaplacian; //Piramide laplaciana da collassare che viene calcolata nella fase 2 dell'algoritmo Pyramid bufferGaussPyramid; //Piramide gaussiana temporanea utilizzata all'interno della fase 2 dell'algoritmo Pyramid bufferLaplacianPyramid; //Piramide laplaciana temporanea utilizzata all'interno della fase 2 dell'algoritmo } WorkingBuffers;

OpenMP

typedef struct { uint32_t end; //Numbero di pixel totali presenti dentro gaussPyramid, ovvero il numero di pixel totali da renderizzare nella fase 2 Kernel filter; //Kernel per il blur già precalcolato Pyramid *bArr; //Array di piramidi gaussiane temporanee utilizzate, una da ogni thread differente, all'interno della fase 2 dell'algoritmo uint32_t *pyrDimensions; //Numero di pixel precalcolato per ogni livello della gaussPyramid, ovvero il numero di pixel per renderizzare ogni singolo livello differente della gaussPyramid nella fase 2 CurrentLevelInfo *cliArr; //Array di CurrentLevelInfo (vedi più avanti nella sezione dedicata), uno per ogni thread Pyramid gaussPyramid; //Piramide gaussiana che viene creata nella fase 1 dell'algoritmo Pyramid outputLaplacian; //Piramide laplaciana da collassare che viene calcolata nella fase 2 dell'algoritmo } WorkingBuffers;

CUDA

typedef struct { Kernel d_filter; //Kernel per il blur già precalcolato salvato sulla memoria globale del device Pyramid d_gaussPyramid; //Piramide gaussiana salvata sulla memoria globale del device che viene creata nella fase 1 dell'algoritmo Pyramid d_outputLaplacian; //Piramide laplaciana da collassare, salvata sulla memoria globale del device, che viene calcolata nella fase 2 dell'algoritmo Image3 *d_img; //Immagine di buffer salvata sulla memoria globale del device usata come input/output per la parte cuda dell'algoritmo } WorkingBuffers;

CurrentLevelInfo

Questa struttura è utilizzata solo per l’implementazione in openmp.
Contiene informazioni sul livello corrente nella fase 2 dell’algoritmo, assieme ad altri dati in cache

typedef struct { uint8_t currentNLevels; //Numero di livelli presenti nelle piramidi buffer che utilizziamo al livello corrente durante la fase 2 uint8_t lev; //Numero del livello corrente durante la fase 2 uint32_t oldY; //coordinata Y attuale nell'elaborazione del livello corrente durante la fase 2 uint32_t width; //larghezza del livello corrente durante la fase 2 uint32_t prevLevelDimension; //Numero di pixel che abbiamo renderizzato fino all'inizio del livello corrente durante la fase 2 uint32_t nextLevelDimension; //Numero di pixel che avremo renderizzato alla fine del livello corrente durante la fase 2. Viene usato per capire quando si deve passare all'elaborazione del livello successivo uint32_t subregionDimension; //Dimensioni della sottoregione R0 nel livello corrente durante la fase 2 uint32_t full_res_roi_yShifted; //Valore precalcolato di full_res_roi_yShifted uint32_t base_y; //Valore precalcolato di base_y uint32_t end_y; //Valore precalcolato di end_y Image3 *currentGaussLevel; //Puntatore al livello corrente della piramide gaussiana in input } CurrentLevelInfo;

Typedef

Pyramid

typedef Image3** Pyramid;

La piramide è un array di immagini le cui dimensioni si riducono quadraticamente ad ogni livello. Il livello zero è quello più grande.

lev[0].width = imgSource.width
lev[n].width = lev[n - 1].width / 2

Stessa cosa per l'height

Esistono due tipi di piramidi, definiti come segue:

Gaussiana

gauss[0] = sourceImg
gauss[n] = downsample(gauss[n - 1])

Dove

downsample(img) è una funzione che dimezza e applica un kernel di blur sull’input.

Laplaciana

lap[nLevels] = gauss[nLevels]
lap[n] = gauss[n] - upsample(gauss[n + 1])

Dove

upsample(img) è una funzione che duplica e applica un kernel di blur all’input.

Pixel

typedef Vec4f Pixel4;

I pixel sono dei vettori, in cui xyzw rappresentano i valori RGBA

AlphaMap

typedef uint8_t* AlphaMap;

Una alpha map è una matrice di uint8_t in cui in ogni posizione viene salvato il canale alpha del corrispettivo pixel di una Image3

Kernel

typedef float* Kernel;

Un Kernel è una matrice di float con valori compresi tra 0 ed 1, utilizzata nella funzione convolve(image, kernel), definita come segue:

Per ogni pixel C dell’immagine in input I:

All’interno dell’algoritmo LLF il kernel viene utilizzato per applicare un filtro di blur nelle fasi di upsampling e downsampling, usate per costruire le piramidi gaussiane e laplaciane

Singlecore

La versione singlecore utilizza l’algoritmo descritto inizialmente, ma con la piccola ottimizzazione dove downsampling e upsampling vengono svolte contemporaneamente al convolving. In questo modo risparmiamo tempo evitando copie extra delle immagini e non abbiamo neanche bisogno di ulteriori buffer durante le fasi intermedie di upsample e downsample.
Le funzioni originali sono comunque riportate nello stesso file.

Risultato dell’implementazione single core

Performance

Image Size Time in ms
480x320 5815
800x533 16612
1920x1279 98914
3840x2558 384384

Misure in ms. Test svolti con processore AMD r5 5600x

OpenMP

Idea algoritmo

Durante il rendering del filtro LLF non abbiamo dipendenze tra i singoli pixel che stiamo renderizzando, bensì abbiamo dipendenze tra i differenti layer delle piramidi.

Durante la fase 1 e 3 dell’algoritmo andiamo a parallelizzare (durante le operazioni di creazione della piramide gaussiana e di collaso di quella laplaciana) solamente il rendering di ogni singolo layer: ogni thread lavorerà sul suo set di pixel indipendentemente dagli altri e tra un layer e l’altro avverrà una sincronizzazione.

Durante la fase 2, invece, non esiste dipendenza né tra layer né tra pixel, e quindi la parallelizzazione avviene renderizzando più pixel contemporaneamente senza alcun bisogno di sincronizzazione. Il rendering di ciascun pixel avverrà analogalmente alla precedente implementazione single core. Ogni thread avrà una propria struttura dati contenente i suoi buffer personali per il rendering, più informazioni e dati “in cache” sul livello corrente che sta renderizzando. Lavorerà indipendentemente da tutti gli altri thread.

All’interno della fase 2 viene applicata una prima ottimizzazione sostanziale a livello di algoritmo: poiché dell’intera piramide laplaciana temporanea ci interessa solamente un pixel (quello che poi verrà copiato nella piramide di output), al posto di renderizzarla tutta, verrá upsampleata e blurrata solo la sezione che lo determina.

Architettura

Rispetto alla versione single core, vengono effettuate delle operazioni in più:

FASE 1:

FASE 2:
Viene lanciata un’unica parallel for che itererà su tutti i pixel di tutti i livelli della piramide gaussiana in input. Poiché solo in questo caso non c’è dipendenza tra i layer della piramide, non serviranno primitive di sincronizzazione tra l’uno e l’altro.

Per ogni pixel:

FASE 3:

Scelte progettuali

La creazione in fase progettuale della struttura dati currentLevelInfo viene dall’idea di voler ridurre al minimo l’overhead della creazione di nuovi thread. Essa è stata creata proprio per fare in modo che ci sia un’unica parallel for, piuttosto che una per ogni layer (che invece sarebbe stata una scelta più “pulita” a livello di codice). La creazione di una tale struttura dati è stata possibile grazie al fatto che non esiste alcun tipo di dipendenza tra pixel nell’elaborazione della piramide laplaciana di output. Questa caratteristica, inoltre, ha semplificato di molto il processo di parallelizzazione e della scrittura di codice, permettendoci di applicare la scelta più ovvia: al posto di parallelizzare il lavoro per renderizzare un singolo pixel della piramide gaussiana di input, lavoriamo su più pixel contemporaneamente. In questo modo, tra l’altro, siamo sempre sicuri che il carico di ogni thread sarà equo rispetto agli altri

Test di correttezza

Abbiamo testato l’algoritmo in openmp con vari parametri e abbiamo sempre riscontrato risultati visibilmente identici alla versione single core.

Risultato dell’implementazione con openMP

Performance

Image Size 24 Threads 16 Threads 12 Threads 8 Threads 4 Threads 2 Threads Single core (Con stessa implementazione)
480x320 70 87 146 168 279 443 872
800x533 205 257 394 473 826 1301 2551
1920x1279 1213 1528 2271 2845 4881 7779 15344
3840x2558 4904 6226 9266 11338 18299 30788 60654

Misure in ms. Test svolti con processore AMD r9 5900x

Possiamo notare come, anche con input piccoli, lo speedup e l’efficenza rimangano costanti a parità di thread.
L’efficenza inversamente proporzionale al numero di thread indica che la gestione della parallelizzazione fa da collo di bottiglia alla velocità di esecuzione del codice

Speedup

Image Size 24 Threads 16 Threads 12 Threads 8 Threads 4 Threads 2 Threads
480x320 12,4571 10,023 5,9726 5,1905 3,1254 1,9684
800x533 12,4439 9,9261 6,4746 5,3932 3,0884 1,9608
1920x1279 12,6496 10,0419 6,7565 5,3933 3,1436 1,9725
3840x2558 12,3683 9,7420 6,5459 5,3496 3,3146 1,9701

Efficienza

Image Size 24 Threads 16 Threads 12 Threads 8 Threads 4 Threads 2 Threads
480x320 0,5190 0,6264 0,4977 0,6488 0,7814 0,9842
800x533 0,5185 0,6204 0,5396 0,6742 0,7721 0,9804
1920x1279 0,5271 0,6276 0,5630 0,6742 0,7859 0,9862
3840x2558 0,5153 0,6089 0,5455 0,6687 0,8287 0,9850

Cuda

Idea algoritmo

Durante la fase 1 e 3 dell’algoritmo i kernel cuda vengono lanciati con un solo blocco e più thread (La motivazione è nelle scelte progettuali) che elaboreranno più pixel, creando così una associazione 1:N tra thread e pixel da elaborare. Per mantenere le dipendenze tra layer, si usa una __syncthreads() alla fine di ogni upsample o downsample. L’unica funzione che fa eccezione è d_clampImage3(), dove invece c’è una associazione 1:1 tra thread e pixel da renderizzare

Durante la fase 2, invece, useremo più blocchi. Ciascun blocco si occuperà di renderizzare più pixel della piramide gaussiana di input, e ciascun thread del blocco si occuperà di upsampleare o downsampleare più pixel della relativa sottoregione. Creiamo così una associazione 1:N tra blocco e pixel della piramide gaussiana di input da renderizzare, ed un’altra associazione 1:N tra ogni thread (di un blocco) e il pixel della sottoregione da renderizzare. Così facendo stiamo applicando una doppia parallelizzazione: Renderizziamo più pixel della piramide gaussiana di input contemporaneamente, e parallelizziamo il rendering di ciascuno di essi.

Durante l’elaborazione del singolo pixel applichiamo una ottimizzazione a livello spaziale del miglioramento sostanziale applicato già nel codice relativo ad OpenMP (Quarto paragrafo dell’“Idea algoritmo”): quando calcoliamo la piramide gaussiana temporanea, ci salviamo solamente gli ultimi due layer più piccoli (Quelli utilizzati successivamente dalla versione già ottimizzata della piramide laplaciana). In questo modo possiamo avere entrambi i layer in shared memory ed ottenere così accessi molto più rapidi in lettura e scrittura durante i vari punti della fase 2, raggiungendo così uno speedup di circa il 60%

A differenza di openmp, in cuda siamo costretti a chiamare più volte il kernel della llf, una volta per ogni layer della piramide gaussiana di input

Architettura

Rispetto all’implementazione single core, sono state riscritte molte funzioni di utility generale appositamente per cuda. In particolare:

Inoltre sono state fatte ulteriori modifiche alle 3 fasi dell’algoritmo:

FASE 1:

FASE 2:
Nella fase 2 l’host looppa per ogni layer della piramide gaussiana di input. Per ciascuno precalcola la dimensione della sottoregione, e chiama il kernel interno sul device utilizzando N blocchi ed M thread, dove N = numero di pixel della piramide gaussiana di input che verranno elaborati contemporaneamente e M = numero di thread usati contemporaneamente per renderizzare un singolo pixel della piramide gaussiana di input
Ciascun blocco:

FASE 3:

Scelte progettuali

Il motivo per cui nelle fasi 1 e 3 lanciamo i kernel (che elaborano piramidi) con un solo blocco è in realtà molto semplice: poiché la dimensione delle piramidi si riduce quadraticamente ci si ritrova sempre con più thread che pixel da elaborare, rendendo inutile il lancio e sincronizzazione di più blocchi. Inoltre, poiché c’è una dipendenza tra i vari layer della piramidi, avremmo dovuto lanciare ogni volta un kernel per elaborare un singolo layer per mantenere la sincronizzazione necessaria, invece di avviarne soltanto uno all’inizio.
d_clampImage3() invece viene lanciato con più blocchi perché lavora su tutti i pixel dell’immagine di partenza, ben più dei thread possibili, ed in più non lavora con layer differenti dove sarebbero state necessarie primitive di sincronizzazione

La ragione per cui, a differenza di openmp, chiamiamo una volta per layer il kernel interno della llf, è perché avere thread e blocchi a lavorare contemporaneamente su layer differenti aumenterebbe la branch divergence: le “grandezze” dei for sono tutte uguali solo all’interno di uno stesso layer. L’esempio più immediato è applicare il downsample sulla sottoregione: su layer differenti la sottoregione avrà dimensioni differenti, di conseguenza i for all’interno del downsample avranno “durate” differenti.

Nella fase 2 NON facciamo una associazione 1:1 tra blocco e pixel da renderizzare per diminuire i tempi di accesso in memoria: Ad inizio del kernel interno della llf viene effettuata una copia dalla global alla shared memory del kernel di blur, in modo che poi ci si possa accedere più rapidamente nelle fasi di upsampling e downsampling. Avere un numero di blocchi troppo alto, vorrebbe dire effettuare questa copia troppe volte, rallentando di fatto l’esecuzione. Il numero di blocchi ideale corrisponde di conseguenza al numero di blocchi massimo che il device riesce a mandare in esecuzione contemporaneamente.

Test di correttezza

Risultato dell’implementazione con CUDA

Performance

Image Size 1Block 32Threads 16Block 64Threads 32Block 128Threads 256Block 256Threads 512Block 256Threads 512Block 512Threads OpenMP 1thread
480x320 1338 73 32 21 17 30 872
800x533 3749 201 83 52 51 84 2551
1920x1279 22331 1091 452 279 285 461 15344
3840x2558 95493 4431 1754 1112 1045 1763 60654

Misure in ms. Test svolti con gpu RTX 3600

Possiamo notare che, come accadeva esattamente con openmp, lo speedup rimane costante a parità di blocchi/threads.
Notiamo anche però come, superato una soglia di blocchi, lo speedup aumenti di pochissimo, poiché avremo più blocchi rispetto a quanti la gpu riesca a eseguirne contemporaneamente e poiché dovremmo effettuare più copie del kernel di blur dalla global memory alla shared memory.
Aumentare di troppo il numero di thread, invece, è proprio controproducente. Questo perché con troppi thread abbiamo un numero esagerato di thread “a non fare niente” per le elaborazioni delle sottoregioni più piccole, ed inoltre riduciamo anche il numero di blocchi che verranno eseguiti contemporaneamente

Speedup

Image Size 1Block 32Threads 16Block 64Threads 32Block 128Threads 256Block 256Threads 512Block 256Threads 512Block 512Threads
480x320 0,6517 11,9452 27,2500 41,5238 51,2941 29,0667
800x533 0,6804 12,6915 30,7349 49,0577 50,0196 30,3690
1920x1279 0,6871 14,0642 33,9469 54,9964 53,8386 33,2842
3840x2558 0,6352 13,6886 34,5804 54,5450 58,0421 34,4039

Implementazione real time con NDI®

Utilizzando NDI® (Nuovo standard industriale per la trasmissione real time di video a bassa compressione over IP, successore di SDI) abbiamo creato un’applicazione capace di ricevere una stream video e di inoltrare il traffico renderizzandoci sopra il filtro llf in tempo reale. La codebase è unica per tutte e tre le versioni; vengono usate flag del preprocessore per chiamare le funzioni corrette.

NDI® is a registered trademark of Vizrt Group.

Algoritmo

Le prime cose che accadono all’avvio del programma sono il caricamento degli argomenti, il set dell’handler per il sigint e l’inizializzazione della libreria NDI®.
Dopo aver trovato la source NDI® col nome corrispondente all’apposito argomento, ci si connette specificando di voler utilizzare RGBX/A come formato colore (Di default viene usato UYVY), si crea il receiver ed il sender NDI®, e si avvia un secondo thread che si occuperà del rendering vero e proprio.

Il secondo thread inizialmente si occuperà di allocare momentaneamente “al minimo” i suoi buffer interni, i semafori ed inizializzare i working buffer della versione di LLF selezionata. A questo punto il programma entra in uno stato di working, riassumibile così:

THREAD PRINCIPALE:

PROCESSING THREAD: