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?
Grayscale bilinear patch extraction - SSE optimization
c++;optimization;image;sse;opencv