目的

OpenCVによるビデオ入力と類似度の測定チュートリアルでは、2つの画像間の類似度を確認するためのPSNRおよびSSIM法をすでに紹介した。ご覧のとおり、その実行処理にはかなりの時間がかかり、特にSSIMの場合は顕著である。しかし、CPU向けのOpenCV実装の性能値に満足できず、かつシステムにNVIDIA CUDA GPUデバイスがある場合でも、まだ望みはある。ビデオカード向けに独自のアルゴリズムを移植したり記述したりしてみるとよい。

このチュートリアルは、OpenCVのGPUモジュールを使ったコーディングへの取り組み方をよく理解するのに役立つ。前提として、core、highgui、imgproc モジュールの扱い方をすでに知っている必要がある。したがって、主な目的は次のとおりである。

CPUと比べて何が違うのか?
PSNRとSSIMのGPUコードを作成する
最大の性能が得られるようにコードを最適化する

ソースコード

ソースコードと映像ファイルは、OpenCVソースライブラリの samples/cpp/tutorial_code/gpu/gpu-basics-similarity/gpu-basics-similarity ディレクトリにもあり、こちらからダウンロードすることもできる。完全なソースコードはかなり長い（コマンドライン引数によるアプリケーションの制御と性能測定のため）。そのため、これらの節がそれらで散らかるのを避けるため、ここでは関数自体のみを示す。

PSNRは浮動小数点数を返し、2つの入力が類似している場合は30から50の値となる（高いほど良い）。

double getPSNR(const Mat& I1, const Mat& I2)
{
 Mat s1;
 absdiff(I1, I2, s1);       // |I1 - I2|
    s1.convertTo(s1, CV_32F);  // cannot make a square on 8 bits
    s1 = s1.mul(s1);           // |I1 - I2|^2
 
 Scalar s = sum(s1);         // sum elements per channel
 
 double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels
 
 double  mse = sse /(double)(I1.channels() * I1.total());
 
 // For very small SSE, add epsilon to approximate infinite PSNR (~361 dB)
 if( sse <= 1e-10) mse+= DBL_EPSILON;
 
 double psnr = 10.0*log10((255*255)/mse);
 return psnr;
}

double getPSNR_CUDA(const Mat& I1, const Mat& I2)
{
 cuda::GpuMat gI1, gI2, gs, t1,t2;
 
    gI1.upload(I1);
    gI2.upload(I2);
 
    gI1.convertTo(t1, CV_32F);
    gI2.convertTo(t2, CV_32F);
 
    cuda::absdiff(t1.reshape(1), t2.reshape(1), gs);
    cuda::multiply(gs, gs, gs);
 
 Scalar s = cuda::sum(gs);
 double sse = s.val[0] + s.val[1] + s.val[2];
 
 double  mse = sse /(double)(gI1.channels() * I1.total());
 
 // For very small SSE, add epsilon to approximate infinite PSNR (~361 dB)
 if( sse <= 1e-10) mse+= DBL_EPSILON;
 
 double psnr = 10.0*log10((255*255)/mse);
 return psnr;
}

struct BufferPSNR                                     // Optimized CUDA versions
{   // Data allocations are very expensive on CUDA. Use a buffer to solve: allocate once reuse later.
 cuda::GpuMat gI1, gI2, gs, t1,t2;
 
 cuda::GpuMat buf;
};

double getPSNR_CUDA_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)
{
    b.gI1.upload(I1);
    b.gI2.upload(I2);
 
    b.gI1.convertTo(b.t1, CV_32F);
    b.gI2.convertTo(b.t2, CV_32F);
 
    cuda::absdiff(b.t1.reshape(1), b.t2.reshape(1), b.gs);
    cuda::multiply(b.gs, b.gs, b.gs);
 
 double sse = cuda::sum(b.gs, b.buf)[0];
 
 double mse = sse /(double)(I1.channels() * I1.total());
 
 // For very small SSE, add epsilon to approximate infinite PSNR (~361 dB)
 if (sse <= 1e-10) mse += DBL_EPSILON;
 
 double psnr = 10.0*log10((255*255)/mse);
 return psnr;
}

