1. 矩的概念
圖像識別的一個核心問題是圖像的特征提取，簡單描述即為用一組簡單的數據（圖像描述量）來描述整個圖像，這組數據越簡單越有代表性越好。良好的特征不受光線、噪點、幾何形變的干擾。圖像識別發展幾十年，不斷有新的特征提出，而圖像不變矩就是其中一個。
矩是概率與統計中的一個概念，是隨機變量的一種數字特征。設$X$為隨機變量，$c$為常數，$k$為正整數。則量$E[(x-c)^k]$稱為$X$關於$c$點的$k$階矩。
比較重要的有兩種情況：
1. $c=0$。這時$a_k=E(X^k)$稱為$X$的$k$階原點矩
2. $c=E(X)$。這時$\mu_k=E[(X-EX)^k]$稱為$X$的$k$階中心矩。
一階原點矩就是期望。一階中心矩$\mu_1=0$，二階中心矩$\mu_2$就是$X$的方差$Var(X)$。在統計學上，高於4階的矩極少使用。$\mu_3$可以去衡量分布是否有偏。$\mu_4$可以去衡量分布（密度）在均值附近的陡峭程度如何。
針對於一幅圖像，我們把像素的坐標看成是一個二維隨機變量$(X,Y)$，那么一幅灰度圖像可以用二維灰度密度函數來表示，因此可以用矩來描述灰度圖像的特征。
不變矩(Invariant Moments)是一處高度濃縮的圖像特征，具有平移、灰度、尺度、旋轉不變性。M.K.Hu在1961年首先提出了不變矩的概念。1979年M.R.Teague根據正交多項式理論提出了Zernike矩。下面主要介紹這兩種矩特征的算法原理與實現。
2. Hu矩
一幅$M\times N$的數字圖像$f(i,j)$，其$p+q$階幾何矩$m_{pq}$和中心矩$\mu_{pq}$為：
$$m_{pq}=\sum_{i=1}^M\sum_{j=1}^Ni^pj^qf(i,j)$$
$$\mu_{pq}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^p(j-\bar{j})^qf(i,j)$$
其中$f(i,j)$為圖像在坐標點$(i,j)$處的灰度值。$\bar{i}=m_{10}/m_{00},\bar{j}=m_{01}/m_{00}$
若將$m_{00}$看作是圖像的灰度質量，則$(\bar{i},\bar{j})$為圖像的質心坐標，那么中心矩$\mu_{pa}$反映的是圖像灰度相對於其灰度質心的分布情況。可以用幾何矩來表示中心矩，0~3階中心矩與幾何矩的關系如下：
$\mu{00}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^0f(i,j)=m_{00}$
$\mu{10}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^0f(i,j)=0$
$\mu{01}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^1f(i,j)=0$
$\mu{11}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^1f(i,j)=m_{11}-\bar{y}m_{10}$
$\mu{20}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^2(j-\bar{j})^0f(i,j)=m_{20}-\bar{y}m_{01}$
$\mu{02}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^2f(i,j)=m_{02}-\bar{y}m_{01}$
$\mu{30}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^3(j-\bar{j})^0f(i,j)=m_{30}-2\bar{x}m_{20}+2\bar{x}^2m_{10}$
$\mu{12}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^2f(i,j)=m_{12}-2\bar{y}m_{11}-\bar{x}m_{02}+2\bar{y}^2m_{10}$
$\mu{21}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^2(j-\bar{j})^1f(i,j)=m_{21}-2\bar{x}m_{11}-\bar{y}m_{20}+2\bar{x}^2m_{01}$
$\mu{03}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^3f(i,j)=m_{03}-2\bar{y}m_{02}+2\bar{y}^2m_{01}$
為了消除圖像比例變化帶來的影響，定義規格化中心矩如下：
$$\eta_{pq}=\frac{\mu_{pa}}{\mu_{00}^\gamma},(\gamma=\frac{p+q}{2},p+q=2,3,\dots)$$
利用二階和三階規格中心矩可以導出下面7個不變矩組$(\Phi_1~\Phi_7)$，它們在圖像平移、旋轉和比例變化時保持不變。
$\Phi_1=\eta_{20}+\eta_{02}$
$\Phi_2=(\eta_{20}-\eta_{02})^2+4\eta_{11}^2$
$\Phi_3=(\eta_{20}-3\eta_{12})^2+3(\eta_{21}-\eta_{03})^2$
$\Phi_4=(\eta_{30}+\eta_{12})^2+(\eta_{21}+\eta_{03})^2$
$\Phi_5=(\eta_{30}+3\eta_{12})(\eta_{30}+\eta_{12})[(\eta_{30}+\eta_{12})^2-3(\eta_{21}+\eta_{03})^2]+(3\eta_{21}-\eta_{03})(\eta_{21}+\eta_{03})[3(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{03})^2]$
$\Phi_6=(\eta_{20}-\eta_{02})[(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{03})^2]+4\eta_{11}(\eta_{30}+\eta_{12})(\eta_{21}+\eta_{03})$
$\Phi_7=(3\eta_{21}-\eta_{03})(\eta_{30}+\eta_{12})[(\eta_{30}+\eta_{12})^2-3(\eta_{21}+\eta_{03})^2]+]+(3\eta_{12}-\eta_{30})(\eta_{21}+\eta_{03})[3(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{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|\Phi_i|$
我們分別計算一幅圖像在，旋轉，噪聲與模糊時的Hu矩。
| 類別 | $log|\Phi_1|$ | $log|\Phi_2|$ | $log|\Phi_3|$ | $log|\Phi_4|$ | $log|\Phi_5|$ | $log|\Phi_6|$ | $log|\Phi_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)\ge 0$，則可作為權函數。
對於一個多項式的序列$f_i$和權函數$W(x)$，定義內積：$<f_m,f_n>=\int_{a}^{b}f_m(x)f_n(x)W(x)dx$
若$n \ne m,<f_m,f_n>=0$，這些多項式則稱為正交多項式。若$f_i$除了正交之外，更有$<f_m,f_n>=1$的話，則稱為規范正交多項式。
那么正交多項式有什么作用呢？答案是：逼近！正交多項式相當於基，任何一個n維多項式函數$f(x)$都可以用一組正交多項式加權求和來逼近。
Zernike在1934年提出了在單位圓上定義的一組正交多項式，即Zernike正交多項式，其定義形式為：
$$R_{nm}(\rho) = \sum_{s=0}^{(n-|m|)/2}\frac{(-1)^s[(n-s)!]\rho^{n-2s}}{s!(\frac{n+|m|}{2}-s)!(\frac{n+|m|}{2}+s)!}$$
$$V_{nm}(x,y)=V_{nm}(\rho,\theta)=R_{nm}(\rho)e^{jm\theta}$$
其中$R_{nm}(\rho)$表示點$(x,y)$的徑向多項式，$V_{nm}(x,y)$為Zernike正交多項式，$n,m$為正交多項式的階數，$n$是非負整數，$n-|m|$是偶數，並且$n\ge|m|$。
Zernike多項式$V_{nm}(x,y)=V_{nm}(\rho,\theta)$是定義在單位圓$x^2+y^2\le 1$上的正交復函數的集合，具有重要的遞推性質，即$R_{nm}$可由$R_{(n-2)m}$和$R_{(n-4)m}$得到，公式如下：
$$R_{nm}(\rho)=\frac{[(K_2^2\rho^2+K_3)R_{(n-2)m}(\rho)+K_4R_{(n-4)m}(\rho)]}{K_1}$$
$$R_{mm}(\rho)=\rho^m$$
式中：$K_1=(n+1)(n-1)(n-2)/2,K_2=2n(n-1)(n-2),K_3=-(n-1)^3,K_4=-n(n-1)(n-3)/2$。
4.2 Zernike矩的定義
由於Zernike多項式的正交完備性，所以在單位圓內的任何圖像$f(x,y)$都可以唯一的用下面式子展開：
$$f(x,y)=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}Z_{nm}V_{n,m}(\rho,\theta)$$
上式中的$Z_{nm}$就是Zernike矩。
對二維函數$f(x,y)$的Zernike矩的定義如下：
$$Z_{nm}=\frac{n+1}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[V_{nm}(\rho,\theta)]f(\rho,\theta)\rho\frac{dy}{dx}d\rho d\theta$$
$$=\frac{n+1}{\pi}\iint R_{nm}(\rho)e^{jm\theta}f(\rho,\theta)d\rho d\theta$$
式中$\rho=\sqrt{x^2+y^2}(-1<x,y<1)$，$\theta$為軸$x$與$\rho$矢量在逆時針方向的夾角；$R_{nm}(\rho)$表示點$(x,y)$的徑向多項式。
4.3 Zernike矩的計算
從Zernike矩的計算公式上來看，對於二維圖像，其Zernike矩$Z_{nm}$為復數，將其實部和虛部分別記為$C_{nm}$和$S_{nm}$，則有：
$$C_{nm}=\frac{2n+2}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[R_{nm}(\rho)cos(m\theta)f(\rho,\theta)\rho d\rho d\theta$$
$$C_{nm}=\frac{2n+2}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[R_{nm}(\rho)sin(m\theta)f(\rho,\theta)\rho d\rho d\theta$$
因為數字圖像是離散形式的點，所以需要將上式離散化，把積分號換為求和號，但是需要作一些坐標變換。
對於$N\times N$的圖像$f(x,y)$，令坐標原點位於圖像的中心，則$-N/2\le x,y\le N/2$，對於像素$(x,y)$，引入2個參數$(r,\sigma)$，唯一對應於像素，其定義為：
$r=max(|x|,|y|)$
如果$|x|=r$，則：
$$\sigma=\frac{2(r-x)y}{|y|}+\frac{xy}{r}$$
如果$|y|=r$，則：
$$\sigma=2y-\frac{xy}{r}$$
我們容易計算出，$r$的取值范圍為$1\sim N/2$，$\sigma$的取值范圍是$1\sim 8r$，再根據參數$(r,\sigma)$定義相應的極坐標：
$$\rho=2r/N,\theta=\pi\sigma(4r)$$
所以，最終我們得到離散化的Zernike矩的計算公式：
$$C_{nm}=\frac{2n+2}{N^2}\sum_{r=1}^{N/2}R_{nm}(2r/N)\sum_{\sigma=1}^{8r}cos\frac{\pi m\sigma}{4r}f(r,\sigma)$$
$$S_{nm}=\frac{2n+2}{N^2}\sum_{r=1}^{N/2}R_{nm}(2r/N)\sum_{\sigma=1}^{8r}sin\frac{\pi m\sigma}{4r}f(r,\sigma)$$
1. 確定圖像的大小$N\times N$，即公式中的$N$；
2. 確定$r,\sigma$的范圍；
3. 利用Zernike多項式的遞推性質計算各階$R_{nm}(\rho)$，並結合上面Zernike矩計算公式，算出$C_{nm},S_{nm}$
4. 對$C_{nm},S_{nm}$求模，進而計算得到$|Z_{nm}|$
現在我們用Zernike矩來計算美女圖像在4種狀態下的值：
| 類別 | $log|Z_{11}|$ | $log|Z_{20}|$ | $log|Z_{22}|$ | $log|Z_{31}|$ | $log|Z_{40}|$ | $log|Z_{42}|$ | $log|Z_{44}|$ |
| 原圖 | 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++版的實現補上。
function [A_nm,zmlist,cidx,V_nm] = zernike(img,n,m) % 功能：計算輸入圖像的Zernike矩 % 輸入：img-灰度圖像 % n-階數 % 輸出：V_nm-n階的Zernike多項式，定義為在極坐標系中p，theta的函數 % cidx-表示虛部值 % A_nm-Zernike矩 if nargin>0 if nargin==1 n=0; end d=size(img); img=double(img); % 取步長 xstep=2/(d(1)-1); ystep=2/(d(2)-1); % 畫方格 [x,y]=meshgrid(-1:xstep:1,-1:ystep:1); circle1= x.^2 + y.^2; % 提取符合circle1<=1的數 inside=find(circle1<=1); % 構造size（d）*size(d)的矩陣 mask=zeros(d); % 構造size(inside)*size(inside)的全為1的矩陣賦值給mask（inside） mask(inside)=ones(size(inside)); % 計算圖像的復數表示 [cimg,cidx]=clipimg(img,mask); % 計算Z的實部和虛部 z=clipimg(x+i*y,mask); % 計算復數的模，sqrt(x,y),z=x+iy p=0.9*abs(z); ; % 計算復數z的輻角值（tanz） theta=angle(z); c=1; for order=1:length(n) n1=n(order); if nargin<3 m=zpossible(n1); end for r=1:length(m) V_nmt=zpoly(n1,m(r),p,theta); % conj是求復數的共軛 zprod=cimg.*conj(V_nmt); % (n1+1)/π*∑∑(zprod); 對於圖像而言求和代替了求積分 A_nm(c)=(n1+1)*sum(sum(zprod))/pi; zmlist(c,1:2)=[n1 m(r)]; if nargout==4 V_nm(:,c)=V_nmt; end c=c+1; end end else end %%%%%子函數%%%%% function [cimg,cindex,dim]=clipimg(img,mask) %功能：計算復數的實部和虛部 dim=size(img); cindex=find(mask~=0); cimg=img(cindex); return; %%%%%子函數%%%%% function [m]=zpossible(n) % 功能：判斷n是偶數還是奇數，是偶數時，m取0,2,4,6等,否則取奇數賦值m if iseven(n) m=0:2:n; else m=1:2:n; end return; %%%%%子函數%%%%% function [V_nm,mag,phase]=zpoly(n,m,p,theta) %功能：計算Zernike矩多項式 R_nm=zeros(size(p)); % 產生size(p)*size(p)的零矩陣賦給R_nm a=(n+abs(m))/2; b=(n-abs(m))/2; total=b; for s=0:total num=((-1)^s)*fac(n-s)*(p.^(n-2*s)); % (-1).-1*(n-s)!r.^(n-2*s) den=fac(s)*fac(a-s)*fac(b-s); % s!*(a-s)!*(b-s)! R_nm=R_nm + num/den; % R_nm是一個實數值的徑向多項式 end mag=R_nm; % 賦值 phase=m*theta; V_nm=mag.*exp(i*phase); % V_nm為n階的Zernike多項式，定義為在極坐標系中p，theta的函數 return; %%%%%子函數%%%%% function [factorial]=fac(n) %功能：求n的階乘 maxno=max(max(n)); zerosi=find(n<=0); %取n小於等於0的數 n(zerosi)=ones(size(zerosi)); factorial=n; findex=n; for i=maxno:-1:2 cand=find(findex>2); candidates=findex(cand); findex(cand)=candidates-1; factorial(cand)=factorial(cand).*findex(cand); end return; function [verdict]=iseven(candy) verdict=zeros(size(candy)); isint=find(isint(candy)==1); divided2=candy(isint)/2; evens=(divided2==floor(divided2)); verdict(isint)=evens; return; function [verdict]=isint(candy) verdict = double(round(candy))==candy; return;
5. 總結
不變矩的應用過程一般包括：
- 選擇合適的不變矩類型；
- 選擇分類器（如神經網絡、最短距離等）；
- 如果是神經網絡分類器，則需要計算學習樣例的不變矩去訓練神經網絡；
- 計算待識別對象的不變矩，輸入神經網絡就可得到待識別對象的類型，或者計算待識別對象不變矩與類別對象不變矩之間的距離，選擇最短距離的類別作為待識別對象的類別。
可以看出，不變矩作用主要目的是描述事物（圖像）的特征。人眼識別圖像的特征往往又表現為“求和”的形式，因此不變矩是對圖像元素進行了積分操作。
不變矩能夠描述圖像整體特征就是因為它具有平移不變形、比例不變性和旋轉不變性等性質。
然而，另一方面圖像的各階不變矩究竟代表的什么特征很難進行直觀的物理解釋。
6. 參考資料
[1] 《現代數字圖像處理》（matlab版）
[2] 正交多項式WIKI
[3] opencv形態描述