圖像特征提取：圖像的矩特征

1. 矩的概念

圖像識別的一個核心問題是圖像的特征提取，簡單描述即為用一組簡單的數據（圖像描述量）來描述整個圖像，這組數據越簡單越有代表性越好。良好的特征不受光線、噪點、幾何形變的干擾。圖像識別發展幾十年，不斷有新的特征提出，而圖像不變矩就是其中一個。

矩是概率與統計中的一個概念，是隨機變量的一種數字特征。設X為隨機變量，c為常數，k為正整數。則量E[(x−c)k]稱為X關於c點的k階矩。

比較重要的有兩種情況：

1. c=0。這時ak=E(Xk)稱為X的k階原點矩

2. c=E(X)。這時μk=E[(X−EX)k]稱為X的k階中心矩。

一階原點矩就是期望。一階中心矩μ1=0，二階中心矩μ2就是X的方差Var(X)。在統計學上，高於4階的矩極少使用。μ3可以去衡量分布是否有偏。μ4可以去衡量分布（密度）在均值附近的陡峭程度如何。

針對於一幅圖像，我們把像素的坐標看成是一個二維隨機變量(X,Y)，那么一幅灰度圖像可以用二維灰度密度函數來表示，因此可以用矩來描述灰度圖像的特征。

不變矩(Invariant Moments)是一處高度濃縮的圖像特征，具有平移、灰度、尺度、旋轉不變性。M.K.Hu在1961年首先提出了不變矩的概念。1979年M.R.Teague根據正交多項式理論提出了Zernike矩。下面主要介紹這兩種矩特征的算法原理與實現。

2. Hu矩

一幅M×N的數字圖像f(i,j)，其p+q階幾何矩mpq和中心矩μpq為：

 

mpq=∑i=1M∑j=1Nipjqf(i,j)

 

 

μpq=∑i=1M∑j=1N(i−i¯)p(j−j¯)qf(i,j)

 

其中f(i,j)為圖像在坐標點(i,j)處的灰度值。i¯=m10/m00,j¯=m01/m00

若將m00看作是圖像的灰度質量，則(i¯,j¯)為圖像的質心坐標，那么中心矩μpa反映的是圖像灰度相對於其灰度質心的分布情況。可以用幾何矩來表示中心矩，0~3階中心矩與幾何矩的關系如下：

μ00=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)0f(i,j)=m00

μ10=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)0f(i,j)=0

μ01=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)1f(i,j)=0

μ11=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)1f(i,j)=m11−y¯m10

μ20=∑Mi=1∑Nj=1(i−i¯)2(j−j¯)0f(i,j)=m20−y¯m01

μ02=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)2f(i,j)=m02−y¯m01

μ30=∑Mi=1∑Nj=1(i−i¯)3(j−j¯)0f(i,j)=m30−2x¯m20+2x¯2m10

μ12=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)2f(i,j)=m12−2y¯m11−x¯m02+2y¯2m10

μ21=∑Mi=1∑Nj=1(i−i¯)2(j−j¯)1f(i,j)=m21−2x¯m11−y¯m20+2x¯2m01

μ03=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)3f(i,j)=m03−2y¯m02+2y¯2m01

為了消除圖像比例變化帶來的影響，定義規格化中心矩如下：

 

ηpq=μpaμγ00,(γ=p+q2,p+q=2,3,…)

 

利用二階和三階規格中心矩可以導出下面7個不變矩組(Φ1 Φ7)，它們在圖像平移、旋轉和比例變化時保持不變。

Φ1=η20+η02

Φ2=(η20−η02)2+4η211

Φ3=(η20−3η12)2+3(η21−η03)2

Φ4=(η30+η12)2+(η21+η03)2

Φ5=(η30+3η12)(η30+η12)[(η30+η12)2−3(η21+η03)2]+(3η21−η03)(η21+η03)[3(η30+η12)2−(η21+η03)2]

Φ6=(η20−η02)[(η30+η12)2−(η21+η03)2]+4η11(η30+η12)(η21+η03)

Φ7=(3η21−η03)(η30+η12)[(η30+η12)2−3(η21+η03)2]+]+(3η12−η30)(η21+η03)[3(η30+η12)2−(η21+η03)2]

3. 利用OpenCV計算Hu矩

opencv里對Hu矩的計算有直接的API，它分為了兩個函數：moments()函數用於計算中心矩，HuMoments函數用於由中心矩計算Hu矩。

Moments moments(InputArray array, bool binaryImage=false )

參數說明

• 輸入參數：array是一幅單通道，8-bits的圖像，或一個二維浮點數組(Point of Point2f)。binaryImage用來指示輸出圖像是否為一幅二值圖像，如果是二值圖像，則圖像中所有非0像素看作為1進行計算。
• 輸出參數：moments是一個類：

class Moments { public: Moments(); Moments(double m00, double m10, double m01, double m20, double m11, double m02, double m30, double m21, double m12, double m03 ); Moments( const CvMoments& moments ); operator CvMoments() const;

里面保存了圖像的2階與3階中心矩的值。

void HuMoments(const Moments& moments, double* hu)

參數說明：

• 輸入參數：moments即為上面一個函數計算得到的moments類型。
• 輸出參數：hu是一個含有7個數的數組。

int main(int argc, char** argv) 
{ 
    Mat image = imread(argv[1]);  
    cvtColor(image, image, CV_BGR2GRAY); 
    Moments mts = moments(image); 
    double hu[7]; 
    HuMoments(mts, hu); 
    for (int i=0; i<7; i++) 
    { 
        cout << log(abs(hu[i])) <<endl; 
    } 
   return 0; 
}

上面代碼中，最終輸出的值為log|Φi|

我們分別計算一幅圖像在，旋轉，噪聲與模糊時的Hu矩。

   

類別	log|Φ1|	log|Φ2|	log|Φ3|	log|Φ4|	log|Φ5|	log|Φ6|	log|Φ7|
原圖	-6.76181	-19.1286	-23.7441	-26.776	-51.7618	-35.8491	-51.534
旋轉	-6.72102	-19.0844	-23.5756	-25.9122	-51.4619	-35.4595	-50.7674
加放噪點	-6.76086	-19.1255	-23.7611	-26.3228	-51.5056	-35.895	-51.6321
模糊	-6.76183	-19.1295	-23.7451	-26.2767	-51.765	-35.8484	-51.5307

4. Zernike矩

Hu矩在圖像描述上有廣泛的應用，但是其低階幾何矩與圖像整體特征有關，不包含太多的圖像細節信息，而高階幾何矩易受噪聲影響，因此很難利用幾何矩恢復圖像。

Zernike矩能夠很容易地構造圖像的任意高階矩，並能夠使用較少的矩來重建圖像。Zernike矩是基於Zernike多項式的正交化函數，雖然其計算比較復雜，但是Zernide矩在圖像旋轉和低噪聲敏感度方面具有較大的優越性。由於Zernike矩具有圖像旋轉不變性，而且可以構造任意高階矩，所以被廣泛應用對目標進行識別中。

4.1 Zernike矩多項式

首先要弄清楚什么是正交多項式。若函數W(x)在區間(a,b)可積，且W(x)≥0，則可作為權函數。

對於一個多項式的序列fi和權函數W(x)，定義內積：<fm,fn>=∫bafm(x)fn(x)W(x)dx

若n≠m,<fm,fn>=0，這些多項式則稱為正交多項式。若fi除了正交之外，更有<fm,fn>=1的話，則稱為規范正交多項式。

那么正交多項式有什么作用呢？答案是：逼近！正交多項式相當於基，任何一個n維多項式函數f(x)都可以用一組正交多項式加權求和來逼近。

 

Zernike在1934年提出了在單位圓上定義的一組正交多項式，即Zernike正交多項式，其定義形式為：

 

Rnm(ρ)=∑s=0(n−|m|)/2(−1)s[(n−s)!]ρn−2ss!(n+|m|2−s)!(n+|m|2+s)!

 

 

Vnm(x,y)=Vnm(ρ,θ)=Rnm(ρ)ejmθ

 

其中Rnm(ρ)表示點(x,y)的徑向多項式，Vnm(x,y)為Zernike正交多項式，n,m為正交多項式的階數，n是非負整數，n−|m|是偶數，並且n≥|m|。

Zernike多項式Vnm(x,y)=Vnm(ρ,θ)是定義在單位圓x2+y2≤1上的正交復函數的集合，具有重要的遞推性質，即Rnm可由R(n−2)m和R(n−4)m得到，公式如下：

 

Rnm(ρ)=[(K22ρ2+K3)R(n−2)m(ρ)+K4R(n−4)m(ρ)]K1

 

 

Rmm(ρ)=ρm

 

式中：K1=(n+1)(n−1)(n−2)/2,K2=2n(n−1)(n−2),K3=−(n−1)3,K4=−n(n−1)(n−3)/2。

4.2 Zernike矩的定義

由於Zernike多項式的正交完備性，所以在單位圓內的任何圖像f(x,y)都可以唯一的用下面式子展開：

 

f(x,y)=∑n=0∞∑m=0∞ZnmVn,m(ρ,θ)

 

上式中的Znm就是Zernike矩。

對二維函數f(x,y)的Zernike矩的定義如下：

 

Znm=n+1π∫10∫2π0[Vnm(ρ,θ)]f(ρ,θ)ρdydxdρdθ

 

 

=n+1π∬Rnm(ρ)ejmθf(ρ,θ)dρdθ

 

式中ρ=x2+y2−−−−−−√(−1<x,y<1)，θ為軸x與ρ矢量在逆時針方向的夾角；Rnm(ρ)表示點(x,y)的徑向多項式。

4.3 Zernike矩的計算

從Zernike矩的計算公式上來看，對於二維圖像，其Zernike矩Znm為復數，將其實部和虛部分別記為Cnm和Snm，則有：

 

Cnm=2n+2π∫10∫2π0[Rnm(ρ)cos(mθ)f(ρ,θ)ρdρdθ

 

 

Cnm=2n+2π∫10∫2π0[Rnm(ρ)sin(mθ)f(ρ,θ)ρdρdθ

 

因為數字圖像是離散形式的點，所以需要將上式離散化，把積分號換為求和號，但是需要作一些坐標變換。

對於N×N的圖像f(x,y)，令坐標原點位於圖像的中心，則−N/2≤x,y≤N/2，對於像素(x,y)，引入2個參數(r,σ)，唯一對應於像素，其定義為：

r=max(|x|,|y|)

如果|x|=r，則：

 

σ=2(r−x)y|y|+xyr

 

如果|y|=r，則：

 

σ=2y−xyr

 

我們容易計算出，r的取值范圍為1∼N/2，σ的取值范圍是1∼8r，再根據參數(r,σ)定義相應的極坐標：

 

ρ=2r/N,θ=πσ(4r)

 

所以，最終我們得到離散化的Zernike矩的計算公式：

 

Cnm=2n+2N2∑r=1N/2Rnm(2r/N)∑σ=18rcosπmσ4rf(r,σ)

 

 

Snm=2n+2N2∑r=1N/2Rnm(2r/N)∑σ=18rsinπmσ4rf(r,σ)

 

1. 確定圖像的大小N×N，即公式中的N；

2. 確定r,σ的范圍；

3. 利用Zernike多項式的遞推性質計算各階Rnm(ρ)，並結合上面Zernike矩計算公式，算出Cnm,Snm

4. 對Cnm,Snm求模，進而計算得到|Znm|

現在我們用Zernike矩來計算美女圖像在4種狀態下的值：

類別	log|Z11|	log|Z20|	log|Z22|	log|Z31|	log|Z40|	log|Z42|	log|Z44|
原圖	11.1732	13.8469	12.3515	12.4391	14.2782	12.6137	11.5745
旋轉	12.3036	13.8309	13.5861	12.0467	13.1320	13.8396	12.7862
加放噪點	11.1538	13.8490	12.3315	12.4316	14.2730	12.5925	11.5591
模糊	11.1636	13.8465	12.3480	12.4367	14.2799	12.6130	11.5752

通過表中，可以看出，Zernike在總體上效果比Hu矩更好（PS:感覺在旋轉上好像差強人意！）

下面是Zernike矩的matlab實現[來自《現代數字圖像-處理技術提高及應用案例詳解》]，這里偷懶了，有機會的話會把C++版的實現補上。

 

5. 總結

不變矩的應用過程一般包括：

1. 選擇合適的不變矩類型；
2. 選擇分類器（如神經網絡、最短距離等）；
3. 如果是神經網絡分類器，則需要計算學習樣例的不變矩去訓練神經網絡；
4. 計算待識別對象的不變矩，輸入神經網絡就可得到待識別對象的類型，或者計算待識別對象不變矩與類別對象不變矩之間的距離，選擇最短距離的類別作為待識別對象的類別。

可以看出，不變矩作用主要目的是描述事物（圖像）的特征。人眼識別圖像的特征往往又表現為“求和”的形式，因此不變矩是對圖像元素進行了積分操作。

不變矩能夠描述圖像整體特征就是因為它具有平移不變形、比例不變性和旋轉不變性等性質。

然而，另一方面圖像的各階不變矩究竟代表的什么特征很難進行直觀的物理解釋。

6. 參考資料

[1] 《現代數字圖像處理》（matlab版）

[2] 正交多項式WIKI

[3] opencv形態描述