优化旋转矩阵 - 关于矩阵

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度旋转要慢很多。