基于暗通道优先算法的去雾应用(Matlab/C++)

基于暗通道优先的单幅图像去雾算法（Matlab/C++）

算法原理：

            参见论文：Single Image Haze Removal Using Dark Channel Prior  [1]

       ① 暗通道定义

     何恺明 通过对大量在户外拍摄的自然景物图片进行统计分析得出一个结论：在绝大多数非天空的局部区域里，某一些像素总会（至少一个颜色通道）具有很低的值。换言之，该区域光强度的最小值是个很小的数（趋于0）。

基于上述结论，我们定义暗通道，用公式描述，对于一幅图像J有如下式子：

 也就是说以像素点x为中心，分别取三个通道内窗口Ω内的最小值，然后再取三个通道的最小值作为像素点x的暗通道的值，如下图所示：

   

 Jc代表J的某一个颜色通道，而Ω(x)是以x为中心的一块方形区域。我们观察得出，除了天空方位，Jdark的强度总是很低并且趋近于0。如果J是户外的无雾图像，我们把Jdark称为J的暗原色,并且把以上观察得出的经验性规律称为暗原色先验。

②大气物理模型

要想从物理模型角度对有雾图像进行清晰化处理，就要了解有雾图像的物理成因，那么就要了解雾天的大气散射模型。

大气散射物理模型包含两部分，第一部分称为直接衰减项(Direct Attenuation)也称直接传播，第二部分称为大气光照(Airlight)

用公式表示如下:

   I是观测到的有雾图像，J是景物反射光强度（也就是清晰的无雾图像），A是全局大气光照强度，t用来描述光线通过介质透射到成像设备过程中没有被散射的部分，去雾的目标就是从I中复原J。那么也就是要通过I求A和t。

    方程右边的第一项J(x)t(x) 叫做直接衰减项，第二项A(1-t(x))则是大气光照。直接衰减项描述的是景物光线在透射媒介中经衰减后的部分，而大气光则是由前方散射引起的，会导致景物颜色的偏移。因为大气层可看成各向同性的，透射率t可表示为：

β为大气的散射系数，该式表明景物光线是随着景物深度d按指数衰减的。

 ③求解透射率t

在论文[1]中，作者给出了推导过程，这里就不再重复，其最后得到透射率t的公式如下：

Ic为输入的有雾图像，对其除以全局大气光照Ac后在利用暗通道定义公式进行求解暗通道。w（0<w<1）是雾的保留系数通常取0.95。

这里需要值得注意的是，求得的t是粗透射率图，并不能直接带入大气模型公式求解，所以需要进行细化后再处理。细化过程见⑤，Ac为全局大气光照，其求法见④。

④求解全局大气光照Ac

论文[1]中作者给出求解全局大气光照的过程如下：

1.首先对输入的有雾图像I求解其暗通道图像Jdark。

2.选择暗通道Jdark内图像总像素点个数（N_imagesize）千分之一（N=N_imagesize/1000）个最亮的像素点，并记录这些像素点（x,y）坐标。

3.再根据这些点的坐标分别在原图像I的三个通道（r,g,b）内找到这些像素点并加和得到（sum_r,sum_g,sum_b）.

4.Ac=[Ar,Ag,Ab]. 其中Ar=sum_r/N;   Ag=sum_g/N;   Ab=sum_b/N.

⑤细化透射率t

    作者在论文[1]中使用了软抠图（soft matting）的方法，详见论文如下：

A Closed-Form Solution to Natural Image Matting[2], 作者：Anat Levin

使用软抠图法对得到的粗透射率t~进行细化。由于这个方法时间和内存花费比较大，后来作者又使用指导性滤波器进行细化粗透射图。效果上指导性滤波要稍差于软抠图法，但在时间和内存花费上具有明显优势，因此这里我们使用指导性滤波器进行细化粗透射率t~。（注：Matlab代码中也附带软抠图法细化透射率的代码。内容见文件包）

关与指导性滤波的详细内容见论文：Guided Image Filtering [3] 作者：何恺明 

⑥求解最后清晰图像

现在，我们得到了A和t，那么带入大气模型公式：

这里，t0参数用来限定透射率t的下限值，其作用也就是在输入图像的浓雾区域保留一定的雾。

附录：

     

ReadMe.txt：

感谢论文Single Image Haze Removal Using Dark Channel Prior作者何凯明。
基于暗通道优先去雾算法的Matlab /C++ 源码为免费开源。只供研究学习使用，如果引用请注明原开发者，及出处。

Matlab版中指导性滤波的源码为原作者何凯明提供。

开发者：赵常凯（不包括Matlab版中 指导性滤波的源码）
时间：2013.8.19

       Matlab 源代码  下载 5.69MB

       C++ 源代码     下载 1.68MB

       说明:C++ 是基于 OpenCV2.2 + Qt 4.8.5 编写的，如果不匹配需要搭建 Qt+OpenCV开发环境。 如果只想查看核心源码文件，请查看Dehazor.cpp文件。

  1 /*--------------------------------------------------------------------------------*\
  2 
  3    This program is free software; permission is hereby granted to use, copy, modify,
  4    and distribute this source code, or portions thereof, for any purpose, without fee,
  5    subject to the restriction that the copyright notice may not be removed
  6    or altered from any source or altered source distribution.
  7    The software is released on an as-is basis and without any warranties of any kind.
  8    In particular, the software is not guaranteed to be fault-tolerant or free from
  9    failure.
 10 
 11    The author disclaims all warranties with regard to this software, any use,
 12    and any consequent failure, is purely the responsibility of the user.
 13 
 14    Copyright (C) 2013-2016 Changkai Zhao, www.cnblogs.com/changkaizhao
 15 \*-------------------------------------------------------------------------------*/
 16 #include "dehazor.h"
 17 
 18 cv::Mat Dehazor::process(const cv::Mat &image)
 19 {
 20     int dimr=image.rows;
 21     int dimc=image.cols;
 22 
 23     cv::Mat rawtemp;
 24     cv::Mat refinedImage_temp;
 25     cv::Mat output_b_temp(image.rows,image.cols,CV_32F);
 26     cv::Mat output_g_temp(image.rows,image.cols,CV_32F);
 27     cv::Mat output_r_temp(image.rows,image.cols,CV_32F);
 28     cv::Mat output_temp(dimr,dimc,image.type());
 29     rawImage.create(image.rows,image.cols,CV_8U); 
 30     refinedImage_temp.create(image.rows,image.cols,CV_32F);
 31     rawtemp.create(image.rows,image.cols,CV_8U);
 32 
 33 
 34     float sumb=0;float sumg=0;float sumr=0;
 35     float Air_b;  float Air_g; float Air_r;
 36     cv::Point2i pt;
 37 
 38     int dimone=floor(dimr*dimc*0.001);
 39     cv::Mat quantile;
 40     quantile=cv::Mat::zeros(2,dimone,CV_8U);
 41 
 42     int dx=floor(windowsize/2);
 43 //    cv::Mat imagetemp;
 44 //    imagetemp.create(image.rows,image.cols,CV_32FC3);
 45 
 46     for(int j=0;j<dimr;j++){
 47         for(int i=0;i<dimc;i++){
 48             int min;
 49             image.at<cv::Vec3b>(j,i)[0] >= image.at<cv::Vec3b>(j,i)[1]?
 50                min=image.at<cv::Vec3b>(j,i)[1]:
 51                min=image.at<cv::Vec3b>(j,i)[0];
 52             if(min>=image.at<cv::Vec3b>(j,i)[2])
 53                 min=image.at<cv::Vec3b>(j,i)[2];
 54             rawImage.at<uchar>(j,i)=min;
 55         }
 56     }
 57 
 58     for(int j=0;j<dimr;j++){
 59         for(int i=0;i<dimc;i++){
 60             int min=255;
 61 
 62             int jlow=j-dx;int jhigh=j+dx;
 63             int ilow=i-dx;int ihigh=i+dx;
 64 
 65             if(ilow<=0)
 66                 ilow=0;
 67             if(ihigh>=dimc)
 68                 ihigh=dimc-1;
 69             if(jlow<=0)
 70                 jlow=0;
 71             if(jhigh>=dimr)
 72                 jhigh=dimr-1;
 73 
 74             for(int m=jlow;m<=jhigh;m++){
 75                 for(int n=ilow;n<=ihigh;n++){
 76                     if(min>=rawImage.at<uchar>(m,n))
 77                         min=rawImage.at<uchar>(m,n);
 78                 }
 79             }
 80            rawtemp.at<uchar>(j,i)=min;
 81         }
 82     }
 83 
 84     for(int i=0;i<dimone;i++){
 85         cv::minMaxLoc(rawtemp,0,0,0,&pt);
 86         sumb+= image.at<cv::Vec3b>(pt.y,pt.x)[0];
 87         sumg+= image.at<cv::Vec3b>(pt.y,pt.x)[1];
 88         sumr+= image.at<cv::Vec3b>(pt.y,pt.x)[2];
 89         rawtemp.at<uchar>(pt.y,pt.x)=0;
 90     }
 91     Air_b=sumb/dimone;
 92     Air_g=sumg/dimone;
 93     Air_r=sumr/dimone;
 94 
 95 
 96     cv::Mat layb;  cv::Mat Im_b;
 97     cv::Mat layg;  cv::Mat Im_g;
 98     cv::Mat layr;  cv::Mat Im_r;
 99 
100     // create vector of 3 images
101     std::vector<cv::Mat> planes;
102     // split 1 3-channel image into 3 1-channel images
103     cv::split(image,planes);
104 
105     layb=planes[0];
106     layg=planes[1];
107     layr=planes[2];
108     Im_b=planes[0];
109     Im_g=planes[1];
110     Im_r=planes[2];
111 
112     layb.convertTo(layb, CV_32F);
113     layg.convertTo(layg, CV_32F);
114     layr.convertTo(layr, CV_32F);
115 
116     Im_b.convertTo(Im_b, CV_32F);
117     Im_g.convertTo(Im_g, CV_32F);
118     Im_r.convertTo(Im_r, CV_32F);
119 
120 
121     for (int j=0; j<dimr; j++) {
122         for (int i=0; i<dimc; i++) {
123         // process each pixel ---------------------
124         layb.at<float>(j,i)=layb.at<float>(j,i)/Air_b;
125 
126         layg.at<float>(j,i)=layg.at<float>(j,i)/Air_g;
127 
128         layr.at<float>(j,i)=layr.at<float>(j,i)/Air_r;
129         // end of pixel processing ----------------
130         } // end of line
131     }
132 
133 
134     rawtemp.convertTo(rawtemp,CV_32F);
135 
136 
137     for(int j=0;j<dimr;j++){
138         for(int i=0;i<dimc;i++){
139             float min;
140             layb.at<float>(j,i) >= layg.at<float>(j,i)?
141                min=layg.at<float>(j,i):
142                min=layb.at<float>(j,i);
143             if(min>=layr.at<float>(j,i))
144                 min=layr.at<float>(j,i);
145             rawtemp.at<float>(j,i)=min;
146         }
147     }
148     for(int j=0;j<dimr;j++){
149         for(int i=0;i<dimc;i++){
150             float min=1;
151 
152             int jlow=j-dx;int jhigh=j+dx;
153             int ilow=i-dx;int ihigh=i+dx;
154 
155             if(ilow<=0)
156                 ilow=0;
157             if(ihigh>=dimc)
158                 ihigh=dimc-1;
159             if(jlow<=0)
160                 jlow=0;
161             if(jhigh>=dimr)
162                 jhigh=dimr-1;
163 
164             for(int m=jlow;m<=jhigh;m++){
165                 for(int n=ilow;n<=ihigh;n++){
166                     if(min>=rawtemp.at<float>(m,n))
167                         min=rawtemp.at<float>(m,n);
168                 }
169             }
170            rawImage.at<uchar>(j,i)=(1-(float)fog_reservation_factor*min)*255;
171         }
172     }
173 
174    refinedImage_temp=guildedfilter_color(image,rawImage,localwindowsize,eps);
175 
176    for(int j=0;j<dimr;j++){
177        for(int i=0;i<dimc;i++){
178            if(refinedImage_temp.at<float>(j,i)<0.1)
179                refinedImage_temp.at<float>(j,i)=0.1;
180        }
181    }
182 
183    cv::Mat onemat(dimr,dimc,CV_32F,cv::Scalar(1));
184 
185    cv::Mat air_bmat(dimr,dimc,CV_32F);
186    cv::Mat air_gmat(dimr,dimc,CV_32F);
187    cv::Mat air_rmat(dimr,dimc,CV_32F);
188 
189    cv::addWeighted(onemat,Air_b,onemat,0,0,air_bmat);
190    cv::addWeighted(onemat,Air_g,onemat,0,0,air_gmat);
191    cv::addWeighted(onemat,Air_r,onemat,0,0,air_rmat);
192 
193 
194    output_b_temp=Im_b-air_bmat;
195    output_g_temp=Im_g-air_gmat;
196    output_r_temp=Im_r-air_rmat;
197 
198    output_b_temp=output_b_temp.mul(1/refinedImage_temp,1)+air_bmat;
199    output_g_temp=output_g_temp.mul(1/refinedImage_temp,1)+air_gmat;
200    output_r_temp=output_r_temp.mul(1/refinedImage_temp,1)+air_rmat;
201 
202 
203    output_b_temp.convertTo(output_b_temp,CV_8U);
204    output_g_temp.convertTo(output_g_temp,CV_8U);
205    output_r_temp.convertTo(output_r_temp,CV_8U);
206 
207    for(int j=0;j<dimr;j++){
208        for(int i=0;i<dimc;i++){
209            output_temp.at<cv::Vec3b>(j,i)[0]=output_b_temp.at<uchar>(j,i);
210            output_temp.at<cv::Vec3b>(j,i)[1]=output_g_temp.at<uchar>(j,i);
211            output_temp.at<cv::Vec3b>(j,i)[2]=output_r_temp.at<uchar>(j,i);
212        }
213    }
214 
215    cv::Mat nom_255(dimr,dimc,CV_32F,cv::Scalar(255));
216 
217    cv::Mat ref_temp(dimr,dimc,CV_32F);
218    ref_temp=refinedImage_temp;
219    ref_temp=ref_temp.mul(nom_255,1);
220    ref_temp.convertTo(refinedImage,CV_8U);
221 
222 
223     return output_temp;
224 }
225 cv::Mat Dehazor::boxfilter(cv::Mat &im, int r)
226 {
227     //im is a CV_32F type mat [0,1] (normalized)
228     //output is the same size to im;
229 
230     int hei=im.rows;
231     int wid=im.cols;
232     cv::Mat imDst;
233     cv::Mat imCum;
234 
235 
236     imDst=cv::Mat::zeros(hei,wid,CV_32F);
237     imCum.create(hei,wid,CV_32F);
238 
239     //cumulative sum over Y axis
240     for(int i=0;i<wid;i++){
241         for(int j=0;j<hei;j++){
242             if(j==0)
243                 imCum.at<float>(j,i)=im.at<float>(j,i);
244             else
245                 imCum.at<float>(j,i)=im.at<float>(j,i)+imCum.at<float>(j-1,i);
246         }
247     }
248 
249 
250     //difference over Y axis
251     for(int j=0;j<=r;j++){
252         for(int i=0;i<wid;i++){
253             imDst.at<float>(j,i)=imCum.at<float>(j+r,i);
254         }
255     }
256     for(int j=r+1;j<=hei-r-1;j++){
257         for(int i=0;i<wid;i++){
258             imDst.at<float>(j,i)=imCum.at<float>(j+r,i)-imCum.at<float>(j-r-1,i);
259         }
260     }
261     for(int j=hei-r;j<hei;j++){
262         for(int i=0;i<wid;i++){
263             imDst.at<float>(j,i)=imCum.at<float>(hei-1,i)-imCum.at<float>(j-r-1,i);
264         }
265     }
266 
267 
268     //cumulative sum over X axis
269       for(int j=0;j<hei;j++){
270           for(int i=0;i<wid;i++){
271               if(i==0)
272                   imCum.at<float>(j,i)=imDst.at<float>(j,i);
273               else
274                   imCum.at<float>(j,i)=imDst.at<float>(j,i)+imCum.at<float>(j,i-1);
275           }
276       }
277       //difference over X axis
278       for(int j=0;j<hei;j++){
279           for(int i=0;i<=r;i++){
280               imDst.at<float>(j,i)=imCum.at<float>(j,i+r);
281           }
282       }
283       for(int j=0;j<hei;j++){
284           for(int i=r+1;i<=wid-r-1;i++){
285               imDst.at<float>(j,i)=imCum.at<float>(j,i+r)-imCum.at<float>(j,i-r-1);
286           }
287       }
288       for(int j=0;j<hei;j++){
289           for(int i=wid-r;i<wid;i++){
290               imDst.at<float>(j,i)=imCum.at<float>(j,wid-1)-imCum.at<float>(j,i-r-1);
291           }
292       }
293 
294     return imDst;
295 }
296 cv::Mat Dehazor::guildedfilter_color(const cv::Mat &Img, cv::Mat &p, int r, float &epsi)
297 {
298 
299     int hei=p.rows;
300     int wid=p.cols;
301 
302     cv::Mat matOne(hei,wid,CV_32F,cv::Scalar(1));
303     cv::Mat N;
304 
305     N=boxfilter(matOne,r);
306 
307 
308 
309     cv::Mat mean_I_b(hei,wid,CV_32F);
310     cv::Mat mean_I_g(hei,wid,CV_32F);
311     cv::Mat mean_I_r(hei,wid,CV_32F);
312     cv::Mat mean_p(hei,wid,CV_32F);
313 
314 
315     cv::Mat Ip_b(hei,wid,CV_32F);
316     cv::Mat Ip_g(hei,wid,CV_32F);
317     cv::Mat Ip_r(hei,wid,CV_32F);
318     cv::Mat mean_Ip_b(hei,wid,CV_32F);
319     cv::Mat mean_Ip_g(hei,wid,CV_32F);
320     cv::Mat mean_Ip_r(hei,wid,CV_32F);
321     cv::Mat cov_Ip_b(hei,wid,CV_32F);
322     cv::Mat cov_Ip_g(hei,wid,CV_32F);
323     cv::Mat cov_Ip_r(hei,wid,CV_32F);
324 
325     cv::Mat II_bb(hei,wid,CV_32F);
326     cv::Mat II_gg(hei,wid,CV_32F);
327     cv::Mat II_rr(hei,wid,CV_32F);
328     cv::Mat II_bg(hei,wid,CV_32F);
329     cv::Mat II_br(hei,wid,CV_32F);
330     cv::Mat II_gr(hei,wid,CV_32F);
331 
332     cv::Mat var_I_bb(hei,wid,CV_32F);
333     cv::Mat var_I_gg(hei,wid,CV_32F);
334     cv::Mat var_I_rr(hei,wid,CV_32F);
335     cv::Mat var_I_bg(hei,wid,CV_32F);
336     cv::Mat var_I_br(hei,wid,CV_32F);
337     cv::Mat var_I_gr(hei,wid,CV_32F);
338 
339     cv::Mat layb;
340     cv::Mat layg;
341     cv::Mat layr;
342     cv::Mat P_32;
343 
344     // create vector of 3 images
345     std::vector<cv::Mat> planes;
346     // split 1 3-channel image into 3 1-channel images
347     cv::split(Img,planes);
348 
349     layb=planes[0];
350     layg=planes[1];
351     layr=planes[2];
352 
353     layb.convertTo(layb, CV_32F);
354     layg.convertTo(layg, CV_32F);
355     layr.convertTo(layr, CV_32F);
356 
357     p.convertTo(P_32,CV_32F);
358     cv::Mat nom_255(hei,wid,CV_32F,cv::Scalar(255));
359 
360 
361 
362     layb=layb.mul(1/nom_255,1);
363     layg=layg.mul(1/nom_255,1);
364     layr=layr.mul(1/nom_255,1);
365     P_32=P_32.mul(1/nom_255,1);
366 
367     cv::Mat mean_I_b_temp=boxfilter(layb,r);
368     cv::Mat mean_I_g_temp=boxfilter(layg,r);
369     cv::Mat mean_I_r_temp=boxfilter(layr,r);
370     cv::Mat mean_p_temp=boxfilter(P_32,r);
371 
372     mean_I_b=mean_I_b_temp.mul(1/N,1);
373     mean_I_g=mean_I_g_temp.mul(1/N,1);
374     mean_I_r=mean_I_r_temp.mul(1/N,1);
375     mean_p=mean_p_temp.mul(1/N,1);
376 
377     Ip_b=layb.mul(P_32,1);
378     Ip_g=layg.mul(P_32,1);
379     Ip_r=layr.mul(P_32,1);
380 
381     cv::Mat mean_Ip_b_temp=boxfilter(Ip_b,r);
382     cv::Mat mean_Ip_g_temp=boxfilter(Ip_g,r);
383     cv::Mat mean_Ip_r_temp=boxfilter(Ip_r,r);
384 
385     mean_Ip_b=mean_Ip_b_temp.mul(1/N,1);
386     mean_Ip_g=mean_Ip_g_temp.mul(1/N,1);
387     mean_Ip_r=mean_Ip_r_temp.mul(1/N,1);
388 
389     cov_Ip_b=mean_Ip_b-mean_I_b.mul(mean_p,1);
390     cov_Ip_g=mean_Ip_g-mean_I_g.mul(mean_p,1);
391     cov_Ip_r=mean_Ip_r-mean_I_r.mul(mean_p,1);
392 
393 
394 //     variance of I in each local patch: the matrix Sigma in Eqn (14).
395 //     Note the variance in each local patch is a 3x3 symmetric matrix:
396 //                   bb, bg, br
397 //           Sigma = bg, gg, gr
398 //                   br, gr, rr
399     II_bb=layb.mul(layb,1);
400     II_gg=layg.mul(layg,1);
401     II_rr=layr.mul(layr,1);
402     II_bg=layb.mul(layg,1);
403     II_br=layb.mul(layr,1);
404     II_gr=layg.mul(layr,1);
405 
406     cv::Mat bb_box=boxfilter(II_bb,r);
407     cv::Mat gg_box=boxfilter(II_gg,r);
408     cv::Mat rr_box=boxfilter(II_rr,r);
409     cv::Mat bg_box=boxfilter(II_bg,r);
410     cv::Mat br_box=boxfilter(II_br,r);
411     cv::Mat gr_box=boxfilter(II_gr,r);
412 
413     var_I_bb=bb_box.mul(1/N,1)-mean_I_b.mul(mean_I_b);
414     var_I_gg=gg_box.mul(1/N,1)-mean_I_g.mul(mean_I_g);
415     var_I_rr=rr_box.mul(1/N,1)-mean_I_r.mul(mean_I_r);
416     var_I_bg=bg_box.mul(1/N,1)-mean_I_b.mul(mean_I_g);
417     var_I_br=br_box.mul(1/N,1)-mean_I_b.mul(mean_I_r);
418     var_I_gr=gr_box.mul(1/N,1)-mean_I_g.mul(mean_I_r);
419 
420     cv::Mat a_b(hei,wid,CV_32F);
421     cv::Mat a_g(hei,wid,CV_32F);
422     cv::Mat a_r(hei,wid,CV_32F);
423 
424     cv::Mat b(hei,wid,CV_32F);
425     cv::Mat sigma(3,3,CV_32F,cv::Scalar(0));
426     cv::Mat inv_sigma(3,3,CV_32F);
427 
428     for(int j=0;j<hei;j++){
429         for(int i=0;i<wid;i++){
430             sigma.at<float>(0,0)=var_I_rr.at<float>(j,i)+epsi;
431             sigma.at<float>(0,1)=var_I_gr.at<float>(j,i);
432             sigma.at<float>(0,2)=var_I_br.at<float>(j,i);
433             sigma.at<float>(1,0)=var_I_gr.at<float>(j,i);
434             sigma.at<float>(2,0)=var_I_br.at<float>(j,i);
435             sigma.at<float>(1,1)=var_I_gg.at<float>(j,i)+epsi;
436             sigma.at<float>(2,2)=var_I_bb.at<float>(j,i)+epsi;
437             sigma.at<float>(1,2)=var_I_bg.at<float>(j,i);
438             sigma.at<float>(2,1)=var_I_bg.at<float>(j,i);
439             inv_sigma=sigma.inv(cv::DECOMP_LU);
440 
441             a_r.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,0)+
442                                cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,0)+
443                                cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,0);
444             a_g.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,1)+
445                                cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,1)+
446                                cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,1);
447             a_b.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,2)+
448                                cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,2)+
449                                cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,2);
450 
451 
452         }
453     }
454     b=mean_p-a_b.mul(mean_I_b,1)-a_g.mul(mean_I_g,1)-a_r.mul(mean_I_r,1);
455 
456     cv::Mat box_ab=boxfilter(a_b,r);
457     cv::Mat box_ag=boxfilter(a_g,r);
458     cv::Mat box_ar=boxfilter(a_r,r);
459     cv::Mat box_b=boxfilter(b,r);
460     cv::Mat q(hei,wid,CV_32F);
461 
462 
463     q=box_ab.mul(layb,1)+box_ag.mul(layg,1)+box_ar.mul(layr,1)+box_b;
464     q=q.mul(1/N,1);
465 
466 
467     return q;
468 
469 }

       截图：