Исчерпывающая билинейная обработка патчей - оптимизация 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 патчей, и кажется хорошей идеей обрабатывать пучки по 8 пикселей за один раз с использованием SSE.

Однако, я не знаю, как иметь дело с умножением на float интерполяции весов (т.е. w00, w01, w10 и w11. Эти веса обязательно положительно и меньше 1, следовательно, умножение не может переполнить unsigned char тип данных.

кто-нибудь знает, как поступить

---

EDIT:

я попытался сделать это следующим образом (при условии, 16х16 патчей), но нет никакой существенной скорости до:

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. (После соответствующего чередования данных и репликации значений V0 и V1 с помощью PUNPCKLBW и т.п.).

В конце разделить на общий вес (PSRLW), перемасштабировать до байтов (PACKUSWB). (Этот шаг может быть выполнен только один раз, объединяя две интерполяции.)

Вы могли бы подумать о удвоении всех весов, чтобы окончательное масштабирование было на 8 бит, а PACKUSWB было бы достаточно, но, к сожалению, оно насыщает значения и там не является ненасыщенным эквивалентом.

Может быть, лучше всего вычислить все 64 интерполяционных веса и суммировать четыре билинейных члена.

UPDATE:

Если цель состоит в том, чтобы интерполировать с фиксированными коэффициентами для всех пикселей каре (на самом деле достигающих субпикселей перевод), стратегия отличается.

Вы будете загружать пробег 8 (16?) пикселей, соответствующих верхним левым углам, пробег 8 сместил один пиксель вправо (соответствующий верхним правым углам) и аналогично для следующего ряда (нижние конусы); умножить и добавить пары (PMADDUBSW) значения пикселей в соответствующие значения интерполяции и объединить пары (PADDW). Храните веса с репликацией.

Другим вариантом является предотвращение (PMADD) и выполнение отдельных умножений (PMULLW) и добавление (PADDW). Это упростит схему реорганизации.

После масштабирования (как указано выше) вы закончите с 8 интерполированными значениями.

Это может работать и для переменных интерполяционных весов, если вы интерполируете ровно один пиксель на квадрат.