優化旋轉矩陣 - 關於矩陣

Question

我試圖優化非常大的圖像的旋轉中心任意角度，最小的是4096×4096或〜1600萬像素。

旋轉總是關於圖像的中心和圖像不一定總是正方形的，但永遠是2

我有機會獲得MKL/TBB，其中MKL是一個優化的BLAS我的目標動力平臺。我不熟悉這個操作是否在BLAS中。所以，裸露在我面前的答案是顯而易見的，只是我不熟悉的BLAS功能。

我盡了最大努力，到目前爲止都是圍繞17-25ms爲4096×4096的圖像（相同的圖像尺寸，這意味着我可能踩遍緩存很不一致）。矩陣是16字節對齊的。現在

，目標不能調整大小。所以，裁剪應該會發生。例如，以45度角旋轉的矩形矩陣必定會夾在角上，該位置的值應爲零。

現在，我盡了最大努力用瓷磚的做法 - 不優雅是尚未被放入片大小也不循環展開。

這裏是我的算法，因爲它代表利用TBB - http://threadingbuildingblocks.org/：

//- cosa = cos of the angle 
//- sina = sin of angle 
//- for those unfamiliar with TBB, this is giving me blocks of rows or cols that 
//- are unique per thread 
void operator() (const tbb::blocked_range2d<size_t, size_t> r) const 
{ 
    double xOffset; 
    double yOffset; 
    int lineOffset; 

    int srcX; 
    int srcY; 

    for (size_t row = r.rows().begin(); row != r.rows().end(); ++row) 
    { 
     const size_t colBegin = r.cols().begin(); 
     xOffset = -(row * sina) + xHelper + (cosa * colBegin); 
     yOffset = (row * cosa) + yHelper + (sina * colBegin); 
     lineOffset = (row * rowSpan); //- all col values are offsets of this row 
     for(size_t col = colBegin; col != r.cols().end(); ++col, xOffset += cosa, yOffset += sina) 
     { 
      srcX = xOffset; 
      srcY = yOffset; 

      if(srcX >= 0 && srcX < colSpan && srcY >= 0 && srcY < rowSpan) 
      { 
       destData[col + lineOffset] = srcData[srcX + (srcY * rowSpan)]; 
      } 
     } 
    } 
} 

我對這個函數的調用是這樣的：

double sina = sin(angle); 
double cosa = cos(angle); 
double centerX = (colSpan)/2; 
double centerY = (rowSpan)/2; 

//- Adding .5 for rounding 
const double xHelper = centerX - (centerX * cosa) + (centerY * sina) + .5; 
const double yHelper = centerY - (centerX * sina) - (centerY * cosa) + .5; 
tbb::parallel_for(tbb::blocked_range2d<size_t, size_t>(0, rowSpan, 0, colSpan), DoRotate(sina, cosa, xHelper, yHelper, rowSpan, colSpan, (fcomplex *)pDestData, (fcomplex *)pSrcData)); 

浮點型複數僅僅是在複數的房子表示。它被定義爲：

struct fcomplex 
{ 
    float real; 
    float imag; 
}; 

所以，我想要做複數值矩陣的繞它的中心在任意角度非常大的圖像儘可能快。

更新：

基於奇妙的反饋，我已經更新到這一點：這是增加約40％。我想知道，但如果什麼都可以做：

