4

私は非常に大きな画像の回転を最適化しようとしています。最小は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:回答として得た提案を考慮し、長方形を回転するときのバグを修正して、以下に解決策を示します。

4

3 に答える 3

3

最初に単純な近似回転(90/190/270)度を実行し、次に0〜90度の間の最終回転を実行すると、かなりの数のことを最適化できる可能性があります。たとえば、if( srcX >= 0 && srcX < colSpan && srcY >= 0 && srcY < rowSpan )テストを最適化すると、キャッシュがより使いやすくなります。あなたのプロファイリングは、91度の回転が1度の回転よりもはるかに遅いことを示しているに違いありません。

于 2011-08-04T15:00:27.833 に答える
1

最適化することはあまりありません。あなたのアルゴリズムは健全です。行ごとに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命令は、さらに適しています。

于 2011-08-04T14:52:54.747 に答える
1

提供された提案を考慮に入れて、私はこの解決策にたどり着きました。また、長方形の画像で問題が発生する元の実装のバグを修正しました。

最初に90度回転するという提案(アフィン変換を使用し、それをスレッド化してから、少しだけ回転させると、マトリックスを2回繰り返す必要があるだけで遅くなることがわかりました)。もちろん、それには多くの要因が関係しており、おそらくメモリ帯域幅が原因で物事がより歪んでいます。したがって、私がテストおよび最適化するマシンの場合、このソリューションは私が提供できる最高のものであることが証明されました。

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%ポジティブではありません。これが最良の答えです。しかし、これは私が提供できた最高のものです。それで、私は自分の解決策をコミュニティに渡すと思いました。

于 2011-08-11T19:11:56.373 に答える