簡単に言うと、私は C++ でコンピューティング集約型の画像処理アプリケーションを開発しています。大きな画像から抽出されたピクセルの小さなブロックで、さまざまな画像ワープを計算する必要があります。プログラムが思ったほど速く実行されません。プロファイリング (OProfile) は、ワーピング/補間関数が CPU 時間の 70% 以上を消費することを示したので、それを試して最適化することは明らかでした。
今まで、タスクに OpenCV 画像処理ライブラリを使用していました。
// some parameters for the image warps (position, stretch, skew)
struct WarpParams;
void Image::get(const WarpParams ¶ms)
{
// fills matrices mapX_ and mapY_ with x and y coordinates of points to be
// inteprolated.
updateCoordMaps(params);
// perform interpolation to obtain pixels at point locations
// P(mapX_[i], mapY_[i]) based on the original data and put the
// result in pixels_. Use bicubic inteprolation.
cv::remap(image_->data(), pixels_, mapX_, mapY_, CV_INTER_CUBIC);
}
私は独自の補間関数を作成し、それをテスト ハーネスに入れて、実験中に正確性を確認し、古いものと比較してベンチマークしました。
私の機能は非常に遅くなりましたが、これは予想通りでした。一般的に、アイデアは次のとおりです。
- mapX_、mapY_ 座標マップを反復処理し、補間する次のピクセルの (実数値) 座標を抽出します。
- 補間されたピクセルを囲む元のイメージからピクセル値 (整数座標) の 4x4 近傍を取得します。
- これらの 16 ピクセルのそれぞれについて畳み込みカーネルの係数を計算します。
- 補間されたピクセルの値を、16 個のピクセル値とカーネル係数の線形結合として計算します。
Wolfdale Core2 Duo で 25us の古い機能を実行しました。新しいものは 587us (!) かかりました。私は熱心に魔法使いの帽子をかぶり、コードのハッキングを開始しました。すべてのブランチを削除し、重複する計算を省略し、3 つのネストされたループを座標マップ上で 1 つのループに変換することができました。これは私が思いついたものです：
void Image::getConvolve(const WarpParams ¶ms)
{
__declspec(align(16)) static float kernelX[4], kernelY[4];
// grab pointers to coordinate map matrices and the original image
const float
*const mapX = mapX_.ptr<float>(),
*const mapY = mapY_.ptr<float>(),
*const img = image_->data().ptr<float>();
// grab pointer to the output image
float *const subset = pixels_.ptr<float>(),
x, y, xint, yint;
const ptrdiff_t imgw = image_->width();
ptrdiff_t imgoffs;
__m128 v_px, v_kernX, v_kernY, v_val;
// iterate over the coordinate matrices as linear buffers
for (size_t idx = 0; idx < pxCount; ++idx)
{
// retrieve coordinates of next pixel from precalculated maps,
// break up each into fractional and integer part
x = modf(mapX[idx], &xint);
y = modf(mapY[idx], &yint);
// obtain offset of the top left pixel from the required 4x4
// neighborhood of the current pixel in the image's
// buffer (sadly, the position will be unaligned)
imgoffs = (((ptrdiff_t)yint - 1) * imgw) + (ptrdiff_t)xint - 1;
// calculate all 4 convolution kernel values for every row and
// every column
tap4Kernel(x, kernelX);
tap4Kernel(y, kernelY);
// load the kernel values for the columns, these don't change
v_kernX = _mm_load_ps(kernelX);
// process a row of the 4x4 neighborhood
// get set of convolution kernel values for the current row
v_kernY = _mm_set_ps1(kernelY[0]);
v_px = _mm_loadu_ps(img + imgoffs); // load the pixel values
// calculate the linear combination of the pixels with kernelX
v_px = _mm_mul_ps(v_px, v_kernX);
v_px = _mm_mul_ps(v_px, v_kernY); // and kernel Y
v_val = v_px; // add result to the final value
imgoffs += imgw;
// offset points now to next row of the 4x4 neighborhood
v_kernY = _mm_set_ps1(kernelY[1]);
v_px = _mm_loadu_ps(img + imgoffs);
v_px = _mm_mul_ps(v_px, v_kernX);
v_px = _mm_mul_ps(v_px, v_kernY);
v_val = _mm_add_ps(v_val, v_px);
imgoffs += imgw;
/*... same for kernelY[2] and kernelY[3]... */
// store resulting interpolated pixel value in the subset's
// pixel matrix
subset[idx] = horizSum(v_val);
}
}
// Calculate all 4 values of the 4-tap convolution kernel for 4 neighbors
// of a pixel and store them in an array. Ugly but fast.
// The "arg" parameter is the fractional part of a pixel's coordinate, i.e.
// a number in the range <0,1)
void Image::tap4Kernel(const float arg, float *out)
{
// chaining intrinsics was slower, so this is done in separate steps
// load the argument into 4 cells of a XMM register
__m128
v_arg = _mm_set_ps1(arg),
v_coeff = _mm_set_ps(2.0f, 1.0f, 0.0f, -1.0f);
// subtract vector of [-1, 0, 1, 2] to obtain coorinates of 4 neighbors
// for kernel calculation
v_arg = _mm_sub_ps(v_arg, v_coeff);
// clear sign bits, this is equivalent to fabs() on all 4
v_coeff = _mm_set_ps1(-0.f);
v_arg = _mm_andnot_ps(v_coeff, v_arg);
// calculate values of abs(argument)^3 and ^2
__m128
v_arg2 = _mm_mul_ps(v_arg, v_arg),
v_arg3 = _mm_mul_ps(v_arg2, v_arg),
v_val, v_temp;
// calculate the 4 kernel values as
// arg^3 * A + arg^2 * B + arg * C + D, using
// (A,B,C,D) = (-0.5, 2.5, -4, 2) for the outside pixels and
// (1.5, -2.5, 0, 1) for inside
v_coeff = _mm_set_ps(-0.5f, 1.5f, 1.5f, -0.5f);
v_val = _mm_mul_ps(v_coeff, v_arg3);
v_coeff = _mm_set_ps(2.5f, -2.5f, -2.5f, 2.5f);
v_temp = _mm_mul_ps(v_coeff, v_arg2);
v_val = _mm_add_ps(v_val, v_temp);
v_coeff = _mm_set_ps(-4.0f, 0.0f, 0.0f, -4.0f),
v_temp = _mm_mul_ps(v_coeff, v_arg);
v_val = _mm_add_ps(v_val, v_temp);
v_coeff = _mm_set_ps(2.0f, 1.0f, 1.0f, 2.0f);
v_val = _mm_add_ps(v_val, v_coeff);
_mm_store_ps(out, v_val);
}
最後に取っておいたメインループに SSE を導入する前でさえ、実行時間を 40us 未満にすることができてうれしかったです。少なくとも 3 倍のスピードアップを期待していましたが、36us でかろうじて速くなり、改善しようとしていた古い get() よりも遅くなりました。さらに悪いことに、ベンチマーク ループを変更してより多くの実行を行ったところ、古い関数の平均実行時間は同じでしたが、私のものは 127us を超えました。ワープとは、結果を計算するために元の画像から広く分散したピクセル値に到達する必要があることを意味します)。
その理由はロードが整列されていないことに違いないと考えましたが、それは仕方がありません (予測できないピクセル値に到達する必要があります)。最適化部門でこれ以上できることは何もなかったので、cv::remap() 関数を調べて、彼らがどのようにそれを行っているかを確認することにしました。入れ子になったループとたくさんの分岐が含まれていることに驚いたことを想像してみてください。また、私が気にする必要のない多くの引数の検証も行います。私が知る限り（コードにコメントはありません）、SSE（アライメントされていないロードも同様です！）は、座標マップから値を抽出して整数に丸めるためにのみ使用され、その後、実際の補間を行う関数が呼び出されます通常の float 演算を使用します。
私の質問は、なぜ私は惨めに失敗したのですか (なぜ私のコードはとても遅く、混乱しているように見えてもなぜ彼らのコードは速いのですか)、コードを改善するにはどうすればよいでしょうか?
ここでは OpenCV コードを貼り付けません。これは既に長すぎるためです。pastebinで確認できます。
VC++2010 のリリース モードでコードをテストし、コンパイルしました。使用した OpenCV は、v2.3.1 のコンパイル済みバイナリ バンドルです。
編集:ピクセル値は 0..1 の範囲の浮動小数点数です。プロファイリングは、tap4Kernel() 関数が関連していないことを示しました。ほとんどの時間は getConvolve() 内で費やされます。
EDIT2:生成されたコードの逆アセンブリをペーストビンに貼り付けました。これは、古い Banias Celeron プロセッサ (SSE2 を搭載) でコンパイルされていますが、外観はほぼ同じです。
EDIT3：すべてのプログラマーがメモリについて知っておくべきことを読んだ後、OpenCV関数が多かれ少なかれ同じアルゴリズムを実装していると誤って想定していることに気付きました。補間するピクセルごとに、その 4x4 の近傍を取得する必要があります。この近傍のピクセルは、画像バッファー内に不連続に配置されています。私は CPU キャッシュを悪用していますが、OpenCV はおそらくそうではありません。私の関数には 5,800,000 回のメモリ アクセスがあるのに対し、OpenCV は 400,000 回しかアクセスしないため、VTune プロファイリングは一致しているように見えます。それらの機能はごちゃごちゃしており、おそらくさらに最適化される可能性がありますが、おそらくメモリとキャッシュの使用に対するよりスマートなアプローチのおかげで、私よりも優位に立つことができます.
アップデート：ピクセル値を XMM レジスタにロードする方法を改善することができました。画像のすべてのピクセルに対して 16 要素のセルを保持するオブジェクトにバッファーを割り当てます。画像の読み込み時に、このセル バッファーに、ピクセルごとに事前に配置された 4x4 近傍のシーケンスを入力します。スペース効率はあまり良くありませんが (画像は 16 倍のスペースを必要とします)、ロードは常に整列され (_mm_loadu_ps() はもう必要ありません)、画像バッファーからピクセルを分散して読み取る必要がなくなります。順次格納されます。驚いたことに、ほとんど改善がありませんでした。アラインされていないロードは 10 倍遅くなる可能性があると聞きましたが、明らかにこれは問題ではありません。しかし、コードの一部をコメントアウトすると、modf() 呼び出しが実行時間の 75% を占めていることがわかりました。それらを排除することに焦点を当て、回答を投稿します。