目標

Thrustは各種CUDAアクセラレーテッドアルゴリズムのための非常に強力なライブラリである。しかしthrustはベクトルを対象に設計されており、ピッチ付き行列（pitched matrix）を扱うようには作られていない。以下のチュートリアルでは、thrustアルゴリズムで利用できるthrustイテレータへ cv::cuda::GpuMat をラップする方法を解説する。

このチュートリアルでは次の方法を示す。

GpuMatをthrustイテレータにラップする
GpuMatを乱数で埋める
GpuMatの1列をインプレースでソートする
0より大きい値を新しいGPU行列にコピーする
thrustでストリームを使う

GpuMatをthrustイテレータにラップする

以下のコードはGpuMatのイテレータを生成する。

/*
    @Brief GpuMatBeginItr returns a thrust compatible iterator to the beginning of a GPU mat's memory.
    @Param mat is the input matrix
    @Param channel is the channel of the matrix that the iterator is accessing.  If set to -1, the iterator will access every element in sequential order
*/
template<typename T>
thrust::permutation_iterator<thrust::device_ptr<T>, thrust::transform_iterator<step_functor<T>, thrust::counting_iterator<int>>>  GpuMatBeginItr(cv::cuda::GpuMat mat, int channel = 0)
{
 if (channel == -1)
    {
        mat = mat.reshape(1);
        channel = 0;
    }
 CV_Assert(mat.depth() == cv::DataType<T>::depth);
 CV_Assert(channel < mat.channels());
 return thrust::make_permutation_iterator(thrust::device_pointer_cast(mat.ptr<T>(0) + channel),
        thrust::make_transform_iterator(thrust::make_counting_iterator(0), step_functor<T>(mat.cols, mat.step / sizeof(T), mat.channels())));
}

/*
@Brief GpuMatEndItr returns a thrust compatible iterator to the end of a GPU mat's memory.
@Param mat is the input matrix
@Param channel is the channel of the matrix that the iterator is accessing.  If set to -1, the iterator will access every element in sequential order
*/
template<typename T>
thrust::permutation_iterator<thrust::device_ptr<T>, thrust::transform_iterator<step_functor<T>, thrust::counting_iterator<int>>>  GpuMatEndItr(cv::cuda::GpuMat mat, int channel = 0)
{
 if (channel == -1)
    {
        mat = mat.reshape(1);
        channel = 0;
    }
 CV_Assert(mat.depth() == cv::DataType<T>::depth);
 CV_Assert(channel < mat.channels());
 return thrust::make_permutation_iterator(thrust::device_pointer_cast(mat.ptr<T>(0) + channel),
        thrust::make_transform_iterator(thrust::make_counting_iterator(mat.rows*mat.cols), step_functor<T>(mat.cols, mat.step / sizeof(T), mat.channels())));
}

ここでの目標は、行列の先頭から始まり、連続する行列要素にアクセスできるよう正しくインクリメントするイテレータを得ることである。連続した行についてはこれは自明だが、ピッチ付き行列の1列についてはどうだろうか。これを実現するには、イテレータが行列の次元とステップを把握している必要がある。この情報はstep_functorに埋め込まれている。

template<typename T> struct step_functor : public thrust::unary_function<int, int>
{
 int columns;
 int step;
 int channels;
    __host__ __device__ step_functor(int columns_, int step_, int channels_ = 1) : columns(columns_), step(step_), channels(channels_)  {   };
    __host__ step_functor(cv::cuda::GpuMat& mat)
    {
 CV_Assert(mat.depth() == cv::DataType<T>::depth);
        columns = mat.cols;
        step = mat.step / sizeof(T);
        channels = mat.channels();
    }
    __host__ __device__
 int operator()(int x) const
 {
 int row = x / columns;
 int idx = (row * step) + (x % columns)*channels;
 return idx;
    }
};

ステップファンクタはインデックス値を受け取り、行列の先頭からの適切なオフセットを返す。カウンティングイテレータは単純にピクセル要素の範囲をインクリメントする。これらをtransform_iteratorに組み合わせると、0からM*Nまでカウントし、GpuMatのピッチ付きメモリを考慮して正しくインクリメントするイテレータが得られる。残念ながらこれにはメモリ位置の情報が一切含まれていないため、それにはthrust::device_ptrが必要になる。デバイスポインタをtransform_iteratorと組み合わせることで、thrustに行列の先頭要素を指し示し、それに応じてステップさせることができる。

GpuMatを乱数で埋める

thrust用のイテレータを作る便利な関数が揃ったので、これらを使ってOpenCV単体ではできないことをやってみよう。残念ながら本稿執筆時点で、OpenCVにはGPUでの乱数生成機能が一切ない。幸いthrustにはそれがあり、両者の相互運用は今や容易である。例は http://stackoverflow.com/questions/12614164/generating-a-random-number-vector-between-0-and-1-0-using-thrust より引用。

まず乱数値を生成するファンクタを書く必要がある。

struct prg
{
  float a, b;
 
  __host__ __device__
    prg(float _a = 0.f, float _b = 1.f) : a(_a), b(_b) {};
 
  __host__ __device__
    float operator()(const unsigned int n) const
  {
    thrust::default_random_engine rng;
    thrust::uniform_real_distribution<float> dist(a, b);
    rng.discard(n);
    return dist(rng);
  }
};

これは整数値を受け取り、aとbの間の値を出力する。次に、thrustのtransformを使って行列を0から10の間の値で埋める。

  {
    cv::cuda::GpuMat d_value(1, 100, CV_32F);
    auto valueBegin = GpuMatBeginItr<float>(d_value);
    auto valueEnd = GpuMatEndItr<float>(d_value);
    thrust::transform(thrust::make_counting_iterator(0), thrust::make_counting_iterator(d_value.cols), valueBegin, prg(-1, 1));
 
    cv::Mat h_value(d_value);
  }

GpuMatの1列をインプレースでソートする

行列の要素を乱数値とインデックスで埋めてみよう。その後、乱数とインデックスをソートする。

  {
    cv::cuda::GpuMat d_data(1, 100, CV_32SC2);
    // Thrust compatible begin and end iterators to channel 1 of this matrix
    auto keyBegin = GpuMatBeginItr<int>(d_data, 1);
    auto keyEnd = GpuMatEndItr<int>(d_data, 1);
    // Thrust compatible begin and end iterators to channel 0 of this matrix
    auto idxBegin = GpuMatBeginItr<int>(d_data, 0);
    auto idxEnd = GpuMatEndItr<int>(d_data, 0);
    // Fill the index channel with a sequence of numbers from 0 to 100
    thrust::sequence(idxBegin, idxEnd);
    // Fill the key channel with random numbers between 0 and 10.  A counting iterator is used here to give an integer value for each location as an input to prg::operator()
    thrust::transform(thrust::make_counting_iterator(0), thrust::make_counting_iterator(d_data.cols), keyBegin, prg(0, 10));
    // Sort the key channel and index channel such that the keys and indices stay together
    thrust::sort_by_key(keyBegin, keyEnd, idxBegin);
 
    cv::Mat h_idx(d_data);
  }

ストリームを使いながら0より大きい値を新しいGPU行列にコピーする

この例では、cv::cuda::Streamsをthrustとどう併用できるかを見ていく。残念ながらこの特定の例では、結果をCPUに返さなければならない関数を使っているため、ストリームの最適な使い方ではない。

  {
    cv::cuda::GpuMat d_value(1, 100, CV_32F);
    auto valueBegin = GpuMatBeginItr<float>(d_value);
    auto valueEnd = GpuMatEndItr<float>(d_value);
    cv::cuda::Stream stream;
    //! [random_gen_stream]
    // Same as the random generation code from before except now the transformation is being performed on a stream
    thrust::transform(thrust::system::cuda::par.on(cv::cuda::StreamAccessor::getStream(stream)), thrust::make_counting_iterator(0), thrust::make_counting_iterator(d_value.cols), valueBegin, prg(-1, 1));
    //! [random_gen_stream]
    // Count the number of values we are going to copy
    int count = thrust::count_if(thrust::system::cuda::par.on(cv::cuda::StreamAccessor::getStream(stream)), valueBegin, valueEnd, pred_greater<float>(0.0));
    // Allocate a destination for copied values
    cv::cuda::GpuMat d_valueGreater(1, count, CV_32F);
    // Copy values that satisfy the predicate.
    thrust::copy_if(thrust::system::cuda::par.on(cv::cuda::StreamAccessor::getStream(stream)), valueBegin, valueEnd, GpuMatBeginItr<float>(d_valueGreater), pred_greater<float>(0.0));
    cv::Mat h_greater(d_valueGreater);
  }

まず、ストリーム上で-1から1の間でランダムに生成したデータをGPU matに埋める。

// Same as the random generation code from before except now the transformation is being performed on a stream

thrust::transform(thrust::system::cuda::par.on(cv::cuda::StreamAccessor::getStream(stream)), thrust::make_counting_iterator(0), thrust::make_counting_iterator(d_value.cols), valueBegin, prg(-1, 1));

thrust::system::cuda::par.on(...)の使用に注目してほしい。これはthrustコードをストリーム上で実行するための実行ポリシーを作成する。cudaツールキットに同梱されているバージョンのthrustにはバグがあり、バージョン7.5の時点でこれは修正されていない。このバグによりコードがストリーム上で実行されない。ただしこのバグは、gitリポジトリ（http://github.com/thrust/thrust.git）から最新版のthrustを使うことで修正できる。次に、以下の述語を使ったthrust::count_ifにより、0より大きい値がいくつあるかを判定する。

template<typename T> struct pred_greater
{
  T value;
  __host__ __device__ pred_greater(T value_) : value(value_){}
  __host__ __device__ bool operator()(const T& val) const
  {
    return val > value;
  }
};

その結果を使って、コピーした値を格納する出力バッファを作成し、続いて同じ述語を用いたcopy_ifで出力バッファを埋める。最後に、確認のためその値をCPU matにダウンロードする。

目次

目標

GpuMatをthrustイテレータにラップする

GpuMatを乱数で埋める

GpuMatの1列をインプレースでソートする

ストリームを使いながら0より大きい値を新しいGPU行列にコピーする