MATLAB代码 求分析 何的去雾算法里面的暗通道算法 每一句都是啥意思啊？又分别对应论文里的啥原理？
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了MATLAB代码 求分析 何的去雾算法里面的暗通道算法 每一句都是啥意思啊？又分别对应论文里的啥原理？相关的知识，希望对你有一定的参考价值。
function dark = darkChannel(imRGB)
r=imRGB(:,:,1);
g=imRGB(:,:,2);
b=imRGB(:,:,3);
[m n] = size(r);
a = zeros(m,n);
for i = 1: m
for j = 1: n
a(i,j) = min(r(i,j), g(i,j));
a(i,j)= min(a(i,j), b(i,j));
end
end
d = ones(15,15);
fun = @(block_struct)min(min(block_struct.data))*d;
dark = blockproc(a, [15 15], fun);
dark = dark(1:m, 1:n);
r=imRGB(:,:,1);
g=imRGB(:,:,2);
b=imRGB(:,:,3); 分别提取三色的灰度图
[m n] = size(r); 提取单色图矩阵的宽度和长度
a = zeros(m,n); 创建m*n的零矩阵a
for i = 1: m
for j = 1: n
a(i,j) = min(r(i,j), g(i,j));
a(i,j)= min(a(i,j), b(i,j)); 依次比较三色分量的最小值提取为暗通道图
end
end
d = ones(15,15); 创建15*15的单位矩阵
fun = @(block_struct)min(min(block_struct.data))*d;
dark = blockproc(a, [15 15], fun); 将图片分成15*15的小块并将每一块变成其中的最小值
dark = dark(1:m, 1:n);
我也是新手啊兄弟只能帮你到这儿了 参考技术A %%%%暗原色去雾算法是建立在户外自然场景暗通道优先法则的基础上的去雾方法，其实就是解一个方程
%%I(x)=J(x)t(x)+A(1-t(x));其中I(x)是受到雾气污染的图像。J(x)是我们需要求的去雾后的图像
%%t(x)是天空中云层的透射分布率，A是天空的亮度
%暗影去雾算法
% 原始图像
img_name = imread('C:\Users\Administrator\Desktop\桌面文件\专业课\设计\2.bmp');
I = double(img_name) / 255;
% 获取图像大小
[h,w,c] = size(I);
%去雾系数
w0 = 0.95;
img_size = w * h;
%初始化结果图像
dehaze = zeros(h,w,c);
%初始化暗影通道图像
win_dark = zeros(h,w);
for i=1:h
for j=1:w
win_dark(i,j)=min(I(i,j,:));%将三个通道中最暗的值赋给dark_I(i,j),显然，三维图变成了二维图
end
end
win_dark = ordfilt2(win_dark,1,ones(9,9),'symmetric');
%计算大气亮度A，相关原理详见论文“Single Image Haze Removal Using Dark Channel Prior”
dark_channel = win_dark;
A = max(max(dark_channel));
[i,j] = find(dark_channel==A);
i = i(1);
j = j(1);
A = mean(I(i,j,:));
%计算初始的transmission map
transmission = 1 - w0 * win_dark / A;
%用guided filter对trasmission map做soft matting
gray_I = I(:,:,1);%这里gray_I 可以是RGB图像中任何一个通道
p = transmission;
r = 80;
eps = 10^-3;
transmission_filter = guidedfilter(gray_I, p, r, eps);
t0=0.1;
t1 = max(t0,transmission_filter);
for i=1:c
for j=1:h
for l=1:w
dehaze(j,l,i)=(I(j,l,i)-A)/t1(j,l)+A;
end
end
end
figure,
imshow(I);title('去雾前')
figure,
imshow(dehaze);title('去雾后')
%里面用到的两个函数如下：（这两个函数是何凯明guided filter中的，详细原理可以参考博客园中laviewpb的分析）
function q = guidedfilter(I, p, r, eps)
% GUIDEDFILTER O(1) time implementation of guided filter.
%
% - guidance image: I (should be a gray-scale/single channel image)
% - filtering input image: p (should be a gray-scale/single channel image)
% - local window radius: r
% - regularization parameter: eps
[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;
a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;
q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
end
function imDst = boxfilter(imSrc, r)
% BOXFILTER O(1) time box filtering using cumulative sum
%
% - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
% - Running time independent of r;
% - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
% - But much faster.
[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));
%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);
%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end 参考技术B 这个哈解决的，要求告诉我
基于暗通道优先算法的去雾应用(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 }
截图：
以上是关于MATLAB代码 求分析 何的去雾算法里面的暗通道算法 每一句都是啥意思啊？又分别对应论文里的啥原理？的主要内容，如果未能解决你的问题，请参考以下文章