精华内容
下载资源
问答
  • 暗通道先验算法

    万次阅读 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    根据大气散射模型复原图像
    

    复原效果
    有雾图像

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

    展开全文
  • 在http://write.blog.csdn.net/postedit/78301999中介绍了图像去雾的相关研究方法,发现目前为止在图像去雾方面,何凯明博士基于暗通道先验算法具有很好的效果,关于该方法的介绍也很多,本篇作下学习笔记和个人...

    http://write.blog.csdn.net/postedit/78301999中介绍了图像去雾的相关研究方法,发现目前为止在图像去雾方面,何凯明博士基于暗通道先验的算法具有很好的效果,关于该方法的介绍也很多,本篇作下学习笔记和个人理解:

    Retinex背景知识

    Atmospheric Scattering Model

    (1)

    图像去雾过程就是根据 I(x) 求解 J(x) 的过程。从上面的公式可以看出,基于物理模型的去雾算法本质是根据已知的有雾图像 I(x) 求出透射率 t(x) 和全局大气光成分A (A>0),进而得到复原图像J(X)(场景光成分)。在参考【2】里有关于该模型的详细解释,其中J(x)t(x)是直接衰减(direct attenuation),A(1-t(x)) 是空气光(airlight)。当大气是同质(说的应该是均匀分布吧)时,传输率t表示为:

                                                                                                                                               (2)

    其中,β是大气的散射系数,d是场景光线的景深。

    有雾图像特征:



    上面几幅图是图像的彩色分布图和直方图,一般的,有雾图像的RGB色彩分布图比较集中,直方图范围也窄,相对集中;清晰图像的RGB色彩分布图比较广泛,直方图分布均匀,色彩感较强。

    基于暗通道先验去雾算法原理:

    下面是原作者的一些感言:

    简单有效的图像去雾技术

    这篇论文研究的问题是图像的去雾技术,它可以还原图像的颜色和能见度,同时也能利用雾的浓度来估计物体的距离,这些在计算机视觉上都有重要应用(例如三维重建,物体识别)。但是之前人们还没找到简单有效的方法来达到这个目的。在这篇论文里,我们找到了一个非常简单的,甚至说令人惊讶统计规律,并提出了有效的去雾方法。
    与之前的方法不同,我们把注意力放到了无雾图像的统计特征上。我们发现,在无雾图像中,每一个局部区域都很有可能会有阴影,或者是纯颜色的东西,又或者是黑色的东西。因此,每一个局部区域都很有可能有至少一个颜色通道会有很低的值。我们把这个统计规律叫做Dark Channel Prior。直观来说,Dark Channel Prior认为每一个局部区域都总有一些很暗的东西。这个规律很简单,但在我们研究的去雾问题上却是本质的基本规律。
    由于雾总是灰白色的,因此一旦图像受到雾的影响,那么这些本来应该很暗的东西就会变得灰白。不仅如此,根据物理上雾的形成公式,我们还能根据这些东西的灰白程度来判断雾的浓度。因此,我们提出的Dark Channel Prior能很有效地去除雾的影响,同时利用物的浓度来估算物体的距离。

    大道之行在于简

    我们这篇文章的三个审稿人都给出了最高的评分。他们认为我们的方法简单而有效。其中一位评委说,Dark Channel Prior的想法听起来很不可思议,但我们却证明了其真实性。另一位评委认为很少有文章能够用如此简单的方法使实验结果获得如此大的提升。还有一位评委甚至亲自实现了我们的方法并确认其可行。孙剑说阅读这样的评审结果是一件让人快乐的事情。而汤老师认为,这篇文章的成功在于三个方面。第一,方法非常简单;第二,对于一个很困难的问题,给出了很好的结果;第三,发现了一个基本的自然规律并且应用在实际的问题中。在迈阿密的演讲结束后,观众也给予了很高的评价。他们跟我说,这是这次CVPR上最有趣的一个演讲。
    一位与会的研究员说,最好的idea,往往就是那些看起来很简单,但说出来大家都会觉得怎么没有人想到过的idea。而我们的idea正好就符合了这一点。我们论文摘要的第一句话是这么说的,“我们提出了一个简单而有效的方法”。或许,这就是对我们这次工作最好的概括——简单的,就是美的。

    下面正式展开基于暗通道去雾算法的原理部分:

    从上面的式(1)可以看出,A, I, J 三个向量是共面的,并且末端点共线,并可以推出t的表达式(3):


    暗通道先验:

    关于暗通道先验的概念,在作者的感言部分有很清楚的解释,就不再说了。其求法如下:

                        (4)

    Ω(x)表示以x为中心的局部块,c为三个通道,Jc表示图像J的彩色通道图像。即先求出rgb通道中像素的最小值,然后将其存进和原始图像尺寸相同的灰度图中,再对该灰度图进行最小值滤波。

    J 是无雾霾图像,除了天空区域,式(4)的值趋近与0,另外作者给出了暗通道存在的几个原因:
    (1)阴影。比如:在城市景观中的车辆、建筑物和窗户内部的阴影;在自然风景中的树叶、树木和岩石的阴影。
    (2)彩色物体或表面。比如:低反射率彩色通道的物体(如绿草、树,红花或叶子,蓝色的水面)容易在暗通道产生亮度低的值。
    (3)暗的对象或表面。比如:暗的树和石头。

    由于户外的图像通常都有阴影并且色彩斑斓,所以这些图像的暗通道确实很暗。

    估计传输率 t:

    假设大气光 A 已知,式(1)可以推出:

        ,                                (5)

    进一步假设局部块Ω(x)是连续的,块的传输率为 t~,式(5)对每个通道去最小值运算,

    t~(x)

      ,   (6)

    由于式(4)趋近于0,A>0,所以有:

    ,                                                  (7)

    把式(7)带入(6),得到:

    ,                                      (8)


    由于天空在无限远处且传输率趋于0,所以式(8)此时也是可以满足的。同时考虑到大气透视(aerial perspective)现象的存在,如果完全移除雾霾,图像看起来会不自然并且可能丢失深度信息,所以需要对远处的对象保留少量的雾霾,于是在式(8)中引入了恒参w,一般取值为0.95,式(8)变为:

                                       (9)

    基于暗通道先验去雾算法的缺陷:

    暗原色先验是一种统计的结果,是对大量户外无雾照片的统计结果,如果目标场景内存在和大气光类似,比如雪地、白色背景墙、大海等,那么由于前提条件就不成立,此时将无法获得满意的效果,但是对于一般的风景照片该算法处理效果会不错。

    基于暗通道先验去雾算法实现:

    式(4)求暗通道先验时,进行两部取最小值操作,该算法的快速实现可以参见【3】 ,这个算法的时间复杂度是O(1)。后来,作者有针对算法效率进行改进,使用了引导滤波算法,

    具体代码请参见参考附录里的博客。我这里只附录下Richal_zhang的matlab实现,她的是按照论文里最基本的实现:

    %=========================================================%  
    %调用规则:(有雾时调用,否则不调用)  
    %实际操作时,according to experiments:  
    %percent=under_50/total  
    %percent<0.1%,取w=0.6  
    %percent>0.1%&&percent<1%,取w=0.45  
    %percenet>1%&&percent<2%,取w=0.3  
    %else not use haze-free-adjust  
    %有雾:绘制出的直方图<50的部分<1%  
    %最后控制台还会输出原图中under_50像素点所占比例  
    %=========================================================%  
    close all  
    clear all  
    clc  
    blockSize=15;               %每个block为15个像素  
    w0=0.6;                      
    t0=0.1;  
    % A=200;  
    I=imread('D:\fcq_proMatlab\test_image\fog\5.jpg');  
    h = figure;  
    %set(gcf,'outerposition',get(0,'screensize'));%获得SystemScreenSize 传递给当前图像句柄gcf的outerposition属性  
    subplot(321)%表示3(行数)*2(列数)的图像,1代表所画图形的序号  
    imshow(I);  
    title('Original Image');  
    
    
    subplot(323);  
    grayI=rgb2gray(I);  
    imshow(grayI,[]);  
    title('原图像灰度图')  
    
    subplot(324);  
    imhist(grayI,64);  
    
    %统计<50的像素所占的比例  
    %%%%%%%%%%%%%%%%%%%%%%  
    [COUNT x]=imhist(grayI);  
    under_50=0;  
    for i=0:50  
        under_50=under_50+COUNT(x==i);  
    end  
    under_50  
    total=size(I,1)*size(I,2)*size(I,3);  
    percent=under_50/total;
    %%%%%%%%%%%%%%%%%%%%%%  
    
    if(percent>0.02)  
        error('This image need not Haze-Free-Proprocessing.');  
    else if(percent<0.001)  
            w=0.6;  
        else if (percent>0.01)  
                w=0.3;  
            else  
                w=0.45;  
            end  
        end  
    end  
    
    [h,w,s]=size(I);  
    min_I=zeros(h,w);  
    
    for i=1:h                   
        for j=1:w  
            dark_I(i,j)=min(I(i,j,:));%取每个点的像素为RGB分量中最低的那个通道的值  
        end  
    end  
    
    Max_dark_channel=double(max(max(dark_I)));  
    dark_channel=double(dark_I);  
    t=1-w0*(dark_channel/Max_dark_channel);  
    
    T=uint8(t*255);  
    
    t=max(t,t0);  
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    I1=double(I);  
    J(:,:,1) = uint8((I1(:,:,1) - (1-t)*Max_dark_channel)./t);  
    
    J(:,:,2) = uint8((I1(:,:,2) - (1-t)*Max_dark_channel)./t);  
    
    J(:,:,3) =uint8((I1(:,:,3) - (1-t)*Max_dark_channel)./t);  
    subplot(322)  
    imshow(J);  
    title('Haze-Free Image:');  
    
    subplot(325);  
    grayJ=rgb2gray(J);  
    imshow(grayJ,[]);  
    title('去雾后灰度图')  
    
    subplot(326);  
    imhist(grayJ,64);  

    参考:

    • 《图像去雾的最新研究进展》[J]. 自动化学报
    • 《Visibility in Bad Weather from a Single Image》[J].IEEE
    • 《STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN THREE COMPARISONS PER ELEMENT》
    • http://qianjiye.de/2015/09/haze-removal-kaiming#outline_0
    • http://blog.csdn.net/u011691310/article/details/16827695
    • http://blog.csdn.net/baimafujinji/article/details/27206237
    • http://blog.csdn.net/baimafujinji/article/details/30060161
    • http://www.cnblogs.com/Imageshop/p/3281703.html
    • http://blog.csdn.net/s12244315/article/details/50292049
    展开全文
  • 基于暗通道先验算法的图像去雾技术已经日益成熟,但是其处理速度慢、天空区域过曝、处理完的图像色彩变暗等缺点也很明显。本文针对这几方面分别提出了求透射率时的优化、纠正天空等明亮区域的错误估计的透射率、采用...
  • 针对暗通道先验去雾算法受限于天空区域,提出了一种天空优化的暗通道去雾算法。在原有公式的基础上,根据天空大气光中值及天空区域的占比添加天空区域约束因子,使其适用于天空区域的处理。同时,在原含雾图像对透射...
  • 基于暗通道先验的去雾算法:2009CVPR+导向滤波+容差控制 基于CLAHE限制自适应直方图去雾算法
  • Opencv实现暗通道先验去雾算法

    千次阅读 2017-10-03 10:58:00
    Opencv实现暗通道先验去雾算法今天读了何凯明博士的《Single Image Haze Removal Using Dark Channel Prior》,用opencv实现了一遍。 其中,暗通道及T(x)采用腐蚀erode算法。 具体代码如下:class ...

    今天读了何凯明博士的《Single Image Haze Removal Using Dark Channel Prior》,用opencv实现了一遍。
    其中,暗通道及T(x)采用腐蚀erode算法。
    具体代码如下:

    class hazeOutProcesser
    {
        public:
    
        hazeOutProcesser(const Mat& srcImg)
        {
            this->srcImg = srcImg.clone();
        }
    
        hazeOutProcesser(){}
    
        void getDarkChannel(const Mat& src, Mat& dst)
        {
            Mat element = getStructuringElement(MORPH_RECT, Size(15, 15));
            // erode(src, dst, element);
            morphologyEx(src, dst, MORPH_ERODE, element);
            // Mat temp = dst.clone();
            // bilateralFilter(temp, dst, 25, 25*2, 25/2);
            cvtColor(dst, dst, COLOR_BGR2GRAY);
        }
    
        std::pair<Mat, Mat> getAtmosphere(const Mat& sortedIndex,
                             const Mat& srcImg,
                             double rate = 0.001)
        {
            long n_pixel = srcImg.rows * srcImg.cols;
            long n_search = n_pixel * rate;
            Mat srcImgVector = srcImg.reshape(1, n_pixel).clone(); //1channel, 3rows(b, g, r)
                                                           //(n_pixel, 3)
            cout << srcImgVector.rows << ", " << srcImgVector.cols << endl;
            Mat accumulator = Mat::zeros(1, 3, CV_32FC1); //[(32f)0, (32f)0, (32f)0]
            for(int i = 0; i < n_search; ++i) {
                int index = sortedIndex.at<int>(i);
                Mat temp = srcImgVector.rowRange(index, index + 1).clone();   //row:index
                accumulator += temp;
            }
            Mat atmosphere = accumulator / n_search;
            for(int i = 0; i < n_pixel; ++i) {
                srcImgVector.rowRange(i, i + 1) = srcImgVector.rowRange(i, i + 1).mul(1 / atmosphere);
            }
            srcImgVector = srcImgVector.reshape(3).reshape(3, srcImg.rows);
            // imshow("srcImgVector", srcImgVector);
            // cout << srcImgVector.rows << ", " << srcImgVector.cols << endl;
            return std::make_pair(atmosphere, srcImgVector);
        }
    
        Mat getTransmission(const Mat& srcImg_divide_atmosphere) 
        {
            float omega = 0.95;
            Mat darkchannel_SDA;
            getDarkChannel(srcImg_divide_atmosphere, darkchannel_SDA);
            darkchannel_SDA *= omega; darkchannel_SDA = 1 - darkchannel_SDA;
            // imshow("sda", darkchannel_SDA);
            return darkchannel_SDA;
        }
    
        Mat hazeOut(const Mat& srcImg, const Mat& T, const Mat& A)
        {
            const float T0 = 0.1;
            // long n_pixel = srcImg.rows * srcImg.cols;
            // Mat srcImgVector = srcImg.reshape(1, n_pixel);
    
            // for(int i = 0; i < n_pixel; ++i) {
            //  srcImgVector.rowRange(i, i + 1) -= A;
            // }
            /*-------------------need improve-------------------
                I use a loop to calculate the result matrix
                but using maxtrix directly might speed up the procession
            */
    
            Mat hazeout(srcImg.size(), CV_32FC3);
            for(size_t r = 0; r < hazeout.rows; ++r) {
                for(size_t c = 0; c < hazeout.cols; ++c) {
                    hazeout.at<Vec3f>(r, c) = 
                        Vec3f((srcImg.at<Vec3f>(r, c)[0] - A.at<float_t>(0, 0)) / max(T.at<float>(r, c), T0) + A.at<float_t>(0, 0),
                              (srcImg.at<Vec3f>(r, c)[1] - A.at<float_t>(0, 1)) / max(T.at<float>(r, c), T0) + A.at<float_t>(0, 1),
                              (srcImg.at<Vec3f>(r, c)[2] - A.at<float_t>(0, 2)) / max(T.at<float>(r, c), T0) + A.at<float_t>(0, 2));
                }
            }
    
    
            return hazeout;
        }
    
        void process()
        {
            getDarkChannel(srcImg, darkchannel);
    
    
            Mat matVector = darkchannel.reshape(1, 1);                       //(1, 116749)
            // cout << matVector.rows << ", " << matVector.cols << endl;
            Mat_<int> sortIndex;
            sortIdx(darkchannel, 
                    sortIndex, 
                    CV_SORT_EVERY_ROW | SORT_DESCENDING);
            std::pair<Mat, Mat> p = getAtmosphere(sortIndex, srcImg);
            /*-------check minimum-------------
                double min = 0;
                minMaxLoc(transmission, &min);
                cout << "min:" << min << endl;
            */
            Mat transmission = getTransmission(p.second);
            hazeOut(srcImg, transmission, p.first).convertTo(dstImg, CV_8UC3, 255, 0);
        }
    
        void setInput(const Mat& input)
        {
            srcImg = input.clone();
        }
    
        Mat getDarkChannelOutput()
        {
            return darkchannel;
        }
    
        Mat getOutput()
        {
            return dstImg;
        }
    
        private:
    
        Mat srcImg;
        Mat dstImg;
        Mat darkchannel;
    };
    #include <opencv2/opencv.hpp>
    using namespace cv;
    using namespace std;
    int main(int argc, char const *argv[])
    {
        Mat srcImg = imread("../../../data/sea1.jpg");
        assert(srcImg.data && srcImg.channels() == 3);
        srcImg.convertTo(srcImg, CV_32FC3, 1.0 / 255, 0);
        hazeOutProcesser processer;
        processer.setInput(srcImg);
        processer.process();
        imshow("darkchannel1", processer.getDarkChannelOutput());
        imshow("src1", srcImg);
        imshow("hazeout1", processer.getOutput());
    
    
        srcImg = imread("../../../data/sea2.png");
        assert(srcImg.data && srcImg.channels() == 3);
        srcImg.convertTo(srcImg, CV_32FC3, 1.0 / 255, 0);
        processer.setInput(srcImg);
        processer.process();
        // imshow("darkchannel2", processer.getDarkChannelOutput());
        imshow("src2", srcImg);
        imshow("hazeout2", processer.getOutput());
    
        srcImg = imread("../../../data/sea4.jpg");
        assert(srcImg.data && srcImg.channels() == 3);
        resize(srcImg, srcImg, Size(), 0.5, 0.5);
        srcImg.convertTo(srcImg, CV_32FC3, 1.0 / 255, 0);
        processer.setInput(srcImg);
        processer.process();
        // imshow("darkchannel2", processer.getDarkChannelOutput());
        imshow("src3", srcImg);
        imshow("hazeout3", processer.getOutput());
        waitKey();
        return 0;
    }

    后记

    为了提高运算速度,尽量使用矩阵运算,当然以上还有需要改善的地方。
    当然,本人水平有限,若有地方存在问题,还欢迎留言一起探讨!

    展开全文
  • 何恺明的暗通道先验(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)

    展开全文
  • 针对暗通道先验在天空区域的失效问题,提出了一种基于亮度模型融合的改进暗通道先验图像去雾算法。首先通过Canny算子分割得到天空区域与非天空区域;其次,利用亮度模拟景深,重构亮度透射率,并通过与暗通道透射率的融合...
  • 传统基于暗通道先验的图像去雾算法不能有效去除有雾图像在景深突变处的雾点,边界处容易引起光晕效应,对此提出一种基于暗通道先验的自适应超像素去雾算法.首先,在暗通道的获取过程中引入自适应方法判断当前像素邻域内...
  • 基于导向滤波的暗通道先验去雾算法(Python语言,可直接运行)1 编译环境2 原理介绍2.1 暗通道先验2.1.1 暗通道先验理论与去雾模型2.1.2 处理步骤2.2 导向滤波求t(x)3 代码4 处理结果4.1 第一组结果及其比较,暗通道...
  • 基于暗通道先验条件图像去雾算法

    千次阅读 2019-08-30 17:41:37
    基于暗通道先验条件图像去雾算法 香港大学何凯明博士于2009发表了一篇论文《Single Image Haze Removal Using Dark Channel Prior 》。在文章中,何凯明博士提出了一种简单而有效的图像先验暗通道消除单输入图像雾...
  • 基于何凯文博士的参考文献He K, Sun J, Tang X. Single image haze removal using dark channel prior[J]. IEEE CVPR, 2009.所写的基于暗通道先验的去雾算法matlab代码,内有代码,论文,及测试图片
  • 基于暗通道先验条件图像去雾算法香港大学何凯明博士于2009发表了一篇论文《Single Image Haze Removal Using Dark Channel Prior 》。在文章中,何凯明博士提出了一种简单而有效的图像先验暗通道消除单输入图像雾的...
  • 暗通道先验去雾算法,可以较好的去除大量雾。
  • 何恺明的暗通道先验(dark channel prior)去雾算法是CV界去雾领域很有名的算法,关于该算法的论文"Single Image Haze Removal Using Dark Channel Prior"一举获得2009年CVPR最佳论文。作者统计了大量的无雾图像,...
  • 北航提出FFA-Net去雾新网络和何凯明博士提出的暗通道去雾算法,现所有源码已开源。其论文链接:https://arxiv.org/abs/1911.07559。而今天我们就将针对这两个项目进行实践。其中得到的去雾效果如下:作者 | 李秋键...
  • 数字图像处理、去雾算法
  • 在写作和发表论文前,这篇博客不会介绍任何实质性的内容,仅为展示对当前基于大气光传播模型和暗通道先验的主流去雾算法中的几处误区或是不合理的地方进行修正后(并加了简单的预处理)获得的去雾效果图(算法比何...
  • 针对以上情况,结合大气散射模型和夜间雾天图像成像特点,提出基于Retinex理论和暗通道先验的去雾算法。首先,根据Retinex理论求得夜间场景的有雾入射光图像和有雾反射光图像;其次,利用暗通道先验得到场景的无雾...
  • 研究暗通道去雾是去年的事情了,当时正值...首先还是简单说下暗通道先验去雾算法,这里只是记录性质,详细推荐看论文,写的非常好: 雾图模型 I(x) ——待去雾的图像 J(x)——无雾图像 A——全球大气光成分 t——折
  • 针对目前检测方法适应性不强、在复杂环境下检测性能不高的问题,提出了一种基于背景动态更新和暗通道先验的烟雾检测算法算法首先通过改进的背景动态更新算法提取运动前景;然后,结合暗通道先验知识确定前景中的...
  • 暗通道先验的去雾算法实现

    千次阅读 2015-11-01 18:37:00
    暗通道优先的图像去雾算法 http://blog.csdn.net/baimafujinji/article/details/30060161 高级图像去雾算法的快速实现 http://blog.csdn.net/huixingshao/article/details/42834939 基于暗通道优先算法...
  • 一、雾图形成模型: —待去雾图像;—恢复的无雾的图像;—全球大气光成分;—透射率 变形: ...c为R、G、B三通道。...经过大量实验,局部找最点进行均匀去雾有很好的效果,由此推出Dark Channel Prior将全局改为
  • 针对暗原色先验去雾算法中雾霾图像明亮区域透射率估计过小,造成图像色彩失真的问题,提出一种新的基于...实验结果表明,比值重估透射率去雾算法暗通道和容差机制去雾算法相比,该算法恢复的图像更接近于真实图像。
  • 经典去雾算法--暗通道先验去雾(DCP)

    千次阅读 2020-04-04 14:03:54
    基于暗通道先验的单幅图像去雾算法来自于何凯明博士2009年的CVPR论文:《Single Image Haze Removal Using Dark Channel Prior》,2009年的CVPR只有一篇文章被选为那年的最佳论文。这是CVPR创立25年以来首次由中国人...
  •  首先看看暗通道先验是什么:  在绝大多数非天空的局部区域里,某一些像素总会有至少一个颜色通道具有很低的值。换言之,该区域光强度的最小值是个很小的数。  我们给暗通道一个数学定义,对于任意的...

空空如也

空空如也

1 2 3 4 5 ... 8
收藏数 155
精华内容 62
关键字:

暗通道先验算法