Chuyển vị nhanh của hình ảnh và tối ưu hóa Sobel Filter thành C (SIMD)

Question

Tôi muốn triển khai thực hiện (thực sự) nhanh Sobel operator cho một người theo dõi tia một người bạn của tôi và tôi đã viết (nguồn có thể được tìm thấy here). Sau đây là những gì tôi tìm ra cho đến nay ...

Trước tiên, giả sử hình ảnh là một chuỗi cửa hàng hình ảnh có kích thước màu xám theo dòng trong mảng số nguyên không dấu 8 bit.

Để viết bộ lọc Sobel thực, tôi cần tính Gx và Gy cho mỗi pixel. Mỗi số này được tính toán nhờ 6 pixel bên cạnh nguồn gốc. Nhưng lệnh SIMD cho phép tôi xử lý 16 hoặc thậm chí 32 (AVX) pixel. Hy vọng rằng, các hạt nhân của các nhà điều hành có một số tài sản tốt đẹp, vì vậy mà tôi có thể tính toán Gy bởi:

• trừ mỗi i và i + 2 hàng và lưu trữ kết quả trong một i + 1 hàng của một số hình ảnh khác (mảng)
• thêm i, gấp đôi so với i + 1 và i + 2 cột cho i + 1 cột của hình ảnh chính thức

tôi sẽ làm điều tương tự (nhưng hoán) để tính toán Gx sau đó thêm hai hình ảnh.

Một số lưu ý:

• Tôi không quan tâm về cấp phát bộ nhớ từ tất cả mọi thứ sẽ được phân bổ ngay từ đầu.
• tôi có thể đối phó với sự tràn và ký vấn đề phân chia các giá trị bằng bốn (nhờ _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

Vấn đề thực sự tôi phải đối mặt là để đi từ các hàng để các cột. Vì tôi không thể tải ảnh trong thanh ghi SIMD theo cách khác. Tôi phải lật hình ảnh ba lần ít nhất phải không?

Khi ảnh gốc. Sau đó, tôi có thể tính toán bước đầu tiên cho Gx và Gy và sau đó lật các ảnh kết quả để tính bước thứ hai.

Vì vậy, đây là câu hỏi của tôi:

• là loại triển khai một ý tưởng tốt?
• Có cách nào để chuyển đổi mảng nhanh hơn thuật toán câm không? (Tôi không nghĩ vậy)
• Trường hợp bị tắc nghẽn ở đâu? (bất kỳ phỏng đoán nào?: P)

Answer 1

Tôi nghĩ rằng chuyển vị/2-pass không tốt để tối ưu hóa mã Toán tử Sobel. Toán tử Sobel không phải là hàm tính toán, do đó việc lãng phí truy cập bộ nhớ cho truy cập transpose/2-pass không tốt cho trường hợp này. Tôi đã viết một số mã kiểm tra Sobel Operator để xem SSE có thể nhận được nhanh như thế nào. mã này không xử lý các pixel cạnh đầu tiên và cuối cùng, và sử dụng các FPU để tính giá trị sqrt().

Toán tử Sobel cần 14 nhân, 1 căn bậc hai, 11 bổ sung, 2 phút/tối đa, 12 truy cập đọc và 1 toán tử truy cập ghi. Điều này có nghĩa là bạn có thể xử lý một thành phần trong chu kỳ 20 ~ 30 nếu bạn tối ưu hóa mã tốt.

FloatSobel() đã mất 2113044 chu kỳ CPU để xử lý xử lý hình ảnh 256 * 256 32.76 chu kỳ/thành phần. Tôi sẽ chuyển đổi mã mẫu này thành 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; 
    } 
} 

Hàm SseSobel() mất 613220 Chu kỳ CPU để xử lý cùng một hình ảnh 256 * 256. Nó mất 9.51 chu kỳ/thành phần và 3.4 thời gian nhanh hơn FPUSobel(). Có một số không gian để tối ưu hóa nhưng nó sẽ không nhanh hơn 4 lần vì nó sử dụng SIMD 4 chiều.

Hàm này sử dụng phương pháp SoA để xử lý 4 pixel cùng một lúc. SoA là tốt hơn so với AoS trong hầu hết các mảng hoặc dữ liệu hình ảnh bởi vì bạn phải transpose/shuffle để sử dụng AoS. Và SoA dễ dàng thay đổi mã C phổ biến thành mã 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; 
    } 
}