OpenCVのHALモジュールと高速な数値計算ライブラリ

はじめに

これは，OpenCV Advent Calendar 2015の記事です．関連記事は，リンク先に目次としてまとめられています．

HALモジュール

OpenCVはhal(Hardware Acceleration Layer)モジュールを提供しています．将来的には，様々なハードウェアに最適化が関数を集めたモジュールにして提供していく予定のようですが，現在では，SIMDベクトル化を行った少量の関数群が提供されています．（3.1のアップデートで拡大したいということがmeeting noteにありました．）

これは，実装したいアルゴリズムの中で良く使われる非常にコアな関数を，各ハードウェアに依存しながらもより高速に実装したモジュールへと分割することで，高速化と汎用性を両立すること目指してして作ったものと想定されます．例えば，精度をほんのわずかに落とせば，Mat全体をexpやlog，powをするような機能を高速化することが出来ます．他にもcos, sin, tan, arctanなど，様々な関数が高速化されるとメリットがあります．

そこで，（かなり未完成ですが）現在，exp，log，Atan，sqrt，sqrtの逆数計算やノルム計算，LU分解，チェロスキー分解が高速化されています．以下はHALモジュール内で実体がある関数の抜粋です．これ以外にもファイルはありますが，ほとんどがヘッダと空のcppだけ用意された状況です．しかし将来的には，色変換やフィルタ，統計処理，行列演算などを将来的には拡張していくことが予想されます．

int LU(float* A, size_t astep, int m, float* b, size_t bstep, int n);
int LU(double* A, size_t astep, int m, double* b, size_t bstep, int n);
bool Cholesky(float* A, size_t astep, int m, float* b, size_t bstep, int n);
bool Cholesky(double* A, size_t astep, int m, double* b, size_t bstep, int n);

int normL1_(const uchar* a, const uchar* b, int n);
float normL1_(const float* a, const float* b, int n);
float normL2Sqr_(const float* a, const float* b, int n);

void exp(const float* src, float* dst, int n);
void exp(const double* src, double* dst, int n);
void log(const float* src, float* dst, int n);
void log(const double* src, double* dst, int n);

void fastAtan2(const float* y, const float* x, float* dst, int n, bool angleInDegrees);
void magnitude(const float* x, const float* y, float* dst, int n);
void magnitude(const double* x, const double* y, double* dst, int n);
void sqrt(const float* src, float* dst, int len);
void sqrt(const double* src, double* dst, int len);
void invSqrt(const float* src, float* dst, int len);
void invSqrt(const double* src, double* dst, int len);

高速なexp計算

512x512の画像をexpするのに，HALモジュールを使うとどれくらい高速化されるのかテストしてみましょう．

まず，coreモジュールの普通のexpを使ってみます．
cv::exp(srcf, dest);
これを使うと，1.35msで計算することが出来ました．

これを，halモジュールのexpに切り替えてみましょう．
cv::hal::exp(srcf.ptr<float>(0), dest.ptr<float>(0), srcf.size().area());
Matを受け付けず，ポインタで渡す必要があるなど使い方が少し異なりますが，これで計算すると，なんと0.66msと半分以下で計算できます．

おまけ
さらに，expやlogでより速いものが欲しいひとは，下記リンク先のfmath.hppを使ってみて下さい．非常に高速です．上記と同じ画像を0.593msで計算できます．以下はサンプルコードです．SIMD等を使っていますが，詳しくは，前回の記事のリンクでもある，組み込み関数（intrinsic）によるSIMD入門を参照してください．

float* s = srcf.ptr<float>(0);
float* d = dest.ptr<float>(0);
int size = srcf.size().area() / 4;
for (int i = 0; i < size; i++)
    {
        __m128 ms = _mm_load_ps(s);
        _mm_store_ps(d, fmath::exp_ps(ms));
        s += 4, d += 4;
    }

行列とEigen

現在，HALモジュールのLU分解やCholesky分解は，C言語で余計なオーバーヘッドがおきないように書かれている程度で，コード自体の最適化はされているわけではありません．試してはいませんが，大きな行列では，行列ライブラリであるEigenを使ったほうが速いでしょう．このライブラリは，フリーで公開されているものの中ではかなり速いものでありかつ使いやすいです．しかし，行列が小さい場合は，OpenCVからEigenの形式への変換のオーバーヘッドほうが大きくなるため，このような小幅な改善でも効果があると思います（コードにもfor small matrixとあります．）．

WITH_EIGENでOpenCVをコンパイルするとcv::cv2eigenとcv::eigen2cvの各関数でOpenCVのMat形式からEigenの行列形式(Eigen::Matrix)に型変換できるようになります．もし巨大な行列の演算が主役のプログラムであればOpenCVよりもEigenを使うべきでしょう．

IntelのマルチメディアライブラリIPP

IPPはIntelから公開されている商用の高速化ライブラリです．OpenCVはもともとIntelから提供されていた経緯もあり最近一部機能を無料でOpenCVに提供してくれています．公開されているライブラリでは，デフォルトでIPPの機能は有効化されています．下記にIPPのオンオフ，状態の取得をする関数を示します．IPPの関数をoffにした状態で速度を計測してみるとその恩恵が分かるでしょう．

namespace ipp
{
CV_EXPORTS void setIppStatus(int status, const char * const funcname = NULL, const char * const filename = NULL, int line = 0);
CV_EXPORTS int getIppStatus();
CV_EXPORTS String getIppErrorLocation();
CV_EXPORTS bool useIPP();
CV_EXPORTS void setUseIPP(bool flag);
} // ipp

このIPPはFFTの高速化も実現してくれています．IPPが導入されるまでのOpenCVではFFTが非常に遅く，IPP導入のおかげで劇的にかわったポイントの一つです．ただし，入力する信号によっては，別の高速なFFTライブラリであるFFTW を使ったほうが高速な場合もあります．またDCTにいたっては，昔のOpenCVでは高速DCTではなく離散コサイン変換をベタで書いてあったため，IPP導入の速度向上が破格です．（それでも，自前で作ったほうが速いですが．．．（後述））

OpenCVの開発チームは，時期3.1ではこのIPPのサポートを拡大すると述べています．現に，ライブラリ関数をincludeするヘッダであるipp.hは，最新版の3.0時のヘッダでは6.5K行だったのが，現在の開発版のヘッダでは8.5K行程度に成長しており，今後の拡大が期待できます．
例えば，3.0では，バイラテラルフィルタの実装ippiFilterBilateral_8u_C1Rが宣言だけされて無効化されていたのですが，ipp.hにbilateralFilterが定義されているため，公開されるようになるかもしれません．

なお，IPPがサポートされている関数は書きの#defineマクロによって，IPPが使えるか？をチェックして関数を呼び出すかどうかをスイッチしています．

#define CV_IPP_CHECK_COND (cv::ipp::useIPP())
#define CV_IPP_CHECK() if(CV_IPP_CHECK_COND)

その他のライブラリ

OpenCV3.0は，線形計画法のシンプレックス法や，最適化をしたいときに使えるDownhill simplex法，共役勾配法などの最適化アルゴリズムを提供し始めるなど新たな数値計算に関連する関数の検討も始めています．しかし，まだライブラリとしては成長段階であり出来があまりよくありません．前述の行列の章で述べたEigenライブラリなどのように，現在では，行列演算や微積分，最適化などで性能が欲しい場合は，外部ライブラリを使うことを検討するべきです．

現状，コンピュータビジョンのライブラリであるOpenCVでは，固有値分解や主成分分析といったコンピュータビジョンで多用する演算はライブラリ内に用意されているもののあまり優れた実装ではありません．またDCTや基本的なものも遅いです．例えば，符号化に使うDCTを例に挙げて説明します（筆者が符号化の研究をしているのでその他の話よりは詳しいので選びます．あまり数値計算の主役とは言いがたいですが．．．）

DCTは静止画像の符号化であるJPEGや動画像の符号化であるH.264/AVC，H.265/HEVCで使われており，世界で最も使われている数値計算アルゴリズムと言っても過言ではありません．このDCTは，4x4の場合，8x8の場合,16x16の場合，32x32の場合，それらを固定小数点であらわした整数で演算した場合，精度を高めた整数演算をした場合，整数の演算で大体直交かさせた場合など各場合において，最低何回の乗算，加算で実現できるかを理論的に解析したり，実践的に実装したものが論文誌として発表されるくらいに高度な議論が積み重なっています．例えば，2005年に公開された論文を，私がSSEを使って実装した4x4,8x8,16x16の浮動小数点DCTの実装は，OpenCVやIPPのDCTよりも圧倒的な速度を持っています．(このDCTを使った，高速なノイズ除去法:乱択冗長DCTもコードを公開しているので，良かったら使ってみてください．)

まとめ

この数値計算自が主要な研究分野であり，この分野の成果をコンピュータビジョンのライブラリひとつで吸収していくことを期待するのは非現実的なことかもしれません．

そこで，必要なときに必要なライブラリを使っていきましょう．有料でいいならIntelから出ているMKLやIPPは非常に優秀なライブラリです．もし高速なコードが必要になったと思ったら，まずIPPやMKLでサポートしていないか調べましょう．また，Wikipediaの数値解析ソフトウェアの項にもリンクが充実しています．

宣伝: もし，個別に最適化して欲しいものがあれば @fukushima1981 までご相談ください．勤務先の連絡先でOKです．共同研究等の形が取れれば可能かもしれません．