OpenCV学習ノート（27）KAZEアルゴリズム原理とソース分析（一）非線形拡散フィルタリング
19732 ワード
KAZEアルゴリズムリソース：
1. 論文: http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla12eccv.pdf
2. プロジェクトホームページ：http://www.robesafe.com/personal/pablo.alcantarilla/kaze.html
3. 作者コード：http://www.robesafe.com/personal/pablo.alcantarilla/code/kaze_フィーチャー1_4.(bootライブラリが必要で、また、その計時関数の使用は複雑で、OpenCVのcv:getTickCountで代替できます。)
4. Computter Vision Talksの評価：http://computer-vision-talks.com/2013/03/porting-kaze-features-to-opencv/
5. Computter Vision TalksブロガーIevgen KhvedcheniaはKAZEをOpenCVのcvに統合します。Feature 2 Dクラスですが、OpenCVを再コンパイルする必要があります。また、アルゴリズムパラメータの調整とMaskフィルタの特徴点を実現する機能がありません。https://github.com/BloodAxe/opencv/tree/kaze-features
6. IevgenのプロジェクトライブラリからKAZEを抽出し、cvを継承するようにパッケージ化しました。Feature 2 Dのクラスは、OpenCVを再コンパイルする必要がなく、パラメータ調整とMaskフィルタリングの機能を実現しました。https://github.com/yuhuazou/kaze_opencv
7. Matlab版のインターフェースプログラムは、1.0版のKAZEコードを実装しています。https://github.com/vlfeat/vlbenchmarks/blob/unstable/%2BlocalFeatures/Kaze.m
概要
ECCV 2012にはSIFTよりも安定した特徴検出アルゴリズムKAZE（［1］）が現れた。KAZEの名前は尺度空間分析の創始者である日本の学者イジマを記念したものです。KAZEは日本語の「風」の語呂合わせで、風の形成が空間における非線形の流動過程であるかのように、KAZE特徴検出は画像領域における非線形拡散処理の過程であることを意味する。
従来のSIFT、SURFなどの特徴検出アルゴリズムは、線形Gaussピラミッドに基づいてマルチスケール分解し、ノイズを除去し、顕著な特徴点を抽出する。しかしGauss分解は局所精度を犠牲にしたものであり，境界の曖昧さと詳細の損失をもたらしやすい。非線形スケール分解はこの問題を解決するために有望であるが、従来の方法は正のオーロラ法に基づいて非線形拡散（Non-linear diffusion）方程式を解く際に反復収束するステップが短すぎて、時間がかかり、計算が複雑度が高い。これにより、KAZEアルゴリズムの著者らは、加性演算子分割アルゴリズム（Additive Operator Spliting，AOS）を用いた非線形拡散フィルタリングを提案し、任意のステップサイズで安定した非線形スケール空間を構築することができる。
注：KAZEアルゴリズムの原理はSIFTとSURFとよく似ています。KAZEを深く知る前に、以下のブログ記事を参考にしてSIFTとSURFについて紹介します。
1. 【OpenCV】SIFT原理とソース分析-魏さんの修行路（いいブログです。心を込めて、真面目で、文脈も豊かです。）
2. Surfアルゴリズム学習心得(一)——アルゴリズム原理-yang_yong'blog
3. Opencv学習ノート（六）SURF学習ノート-crazy_sparrowのコラム
1.1非線形拡散フィルタリング
1.1.1 Perona-Maalik拡散方程式
具体的には、非線形拡散フィルタリング方法は、画像輝度（L）の異なるスケールでの変化をある形式の流動関数（flow function）の発散と見なすものであり、非線形偏微分方程式により記述することができる。
適切な伝導関数c（x，y，t）を設定することで，画像の局所構造に拡散を適応させることができる。伝導関数はスカラーでも良いし、テンソルでもいいです。時間tはスケールパラメータとして、その値が大きいほど画像の表現形式が簡単になる。PeronaとMalikは伝導関数の構造を提案した。
そのうちの▽Lσガウス平滑化後の画像Lです。σの勾配です。関数g()の形は以下の通りです。
ここで、関数g 1は高コントラストのエッジを優先的に保持し、g 2は幅の大きい領域を優先的に保持し、g 3は領域内部を効果的に平滑化して境界情報を保持することができる（KAZEコードではデフォルトでは関数g 2を使用する）。
関数g 1、g 2、g 3の実現コードは以下の通りです。nldifusion_functions.cpp中)：
//*************************************************************************************
//*************************************************************************************
/**
* @brief This function computes the Perona and Malik conductivity coefficient g1
* g1 = exp(-|dL|^2/k^2)
* @param src Input image
* @param dst Output image
* @param Lx First order image derivative in X-direction (horizontal)
* @param Ly First order image derivative in Y-direction (vertical)
* @param k Contrast factor parameter
*/
void PM_G1(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
{
cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly))/(k*k),dst);
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This function computes the Perona and Malik conductivity coefficient g2
* g2 = 1 / (1 + dL^2 / k^2)
* @param src Input image
* @param dst Output image
* @param Lx First order image derivative in X-direction (horizontal)
* @param Ly First order image derivative in Y-direction (vertical)
* @param k Contrast factor parameter
*/
void PM_G2(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
{
dst = 1./(1. + (Lx.mul(Lx) + Ly.mul(Ly))/(k*k));
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This function computes Weickert conductivity coefficient g3
* @param src Input image
* @param dst Output image
* @param Lx First order image derivative in X-direction (horizontal)
* @param Ly First order image derivative in Y-direction (vertical)
* @param k Contrast factor parameter
* @note For more information check the following paper: J. Weickert
* Applications of nonlinear diffusion in image processing and computer vision,
* Proceedings of Algorithmy 2000
*/
void Weickert_Diffusivity(const cv::Mat &src, cv::Mat &dst, cv::Mat &Lx, cv::Mat &Ly, float k )
{
cv::Mat modg;
cv::pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg);
cv::exp(-3.315/modg, dst);
dst = 1.0 - dst;
}
パラメータkは拡散レベルを制御するコントラスト因子であり、どれぐらいのエッジ情報を保持するかを決定することができ、その値が大きいほど、保持するエッジ情報が少ない。KAZEアルゴリズムでは、パラメータkの値は勾配画像▽Lである。σのヒストグラムの70%のビットの値:
計算パラメータkのソースコードは以下の通りです。nldifusion_functions.cpp中)：
//*************************************************************************************
//*************************************************************************************
/**
* @brief This function computes a good empirical value for the k contrast factor
* given an input image, the percentile (0-1), the gradient scale and the number of
* bins in the histogram
* @param img Input image
* @param perc Percentile of the image gradient histogram (0-1)
* @param gscale Scale for computing the image gradient histogram
* @param nbins Number of histogram bins
* @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel
* @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel
* @return k contrast factor
*/
float Compute_K_Percentile(const cv::Mat &img, float perc, float gscale, unsigned int nbins, unsigned int ksize_x, unsigned int ksize_y)
{
float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0;
unsigned int nbin = 0, nelements = 0, nthreshold = 0, k = 0;
float npoints = 0.0; // number of points of which gradient greater than zero
float hmax = 0.0; // maximum gradient
// Create the array for the histogram
float *hist = new float[nbins];
// Create the matrices
cv::Mat gaussian = cv::Mat::zeros(img.rows,img.cols,CV_32F);
cv::Mat Lx = cv::Mat::zeros(img.rows,img.cols,CV_32F);
cv::Mat Ly = cv::Mat::zeros(img.rows,img.cols,CV_32F);
// Set the histogram to zero, just in case
for( unsigned int i = 0; i < nbins; i++ )
{
hist[i] = 0.0;
}
// Perform the Gaussian convolution
Gaussian_2D_Convolution(img,gaussian,ksize_x,ksize_y,gscale);
// Compute the Gaussian derivatives Lx and Ly
Image_Derivatives_Scharr(gaussian,Lx,1,0);
Image_Derivatives_Scharr(gaussian,Ly,0,1);
// Skip the borders for computing the histogram
for( int i = 1; i < gaussian.rows-1; i++ )
{
for( int j = 1; j < gaussian.cols-1; j++ )
{
lx = *(Lx.ptr<float>(i)+j);
ly = *(Ly.ptr<float>(i)+j);
modg = sqrt(lx*lx + ly*ly);
// Get the maximum
if( modg > hmax )
{
hmax = modg;
}
}
}
// Skip the borders for computing the histogram
for( int i = 1; i < gaussian.rows-1; i++ )
{
for( int j = 1; j < gaussian.cols-1; j++ )
{
lx = *(Lx.ptr<float>(i)+j);
ly = *(Ly.ptr<float>(i)+j);
modg = sqrt(lx*lx + ly*ly); // modg can be stored in a matrix in last step in oder to avoid re-computation
// Find the correspondent bin
if( modg != 0.0 )
{
nbin = floor(nbins*(modg/hmax));
if( nbin == nbins )
{
nbin--;
}
hist[nbin]++;
npoints++;
}
}
}
// Now find the perc of the histogram percentile
nthreshold = (unsigned int)(npoints*perc);
// find the bin (k) in which accumulated points are greater than 70% (perc) of total valid points (npoints)
for( k = 0; nelements < nthreshold && k < nbins; k++)
{
nelements = nelements + hist[k];
}
if( nelements < nthreshold )
{
kperc = 0.03;
}
else
{
kperc = hmax*((float)(k)/(float)nbins);
}
delete hist;
return kperc;
}
注：非線形拡散フィルタの適用については、［2］を参照してください。
1.1.2 AOSアルゴリズム
非線形偏微分方程式は解析解がないので，一般的に数値解析の方法で反復的に解いた。従来，明示的差分形式を用いた解法は，小さなステップだけを用いて収束が緩やかである。このために、方程式を以下のように離散化しました。
ここでAlは、各次元（l）に画像をアップロードするガイドを示すマトリクスである。この方程式の解は以下の通りである。
この解法は任意の時間ステップに対して（τ）すべて有効です。上式マトリックスAlは三対角行列で対角優位（tridiagonal and diagonally dominant marix）であり、このような線形システムはThomasアルゴリズムによって迅速に解くことができる。AOSの応用については、［3］を参照してください。
アルゴリズムの実現ソースは以下の通りです。
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs a scalar non-linear diffusion step using AOS schemes
* @param Ld Image at a given evolution step
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
* @note If c is constant, the diffusion will be linear
* If c is a matrix of the same size as Ld, the diffusion will be nonlinear
* The stepsize can be arbitrarilly large
*/
void KAZE::AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
{
AOS_Rows(Ldprev,c,stepsize);
AOS_Columns(Ldprev,c,stepsize);
Ld = 0.5*(Lty + Ltx.t());
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs a scalar non-linear diffusion step using AOS schemes
* Diffusion in each dimension is computed independently in a different thread
* @param Ld Image at a given evolution step
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
* @note If c is constant, the diffusion will be linear
* If c is a matrix of the same size as Ld, the diffusion will be nonlinear
* The stepsize can be arbitrarilly large
*/
#if HAVE_THREADING_SUPPORT
void KAZE::AOS_Step_Scalar_Parallel(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
{
boost::thread *AOSth1 = new boost::thread(&KAZE::AOS_Rows,this,Ldprev,c,stepsize);
boost::thread *AOSth2 = new boost::thread(&KAZE::AOS_Columns,this,Ldprev,c,stepsize);
AOSth1->join();
AOSth2->join();
Ld = 0.5*(Lty + Ltx.t());
delete AOSth1;
delete AOSth2;
}
#endif
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs performs 1D-AOS for the image rows
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
*/
void KAZE::AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
{
// Operate on rows
for( int i = 0; i < qr.rows; i++ )
{
for( int j = 0; j < qr.cols; j++ )
{
*(qr.ptr<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i+1)+j);
}
}
for( int j = 0; j < py.cols; j++ )
{
*(py.ptr<float>(0)+j) = *(qr.ptr<float>(0)+j);
}
for( int j = 0; j < py.cols; j++ )
{
*(py.ptr<float>(py.rows-1)+j) = *(qr.ptr<float>(qr.rows-1)+j);
}
for( int i = 1; i < py.rows-1; i++ )
{
for( int j = 0; j < py.cols; j++ )
{
*(py.ptr<float>(i)+j) = *(qr.ptr<float>(i-1)+j) + *(qr.ptr<float>(i)+j);
}
}
// a = 1 + t.*p; (p is -1*p)
// b = -t.*q;
ay = 1.0 + stepsize*py; // p is -1*p
by = -stepsize*qr;
// Call to Thomas algorithm now
Thomas(ay,by,Ldprev,Lty);
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs performs 1D-AOS for the image columns
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
*/
void KAZE::AOS_Columns(const cv::Mat &Ldprev, const cv::Mat &c, const float stepsize)
{
// Operate on columns
for( int j = 0; j < qc.cols; j++ )
{
for( int i = 0; i < qc.rows; i++ )
{
*(qc.ptr<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i)+j+1);
}
}
for( int i = 0; i < px.rows; i++ )
{
*(px.ptr<float>(i)) = *(qc.ptr<float>(i));
}
for( int i = 0; i < px.rows; i++ )
{
*(px.ptr<float>(i)+px.cols-1) = *(qc.ptr<float>(i)+qc.cols-1);
}
for( int j = 1; j < px.cols-1; j++ )
{
for( int i = 0; i < px.rows; i++ )
{
*(px.ptr<float>(i)+j) = *(qc.ptr<float>(i)+j-1) + *(qc.ptr<float>(i)+j);
}
}
// a = 1 + t.*p';
ax = 1.0 + stepsize*px.t();
// b = -t.*q';
bx = -stepsize*qc.t();
// Call Thomas algorithm again
// But take care since we need to transpose the solution!!
Thomas(ax,bx,Ldprev.t(),Ltx);
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method does the Thomas algorithm for solving a tridiagonal linear system
* @note The matrix A must be strictly diagonally dominant for a stable solution
*/
void KAZE::Thomas(cv::Mat a, cv::Mat b, cv::Mat Ld, cv::Mat x)
{
// Auxiliary variables
int n = a.rows;
cv::Mat m = cv::Mat::zeros(a.rows,a.cols,CV_32F);
cv::Mat l = cv::Mat::zeros(b.rows,b.cols,CV_32F);
cv::Mat y = cv::Mat::zeros(Ld.rows,Ld.cols,CV_32F);
/** A*x = d; */
/** / a1 b1 0 0 0 ... 0 \ / x1 \ = / d1 \ */
/** | c1 a2 b2 0 0 ... 0 | | x2 | = | d2 | */
/** | 0 c2 a3 b3 0 ... 0 | | x3 | = | d3 | */
/** | : : : : 0 ... 0 | | : | = | : | */
/** | : : : : 0 cn-1 an | | xn | = | dn | */
/** 1. LU decomposition
/ L = / 1 \ U = / m1 r1 \
/ | l1 1 | | m2 r2 |
/ | l2 1 | | m3 r3 |
/ | : : : | | : : : |
/ \ ln-1 1 / \ mn / */
for( int j = 0; j < m.cols; j++ )
{
*(m.ptr<float>(0)+j) = *(a.ptr<float>(0)+j);
}
for( int j = 0; j < y.cols; j++ )
{
*(y.ptr<float>(0)+j) = *(Ld.ptr<float>(0)+j);
}
// 2. Forward substitution L*y = d for y
for( int k = 1; k < n; k++ )
{
for( int j=0; j < l.cols; j++ )
{
*(l.ptr<float>(k-1)+j) = *(b.ptr<float>(k-1)+j) / *(m.ptr<float>(k-1)+j);
}
for( int j=0; j < m.cols; j++ )
{
*(m.ptr<float>(k)+j) = *(a.ptr<float>(k)+j) - *(l.ptr<float>(k-1)+j)*(*(b.ptr<float>(k-1)+j));
}
for( int j=0; j < y.cols; j++ )
{
*(y.ptr<float>(k)+j) = *(Ld.ptr<float>(k)+j) - *(l.ptr<float>(k-1)+j)*(*(y.ptr<float>(k-1)+j));
}
}
// 3. Backward substitution U*x = y
for( int j=0; j < y.cols; j++ )
{
*(x.ptr<float>(n-1)+j) = (*(y.ptr<float>(n-1)+j))/(*(m.ptr<float>(n-1)+j));
}
for( int i = n-2; i >= 0; i-- )
{
for( int j = 0; j < x.cols; j++ )
{
*(x.ptr<float>(i)+j) = (*(y.ptr<float>(i)+j) - (*(b.ptr<float>(i)+j))*(*(x.ptr<float>(i+1)+j)))/(*(m.ptr<float>(i)+j));
}
}
}
引き続き…
Ref:
[1]http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla12eccv.pdf
[2]http://manu16.magtech.com.cn/geoprog/CN/article/downloadArticleFile.do?attachType=PDF&id=3146
[3]http://file.lw23.com/8/8e/8ec/8ecd21e4-b030-4e05-9333-40cc2d97bde4.pdf