优化矩阵旋转 - 关于矩阵中心的任意角度
Posted
技术标签:
【中文标题】优化矩阵旋转 - 关于矩阵中心的任意角度【英文标题】：Optimized Matrix Rotation - arbitrary angle about center of matrix 【发布时间】：2011-08-04 14:08:16 【问题描述】：我正在尝试优化超大图像的旋转，最小的是 4096x4096 或约 1600 万像素。
旋转总是围绕图像的中心，图像不一定总是正方形，但总是 2 的幂。
我可以访问 MKL/TBB，其中 MKL 是针对我的目标平台优化的 BLAS。我不熟悉这个操作是否在 BLAS 中。
到目前为止，对于 4096x4096 图像，我的最佳尝试大约是 17 到 25 毫秒（对于相同的图像大小非常不一致，这意味着我可能会在整个缓存中踩踏）。矩阵是 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 ) );
fcomplex 只是复数的内部表示。定义为：
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： 我在下面提出了一个解决方案，考虑了我作为答案得到的建议以及在旋转矩形时修复了一个错误。
【问题讨论】：
我edited out the clutter 第二次了。 这为 GPU 计算而尖叫。也许这是你的一个选择。 我记得以前使用Bresenham's algorithm 来避免在这种情况下进行浮点运算。这可能是一个想法，因为您的 xOffset 和 yOffset 修改似乎是 Bresenham 的不错选择。 @Axel Gneiting 不幸的是，我还不能把它带到 OpenCL 或 CUDA 世界。我们不保证客户的硬件。但是，我已经获得了下一次迭代的许可，可以为具有兼容硬件的客户编写 gpu 版本。到那时！ 【参考方案1】：如果您首先执行一个简单的近似旋转 (90/190/270) 度，然后在 0-90 度之间进行最终旋转，您也许可以优化很多东西。例如。然后您可以优化if( srcX >= 0 && srcX < colSpan && srcY >= 0 && srcY < rowSpan ) 测试，它会更加缓存友好。我敢打赌，您的分析表明 91 度旋转比 1 度旋转慢得多。
【讨论】：
它确实显示了这一点，这里有 20 个输出时序： 每张幻灯片滑动图像大约 7.5 度，时间以毫秒为单位：0.0274989 - 7.5 deg 0.0265466 - 15 deg 0.0246373 - 22.5 deg 0.0245368 - 30 deg 0.0246203 - 37.5 deg 0.0241758 - 45 deg 0.0232378 - 52.5 deg 0.0216996 - 60 deg 0.0233684 - 67.5 deg 0.021255 - 75 deg 0.0179 - 82.5 deg 0.0165404 - 90 deg 0.0179594 - 97.5 deg 0.0185623 - 105 deg 0.0203935 - 112.5 deg 0.0218039 - 120 deg 0.0243723 - 127.5 deg 0.0228118 - 135 deg 0.023037 - 142.5 deg
顺便说一句，要将大图像旋转 90 度以上，一次旋转 16x16 个子块。这应该让缓存满意。【参考方案2】：
没有太多需要优化的地方。你的算法是合理的。您正在逐行写入 dstData（这有利于缓存/内存）强制每个线程顺序写入。
唯一剩下的就是循环展开你的内部 for...loop ~4x（对于 32 位系统）或 8x（对于 64 位系统）。它可能会为您带来大约 10-20% 的改进。由于问题的性质（强制从 srcData 进行随机读取），您总是会有时间上的差异。
我会进一步思考......
您的内部 for...循环是矢量化的强大目标。 考虑静态向量：
// 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
x86 的 SSE 指令非常适合处理此类数据。甚至从双精度数到整数的转换。允许 256 位向量（4 个双精度）的 AVX 指令更适合。
【讨论】：
我最终转换为浮点数，因为无论如何我都会使用整数。然后我在花车上做了 SIMD（一次做 4 个很好！）我的时间大约是原来的 60%！时间以毫秒为单位，20转：0.0197117 987654323 0.0198210 987654325 0.0170747 987654327 0.0171950 987654329 0.0165822 987654331 0.0153101 987654333 0.0146376 987654335 @ 0.01569760.01666880.01604940.01719730.01820940.0184469【参考方案3】：
考虑到所提供的建议，我得出了这个解决方案。此外，我修复了原始实现中的一个错误，该错误导致矩形图像出现问题。
首先旋转 90 度的建议（使用仿射变换和线程化，然后旋转较小的度数被证明是因为必须迭代矩阵两次而更慢）。当然，很多因素都会影响到这一点，而且很可能内存带宽会导致事情变得更加扭曲。因此，对于我正在测试和优化的机器，这个解决方案被证明是我能提供的最好的。
使用 16x16 的图块：
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% 肯定这是最好的答案。但是，这是我能提供的最好的。所以，我想我会将我的解决方案传递给社区。p>
【讨论】：
以上是关于优化矩阵旋转 - 关于矩阵中心的任意角度的主要内容，如果未能解决你的问题，请参考以下文章