はじめに

画像処理プログラムを書くにあたって，並列処理ができると上限はあるものの，基本的に高速化できます。そんな並列処理を主に担うプロセッサーの代表格的存在としてGPUがあります。この記事では，そんなGPUを駆使した画像処理プログラムの一例を書いていきたいと思います。なお，OpenCVやCUDAに関する基礎知識はある前提で書きますので，ご了承ください。
また，cv::cuda::GpuMatなるクラスも存在するのですが，あまり汎用性が高くないと判断したので今回は用いらないこととします。

並列化前

まず，今回書きたい並列画像処理と同様の結果を得られる可読性重視のコードを載せておきます。つまり，シンプルにfor文の2重ループで実装した場合です。

source.cpp

#include <opencv2/opencv.hpp>

using namespace cv;

void processing(uchar* src, uchar* dst, size_t src_pitch, size_t dst_pitch, int w, int h, int ch) {

	for (int x = 0; x < w; x++) {
		for (int y = 0; y < h; y++) {

			uchar* src_row = src + y * src_pitch;
			uchar* src_pix = src_row + x * ch;

			uchar* dst_row = dst + y * dst_pitch;
			uchar* dst_pix = dst_row + x * ch;

			for (int c = 0; c < ch; c++) {

				*(dst_pix + c) = *(src_pix + (ch - 1 - c));

			}

		}
	}

}

int main(void) {

	Mat input = imread("image.jpg");
	Mat output = Mat(input.rows, input.cols, CV_8UC3);

	processing(input.data, output.data, input.step, output.step, input.cols, input.rows, input.channels());

	imshow("output", output);
	waitKey(0);

    input.release();
    output.release();

	return 0;

}

コードを見てもらえばお分かりいただけるかと思いますが，まあつまりは各ピクセルにおけるR成分とB成分を入れ替えて，色を反転させるようなエフェクトになります。
基礎と言えば基礎なのか，判断が難しいところですが，画像処理というよりもメモリ管理に関するお話なのでpitchってなに？ という話とかは後に追いかけながらすることにします。

CUDAの世界へ

CUDAというかNVIDIA製のグラフィックボードで各種演算を行うためには，まずグラフィックボード側でメモリを確保する必要がありますね。ちなみに，CUDAにおいてはCPU側のことをホスト(Host)，グラフィックボード側のことをデバイス(Device)と呼びますので，これ以降ではそのように名前を統一して用いることにします。

メモリの確保について

デバイス側のメモリを確保する手段はいくつかありますが，今回画像処理ということで2次元的なメモリの扱い方をするのに最適な関数の一つとして，cudaMallocPitchを用いることにします。というのも，OpenCVのMatクラスにおけるピクセルデータをメモリに保管する扱い方として，行と行の間のメモリを余分に確保されている場合があります。これを考慮した各行の長さをピッチ(pitch) と呼びます。故に，余分に確保されたメモリを無視して，例えば先ほどの2重forループ内で次のように画像処理プログラムを書いてしまうとメモリサイズが合わず，画像全体に対して処理できない可能性があります。

example.cpp

...

for (int c = 0; c < ch; c++) {

    dst[ c + (x + y * w) * ch ] = src[ c + (x + y * w) * ch ];

}

...

そこでピッチを考慮します。OpenCVのMatクラスではすべての行において，ピッチが等しく確保されるようになっており，それはMatクラスのstep というメンバ変数に保管されています。これを用いてピッチを考慮したのが最初に挙げたコードでした。
しかるに，cudaMallocPitch関数ではピッチを考慮するようにデバイス側でメモリを確保してくれるのでMatクラスと相性が良いのです。ということで，まずメモリを確保しましょう。

source.cu

...

Mat input = imread("image.jpg");
Mat output = Mat(input.rows, input.cols, CV_8UC3);

int w = input.cols;
int h = input.rows;
int c = input.channels();

uchar* cuda_src = nullptr;
uchar* cuda_dst = nullptr;

size_t cuda_src_pitch = 0;
size_t cuda_dst_pitch = 0;
cudaMallocPitch(&cuda_src, &cuda_src_pitch, w * c * sizeof(uchar), h);
cudaMallocPitch(&cuda_dst, &cuda_dst_pitch, w * c * sizeof(uchar), h);

...

cudaMallocPitch関数の第二引数に渡したポインタには，確保したメモリのピッチが入ります。

デバイスで確保したメモリへ画像データを転送

ピッチを考慮したデータ転送に適する関数として，cudaMemcpy2D関数があります。この関数は，第二引数に転送先のピッチ，また第4引数に転送元のピッチを指定することができます。

source.cu

...

cudaMemcpy2D(cuda_src, cuda_src_pitch, input.data, input.step, w * c * sizeof(uchar), h, cudaMemcpyHostToDevice);

...

カーネル関数の実行に関するブロックとグリッドの設定

デバイス側での演算を行うカーネル関数を実行するにあたり，設定しなければならないスレッドやブロック，グリッドですが，まず前提としてCUDAでは，全グリッド・ブロック共通で一度に動かすことのできるスレッド数に最大値が存在します。それを取得するのは，cudaGetDeviceProperties関数よりデバイスのプロパティを取得することにより可能です。

getProp.cu

...

cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, 0); // 第二引数でプロパティを取得するデバイスを指定

// 最大スレッド数をコンソールに表示
printf("max threads per block:%d\n", prop.maxThreadsPerBlock);

...