SSIMは画像のMSSIMを返す。これも0から1の間の浮動小数点数だが（高いほど良い）、チャンネルごとに1つの値を持つ。そのため、OpenCVのデータ構造 Scalar を返す。

Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
 const double C1 = 6.5025, C2 = 58.5225;
    /***************************** INITS **********************************/
 int d     = CV_32F;
 
 Mat I1, I2;
    i1.convertTo(I1, d);           // cannot calculate on one byte large values
    i2.convertTo(I2, d);
 
 Mat I2_2   = I2.mul(I2);        // I2^2
 Mat I1_2   = I1.mul(I1);        // I1^2
 Mat I1_I2  = I1.mul(I2);        // I1 * I2
 
    /*************************** END INITS **********************************/
 
 Mat mu1, mu2;   // PRELIMINARY COMPUTING
 GaussianBlur(I1, mu1, Size(11, 11), 1.5);
 GaussianBlur(I2, mu2, Size(11, 11), 1.5);
 
 Mat mu1_2   =   mu1.mul(mu1);
 Mat mu2_2   =   mu2.mul(mu2);
 Mat mu1_mu2 =   mu1.mul(mu2);
 
 Mat sigma1_2, sigma2_2, sigma12;
 
 GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= mu1_2;
 
 GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= mu2_2;
 
 GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
    sigma12 -= mu1_mu2;
 
 Mat t1, t2, t3;
 
    t1 = 2 * mu1_mu2 + C1;
    t2 = 2 * sigma12 + C2;
    t3 = t1.mul(t2);              // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
    t1 = mu1_2 + mu2_2 + C1;
    t2 = sigma1_2 + sigma2_2 + C2;
    t1 = t1.mul(t2);               // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
 Mat ssim_map;
 divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;
 
 Scalar mssim = mean( ssim_map ); // mssim = average of ssim map
 return mssim;
}

Scalar getMSSIM_CUDA( const Mat& i1, const Mat& i2)
{
 const float C1 = 6.5025f, C2 = 58.5225f;
    /***************************** INITS **********************************/
 cuda::GpuMat gI1, gI2, gs1, tmp1,tmp2;
 
    gI1.upload(i1);
    gI2.upload(i2);
 
    gI1.convertTo(tmp1, CV_MAKE_TYPE(CV_32F, gI1.channels()));
    gI2.convertTo(tmp2, CV_MAKE_TYPE(CV_32F, gI2.channels()));
 
    vector<cuda::GpuMat> vI1, vI2;
    cuda::split(tmp1, vI1);
    cuda::split(tmp2, vI2);
 Scalar mssim;
 
 Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(vI2[0].type(), -1, Size(11, 11), 1.5);
 
 for( int i = 0; i < gI1.channels(); ++i )
    {
 cuda::GpuMat I2_2, I1_2, I1_I2;
 
        cuda::multiply(vI2[i], vI2[i], I2_2);        // I2^2
        cuda::multiply(vI1[i], vI1[i], I1_2);        // I1^2
        cuda::multiply(vI1[i], vI2[i], I1_I2);       // I1 * I2
 
        /*************************** END INITS **********************************/
 cuda::GpuMat mu1, mu2;   // PRELIMINARY COMPUTING
        gauss->apply(vI1[i], mu1);
        gauss->apply(vI2[i], mu2);
 
 cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
        cuda::multiply(mu1, mu1, mu1_2);
        cuda::multiply(mu2, mu2, mu2_2);
        cuda::multiply(mu1, mu2, mu1_mu2);
 
 cuda::GpuMat sigma1_2, sigma2_2, sigma12;
 
        gauss->apply(I1_2, sigma1_2);
        cuda::subtract(sigma1_2, mu1_2, sigma1_2); // sigma1_2 -= mu1_2;
 
        gauss->apply(I2_2, sigma2_2);
        cuda::subtract(sigma2_2, mu2_2, sigma2_2); // sigma2_2 -= mu2_2;
 
        gauss->apply(I1_I2, sigma12);
        cuda::subtract(sigma12, mu1_mu2, sigma12); // sigma12 -= mu1_mu2;
 
 cuda::GpuMat t1, t2, t3;
 
        mu1_mu2.convertTo(t1, -1, 2, C1); // t1 = 2 * mu1_mu2 + C1;
        sigma12.convertTo(t2, -1, 2, C2); // t2 = 2 * sigma12 + C2;
        cuda::multiply(t1, t2, t3);        // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
        cuda::addWeighted(mu1_2, 1.0, mu2_2, 1.0, C1, t1);       // t1 = mu1_2 + mu2_2 + C1;
        cuda::addWeighted(sigma1_2, 1.0, sigma2_2, 1.0, C2, t2); // t2 = sigma1_2 + sigma2_2 + C2;
        cuda::multiply(t1, t2, t1);                              // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
 cuda::GpuMat ssim_map;
        cuda::divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;
 
 Scalar s = cuda::sum(ssim_map);
        mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols);
 
    }
 return mssim;
}

struct BufferMSSIM                                     // Optimized CUDA versions
{   // Data allocations are very expensive on CUDA. Use a buffer to solve: allocate once reuse later.
 cuda::GpuMat gI1, gI2, gs, t1,t2;
 
