あかすくぱるふぇ

同人サークル「あかすくぱるふぇ」のブログです。

CUDA Sampleでよく行われている処理。
ただ、CUDA Sampleはごちゃごちゃしていてよくわからないので、最小構成を抜き出してみました。
int main(int argc, char *argv[])
{
// 画像ファイル入力
cvImg = cv::imread("Lenna.bmp", 0);

// 初期化
init(argc, argv);

glutMainLoop();

// 終了処理
finalize();

return 0;
}
main関数。
画像ファイルの入力と初期化を行った後に、glutMainLoop()を実行します。

void init(int argc, char *argv[])
{
size_t dataSize = cvImg.cols * cvImg.rows * sizeof(unsigned char);

glInit(argc, argv); // GLの初期化
GenTexture(); // テクスチャの生成
GenAndRegisterPBO(dataSize); // PBOの生成とCUDAへの登録
   
// 入力画像のGPU側メモリ確保
  checkCudaErrors(cudaMalloc((unsigned char**)&dImgIn, dataSize));
}
初期化関数です。
以下に各関数を載せます。

void glInit(int argc, char *argv[])
{
glutInit(&argc, argv);
glutInitWindowSize(cvImg.cols, cvImg.rows);
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
glutCreateWindow(argv[0]);
glewInit();
glutDisplayFunc(display);
}

void GenTexture()
{
glGenTextures(1, &tex);
glBindTexture(GL_TEXTURE_2D, tex);
glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, cvImg.cols, cvImg.rows,
      0, GL_LUMINANCE, GL_UNSIGNED_BYTE, NULL);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
}

void GenAndRegisterPBO(const size_t dataSize)
{
glGenBuffers(1, &pbo);
glBindBuffer(GL_PIXEL_UNPACK_BUFFER, pbo);
glBufferData(GL_PIXEL_UNPACK_BUFFER, dataSize, 0, GL_STREAM_DRAW);
glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);
cudaGraphicsGLRegisterBuffer(&cuda_pbo_resource, pbo,
      cudaGraphicsMapFlagsWriteDiscard);
}
glInit()はGLUT, GLEWのおなじみの処理。
GenTexture()はテクスチャ生成のおなじみの処理。
GenAndRegisterPBO()はPBOを生成し、そのPBOをCUDAに登録しています。
CUDAへの登録はcudaGraphicsGLRegisterBuffer()で行います。
cudaGraphicsMapFlagsWriteDiscardはPBOがCUDAからWriteOnlyであるというフラグです。

void display()
{
// カーネル関数で書き換える画像(PBO)のポインタを取得
unsigned char *dImgOut;
size_t dataSizeTmp;
checkCudaErrors(cudaGraphicsMapResources(1, &cuda_pbo_resource, 0));
checkCudaErrors(cudaGraphicsResourceGetMappedPointer(
(void **)&dImgOut, &dataSizeTmp, cuda_pbo_resource));

// 入力画像をGPUのグローバルメモリに転送
size_t dataSize = cvImg.cols * cvImg.rows * sizeof(unsigned char);
checkCudaErrors(cudaMemcpy(dImgIn, cvImg.data, dataSize, cudaMemcpyHostToDevice));

// カーネル関数の実行
dim3 grid(cvImg.rows);
dim3 block(1024);
invertKernel << <grid, block >> >(dImgOut, dImgIn, cvImg.cols);

checkCudaErrors(cudaGraphicsUnmapResources(1, &cuda_pbo_resource, 0));

glBindTexture(GL_TEXTURE_2D, tex);
glBindBuffer(GL_PIXEL_UNPACK_BUFFER, pbo);
glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, cvImg.cols, cvImg.rows,
GL_LUMINANCE, GL_UNSIGNED_BYTE, 0);
glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);

glEnable(GL_TEXTURE_2D);
glBegin(GL_QUADS);
{
glTexCoord2f(0, 1);
glVertex2f(-1, -1);
glTexCoord2f(1, 1);
glVertex2f(1, -1);
glTexCoord2f(1, 0);
glVertex2f(1, 1);
glTexCoord2f(0, 0);
glVertex2f(-1, 1);
}
glEnd();

glBindTexture(GL_TEXTURE_2D, 0);
glutSwapBuffers();
}
display関数です。
PBOのポインタ取得、入力画像のGPU転送、カーネル実行、テクスチャマッピングを行っています。
処理結果をCPUに戻さずに、PBOからテクスチャにコピー(GPU内で閉じてる)して描画しているのがポイントです。

念のため、カーネル関数の一例も載せておきます。
__global__ void invertKernel(unsigned char *out, const unsigned char *in, const int width)
{
for (int i = threadIdx.x; i < width; i += blockDim.x) {
int index = blockIdx.x * width + i;
out[index] = 255 - in[index];
}
}

以上です。

