灰度双线性补丁提取 - SSE 优化

Question

我的程序大量使用通过双线性插值从较大灰度图像中提取的小子图像。

我为此目的使用了以下函数：

bool extract_patch_bilin(const cv::Point2f &patch_ctr, const cv::Mat_<uchar> &img, cv::Mat_<uchar> &patch)
{
    const int hsize = patch.rows/2;

    // ...
    // Precondition checks: patch is a preallocated square matrix and both patch and image have continuous buffers
    // ...

    int floorx=(int)floor(patch_ctr.x)-hsize, floory=(int)floor(patch_ctr.y)-hsize;
    if(floorx<0 || img.cols-1<floorx+patch.cols || floory<0 || img.rows-1<floory+patch.rows)
        return false;

    float x=patch_ctr.x-hsize-floorx;
    float y=patch_ctr.y-hsize-floory;
    float xy = x*y;
    float w00=1-x-y+xy, w01=x-xy, w10=y-xy, w11=xy;
    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)
            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;
}

长话短说，我已经在程序的其他部分使用并行化，我正在考虑使用 SSE 来优化这个函数的执行，因为我主要使用 8x8 补丁，处理束似乎是个好主意使用 SSE 一次 8 个像素。

但是，我不确定如何处理乘以 float 插值权重（即 w00、w01、w10 和 w11。这些权重必须为正且小于 1，因此乘法不能溢出unsigned char 数据类型。

有人知道怎么做吗？

---

编辑：

我尝试按如下方式执行此操作（假设 16x16 补丁），但没有明显的加速：

bool extract_patch_bilin_16x16(const cv::Point2f& 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
    float x=patch_ctr.x-hsize-floorx;
    float  y=patch_ctr.y-hsize-floory;
    float  xy = x*y;
    float  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;
    // Precompute weighting 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
    int ngroups = patch.rows>>4;
    for(int v=0; v<patch.rows; ++v,buff_img0+=img_stride,buff_img1+=img_stride) {
        for(int g=0; g<ngroups; ++g,buff_patch+=16,buff_img0+=16,buff_img1+=16) {
                ////////////////////////////////
                // 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 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 approximate 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);
                ////////////////////////////////
                // Process the higher 8 values
                ////////////////////////////////
                // Unpack into 16-bits integers
                __m128i val00_hi = _mm_unpackhi_epi8(val00,CONST_0);
                __m128i val01_hi = _mm_unpackhi_epi8(val01,CONST_0);
                __m128i val10_hi = _mm_unpackhi_epi8(val10,CONST_0);
                __m128i val11_hi = _mm_unpackhi_epi8(val11,CONST_0);
                // Multiply with the integer weights
                __m128i w256val00_hi = _mm_mullo_epi16(val00_hi,w00x256_16i);
                __m128i w256val01_hi = _mm_mullo_epi16(val01_hi,w01x256_16i);
                __m128i w256val10_hi = _mm_mullo_epi16(val10_hi,w10x256_16i);
                __m128i w256val11_hi = _mm_mullo_epi16(val11_hi,w11x256_16i);
                // Divide by 256 to get the approximate result of the multiplication with floating-point weights
                __m128i wval00_hi = _mm_srli_epi16(w256val00_hi,8);
                __m128i wval01_hi = _mm_srli_epi16(w256val01_hi,8);
                __m128i wval10_hi = _mm_srli_epi16(w256val10_hi,8);
                __m128i wval11_hi = _mm_srli_epi16(w256val11_hi,8);
                // Add pairwise
                __m128i sum0_hi = _mm_add_epi16(wval00_hi,wval01_hi);
                __m128i sum1_hi = _mm_add_epi16(wval10_hi,wval11_hi);
                __m128i final_hi = _mm_add_epi16(sum0_hi,sum1_hi);
                ////////////////////////////////
                // Repack all values
                ////////////////////////////////
                __m128i final_val = _mm_packus_epi16(final_lo,final_hi);
                _mm_storeu_si128((__m128i*)buff_patch,final_val);
        }
    }
}

知道可以做些什么来提高速度吗？

Answer 1

我会考虑坚持使用整数：您的权重是 1/64 的倍数，因此使用定点 8.6 就足够了，并且适合 16 位数字。

双线性插值最好作为三个线性插值（两个在 Y 上，然后在 X 上一个；您可以将第二个 Y 插值重复用于相邻的补丁）。

要在两个值之间执行线性插值，您将为所有插值权重 P 和 Q（8 到 1 和 0 到 7）预存储一次，然后将它们成对相乘和相加，例如 V0.P[i] +V1.Q[i]。这可以使用 PMADDUBSW 指令有效地完成。 （经过适当的数据交织，并使用 PUNPCKLBW 等复制值 V0 和 V1）。

最后，除以总权重（PSRLW），重新缩放为字节（PACKUSWB）。 （这一步只能执行一次，结合两次插值。）

您可以考虑将所有权重加倍，这样最终的缩放比例为 8 位，PACKUSWB 就足够了，但不幸的是，它使值饱和并且没有不饱和的等价物。

可能是预先计算所有 64 个插值权重并对四个双线性项求和更好。

更新：

如果目标是对所有像素四边形使用固定系数进行插值（实际上实现亚像素平移），则策略不同。

您将加载与左上角对应的 8 (16 ?) 个像素的运行，将 8 个像素向右移动一个像素（对应于右上角）的运行，对于下一行 (底角）；将像素值成对相乘和相加 (PMADDUBSW) 到相应的插值权重，然后组合成对 (PADDW)。通过复制存储权重。

另一个选项是避免 (PMADD) 并执行单独的乘法 (PMULLW) 和加法 (PADDW)。这将简化重组方案。

缩放后（如上），您最终会得到 8 个插值。

这也适用于可变插值权重，只要每个四边形恰好插值一个像素。