 cuda::GpuMat I1_2, I2_2, I1_I2;
    vector<cuda::GpuMat> vI1, vI2;
 
 cuda::GpuMat mu1, mu2;
 cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
 
 cuda::GpuMat sigma1_2, sigma2_2, sigma12;
 cuda::GpuMat t3;
 
 cuda::GpuMat ssim_map;
 
 cuda::GpuMat buf;
};

Scalar getMSSIM_CUDA_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b)
{
 const float C1 = 6.5025f, C2 = 58.5225f;
    /***************************** INITS **********************************/
 
    b.gI1.upload(i1);
    b.gI2.upload(i2);
 
 cuda::Stream stream;
 
    b.gI1.convertTo(b.t1, CV_32F, stream);
    b.gI2.convertTo(b.t2, CV_32F, stream);
 
    cuda::split(b.t1, b.vI1, stream);
    cuda::split(b.t2, b.vI2, stream);
 Scalar mssim;
 
 Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(b.vI1[0].type(), -1, Size(11, 11), 1.5);
 
 for( int i = 0; i < b.gI1.channels(); ++i )
    {
        cuda::multiply(b.vI2[i], b.vI2[i], b.I2_2, 1, -1, stream);        // I2^2
        cuda::multiply(b.vI1[i], b.vI1[i], b.I1_2, 1, -1, stream);        // I1^2
        cuda::multiply(b.vI1[i], b.vI2[i], b.I1_I2, 1, -1, stream);       // I1 * I2
 
        gauss->apply(b.vI1[i], b.mu1, stream);
        gauss->apply(b.vI2[i], b.mu2, stream);
 
        cuda::multiply(b.mu1, b.mu1, b.mu1_2, 1, -1, stream);
        cuda::multiply(b.mu2, b.mu2, b.mu2_2, 1, -1, stream);
        cuda::multiply(b.mu1, b.mu2, b.mu1_mu2, 1, -1, stream);
 
        gauss->apply(b.I1_2, b.sigma1_2, stream);
        cuda::subtract(b.sigma1_2, b.mu1_2, b.sigma1_2, cuda::GpuMat(), -1, stream);
 //b.sigma1_2 -= b.mu1_2;  - This would result in an extra data transfer operation
 
        gauss->apply(b.I2_2, b.sigma2_2, stream);
        cuda::subtract(b.sigma2_2, b.mu2_2, b.sigma2_2, cuda::GpuMat(), -1, stream);
 //b.sigma2_2 -= b.mu2_2;
 
        gauss->apply(b.I1_I2, b.sigma12, stream);
        cuda::subtract(b.sigma12, b.mu1_mu2, b.sigma12, cuda::GpuMat(), -1, stream);
 //b.sigma12 -= b.mu1_mu2;
 
 //here too it would be an extra data transfer due to call of operator*(Scalar, Mat)
        cuda::multiply(b.mu1_mu2, 2, b.t1, 1, -1, stream); //b.t1 = 2 * b.mu1_mu2 + C1;
        cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
        cuda::multiply(b.sigma12, 2, b.t2, 1, -1, stream); //b.t2 = 2 * b.sigma12 + C2;
        cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -12, stream);
 
        cuda::multiply(b.t1, b.t2, b.t3, 1, -1, stream);     // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
        cuda::add(b.mu1_2, b.mu2_2, b.t1, cuda::GpuMat(), -1, stream);
        cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
 
        cuda::add(b.sigma1_2, b.sigma2_2, b.t2, cuda::GpuMat(), -1, stream);
        cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -1, stream);
 
 
        cuda::multiply(b.t1, b.t2, b.t1, 1, -1, stream);     // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
        cuda::divide(b.t3, b.t1, b.ssim_map, 1, -1, stream);      // ssim_map =  t3./t1;
 
        stream.waitForCompletion();
 
 Scalar s = cuda::sum(b.ssim_map, b.buf);
        mssim.val[i] = s.val[0] / (b.ssim_map.rows * b.ssim_map.cols);
 
    }
 return mssim;
}

どうやって行うか? - GPU

上で見たように、各処理に対して3種類の関数がある。1つはCPU用、2つはGPU用である。GPU用を2つ作った理由は、CPUからGPUへ単純に移植するだけでは実際にはかえって遅くなることが多い、という点を示すためである。性能向上を得たいなら、いくつかの規則を覚えておく必要があり、それらについては後ほど詳しく説明する。

