画像ピクセル値へのアクセスと計算速度

03 3月 2010 Under: opencv2.x-samples

C

#include <cv.h>
#include <highgui.h>
#include <stdio.h>

int
main (int argc, char **argv)
{
  int x, y;
  int width=1025, height=978;
  IplImage *img=0, *hsv_img=0;
  double c, f;
  
  // (1)allocate and initialize an image
  img = cvCreateImage(cvSize(width, height), IPL_DEPTH_8U, 3);
  if(img==0) return -1;
  cvZero(img);
  hsv_img = cvCloneImage(img);
  cvCvtColor(img, hsv_img, CV_BGR2HSV);
  f = cvGetTickFrequency()*1000;

  // (2)create hue-value gradation image
  c = cvGetTickCount();
  for(y=0; y<height; y++) {
    for(x=0; x<width; x++) {
      int a = img->widthStep*y+(x*3);
      hsv_img->imageData[a+0] = (x*180/width);
      hsv_img->imageData[a+1] = 255;
      hsv_img->imageData[a+2] = ((height-y)*255/height);
    }
  }
  printf("%f\n", (cvGetTickCount()-c)/f);

  cvCvtColor(hsv_img, img, CV_HSV2BGR);

  // (3)show the iamge, and quit when any key pressed
  cvNamedWindow ("Gradation", CV_WINDOW_AUTOSIZE);
  cvShowImage ("Gradation", img);
  cvWaitKey (0);

  cvDestroyWindow("Gradation");
  cvReleaseImage(&img);
  
  return 0;
}

// (1)指定サイズの画像を作成します.

関数 cvCreateImage によって,指定サイズの画像を作成します.確保されるデータ領域は,4バイト境界にアラインメントが揃えられます.このようにして作成された画像の色空間を,関数 cvCvtColor によって RGBからHSVに変換します.
また,処理時間を計測するために,1ミリ秒毎の tick 数を求めます. C++ インタフェースの getTickFrequency が 1 秒毎の tick 数を返すのに対して, cvGetTickFrequency は 1 マイクロ秒毎の tick 数を返すことに注意してください.

// (2)色相-明度のグラデーション画像を作成します.

IplImage 構造体の imageData へのポインタを取得し,色相と明度のグラデーション画像を作成します. [0] が色相 “Hue”,[1] が彩度 “Saturation”,[2] が明度 “Value” になります.

// (3)画像を表示します.

作成されたグラデーション画像を表示し,何かキーが押されるとプログラムを終了します.

C++

#include <cv.h>
#include <highgui.h>
#include <iostream>

using namespace cv;

#define NO_GAP

int
main (int argc, char **argv)
{
  const int w=1025, h=978;

  // (1)create an image with specified size
#ifdef NO_GAP
  Mat img = Mat::zeros(Size(w,h), CV_8UC3);
  Mat hsv_img;
  cvtColor(img, hsv_img, CV_BGR2HSV);
#else
  Mat img = Mat::zeros(Size(w+1,h), CV_8UC3);
  cvtColor(img, img, CV_BGR2HSV);
  Rect roi_rect(0,0,w,h);
  Mat hsv_img(img, roi_rect);
#endif
  int ch = img.channels(); // =3
  std::cout << "step:" << hsv_img.step << ", size:" << w*ch << std::endl;
  double c, f = getTickFrequency()/1000;
  int width = hsv_img.cols, height = hsv_img.rows;

  // (2a)create hue-value gradation image
  c = getTickCount();
  for(int y=0; y<height; ++y) {
    for(int x=0; x<width; ++x) {
      int a = hsv_img.step*y+(x*ch);
      hsv_img.data[a+0] = (x*180/width);
      hsv_img.data[a+1] = 255;
      hsv_img.data[a+2] = ((height-y)*255/height);
    }
  }
  std::cout << (getTickCount()-c)/f << "\t";

  // (2b)create hue-value gradation image
  c = getTickCount();
  for(int y=0; y<height; ++y) {
    uchar *p = hsv_img.ptr(y);
    for(int x=0; x<width; ++x) {
      p[x*ch+0] = (x*180/width);
      p[x*ch+1] = 255;
      p[x*ch+2] = ((height-y)*255/height);
    }
  }
  std::cout << (getTickCount()-c)/f << "\t";

  // (2c)create hue-value gradation image
  c = getTickCount();
  if(hsv_img.isContinuous()) {
    width *= height;
    height = 1;
  }
  for(int y=0; y<height; ++y) {
    uchar *p = hsv_img.ptr(y);
    for(int x=0; x<width; ++x) {
      p[x*ch+0] = (x*180/width);
      p[x*ch+1] = 255;
      p[x*ch+2] = ((height-y)*255/height);
    }
  }
  std::cout << (getTickCount()-c)/f << "\t";
  width = hsv_img.cols;
  height = hsv_img.rows;
  
  // (2d)create hue-value gradation image
  c = getTickCount();
  for(int y=0; y<height; ++y) {
    for(int x=0; x<width; ++x) {
      Vec3b &p = hsv_img.at<Vec3b>(y,x);
      p[0] = (x*180/width);
      p[1] = 255;
      p[2] = ((height-y)*255/height);
    }
  }
  std::cout << (getTickCount()-c)/f << "\t";

  // (2e)create hue-value gradation image
  c = getTickCount();
  Mat_<Vec3b>& ref_img = (Mat_<Vec3b>&)hsv_img;
  for(int y=0; y<height; ++y) {
    for(int x=0; x<width; ++x) {
      ref_img(y, x)[0] = (x*180/width);
      ref_img(y, x)[1] = 255;
      ref_img(y, x)[2] = ((height-y)*255/height);
    }
  }
  std::cout << (getTickCount()-c)/f << "\t";

  // (2f)create hue-value gradation image
  c = getTickCount();
  vector<Mat> planes;
  split(hsv_img, planes);
  MatIterator_<uchar> it_h = planes[0].begin<uchar>();
  MatIterator_<uchar> it_s = planes[1].begin<uchar>();
  MatIterator_<uchar> it_v = planes[2].begin<uchar>();
  for(int c=0; it_h!=planes[0].end<uchar>(); ++it_h, ++it_s, ++it_v, ++c) {
    *it_h = ((c%width)*180/width);
    *it_s = 255;
    *it_v = ((height-c/width)*255/height);
  }
  merge(planes, hsv_img);
  std::cout << (getTickCount()-c)/f << "\t";

  // (2g)create hue-value gradation image (when no-gap)
#ifdef NO_GAP
  c = getTickCount();
  MatIterator_<Vec3b> it = hsv_img.begin<Vec3b>();
  for(int c=0; it!=hsv_img.end<Vec3b>(); ++it,++c) {
    int x=c%width;
    int y=c/width;
    (*it)[0] = (x*180/width);
    (*it)[1] = 255;
    (*it)[2] = ((height-y)*255/height);
  }
  std::cout << (getTickCount()-c)/f << "\t";
#endif
  std::cout << std::endl;
  
  // (3)show created gradation image
  cvtColor(hsv_img, img, CV_HSV2BGR);
  namedWindow("Garadation", CV_WINDOW_AUTOSIZE);
  imshow("Gradation", img);
  waitKey(0);
  
  return 0;
}

// (1)指定サイズの画像を作成します.

行間にギャップが存在しない領域を持つMatインスタンスを作成するには,単に Mat コンストラクタ, Mat::zeros などを利用します.逆に,ギャップを持つ画像を作成する方法の一つには,対象画像よりも横幅が大きい画像を作成し,その部分領域(部分行列,ROI)として画像を作成する方法があります.実際に,Mat.stepと width の3倍(チャンネル数倍)を比較してみると,行間にギャップが存在する事が確認できます.また, cvCreateImagecvLoadImage などの IplImage 構造体を対象とする関数が確保するデータ領域は4バイト境界にアラインメントが揃えられるので,次の用にするとギャップを持つ画像領域が確保されることになります.

Ptr<IplImage> ipl_img = cvCreateImage(cvSize(w,h), 8, 3);
cvCvtColor(ipl_img, ipl_img, CV_BGR2HSV);
Mat img, hsv_img(ipl_img);

このようにして作成された画像の色空間を,関数 cvtColor によって RGBからHSVに変換します.
また,処理時間を計測するために,1ミリ秒毎の tick 数を求めます.C インタフェースの cvGetTickFrequency が 1 マイクロ秒毎の tick 数を返すのに対して, getTickFrequency は 1 秒毎の tick 数を返すことに注意してください.

// (2a)色相-明度のグラデーション画像を作成します.

この方法は,OpenCV 1.x での処理とほぼ同じです.Mat.step を IplImage の widthStep と読み替えると,容易く理解できるでしょう.

// (2b)色相-明度のグラデーション画像を作成します.

Mat クラスには,指定行の先頭へのポインタを返すメソッド Mat.ptr() が存在します.これの引数として y (行のインデックス)を与えることで,(2a) と同様の処理が可能です.実際,OpenCV 内部でもこれらの処理はほぼ等価なものになります.

// (2c)色相-明度のグラデーション画像を作成します.

これは (2b) の処理を少し最適化したものです.あらかじめギャップの有無を調べ,ギャップが存在しない場合は,MxN の行列を 1xMN の 1行の行列として扱うことで,列のループ無くします.

// (2d)色相-明度のグラデーション画像を作成します.

これは,Mat.at メソッドを用いて,指定座標 (x,y) の参照を取得する方法です( at の引数は y, x の順になることに注意してください).

// (2e)色相-明度のグラデーション画像を作成します.

(2d) のように Mat.at を用いる代わりに,あらかじめ Mat_ クラスへキャストしたインスタンスの参照を取得し,そこに定義された()演算子を用いて指定座標 (x, y) への参照を取得します.

// (2f)色相-明度のグラデーション画像を作成します.

イテレータを利用する例です.Mat クラスのメソッドが返すイテレータは,チャンネル毎のイテレータなので,まず画像をチャンネル毎に分離してから処理を行います.また,最後に処理結果をマージする必要があります.

// (2g)色相-明度のグラデーション画像を作成します.

(無理に)イテレータを利用する例です.3チャンネルを1平面と見なし,イテレータにより処理を行いますが,当然ギャップを飛び越すことができないので,ギャップなしの画像に対してのみ用いることができます.

// (3)画像を表示します.

作成されたグラデーション画像を表示し,何かキーが押されるとプログラムを終了します.

———-

(2a) (2b) (2c) (2d) (2e) (2f) (2g)
without gap 10.8366 10.524 10.164 18.940 18.973 73.622 60.040
with gap 10.000 10.261 10.221 19.191 19.108 73.484

上記の表は,処理時間を計測した結果の一例です(100試行の平均 [ms] ).(2a),(2b),(2c) は,有効桁数を考えるとほぼで同じであり,(2d)と(2e)も同様に同じ処理時間である一方,(2f),(2g)の例は極端に遅くなっています.

(追記→)
gcc -O0

(2a) (2b) (2c) (2d) (2e) (2f) (2g)
without gap 10.837 8.613 8.552 22.531 30.331 136.910 116.788
with gap 10.930 9.260 8.711 23.410 28.617 136.247

gcc -O3

(2a) (2b) (2c) (2d) (2e) (2f) (2g)
without gap 3.620 3.390 3.488 3.670 4.909 57.061 46.872
with gap 3.658 3.408 3.392 3.690 4.361 57.937

上記の表は,処理時間を計測した結果の一例です(100試行の平均 [ms] ).また,最適化オプションによって,(2a),(2b),(2c),(2d),(2e) まで,ほぼ変わらない速度で実行されていることが分かります.

(←追記)

画像内のピクセルを横断するような処理の多くは,処理対象となる画素の空間情報(どこにその画素が位置するか)を利用するため,たとえ1チャンネル画像であったとしても,MatIterator_ を用いて処理するには向いていません.また,今回の場合のようにイテレータがクラスとして実装されている場合,ポインタを直接操作する処理と比較して時間がかかります.Mat.dataはメモリ上に連続している必要は無いため,イテレータがインクリメント/デクリメントされる度に,その指す場所がスライス(2次元画像ならば「行」)から外れていないかを確認し,外れているならば適切な場所を指すようにアドレスを再計算する必要があるためです.また,OpenCV2.2以降,MatNDがMatと統合されたため,行列の次元(Mat.dims>=2)が高くなることもあるので,その場合はこの計算によりコストがかかります.

しかし,例えば以下のような場合では,イテレータを用いて処理を簡単に記述できます.

ある行列において,値が正である要素の合計値を求める処理を考えてみます.

・行と列のループで処理する場合

double sum=0;
int cols = M.cols, rows = M.rows;
for(int i = 0; i < rows; i++) {
  const double* Mi = M.ptr<double>(i);
  for(int j = 0; j < cols; j++) {
    sum += std::max(Mi[j], 0.);
  }
}

・イテレータを利用して処理する場合

double sum=0;
MatConstIterator_<double> it = M.begin<double>(), it_end = M.end<double>();
for(; it != it_end; ++it) {
  sum += std::max(*it, 0.);
}

このように,処理対象のインデックスを必要としない処理の場合は特に,イテレータを利用して簡潔かつ安全に処理を行うことが可能です.また,STLの任意のアルゴリズムに適用できるという利点もあります(MatIterator_ はランダムアクセスイテレータです).

実行結果例

色相-明度のグラデーション画像.色相は左から右(0-180),明度は上から下(255-0)に変化します.