Beispiel: Sobel-Operator

Die Bilderkennung identifiziert häufig im ersten Schritt die Umrisse von möglichen Objekten. Der Sobel-Operator [2] detektiert Kanten und findet sich auch in gängigen Bildverarbeitungsprogrammen wie Gimp. Für Grauwert-Bilder funktioniert der Sobel-Operator so: Die Farbe Schwarz bekommt den Grauwert 0, Weiß wird mit dem Wert 255 kodiert. Grautöne liegen dazwischen. Eine Kante ist eine abrupte Änderung des Grauwerts, die es nun zu detektieren gilt.

Der Sobel-Operator sieht sich für jeden Bildpunkt die Nachbarpixel an und berechnet Differenzen in horizontaler und vertikaler Richtung. Die Differenzbildung hebt Unterschiede hervor (Differenz von Grauwerten ist groß) und unterdrückt Flächen ähnlichen oder gleichen Grauwerts (Differenz ist klein oder null). Das Ergebnis ist oft sehr dunkel, wobei die markanten Umrisse hervorstechen.

Das Berechnungsgitter baut sich aus einer zweidimensionalen Struktur von Blöcken auf. Jeder Thread eines Blocks ermittelt den Sobel-Wert für genau einen Bildpunkt aus einer Gewichtung der Grauwerte der Nachbarpunkte. Der Operator berechnet dazu Differenzen in horizontaler und vertikaler Richtung getrennt.

Abbildung 4 zeigt das Prinzip für die horizontale Richtung. Die Zahlen der Nachbarpunkte geben die Gewichtung an, wie diese in der Differenzbildung berücksichtigt werden:

Abbildung 4: Der Sobel-Wert für einen Bildpunkt berechnet sich aus den Grauwert-Differenzen der Nachbarpunkte. Neben dem dargestellten horizontalen benötigt die Berechnung des Gradienten auch den vertikalen Sobel-Operator.

sh(x,y)=1*g(x-1,y+1)-1*g(x-1,y-1)+2*g(x,y+1) -2*g(x,y-1)+1*g(x+1,y+1)-1*g(x+1,y-1)

In analoger Weise wird sv als Differenz in vertikaler Richtung berechnet. Der zweite Schritt ermittelt den Betrag des Gradienten:

s(x,y)=sqrt(sh*sh+sv*sv)

Die Grundstruktur eines Cuda-Programms ist in den meisten Fällen ähnlich zu der »main()«-Funktion in Listing 2. Die »main()«-Funktion läuft auf dem Host (CPU des PC) und kümmert sich um das Speichermanagement zum Device (GPU der Grafikkarte). Die Funktion »readpgm()« liest die Bild-Pixel und die -Dimensionen aus einer Datei.

01 int main() {
02 
03 cout << "Beispielprogramm fuer Sobel-Operator mit CUDA-Umgebung" << endl;
04 
05 // --Grauwertbild einlesen und Speicher auf dem Host allokieren
06 unsigned char* ipicture=0;
07 int length=0,width=0;
08 readpgm("pic/starbucks.pgm",ipicture,length,width);
09 
10 // --Speicher auf dem Device (GPU) allokieren
11 int memsize=length*width;
12 unsigned char* device_ipicture=0,*device_opicture=0;
13 // --Speicher fuer Ein- sowie Ausgabebild
14 cudaMalloc((void**)&device_ipicture,memsize);
15 cudaMalloc((void**)&device_opicture,memsize);
16 
17 // --Speicher vom Host auf das Device kopieren
18 cudaMemcpy(device_ipicture,ipicture,memsize,cudaMemcpyHostToDevice);
19 
20 // --Kernelfunktion starten
21 sobel(device_ipicture,device_opicture,length,width);
22 
23 // --Speicher vom Device auf den Host kopieren
24 unsigned char* opicture=new unsigned char[memsize];
25 cudaMemcpy(opicture,device_opicture,memsize,cudaMemcpyDeviceToHost);
26 
27 // --Speicher auf dem Device freigeben
28 cudaFree(device_ipicture);
29 cudaFree(device_opicture);
30 
31 savepgm("result.pgm",opicture,length,width);
32 
33 // --Speicher auf dem Host freigeben
34 delete ipicture;
35 delete opicture;
36 
37 return 0;
38 }

Hin und her kopieren

CPU und GPU verwenden unterschiedlichen Speicher. Analog zu den Standard-C-Funktionen gibt es daher »cudaMalloc()« und »cudaFree()«, um Speicher auf der Grafikkarte zu allokieren und freizugeben. Das Beispiel legt sowohl auf dem Host als auch auf dem Device jeweils zwei Speicherbereiche an, für das Eingangs- und das Sobel-Bild.

Der nächste Schritt kopiert den Speicherinhalt des Eingangsbilds (Abbildung 5) vom Host auf das Device. Hierzu stellt Cuda die Funktion »cudaMemcpy()« bereit. Diese arbeitet ähnlich dem Standard-C-Pendant, mit einer Ausnahme: Der zusätzliche Parameter gibt an, ob sie vom Host zum Device oder in die andere Richtung kopiert.

Abbildung 5: Ein Anwendungsfall für den Sobel-Operator und seine parallelisierte Umsetzung: Gesucht sind die markanten Kanten an einer New Yorker Straßenecke. Diese Bilddatei dient als Eingabe.