GPUモジュールは、CPU側の対応物にできるだけ似るように開発された。これにより移植作業が容易になる。コードを書く前にまずやるべきことは、GPUモジュールをプロジェクトにリンクし、そのモジュールのヘッダファイルをインクルードすることである。GPUのすべての関数とデータ構造は、cv名前空間のgpuサブ名前空間に含まれている。use namespaceキーワードでデフォルト名前空間に加えてもよいし、混乱を避けるためにcv::を付けてどこでも明示的に書いてもよい。ここでは後者を採用する。

#include <opencv2/gpu.hpp> // GPU structures and methods

GPUは「graphics processing unit（グラフィックス処理装置）」の略である。もともとはグラフィカルなシーンを描画するために作られた。こうしたシーンは何らかの形で大量のデータの上に成り立っている。とはいえ、それらはすべてが逐次的に互いに依存しているわけではなく、並列処理が可能である。このためGPUは複数の小さな処理ユニットを内蔵している。これらは最先端のプロセッサではなく、CPUとの1対1の比較では後れを取る。しかし、その強みは数の多さにある。近年、GPUのこの大規模な並列処理能力を非グラフィカルなシーンでも活用しようという傾向が高まっている。描画も同様である。これにより、GPGPU（general-purpose computation on graphics processing units、GPUによる汎用計算）が生まれた。

GPUは独自のメモリを持つ。OpenCVでハードディスクからデータをMatオブジェクトに読み込むと、それはシステムメモリ上に置かれる。CPUは（キャッシュを介して）これに何らかの形で直接作用するが、GPUはそうできない。GPUは計算に必要な情報をシステムメモリから自身のメモリへ転送しなければならない。これはアップロード処理によって行われ、時間がかかる。最終的に、結果はCPUが見て使えるようにシステムメモリへダウンロードし直す必要がある。小さな関数をGPUに移植することは推奨されない。アップロード/ダウンロード時間が、並列実行で得られる時間より大きくなるからである。

Matオブジェクトはシステムメモリ（またはCPUキャッシュ）にのみ格納される。OpenCVの行列をGPUへ持って行くには、そのGPU側の対応物であるcv::cuda::GpuMatを使う必要がある。これはMatと同様に動作するが、2Dに限定され、関数が参照を返さない（GPU参照とCPU参照を混在させられない）という制約がある。MatオブジェクトをGPUへアップロードするには、クラスのインスタンスを作成した後にupload関数を呼ぶ。ダウンロードには、Matオブジェクトへの単純な代入を使ってもよいし、download関数を使ってもよい。

Mat I1;         // Main memory item - read image into with imread for example
gpu::GpuMat gI; // GPU matrix - for now empty
gI1.upload(I1); // Upload a data from the system memory to the GPU memory
 
I1 = gI1;       // Download, gI1.download(I1) will work too

データがGPUメモリ上に揃ったら、OpenCVのGPU対応関数を呼べる。ほとんどの関数はCPU上と同じ名前のままで、GpuMat入力のみを受け付ける点だけが異なる。

もう1つ覚えておくべきことは、すべてのチャンネル数で効率的なアルゴリズムをGPU上で作れるわけではない、という点である。一般に、GPU画像の入力画像はチャンネルが1または4のいずれかで、要素サイズはcharまたはfloat型のいずれかである必要があることがわかった。GPUではdoubleはサポートされていない、残念ながら。一部の関数に他の型のオブジェクトを渡すと例外がスローされ、エラー出力にエラーメッセージが出る。ドキュメントはほとんどの箇所で入力に受け付けられる型を詳述している。入力として3チャンネル画像がある場合は、2つの選択肢がある。新しいチャンネルを追加する（そしてchar要素を使う）か、画像を分割して各画像に対して関数を呼ぶかである。前者はメモリを無駄にするため、あまり推奨されない。

要素の位置（隣接要素）が問題にならない一部の関数では、手っ取り早い解決策はシングルチャンネル画像へ整形し直すことである。これはPSNRの実装の場合に当てはまり、そこではabsdiffメソッドにおいて隣接要素の値は重要ではない。しかしGaussianBlurではこれは選択肢にならず、SSIMではsplitメソッドを使う必要がある。この知識があれば、GPUで動くコード（私のGPU版のような）を作って実行できる。それがCPU実装より遅くなりうることを知って驚くだろう。

最適化

この理由は、メモリ確保とデータ転送のコストを無駄にしているからである。そしてGPUではこれが極めて高い。最適化のもう1つの可能性は、cv::cuda::Streamの助けを借りて非同期のOpenCV GPU呼び出しも導入することである。

