Я пытаюсь оптимизировать вращение очень больших изображений, наименьшее - 4096x4096 или ~ 16 миллионов пикселей.
Вращение всегда происходит вокруг центра изображения, и изображения не обязательно всегда квадратные, но всегда будут степенью двойки.
У меня есть доступ к MKL/TBB, где MKL — это оптимизированный BLAS для моих целевых платформ. Я не знаю, есть ли эта операция вообще в BLAS.
Мои лучшие попытки на данный момент составляют около 17-25 мс (очень непоследовательно для одного и того же размера изображения, что означает, что я, вероятно, топчу весь кеш) для изображений 4096x4096. Матрицы выровнены по 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: я разместил решение ниже, принимая во внимание предложения, которые я получил в качестве ответов, а также исправление ошибки при вращении прямоугольников.