Schnellere oder tiefere Optimierung der Medianfilterung für Nvidia-GPUs

    Einleitung


    In einem früheren Beitrag habe ich versucht zu beschreiben, wie einfach es ist, die GPU für die Bildverarbeitung zu nutzen. Das Schicksal war für mich eine Gelegenheit, die Medianfilterung für die GPU zu verbessern. In diesem Beitrag werde ich versuchen, anhand des Beispiels der Medianfilterung zu erläutern, wie Sie die GPU bei der Bildverarbeitung noch leistungsfähiger machen können. Wir werden die GPU GTX 780 ti mit dem optimierten Code vergleichen , der auf dem modernen Intel Core i7 Skylake 4.0 GHz-Prozessor mit einem Satz von Vektorregistern AVX2 ausgeführt wird. Die erreichte Filtergeschwindigkeit von 3x3 Quadrat bei 51 GPixeln / Sek. Für die GPU GTX 780Ti und die spezifische Filtergeschwindigkeit von 3x3 Quadrat bei 10,2 GPixeln / Sek. Pro 1 TFlops für eine einzelne Genauigkeit zu diesem Zeitpunkt sinddas höchste bekannte in der Welt .


    Um zu verdeutlichen, was gerade passiert, wird dringend empfohlen, den vorherigen Beitrag zu lesen , da dieser Beitrag Tipps, Empfehlungen und Vorgehensweisen zur Verbesserung der bereits geschriebenen Version enthält. Um Sie daran zu erinnern, was wir tun:


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

    Medianfilter 3x3 Quadrat


    Eine Möglichkeit zur Optimierung des Bildverarbeitungscodes auf einer GPU besteht darin, wiederholte Vorgänge zu suchen und hervorzuheben und sie nur einmal in einem Warp auszuführen. In einigen Fällen müssen Sie den Algorithmus geringfügig ändern. Wenn Sie die Filterung in Betracht ziehen, können Sie feststellen, dass benachbarte Threads beim Laden eines 3x3-Quadrats um ein Pixel die gleichen logischen Operationen ausführen und den gleichen Speicher laden:


    __global__ __launch_bounds__(256, 8) 
    void mFilter_3х3(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)
            {
                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;  // <----- Здесь каждый warp делает лишние
                                                                                          // 6 операций загрузки и 12 логических 
                                                                                         // операций
                idy += BL_Y;
            }
            /*  Остальной код далее  */
        }
    }

    Mit ein wenig Überlegung können Sie die Sortierung verbessern, indem Sie dieselben Vorgänge nur einmal ausführen. Zur Verdeutlichung werde ich zeigen, was Sortieren war und was wurde und welche Vorteile dieser Ansatz bietet.
    In einem früheren Beitrag habe ich die bereits erfundene unvollständige Metzgersortierung für 9 Elemente verwendet, mit der Sie den Median ermitteln können, ohne bedingte Operatoren für die minimale Anzahl von Operationen zu verwenden. Wir stellen diese Sortierung in Form von Iterationen dar, nach denen jeweils unabhängige Elemente sortiert werden, wobei jedoch Abhängigkeiten zwischen den Iterationen bestehen. Für jede Iteration werden Paare von Elementen zum Sortieren ausgewählt.
    Alte Sorte:



    Neue Sortierung:



    Anscheinend können beide Ansätze als gleich bezeichnet werden. Betrachten wir dies jedoch unter dem Gesichtspunkt der Optimierung in Bezug auf die Anzahl der Operationen und die Möglichkeit, die berechneten Werte von anderen Threads zu verwenden. Sie können Spalten nur einmal laden und sortieren, zwischen Threads austauschen und den Median zählen. Das sieht dann optisch so aus:



    Der folgende Kernel entspricht einem solchen Algorithmus. Für den Austausch habe ich mich entschieden, keinen gemeinsamen Speicher zu verwenden, sondern neue __shfl-Operationen, mit denen Sie Werte innerhalb eines Warps ohne Synchronisation austauschen können. Somit ist es möglich, wiederholte Berechnungen teilweise zu vermeiden.


    __global__ BOUNDS(BL_X, 2048 / BL_X)
    void mFilter_3x3(const int H, const int W, const unsigned char * __restrict in, unsigned char *out)
    {
        // индекс по Оси Х
        int idx = threadIdx.x + blockIdx.x * blockDim.x;
        // смещаемся немного назад, так как в каждом warp'е первая и последняя нить не участвует в итоговом вычислении
        idx = idx - 2 * (idx / 32);
        // умножаем на 4, так как мы обрабатываем по 4 точки сразу
        int idy = threadIdx.y + blockIdx.y * blockDim.y * 4 + 1;
        unsigned a[9];
        unsigned RE[3];
        // первая и последняя нить не участвует в итоговом вычислении
        bool bound = (idxW > 0 && idxW < 31 && idx < W - 1);
        // считываем и сортируем столбцы
        RE[0] = in[(idy - 1) * W + idx] | (in[(idy + 0) * W + idx] << 8) | (in[(idy + 1) * W + idx] << 16) | (in[(idy + 2) * W + idx] << 24);
        RE[1] = in[(idy + 0) * W + idx] | (in[(idy + 1) * W + idx] << 8) | (in[(idy + 2) * W + idx] << 16) | (in[(idy + 3) * W + idx] << 24);
        RE[2] = in[(idy + 1) * W + idx] | (in[(idy + 2) * W + idx] << 8) | (in[(idy + 3) * W + idx] << 16) | (in[(idy + 4) * W + idx] << 24);
        Sort(RE[0], RE[1]);
        Sort(RE[1], RE[2]);
        Sort(RE[0], RE[1]);
        // делаем обмены 
        a[0] = __shfl_down(RE[0], 1);
        a[1] = RE[0];
        a[2] = __shfl_up(RE[0], 1);
        a[3] = __shfl_down(RE[1], 1);
        a[4] = RE[1];
        a[5] = __shfl_up(RE[1], 1);
        a[6] = __shfl_down(RE[2], 1);
        a[7] = RE[2];
        a[8] = __shfl_up(RE[2], 1);
        // ищем максимум 
        a[2] = __vmaxu4(a[0], __vmaxu4(a[1], a[2]));
        // ищем медиану
        Sort(a[3], a[4]);
        a[4] = __vmaxu4(__vminu4(a[4], a[5]), a[3]);
        // ищем минимум
        a[6] = __vminu4(a[6], __vminu4(a[7], a[8]));
        // ищем финально медиану
        Sort(a[2], a[4]);
        a[4] = __vmaxu4(__vminu4(a[4], a[6]), a[2]);
        // записываем результат, если не вышли за границу картинки
        if (idy + 0 < H - 1 && bound)
            out[(idy) * W + idx] = a[4] & 255;
        if (idy + 1 < H - 1 && bound)
            out[(idy + 1) * W + idx] = (a[4] >> 8) & 255;
        if (idy + 2 < H - 1 && bound)
            out[(idy + 2) * W + idx] = (a[4] >> 16) & 255;
        if (idy + 3 < H - 1 && bound)
            out[(idy + 3) * W + idx] = (a[4] >> 24) & 255;    
    }

    Die zweite Optimierung hängt mit der Tatsache zusammen, dass jeder Thread nicht 4 Punkte gleichzeitig verarbeiten kann, sondern beispielsweise 2 Sätze mit 4 Punkten oder N Sätze mit 4 Punkten. Dies ist erforderlich, um die Belastung des Geräts zu maximieren, da nur für einen 3x3-Filter ein Satz von 4 Punkten nicht ausreicht.


    Die dritte Optimierung, die noch durchgeführt werden kann, besteht darin, die Breite und Höhe des Bildes zu ersetzen. Dafür gibt es mehrere Argumente:


    1. In der Regel haben Bildgrößen feste Werte, z. B. die gängigen Formate Full HD, Ultra HD und 720p. Sie können nur einen Satz vorkompilierter Kernel haben. Diese Optimierung führt zu einer Produktivitätssteigerung von 10-15%.
    2. Mit dem Cuda Toolkit 7.5 ist es möglich, eine dynamische Kompilierung durchzuführen, mit der Sie Werte zur Laufzeit ersetzen können.

    Vollständig optimierter Code sieht dann so aus. Indem Sie die Anzahl der Punkte variieren, können Sie maximale Leistung erzielen. In meinem Fall wurde die maximale Leistung mit numP_char = 3 erreicht, dh drei Sätze mit 4 Punkten oder 12 Punkten pro Thread.


    __global__ BOUNDS(BL_X, 2048 / BL_X)
    void mFilter_3x3(const unsigned char * __restrict in, unsigned char *out)
    {
        int idx = threadIdx.x + blockIdx.x * blockDim.x;
        idx = idx - 2 * (idx / 32);
        int idxW = threadIdx.x % 32;
        int idy = threadIdx.y + blockIdx.y * blockDim.y * (4 * numP_char) + 1;
        unsigned a[numP_char][9];
        unsigned RE[numP_char][3];
        bool bound = (idxW > 0 && idxW < 31 && idx < W - 1);
    #pragma unroll
        for (int z = 0; z < numP_char; ++z)
        {
            RE[z][0] = in[(idy - 1 + z * 4) * W + idx] | (in[(idy - 1 + 1 + z * 4) * W + idx] << 8) | (in[(idy - 1 + 2 + z * 4) * W + idx] << 16) | (in[(idy - 1 + 3 + z * 4) * W + idx] << 24);
            RE[z][1] = in[(idy + z * 4) * W + idx] | (in[(idy + 1 + z * 4) * W + idx] << 8) | (in[(idy + 2 + z * 4) * W + idx] << 16) | (in[(idy + 3 + z * 4) * W + idx] << 24);
            RE[z][2] = in[(idy + 1 + z * 4) * W + idx] | (in[(idy + 1 + 1 + z * 4) * W + idx] << 8) | (in[(idy + 1 + 2 + z * 4) * W + idx] << 16) | (in[(idy + 1 + 3 + z * 4) * W + idx] << 24);
            Sort(RE[z][0], RE[z][1]);
            Sort(RE[z][1], RE[z][2]);
            Sort(RE[z][0], RE[z][1]);
            a[z][0] = __shfl_down(RE[z][0], 1);
            a[z][1] = RE[z][0];
            a[z][2] = __shfl_up(RE[z][0], 1);
            a[z][3] = __shfl_down(RE[z][1], 1);
            a[z][4] = RE[z][1];
            a[z][5] = __shfl_up(RE[z][1], 1);
            a[z][6] = __shfl_down(RE[z][2], 1);
            a[z][7] = RE[z][2];
            a[z][8] = __shfl_up(RE[z][2], 1);
            a[z][2] = __vmaxu4(a[z][0], __vmaxu4(a[z][1], a[z][2]));
            Sort(a[z][3], a[z][4]);
            a[z][4] = __vmaxu4(__vminu4(a[z][4], a[z][5]), a[z][3]);
            a[z][6] = __vminu4(a[z][6], __vminu4(a[z][7], a[z][8]));
            Sort(a[z][2], a[z][4]);
            a[z][4] = __vmaxu4(__vminu4(a[z][4], a[z][6]), a[z][2]);
        }
    #pragma unroll
        for (int z = 0; z < numP_char; ++z)
        {
            if (idy + z * 4 < H - 1 && bound)
                out[(idy + z * 4) * W + idx] = a[z][4] & 255;
            if (idy + 1 + z * 4 < H - 1 && bound)
                out[(idy + 1 + z * 4) * W + idx] = (a[z][4] >> 8) & 255;
            if (idy + 2 + z * 4 < H - 1 && bound)
                out[(idy + 2 + z * 4) * W + idx] = (a[z][4] >> 16) & 255;
            if (idy + 3 + z * 4 < H - 1 && bound)
                out[(idy + 3 + z * 4) * W + idx] = (a[z][4] >> 24) & 255;
        }
    }

    5x5 Median Filter


    Leider konnte ich keine andere Sortierung für das 5x5-Quadrat finden. Das einzige, was Sie sparen können, ist das Herunterladen und Kombinieren von 4 Punkten in einem vorzeichenlosen int. Ich sehe keinen Sinn darin, einen noch längeren Code mitzubringen, da alle Transformationen analog zu einem 3x3-Quadrat durchgeführt werden können.
    In diesem Artikel beschreiben wir einige der Ideen für den kombinierten Betrieb der beiden Quadrate mit der Einführung eines 20 - Elemente. Die von den Autoren vorgeschlagene Methode der vergesslichen Auswahlsortierung führt jedoch immer noch mehr Operationen aus als das unvollständige Metzgersortierungsnetz für 25 Elemente, selbst wenn zwei benachbarte 5 × 5-Quadrate kombiniert werden.


    Leistung


    Abschließend gebe ich Ihnen einige Zahlen, die ich während des Optimierungsprozesses erhalten habe.


    Vorlaufzeit, msBeschleunigung
    Skalare CPUAVX2-CPUGPU Skalar / GPU AVX2 / GPU
    3x3 1920x108022.90,2550,044 (47,1 GP / s)5205.8
    3x3 4096x216097,91,330,172 (51,4 GP / s)5697,73
    5x5 1920x1080134,31,470,260 (7,9 GP / s)5165.6
    5x5 4096x2160569,26,691.000 (8,8 GP / s)5696,69

    Fazit


    Abschließend möchte ich festhalten, dass unter allen gefundenen Artikeln und Implementierungen der Medianfilterung für die GPU der bereits erwähnte Artikel von Interesse ist"Optimierte Hochgeschwindigkeitsimplementierung eines GPU-basierten Medianfilters", veröffentlicht im Jahr 2013. In diesem Artikel schlugen die Autoren einen völlig anderen Sortierungsansatz vor, nämlich die Methode der vergesslichen Auswahl. Das Wesentliche dieser Methode ist, dass wir roundUp (N / 2) + 1 Elemente nehmen, von links nach rechts und zurück gehen und so die minimalen und maximalen Elemente an den Kanten erhalten. Vergessen Sie als Nächstes diese Elemente, fügen Sie am Ende eines der verbleibenden Elemente hinzu und wiederholen Sie den Vorgang. Wenn es bereits nichts hinzuzufügen gibt, erhalten wir ein Array mit 3 Elementen, aus denen der Median leicht ausgewählt werden kann. Einer der Vorteile dieses Ansatzes ist die geringere Anzahl verwendeter Register im Vergleich zu den bekannten Sortierungen.


    Der Artikel besagt, dass die Autoren beim Tesla C2050 ein Ergebnis von 1,8 GPixel / Sek. Erzielt haben. Die Leistung dieser Karte in einfacher Genauigkeit wird auf 1 TFlops geschätzt. Die Leistung des am Test teilnehmenden 780Ti wird auf 5 TFlops geschätzt. Somit ist die spezifische Rechengeschwindigkeit pro 1 TFlops des von mir vorgeschlagenen Algorithmus für ein 3 × 3-Quadrat ungefähr 5,5-mal höher und für ein 5 × 5-Quadrat 2-mal höher als die von den Autoren des Artikels vorgeschlagene. Dieser Vergleich ist nicht ganz richtig, aber näher an der Wahrheit. Auch in diesem Artikel wurde erwähnt, dass zu dieser Zeit ihre Implementierung die schnellste von allen bekannten war.


    Die erzielte Beschleunigung gegenüber der AVX2-Version betrug durchschnittlich das 6-fache. Wenn Sie neue Karten verwenden, die auf der Pascal-Architektur basieren, kann sich diese Beschleunigung um mindestens das Zweifache erhöhen, was ungefähr 100 GPixel / s entspricht.


    Jetzt auch beliebt: