NPP : Canny変換(そのに)

Canny変換は明るさの変化点を見つけることで輪郭を検出します。 カラー画像で色は違うけど明るさの同じ領域が接しているとモノクロ化したときに明るさに変化がないため輪郭が検出できなくなるんですね。

カラー画像をRGB3枚の画像にバラし、それぞれにCanny変換をかけて再合成してみました。

/*
 * DO NOT FORGET nvcc option : --expt-extended-lambda
 */

// std
#include <iostream>

// OpenCV
#include <opencv2/opencv.hpp>

// CUDA
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <npp.h>

// カーネル関数 二次元のtransform 
//    dst[y][x] = fun(src[y][x]) 
//       where : 0 <= x < width, 0 <= y < height
template<typename T, typename U, typename Function>
__global__ void kernel_transform2D(unsigned int  width, unsigned int height, 
                                        const T* src,         size_t src_pitch,
                                              U* dst,         size_t dst_pitch,
                                       Function  fun) {
  unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
  if ( x < width && y < height ) {
    U* dst_ptr = ((U*)((char*)dst + dst_pitch*y)) + x;
    *dst_ptr = fun(((const T*)((const char*)src + src_pitch*y))[x], *dst_ptr);
  }
}

void color2gray(unsigned int  width, unsigned int height, 
                      uchar3* src,         size_t src_pitch,
                       uchar* dst,         size_t dst_pitch) {
  kernel_transform2D<<<dim3((width+31)/32, (height+7)/8), dim3(32,8)>>>(
    width, height, 
    src, src_pitch, 
    dst, dst_pitch,
    [] __device__ (const uchar3 v, uchar) -> uchar { 
       int t = (v.x + v.y*7 + v.z*2)/10; 
       if ( t <   0 ) t = 0; 
       if ( t > 255 ) t = 255; 
       return (uchar)t; 
    }
  );
}

void color2gray_channel(unsigned int  width, unsigned int height, 
                        uchar3* src, size_t src_pitch,
                        uchar*  dst, size_t dst_pitch,
                       int    channel) {
  kernel_transform2D<<<dim3((width+31)/32, (height+7)/8), dim3(32,8)>>>(
    width, height, 
    src, src_pitch, 
    dst, dst_pitch,
    [=] __device__ (const uchar3 v, uchar) -> uchar 
      { return ((const uchar*)&v)[channel]; }
  );
}

void gray2color_channel(unsigned int  width, unsigned int height, 
                        uchar*  src, size_t src_pitch,
                        uchar3* dst, size_t dst_pitch,
                        int channel) {
  kernel_transform2D<<<dim3((width+31)/32, (height+7)/8), dim3(32,8)>>>(
    width, height, 
    src, src_pitch, 
    dst, dst_pitch,
    [=] __device__ (const uchar v, uchar3 c) -> uchar3 
      { uchar3 t = c; ((uchar*)&t)[channel] = v; return t; }
  );
}

int main(int argc, char *argv[]) {
  cv::VideoCapture camera(0);

  cv::namedWindow("original", CV_WINDOW_AUTOSIZE);
  cv::namedWindow("canny", CV_WINDOW_AUTOSIZE);

  cv::Mat frame;
  cv::Mat canny;

  uchar3* d_frame;
  uchar*  d_gray_base;
  uchar*  d_canny_base;

  uchar*  d_gray[3];
  uchar*  d_canny[3];
  size_t  d_frame_pitch;
  size_t  d_gray_pitch;
  size_t  d_canny_pitch;
  Npp8u*  d_buffer;
  NppiSize size;

 
  // 一発目のキャプチャでフレームのサイズがわかるから
  // (そして多分その後ずっと変わらんだろから)
  // それを基にdevice-memoryを確保
  camera >> frame;

  size.width = (int)frame.size().width;
  size.height = (int)frame.size().height;

  cudaMallocPitch(&d_frame,      &d_frame_pitch, size.width*sizeof(uchar3), size.height);
  cudaMallocPitch(&d_gray_base,  &d_gray_pitch,  size.width,                size.height*3);
  cudaMallocPitch(&d_canny_base, &d_canny_pitch, size.width,                size.height*3);

  for ( size_t i = 0; i < 3; ++i ) {
    d_gray[i]  = d_gray_base  + d_gray_pitch *size.height*i;
    d_canny[i] = d_canny_base + d_canny_pitch*size.height*i;
  }

  // Cannyに引き渡すパラメータ
  NppiSize  nroi    = size;
  NppiPoint noffset = { 0, 0 };
  // 以下のパラメータはイイカンジになるよう適宜調整。
  Npp16s                 nlow_threshold  = 50;
  Npp16s                 nhigh_threshold = 150;
  NppiDifferentialKernel nkernel   = NPP_FILTER_SOBEL;
  NppiMaskSize           nmasksize = NPP_MASK_SIZE_3_X_3;


  // Canny変換に必要なバッファを確保
  {
  int buffer_size;
  nppiFilterCannyBorderGetBufferSize(size, &buffer_size);
  cudaMalloc(&d_buffer, buffer_size);
  }

  canny = frame.clone();
  std::cout 
    << "width,height   = " << size.width << ',' << size.height  
    << "\nstep           = " << frame.step 
    << "\ndepth, channel = " << frame.depth() << ',' << frame.channels()
    << "\n***** [ESC] to exit. *****\n";

  while ( cv::waitKey(10) != 0x1b ) {
    // [1] 画像を frame にキャプチャ
    camera >> frame;
    cv::imshow("original", frame);

    // [2] frame から d_frame へコピー
    cudaMemcpy2D(d_frame, d_frame_pitch, frame.data, frame.step, 
                 size.width*sizeof(uchar3), size.height, cudaMemcpyDefault);

    for ( int i = 0; i < 3; ++i ) {
      // [3] d_frame をモノクロ化して d_gray へ
      color2gray_channel(size.width, size.height, d_frame, d_frame_pitch, d_gray[i], d_gray_pitch,i);

      // [4] d_gray に Canny変換カマして d_canny へ
      nppiFilterCannyBorder_8u_C1R(d_gray[i],  (int)d_gray_pitch,  size, noffset,
                                   d_canny[i], (int)d_canny_pitch, nroi,
                                   nkernel, nmasksize,
                                   nlow_threshold, nhigh_threshold,
                                   nppiNormL2, NPP_BORDER_REPLICATE, 
                                   d_buffer);

      // [5] d_canny をカラー化(RGBを同じ値にするだけ)して d_frame へ
      gray2color_channel(size.width, size.height, d_canny[i], d_canny_pitch, d_frame, d_frame_pitch, i);
    }

    // [6] d_frame を canny へコピー
    cudaMemcpy2D(canny.data, canny.step, d_frame, d_frame_pitch, 
                 size.width*sizeof(uchar3), size.height, cudaMemcpyDefault);

    // [7] 描画!
    cv::imshow("canny", canny);
  }

  // あとしまつ
  cudaFree(d_frame);
  cudaFree(d_gray_base);
  cudaFree(d_canny_base);
  cudaFree(d_buffer);
}

