Szybka transpozycja obrazu i optymalizacja filtra Sobela w C (SIMD)

Question

Chcę zaimplementować naprawdę (naprawdę) szybki Sobel operator dla ray-tracera, mojego przyjaciela i napisałem (źródła można znaleźć here). Oto co do tej pory wymyśliłem ...

Po pierwsze, załóżmy, że obraz jest linią obrazów w skali szarości, wiersz po linii, w 8-bitowej tablicy liczb całkowitych bez znaku.

Aby napisać prawdziwy filtr Sobela, muszę obliczyć Gx i Gy dla każdego piksela. Każda z tych liczb jest obliczana za pomocą 6 pikseli obok miejsca pochodzenia. Ale instrukcja SIMD pozwala mi radzić sobie z 16 lub nawet 32 ​​(AVX) pikselami. Mamy nadzieję, że jądro operatora posiada kilka nieruchomości, tak, że można obliczyć Gy przez:

• odjęcie każdego i oraz i + 2 rzędy i zapisać wynik w I + 1 rzędzie innego obrazu (tablicę)
• dodawania i, dwukrotnie z i + 1 oraz i + 2 kolumny dają i + 1 kolumnę końcowego obrazu

zrobiłbym to samo (ale transponowane) obliczyć Gx następnie dodać dwa obrazy.

Kilka uwag:

• Nie dbam o alokacji pamięci, ponieważ wszystko zostanie przydzielone na początku.
• mogę poradzić sobie z przepełnienia i podpisać problemu podzielenie wartości przez cztery (dzięki _mm_srli_epi8) (uint8_t >> 2 - uint8_t >> 2) = int7_t //really store as int8_t
int7_t + uint8_t << 1 >> 2 + int7_t = uint8_t
//some precision is lost but I don't care

Prawdziwy problem mam skierowane jest pójść z wierszy do kolumn. Ponieważ nie mogłem załadować obrazu w rejestrze SIMD w inny sposób. Muszę odwrócić obraz co najmniej trzy razy, prawda?

Po oryginalnym obrazie. Następnie mogę obliczyć pierwszy krok dla Gx i Gy, a następnie odwrócić wynikowe obrazy, aby obliczyć drugi krok.

Tak, tutaj jest moje pytania:

• Czy tego rodzaju realizacji to dobry pomysł?
• Czy istnieje sposób na przetransponowanie tablicy szybciej niż głupi algorytm? (Nie sądzę)
• Gdzie będą wąskie gardła? (jak sądzę?): P)

Answer 1

Myślę, że transpozycja/2-przebieg nie jest dobra do optymalizacji kodu operatora Sobela. Sobel Operator nie jest funkcją obliczeniową, więc marnowanie dostępu do pamięci w przypadku transpozycji/2-pasmowego dostępu nie jest dobre w tym przypadku. Napisałem kilka kodów testowych Sobel Operator, aby zobaczyć, jak szybko może uzyskać SSE. ten kod nie obsługuje pikseli o pierwszej i ostatniej krawędzi, a do obliczenia wartości sqrt() można użyć jednostek FPU.

Obsługa operatora Sobela: 14 mnożenie, 1 pierwiastek kwadratowy, 11 dodawania, 2 min/maks., 12 uprawnień do odczytu i 1 operator dostępu do zapisu. Oznacza to, że możesz przetworzyć komponent w cyklu 20 ~ 30, jeśli dobrze zoptymalizujesz kod.

Funkcja FloatSobel() zajęła 2113044 cykli procesora w celu przetworzenia obrazu 256 * 256 przetwarzania 32,76 cyklu/komponentu. Przekształcę ten przykładowy kod w SSE.

void FPUSobel() 
{ 
    BYTE* image_0 = g_image + g_image_width * 0; 
    BYTE* image_1 = g_image + g_image_width * 1; 
    BYTE* image_2 = g_image + g_image_width * 2; 
    DWORD* screen = g_screen + g_screen_width*1; 

    for(int y=1; y<g_image_height-1; ++y) 
    { 
     for(int x=1; x<g_image_width-1; ++x) 
     { 
      float gx = image_0[x-1] * (+1.0f) + 
         image_0[x+1] * (-1.0f) + 
         image_1[x-1] * (+2.0f) + 
         image_1[x+1] * (-2.0f) + 
         image_2[x-1] * (+1.0f) + 
         image_2[x+1] * (-1.0f); 

      float gy = image_0[x-1] * (+1.0f) + 
         image_0[x+0] * (+2.0f) + 
         image_0[x+1] * (+1.0f) + 
         image_2[x-1] * (-1.0f) + 
         image_2[x+0] * (-2.0f) + 
         image_2[x+1] * (-1.0f); 


      int result = (int)min(255.0f, max(0.0f, sqrtf(gx * gx + gy * gy))); 

      screen[x] = 0x01010101 * result; 
     } 
     image_0 += g_image_width; 
     image_1 += g_image_width; 
     image_2 += g_image_width; 
     screen += g_screen_width; 
    } 
} 

Funkcja SseSobel() zajęła procesorowi procesorowi 613220 przetwarzanie tego samego obrazu 256 * 256. Wymagało to 9,51 cyklu/komponentu i 3,4 razy szybciej niż FPUSobel(). Jest kilka miejsc do optymalizacji, ale nie będzie to więcej niż 4 razy, ponieważ wykorzystano 4-drożny SIMD.

Ta funkcja wykorzystała podejście SoA do przetwarzania 4 pikseli naraz. SoA jest lepsza niż AoS w większości danych tablicowych lub obrazowych, ponieważ musisz przetransponować/przetasować, aby użyć AoS. A SoA znacznie łatwiej zmienia wspólny kod C na kody SSE.

void SseSobel() 
{ 
    BYTE* image_0 = g_image + g_image_width * 0; 
    BYTE* image_1 = g_image + g_image_width * 1; 
    BYTE* image_2 = g_image + g_image_width * 2; 
    DWORD* screen = g_screen + g_screen_width*1; 

    __m128 const_p_one = _mm_set1_ps(+1.0f); 
    __m128 const_p_two = _mm_set1_ps(+2.0f); 
    __m128 const_n_one = _mm_set1_ps(-1.0f); 
    __m128 const_n_two = _mm_set1_ps(-2.0f); 

    for(int y=1; y<g_image_height-1; ++y) 
    { 
     for(int x=1; x<g_image_width-1; x+=4) 
     { 
      // load 16 components. (0~6 will be used) 
      __m128i current_0 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_0+x-1)), _mm_setzero_si128()); 
      __m128i current_1 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_1+x-1)), _mm_setzero_si128()); 
      __m128i current_2 = _mm_unpacklo_epi8(_mm_loadu_si128((__m128i*)(image_2+x-1)), _mm_setzero_si128()); 

      // image_00 = { image_0[x-1], image_0[x+0], image_0[x+1], image_0[x+2] } 
      __m128 image_00 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_0, _mm_setzero_si128())); 
      // image_01 = { image_0[x+0], image_0[x+1], image_0[x+2], image_0[x+3] } 
      __m128 image_01 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_0, 2), _mm_setzero_si128())); 
      // image_02 = { image_0[x+1], image_0[x+2], image_0[x+3], image_0[x+4] } 
      __m128 image_02 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_0, 4), _mm_setzero_si128())); 
      __m128 image_10 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_1, _mm_setzero_si128())); 
      __m128 image_12 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_1, 4), _mm_setzero_si128())); 
      __m128 image_20 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(current_2, _mm_setzero_si128())); 
      __m128 image_21 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_2, 2), _mm_setzero_si128())); 
      __m128 image_22 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_srli_si128(current_2, 4), _mm_setzero_si128())); 

      __m128 gx = _mm_add_ps(_mm_mul_ps(image_00,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_02,const_n_one), 
         _mm_add_ps(_mm_mul_ps(image_10,const_p_two), 
         _mm_add_ps(_mm_mul_ps(image_12,const_n_two), 
         _mm_add_ps(_mm_mul_ps(image_20,const_p_one), 
            _mm_mul_ps(image_22,const_n_one)))))); 

      __m128 gy = _mm_add_ps(_mm_mul_ps(image_00,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_01,const_p_two), 
         _mm_add_ps(_mm_mul_ps(image_02,const_p_one), 
         _mm_add_ps(_mm_mul_ps(image_20,const_n_one), 
         _mm_add_ps(_mm_mul_ps(image_21,const_n_two), 
            _mm_mul_ps(image_22,const_n_one)))))); 

      __m128 result = _mm_min_ps(_mm_set1_ps(255.0f), 
          _mm_max_ps(_mm_set1_ps(0.0f), 
             _mm_sqrt_ps(_mm_add_ps(_mm_mul_ps(gx, gx), _mm_mul_ps(gy,gy))))); 

      __m128i pack_32 = _mm_cvtps_epi32(result); //R32,G32,B32,A32 
      __m128i pack_16 = _mm_packs_epi32(pack_32, pack_32); //R16,G16,B16,A16,R16,G16,B16,A16 
      __m128i pack_8 = _mm_packus_epi16(pack_16, pack_16); //RGBA,RGBA,RGBA,RGBA 
      __m128i unpack_2 = _mm_unpacklo_epi8(pack_8, pack_8); //RRGG,BBAA,RRGG,BBAA 
      __m128i unpack_4 = _mm_unpacklo_epi8(unpack_2, unpack_2); //RRRR,GGGG,BBBB,AAAA 

      _mm_storeu_si128((__m128i*)(screen+x),unpack_4); 
     } 
     image_0 += g_image_width; 
     image_1 += g_image_width; 
     image_2 += g_image_width; 
     screen += g_screen_width; 
    } 
}