Wstęp

Rozmycie obrazu (w najprostszej wersji) realizuje się jako sumę ośmiu sąsiednich pikseli 'o' piksela 'x', lub żeby być precyzyjnym: blur(x) = (a+b+c + d+x+e + f+g+h)/9.

o o o   a b c
o x o   d x e
o o o   f g h

Rozmycie przeprowadza się w podwójnej pętli:

for y := y_min to y_max do
        for x := x_min to x_max do
                begin
                        a = getpixel(x-1, y-1);
                        b = getpixel(x  , y-1);
                        c = getpixel(x+1, y-1);

                        d = getpixel(x-1, y);
                        x = getpixel(x  , y);
                        e = getpixel(x+1, y);

                        f = getpixel(x-1, y+1);
                        g = getpixel(x  , y+1);
                        h = getpixel(x+1, y+1);

                        avg = (a+b+c + d+x+e + f+g+h)/9;

                        putpixel(x, y, avg);  # zapis do *nowego* obrazu
                end;

Jak łatwo policzyć na każdy rozmywany piksel należy pobrać (getpixel) aż 9 pikseli, co dla całego obrazu daje 9 ⋅ (y_max − y_min) ⋅ (x_max − x_min). Np. dla obrazu o rozdzielczości 320x200 daje to aż 576000 odwołań do obrazu.

Ponadto należy, dla każdego piksela, wykonać 8 operacji dodawania.

W artykule pokażę jak te liczby zmniejszyć — i w efekcie przyspieszyć algorytm — oraz jak użyć rozkazów MMX/SSE.

Uwaga 1: Ze względu na wydajność należy rozmywać obraz obcięty z każdej strony po jednym pikselu, tj. np. zamiast (0,0)-(320,20) rozmywać (1,1)-(319,199). Problem sprawiają bowiem piksele brzegowe, np. w punkcie (0,0) musielibyśmy pobrać 5 nieistniejących w obrazie pikseli. Pół biedy z przyjęciem ich wartości (czarne, białe, powtórzone, losowe, itp.), największe zmartwienie to właśnie sprawdzenie czy piksel „istnieje” czy nie.

Dopisek 14.07.2007: jednak w przedstawionej tutaj modyfikacji, dość łatwo poradzić sobie z tymi „brakującymi” pikselami. W przykładowym programie zostało to rozwiązane, można ustalić kolor obrzeża.

Uwaga 2: W przykładowym kodzie wykorzystuję funkcje getpixel/putpixel, które stawiają, bądź pobierają piksel na zadanym, domyślnym obrazie. Faktycznie te operacje wykonuje się na wskaźnikach, ale dla czytelności użyłem postaci funkcyjnej.

Optymalizacja

Rozpatrzmy rozmycie dwóch kolejnych pikseli w jednej linii, tzn. najpierw rozmywamy piksel (x,y) ozn. X, a następnie (x+1,y) ozn. Y.

a b c d
e X Y f
g h i j

blur(X) = (a+b+c + e+X+Y + g+h+i)/9
blur(Y) = (b+c+d + X+Y+f + h+i+j)/9

Zobaczmy, że w wyrażeniu na blur(Y) część pikseli powtarza się, doszły zaledwie trzy nowe: d, f oraz j. Teraz jeszcze jedna rzecz: celowo grupuję w wyrażeniach po 3 piksele, bowiem każda „grupa” (suma cząstkowa) pochodzi z innej linii. Proszę teraz zobaczyć jak wyglądałby wzór na rozmycie piksela h, znajdującego się pod pikselem X, tj. (x,y+1):

blur(h) = (e+X+Y + g+h+i + .+.+.)/9

Aż dwie „grupy” występujące w wyrażeniu na blur(X) powtarzają się w blur(h). Zatem wystarczy zapamiętać wszystkie sumy cząstkowe z dwóch wcześniejszych linii y, co oszczędzi 4 dodawania i 6 odczytów pikseli. Do efektywnego liczenia sum poszczególnych grup wykorzystamy odnotowany wcześniej fakt, że przejście do następnego piksela w jednej linii y powoduje pojawienie się 3 nowych pikseli w wyrażeniu: po jednym na każdą grupę.

Na początku kod w C który dla zadanej linii y liczy wszystkie sumy w obrazie grayscale. Ponieważ na sumę cząstkową składają się wartości trzech pikseli (bajtów), więc musi być ona przechowywana w słowie.