こんなんができましたよ。

NPP : Canny変換

NPP(NVIDIA Performance Primitive) の中に Canny変換 を見つけました。
どうやら CUDA 8.0 で新たに追加されたみたいです。

Canny変換は画像の輪郭を抽出するもので、Sobel/Scharr変換よりシャープな輪郭線を描いてくれます。 Sobel/Scharr変換で得られた輝度勾配の稜線を見つけてくれるってゆーか。

早速試してみました。 OpenCV 3.1 を使ってWeb-cameraからの画像のキャプチャと描画を行います。ダンドリはこんな。

1. Web-cameraからキャプチャした画像を
2. Device-memoryにコピー
3. モノクロ化し
4. Canny変換を施します。
5. 変換後のモノクロ画像をカラー化(RGBに同じ値を入れるだけ)し
6. Hostに書き戻して
7. 描画!

/*
 * DO NOT FORGET nvcc option : --expt-extended-lambda
 */

// std
#include <iostream>

// OpenCV
#include <opencv2/opencv.hpp>

// CUDA
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <npp.h>

// カーネル関数 二次元のtransform 
//    dst[y][x] = fun(src[y][x]) 
//       where : 0 <= x < width, 0 <= y < height
template<typename T, typename U, typename Function>
__global__ void kernel_transform2D(unsigned int  width, unsigned int height, 
                                        const T* src,         size_t src_pitch,
                                              U* dst,         size_t dst_pitch,
                                       Function  fun) {
  unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
  if ( x < width && y < height ) {
    ((U*)((char*)dst + dst_pitch*y))[x] = fun(((const T*)((const char*)src + src_pitch*y))[x]);
  }
}

void color2gray(unsigned int  width, unsigned int height, 
                      uchar3* src,         size_t src_pitch,
                       uchar* dst,         size_t dst_pitch) {
  kernel_transform2D<<<dim3((width+31)/32, (height+7)/8), dim3(32,8)>>>(
    width, height, 
    src, src_pitch, 
    dst, dst_pitch,
    [] __device__ (const uchar3 v) -> uchar { 
       int t = (v.x + v.y*7 + v.z*2)/10; 
       if ( t <   0 ) t = 0; 
       if ( t > 255 ) t = 255; 
       return (uchar)t; 
    }
  );
}

void gray2color(unsigned int  width, unsigned int height, 
                       uchar* src,         size_t src_pitch,
                      uchar3* dst,         size_t dst_pitch) {
  kernel_transform2D<<<dim3((width+31)/32, (height+7)/8), dim3(32,8)>>>(
    width, height, 
    src, src_pitch, 
    dst, dst_pitch,
    [] __device__ (const uchar v) -> uchar3 { return make_uchar3(v,v,v); }
  );
}

int main(int argc, char *argv[]) {
  cv::VideoCapture camera(0);

  cv::namedWindow("original", CV_WINDOW_AUTOSIZE);
  cv::namedWindow("canny", CV_WINDOW_AUTOSIZE);

  cv::Mat frame;
  cv::Mat canny;

  uchar3* d_frame;
  uchar*  d_gray;
  uchar*  d_canny;
  size_t  d_frame_pitch;
  size_t  d_gray_pitch;
  size_t  d_canny_pitch;
  Npp8u*  d_buffer;
  NppiSize size;

 
  // 一発目のキャプチャでフレームのサイズがわかるから
  // (そして多分その後ずっと変わらんだろから)
  // それを基にdevice-memoryを確保
  camera >> frame;
  size.width = (int)frame.size().width;
  size.height = (int)frame.size().height;

  cudaMallocPitch(&d_frame, &d_frame_pitch, size.width*sizeof(uchar3), size.height);
  cudaMallocPitch(&d_gray,  &d_gray_pitch,  size.width,                size.height);
  cudaMallocPitch(&d_canny, &d_canny_pitch, size.width,                size.height);

  // Cannyに引き渡すパラメータ
  NppiSize  nroi    = size;
  NppiPoint noffset = { 0, 0 };
  Npp16s    nlow_threshold  = 50;  // これと
  Npp16s    nhigh_threshold = 150; // これは適宜調整。

  // Canny変換に必要なバッファを確保
  {
  int buffer_size;
  nppiFilterCannyBorderGetBufferSize(size, &buffer_size);
  cudaMalloc(&d_buffer, buffer_size);
  }

  canny = frame.clone();
  std::cout 
    << "width,height   = " << size.width << ',' << size.height  
    << "\nstep           = " << frame.step 
    << "\ndepth, channel = " << frame.depth() << ',' << frame.channels()
    << "\n***** [ESC] to exit. *****\n";

  while ( cv::waitKey(10) != 0x1b ) {

    // [1] 画像を frame にキャプチャ
    camera >> frame;
    cv::imshow("original", frame);

    // [2] frame から d_frame へコピー
    cudaMemcpy2D(d_frame, d_frame_pitch, frame.data, frame.step, 
                 size.width*sizeof(uchar3), size.height, cudaMemcpyDefault);
    // [3] d_frame をモノクロ化して d_gray へ
    color2gray(size.width, size.height, d_frame, d_frame_pitch, d_gray, d_gray_pitch);

    // [4] d_gray に Canny変換カマして d_canny へ
    nppiFilterCannyBorder_8u_C1R(d_gray,  (int)d_gray_pitch,  size, noffset,
                                 d_canny, (int)d_canny_pitch, nroi,
                                 NPP_FILTER_SOBEL, NPP_MASK_SIZE_3_X_3,
                                 nlow_threshold, nhigh_threshold,
                                 nppiNormL2, NPP_BORDER_REPLICATE, 
                                 d_buffer);

    // [5] d_canny をカラー化(RGBを同じ値にするだけ)して d_frame へ
    gray2color(size.width, size.height, d_canny, d_canny_pitch, d_frame, d_frame_pitch);

    // [6] d_frame を canny へコピー
    cudaMemcpy2D(canny.data, canny.step, d_frame, d_frame_pitch, 
                 size.width*sizeof(uchar3), size.height, cudaMemcpyDefault);

    // [7] 描画!
    cv::imshow("canny", canny);
  }

  // あとしまつ
  cudaFree(d_frame);
  cudaFree(d_gray);
  cudaFree(d_canny);
  cudaFree(d_buffer);
}

こんな輪郭線を描いてくれます;

カスタム・アロケータ

CUDAでは目的/用途に応じて様々なメモリの確保/解放APIが用意されています。

1. ピン留めされたHost-memory : cudaMallocHost / cudaFreeHost
2. Host/Device双方で共用できるManaged-memory : cudaMallocManaged / cudaFree
3. Device-memory : cudaMalloc : cudaFree

Device側はさておき、Host側は上記1.2および通常のnew[] / delete[] の３種のメモリ確保/解放を使い分けることになります。

C++屋が(可変長)配列を扱う際。日常的にstd::vectorのお世話になるのですが、std::vector<T>は(デフォルトで)new[]/delete[]が内部的なメモリ管理に使われます。

このメモリ管理をcudaMallocHost/cudaFreeHostあるいはcudaMallocManaged/cudaFreeに差し替えることができれば ピン留めされたvector と Host/Deviceの双方でで共用できるvector が使えて便利。

ってわけで実装しました。ついでに unnique_device_ptr と cuda runtime 例外も。

/* cuda_except.h */
#ifndef CUDA_EXCEPT_H_
#define CUDA_EXCEPT_H_

#include <cuda_runtime.h>
#include <stdexcept>

namespace cu {

class cuda_error : public std::runtime_error {
  cudaError_t err_;
public:
  cuda_error(cudaError_t error) : std::runtime_error(cudaGetErrorString(error)), err_(error) {}
  cudaError_t code() const { return err_; }
  const char* name() const { return cudaGetErrorName(err_); }
};

}
#endif

/* cuda_allocator.h */
#ifndef CUDA_ALLOCATOR_H_
#define CUDA_ALLOCATOR_H_

#include "cuda_except.h"