GPU上でのメモリ確保はかなりのコストがかかる。したがって、可能であれば新規メモリの確保はできるだけ少ない回数にする。複数回呼び出すつもりの関数を作るなら、その関数のローカルな引数を最初の呼び出し時に一度だけ確保するのがよい。そのためには、使用するすべてのローカル変数を含むデータ構造を作る。例えばPSNRの場合は次のようになる:
struct BufferPSNR // Optimized GPU versions

{ // Data allocations are very expensive on GPU. Use a buffer to solve: allocate once reuse later.

gpu::GpuMat gI1, gI2, gs, t1,t2;

gpu::GpuMat buf;

};

次に、メインプログラムでこれのインスタンスを作る:
BufferPSNR bufferPSNR;

そして最後に、関数を呼ぶたびにこれを渡す:
double getPSNR_GPU_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)

これで、これらのローカルな引数にはb.gI1、b.bufなどとしてアクセスできる。GpuMatは、新しい行列サイズが前回と異なる場合にのみ、新たな呼び出しで自身を再確保する。
不要な関数間のデータ転送を避ける。GPUに移ると、どんな小さなデータ転送でも無視できなくなる。したがって、可能であればすべての計算をその場（in-place）で行う（言い換えれば新しいメモリオブジェクトを作らない。理由は前の項目で説明したとおりである）。例えば、算術演算は1行の式で表す方が書きやすいかもしれないが、それは遅くなる。SSIMの場合、ある時点で次を計算する必要がある:
b.t1 = 2 * b.mu1_mu2 + C1;

上記の呼び出しは成功するが、そこに隠れたデータ転送が存在することに注意せよ。加算を行う前に、乗算の結果をどこかに格納する必要がある。したがって、裏でローカルな行列を作り、それにC1の値を加え、最後にそれをt1へ代入する。これを避けるために、算術演算子の代わりにgpu関数を使う:
gpu::multiply(b.mu1_mu2, 2, b.t1); //b.t1 = 2 * b.mu1_mu2 + C1;

gpu::add(b.t1, C1, b.t1);
非同期呼び出しを使用する（cv::cuda::Stream）。デフォルトでは、GPU関数を呼び出すたびに、その呼び出しが完了するのを待ってから結果とともに戻る。しかし、非同期呼び出しを行うことも可能であり、これは演算の実行を呼び出してアルゴリズムのためのコストの高いデータ割り当てを行い、すぐに戻ってくることを意味する。これで、必要であれば別の関数を呼び出せる。MSSIMにとって、これは小さな最適化ポイントである。デフォルト実装では、画像をチャンネルに分割し、各チャンネルに対してGPU関数を呼び出す。ストリームを使うとわずかながら並列化が可能である。ストリームを使うことで、GPUが既に与えられたメソッドを実行している間に、データ割り当てとアップロード操作を行える。例えば、2枚の画像をアップロードする必要があるとする。これらを次々にキューに入れ、それを処理する関数を呼び出す。関数はアップロードの完了を待つが、その間に次に実行される関数のための出力バッファの割り当てを行う。
gpu::Stream stream;

stream.enqueueConvert(b.gI1, b.t1, CV_32F); // Upload

gpu::split(b.t1, b.vI1, stream); // Methods (pass the stream as final parameter).

gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2

CV_32F
#define CV_32F
Definition interface.h:59

結果と結論

Intel P8700ノートPC用CPUとローエンドのNVIDIA GT220Mの組み合わせでは、性能数値は次のとおりである:

Time of PSNR CPU (averaged for 10 runs): 41.4122 milliseconds. With result of: 19.2506
Time of PSNR GPU (averaged for 10 runs): 158.977 milliseconds. With result of: 19.2506
Initial call GPU optimized:              31.3418 milliseconds. With result of: 19.2506
Time of PSNR GPU OPTIMIZED ( / 10 runs): 24.8171 milliseconds. With result of: 19.2506
 
Time of MSSIM CPU (averaged for 10 runs): 484.343 milliseconds. With result of B0.890964 G0.903845 R0.936934
Time of MSSIM GPU (averaged for 10 runs): 745.105 milliseconds. With result of B0.89922 G0.909051 R0.968223
Time of MSSIM GPU Initial Call            357.746 milliseconds. With result of B0.890964 G0.903845 R0.936934
Time of MSSIM GPU OPTIMIZED ( / 10 runs): 203.091 milliseconds. With result of B0.890964 G0.903845 R0.936934

どちらの場合も、CPU実装と比べてほぼ100%の性能向上を達成できた。これがあなたのアプリケーションを動作させるのに必要なちょうどの改善かもしれない。これの実行中の様子はこちらのYouTubeで見られる。

目次

目的

ソースコード

どうやって行うか? - GPU

最適化

結果と結論