#define SIZE_X  ...     // rozmiary
#define SIZE_Y  ...     // obrazu XxY

typedef byte            Image_line[SIZE_X];     // linia obrazu
typedef Image_line      Image[SIZE_Y];          // obraz
typedef word            Sum[SIZE_X-2];          // tablica sum

/*
 *      image   - wskaźnik do początku y-owej linii
 *      count   - szerokość obrazu
 *      sum     - tablica sum (rozmiar count-2)
 */
void calc_sum(Image_line L, Sum S) {
        byte a, b, c;
        int  x, i;

        a = L[0];       // pobranie dwóch początkowych
        b = L[1];       // pikseli
        i = 0;
        for (x=2; x<count; x++) {
                c = L[x];               // w każdej pętli pobierany
                sum[i++] = a+b+c;       // jest jeden nowy piksel

                a = b;                  // natomiast dwa pozostają
                b = c;                  // po wcześniejszej iteracji
        }
}

Funkcja rozmywająca obraz.

/*
 *      image_in        - obraz rozmywany
 *      image_out       - wynik rozmycia
 */
void blur_image(Image image_in, Image image_out) {
        Sum  a,  b,  c;
        Sum *A, *B, *C, *T;
        word s;
        int x, y;

        A = &a; // używane są wskaźniki do tablic sum
        B = &b; // żeby uniknąć przesłań pamięci;
        C = &c; // wystarczy "rotacja" wskaźników

        calc_sum(image_in[0], A);
        calc_sum(image_on[1], B);
        for (y=1; y < SIZE_Y-1; y++) {

                // w każdej pętli liczona są wyłącznie sumy cząstkowe
                // dla linii y+1; sumy cząstkowe z linii y i y-1
                //pochodzą z wcześniejszej iteracji
                calc_sum(image_in[y+1], C);

                // rozmycie
                for (x=0; x < SIZE_X-2; x++) {
                        s = (A[x] + B[x] + C[x])/9;
                        putpixel(image_out, x, y, s);
                }

                // "rotacja" wskaźników
                T = A;
                A = B;
                B = C;
                C = T;
        }
}

W implementacji asemblerowej powyższej funkcji problemem może być przechowywanie trzech wskaźników, będą bowiem potrzebne 3 rejestry. Można połączyć tablice A, B i C w jedną i operować na przesunięciu względem jej początku.

void blur_image(Image image_in, Image image_out) {
#define SUM_SIZE (SIZE_X - 2)
        word abc[SUM_SIZE*3];
        int  offset = 0;
        word s;
        int x, y;

        calc_sum(image_in[0], (Sum*)abc[0*SUM_SIZE]);
        calc_sum(image_in[1], (Sum*)abc[1*SUM_SIZE]);
        for (y=1; y<SIZE_Y-1; y++) {

                calc_sum(image[y+1], (Sum*)t[offset]);

                // rozmycie
                for (x=0; x < SIZE_X-2; x++) {
                        s = (abc[0*SUM_SIZE] +
                             abc[1*SUM_SIZE] +
                             abc[2*SUM_SIZE])/9;
                        putpixel(image_out, x, y, s);
                }

                // rotacja bufora
                if (offset == 2*SUM_SIZE)
                        offset = 0;
                else
                        offset += SUM_SIZE;
        }
#undef SUM_SIZE
}

calc_sum:
; edi - tablica 'ImageLine'
; esi - tablica 'Sum'

mov ecx, SIZE_X-2 ; liczba iteracji

xor ah, ah
xor bh, bh
xor dh, dh

mov al, [edi  ] ; pobierz a (=[edi-2])
mov bl, [edi+1] ; pobierz b (=[edi-1])
                ;         c (=[edi])

add edi, 2

_sum:
        mov dl, [edi] ; pobierz c, tylko jedną nowa wartość

        add ax, bx    ; a jest tracone w następnej iteracji
        add ax, dx    ; więc zostanie użyte do obliczeń

        mov [esi], ax ; zapisz sumę do tablicy 'Sum'

        mov ax, bx    ; a=b
        mov bx, dx    ; b=c

        inc edi
        add esi, byte 2

        loop _sum
ret

calc_sum:
mov ecx, SIZE_X/8

pxor mm7, mm7

_sum:
        ; edi -> |h|g|f|e|d|c|b|a|...
        movq mm0, [edi+0]  ; mm0 = |h|g|f|e|d|c|b|a|
        movq mm1, [edi+1]  ; mm1 = |i|h|g|f|e|d|c|b|
        movq mm2, [edi+2]  ; mm2 = |j|i|h|g|f|e|d|c|

        movq mm3, mm0
        movq mm4, mm1
        movq mm5, mm2

        punpcklbw mm0, mm7 ; mm0 = | d | c | b | a |
        punpcklbw mm1, mm7 ; mm1 = | e | d | c | b |
        punpcklbw mm2, mm7 ; mm2 = | f | e | d | c |

        punpckhbw mm3, mm7 ; mm3 = | h | g | f | e |
        punpckhbw mm4, mm7 ; mm3 = | i | h | g | f |
        punpckhbw mm5, mm7 ; mm3 = | j | i | h | g |

        paddw mm0, mm1     ;
        paddw mm0, mm2     ; mm0 = |def|cde|bcd|abc|

        paddw mm3, mm4     ;
        paddw mm3, mm5     ; mm3 = |hij|ghi|fgh|efg|

        movq [esi+0], mm0
        movq [esi+8], mm1

        add edi, 8
        add esi, 16
        loop _sum
ret

Dzielenie

Nieporuszanym jeszcze problemem jest dzielenie przez 9. W przypadku kodu x86 (a także kodu w C) w zupełności wystarczy tablica z prekalkulowanymi wartościami.

byte div_tbl[9*255];
...
for (i=0; i<9*255; i++)
        div_tbl[i] = i/9;

...
void blur(Image image)
...
                s = div_tbl[A[x] + B[x] + C[x]];
...

Natomiast MMX dostarcza tylko rozkazów mnożenia, wykorzystamy jednak pewien trik. Jak być może wiesz, przesunięcie bitowe w prawo, to dzielenie przez pewną potęgę dwójki. Można tę własność wykorzystać — wzór na dzielenie całkowitoliczbowe przez dowolną liczbę k jest następujący:

x/k = (x*(2^n/k)) >> n

Czynnik 2^n/k jest liczony w trakcie kompilacji; w rozważanym przypadku (k=9, n=16) wyniesie on d = 65536/9 = 7281 = 0x1c71. Dlaczego n=16? Rozkazy MMX pmulhw oraz pmuluhw obliczają starsze słowo z mnożenia, a więc faktycznie, przesunięty w prawo o 16 bitów, (32-bitowy) wynik mnożenia.

Dopisek 14.07.2007: żeby naprawdę mnożenie, a potem przesunięcie dało identyczny wynik co dzielenie, należy przed przesunięciem dodać do wyniku mnożenia wartość d, czyli (x * d + d) >> 16. Bez dodawania różnica wartości pomiędzy (x * d) >> 16 a zwykłym dzieleniem całkowitoliczbowym jest co najwyżej 1. Dlatego wolałem zostać przy tym prostszym wyrażeniu, bo z jednej strony różnica jest niezauważalna, z drugiej odpada jedno działanie.

Co prawda rozkaz pmulhw mnoży liczby ze znakiem, lecz 0x1c71 w zapisie U2 jest liczbą dodatnią, więc nie będzie problemu z jego wykorzystaniem.

segment .data

mmx_coef dd 0x1c711c71, 0x1c711c71

segment .code

calc_average_mmx:
; edi - A
; esi - B
; ebx - C
; edx - linia obrazu do której zapisywany jest wynik

movq mm2, [mmx_coef]
mov  ecx, SIZE_X-2

_sum:
        movq mm0, [edi]
        addw mm0, [esi]
        addw mm0, [ebx]   ; sumuj 4 młodsze

        movq mm1, [edi+8]
        addw mm1, [esi+8]
        addw mm1, [ebx+8] ; i 4 starsze grupy pikseli

        pmulhw mm0, mm2   ; mm0/9, mm1/9
        pmulhw mm1, mm2   ; packed_word będą w zakresie 0..255

        packuswb mm0, mm1 ; połącz wyniki

        movq [edx], mm0   ; zapisz

        add edi, 16 ; przesuń wskaźniki
        add esi, 16 ;
        add ebx, 16 ;
        add edx, 8  ;
        loop _sum
ret

Przy założeniu, że przedstawiona powyżej funkcja calc_average_mmx, została zadeklarowana w C:

calc_average_mmx(Image_line, Sum, Sum, Sum);