Die »sobel()«-Funktion (Listing 3) dimensioniert das Gitter. Sie berücksichtigt dabei, dass sich auf dem Bildrand wegen der fehlenden Nachbarpunkte kein Sobel-Wert berechnen lässt. Eine Faustregel besagt, dass jeder Block mindestens 32 Threads hat und dass es mindestens genauso viele Blöcke wie Prozessoren gibt. Ansonsten würde die zur Verfügung stehende Rechenleistung der GPU nicht effizient genutzt, da einige Prozessoren sich langweilen. Danach startet der Kernel. Die Funktion wird neben den Funktionsparametern mit den Gitter- und Block-Dimensionen in dreifach spitzen Klammern aufgerufen.

01 // --Hier werden die Gitter- und Blockdimensionen definiert und der Kernel gestartet.
02 
03 const int numThreadsPerDim = 16;
04 
05 extern "C" void sobel(unsigned char* in, unsigned char* out, int length, int width) {
06 
07 // --Dimension eines Blocks ist quadratisch (bspw. 16x16 = 256)
08 dim3 dimBlock(numThreadsPerDim,numThreadsPerDim);
09 
10 // --Berechnung der Anzahl der Bloecke (Dimension des Gitters). Das Bildformat wird
11 // --um 2 Punkte jeweils reduziert, da auf dem Bildrand nicht gerechnet wird.
12 dim3 dimGrid((length-2+dimBlock.x-1)/dimBlock.x,(width-2+dimBlock.y-1)/dimBlock.y);
13 
14 // --Starte die Kernelfunktion auf dem Gitter
15 sobel_operator<<<dimGrid,dimBlock>>>(in,out,length,width);
16 
17 }

Haben alle Threads ihre Arbeit zu Ende ausgeführt, kopiert der letzte Schritt mit »cudaMemcpy()« das Ergebnis vom Device zurück auf den Host. Die Funktion »savepgm()« speichert das Ergebnisbild in einer Datei. Im Gegensatz zu den schwergewichtigen Threads in Multicore-CPU-Systemen besteht der Kernel in vielen Fällen nur aus wenigen Zeilen Code, meist unter 100. Die offene Frage bleibt, woher der Kernel weiß, auf welchen Daten er arbeiten soll, wenn alle Threads mit den gleichen Parametern starten. Daher ist eine Kernelfunktion (Listing 4) so aufgebaut, dass sie anhand von »blockDim«, »blockIdx« und »threadIdx« erfahren kann, wo sie gerade läuft. Danach erfolgt die eigentliche Bearbeitung. In Abbildung 6 ist das Ergebnis nach Anwendung des Sobel-Operators dargestellt, das ursprünglich sehr dunkle Bild ist invertiert.

Abbildung 6: Das Bild nach Anwendung des Sobel-Operators. Flächen gleichen Grauwerts werden unterdrückt, Kanten werden hingegen hervorgehoben. Zur besseren Darstellung ist das Bild invertiert.

01 // --Kernel-Funktion fuer Sobel-Operator
02 
03 #define AT(x,y)         ((y)*length+(x))
04 
05 __global__ void sobel_operator(unsigned char* in, unsigned char* out, int length, int width) {
06 
07 // --Wer bin ich? Bestimmung des Bildpunktes
08 // --Es wird jeweils um 1 inkrementiert, da auf dem Rand kein Sobel-Wert berechnet wird.
09 int x=blockDim.x*blockIdx.x+threadIdx.x+1;
10 int y=blockDim.y*blockIdx.y+threadIdx.y+1;
11 
12 // --Anwendung des Operators (Achte auf Rand!)
13 if (x<length-2 && y<width-2) {
14 
15     // --Sobel-Operator in horizontaler Richtung
16     int sh=1*in[AT(x-1,y+1)]-1*in[AT(x-1,y-1)]+
17            2*in[AT(x,y+1)]-2*in[AT(x,y-1)]+
18            1*in[AT(x+1,y+1)]-1*in[AT(x+1,y-1)];
19 
20     // --Sobel-Operator in vertikaler Richtung
21     int sv=1*in[AT(x+1,y+1)]-1*in[AT(x-1,y+1)]+
22            2*in[AT(x+1,y)]-2*in[AT(x-1,y)]+
23            1*in[AT(x+1,y-1)]-1*in[AT(x-1,y-1)];
24 
25     //  --Zusaetzlich invertieren zur besseren Darstellung
26     out[AT(x,y)]=255-(unsigned char)sqrt(float(sh)*float(sh)+float(sv)*float(sv));
27     }
28 }

Nicht alle Probleme lassen sich so offensichtlich parallelisieren und lösen wie bei der Bildbearbeitung mit dem Sobel-Operator. Die GPU ist mit dem Plan entstanden, ein zwei- oder dreidimensionales Bild in Bereiche zu zerlegen und den Grafikprozessor dann jeden Bereich, jedes Pixel fast unabhängig von allen anderen berechnen zu lassen.

Die Welt ist aber manchmal komplexer. Berechnungen, die einen Wert ermitteln, der vom Vorgänger abhängig ist, lassen sich sequenziell mit einer For-Schleife schnell lösen. Die Rückführung eines Vektors von Werten auf einen einzelnen wird als Reduktion bezeichnet. Schon eine Operation wie die Bildung der Summe einer Zahlenfolge ist auf den ersten Blick schwer parallel auszuführen. Aber hier hilft die Mathematik weiter: Es gibt einige Operationen, bei denen es egal ist, in welcher Reihenfolge sie ablaufen. Sind etwa die Zahlen 12, 4 und 23 zu addieren, spielt es keine Rolle, ob man erst 12 und 4 addiert und dann anschließend 23 oder ob man zuerst 4 und 23 und dann 12 addiert. Diese nützliche Eigenschaft wird als Assoziativität bezeichnet.