GPU-Bildverarbeitungsoptimierung am Beispiel der Medianfilterung

  • Tutorial

Einleitung


Grafikbeschleuniger (GPUs) sind seit langem für die Bild- und Videoverarbeitung entwickelt worden. Ab einem bestimmten Zeitpunkt wurden GPUs für allgemeine Zwecke eingesetzt. Aber auch die Entwicklung von Zentralprozessoren stand nicht still: Intel entwickelt aktiv Vektorerweiterungen (AVX256, AVX512, AVX1024). Als Ergebnis erscheinen verschiedene Prozessoren - Core, Xeon, Xeon Phi. Die Bildverarbeitung kann einer solchen Klasse von Algorithmen zugeordnet werden, die leicht vektorisiert werden können.
Die Praxis zeigt jedoch, dass es trotz des hohen Compiler-Niveaus und der Anpassungsfähigkeit der Zentral- und Coprozessoren von Xeon Phi nicht so einfach ist, Bilder mit Vektoranweisungen zu verarbeiten, da moderne Compiler mit automatischer Vektorisierung schlecht zurechtkommen und die Verwendung von vektoreigenen Funktionen recht schwierig ist. Es stellt sich auch die Frage, manuell vektorisierte Code- und Skalarabschnitte zu kombinieren.


In diesem Artikel wurde untersucht, wie mit Vektorerweiterungen die Medianfilterung für Zentraleinheiten optimiert werden kann. Für die Reinheit des Experiments habe ich dieselben Bedingungen wie der Autor verwendet: Ich habe einen 3x3-Quadrat-Median-Filteralgorithmus verwendet. Diese Bilder belegen 8 Bit pro Pixel. Die Frame-Verarbeitung (abhängig von der Größe des Filterquadrats) wird nicht durchgeführt. und versuchte, diesen Kern auf der GPU zu optimieren. Ein Algorithmus zur Medianfilterung von 5x5-Quadraten wurde ebenfalls in Betracht gezogen. Das Wesen des Algorithmus ist wie folgt:


  • Für jeden Punkt des Quellbildes wird eine bestimmte Nachbarschaft aufgenommen (in unserem Fall 3x3 / 5x5).
  • Die Punkte in dieser Nachbarschaft werden nach zunehmender Helligkeit sortiert.
  • Der Mittelpunkt der sortierten Nachbarschaft wird im endgültigen Bild aufgezeichnet.

Erste GPU-Kernel-Implementierung


Um keine "unnötigen" Optimierungen vorzunehmen, verwenden wir einige Berechnungen des Autors des genannten Artikels. Am Ende erhalten Sie einen einfachen Code für den Rechenkern:


__device__ __inline__ void Sort(int &a, int &b)
{
    const int d = a - b;
    const int m = ~(d >> 8);
    b += d & m;
    a -= d & m;
}
__global__ void mFilter(const int H, const int W, const unsigned char *in, unsigned char *out)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x + 1;
    int idy = threadIdx.y + blockIdx.y * blockDim.y + 1;
    int a[9] = {};
    if (idy < H - 1 && idx < W - 1)
    {        
        for (int z2 = 0; z2 < 3; ++z2)
            for (int z1 = 0; z1 < 3; ++z1)
                a[3 * z2 + z1] = in[(idy - 1 + z2) * W + idx - 1 + z1];
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
        Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
        Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
        Sort(a[4], a[2]);
        out[idy * W + idx] = (unsigned char)a[4];
    }
}

Für jeden Punkt wird die Nachbarschaft gemäß diesem Bild geladen (das Bildpixel, für das das Ergebnis berechnet werden soll, ist gelb hervorgehoben):



In dieser Implementierung gibt es ein Minus - das Bild wird im 8-Bit-Format gespeichert und die Berechnungen erfolgen in 32-Bit. Wenn Sie sich den Kampfcode einer Vektorimplementierung auf der CPU ansehen, sehen Sie, dass viele Optimierungen angewendet wurden, einschließlich des Ausschlusses unnötiger Operationen während des Sortierens, da je nach Compiler Optimierungen möglicherweise nicht funktionieren. Der Nvidia-Compiler hat keine derartigen Nachteile - er kann problemlos zwei Schleifen in z2 und z1 implementieren, selbst wenn #pragma unroll nicht angegeben ist (in diesem Fall entscheidet der Compiler, die Schleife selbst zu implementieren, und wenn nicht genügend Register vorhanden sind, werden die Schleifen nicht implementiert). und auch unnötige Operationen von Sortierfunktionen nach deren Ersetzung werden ausgeschlossen ( Inline - Anzeige)ist eine Empfehlung, wenn Sie die Funktion unabhängig davon , wie der Compiler entscheidet , es ist notwendig , um den Einsatz ersetzen wollen forceinline ).
Ich werde auf einer Nvidia GTX 780 Ti-GPU testen, die ungefähr im selben Veröffentlichungsjahr wie der Zentralprozessor im obigen Artikel (Intel Core-i7 4770, 3,4 GHz) veröffentlicht wurde. Bildgröße - Schwarzweiß 1920 x 1080 Pixel. Wenn Sie die kompilierte erste Version auf einer GPU ausführen, erhalten wir die folgenden Zeiten:


Vorlaufzeit, msRelative Beschleunigung
Skalare CPUAVX2-CPUGPU Version1Skalar / AVX2Skalar / Version1 AVX2 / Version1
24.8140,4240,37258,566,71.13

Wie Sie aus den oben genannten Zeiten ersehen können, erhalten wir eine gute Beschleunigung, wenn Sie den Code „nativ“ umschreiben und auf der GPU ausführen. Es gibt keine Probleme, das Bild an der Größe auszurichten - der Code auf der GPU funktioniert bei jeder Bildgröße (bzw. ab einer bestimmten Größe) gleich gut wenn die GPU voll geladen ist). Sie müssen auch nicht über die Speicherausrichtung nachdenken - der Compiler wird alles für uns tun, es reicht aus, Speicher mit den speziellen CUDA-RunTime-Funktionen zuzuweisen.
Viele, die auf Grafikbeschleuniger gestoßen sind, werden sofort sagen, dass es notwendig ist, Daten herunterzuladen und wieder hochzuladen. Mit einer PCIe 3.0-Geschwindigkeit von 15 GB / s für 2 MB Bild werden etwa 130 Mikrosekunden für die Einwegübertragung erzielt. Wenn also ein Bild-Streaming stattfindet, dauert das Laden und Entladen ca. 260 Mikrosekunden, während das Konto ca. 372 Mikrosekunden benötigt. Dies bedeutet, dass selbst in der „nativsten“ Implementierung Übertragungen an die GPU abgedeckt sind. Wir können auch den Schluss ziehen, dass die Verarbeitung des Bilderstroms schneller als 260 Mikrosekunden pro Bild fehlschlägt. Wenn dieser Filter jedoch zusammen mit weiteren zehn Filtern während der Verarbeitung eines Bildes verwendet wird, ist eine weitere Optimierung dieses Kernels nützlich.


Zweite GPU-Core-Implementierung


Um die Effizienz der Sortiervorgänge zu erhöhen, können Sie die kürzlich hinzugefügten SIMD-Anweisungen verwenden. Diese Anweisungen sind intrinsische Funktionen, mit denen Sie in einer Operation eine mathematische Funktion für zwei nicht signierte / signierte Kurzzeichen oder für vier nicht signierte / signierte Zeichen berechnen können. Eine Liste der verfügbaren Funktionen finden Sie in der Cuda ToolKit- Dokumentation . Wir brauchen zwei Funktionen: vminu4 (ohne Vorzeichen a, ohne Vorzeichen b) und vmaxu4 (ohne Vorzeichen a, ohne Vorzeichen b); Berechnen des Minimums und Maximums für 4x ohne vorzeichenbehaftete 8-Bit-Werte. Um vier Werte zu laden, teilen Sie das gesamte Bild horizontal in vier Teile. Jeder Block lädt Daten aus vier Bändern und kombiniert vier geladene Punkte mit einer Bitverschiebung im vorzeichenlosen int:



Die aktive Verarbeitungszone ist das oberste Quadrat mit einer Größe von 1920 x 270. Dementsprechend wird ein Block mit 128 x 1 Threads in Orange angezeigt, in den 128 * 4 Elemente geladen werden. Wie Sie wissen, besteht die Größe einer Verzerrung auf einer GPU aus 32 Threads. Alle Threads eines Warps führen jeweils einen Befehl pro Zyklus aus. Da jeder Thread 4 Elemente des Typs "unsigned char" hochlädt, führt ein Warp in einem Befehl eine Verarbeitung von insgesamt 128 Elementen aus, dh es wird eine Verarbeitung von 1024 Bit erhalten. Für die CPU stehen derzeit 256-Bit-AVX2-Register zur Verfügung (die sich in Zukunft in 512 und 1024 Bit verwandeln werden). Mit diesen Optimierungen wird der Code des Rechenkerns in die folgende Form konvertiert:


__device__ __inline__ void Sort(unsigned int &a, unsigned int &b)
{
    unsigned t = a;
    a = __vminu4(t, b);
    b = __vmaxu4(t, b);
}
__global__ __launch_bounds__(256, 8) 
void mFilter_3(const int H, const int W, const unsigned char * __restrict in, unsigned char *out, const unsigned stride_Y)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x + 1;
    int idy = threadIdx.y + blockIdx.y * blockDim.y + 1;
    unsigned a[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
    if (idx < W - 1)
    {  
        for (int z3 = 0; z3 < 4; ++z3)
        {
            if (idy < H - 1)
            {
                const int shift = 8 * (3 - z3);
                for (int z2 = 0; z2 < 3; ++z2)
                    for (int z1 = 0; z1 < 3; ++z1)            
                        a[3 * z2 + z1] |= in[(idy - 1 + z2) * W + idx - 1 + z1] << shift;               
            }
            idy += stride_Y;
        }      
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
        Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
        Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
        Sort(a[4], a[2]);
        idy = threadIdx.y + blockIdx.y * blockDim.y + 1;
        for (int z = 0; z < 4; ++z)
        {
            if (idy < H - 1)
                out[idy * W + idx] = (unsigned char)((a[4] >> (8 * (3 - z))) & 255);
            idy += stride_Y;
        }
    }
}

Infolge kleiner Änderungen haben wir die Anzahl der verarbeiteten Elemente auf 4x erhöht (oder 128 Elemente pro Anweisung). Die Komplexität des Codes hat zugenommen und der Compiler generiert zu viele Register, sodass die GPU nicht zu 100% geladen werden kann. Da wir noch einfache Berechnungen haben, teilen wir dem Compiler mit, dass wir die maximal mögliche Anzahl von Threads auf jedem Stream-Prozessor ausführen möchten (die maximale Anzahl für diese GPU beträgt 2048 Threads), wobei __launch_bounds (128, 16) verwendet werden Eine Konfiguration von 16 Blöcken mit jeweils 128 Threads wurde gestartet. Eine weitere Neuerung: Wir fügen das Schlüsselwort "restricted" hinzu, um den Textur-Cache zum Laden doppelter Daten zu verwenden (wenn Sie analysieren, gibt es viele, die meisten Punkte werden dreimal gelesen). Damit der Textur-Cache aktiviert wird, müssen Sie ebenfalls const angeben.


const unsigned char * __restrict in.

Relevante Zeiten der zweiten Version:


Vorlaufzeit, msRelative Beschleunigung
Skalare CPUAVX2-CPUGPU Version1Version2 GPUSkalar / AVX2Skalar / Version2 AVX2 / Version2 Version1 / Version2
24.8140,4240,3720,20758,5119,82,041,79

Dritte GPU-Core-Implementierung


Bereits erhalten keine schlechte Beschleunigung. Sie können jedoch feststellen, dass beim Laden von 4 Elementen mit dieser Entfernung zu viele Downloads wiederholt werden. Dies verringert die Leistung erheblich. Gibt es eine Möglichkeit, das, was wir bereits direkt in die Register geladen haben, zu verwenden und die Lokalität der Daten zu verbessern? Ja das kannst du Lassen Sie diese 4 Zeilen nun nahe beieinander liegen, und der Anfang der Blöcke wird über die gesamte Höhe mit einem Abstand von 4 "verschmiert". Angesichts der Größe des Fadenblocks entlang Y erhalten wir also den folgenden Code:


#define BL_X 128
#define BL_Y 1
__device__ __inline__ void Sort(unsigned int &a, unsigned int &b)
{
    unsigned t = a;
    a = __vminu4(t, b);
    b = __vmaxu4(t, b);
}
__global__ __launch_bounds__(256, 8) 
void mFilter_3a(const int H, const int W, const unsigned char * __restrict in, unsigned char *out)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x + 1;
    int idy = threadIdx.y + blockIdx.y * blockDim.y * 4 + 1;
    unsigned a[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
    if (idx < W - 1)
    {        
        for (int z3 = 0; z3 < 4; ++z3)
        {
            //if (idy < H - 1) <---- this is barrier for compiler optimization
            {
                const int shift = 8 * (3 - z3);
                for (int z2 = 0; z2 < 3; ++z2)
                    for (int z1 = 0; z1 < 3; ++z1)            
                        a[3 * z2 + z1] |= in[(idy - 1 + z2) * W + idx - 1 + z1] << shift;               
            }
            idy += BL_Y;
        }
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
        Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
        Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
        Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
        Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
        Sort(a[4], a[2]);
        idy = threadIdx.y + blockIdx.y * blockDim.y * 4 + 1;
        for (int z = 0; z < 4; ++z)
        {
            if (idy < H - 1)
                out[idy * W + idx] = (unsigned char)((a[4] >> (8 * (3 - z))) & 255);
            idy += BL_Y;
        }
    }
}

Durch diesen optimierten Code kann der Compiler die bereits geladenen Daten verwenden. Das einzige Problem ist die Überprüfung des Zustands (idy <H - 1). Um dieses Häkchen zu entfernen, müssen Sie die erforderliche Anzahl von Zeilen hinzufügen, dh das Bild vertikal um 4 ausrichten und alle zusätzlichen Elemente mit der maximalen Anzahl (255) ausfüllen. Nach der Konvertierung des Codes in die dritte Option lauten die Zeiten wie folgt:


Vorlaufzeit, msRelative Beschleunigung
Skalare CPUAVX2-CPUGPU Version1Version2 GPUVersion2 GPUSkalar / AVX2Skalar / Version3 AVX2 / Version3
24.8140,4240,3720,2070,13058,5190,83.26

Darüber hinaus ist es bereits recht schwierig zu optimieren, da es nicht so viele Berechnungen gibt. Das Problem der Verwendung von Shared Memory wurde hier nicht angesprochen. Es gab einige Versuche, Shared Memory zu verwenden, aber ich konnte nicht die beste Zeit bekommen. Wenn wir die Bildverarbeitungsgeschwindigkeit mit seiner Größe verknüpfen, können Sie unabhängig davon, ob das Bild in der Größe ausgerichtet ist oder nicht, ungefähr 15 Gigapixel / Sekunde erhalten.


Der Kern des 5x5-Medianfilters auf der GPU


Wenn Sie alle Optimierungen anwenden, können Sie einen Filter mit einem 5x5-Quadrat mit geringfügigen Änderungen implementieren. Die Komplexität dieses Codes liegt in der korrekten Implementierung einer partiellen bitonischen Sortierung für 25 Elemente. Analog dazu können Sie Filter beliebiger Größe implementieren (beginnend mit einer bestimmten Filtergröße wird der allgemeine Sortieralgorithmus wahrscheinlich schnell sein).