・分析、設計、実装、テストといった開発工程がどれも、途切れることなく続く連続的な取り組みになる。
・イテレーションは全力疾走する短距離走で、この期間を通じて、ユーザーストーリーをリリース可能な動くソフトウェアに変換する。
・リリースはMinimalかつMarketableであるべき。システムの価値の80%が機能の20%からもたらされている。
・コードはいつもばっちり結合した状態にしておく。最後まで全体として動作するかわからない状況に甘んじてはいけない。
・メンバーの役割は気にしない。気にかけるのは、その役割によって果たされる仕事が、チームの誰かによってちゃんとこなされるようになっているかどうかだけ。
・顧客に積極的にかかわってもらう。
・「自分がいる理由」、「目的」、「顧客」、「主要課題」、「判断者」は少なくとも明確にする。
・見積もりを行うのは、ターゲットが達成可能な程度に現実的なものかどうか判断するためである。
・見積もりは理想日単位ではなく、相対的なポイントで表現した方がよい。
・リファクタリングは常に細かくやり続ける。
・カバレッジ100%を目指すべきではない。網羅性ではなく、リリース準備が整っていると確信を持つことが大事。
・ビルド時間は10分以内が望ましい。
・参考書:レガシーコード改善ガイド、テスト駆動開発入門

以下の記事のコードに対して、ループ展開とシェアードメモリによる高速化を試みました。
【CUDA】グローバルメモリとテクスチャメモリの速度比較

以下の修士論文を参考にさせていただきました。
GPUを用いた画像処理の高速化

コードは以下の通り。
const int blockSize = 32;
const int blockSizeWithoutEdge = blockSize - 2;
const int blockSize1d = blockSize * blockSize;

// 平均値フィルタ(GlobalMemory)
__global__ void meanFilterKernel(unsigned char *out, const unsigned char *in,
const int width, const int height, const int filterRadius, const float fscale)
{
if (blockIdx.x - filterRadius < 0 || blockIdx.x + filterRadius >= height) return;

for (int i = threadIdx.x; i < width; i += blockDim.x) {
int sum = 0;
for (int fy = -filterRadius; fy <= filterRadius; ++fy) {
int y = blockIdx.x + fy;
for (int fx = -filterRadius; fx <= filterRadius; ++fx) {
sum += in[y * width + i + fx];
}
}
out[blockIdx.x * width + i] = sum * fscale;
}
}

// 平均値フィルタ(GlobalMemory, Unwind)
__global__ void meanFilterKernelUnwind(unsigned char *out, const unsigned char *in,
const int width, const int height)
{
if (blockIdx.x == 0 || blockIdx.x == height - 1) return;

const float scale = 0.1111111111;
for (int i = threadIdx.x; i < width; i += blockDim.x) {
int sum = 0;
const unsigned char* tmp = in + (blockIdx.x - 1) * width + i - 1;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
tmp += width - 2;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
tmp += width - 2;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
out[blockIdx.x * width + i] = sum * scale;
}
}

// 平均値フィルタ(SharedMemory, Unwind)
__constant__ int dTopLeftIndex[8000];
__global__ void meanFilterKernelUnwindTileShared(unsigned char *out,
const unsigned char *in, const int width)
{
__shared__ unsigned char tile[blockSize1d];
__shared__ int topLeftIndex;
int indexInTile = threadIdx.y * blockDim.x + threadIdx.x;

if (threadIdx.x == 0 && threadIdx.y == 0) topLeftIndex = dTopLeftIndex[blockIdx.x];
__syncthreads();

int indexInImage = topLeftIndex + threadIdx.y * width + threadIdx.x;
tile[indexInTile] = in[indexInImage];
__syncthreads();

if (threadIdx.x == 0 || threadIdx.x == blockDim.x - 1 ||
threadIdx.y == 0 || threadIdx.y == blockDim.y - 1) return;

const float scale = 0.1111111111;
int sum = 0;
const unsigned char* tmp = tile + indexInTile - blockDim.x - 1;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
tmp += blockDim.x - 2;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
tmp += blockDim.x - 2;
sum += *tmp++;
sum += *tmp++;
sum += *tmp;
out[indexInImage] = sum * scale;
}