namespace cu {

template <class T>
struct host_allocator {
  typedef T value_type;
  host_allocator() noexcept {} //default ctor not required by STL
  template<class U> host_allocator(const host_allocator<U>&) noexcept {}
  template<class U> bool operator==(const host_allocator<U>&) const noexcept { return true; }
  template<class U> bool operator!=(const host_allocator<U>&) const noexcept { return false; }
  T* allocate(const size_t n) const;
  void deallocate(T* const p) const noexcept;
  void deallocate(T* const p, size_t) const noexcept { deallocate(p); }
};

template <class T>
struct managed_allocator {
  typedef T value_type;
  managed_allocator() noexcept {} //default ctor not required by STL
  template<class U> managed_allocator(const managed_allocator<U>&) noexcept {}
  template<class U> bool operator==(const managed_allocator<U>&) const noexcept { return true; }
  template<class U> bool operator!=(const managed_allocator<U>&) const noexcept { return false; }
  T* allocate(const size_t n) const;
  void deallocate(T* const p) const noexcept;
  void deallocate(T* const p, size_t) const noexcept { deallocate(p); }
};

template <class T>
struct device_allocator {
  typedef T value_type;
  device_allocator() noexcept {} //default ctor not required by STL
  template<class U> device_allocator(const device_allocator<U>&) noexcept {}
  template<class U> bool operator==(const device_allocator<U>&) const noexcept { return true; }
  template<class U> bool operator!=(const device_allocator<U>&) const noexcept { return false; }
  T* allocate(const size_t n) const;
  void deallocate(T* const p) const noexcept;
  void deallocate(T* const p, size_t) const noexcept { deallocate(p); }
};

#include <cuda_runtime.h>

template <class T>
T* host_allocator<T>::allocate(const size_t n) const {
  if ( n == 0 ) return nullptr;
  if ( n > static_cast<size_t>(-1) / sizeof(T) ) throw std::bad_array_new_length();
  void* pv = nullptr;
  cudaError_t err = cudaMallocHost(&pv, n*sizeof(T));
  if ( err != cudaSuccess ) throw cuda_error(err);
  return static_cast<T*>(pv);
}

template<class T> 
void host_allocator<T>::deallocate(T * const p) const noexcept{
  cudaError_t err = cudaFreeHost(p);
//if ( err != cudaSuccess ) throw cuda_error(err);
}

template <class T>
T* managed_allocator<T>::allocate(const size_t n) const {
  if ( n == 0 ) return nullptr;
  if ( n > static_cast<size_t>(-1) / sizeof(T) ) throw std::bad_array_new_length();
  void* pv = nullptr;
  cudaError_t err = cudaMallocManaged(&pv, n*sizeof(T));
  if ( err != cudaSuccess ) throw cuda_error(err);
  return static_cast<T*>(pv);
}

template<class T> 
void managed_allocator<T>::deallocate(T * const p) const noexcept {
  cudaError_t err = cudaFree(p);
//if ( err != cudaSuccess ) throw cuda_error(err);
}

template <class T>
T* device_allocator<T>::allocate(const size_t n) const {
  if ( n == 0 ) return nullptr;
  if ( n > static_cast<size_t>(-1) / sizeof(T) ) throw std::bad_array_new_length();
  void* pv = nullptr;
  cudaError_t err = cudaMalloc(&pv, n*sizeof(T));
  if ( err != cudaSuccess ) throw cuda_error(err);
  return static_cast<T*>(pv);
}

template<class T> 
void device_allocator<T>::deallocate(T * const p) const noexcept {
  cudaError_t err = cudaFree(p);
//if ( err != cudaSuccess ) throw cuda_error(err);
}

}
#endif

/* unique_device_ptr.h */
#ifndef UNIQUE_DEVICE_PTR_H_
#define UNIQUE_DEVICE_PTR_H_

#include <cuda_runtime.h>
#include <memory>

namespace cu {

template<typename T> struct device_delete {
  device_delete() noexcept = default;
  void operator()(T* ptr) const { cudaFree(ptr); }
};

template<typename T> struct device_delete<T[]> {
  device_delete() noexcept = default;
  void operator()(T* ptr) const { cudaFree(ptr); }
};

template<typename T>
using device_unique_ptr = std::unique_ptr<T,device_delete<T>>;

}

#endif

Windows版CUDA Toolkitではおなじみの配列の足し算をre-writeしてみました。メモリ管理とエラー処理がぐっと楽になります。

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

#include <stdio.h>

#include "cuda_allocator.h"
#include "unique_device_ptr.h"

#include <vector>

template<typename T> using host_vector = std::vector<T, cu::host_allocator<T>>;

void addWithCuda(int *c, const int *a, const int *b, unsigned int size);

__global__ void addKernel(int *c, const int *a, const int *b) {
    int i = threadIdx.x;
    c[i] = a[i] + b[i];
}

inline void cuda_check(cudaError_t status) {
  if (status != cudaSuccess) throw cu::cuda_error(status);
}

int main() {
  try {
    const int arraySize = 5;
    host_vector<int> a = {  1,  2,  3,  4,  5 };
    host_vector<int> b = { 10, 20, 30, 40, 50 };
    host_vector<int> c(arraySize, 0);

    // Add vectors in parallel.
    addWithCuda(c.data(), a.data(), b.data(), arraySize);

    printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
        c[0], c[1], c[2], c[3], c[4]);
  } catch ( const cu::cuda_error& er ) {
    fprintf(stderr, "%s : %s\n", er.name(), er.what());
  }

   // cudaDeviceReset must be called before exiting in order for profiling and
  // tracing tools such as Nsight and Visual Profiler to show complete traces.
  cudaError_t cudaStatus = cudaDeviceReset();
  if (cudaStatus != cudaSuccess) {
    fprintf(stderr, "cudaDeviceReset failed!");
    return 1;
  }

}

// Helper function for using CUDA to add vectors in parallel.
void addWithCuda(int *c, const int *a, const int *b, unsigned int size) {
  cudaError_t cudaStatus;

  cu::device_allocator<int> alloc;

  // Choose which GPU to run on, change this on a multi-GPU system.
  cudaStatus = cudaSetDevice(0);
  cuda_check(cudaStatus);

  cu::device_unique_ptr<int[]> dev_a(alloc.allocate(size));
  cu::device_unique_ptr<int[]> dev_b(alloc.allocate(size));
  cu::device_unique_ptr<int[]> dev_c(alloc.allocate(size));

  cudaStatus = cudaMemcpyAsync(dev_a.get(), a, size*sizeof(int), cudaMemcpyHostToDevice);
  cuda_check(cudaStatus);
  cudaStatus = cudaMemcpyAsync(dev_b.get(), b, size*sizeof(int), cudaMemcpyHostToDevice);
  cuda_check(cudaStatus);

  // Launch a kernel on the GPU with one thread for each element.
  addKernel<<<1, size>>>(dev_c.get(), dev_a.get(), dev_b.get());

  // Check for any errors launching the kernel
  cudaStatus = cudaGetLastError();
  cuda_check(cudaStatus);
    
  cudaStatus = cudaMemcpyAsync(c, dev_c.get(), size*sizeof(int), cudaMemcpyDeviceToHost);
  cuda_check(cudaStatus);

  // cudaDeviceSynchronize waits for the kernel to finish, and returns
  // any errors encountered during the launch.
  cudaStatus = cudaDeviceSynchronize();
  cuda_check(cudaStatus);

}

cuFFT: フーリエ変換で雑音を消す

フーリエ変換(Fourier Transform)は信号処理/解析のド定番。

周期を持ったあらゆる波は異なる周波数のサイン波の重ね合わせで作り出すことができ、与えられた波形から、その成分(サイン波)を割り出すのがフーリエ変換、CUDAにはフーリエ変換ライブラリ: cuFFT が入ってます。軽く使ってみましょうね。

まず入力波を作ります。350,400,450Hzのサイン波を重ね合わせ、さらに一様乱数で生成したノイズを乗せましょう。

  const size_t N = 44100U; // データ数(44.1kHzサンプリングでの1秒分)

  vector<float> signal(N);
  vector<float> h_in(N);

  // 振幅 ±2 のホワイト・ノイズ
  mt19937 mt;
  uniform_real_distribution<float> rnd(-2.0f, 2.0f);

  // 350,400,450Hzのサイン波にノイズを乗せる
  float omega = 2.0f * 3.1416f / N;
  for ( unsigned int i = 0; i < N; ++i ) {
    signal[i] = 
       sinf(omega * 350.0f * (float)i) * 1.0f +
       sinf(omega * 400.0f * (float)i) * 0.8f +
       sinf(omega * 450.0f * (float)i) * 0.6f ;
    h_in[i]   = signal[i] + rnd(mt);
  }

こんなのができました。'赤'はサイン波の合成、それにノイズを乗せたのが'青'です。

この(ノイズまみれの)信号にフーリエ変換を施します。

  // device-memoryの確保(入/出力兼用)
  float* d_real = nullptr;
  cudaMalloc(&d_real, N*sizeof(float));
  float2* d_cplx = reinterpret_cast<float2*>(d_real);

  // フーリエ変換
  cudaMemcpy(d_real, h_in.data(), N*sizeof(float), cudaMemcpyHostToDevice);

  cufftHandle plan_f;
  cufftPlan1d(&plan_f, N, CUFFT_R2C, 1); // Real to Complex (forward)
  cufftExecR2C(plan_f, d_real, d_cplx);

  vector<float2> h_mid(N/2); // スペクトル(フーリエ変換の結果)
  cudaMemcpy(h_mid.data(), d_cplx, N*sizeof(float), cudaMemcpyDeviceToHost);

変換結果がコレ。

350,400,450に大きなピークが見られますね。ノイズは様々な周波数の波がちょっとずつ重なったものなのでグラフの底に貼りつく'モジョモジョ'した部分に現れます。

で、このデータから300Hz以下と500Hz以上の部分をばっさり削ってしまいます。帯域フィルタ(band-pass filter)ってやつです。

  // band-pass filter
  // 300Hz以下/500Hz以上の信号をカットする
  cudaMemset(d_cplx     , 0,      300U  * sizeof(float2));
  cudaMemset(d_cplx+500U, 0, (N/2-500U) * sizeof(float2));

よーするに邪魔なノイズ成分の多くを削り取ったことになります。

しかるのち逆フーリエ変換をかけて、周波数軸から時間軸に戻します。

  // 逆フーリエ変換
  cufftHandle plan_i;
  cufftPlan1d(&plan_i, N, CUFFT_C2R, 1); // Complex to Real (inverse)
  cufftExecC2R(plan_i, d_cplx, d_real);

  // 結果の出力
  vector<float>  h_out(N);
  cudaMemcpy(h_out.data(), d_real, N*sizeof(float), cudaMemcpyDeviceToHost);

結果がコレ。ノイズが消えました♪

コチラ↓が全コード:

/*
 * Noise Reduction with cuFFT
 */
#include <cuda_runtime.h>
#include <cufft.h>

#include <iostream>
#include <random>
#include <vector>
#include <cmath>

using namespace std;

int main() {

  const size_t N = 44100U; // データ数(44.1kHzサンプリングでの1秒分)

  vector<float> signal(N);
  vector<float> h_in(N);

  // 振幅 ±2 のホワイト・ノイズ
  mt19937 mt;
  uniform_real_distribution<float> rnd(-2.0f, 2.0f);

  // 350,400,450Hzのサイン波にノイズを乗せる
  float omega = 2.0f * 3.1416f / N;
  for ( unsigned int i = 0; i < N; ++i ) {
    signal[i] = 
       sinf(omega * 350.0f * (float)i) * 1.0f +
       sinf(omega * 400.0f * (float)i) * 0.8f +
       sinf(omega * 450.0f * (float)i) * 0.6f ;
    h_in[i]   = signal[i] + rnd(mt);
  }

  // device-memoryの確保(入/出力兼用)
  float* d_real = nullptr;
  cudaMalloc(&d_real, N*sizeof(float));
  float2* d_cplx = reinterpret_cast<float2*>(d_real);

  // フーリエ変換
  cudaMemcpy(d_real, h_in.data(), N*sizeof(float), cudaMemcpyHostToDevice);

  cufftHandle plan_f;
  cufftPlan1d(&plan_f, N, CUFFT_R2C, 1); // Real to Complex (forward)
  cufftExecR2C(plan_f, d_real, d_cplx);

  vector<float2> h_mid(N/2); // スペクトル(フーリエ変換の結果)
  cudaMemcpy(h_mid.data(), d_cplx, N*sizeof(float), cudaMemcpyDeviceToHost);

  // band-pass filter
  // 300Hz以下/500Hz以上の信号をカットする
  cudaMemset(d_cplx     , 0,      300U  * sizeof(float2));
  cudaMemset(d_cplx+500U, 0, (N/2-500U) * sizeof(float2));

  // 逆フーリエ変換
  cufftHandle plan_i;
  cufftPlan1d(&plan_i, N, CUFFT_C2R, 1); // Complex to Real (inverse)
  cufftExecC2R(plan_i, d_cplx, d_real);

  // 結果の出力
  vector<float>  h_out(N);
  cudaMemcpy(h_out.data(), d_real, N*sizeof(float), cudaMemcpyDeviceToHost);

  cout << "signal, noised, processed, spectrum" << endl;
  for ( unsigned int i = 0; i < 500; ++i ) {
    cout << signal[i] << ',' 
         << h_in[i] << ',' 
         << h_out[i]/N << ',' 
         << cuCabsf(h_mid[i]) << endl;
  } 

  cudaFree(d_real);
  cufftDestroy(plan_f);
  cufftDestroy(plan_i);

}

(original: 2015-05-21 #36 #37)

cuSOLVER: 鶴亀算を解く

「つる と かめ があわせて３匹います。
　足の数はあわせて10本でした。
　問1: つる と かめ はそれぞれ何匹ですか?
　問2: 足の数があわせて8本なら、それぞれ何匹ですか?」

小学校の算数でやりましたよね、鶴亀算。
コ難しくいえば 連立一次方程式の解を求める問題です。

行列で表現するとこんな。

与えられた A と B から X を求めよ、ってわけ。

CUDA 7.0からこんな問題をCUDAで解いてくれるライブラリ: cuSOLVER がついてます。早速使ってみましたよ。

#include <cuda_runtime.h>
#include <cusolverDn.h> // dense LAPACK

#include <cassert>
#include <iostream>
using namespace std;

template<typename T>
inline size_t bytesof(unsigned int s) { return s * sizeof(T); }

template<typename T>
T* allocate(unsigned int size) {
  T* result = nullptr;
  cudaError_t status = cudaMalloc(&result, bytesof<T>(size));
  assert( status == cudaSuccess );
  return result;
}

int main() {
  cusolverStatus_t status;

  // dense LAPACK
  cusolverDnHandle_t handle;
  status = cusolverDnCreate(&handle);
  assert( status == CUSOLVER_STATUS_SUCCESS );

  int n = 2; // 2x2 正方行列

  float A[] = {  1.0f,  1.0f, 
                 2.0f,  4.0f };
  float* dA = allocate<float>(n*n);
  cudaMemcpy(dA, A, bytesof<float>(n*n), cudaMemcpyHostToDevice);
  int lda = 2;

  // 必要なバッファ量を求め、確保する
  int worksize;
  status = cusolverDnSgetrf_bufferSize(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             &worksize);
  assert( status == CUSOLVER_STATUS_SUCCESS );
#ifdef _DEBUG
  cout << "worksize = " << worksize << endl;
#endif

  float* workspace = allocate<float>(worksize);

  // 計算結果に関する情報
  int* devInfo = allocate<int>(1);
  // ピボット
  int* pivot = allocate<int>(n);

  // LU分解 : dAに結果が求まる(それとpivot)
  status = cusolverDnSgetrf(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             workspace,
             pivot,
             devInfo);
#ifdef _DEBUG
  int info;
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  //         鶴と亀の総数, 足の総数
  float B[] = {  3.0f,       10.0f , 
                 3.0f,        8.0f };
  int nrhs = 2; // 問題数
  float* dB = allocate<float>(n*nrhs);
  cudaMemcpy(dB, B, bytesof<float>(n*nrhs), cudaMemcpyHostToDevice);
  int ldb = 2;

  // AX = B を解く (解XはBをoverrideする)
  status = cusolverDnSgetrs(
             handle,
             CUBLAS_OP_T,
             n,     // 行(=列)
             nrhs,  // 問題数
             dA,    // A
             lda,   // Aのヨコハバ
             pivot, // LU分解で得られたピボット
             dB,    // B
             ldb,   // Bのヨコハバ
             devInfo);
#ifdef _DEBUG
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  // 結果を取得し、出力する
  float X[16];
  cudaMemcpy(X, dB, bytesof<float>(n*nrhs), cudaMemcpyDeviceToHost);
  for ( int i = 0; i < nrhs; ++i ) {
    float* q = B + i*2;
    float* a = X + i*2;
    cout << "総数 = " << q[0] << ", 足の数 = " << q[1]
         << "\t 解 :　"
         <<"鶴 = " << a[0] << " , 亀 = " << a[1] <<endl;
  }

  cudaFree(workspace);
  cudaFree(dA);
  cudaFree(dB);
  cudaFree(devInfo);
  cudaFree(pivot);

  cusolverDnDestroy(handle);
}

実行結果:

(original: 2016-06-03 #40)

streamで速くする(3) ～からくり

CPUとGPUはPCI-busに隔てられてそれぞれが勝手に動くことができます。
CPUはGPUの仕事の完了を待って次の仕事を依頼する必要はないんですわ(そうでないとGPUが重たい仕事してる間CPUがぼーっと待ってにゃならんですから)。仕事の完了を待たずに次の仕事を叩き込むことができるってことは、仕事の待ち行列が用意されてるってことで、それがstreamです。

streamに溜まった仕事(メモリ・コピー と kernel実行)はふたつのエンジンが捌きます。kernel実行をCUDA-core群に割り当てるGigaThread Engineと、メモリ・コピーの依頼を受けてPCI-bus経由でCPU/GPU間のデータ転送を司るCopy Engine、この二つは独立して動けます。

んだから、kernel実行中には次のkernel実行は待たされ、メモリ・コピー中は次のメモリ・コピーは待たされるけど、kernel実行とメモリ・コピーとはそれぞれ異なるEngineで処理されるために同時にやれるってわけですわ。

GeForce GTX9xx, GTX10xx になるとCPU→GPU用 と GPU→CPU用のふたつのCopy Engineを持ってます。PCI-busは上り/下りを同時に転送できるのでより速くなります。

GTX980でのtimelineはこんなカンジ:

緑の帯がCPU→GPU、紫の帯がGPU→CPU 両者がoverlapしてる様子が見て取れます。

streamで速くする(2)

streamはその名の通り"一連の(処理の)流れ"です。 streamに乗せられたjob(仕事)は並んだ順に実行されます。 複数のstreamにjobを乗せておけば、GPUはその中から実行可能なものを一つ選んで処理します。 ホスト・メモリ間のデータ転送にGPUは関与しないので、データ転送の間でも他のstreamに乗ったjobを処理できます。

ここまでにいくつかサンプル・コードを紹介してきましたが、それらの中には stream を明示的に使ったものはありません。 けども実は暗黙の default stream が一本だけあって、その default stream に cudaMemcpy やら kernel呼び出しやらのjobを乗っけていたんです。

３本のstream(ssin, scos, stan)を用意し、sin,cos,tanの表作成をそれぞれの stream に乗せて処理しましょう。

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <math_functions.h>

// 時間稼ぎ
__device__ float time_waste() {
  float result = 0.0f;
  const int TIMES = 10U;
  for ( int i = 0; i < TIMES; ++i ) {
    result += 1.0f / TIMES;
  }
  return result;
}

// 度→ラジアン
__device__ __forceinline__ float deg2rad(float angle) {
  return angle / 180.0f * 3.01416f * time_waste();
}

// sine
__global__ void kernel_sin(float* out, const float* angle, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = sinf(deg2rad(angle[i]));
  }
}

// cosine
__global__ void kernel_cos(float* out, const float* angle, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = cosf(deg2rad(angle[i]));
  }
}

// tangent
__global__ void kernel_tan(float* out, const float* angle, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = tanf(deg2rad(angle[i]));
  }
}

// helper funcs
template<typename T> inline T* allocate_host(size_t size)
  { T* p = nullptr; cudaMallocHost(&p, size*sizeof(T)); return p; }

template<typename T> inline void free_host(T* host)
  { cudaFreeHost(host); }

template<typename T> inline T* allocate_device(size_t size) 
  { T* p = nullptr; cudaMalloc(&p, size*sizeof(T)); return p; }

template<typename T> inline void free_device(T* device)
  { cudaFree(device); }

inline cudaEvent_t create_event() 
  { cudaEvent_t e = nullptr; cudaEventCreate(&e); return e; };

inline void destroy_event(cudaEvent_t event)
  { cudaEventDestroy(event); }

inline cudaStream_t create_stream() 
  { cudaStream_t s = nullptr; cudaStreamCreate(&s); return s; };

inline void destroy_stream(cudaStream_t stream)
  { cudaStreamDestroy(stream); }

#include <iostream>

using namespace std;

const unsigned int N = 1000U * 1000U;
const unsigned int BN = N*sizeof(float);

void multiple_stream() {
  cout << "--- multiple stream ---" << endl;

  float* angle = allocate_host<float>(N); // in
  float* sin   = allocate_host<float>(N); // out
  float* cos   = allocate_host<float>(N); // out
  float* tan   = allocate_host<float>(N); // out

  for ( unsigned int i = 0; i < N; ++i ) {
    angle[i] = 360.0f * i / N;
  }

  float* dangle = allocate_device<float>(N);
  float* dsin   = allocate_device<float>(N);
  float* dcos   = allocate_device<float>(N);
  float* dtan   = allocate_device<float>(N);

  unsigned int size  = N;
  unsigned int block = 128U;
  unsigned int grid  = (size + block -1U) / block;

  cudaEvent_t start = create_event();
  cudaEvent_t mid   = create_event();
  cudaEvent_t stop  = create_event();

  cudaStream_t ssin = create_stream();
  cudaStream_t scos = create_stream();
  cudaStream_t stan = create_stream();

  cudaEventRecord(start);
  cudaMemcpy(dangle, angle, BN, cudaMemcpyHostToDevice);
  cudaEventRecord(mid);

  kernel_sin<<<grid,block,0,ssin>>>(dsin, dangle, size);
  cudaMemcpyAsync(sin, dsin, BN, cudaMemcpyDeviceToHost, ssin);

  kernel_cos<<<grid,block,0,scos>>>(dcos, dangle, size);
  cudaMemcpyAsync(cos, dcos, BN, cudaMemcpyDeviceToHost, scos);

  kernel_tan<<<grid,block,0,stan>>>(dtan, dangle, size);
  cudaMemcpyAsync(tan, dtan, BN, cudaMemcpyDeviceToHost, stan);

  cudaEventRecord(stop);
  cudaDeviceSynchronize();
  float elapsed;
  cudaEventElapsedTime(&elapsed, start, mid);
  cout << elapsed << " + ";
  cudaEventElapsedTime(&elapsed, mid, stop);
  cout << elapsed << " = ";
  cudaEventElapsedTime(&elapsed, start, stop);
  cout << elapsed << " [ms]" << endl;
  for ( unsigned int i = 0; i < 5; ++i ) {
    cout << angle[i] << " :"
         << " sin=" << sin[i] 
         << " cos=" << cos[i] 
         << " tan=" << tan[i] 
         << endl; 
  }
  cout << "..." << endl;

  for ( auto event  : { start, mid,  stop }) 
    destroy_event(event);
  for ( auto stream : { ssin,  scos, stan }) 
    destroy_stream(stream);
  for ( auto host   : { angle,  sin,  cos,  tan  }) 
    free_host(host);
  for ( auto device : { dangle, dsin, dcos, dtan }) 
    free_device(device);

}

int main() {
  multiple_stream();
}

default stream １本で処理した時のタイムラインがコレ。

対して３本(+default)のstreamでやったのがコチラ。

ね、データ転送とkernel実行が同時に行われてるでしょ。重なった分だけ速くなるんです。

実装上のキモはふたつ:

• cudaMemcpyじゃなく cudaMemcpyAsyncを使うべし。 後者はコピー完了を待たずに(job投入後ただちに)返って来るので、すぐさま次のjobを投入できます。

• ホスト側メモリを malloc/free や new/delete じゃなく cudaMallocHost/cudaFreeHost で確保/解放すべし。 そうしないとデータ転送中に kernelが動いてくれないのですよ。

... streamのおハナシはまだまだ続きます。

(original: 2015-04-16)