void operator() (const tbb::blocked_range2d<size_t, size_t> r) const 
{ 
    float xOffset; 
    float yOffset; 
    int lineOffset; 

    __m128i srcXints; 
    __m128i srcYints; 

    __m128 dupXOffset; 
    __m128 dupYOffset; 

    for (size_t row = r.rows().begin(); row != r.rows().end(); ++row) 
    { 
     const size_t colBegin = r.cols().begin(); 
     xOffset = -(row * sina) + xHelper + (cosa * colBegin); 
     yOffset = (row * cosa) + yHelper + (sina * colBegin); 
     lineOffset = (row * rowSpan); //- all col values are offsets of this row 

     for(size_t col = colBegin; col != r.cols().end(); col+=4, xOffset += dupOffsetsX.m128_f32[3], yOffset += dupOffsetsY.m128_f32[3]) 
     { 
      dupXOffset = _mm_load1_ps(&xOffset); //- duplicate the x offset 4 times into a 4 float field 
      dupYOffset = _mm_load1_ps(&yOffset); //- duplicate the y offset 4 times into a 4 float field 

      srcXints = _mm_cvtps_epi32(_mm_add_ps(dupOffsetsX, dupXOffset)); 
      srcYints = _mm_cvtps_epi32(_mm_add_ps(dupOffsetsY, dupYOffset)); 

      if(srcXints.m128i_i32[0] >= 0 && srcXints.m128i_i32[0] < colSpan && srcYints.m128i_i32[0] >= 0 && srcYints.m128i_i32[0] < rowSpan) 
      { 
       destData[col + lineOffset] = srcData[srcXints.m128i_i32[0] + (srcYints.m128i_i32[0] * rowSpan)]; 
      } 

      if(srcXints.m128i_i32[1] >= 0 && srcXints.m128i_i32[1] < colSpan && srcYints.m128i_i32[1] >= 0 && srcYints.m128i_i32[1] < rowSpan) 
      { 
       destData[col + 1 + lineOffset] = srcData[srcXints.m128i_i32[1] + (srcYints.m128i_i32[1] * rowSpan)]; 
      } 

      if(srcXints.m128i_i32[2] >= 0 && srcXints.m128i_i32[2] < colSpan && srcYints.m128i_i32[2] >= 0 && srcYints.m128i_i32[2] < rowSpan) 
      { 
       destData[col + 2 + lineOffset] = srcData[srcXints.m128i_i32[2] + (srcYints.m128i_i32[2] * rowSpan)]; 
      } 

      if(srcXints.m128i_i32[3] >= 0 && srcXints.m128i_i32[3] < colSpan && srcYints.m128i_i32[3] >= 0 && srcYints.m128i_i32[3] < rowSpan) 
      { 
       destData[col + 3 + lineOffset] = srcData[srcXints.m128i_i32[3] + (srcYints.m128i_i32[3] * rowSpan)]; 
      } 
     } 
    } 
} 

更新2： 我把下面的解決方案，同時考慮到建議我得到的答案以及旋轉矩形時修復bug。

Answer 1

考慮到所提供的建議，我在這個解決方案已經到達。另外，我修復了原始實現中的一個錯誤，導致矩形圖像出現問題。

（只是從具有遍歷矩陣兩次使用仿射變換和線程的是，再經證明是較慢的更小的度旋轉）90度旋轉首先的建議。當然，有很多因素會影響到這一點，而最有可能的內存帶寬會導致更多的偏差。所以，對於我正在測試和優化的機器來說，這個解決方案被證明是我能提供的最好的解決方案。

使用16×16瓦：

class DoRotate 
{ 
const double sina; 
const double cosa; 

const double xHelper; 
const double yHelper; 

const int rowSpan; 
const int colSpan; 

mutable fcomplex *destData; 
const fcomplex *srcData; 

const float *offsetsX; 
const float *offsetsY; 

__m128 dupOffsetsX; 
__m128 dupOffsetsY; 

public: 
void operator() (const tbb::blocked_range2d<size_t, size_t> r) const 
{ 
    float xOffset; 
    float yOffset; 
    int lineOffset; 

    __m128i srcXints; 
    __m128i srcYints; 

    __m128 dupXOffset; 
    __m128 dupYOffset; 

    for (size_t row = r.rows().begin(); row != r.rows().end(); ++row) 
    { 
     const size_t colBegin = r.cols().begin(); 
     xOffset = -(row * sina) + xHelper + (cosa * colBegin); 
     yOffset = (row * cosa) + yHelper + (sina * colBegin); 
     lineOffset = (row * colSpan); //- all col values are offsets of this row 

     for(size_t col = colBegin; col != r.cols().end(); col+=4, xOffset += (4 * cosa), yOffset += (4 * sina)) 
     { 
      dupXOffset = _mm_load1_ps(&xOffset); //- duplicate the x offset 4 times into a 4 float field 
      dupYOffset = _mm_load1_ps(&yOffset); //- duplicate the y offset 4 times into a 4 float field 

      srcXints = _mm_cvttps_epi32(_mm_add_ps(dupOffsetsX, dupXOffset)); 
      srcYints = _mm_cvttps_epi32(_mm_add_ps(dupOffsetsY, dupYOffset)); 

      if(srcXints.m128i_i32[0] >= 0 && srcXints.m128i_i32[0] < colSpan && srcYints.m128i_i32[0] >= 0 && srcYints.m128i_i32[0] < rowSpan) 
      { 
       destData[col + lineOffset] = srcData[srcXints.m128i_i32[0] + (srcYints.m128i_i32[0] * colSpan)]; 
      } 

      if(srcXints.m128i_i32[1] >= 0 && srcXints.m128i_i32[1] < colSpan && srcYints.m128i_i32[1] >= 0 && srcYints.m128i_i32[1] < rowSpan) 
      { 
       destData[col + 1 + lineOffset] = srcData[srcXints.m128i_i32[1] + (srcYints.m128i_i32[1] * colSpan)]; 
      } 

      if(srcXints.m128i_i32[2] >= 0 && srcXints.m128i_i32[2] < colSpan && srcYints.m128i_i32[2] >= 0 && srcYints.m128i_i32[2] < rowSpan) 
      { 
       destData[col + 2 + lineOffset] = srcData[srcXints.m128i_i32[2] + (srcYints.m128i_i32[2] * colSpan)]; 
      } 

      if(srcXints.m128i_i32[3] >= 0 && srcXints.m128i_i32[3] < colSpan && srcYints.m128i_i32[3] >= 0 && srcYints.m128i_i32[3] < rowSpan) 
      { 
       destData[col + 3 + lineOffset] = srcData[srcXints.m128i_i32[3] + (srcYints.m128i_i32[3] * colSpan)]; 
      } 
     } 
    } 
} 

DoRotate(const double pass_sina, const double pass_cosa, const double pass_xHelper, const double pass_yHelper, 
      const int pass_rowSpan, const int pass_colSpan, const float *pass_offsetsX, const float *pass_offsetsY, 
      fcomplex *pass_destData, const fcomplex *pass_srcData) : 
    sina(pass_sina), cosa(pass_cosa), xHelper(pass_xHelper), yHelper(pass_yHelper), 
    rowSpan(pass_rowSpan), colSpan(pass_colSpan), 
    destData(pass_destData), srcData(pass_srcData) 
{ 
    dupOffsetsX = _mm_load_ps(pass_offsetsX); //- load the offset X array into one aligned 4 float field 
    dupOffsetsY = _mm_load_ps(pass_offsetsY); //- load the offset X array into one aligned 4 float field 
} 
}; 

