3
\$\begingroup\$

My program makes an intensive use of small sub-images extracted using bilinear interpolation from larger grayscale images (some call this problem subpixel translation).

I am using the following C++ function for this purpose:

template<typename T>
bool extract_patch_bilin(const cv::Point_<T>& patch_ctr, const cv::Mat_<uchar> &img, cv::Mat_<uchar> &patch)
{
    // ...
    // Precondition checks
    // ...

    const int hsize = patch.rows/2;
    int floorx=(int)floor(patch_ctr.x)-hsize, floory=(int)floor(patch_ctr.y)-hsize;
    // Check that the full extracted patch is inside the image
    if(floorx<0 || img.cols-1<floorx+patch.cols || floory<0 || img.rows-1<floory+patch.rows)
        return false;

    // Compute the constant bilinear weights
    T x=patch_ctr.x-hsize-floorx;
    T y=patch_ctr.y-hsize-floory;
    T xy = x*y;
    T w00=1-x-y+xy, w01=x-xy, w10=y-xy, w11=xy;
    // Prepare image resampling loop
    int img_stride = img.cols-patch.cols;
    uchar* buff_img0 = (uchar*)img.data+img.cols*floory+floorx;
    uchar* buff_img1 = buff_img0+img.cols;
    uchar* buff_patch = (uchar*)patch.data;
    for(int v=0; v<patch.rows; ++v,buff_img0+=img_stride,buff_img1+=img_stride) {
        for(int u=0; u<patch.cols; ++u,++buff_patch,++buff_img0,++buff_img1) {
            // Interpolate pixel value and clamp to [0,255] range
            buff_patch[0] = cv::saturate_cast<uchar>(buff_img0[0]*w00+buff_img0[1]*w01+buff_img1[0]*w10+buff_img1[1]*w11);
        }
    }
    return true;
}

Since this function currently is the bottleneck of my program (it is called several millions times) and since parallelisation is already used for other parts of my program, I tried to optimize it using SSE. This is the implementation I wrote to extract 8x8 patches:

template<typename T>
bool extract_patch_bilin_8x8(const cv::Point_<T>& patch_ctr, const cv::Mat_<uchar> &img, cv::Mat_<uchar> &patch)
{
    // ...
    // Precondition checks
    // ...

    const int hsize = patch.rows/2;
    int floorx=(int)floor(patch_ctr.x)-hsize, floory=(int)floor(patch_ctr.y)-hsize;
    // Check that the full extracted patch is inside the image
    if(floorx<0 || img.cols-1<floorx+patch.cols || floory<0 || img.rows-1<floory+patch.rows)
        return false;
    // Also check that we are not too close to the end of the image buffer, to avoid illegal access when loading pixel values
    if(floory+patch.rows==img.rows-1 && img.cols-9<floorx+patch.cols)
        return false;

    // Compute the constant bilinear weights
    T x=patch_ctr.x-hsize-floorx;
    T y=patch_ctr.y-hsize-floory;
    T xy = x*y;
    T w00=1-x-y+xy, w01=x-xy, w10=y-xy, w11=xy;
    // Prepare image resampling loop
    uchar* buff_img0 = (uchar*)img.data+img.cols*floory+floorx;
    uchar* buff_img1 = buff_img0+img.cols;
    uchar* buff_patch = (uchar*)patch.data;
    // Precompute weighting SIMD variables
    const __m128i CONST_0 = _mm_setzero_si128();
    __m128i w00x256_32i = _mm_set1_epi32(cvRound(w00*256));
    __m128i w01x256_32i = _mm_set1_epi32(cvRound(w01*256));
    __m128i w10x256_32i = _mm_set1_epi32(cvRound(w10*256));
    __m128i w11x256_32i = _mm_set1_epi32(cvRound(w11*256));
    __m128i w00x256_16i = _mm_packs_epi32(w00x256_32i,w00x256_32i);
    __m128i w01x256_16i = _mm_packs_epi32(w01x256_32i,w01x256_32i);
    __m128i w10x256_16i = _mm_packs_epi32(w10x256_32i,w10x256_32i);
    __m128i w11x256_16i = _mm_packs_epi32(w11x256_32i,w11x256_32i);
    // Process pixels
    for(int v=0; v<patch.rows; ++v,buff_img0+=img.cols,buff_img1+=img.cols,buff_patch+=8) {
        ////////////////////////////////
        // Load the data (16 pixels in one load)
        ////////////////////////////////
        __m128i val00 = _mm_loadu_si128((__m128i*)buff_img0);
        __m128i val01 = _mm_loadu_si128((__m128i*)(buff_img0+1));
        __m128i val10 = _mm_loadu_si128((__m128i*)buff_img1);
        __m128i val11 = _mm_loadu_si128((__m128i*)(buff_img1+1));
        ////////////////////////////////
        // Process only the lower 8 values
        ////////////////////////////////
        // Unpack into 16-bits integers
        __m128i val00_lo = _mm_unpacklo_epi8(val00,CONST_0);
        __m128i val01_lo = _mm_unpacklo_epi8(val01,CONST_0);
        __m128i val10_lo = _mm_unpacklo_epi8(val10,CONST_0);
        __m128i val11_lo = _mm_unpacklo_epi8(val11,CONST_0);
        // Multiply with the integer weights
        __m128i w256val00_lo = _mm_mullo_epi16(val00_lo,w00x256_16i);
        __m128i w256val01_lo = _mm_mullo_epi16(val01_lo,w01x256_16i);
        __m128i w256val10_lo = _mm_mullo_epi16(val10_lo,w10x256_16i);
        __m128i w256val11_lo = _mm_mullo_epi16(val11_lo,w11x256_16i);
        // Divide by 256 to get the result of the multiplication with floating-point weights
        __m128i wval00_lo = _mm_srli_epi16(w256val00_lo,8);
        __m128i wval01_lo = _mm_srli_epi16(w256val01_lo,8);
        __m128i wval10_lo = _mm_srli_epi16(w256val10_lo,8);
        __m128i wval11_lo = _mm_srli_epi16(w256val11_lo,8);
        // Add pairwise
        __m128i sum0_lo = _mm_add_epi16(wval00_lo,wval01_lo);
        __m128i sum1_lo = _mm_add_epi16(wval10_lo,wval11_lo);
        __m128i final_lo = _mm_add_epi16(sum0_lo,sum1_lo);
        ////////////////////////////////
        // Repack all values
        ////////////////////////////////
        __m128i final_val = _mm_packus_epi16(final_lo,CONST_0);
        if(v<patch.rows-1)
            _mm_storeu_si128((__m128i*)buff_patch,final_val);
        else {
            // Careful for the last row, we have to process the 8 pixels manually to avoid illegal access outside the patch image buffer
            for(int u=0; u<patch.rows; ++u)
                buff_patch[u] = final_val.m128i_u8[u];
        }
    }
}

This implementation only compute an approximation of the requested patch, since it uses integer operations instead of floating-point ones, but this is enough for my purposes.

My real problem is that there is no significant speed-up using this SSE implementation. I also tried a similar implementation for 16x16 patches, in order to avoid wasting half the computational power, but in this case I have to process twice as many rows so I'm not sure that this is really a better idea. And even with this 16x16 implementation, there is no speedup either.

Any idea on how I could accelerate my initial function using SSE (or any other way)?

Here is a Visual Studio (2010) solution to quickly set up a benchmark environment.

I am getting a speed-up of 3.5 with the SSE function over the regular function, but I was expecting a speed-up closer to 8. Would 3.5 be considered a reasonable speed-up for SSE optimization? Is it possible to improve this?

\$\endgroup\$
  • \$\begingroup\$ What's your target instruction set? SSE2 is a safe minimum these days but you might have more options if you can use AVX or AVX2. \$\endgroup\$ – mattnewport Oct 30 '14 at 22:20
  • \$\begingroup\$ @mattnewport Here's the list of supported instruction sets: MMX, SSE, SSE2, SSE3, SSSE3, SSE4.1, SSE4.2, EM64T, VT-x, AES, AVX. Maybe I should also mention that I'm working on Windows 7 with VS2010. \$\endgroup\$ – BConic Oct 31 '14 at 8:20
  • 2
    \$\begingroup\$ Just a suggestion if you want to get some meaningful answers to this question: when it comes to optimization you really need to test and profile. I'd be much more likely to take a stab at this if you linked to a VS project that set up a test environment for your function with some representative data. As it stands, it would be too much work for me to set up a test case before I could test some optimizations. \$\endgroup\$ – mattnewport Oct 31 '14 at 18:35
  • \$\begingroup\$ @mattnewport I just posted the link to the project ZIP file, thanks for the suggestion. \$\endgroup\$ – BConic Nov 1 '14 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.