void main()
{
cv::Mat img = cv::imread("Room_02000.bmp", 0);
int pixelNum = img.cols * img.rows;
int dataSize = pixelNum * sizeof(unsigned char);

unsigned char *hOut, *dIn, *dOut;
hOut = new unsigned char[pixelNum];
checkCudaErrors(cudaMalloc((unsigned char**)&dIn, dataSize));
checkCudaErrors(cudaMalloc((unsigned char**)&dOut, dataSize));
checkCudaErrors(cudaMemcpy(dIn, img.data, dataSize, cudaMemcpyHostToDevice));

cudaArray *cu_array;
cudaChannelFormatDesc desc = cudaCreateChannelDesc<unsigned char>();
checkCudaErrors(cudaMallocArray(&cu_array, &desc, img.cols, img.rows));
checkCudaErrors(cudaMemcpyToArray(cu_array, 0, 0, img.data, dataSize,
cudaMemcpyHostToDevice));
checkCudaErrors(cudaBindTextureToArray(tex, cu_array));

LARGE_INTEGER freq;
QueryPerformanceFrequency(&freq);

int filterRadius = 1;
int filterDiameter = filterRadius * 2 + 1;
float fscale = 1.0f / filterDiameter / filterDiameter;

// CPU
LARGE_INTEGER start, end;
QueryPerformanceCounter(&start);
for (int h = filterRadius; h < img.rows - filterRadius; ++h) {
for (int w = filterRadius; w < img.cols - filterRadius; ++w) {
int sum = 0;
for (int fy = -filterRadius; fy <= filterRadius; ++fy) {
int y = h + fy;
for (int x = -filterRadius; x <= filterRadius; ++x) {
sum += img.data[y * img.cols + w + x];
}
}
hOut[h * img.cols + w] = sum * fscale;
}
}
QueryPerformanceCounter(&end);
printf("CPU:%lf[ms]\n", (double)(end.QuadPart - start.QuadPart)
/ freq.QuadPart * 1000.0);
//memcpy(img.data, hOut, dataSize);

// GPU(GlobalMemory)
dim3 block(1024);
dim3 grid(img.rows);
int n = 100;
QueryPerformanceCounter(&start);
for (int i = 0; i < n; ++i) {
meanFilterKernel << <grid, block >> >(dOut, dIn, img.cols, img.rows,
filterRadius, fscale);
checkCudaErrors(cudaDeviceSynchronize());
}
QueryPerformanceCounter(&end);
printf("GPU(GlobalMemory):%lf[ms]\n", (double)(end.QuadPart - start.QuadPart)
/ freq.QuadPart * 1000.0 / n);
//checkCudaErrors(cudaMemcpy(img.data, dOut, dataSize, cudaMemcpyDeviceToHost));

// GPU(GlobalMemory, Unwind)
QueryPerformanceCounter(&start);
for (int i = 0; i < n; ++i) {
meanFilterKernelUnwind << <grid, block >> >(dOut, dIn, img.cols, img.rows);
checkCudaErrors(cudaDeviceSynchronize());
}
QueryPerformanceCounter(&end);
printf("GPU(GlobalMemory, Unwind):%lf[ms]\n", (double)(end.QuadPart - start.QuadPart)
/ freq.QuadPart * 1000.0 / n);
//checkCudaErrors(cudaMemcpy(img.data, dOut, dataSize, cudaMemcpyDeviceToHost));

//
// GPU(SharedMemory, Unwind)
//

// ブロックの左上端の座標を記録
int blockNumX = img.cols / blockSizeWithoutEdge;
int blockNumY = img.rows / blockSizeWithoutEdge;

vector<int> topLeftIndex, topLeftX, topLeftY;
for (int by = 0; by < blockNumY; ++by) {
int y = by * blockSizeWithoutEdge;
for (int bx = 0; bx < blockNumX; ++bx) {
int x = bx * blockSizeWithoutEdge;
topLeftIndex.push_back(y * img.cols + x);
}
}

// 右端ブロック
if (img.cols % blockSizeWithoutEdge != 0) {
int x = img.cols - blockSize;
for (int by = 0; by < blockNumY; ++by) {
int y = by * blockSizeWithoutEdge;
topLeftIndex.push_back(y * img.cols + x);
}
}
// 下端ブロック
if (img.rows % blockSizeWithoutEdge != 0) {
int y = img.rows - blockSize;
for (int bx = 0; bx < blockNumX; ++bx) {
int x = bx * blockSizeWithoutEdge;
topLeftIndex.push_back(y * img.cols + x);
}
}
// 右下端ブロック
if (img.cols % blockSizeWithoutEdge != 0 && img.rows % blockSizeWithoutEdge != 0) {
int x = img.cols - blockSize;
int y = img.rows - blockSize;
topLeftIndex.push_back(y * img.cols + x);
}

checkCudaErrors(cudaMemcpyToSymbol(dTopLeftIndex, &topLeftIndex[0], sizeof(int)
* topLeftIndex.size()));

dim3 blockTile(blockSize, blockSize);
dim3 gridTile(topLeftIndex.size());
QueryPerformanceCounter(&start);
for (int i = 0; i < n; ++i) {
meanFilterKernelUnwindTileShared << <gridTile, blockTile >> >(dOut, dIn, img.cols);
checkCudaErrors(cudaDeviceSynchronize());
}
QueryPerformanceCounter(&end);
printf("GPU(SharedMemory, Unwind):%lf[ms]\n", (double)(end.QuadPart - start.QuadPart)
/ freq.QuadPart * 1000.0 / n);
//checkCudaErrors(cudaMemcpy(img.data, dOut, dataSize, cudaMemcpyDeviceToHost));
}
結果は以下の通り。
無題

ループ展開によって処理時間が約2/3になりましたが、シェアードメモリを利用しても
それ以上は速くなりませんでした。

以上です。

※追記
フィルタサイズを大きくすると、シェアードメモリが相対的に速くなりました。
3x3の場合は元々グローバルメモリへのアクセス回数が少ないので、シェアードメモリの効果が見えにくかったということですね。

↑このページのトップヘ