，然後調用代碼：

double sina = sin(radians); 
double cosa = cos(radians); 

double centerX = (colSpan)/2; 
double centerY = (rowSpan)/2; 

//- Adding .5 for rounding to avoid periodicity 
const double xHelper = centerX - (centerX * cosa) + (centerY * sina) + .5; 
const double yHelper = centerY - (centerX * sina) - (centerY * cosa) + .5; 

float *offsetsX = (float *)_aligned_malloc(sizeof(float) * 4, 16); 
offsetsX[0] = 0.0f; 
offsetsX[1] = cosa; 
offsetsX[2] = cosa * 2.0; 
offsetsX[3] = cosa * 3.0; 
float *offsetsY = (float *)_aligned_malloc(sizeof(float) * 4, 16); 
offsetsY[0] = 0.0f; 
offsetsY[1] = sina; 
offsetsY[2] = sina * 2.0; 
offsetsY[3] = sina * 3.0; 

//- tiled approach. Works better, but not by much. A little more stays in cache 
tbb::parallel_for(tbb::blocked_range2d<size_t, size_t>(0, rowSpan, 16, 0, colSpan, 16), DoRotate(sina, cosa, xHelper, yHelper, rowSpan, colSpan, offsetsX, offsetsY, (fcomplex *)pDestData, (fcomplex *)pSrcData)); 

_aligned_free(offsetsX); 
_aligned_free(offsetsY); 

我在沒有辦法100％肯定的，這是最好的答案。但是，這是我能夠提供的最好的。所以，我想我會把我的解決方案傳遞給社區。

Answer 2

沒有太多優化。你的算法是健全的。您正在按行寫入dstData（這對緩存/內存有好處），從而強制每個線程進行順序寫入。

的唯一剩下的就是循環展開你內心的......循環〜4倍（對於32位系統）或8倍（對於64位系統）。它可能會爲你提供10-20％的改善。由於問題的性質（強制從srcData中隨機讀取），您總是會有時間差異。

我將進一步深思......

你內心的...循環是量化很強的目標。 考慮靜態載體：

// SSE instructions MOVDDUP (move and duplicate) MULPD (multiply packed double) 
double[] vcosa = [cosa, cosa, cosa, cosa] * [1.0, 2.0, 3.0, 4.0] 
double[] vsina = [sina, sina, sina, sina] * [1.0, 2.0, 3.0, 4.0] 

vxOffset = [xOffset, xOffset, xOffset, xOffset] 
vyOffset = [yOffset, yOffset, yOffset, yOffset] 

// SSE instructions ADDPD (add packed double) and CVTPD2DQ (convert packed double to signed integer) 
vsrcX = vcosa + vxOffset 
vsrcY = vsina + vyOffset 

86的SSE指令集非常適合用於處理這種類型的數據。即使從雙打轉換爲整數。 AVX指令允許256位向量（4倍）更適合。

Answer 3

您可能能夠優化不少東西，如果你先進行0-90度之間簡單的近似旋轉（90/190/270）度，然後最後的旋轉。例如。那麼您可以優化if(srcX >= 0 && srcX < colSpan && srcY >= 0 && srcY < rowSpan)測試，並且它將更加緩存友好。我敢打賭，你的分析表明91度旋轉比1度旋轉要慢很多。