精华内容
下载资源
问答
  • 基于暗通道先验的去雾算法:2009CVPR+导向滤波+容差控制 基于CLAHE限制自适应直方图去雾算法
  • 本资源是基于导向滤波的暗通道先验图像去雾代码,能够完整的运行。
  • 暗通道先验实现图像去雾,darkchannel,基于透射率大气光,暗通道、
  • 针对暗通道先验去雾算法受限于天空区域,提出了一种天空优化的暗通道去雾算法。在原有公式的基础上,根据天空大气光中值及天空区域的占比添加天空区域约束因子,使其适用于天空区域的处理。同时,在原含雾图像对透射...
  • 暗通道先验图像去雾

    2014-11-10 16:24:39
    本程序实现的是CVPR2009 best paper 利用暗通道先验的单幅图像去雾
  • 暗通道先验图像去雾程序
  • 暗通道先验去雾matlab代码,包括导向滤波代码
  • 传统基于暗通道先验的图像去雾算法不能有效去除有雾图像在景深突变处的雾点,边界处容易引起光晕效应,对此提出一种基于暗通道先验的自适应超像素去雾算法.首先,在暗通道的获取过程中引入自适应方法判断当前像素邻域内...
  • 暗通道先验算法

    万次阅读 2018-10-12 17:32:04
    本文创新点: 作者提出暗通道先验算法,即:在户外无雾图像中(非天空区域)的大部分局部区域,存在一些像素点(暗像素)在至少一个颜色通道中具有非常低的值,趋近于0。(暗通道图可以估计雾浓度信息) 大气散射...

    本文创新点: 作者提出暗通道先验算法,即:在户外无雾图像中(非天空区域)的大部分局部区域,存在一些像素点(暗像素)在至少一个颜色通道中具有非常低的值,趋近于0。(暗通道图可以估计雾浓度信息)
    大气散射模型
    在计算机视觉中,描述有雾图像的模型可以表示为
    I(x)=J(x)t(x)+A(1t(x))I(x)=J(x)t(x)+A(1-t(x)) (1)
    I(x)I(x)为观察到的有雾图像,J(x)J(x)为无雾图像,AA是大气光值,t(x)t(x)为透射率,表示能够到达计算机系统的没有被散射掉的一部分光。
    去雾的目的是从无雾图像 I(x)I(x)中恢复J(x)J(x)AAt(x)t(x)
    暗通道先验
    作者根据对5000多幅无雾图像的暗通道图数据观察发现:约75%的像素值为0,且90%的像素点具有非常低的值,且集中在[0,16]。由此提出暗通道先验理论,即对于一副无雾图像,其暗通道可以表示为:
    Jdark=minyΩ(x)(minc(r,g,b)Jc(y))J^{dark}=min_{y\in \Omega(x)} (min_{c\in (r,g,b)} J^c(y)) (2)
    Jc(y)J^c(y)表示 JJ的任意一个颜色通道,Ω(x)\Omega(x) 表示在像素点xx的窗口。
    根据暗通道先验理论得出:
    Jdark0J^{dark}\rightarrow 0 (3)
    本文思路:
    1 透射率估计
    假设大气光已知,利用大气光对大气散射模型(1)作归一化处理:
    Ic(x)Ac =t(x)Jc(x)Ac +1t(x) \frac{I^c(x)}{A^c} \ =t(x) \frac{J^c(x)}{A^c}\ +1-t(x) (4)
    假设透射率局部区域Ω(x)\Omega(x) 是恒定不变的,计算式(4)的暗通道,即求最小通道和作最小值滤波:
    minyΩ(x)(mincIc(y)Ac)=t~(x)minyΩ(x)(mincJc(y)Ac)min_{y\in \Omega(x)}(min_c \frac{I^c(y)}{A^c})=\widetilde{t}(x) min_{y\in\Omega(x)}(min_c \frac{J^c(y)}{A^c}) (5)
    在无雾图像中, 由于AcA^c是确定的,根据暗通道先验可得JJ的暗通道为0,即:
    Jdark=minyΩ(x)(mincJc(y)Ac)=0J^{dark}=min_{y\in\Omega(x)}(min_c\frac{J^c(y)}{A^c})=0 (6)
    将式(6)带入(5),可得透射率的表达式:
    t~(x)=1minyΩ(x)(mincIc(y)Ac)\widetilde{t}(x)=1-min_{y\in \Omega(x)}(min_c \frac{I^c(y)}{A^c}) (7)
    为了在景深处保持一定的雾感,使图像看起来更真实自然,本文对式(7)引进一个恒定参数ω(0<ω1)\omega(0<\omega \le1)
    t~(x)=1ωminyΩ(x)(mincIc(y)Ac)\widetilde{t}(x)=1-\omega min_{y\in \Omega(x)}(min_c \frac{I^c(y)}{A^c}) (8)
    本文中,ω\omega=0.95。
    透射率优化
    由于暗通道中使用了最小滤波,因此式(8)得到的透射率含有halo效应和块状效应,为了解决这一问题,作者先后提出了soft-matting和导向滤波的优化算法来优化透射率,值得注意的是,soft-matting算法可以很好地消除halo现象和块状现象,但其时间复杂度大大增加;导向滤波算法时间复杂度较小,但其复原后的图像去雾不彻底,在边缘突变区域仍存在一定程度的雾。
    2大气光估计
    在前人的研究中,把大气光取值定在图像中雾最不透明区域。
    在图像中,雾浓度越低,其暗通道图越暗,像素点值越小;雾浓度越高,其暗通道图越亮,像素点值越大,因此,暗通道图可以较好的反映雾浓度信息。
    本文对于大气光值的选取方法为:先在暗通道图中选出图中前0.1%的像素值最大的像素点(这些像素点通常表示的是雾最不透明的点),这些像素点对应到有雾图像中,选取像素值最高的像素点作为大气光A。
    3图像复原
    根据大气散射模型,将大气光A和t透射率t~(x)\widetilde{t}(x)代入式(1),可得出最终的复原场景:
    J(x)=I(x)Amax(t(x),t0)+AJ(x)= \frac{I(x)-A}{max(t(x),t_0)}+A (9)
    其中,t0t_0值为0.1。由于场景点通常不如大气光亮,因此复原后的图像整体较暗,因此我们需要对复原图像J(x)J(x)进行增强。
    本文算法不足
    1 暗通道中使用了最小滤波器,容易产生块状现象和halo效应。
    2 本文对一些特定图像(含有白色景物和天空区域等)失效。
    3 当大气散射模型失效时,本文算法会失败。
    伪代码

    1  对有雾图像作两次最小滤波min(min(I(x))),求暗通道;
    2  求大气光:range=ceil(img_size*0.001);%取暗原色中最亮的0.1%的点数
                radi_pro=zeros(range,1); %用于记录最亮点内对应图片点象素的三个通道的颜色强度
                A=max(radi_pro)/3;%大气光的值(3个颜色通道的平均值)
    3    求透射率: t = 1 - w*(I^dark./A)
    4    根据大气散射模型复原图像
    

    复原效果
    有雾图像

    在这里插入图片描述
    复原图像

    展开全文
  • 传统的暗通道先验已成功地运用于单一图像去模糊问题,但是,当模糊图像具有显著噪声时,暗通道先验无法对模糊核估计起到作用.因此,得益于分数阶计算能够有效地抑制信号的噪声并对信号的低频部分进行增强,将分数阶计算...
  • 针对暗通道先验在天空区域的失效问题,提出了一种基于亮度模型融合的改进暗通道先验图像去雾算法。首先通过Canny算子分割得到天空区域与非天空区域;其次,利用亮度模拟景深,重构亮度透射率,并通过与暗通道透射率的融合...
  • 基于暗通道先验和场景辐射约束相结合的单图像去雾
  • 结合暗通道先验模型和变分模型的单图像去雾和去噪
  • 基于何凯文博士的参考文献He K, Sun J, Tang X. Single image haze removal using dark channel prior[J]. IEEE CVPR, 2009.所写的基于暗通道先验的去雾算法matlab代码,内有代码,论文,及测试图片
  • 何凯明暗通道先验法 原文 翻译 ppt 及大气光模型论文 和图像除雾相关的资料 Single Image Haze Removal Using Dark Channel Prior Optimized contrast enhancement for realtime image and video dehazing modeling...
  • 基于暗通道先验算法的图像去雾技术已经日益成熟,但是其处理速度慢、天空区域过曝、处理完的图像色彩变暗等缺点也很明显。本文针对这几方面分别提出了求透射率时的优化、纠正天空等明亮区域的错误估计的透射率、采用...
  • 针对以上情况,结合大气散射模型和夜间雾天图像成像特点,提出基于Retinex理论和暗通道先验的去雾算法。首先,根据Retinex理论求得夜间场景的有雾入射光图像和有雾反射光图像;其次,利用暗通道先验得到场景的无雾...
  • 基于暗通道先验条件图像去雾算法

    千次阅读 2019-08-30 17:41:37
    基于暗通道先验条件图像去雾算法 香港大学何凯明博士于2009发表了一篇论文《Single Image Haze Removal Using Dark Channel Prior 》。在文章中,何凯明博士提出了一种简单而有效的图像先验暗通道消除单输入图像雾...

    基于暗通道先验条件图像去雾算法

    原文链接:基于暗通道先验条件图像去雾算法

    香港大学何凯明博士于2009发表了一篇论文《Single Image Haze Removal Using Dark Channel Prior 》。在文章中,何凯明博士提出了一种简单而有效的图像先验暗通道消除单输入图像雾的算法。暗通道先验是一种无雾室外图像的统计。这是基于一个关键的观察——大多数无雾室外图像的局部区域包含一些像素,这些像素在至少一个颜色通道(R,G,B)中的强度非常低。利用这一先验模型,可以直接估计图像中薄雾的厚度,并恢复高质量的无薄雾图像。各种室外雾霾图像的结果表明了该方法的有效性。此外,作为除雾的副产物,还可以获得高质量的深度图。
    在这里插入图片描述
    图1 使用单个图像去除薄雾。
    如图1所示,(a)为输入雾度图像,(b)为除雾后的图像,(c)深度图。
    暗通道先验图像去雾计算流程:
    在这里插入图片描述
    图2 暗通道先验图像去雾计算过程

    1 暗通道图像

    暗通道先验是基于以下对无雾室外图像的观察:在大多数非天空斑块中,至少有一个颜色通道在某些像素处的强度非常低。换句话说,这样一个窗口的最小强度值应该很低。形式上,对于图像j,定义:
    在这里插入图片描述
    其中jc是j的颜色通道,Ω(x)是以x为中心的一个局部窗口。观察表明,除了天空区域,jdark的强度较低,如果j是无雾室外图像,则趋向于零(jdark__>0)。称jdark为j的暗通道,上面的统计观测或知识为暗通道先验。

    2 透射率图像

    透射率公式推导过程在此不再赘述,请参看《Single Image Haze Removal Using Dark Channel Prior》。
    实际上,即使在晴朗的日子里,大气也不是绝对没有任何粒子的。所以,当我们看远处的物体时,雾仍然存在。此外,雾的存在是一个基本线索或人类感知深度。这种现象叫做空中透视。如果我们彻底消除图像中的雾,图像可能会显得不自然,深度感也可能会消失。所以可以选择通过在方程中引入一个常数参数ω(0<ω≤1),为遥远的物体保留一小部分雾:
    在这里插入图片描述
    其中t(x)为透射率,ω(0<ω≤1)为常数参数,Ω(x)是以x为中心的一个局部窗口,Ac全球大气光强度,Ic观察到的强度。

    3 带雾图像去雾

    在这里插入图片描述
    图3 去雾效果
    如图3所示,(a)为输入雾化图像,(b)为透射率图,©为Soft Matting后的精细透射图(d)为最终无雾图像。
    在这里插入图片描述
    典型的t0值为0.1,由于场景的亮度通常不如大气光的亮度,所以去雾后的图像看起来很暗淡。因此,增加了J(X)的曝光以达到最佳效果。如图3(d)是最后恢复的无雾图像。
    在这里插入图片描述
    图4 除雾效果
    如图4所示,顶部:输入模糊图像。中间:恢复无雾图像。底部:深度图。顶部行图像中的红色矩形指示的地方为自动获取大气光强的位置。

    4 matlab基于暗通道先验实现图像去雾

    在这里插入图片描述
    图5 实验图1
    在这里插入图片描述
    图6 实验图2
    在这里插入图片描述
    图7 实验图3
    Matlab 源码:

    clear all
    close all
    clc
     
    w0=0.85;   %0.65  乘积因子用来保留一些雾,1时完全去雾     
    t0=0.1;
     
     
    I=imread('test1.jpg');
     
    Ir = I(:,:,1);
    [h,w,s]=size(I);
     
    min_I=zeros(h,w);
    dark_I = zeros(h,w);
    %下面取得暗影通道图像
    for i=1:h                 
        for j=1:w
            dark_I(i,j)=min(I(i,j,:));
        end
    end
    dark_I = uint8(dark_I);
    img_dark = ordfilt2(dark_I,1,ones(5,5));
     
    Max_dark_channel=double(max(max(img_dark)))  %天空亮度
    dark_channel=double(img_dark);
    t1=1-w0*(dark_channel/Max_dark_channel);   %取得透谢分布率图
    t2=max(t1,t0);
    T=uint8(t1*255); 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    I1=double(I);
    J(:,:,1) = uint8((I1(:,:,1) - (1-t2)*Max_dark_channel)./t2);
    J(:,:,2) = uint8((I1(:,:,2) - (1-t2)*Max_dark_channel)./t2);
    J(:,:,3) =uint8((I1(:,:,3) - (1-t2)*Max_dark_channel)./t2);
    figure,
    set(gcf,'outerposition',get(0,'screensize'));
    subplot(221),imshow(I),title('原始图像');
    subplot(222),imshow(J),title('去雾后的图像');
    subplot(223),imshow(img_dark),title('dark channnel的图形');
    subplot(224),imshow(T),title('透射率t的图形');
    imwrite(J,'wu1.jpg');
    

    在这里插入图片描述
    图8 实验结果1
    在这里插入图片描述
    图9 实验结果2
    在这里插入图片描述
    图10 实验结果3
    欢迎关注微信公众号:FPGA开源工作室
    获取更多学习资料。
    FPGA开源工作室

    展开全文
  • 何凯明2009年的IEEE最佳论文翻译,基于暗通道先验的图像去雾
  • 为提升水下图像的视觉效果, 提出了基于红色暗通道先验(RDCP)和逆滤波的水下图像复原算法。该算法首先简化Jaffe-McGlamery水下光学成像模型, 在此基础上, 利用RDCP消除水下成像过程中后向散射引起的图像雾化效果;...
  • 经典去雾算法-何凯明09年提出暗通道先验去雾(Single Image Haze Removal Using Dark Channel Prior) 暗通道去雾公式:I(x) = f(x)*t(x) + (1 – t(x))*A I(x)为待去雾图像,f(x)为去雾后图像,t(x)为透射率(0,1)...

            经典去雾算法-何凯明09年提出暗通道先验去雾(Single Image Haze Removal Using Dark Channel Prior)

            暗通道去雾公式:I(x) = f(x)*t(x) + (1 – t(x))*A

            I(x)为待去雾图像,f(x)为去雾后图像,t(x)为透射率(0,1),A为大气光成分。

            根据公式,去雾算法可解释为:有雾时,相机获取到的图像为两部分组成,一部分为被拍摄物体发射光线穿过雾霾后的光线,另一部分为大气光被雾霾反射后的光线。被拍摄物体发射光线f(x)通过雾霾,雾霾透射率为t(x),那么被摄物体光线到达相机后值为 f(x)*t(x)。原始大气光值为A,大气光的方向可看做与被拍摄物体光线完全相反,大气光一部分穿过雾霾,一部分被雾霾反射,被反射后的光线值即为A - A*t(x))。

            整个去雾流程如下图:

            有了算法模型,接下来对各个部分进行细化实现,首先计算透射率。暗通道先验法指出,根据大量图像统计,无雾图像RGB通道总有一个通道值趋近于0,也就是说在该像素上透射率t(x)趋近于1,f(x)与I(x)几乎相等。而有雾图像,由于雾霾反射的大气光干扰,RGB通道上的最低值,也就是被反射后的大气光值,通过这个值即可计算出该像素上的透射率,当然实际情况下需要考虑其他因素造成的干扰,一般情况下需要对原始图像通过RGB最低值计算的灰度图进行最低值滤波,根据某个像素相邻范围的最低值计算透射率。

            整个计算透射率过程分为三个步骤:第一步为原始图像通过RGB最低值转灰度图;第二步对灰度图进行最低通道滤波,滤波器半径可选,取滤波器内最低值作为中心像素值;第三部对低通滤波后的灰度图进行高斯低通滤波,主要作用为平滑图像,让边缘过度平缓。

            透射率计算流程如下:

            计算大气光成分值则比较简单,大气光成分A可以看做一个独立的RGB像素,包含三个值。一般是对原始图像RGB通道统计0~255概率分布,分别从各个通道最大值开始取占比大于万分之一的值作为A中对应值。

            完成透射率和大气光成分计算后,接下来就是根据公式计算去雾图像。实际上还需要计算透射率图像,透射率图像通过高斯低通滤波后图像进行反转得到,在实际应用开发上,这一步可以省略,直接通过滤波后灰度图实际上可以看做雾霾反射率图T(x),那么透射率即t(x)=1-T(x)。去雾处理后的图像一般会比较暗,可以通过自动对比度或自动色阶进行优化。

            基于暗通道的去雾对于浓度不均匀的雾霾去雾效果比较好,但计算透射率过程由于最低值滤波和高斯低通滤波随着卷积核增大,消耗计算量也越大,处理单张影像还好,如果用于批量处理影像,显然不太实用。当然在实际应用中,可以考虑减少滤波核半径,牺牲一些效果,提升处理速度。把最低值通滤波和高斯低通滤波核半径都设置为1进行测试。

    测试图像:

    去雾效果图:

    附上代码,java版的实现:

    package tools;
    
    import java.awt.image.BufferedImage;
    import java.util.ArrayList;
    import java.util.List;
    
    public class AutoHazeRemoval {
    	
    	//自动去雾
    	public BufferedImage hazeRemoval(BufferedImage image) {
    		BufferedImage tempImage = new BufferedImage(image.getWidth(), image.getHeight(), image.getType());
    		BufferedImage grayImage = new ImageGray().transferGrayImageByLowest(image);
    		grayImage = new LowestFilter().lowestFilter(grayImage, 1);
    		grayImage = new GaussianBlur().gaussBlur(grayImage, 1, 1);
    		int[] A = extractA(image);
    		for (int i = 0; i < image.getWidth(); i++) {
    			for (int j = 0; j < image.getHeight(); j++) {
    				int rgb = image.getRGB(i, j);
    				double R = (rgb >> 16) & 0xff;
    				double G = (rgb >> 8) & 0xff;
    				double B = rgb & 0xff;
    				//计算透射率
    				int transmissionValue = grayImage.getRGB(i, j) & 0xff;
    				double transmission = 0.95 * (double)transmissionValue / 255;
    				R = (R - A[0] * transmission) / (1 - transmission); 
    				G = (G - A[1] * transmission) / (1 - transmission); 
    				B = (B - A[2] * transmission) / (1 - transmission); 
    				rgb = (255 & 0xff) << 24 | (clamp((int)R) & 0xff) << 16
    						| (clamp((int)G) & 0xff) << 8 | (clamp((int)B) & 0xff);
    				tempImage.setRGB(i, j, rgb);
    			}
    		}
    		return tempImage;
    	}
    	/**
    	 * 通过RGB在0~255频率分布计算全局大气光成分
    	 * 
    	 * @param image
    	 * @return
    	 */
    	public int[] extractA(BufferedImage image) {
    		List<Integer[]> list = generateBinary(image);
    		
    		int[] result = new int[3];
    		result[0] = getValue(image,list.get(0));
    		result[1] = getValue(image,list.get(1));
    		result[2] = getValue(image,list.get(2));
    		return result;
    	}
    
    	private int getValue(BufferedImage image,Integer[] list) {
    		double temp = 0.0;
    		int result = 0;
    		for(int i = 255; i > 0; i--) {
    			double num = list[i] / (double)(image.getWidth()*image.getHeight());
    			temp += num;
    			if (temp >= 0.0001) {
    				result = i;
    				break;
    			}
    		}
    		return result;
    	}
    	/**
    	 * 图像二值化 计算图像RGB值从0~255之间分布
    	 * 
    	 * @return
    	 */
    	public List<Integer[]> generateBinary(BufferedImage image) {
    		List<Integer[]> list = new ArrayList<>();
    		Integer[] rlist = new Integer[256];
    		Integer[] glist = new Integer[256];
    		Integer[] blist = new Integer[256];
    		// 通过循环,往集合里面填充0~255个位置,初始值都为0
    		for (int i = 0; i < 256; i++) {
    			rlist[i] = 0;
    			glist[i] = 0;
    			blist[i] = 0;
    		}
    		for (int i = 0; i < image.getWidth(); i++) {
    			for (int j = 0; j < image.getHeight(); j++) {
    				int rgb = image.getRGB(i, j);
    				int r = (rgb >> 16) & 0xff;
    				int g = (rgb >> 8) & 0xff;
    				int b = rgb & 0xff;
    				rlist[r] = rlist[r] + 1;
    				glist[g] = rlist[g] + 1;
    				blist[b] = rlist[b] + 1;
    			}
    		}
    		list.add(rlist);
    		list.add(glist);
    		list.add(blist);
    		return list;
    	}
    
    	// 判断a,r,g,b值,大于256返回256,小于0则返回0,0到256之间则直接返回原始值
    	private int clamp(int rgb) {
    		if (rgb > 255)
    			return 255;
    		if (rgb < 0)
    			return 0;
    		return rgb;
    	}
    }
    

     

    展开全文
  • 何恺明的暗通道先验(dark channel prior)去雾算法是CV界去雾领域很有名的算法,关于该算法的论文"Single Image Haze Removal Using Dark Channel Prior"一举获得2009年CVPR最佳论文。作者统计了大量的无雾图像,...

    何恺明的暗通道先验(dark channel prior)去雾算法是CV界去雾领域很有名的算法,关于该算法的论文"Single Image Haze Removal Using Dark Channel Prior"一举获得2009年CVPR最佳论文。作者统计了大量的无雾图像,发现一条规律:每一幅图像的每一个像素的RGB三个颜色通道中,总有一个通道的灰度值很低。基于这个几乎可以视作是定理的先验知识,作者提出暗通道先验的去雾算法。

           作者首先介绍去雾模型如下:

    如果你不是CV界新手的话,应该对上式倒背如流,在此不再对上式中的各个参量作详细介绍。对于任意一幅输入图像,定义其暗通道的数学表达式为:

    其中c表示rgb三通道中的某一通道。上式表示在一幅输入图像中,先取图像中每一个像素的三通道中的灰度值的最小值,得到一幅灰度图像,再在这幅灰度图像中,以每一个像素为中心取一定大小的矩形窗口,取矩形窗口中灰度值最小值代替中心像素灰度值(最小值滤波),从而得到该雾天图像的暗通道图像。暗通道图像为灰度图像,通过大量统计并观察发现,暗通道图像的灰度值是很低的,所以将整幅暗通道图像中所有像素的灰度值近似为0,即:

            作者在文中假设大气光A为已知量,以便节省篇幅介绍传输函数的求解方法。在此介绍一种简单普遍的求解大气光的方法:对于任意一幅雾天图像,取其暗通道图像灰度值最大的0.1%的像素点对应于原雾天图像的像素位置的每个通道的灰度值的平均值,从而计算出每个通道的大气光值,即大气光值A是一个三元素向量,每一个元素对应于每一个颜色通道。

           对于成像模型,将其归一化,即两边同时除以每个通道的大气光值:

    假设在图像中一定大小的矩形窗口内,传输函数的值为定值,对上式两边用最小化算子(minimum operators)作最小化运算:

     由于在矩形区域内为定值,故将其拿出运算符外部。由于场景辐射(scene radiance)是无雾图像,将暗通道先验应用于J,则有:

    由于总是正值,则有:

     将上式代入到最小化运算的式子中,即可得到传输函数的估计值为:

    为了防止去雾太过彻底,恢复出的景物不自然,应引入参数,重新定义传输函数为:

            对于求解得到的传输函数,应用滤波器对其进行细化,文章中介绍的方法是软抠图的方法,此方法过程复杂,速度缓慢,因此采用导向滤波对传输函数进行滤波。导向滤波的原理此处不再赘述,其伪代码为:

     上述伪代码中,I表示导向图像(guided image),p为输入图像(input image),q为输出图像(output image),表示均值滤波,r为窗口半径。

            将求解得到的去穷远处大气光值和传输函数代入到大气散射模型之中,反解即可得到恢复的目标图像,即场景辐射。本例使用C++实现上述去雾过程,代码如下:

    复制代码

    #include <iostream>
    #include <stdlib.h>
    #include <time.h>
    #include <cmath>
    #include <algorithm>
    #include <opencv2\opencv.hpp>
    
    using namespace cv;
    using namespace std;
    
    typedef struct Pixel 
    {
        int x, y;
        int data;
    }Pixel;
    
    bool structCmp(const Pixel &a, const Pixel &b) 
    {
        return a.data > b.data;//descending降序
    }
    
    Mat minFilter(Mat srcImage, int kernelSize);
    void makeDepth32f(Mat& source, Mat& output);
    void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon);
    Mat getTransmission_dark(Mat& srcimg, Mat& darkimg, int *array, int windowsize);
    Mat recover(Mat& srcimg, Mat& t, float *array, int windowsize);
    
    int main() 
    {
        string loc = "E:\\darkchannel\\firstopencv\\sky.jpg";
        double scale = 1.0;
        string name = "forest";
        clock_t start, finish;
        double duration;
    
        cout << "A defog program" << endl
            << "----------------" << endl;
    
        Mat image = imread(loc);
        Mat resizedImage;
        int originRows = image.rows;
        int originCols = image.cols;
    
        if (scale < 1.0) 
        {
            resize(image, resizedImage, Size(originCols * scale, originRows * scale));
        }
        else 
        {
    
            scale = 1.0;
            resizedImage = image;
        }
    
        int rows = resizedImage.rows;
        int cols = resizedImage.cols;
        Mat convertImage;
        resizedImage.convertTo(convertImage, CV_32FC3, 1 / 255.0, 0);
        int kernelSize = 15 ? max((rows * 0.01), (cols * 0.01)) : 15 < max((rows * 0.01), (cols * 0.01));
        //int kernelSize = 15;
        int parse = kernelSize / 2;
        Mat darkChannel(rows, cols, CV_8UC1);
        Mat normalDark(rows, cols, CV_32FC1);
        int nr = rows;
        int nl = cols;
        float b, g, r;
        start = clock();
        cout << "generating dark channel image." << endl;
        if (resizedImage.isContinuous()) 
        {
            nl = nr * nl;
            nr = 1;
        }
        for (int i = 0; i < nr; i++) 
        {
            float min;
            const uchar* inData = resizedImage.ptr<uchar>(i);
            uchar* outData = darkChannel.ptr<uchar>(i);
            for (int j = 0; j < nl; j++) 
            {
                b = *inData++;
                g = *inData++;
                r = *inData++;
                min = b > g ? g : b;
                min = min > r ? r : min;
                *outData++ = min;
            }
        }
        darkChannel = minFilter(darkChannel, kernelSize);
    
        imshow("darkChannel", darkChannel);
        cout << "dark channel generated." << endl;
    
        //estimate Airlight
        //开一个结构体数组存暗通道,再sort,取最大0.1%,利用结构体内存储的原始坐标在原图中取点
        cout << "estimating airlight." << endl;
        rows = darkChannel.rows, cols = darkChannel.cols;
        int pixelTot = rows * cols * 0.001;
        int *A = new int[3];
        Pixel *toppixels, *allpixels;
        toppixels = new Pixel[pixelTot];
        allpixels = new Pixel[rows * cols];
    
    
        for (unsigned int r = 0; r < rows; r++) 
        {
            const uchar *data = darkChannel.ptr<uchar>(r);
            for (unsigned int c = 0; c < cols; c++) 
            {
                allpixels[r*cols + c].data = *data;
                allpixels[r*cols + c].x = r;
                allpixels[r*cols + c].y = c;
            }
        }
        std::sort(allpixels, allpixels + rows * cols, structCmp);
    
        memcpy(toppixels, allpixels, pixelTot * sizeof(Pixel));
    
        float A_r, A_g, A_b, avg, maximum = 0;
        int idx, idy, max_x, max_y;
        for (int i = 0; i < pixelTot; i++) 
        {
            idx = allpixels[i].x; idy = allpixels[i].y;
            const uchar *data = resizedImage.ptr<uchar>(idx);
            data += 3 * idy;
            A_b = *data++;
            A_g = *data++;
            A_r = *data++;
            //cout << A_r << " " << A_g << " " << A_b << endl;
            avg = (A_r + A_g + A_b) / 3.0;
            if (maximum < avg) 
            {
                maximum = avg;
                max_x = idx;
                max_y = idy;
            }
        }
    
        delete[] toppixels;
        delete[] allpixels;
    
        for (int i = 0; i < 3; i++) 
        {
            A[i] = resizedImage.at<Vec3b>(max_x, max_y)[i];
        }
        cout << "airlight estimated as: " << A[0] << ", " << A[1] << ", " << A[2] << endl;
        //cout << endl;
    
        //暗通道归一化操作(除A)
        //(I / A)
        cout << "start normalization of input image I." << endl;
        float tmp_A[3];
        tmp_A[0] = A[0] / 255.0;
        tmp_A[1] = A[1] / 255.0;
        tmp_A[2] = A[2] / 255.0;
        for (int i = 0; i < nr; i++) 
        {
            float min = 1.0;
            const float* inData = convertImage.ptr<float>(i);
            float* outData = normalDark.ptr<float>(i);
            for (int j = 0; j < nl; j++) 
            {
                b = *inData++ / tmp_A[0];
                g = *inData++ / tmp_A[1];
                r = *inData++ / tmp_A[2];
                min = b > g ? g : b;
                min = min > r ? r : min;
                *outData++ = min;
            }
        }
        cout << "normalization finished." << endl << "generating relative dark channel image." << endl;
        //暗通道最小滤波
        normalDark = minFilter(normalDark, kernelSize);
        cout << "dark channel image generated." << "start estimating transmission and guided image filtering." << endl;
        imshow("normal", normalDark);
        int kernelSizeTrans = std::max(3, kernelSize);
        //求t与将t进行导向滤波
    
        Mat trans = getTransmission_dark(convertImage, normalDark, A, kernelSizeTrans);
        cout << "tansmission estimated and guided filtered." << endl;
        imshow("filtered t", trans);
        cout << "start recovering." << endl;
        Mat finalImage = recover(convertImage, trans, tmp_A, kernelSize);
        cout << "recovering finished." << endl;
        Mat resizedFinal;
        if (scale < 1.0) 
        {
            resize(finalImage, resizedFinal, Size(originCols, originRows));
            imshow("final", resizedFinal);
        }
        else 
        {
            imshow("final", finalImage);
        }
        finish = clock();
        duration = (double)(finish - start);
        cout << "defog used " << duration << "ms time;" << endl;
        waitKey(0);
    
        finalImage.convertTo(finalImage, CV_8UC3, 255);
        imwrite(name + "_refined.png", finalImage);
        destroyAllWindows();
        image.release();
        resizedImage.release();
        convertImage.release();
        darkChannel.release();
        trans.release();
        finalImage.release();
        return 0;
    }
    
    Mat minFilter(Mat srcImage, int kernelSize) 
    {
        int radius = kernelSize / 2;
    
        int srcType = srcImage.type();
        int targetType = 0;
        if (srcType % 8 == 0) 
        {
            targetType = 0;
        }
        else 
        {
            targetType = 5;
        }
        Mat ret(srcImage.rows, srcImage.cols, targetType);
        Mat parseImage;
        copyMakeBorder(srcImage, parseImage, radius, radius, radius, radius, BORDER_REPLICATE);
        for (unsigned int r = 0; r < srcImage.rows; r++) 
        {
            float *fOutData = ret.ptr<float>(r);
            uchar *uOutData = ret.ptr<uchar>(r);
            for (unsigned int c = 0; c < srcImage.cols; c++) 
            {
                Rect ROI(c, r, kernelSize, kernelSize);
                Mat imageROI = parseImage(ROI);
                double minValue = 0, maxValue = 0;
                Point minPt, maxPt;
                minMaxLoc(imageROI, &minValue, &maxValue, &minPt, &maxPt);
                if (!targetType) 
                {
                    *uOutData++ = (uchar)minValue;
                    continue;
                }
                *fOutData++ = minValue;
            }
        }
        return ret;
    }
    
    void makeDepth32f(Mat& source, Mat& output)
    {
        if ((source.depth() != CV_32F) > FLT_EPSILON)
            source.convertTo(output, CV_32F);
        else
            output = source;
    }
    
    void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon)
    {
        CV_Assert(radius >= 2 && epsilon > 0);
        CV_Assert(source.data != NULL && source.channels() == 1);
        CV_Assert(guided_image.channels() == 1);
        CV_Assert(source.rows == guided_image.rows && source.cols == guided_image.cols);
    
        Mat guided;
        if (guided_image.data == source.data)
        {
            //make a copy
            guided_image.copyTo(guided);
        }
        else
        {
            guided = guided_image;
        }
    
        //将输入扩展为32位浮点型,以便以后做乘法
        Mat source_32f, guided_32f;
        makeDepth32f(source, source_32f);
        makeDepth32f(guided, guided_32f);
    
        //计算I*p和I*I
        Mat mat_Ip, mat_I2;
        multiply(guided_32f, source_32f, mat_Ip);
        multiply(guided_32f, guided_32f, mat_I2);
    
        //计算各种均值
        Mat mean_p, mean_I, mean_Ip, mean_I2;
        Size win_size(2 * radius + 1, 2 * radius + 1);
        boxFilter(source_32f, mean_p, CV_32F, win_size);
        boxFilter(guided_32f, mean_I, CV_32F, win_size);
        boxFilter(mat_Ip, mean_Ip, CV_32F, win_size);
        boxFilter(mat_I2, mean_I2, CV_32F, win_size);
    
        //计算Ip的协方差和I的方差
        Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
        Mat var_I = mean_I2 - mean_I.mul(mean_I);
        var_I += epsilon;
    
        //求a和b
        Mat a, b;
        divide(cov_Ip, var_I, a);
        b = mean_p - a.mul(mean_I);
    
        //对包含像素i的所有a、b做平均
        Mat mean_a, mean_b;
        boxFilter(a, mean_a, CV_32F, win_size);
        boxFilter(b, mean_b, CV_32F, win_size);
    
        //计算输出 (depth == CV_32F)
        output = mean_a.mul(guided_32f) + mean_b;
    }
    
    Mat getTransmission_dark(Mat& srcimg, Mat& darkimg, int *array, int windowsize)
    {
        //t = 1 - omega * minfilt(I / A);
        float avg_A;
        //convertImage是一个CV_32FC3的图
        if (srcimg.type() % 8 == 0) {
            avg_A = (array[0] + array[1] + array[2]) / 3.0;
        }
        else {
            avg_A = (array[0] + array[1] + array[2]) / (3.0 * 255.0);
        }
        float w = 0.95;
        int radius = windowsize / 2;
        int nr = srcimg.rows, nl = srcimg.cols;
        Mat transmission(nr, nl, CV_32FC1);
    
        for (int k = 0; k<nr; k++) 
        {
            const float* inData = darkimg.ptr<float>(k);
            float* outData = transmission.ptr<float>(k);
            float pix[3] = { 0 };
            for (int l = 0; l<nl; l++)
            {
                *outData++ = 1.0 - w * *inData++;
            }
        }
        imshow("t", transmission);
    
        Mat trans(nr, nl, CV_32FC1);
        Mat graymat(nr, nl, CV_8UC1);
        Mat graymat_32F(nr, nl, CV_32FC1);
    
        if (srcimg.type() % 8 != 0) 
        {
            cvtColor(srcimg, graymat_32F, CV_BGR2GRAY);
            guidedFilter(transmission, graymat_32F, trans, 6 * windowsize, 0.001);
        }
        else 
        {
            cvtColor(srcimg, graymat, CV_BGR2GRAY);
    
            for (int i = 0; i < nr; i++) 
            {
                const uchar* inData = graymat.ptr<uchar>(i);
                float* outData = graymat_32F.ptr<float>(i);
                for (int j = 0; j < nl; j++)
                    *outData++ = *inData++ / 255.0;
            }
            guidedFilter(transmission, graymat_32F, trans, 6 * windowsize, 0.001);
        }
        return trans;
    }
    
    Mat recover(Mat& srcimg, Mat& t, float *array, int windowsize)
    {
        //J(x) = (I(x) - A) / max(t(x), t0) + A;
        int radius = windowsize / 2;
        int nr = srcimg.rows, nl = srcimg.cols;
        float tnow = t.at<float>(0, 0);
        float t0 = 0.1;
        Mat finalimg = Mat::zeros(nr, nl, CV_32FC3);
        float val = 0;
    
        //Be aware that transmission is a grey image
        //srcImg is a color image
        //finalImg is a color image
        //Mat store color image a pixel per 3 position
        //store grey image a pixel per 1 position
        for (unsigned int r = 0; r < nr; r++) {
            const float* transPtr = t.ptr<float>(r);
            const float* srcPtr = srcimg.ptr<float>(r);
            float* outPtr = finalimg.ptr<float>(r);
            for (unsigned int c = 0; c < nl; c++) {
                //transmission image is grey, so only need 
                //to move once per calculation, using index 
                //c(a.k.a. columns) to move is enough 
                tnow = *transPtr++;
                tnow = std::max(tnow, t0);
                for (int i = 0; i < 3; i++) {
                    //so to calculate every color channel per pixel
                    //move the ptr once after one calculation.
                    //after 3 times move, calculation for a pixel is done
                    val = (*srcPtr++ - array[i]) / tnow + array[i];
                    *outPtr++ = val + 10 / 255.0;
                }
            }
        }
        cout << finalimg.size() << endl;
        return finalimg;
    }

    复制代码

            运行代码的测试结果为:

            可以很明显地观察到,所选素材中含有大量的天空区域。利用暗通道先验规律对含有大量天空区域的图像作暗通道图像运算时,会发现天空区域的暗通道的灰度值并非趋向于零,因此造成在天空区域的传输函数估计值很低,造成误差。为了弥补这个误差,作者对传输函数的值设定了一个下界,这个下界通常取0.1,然而设定这个下界后,会发现天空区域如上述去雾后的结果那样出现光圈(halo)效应,天空部分的去雾效果很不自然。这也是暗通道先验去雾算法的劣处所在。

            

             基于暗通道先验的去雾算法所得到的去雾结果令人惊叹,其理论过程简单易行,至今仍是去雾领域的经典算法。近年来,许多学者对该算法不断从技术上改进,也有学者试图从别的角度解释暗通道先验规律的成立原因,下面从几何的角度解释暗通道先验去雾算法之所以成立的原因。为了与先前介绍的暗通道先验算法中的符号区分,接下来的文中出现的特定的符号的意义如下:表示图像中特定大小的矩阵窗口内的像素集合;表示将对应于中的像素的灰度值在RGB空间中表示的坐标的集合。

           本例介绍的论文为“Single Image Dehazing Using Color Ellipsoid Prior”,作者在这篇文章中介绍了暗通道先验去雾算法的几何数学解释过程。

       无论是何种去雾方法,将去雾前后图像的同一片区域投射在RGB坐标空间中,形成如上图所示的点坐标的分布。将去雾前的坐标和去雾后的坐标相连,形成RGB空间的向量,将不同的去雾方法对应得到的向量定义为先验向量(prior vector)。上图中,通过先验向量P去雾扩大的无雾坐标分布要比Q扩大的要小,这意味着通过先验向量P去雾得到的结果保持了饱和度,却丧失了对比度。反之,通过先验向量Q去雾得到的结果保持了对比度,却产生了过度饱和的现象,所以通过先验向量Q的去雾结果中的树叶部分出现黑点(dark spots)。所有的去雾算法归结为设定合适的先验向量,求解合适的传输函数,满足:尽可能地增大图像对比度,同时抑制过度饱和。

           上图表示雾天图像两个相近区域的像素灰度值分布情况,其中一区域不含有不规则像素值(irregular pixel values),另一区域则含有不规则像素值。不规则像素值的分布是零星的,两个区域的大部分像素值的分布区域几乎一致,所以两个区域的去雾结果大体是相同的。然而在选取先验向量时,若从不规则像素值区域选取的话,则会产生很大的误差,因此设计的先验向量应该具有数据鲁棒性(statistically robust)以克服不规则像素值的干扰。

            如上图所示,雾天图像的局部区域中含有树叶部分和背景部分,树叶部分和背景部分的像素灰度值的RGB空间中如图所示标出。在两个部分的边缘地带,有大量的像素灰度值是由两个区域混合的,这就造成选取混合区域的先验向量进行去雾时,会在边缘地带造成光圈效应(halo effect)。后续算法会解释如何去除这种光圈效应。

              下面介绍颜色椭球先验(color ellipsoid prior):在RGB三维空间中,含雾图像的向量往往聚集在一起并有序的排列着。在三维空间中,向量可以用空间中的坐标表示,所以向量和颜色空间中的点是等价的。向量聚集区域可以统计地看做是一个椭球,椭球的形状由向量的聚集状态决定。这种方法的缺点是无法统计随意分布的向量区域。

       如图(a)和图(b)所示,其为图像区域对应的向量分布区域的示意图。前者为没有随机分布点的规则椭球向量分布,后者是含有不规则向量的椭球向量分布,但是两者可以视为同等椭球向量分布。

           现在根据椭球模型构造数学表达式,设为RGB空间的向量变量,对应于图像中像素区域(i为序数标记)的向量椭球区域可以表示为:

     

     其中,表示雾天图像在区域内的灰度值的均值对应的向量,其表达式为:

    其中为雾天图像的归一化像素值,为三维向量。为无穷远处的大气光值,此处假设其已知:

    方差定义为:

    在RGB空间中,均值决定了球的中心位置,方差决定了球的方向和大小。

            上述各式为椭球模型的基本定义的参量,现在分析CEP值和传输函数的计算过程。文中有一段文字“The proposed method calculates the vector having the smallest color component among the vectors in the ellipsoid  to maximally stretch the hazy vector cluster while preventing the color components of the majority of streched vectors from being negative”,这段话的意思是,通过先验向量将原来的图像雾天的像素灰度值分布的椭球拉伸映射到灰度值很低的椭球中(无雾图像的像素灰度值相对于原雾天图像是比较低的)。故在设计合适的先验向量和先验值(本例中为CEP值)去雾时,应满足如前所述的条件。通过对大气散射模型反解,可以得到目标辐射景象为:

    为了使得去雾图像,即目标辐射场景的灰度值尽可能低,则传输函数的值应该尽可能地大,又因为传输函数与先验值的关系为,所以先验值应该尽可能地小。所以本例中定义的先验值,即CEP值的数学表达式为:

    该式表示,在椭球中的每一个向量取其三颜色通道分量的最小值,再取这些最小值当中的最小值作为CEP值(这与暗通道先验规律具有一致性)。

           下面从几何的角度求解CEP值:首先定义三个坐标轴的单位向量分别为:;将定义为与单位向量垂直且通过坐标原点的平面;如此,则为向量的R,G,B某一分量的值,并且该数值表示向量对应于空间中的坐标(即在椭球面上或者椭球内的点)到平面的距离;设为使得距离最小的向量,即:

    之所以有标记c,是因为表示椭球面上或椭球内部中到平面距离最小的点对应的向量,中的c的取值有三个:,所以对应的中c的取值也有三个且和中c的取值一致,分别为,三者分别表示椭球中到平面距离最小的点坐标对应的向量。

             如上所示的(a),(b),(c)分别为对应在空间中的坐标点和其到平面的距离的示意图。容易观察到,到三个平面距离最近的点均在椭球面上,且在点处椭球的切面与对应的平面平行。由于点在椭球面上,则满足表达式:

    上式表示点在以为中心,为协方差矩阵(该协方差矩阵规定了椭球的形状和方向)逆矩阵的椭球面上。则关于CEP值,有如下推导:

     由于点处椭球的切面与对应的平面平行,所以在点处的外单位法向量(outward normal vector)为上图所示的。根据文中所述定理:The gradient direction of a vector on the ellipse surface is perpendicular to the tangent plane,and so the outward normal vector must have the same direction as the gradient of the ellipse surface,即在球面处的点的梯度方向和在该点处的切平面垂直,在该点处的外单位法向量和梯度方向是一致的。根据该定理,可得下式:

    其中,梯度部分的计算过程为:

    则有:

    所以有:

    其中,。由此重新表达为:

    所以中的向量到平面距离最小值为:

    所以,CEP值的表达式为:

    故对应的传输函数的表达式为:

            下面根据上述推导过程介绍如何设计算法,设是使得CEP值在三通道中最小的那一通道,即:

    则传输函数可以重新表示为:

    若直接计算每一点邻域矩形的均值与方差的差值,再取其三通道中的最小值较为繁琐,作者提出一种近似方法,利用这种近似方法与原原本本的计算结果相差无几,但是大大提高了速度和效率。近似方法为:

    其中,表示输入雾天图像在像素点j处c通道的灰度值。

           设大气光已知,按照上述过程计算每一像素点处的均值和方差,并将其代入传输函数公式中以计算传输函数,由此即可恢复目标图像,处理结果为:

    上述结果在边缘部分出现了很模糊的现象,这是由于在选择椭球区域时,所选区域的像素梯度过大,某一点的像素灰度值与其周围像素的灰度值差别很大。因此有必要对求取的传输函数进行细化,作者分别对均值和方差进行细化(即滤波):对于均值部分,作者采用了与导向滤波十分类似的滤波器,但不完全是导向滤波;对于方差部分,作者采用了均值滤波。滤波过程的伪代码如下:

     

     

           用MATLAB对上述伪代码的仿真为:

    复制代码

    clear all;
    I=imread('E:\picture\t.bmp');
    I=double(I);
    I=I./255;
    dark=darkfunction(I);
    [m,n,o]=size(I);
    imsize=m*n;
    numpx=floor(imsize/1000);
    J_dark_vec=reshape(dark,imsize,1);
    I_vec=reshape(I,imsize,3);
    [J_dark_vec,indices]=sort(J_dark_vec);
    indices=indices(imsize-numpx+1:end);
    atmSum=zeros(1,3);
    for ind=1:numpx
        atmSum=atmSum+I_vec(indices(ind),:);
    end;
    atmospheric=atmSum/numpx;
    A=atmospheric;
    
    omega=0.95;
    r=10;epsilon=0.002;
    epsilon=(1/mean(A))^2*epsilon;
    R=I(:,:,1)/A(1);G=I(:,:,2)/A(2);B=I(:,:,3)/A(3);
    Imin=min(min(R,G),B);
    
    N=boxfilter(ones(m,n),r);
    u_m=boxfilter(Imin,r)./N;
    %对均值滤波,与导向滤波相似,但不是导向滤波
    u_mm=boxfilter(Imin.*Imin,r)./N;
    sigma_m=u_mm-u_m.*u_m;
    a=sigma_m./(sigma_m+epsilon);
    b=(1-a).*u_m;
    u_a=boxfilter(a,r)./N;
    u_b=boxfilter(b,r)./N;
    u_m=u_a.*Imin+u_b;
    
    sigma_mm=boxfilter((Imin-u_m).*(Imin-u_m),r)./N;
    %对方差滤波,此处为均值滤波,也可仿上述滤波过程进行滤波
    sigma_mm=boxfilter(sigma_mm,r)./N;
    
    t=1-omega*(u_m-sqrt(abs(sigma_mm)));
    
    for i=1:m
        for j=1:n
            for k=1:3
                I(i,j,k)=(I(i,j,k)-A(k))/max(t(i,j),0.1)+A(k);
            end;
        end;
    end;
    imshow(I);
    
    function [dark] =darkfunction(I)
    R=I(:,:,1);
    G=I(:,:,2);
    B=I(:,:,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);
    end
    
    function [J] = boxfilter(I,r)
    [m,n]=size(I);
    J=zeros(size(I));
    Icum=cumsum(I,1);
    J(1:r+1,:)=Icum(r+1:2*r+1,:);
    J(r+2:m-r,:)=Icum(2*r+2:m,:)-Icum(1:m-2*r-1,:);
    J(m-r+1:m,:)=repmat(Icum(m,:),[r,1])-Icum(m-2*r:m-r-1,:);
    Icum=cumsum(J,2);
    J(:,1:r+1)=Icum(:,r+1:2*r+1);
    J(:,r+2:n-r)=Icum(:,2*r+2:n)-Icum(:,1:n-2*r-1);
    J(:,n-r+1:n)=repmat(Icum(:,n),[1,r])-Icum(:,n-2*r:n-r-1);
    end

    复制代码

            用C++实现算法为:

    复制代码

    #include <iostream>
    #include <stdlib.h>
    #include <time.h>
    #include <cmath>
    #include <algorithm>
    
    #include <opencv2\opencv.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    #include <opencv2\highgui\highgui.hpp>
    
    using namespace std;
    using namespace cv;
    
    int main()
    {
        Mat GF_smooth(Mat& src, int s, double epsilon);
        Mat staticMin(Mat& I, int s, double eeps, double alpha);
        int est_air(Mat& src, int s, double *A_r, double *A_g, double *A_b);
        Mat est_trans_fast(Mat& src, int s, double eeps, double k, double A_r, double A_g, double A_b);
        Mat rmv_haze(Mat& src, Mat& t, double A_r, double A_g, double A_b);
    
        string loc = "E:\\sphere\\firstopencv\\t.bmp";
        string name = "forest";
        clock_t start, finish;
        double duration;
    
        cout << "A defog program" << endl
             << "---------------" << endl;
    
        start = clock();
            Mat image = imread(loc);
        imshow("hazyimage", image);
        cout << "input hazy image" << endl;
        
        int s = 15;
        double eeps = 0.002, omega = 0.95;
        double air_r, air_g, air_b;
        est_air(image, s, &air_r, &air_g, &air_b);
        cout << "airlight estimation as:" << air_r << ", " << air_g << ", " << air_b << endl;
    
        Mat t;
        t=est_trans_fast(image, s, eeps, omega, air_r, air_g, air_b);
    
        Mat dehazedimage;
        dehazedimage = rmv_haze(image, t, air_r, air_g, air_b);
        imshow("dehaze", dehazedimage);
        finish = clock();
        duration = (double)(finish - start);
        cout << "defog used" << duration << "ms time;" << endl;
        waitKey(0);
    
        imwrite(name + "_refined.png", dehazedimage);
        destroyAllWindows();
        image.release();
        dehazedimage.release();
        return 0;
    }
    
    //---------------------- GUIDED FILTER -------------------//
    Mat GF_smooth(Mat& src, int s, double epsilon)
    {
        src.convertTo(src, CV_64FC1);
    
        Mat mean_I;
        blur(src, mean_I, Size(s, s), Point(-1, -1));
    
        Mat II = src.mul(src);
        Mat var_I;
        blur(II, var_I, Size(s, s), Point(-1, -1));
        var_I = var_I - mean_I.mul(mean_I);
    
        Mat a = var_I / ((var_I + epsilon));
        Mat b = mean_I - a.mul(mean_I);
    
        Mat mean_a;
        blur(a, mean_a, Size(s, s), Point(-1, -1));
        Mat mean_b;
        blur(b, mean_b, Size(s, s), Point(-1, -1));
    
        return mean_a.mul(src) + mean_b;
    }
    
    
    Mat staticMin(Mat& I, int s, double eps, double alpha)
    {
        Mat mean_I = GF_smooth(I, s, eps);
    
        Mat var_I;
        blur((I - mean_I).mul(I - mean_I), var_I, Size(s, s), Point(-1, -1));
    
        Mat mean_var_I;
        blur(var_I, mean_var_I, Size(s, s), Point(-1, -1));
    
        Mat z_I;
        sqrt(mean_var_I, z_I);
    
        return mean_I - alpha*z_I;
    }
    
    
    //---------------------- DEHAZING FUNCTIONS -------------------//
    int est_air(Mat& src, int s, double *A_r, double *A_g, double *A_b)
        {
            /// Split RGB to channels
            src.convertTo(src, CV_64FC3);
            vector<Mat> channels(3);
            split(src, channels);        // separate color channels
    
            Mat R = channels[2];
            Mat G = channels[1];
            Mat B = channels[0];
    
            Mat Im = min(min(R, G), B);
    
            /// Estimate airlight
            Mat blur_Im;
            blur(Im, blur_Im, Size(s, s), Point(-1, -1));
    
            int maxIdx[2] = { 0, 0 };
            minMaxIdx(blur_Im, NULL, NULL, NULL, maxIdx);
    
            int width = R.cols;
            *A_r = ((double*)R.data)[maxIdx[0] * R.cols + maxIdx[1]];
            *A_g = ((double*)G.data)[maxIdx[0] * R.cols + maxIdx[1]];
            *A_b = ((double*)B.data)[maxIdx[0] * R.cols + maxIdx[1]];
    
            return 0;
        }
    
        Mat est_trans_fast(Mat& src, int s, double eeps, double k, double A_r, double A_g, double A_b)
        {
            /// Split RGB to channels
            src.convertTo(src, CV_64FC3);
            vector<Mat> channels(3);
            split(src, channels);        // separate color channels
    
            Mat R = channels[2];
            Mat G = channels[1];
            Mat B = channels[0];
    
            /// Estimate transmission
            Mat R_n = R / A_r;
            Mat G_n = G / A_g;
            Mat B_n = B / A_b;
    
            Mat Im = min(min(R_n, G_n), B_n);
    
            eeps = (3 * 255 / (A_r + A_g + A_b))*(3 * 255 / (A_r + A_g + A_b))*eeps;
            double alpha = 2;
            Mat z_Im = staticMin(Im, s, eeps, alpha);
    
            return min(max(0.001, 1 - k*z_Im), 1);
        }
    
    
        Mat rmv_haze(Mat& src, Mat& t, double A_r, double A_g, double A_b)
        {
            /// Split RGB to channels
            src.convertTo(src, CV_64FC3);
            vector<Mat> channels(3);
            split(src, channels);        // separate color channels
    
            Mat R = channels[2];
            Mat G = channels[1];
            Mat B = channels[0];
    
            /// Remove haze
            channels[2] = (R - A_r) / t + A_r;
            channels[1] = (G - A_g) / t + A_g;
            channels[0] = (B - A_b) / t + A_b;
    
            Mat dst;
            merge(channels, dst);
    
            dst.convertTo(dst, CV_8UC3);
            return dst;
        }

    复制代码

            用MATLAB或者C++运行的实验结果为:

            结果分析:通过该实验结果发现,颜色椭球先验(CEP)的去雾算法的运行结果和暗通道先验(DCP)算法的运行结果基本上是一致的,这也验证了CEP是DCP在几何上的解释,两者具有内在一致性,在本质上是相同的。但是,两者也有细微的差别。由于两种算法的实现结果一致,不好比较,现在将两者在传输函数优化之前的去雾效果进行比较,下图为DCP所求的传输函数在导向滤波之前的去雾效果:

    下图为CEP所求传输函数在滤波优化前的去雾效果:

    很明显地观察到:DCP的边缘效应比CEP的边缘效应更明显一些,说明基于CEP计算的传输函数要比基于DCP计算的传输函数更精确一些。这也体现了文章中所说的,与之前的DCP相比较,CEP去雾效果更具有数据鲁棒性(statistically robust)

    展开全文
  • 针对雾霾天气条件下,大气粒子的散射作用导致的图像质量下降问题,提出一种基于暗通道先验知识与局部多项式核回归算法相结合的去雾方法。根据暗通道先验原理估计出大气光强度和初始透射率,采用局部多项式核回归对...
  • 基于暗通道先验条件图像去雾算法香港大学何凯明博士于2009发表了一篇论文《Single Image Haze Removal Using Dark Channel Prior 》。在文章中,何凯明博士提出了一种简单而有效的图像先验暗通道消除单输入图像雾的...
  • Opencv实现暗通道先验去雾算法

    千次阅读 2017-10-03 10:58:00
    Opencv实现暗通道先验去雾算法今天读了何凯明博士的《Single Image Haze Removal Using Dark Channel Prior》,用opencv实现了一遍。 其中,暗通道及T(x)采用腐蚀erode算法。 具体代码如下:class ...
  • 基于导向滤波的暗通道先验去雾算法(Python语言,可直接运行)1 编译环境2 原理介绍2.1 暗通道先验2.1.1 暗通道先验理论与去雾模型2.1.2 处理步骤2.2 导向滤波求t(x)3 代码4 处理结果4.1 第一组结果及其比较,暗通道...
  • 针对目前检测方法适应性不强、在复杂环境下检测性能不高的问题,提出了一种基于背景动态更新和暗通道先验的烟雾检测算法。算法首先通过改进的背景动态更新算法提取运动前景;然后,结合暗通道先验知识确定前景中的...
  • 本文主要讲解,何恺明暗通道先验去雾以及加权最小二乘算法的原理部分,对应代码可在资源处下载。1暗通道先验去雾模型 公式(1)是常用的雾霭形成模型: (1)其中,I(x)为实际拍摄到的有雾视频帧,J(x)为去雾后的原...
  • 暗通道先验去雾 辐射图解 模糊图像模型: 问题? I 是已知的,目标恢复 J ,而 t , A 未知。工作:尽可能得到最优的 t , A 值。 什么是暗通道,暗像素? 在多数非天空的图像局部区域,存在某些像素(暗像素)在...
  • 暗通道先验去雾

    2016-01-16 12:54:34
    经过自己的优化,运行速度较快,含相关报告和文献

空空如也

空空如也

1 2 3 4 5 ... 12
收藏数 224
精华内容 89
关键字:

暗通道先验