można ją wykorzystać w sposób następujący:

void blur_image(Image image_in, Image image_out) {
        Sum  a,  b,  c;
        Sum *A, *B, *C, *T;
        word s;
        int x, y;

        A = &a;
        B = &b;
        C = &c;

        calc_sum(image_in[0], A);
        calc_sum(image_out[1], B);
        for (y=1; y<SIZE_Y-1; y++) {

                calc_sum(image[y+1], C);
                calc_average_mmx(image_out[y], A, B, C);

                T = A;
                A = B;
                B = C;
                C = T;
        }
}

Implementacja SSE2 oraz MMX2 [31.05.2008]

W przykładowym programie zamieściłem również implementację SSE2. Poza wyrównaniem wskaźników i użyciu rozkazu PALIGNR do pobierania niewyrównanych danych w calc_sum cały kod zasadniczo został przełożony z wersji MMX automatycznie.

Natomiast w procedurze nazwanej umownie MMX2 zamiast odwołań do niewyrównanych adresów wyciąga się odpowiednie dane za pomocą dwóch przesunięć bitowych oraz jednej opracji sumy bitowe. Zamiast:

0:
        movq 0(%esi), %mm0
        movq 1(%esi), %mm1
        movq 2(%esi), %mm2

        // ...
        // dalsza część pętli [wykorzystywane dodatkowo mm3, mm4]
        // ...

        jnz 0b

Można czytać w jednej iteracji wyrównane do granicy 8 bajtów dane (a przetwarzać 16 kolejnych bajtów) i tworzyć podzakresy sekwencją dwóch przesunięć bitowych i operacji sumy:

;
        movq 0(%esi), %mm6
0:
        // mm0:mm3 -- 16 kolejnych bajtów
        movq %mm6,    %mm0
        movq 8(%esi), %mm3

        // zapamiętanie mm3 dla następnej iteracji
        movq %mm3,    %mm6

        movq  %mm0,   %mm1
        movq  %mm0,   %mm2
        movq  %mm3,   %mm4

        // utworzenie mm1
        pslrq $8,     %mm1
        psllq $64-8,  %mm3
        por   %mm3,   %mm1

        // utworzenie mm2
        pslrq $16,    %mm2
        psllq $64-16, %mm4
        por   %mm4,   %mm2

        // ...
        // dalsza część pętli [wykorzystywane dodatkowo mm3, mm4]
        // ...

        jnz 0b

Na Pentium II, pomimo sporego rozrostu kodu, to rozwiązanie pozwala przyspieszyć kod o ok. 30 procent; przyspieszenie w stosunku do analogicznego kodu x86 wzrasta z ok. 3.8 do ok. 4.7. Na Celeronie M praktycznie bez zmian, chyba nawet wolniej działała druga wersja.

Podsumowanie

Celem opisanych wyżej algorytmów, jest minimalizacja odwołań do pamięci i liczby dodawań. W przypadku bezpośredniej implementacji wzoru liczba odwołań do pamięci obrazu wynosiła 9*(SIZE_X-2)*(SIZE_Y-2), natomiast liczba operacji dodawania 8*(SIZE_X-2)*(SIZE_Y-2).

Po optymalizacji liczba odwołań do pamięci obrazu zmaleje, wyniesie tylko SIZE_X*SIZE_X. Co prawda liczba odwołań do tablic sum a, b i c (lub abc) nie będzie mniejsza, ale z racji ich relatywnie małego rozmiaru, będzie mniej nietrafień w cache.

Liczba operacji sumowania wynosie natomiast 2*(SIZE_X-2)*(SIZE_Y) + 2*(SIZE_Y)*(SIZE_X-2).

Porównanie metody wejściowej i zoptymalizowanego algorytmu dla obrazu o wielkości 320x200 pikseli

Przy użyciu kodów MMX, zaprezentowanych w tym artykule, przyrost prędkości był naprawdę imponujący — lekko licząc ok. 2-3 razy w stosunku do kodu x86 używającego tej samej, zoptymalizowanej metody.

Poniżej wyniki testów przykładowego programu przeprowadzane na dwóch maszynach działających pod kontrolą Linuksa.

Przykład

• W repozytorium jest przykładowy program zawierający pięć procedur realizujących ulepszoną metodę:
  • procedura zakodowana w C
  • wersja asemblerowa (x86)
  • wersja MMX
  • druga wersja MMX
  • SSE2