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?