__global__ __launch_bounds__(256, 8) 
void mFilter5a(const int H, const int W, const unsigned char * in, unsigned char *out, const unsigned stride_Y)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x + 2;
    int idy = 4 * threadIdx.y + blockIdx.y * blockDim.y * 4 + 2;
    if (idx < W - 2)
    {
        unsigned a[25] = { 0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0, 0, 0 };
        for (int t = 0; t < 4; ++t)
        {
            const int shift = 8 * (3 - t);
            for (int q = 0; q < 5; ++q)
                for (int z = 0; z < 5; ++z)
                    a[q * 5 + z] |= (in[(idy - 2 + q) * W + idx - 2 + z] << shift);
            idy++;
        }
        Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[2], a[4]); Sort(a[2], a[3]); Sort(a[6], a[7]);
        Sort(a[5], a[7]); Sort(a[5], a[6]); Sort(a[9], a[10]); Sort(a[8], a[10]); Sort(a[8], a[9]);
        Sort(a[12], a[13]); Sort(a[11], a[13]); Sort(a[11], a[12]); Sort(a[15], a[16]); Sort(a[14], a[16]);
        Sort(a[14], a[15]); Sort(a[18], a[19]); Sort(a[17], a[19]); Sort(a[17], a[18]); Sort(a[21], a[22]);
        Sort(a[20], a[22]); Sort(a[20], a[21]); Sort(a[23], a[24]); Sort(a[2], a[5]); Sort(a[3], a[6]);
        Sort(a[0], a[6]); Sort(a[0], a[3]); Sort(a[4], a[7]); Sort(a[1], a[7]); Sort(a[1], a[4]);
        Sort(a[11], a[14]); Sort(a[8], a[14]); Sort(a[8], a[11]); Sort(a[12], a[15]); Sort(a[9], a[15]);
        Sort(a[9], a[12]); Sort(a[13], a[16]); Sort(a[10], a[16]); Sort(a[10], a[13]); Sort(a[20], a[23]);
        Sort(a[17], a[23]); Sort(a[17], a[20]); Sort(a[21], a[24]); Sort(a[18], a[24]); Sort(a[18], a[21]);
        Sort(a[19], a[22]); Sort(a[9], a[18]); Sort(a[0], a[18]); Sort(a[8], a[17]); Sort(a[0], a[9]);
        Sort(a[10], a[19]); Sort(a[1], a[19]); Sort(a[1], a[10]); Sort(a[11], a[20]); Sort(a[2], a[20]);
        Sort(a[12], a[21]); Sort(a[2], a[11]); Sort(a[3], a[21]); Sort(a[3], a[12]); Sort(a[13], a[22]);
        Sort(a[4], a[22]); Sort(a[4], a[13]); Sort(a[14], a[23]); Sort(a[5], a[23]); Sort(a[5], a[14]);
        Sort(a[15], a[24]); Sort(a[6], a[24]); Sort(a[6], a[15]); Sort(a[7], a[16]); Sort(a[7], a[19]);
        Sort(a[13], a[21]); Sort(a[15], a[23]); Sort(a[7], a[13]); Sort(a[7], a[15]); Sort(a[1], a[9]);
        Sort(a[3], a[11]); Sort(a[5], a[17]); Sort(a[11], a[17]); Sort(a[9], a[17]); Sort(a[4], a[10]);
        Sort(a[6], a[12]); Sort(a[7], a[14]); Sort(a[4], a[6]); Sort(a[4], a[7]); Sort(a[12], a[14]);
        Sort(a[10], a[14]); Sort(a[6], a[7]); Sort(a[10], a[12]); Sort(a[6], a[10]); Sort(a[6], a[17]);
        Sort(a[12], a[17]); Sort(a[7], a[17]); Sort(a[7], a[10]); Sort(a[12], a[18]); Sort(a[7], a[12]);
        Sort(a[10], a[18]); Sort(a[12], a[20]);  Sort(a[10], a[20]); Sort(a[10], a[12]);
        idy = 4 * threadIdx.y + blockIdx.y * blockDim.y * 4 + 2;
        for (int z = 0; z < 4; ++z)
        {
            if (idy < H - 2)
                out[idy * W + idx] = (unsigned char)((a[12] >> (8 * (3 - z))) & 255);
            idy++;
        }
    }
}

Die Verarbeitungsgeschwindigkeit eines Filters dieser Größe ist bei Verwendung von AVX2 ungefähr 5,7-mal schneller, während die GPU-Version nur 3,8-mal langsamer wird (ca. 0,5 ms). Die endgültige Beschleunigung bei der GPU mit einer mittleren Filterung von 5 × 5 Quadrat kann im Vergleich zur AVX2-Version das Fünffache erreichen.


Fazit


Abschließend können wir einige der Ergebnisse zusammenfassen. Durch die Optimierung des sequentiellen Skalarprogramms können Sie mit der GPU eine hervorragende Beschleunigung erzielen. Bei Medianfiltern können die Beschleunigungswerte je nach Radius des Filters das 1000-fache erreichen! Natürlich verringert eine sehr optimale Version, die AVX2-Anweisungen verwendet, die Lücke zwischen der GPU und der CPU dramatisch. Persönlich kann ich jedoch aus eigener Erfahrung den Schluss ziehen, dass die manuelle Vektorisierung des Codes nicht so einfach ist, wie es tatsächlich scheint. Im Gegenteil, selbst die "native" Implementierung der GPU führt zu relativ hohen Beschleunigungen. Und in den meisten Fällen reicht das aus. Und wenn Sie bedenken, dass die Anzahl der Grafikkarten in einem Motherboard mehr als eins sein kann, kann sich die Gesamtbeschleunigung bei der Verarbeitung des Bild- oder Videoflusses linear erhöhen.
Die Frage, ob die GPU verwendet werden soll oder nicht, ist rhetorisch. Natürlich zeigt die GPU bei der Verarbeitung von Bildern und Videos eine bessere Leistung, und mit zunehmender Rechenlast steigt der Unterschied zwischen der Verarbeitungsgeschwindigkeit auf der CPU und der GPU zugunsten der letzteren.


Kampfoptimierungen finden Sie im zweiten Artikel: Schnellere oder tiefere Optimierung der Medianfilterung für Nvidia-GPUs


Jetzt auch beliebt: