基于Retinex的人脸关照不变量的提取

为了解决人脸的关照问题，现在的研究方法主要有一下3类：

一、基于不变特征的方法、关照变化建模的方法和人脸图像关照归一化的方法。基于不变特征的方法是利用人脸的关照不变特征进行人脸识别，一般指利用朗博关照模型从图像中消除关照的影响。

二、光照建模的方法，主要是指在一个在一个子空间中表示不同关照引起的变化，并估计参数模型，但是此类方法计算量大，不能应用与实时的人脸识别系统中。

三、人脸图像归一化的方法，是指利用基本的图像处理技术对光照图像进行预处理。如：直方图均衡化和伽马校正等方法。

Retinex理论是Land提出的从信号分析角度出发，直接从单张图像中提取光照不变量。

步骤：

首先从原始图像中估计光照图像，然后在对数域里从原图中减去光照图像得到增强后的图像

Retinex方法主要包括两部分：关照的估计和归一化。

 

光照模型和对数变换：

根据Retinex理论：图像由入射分量和反射分量构成：

 

其中，表示图像的亮度分量，是缓慢变化的低频信号；表示物体的反射性质，是图像的光照不变量，是图像的高频信息；

[cpp]  view plain copy

1. //INFaceTool.h  
2. #include<opencv2\opencv.hpp>  
3. using namespace cv;  
4. Mat normalize8(Mat &X,int model=1);  
5. // The function applies the adaptive single scale retinex algorithm to an image.  
6. void adative_sigle_scale_retinex(Mat &img,Mat &L,Mat &Gradient,Mat&LocalInhomogenity,int T=10);  
7. //This is an auxialry function for computing the iterative convolution  
8. Mat convolute(Mat &X,Mat&y,Mat&N);  

[cpp]  view plain copy

1. //INFaceTool.cpp  
2. #include<opencv2\opencv.hpp>  
3. #include<core\core.hpp>  
4. #include"INFaceTool.h"  
5. using namespace cv;  
6. //The function adjusts the dynamic range of the grey scale image to the interval [0,255] or [0,1]  
7. Mat normalize8(Mat&X,int model)  
8. {     
9.     double max_v_x;  
10.     double min_v_x;  
11.     Mat Y(X.rows,X.cols,CV_64F);  
12.     Mat Ones=Mat::ones(X.rows,X.cols,CV_64F);  
13.     minMaxLoc(X,&min_v_x,&max_v_x,NULL,NULL);  
14.     if(model==1)  
15.     {  
16.         Y=((X-min_v_x*Ones)/(max_v_x*Ones-min_v_x*Ones))*255;  
17.         int nl=Y.rows;  
18.         int nc=Y.cols;  
19.         if(Y.isContinuous())  
20.     {  
21.         nc=nc*nl;  
22.         nl=1;  
23.     }  
24.     for(int j=0;j<nl;j++)  
25.     {  
26.         double*data=Y.ptr<double>(j);  
27.         for(int i=0;i<nc;i++)  
28.         {  
29.             //std::cout<<"data"<<data[i];  
30.             data[i]=ceil(data[i]);  
31.               
32.         }  
33.     }  
34.         return Y;  
35.     }  
36.     else  
37.     {  
38.         Y=(X-min_v_x*Ones)/(max_v_x*Ones-min_v_x*Ones);  
39.         return Y;  
40.     }  
41.       
42. }  
43. void adative_sigle_scale_retinex(Mat &img,Mat &L,Mat &Gradient,Mat&LocalInhomogenity,int T)  
44. {     
45.     Mat X=img.clone();  
46.     if(X.channels()==3)  
47.     {  
48.         cvtColor(X,X,CV_RGB2GRAY);  
49.     }  
50.   
51.     X.convertTo(X,CV_64F);  
52.     X=normalize8(X);  
53.     // Compute spatial gradient in x and y directions  
54.     int a=X.rows;  
55.     int b=X.cols;  
56.     Mat X1=Mat::zeros(a,b+2,CV_64F);  
57.     Mat X1r=X1(Rect(2,0,b,a)); X.copyTo(X1r);  
58.     Mat X2=Mat::zeros(a,b+2,CV_64F);  
59.     Mat X2r=X2(Rect(0,0,b,a));X.copyTo(X2r);  
60.     Mat Gx=X1(Rect(1,0,b,a))-X2(Rect(1,0,b,a));  
61.     X1=Mat::zeros(a+2,b,CV_64F);  
62.     X1r=X1(Rect(0,2,b,a));X.copyTo(X1r);  
63.     X2=Mat::zeros(a+2,b,CV_64F);  
64.     X2r=X2(Rect(0,0,b,a));X.copyTo(X2r);  
65.     Mat Gy=X1(Rect(0,1,b,a))-X2(Rect(0,1,b,a));  
66.     pow(Gx,2,Gx);  
67.     pow(Gy,2,Gy);  
68.     Mat I;  
69.     pow((Gx+Gy),0.5,I);  
70.       
71.     //Compute local inhomogenity  
72.     Mat tao=Mat::zeros(a,b,CV_64F);  
73.     Mat Xtmp=Mat::zeros(a+2,b+2,CV_64F);  
74.     Mat Xtmpr=Xtmp(Rect(1,1,b,a));X.copyTo(Xtmpr);  
75.     float suma=0;  
76.     for(int i=1;i<a+1;i++)  
77.     {  
78.         for(int j=1;j<b+1;j++)  
79.         {  
80.             suma=0;  
81.             for(int k=-1;k<2;k++)  
82.             { for(int h=-1;h<2;h++)  
83.                 {  
84.                    suma=suma+abs(X.at<double>(i-1,j-1)-Xtmp.at<double>(i+k,j+h));  
85.                 }  
86.             }  
87.             suma=suma/9;  
88.             tao.at<double>(i-1,j-1)=suma;  
89.   
90.         }  
91.     }  
92.     double tao_min,tao_max;  
93.     minMaxLoc(tao,&tao_min,&tao_max,NULL,NULL);  
94.     Mat tao_slash=(tao-tao_min)/(tao_max-tao_min);  
95.     int nl=tao_slash.rows;  
96.     int nc=tao_slash.cols;  
97.     if(tao_slash.isContinuous())  
98.     {  
99.         nc=nc*nl;  
100.         nl=1;  
101.     }  
102.     for(int j=0;j<nl;j++)  
103.     {  
104.         double*data=tao_slash.ptr<double>(j);  
105.         for(int i=0;i<nc;i++)  
106.         {  
107.             //std::cout<<"data"<<data[i];  
108.             data[i]=sin(CV_PI/2*data[i]);  
109.               
110.         }  
111.     }  
112.       
113.     Scalar m1=mean(I);  
114.     Scalar m2=mean(tao_slash);  
115.     double meanI=m1.val[0];  
116.     double meanH=m2.val[0];  
117.     //Set needed parameters if they are not provided as inputs  
118.     double S=10*exp(-(meanI)/10);  
119.     double h=0.1*exp(-(meanH)/0.1);  
120.   
121.     //Determine weight functions  
122.     Mat sqrt_tao_slash;  
123.     Mat sqrt_I;  
124.     pow(tao_slash/h,0.5,sqrt_tao_slash);  
125.     pow(I/S,0.5,sqrt_I);  
126.     Mat alpha=1/(1+sqrt_tao_slash);  
127.     Mat beta=1/(1+sqrt_I);  
128.     Mat weight=alpha.mul(beta);  
129.     //precompute Ns  
130.     Mat w=Mat::zeros(a+2,b+2,CV_64F);  
131.     Mat wr=w(Rect(1,1,b,a));  
132.     weight.copyTo(wr);  
133.     Mat N=Mat::zeros(a,b,CV_64F);  
134.     for(int i=1;i<a+1;i++)  
135.     {  
136.           
137.         for(int j=1;j<b+1;j++)  
138.         {   suma=0;  
139.             for(int k=-1;k<2;k++)  
140.             {  
141.                 for(h=-1;h<2;h++)  
142.                 {  
143.                     suma=suma+(w.at<double>(i+k,j+h));  
144.                 }  
145.             }  
146.              N.at<double>(i-1,j-1)=suma;  
147.         }  
148.     }  
149.     // Start iterative convolution  
150.     Mat L_old;  
151.     X.copyTo(L_old);  
152.     Mat L_new,L_new_s;  
153.     L_new=Mat::zeros(a,b,CV_64F);  
154.     L_new_s=Mat::zeros(a,b,CV_64F);  
155.     for(int i=0;i<T;i++)  
156.     {  
157.         L_new_s=convolute(L_old,weight,N);  
158.         L_new=max(L_new_s,L_old);  
159.         L_new.reshape(0,a);  
160.         L_new.copyTo(L_old);  
161.     }     
162.       
163.     //Produce ilumination invariant representation of input image  
164.      nc=X.cols;  
165.      nl=X.rows;  
166.     if(X.isContinuous())  
167.     {  
168.         nc=nc*nl;  
169.         nl=1;  
170.     }  
171. // R=Mat::zeros(a,b,CV_64F);  
172.  L=Mat::zeros(a,b,CV_64F);;  
173.     for(int j=0;j<nl;j++)  
174.     { //double*dataR=R.ptr<double>(j);  
175.       double*dataX=X.ptr<double>(j);  
176.       double*dataL=L.ptr<double>(j);  
177.       double*dataL_new=L_new.ptr<double>(j);  
178.         for(int i=0;i<nc;i++)  
179.         {  
180.             //dataR[i]=log(dataX[i]+1)-log(dataL_new[i]+1);  
181.             dataL[i]=log(dataL_new[i]+1);  
182.               
183.         }  
184.     }     
185.   
186.     L=normalize8(L);  
187.    // R=normalize8(R);  
188.     L.convertTo(L,CV_8U);  
189.     //R.convertTo(R,CV_8U);  
190.     I.convertTo(I,CV_8U);  
191.     I.copyTo(Gradient);  
192.     tao_slash.copyTo(tao_slash);  
193.     //tao_slash.copyTo(LocalInhomogenity);  
194.   
195.       
196.   
197. }  
198. //This is an auxialry function for computing the iterative convolution  
199. Mat convolute(Mat &X,Mat&y,Mat&N)  
200. {  
201.     int a=X.rows;  
202.     int b=X.cols;  
203.     Mat X1=Mat::zeros(a+2,b+2,CV_64F);  
204.     Mat X1r=X1(Rect(1,1,b,a));X.copyTo(X1r);  
205.   
206.     Mat w=Mat::zeros(a+2,b+2,CV_64F);  
207.     Mat wr=w(Rect(1,1,b,a));y.copyTo(wr);  
208.     Mat Y=Mat::zeros(a,b,CV_64F);  
209.     double tt;  
210.     for(int i=1;i<a;i++)  
211.     {  
212.         for(int j=1;j<b;j++)  
213.         {  
214.               
215.             tt=sum(X1(Rect(j-1,i-1,3,3)).mul(w(Rect(j-1,i-1,3,3)))).val[0];  
216.             Y.at<double>(i-1,j-1)=tt/N.at<double>(i-1,j-1);  
217.       
218.         }  
219.     }  
220.     return Y;  
221.       
222. }