これを考慮したうえで，さらにシェアードメモリの存在を考えたいところです。シェアードメモリとは，簡単に言うと同じブロック内のスレッド同士ではメモリのやり取りが高速に行うことができる，その存在のことだと思ってもらえればよいかと思います。（正直あまり詳しくない、、）
そこで，このシェアードメモリも考慮するとなるとき，1ブロックを正方形になるように確保すると考えやすいと思います。（本当に高速なのかはわからないが、、）
また，調べてみたところmaxThreadPerBlockは基本的に2の累乗になるように設計されているらしく，ビット演算で扱いやすいです。つまり，より高速に計算でき得る。
そこでまず，maxThreadPerBlockが2の何乗なのか，その指数を計算する関数を定義しておきます。

source.cu

...

int log2_bitwise(unsigned int n) {

    int exponent = 0;
    while (n >>= 1) exponent++;
    return exponent;

}

...

ところで，maxThreadPerBlockについて，自分の環境では1024でした。例えばこの値について，次のように2乗の形を作ることができます。

\begin{split}
    1024 &= 2^{10} \\
         &= (2^5)^2 \\
\end{split}

つまり，1ブロックを正方形に考えた根底はここにもあって，横・縦軸方向ともに$\sqrt ($maxThreadPerBlock$)$で計算でき，尚且つ環境によってはこれをビットシフトのみで計算することができます。
よって，これらを元にブロック数・グリッド数を計算するコードは以下のように記述できます。

source.cu

...

int exponent = log2_bitwise(prop.maxThreadPerBlock);
dim3 block(
    1 << (exponent >> 1),
    1 << (exponent >> 1),
    1
);
dim3 grid(
    (w / block.x) + (w ^ block.x),
    (h / block.x) + (h ^ block.y),
    1
);

...

カーネル関数の実装

今指定したスレッド・ブロック・グリッドの設定を元に，カーネル関数内で各ピクセルに対応する座標は次のように計算できます。

source.cu

...

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;

...

よって，最初に考えたコードと組み合わせることで，次のようにカーネル関数が記述できます。

source.cu

...

__global__ void kernel(
	uchar* src, uchar* dst,
	unsigned int src_pitch, unsigned int dst_pitch,
	int w, int h, int ch
) {

	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	if (x >= w || y >= h) return;

	uchar* src_row = src + y * src_pitch;
	uchar* src_pix = src_row + x * ch;

	uchar* dst_row = dst + y * dst_pitch;
	uchar* dst_pix = dst_row + x * ch;

	for (int c = 0; c < ch; c++) {
	
		*(dst_pix + c) = *(src_pix + (ch - 1 - c));
	
	}

}

...

あとは，デバイス側からホスト側のメモリ(Matクラスのdata)に計算結果を転送することで完了です。その際に，cudaMemcpy2D関数に渡す引数が先に書いたコードと逆になるので注意しましょう。

完成した全体のコード

上で書いたことをすべてまとめると次のコードが完成します。

source.cu

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <opencv2/opencv.hpp>

using namespace cv;

__global__ void kernel(
	uchar* src, uchar* dst,
	unsigned int src_pitch, unsigned int dst_pitch,
	int w, int h, int ch
) {

	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;

	if (x >= w || y >= h) return;

	uchar* src_row = src + y * src_pitch;
	uchar* src_pix = src_row + x * ch;

	uchar* dst_row = dst + y * dst_pitch;
	uchar* dst_pix = dst_row + x * ch;

	for (int c = 0; c < ch; c++) {
	
		*(dst_pix + c) = *(src_pix + (ch - 1 - c));
	
	}

}


int log2_bitwise(unsigned int n) {

    int exponent = 0;
    while (n >>= 1) exponent++;
    return exponent;

}

int main(void) {

    Mat input = imread("image.jpg");
	Mat output = Mat(input.rows, input.cols, CV_8UC3);
	
	int w = input.cols;
	int h = input.rows;
	int c = input.channels();

	uchar* cuda_src = nullptr;
	uchar* cuda_dst = nullptr;

	size_t cuda_src_pitch = 0;
	size_t cuda_dst_pitch = 0;
	cudaMallocPitch(&cuda_src, &cuda_src_pitch, w * c * sizeof(uchar), h);
	cudaMallocPitch(&cuda_dst, &cuda_dst_pitch, w * c * sizeof(uchar), h);
	
	cudaMemcpy2D(cuda_src, cuda_src_pitch, input.data, input.step, w * c * sizeof(uchar), h, cudaMemcpyHostToDevice);
	
	cudaDeviceProp prop;
	cudaGetDeviceProperties(&prop, 0);
	
    int exponent = log2_bitwise(prop.maxThreadPerBlock);
    dim3 block(
        1 << (exponent >> 1),
        1 << (exponent >> 1),
        1
    );
    dim3 grid(
        (w / block.x) + (w ^ block.x),
        (h / block.x) + (h ^ block.y),
        1
    );

	kernel <<<grid, block>>> (cuda_src, cuda_dst, cuda_src_pitch, cuda_dst_pitch, w, h, c);
	
	cudaMemcpy2D(output.data, output.step, cuda_dst, cuda_dst_pitch, w * c * sizeof(uchar), h, cudaMemcpyDeviceToHost);

	cudaFree(cuda_src);
	cudaFree(cuda_dst);
	
	imshow("output", output);
	waitKey(0);

	return 0;

}

最後に

途中，ブロック数の計算の段階で，ある程度環境に依存するような高速化を図りました。今回書いたコードを参考にされる場合そのことを念頭にしていただけると良いかと思います。
何か書きミスやアドバイスなどありましたらお気軽にコメントいただければと思います。おつかれさまでした！