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

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

算法原理：

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

       ① 暗通道定义

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

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

 

 

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

   

 J^c代表J的某一个颜色通道，而Ω(x)是以x为中心的一块方形区域。我们观察得出，除了天空方位，J^dark的强度总是很低并且趋近于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的公式如下：

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

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

 

④求解全局大气光照A^c

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

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

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

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

4.A^c=[A^r,A^g,A^b]. 其中A^r=sum_r/N;   A^g=sum_g/N;   A^b=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，那么带入大气模型公式：

这里，t₀参数用来限定透射率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文件。

 

/*--------------------------------------------------------------------------------*\

   This program is free software; permission is hereby granted to use, copy, modify,
   and distribute this source code, or portions thereof, for any purpose, without fee,
   subject to the restriction that the copyright notice may not be removed
   or altered from any source or altered source distribution.
   The software is released on an as-is basis and without any warranties of any kind.
   In particular, the software is not guaranteed to be fault-tolerant or free from
   failure.

   The author disclaims all warranties with regard to this software, any use,
   and any consequent failure, is purely the responsibility of the user.

   Copyright (C) 2013-2016 Changkai Zhao, www.cnblogs.com/changkaizhao
\*-------------------------------------------------------------------------------*/
#include "dehazor.h"

cv::Mat Dehazor::process(const cv::Mat &image)
{
    int dimr=image.rows;
    int dimc=image.cols;

    cv::Mat rawtemp;
    cv::Mat refinedImage_temp;
    cv::Mat output_b_temp(image.rows,image.cols,CV_32F);
    cv::Mat output_g_temp(image.rows,image.cols,CV_32F);
    cv::Mat output_r_temp(image.rows,image.cols,CV_32F);
    cv::Mat output_temp(dimr,dimc,image.type());
    rawImage.create(image.rows,image.cols,CV_8U); 
    refinedImage_temp.create(image.rows,image.cols,CV_32F);
    rawtemp.create(image.rows,image.cols,CV_8U);


    float sumb=0;float sumg=0;float sumr=0;
    float Air_b;  float Air_g; float Air_r;
    cv::Point2i pt;

    int dimone=floor(dimr*dimc*0.001);
    cv::Mat quantile;
    quantile=cv::Mat::zeros(2,dimone,CV_8U);

    int dx=floor(windowsize/2);
//    cv::Mat imagetemp;
//    imagetemp.create(image.rows,image.cols,CV_32FC3);

    for(int j=0;j<dimr;j++){
        for(int i=0;i<dimc;i++){
            int min;
            image.at<cv::Vec3b>(j,i)[0] >= image.at<cv::Vec3b>(j,i)[1]?
               min=image.at<cv::Vec3b>(j,i)[1]:
               min=image.at<cv::Vec3b>(j,i)[0];
            if(min>=image.at<cv::Vec3b>(j,i)[2])
                min=image.at<cv::Vec3b>(j,i)[2];
            rawImage.at<uchar>(j,i)=min;
        }
    }

    for(int j=0;j<dimr;j++){
        for(int i=0;i<dimc;i++){
            int min=255;

            int jlow=j-dx;int jhigh=j+dx;
            int ilow=i-dx;int ihigh=i+dx;

            if(ilow<=0)
                ilow=0;
            if(ihigh>=dimc)
                ihigh=dimc-1;
            if(jlow<=0)
                jlow=0;
            if(jhigh>=dimr)
                jhigh=dimr-1;

            for(int m=jlow;m<=jhigh;m++){
                for(int n=ilow;n<=ihigh;n++){
                    if(min>=rawImage.at<uchar>(m,n))
                        min=rawImage.at<uchar>(m,n);
                }
            }
           rawtemp.at<uchar>(j,i)=min;
        }
    }

    for(int i=0;i<dimone;i++){
        cv::minMaxLoc(rawtemp,0,0,0,&pt);
        sumb+= image.at<cv::Vec3b>(pt.y,pt.x)[0];
        sumg+= image.at<cv::Vec3b>(pt.y,pt.x)[1];
        sumr+= image.at<cv::Vec3b>(pt.y,pt.x)[2];
        rawtemp.at<uchar>(pt.y,pt.x)=0;
    }
    Air_b=sumb/dimone;
    Air_g=sumg/dimone;
    Air_r=sumr/dimone;


    cv::Mat layb;  cv::Mat Im_b;
    cv::Mat layg;  cv::Mat Im_g;
    cv::Mat layr;  cv::Mat Im_r;

    // create vector of 3 images
    std::vector<cv::Mat> planes;
    // split 1 3-channel image into 3 1-channel images
    cv::split(image,planes);

    layb=planes[0];
    layg=planes[1];
    layr=planes[2];
    Im_b=planes[0];
    Im_g=planes[1];
    Im_r=planes[2];

    layb.convertTo(layb, CV_32F);
    layg.convertTo(layg, CV_32F);
    layr.convertTo(layr, CV_32F);

    Im_b.convertTo(Im_b, CV_32F);
    Im_g.convertTo(Im_g, CV_32F);
    Im_r.convertTo(Im_r, CV_32F);


    for (int j=0; j<dimr; j++) {
        for (int i=0; i<dimc; i++) {
        // process each pixel ---------------------
        layb.at<float>(j,i)=layb.at<float>(j,i)/Air_b;

        layg.at<float>(j,i)=layg.at<float>(j,i)/Air_g;

        layr.at<float>(j,i)=layr.at<float>(j,i)/Air_r;
        // end of pixel processing ----------------
        } // end of line
    }


    rawtemp.convertTo(rawtemp,CV_32F);


    for(int j=0;j<dimr;j++){
        for(int i=0;i<dimc;i++){
            float min;
            layb.at<float>(j,i) >= layg.at<float>(j,i)?
               min=layg.at<float>(j,i):
               min=layb.at<float>(j,i);
            if(min>=layr.at<float>(j,i))
                min=layr.at<float>(j,i);
            rawtemp.at<float>(j,i)=min;
        }
    }
    for(int j=0;j<dimr;j++){
        for(int i=0;i<dimc;i++){
            float min=1;

            int jlow=j-dx;int jhigh=j+dx;
            int ilow=i-dx;int ihigh=i+dx;

            if(ilow<=0)
                ilow=0;
            if(ihigh>=dimc)
                ihigh=dimc-1;
            if(jlow<=0)
                jlow=0;
            if(jhigh>=dimr)
                jhigh=dimr-1;

            for(int m=jlow;m<=jhigh;m++){
                for(int n=ilow;n<=ihigh;n++){
                    if(min>=rawtemp.at<float>(m,n))
                        min=rawtemp.at<float>(m,n);
                }
            }
           rawImage.at<uchar>(j,i)=(1-(float)fog_reservation_factor*min)*255;
        }
    }

   refinedImage_temp=guildedfilter_color(image,rawImage,localwindowsize,eps);

   for(int j=0;j<dimr;j++){
       for(int i=0;i<dimc;i++){
           if(refinedImage_temp.at<float>(j,i)<0.1)
               refinedImage_temp.at<float>(j,i)=0.1;
       }
   }

   cv::Mat onemat(dimr,dimc,CV_32F,cv::Scalar(1));

   cv::Mat air_bmat(dimr,dimc,CV_32F);
   cv::Mat air_gmat(dimr,dimc,CV_32F);
   cv::Mat air_rmat(dimr,dimc,CV_32F);

   cv::addWeighted(onemat,Air_b,onemat,0,0,air_bmat);
   cv::addWeighted(onemat,Air_g,onemat,0,0,air_gmat);
   cv::addWeighted(onemat,Air_r,onemat,0,0,air_rmat);


   output_b_temp=Im_b-air_bmat;
   output_g_temp=Im_g-air_gmat;
   output_r_temp=Im_r-air_rmat;

   output_b_temp=output_b_temp.mul(1/refinedImage_temp,1)+air_bmat;
   output_g_temp=output_g_temp.mul(1/refinedImage_temp,1)+air_gmat;
   output_r_temp=output_r_temp.mul(1/refinedImage_temp,1)+air_rmat;


   output_b_temp.convertTo(output_b_temp,CV_8U);
   output_g_temp.convertTo(output_g_temp,CV_8U);
   output_r_temp.convertTo(output_r_temp,CV_8U);

   for(int j=0;j<dimr;j++){
       for(int i=0;i<dimc;i++){
           output_temp.at<cv::Vec3b>(j,i)[0]=output_b_temp.at<uchar>(j,i);
           output_temp.at<cv::Vec3b>(j,i)[1]=output_g_temp.at<uchar>(j,i);
           output_temp.at<cv::Vec3b>(j,i)[2]=output_r_temp.at<uchar>(j,i);
       }
   }

   cv::Mat nom_255(dimr,dimc,CV_32F,cv::Scalar(255));

   cv::Mat ref_temp(dimr,dimc,CV_32F);
   ref_temp=refinedImage_temp;
   ref_temp=ref_temp.mul(nom_255,1);
   ref_temp.convertTo(refinedImage,CV_8U);


    return output_temp;
}
cv::Mat Dehazor::boxfilter(cv::Mat &im, int r)
{
    //im is a CV_32F type mat [0,1] (normalized)
    //output is the same size to im;

    int hei=im.rows;
    int wid=im.cols;
    cv::Mat imDst;
    cv::Mat imCum;


    imDst=cv::Mat::zeros(hei,wid,CV_32F);
    imCum.create(hei,wid,CV_32F);

    //cumulative sum over Y axis
    for(int i=0;i<wid;i++){
        for(int j=0;j<hei;j++){
            if(j==0)
                imCum.at<float>(j,i)=im.at<float>(j,i);
            else
                imCum.at<float>(j,i)=im.at<float>(j,i)+imCum.at<float>(j-1,i);
        }
    }


    //difference over Y axis
    for(int j=0;j<=r;j++){
        for(int i=0;i<wid;i++){
            imDst.at<float>(j,i)=imCum.at<float>(j+r,i);
        }
    }
    for(int j=r+1;j<=hei-r-1;j++){
        for(int i=0;i<wid;i++){
            imDst.at<float>(j,i)=imCum.at<float>(j+r,i)-imCum.at<float>(j-r-1,i);
        }
    }
    for(int j=hei-r;j<hei;j++){
        for(int i=0;i<wid;i++){
            imDst.at<float>(j,i)=imCum.at<float>(hei-1,i)-imCum.at<float>(j-r-1,i);
        }
    }


    //cumulative sum over X axis
      for(int j=0;j<hei;j++){
          for(int i=0;i<wid;i++){
              if(i==0)
                  imCum.at<float>(j,i)=imDst.at<float>(j,i);
              else
                  imCum.at<float>(j,i)=imDst.at<float>(j,i)+imCum.at<float>(j,i-1);
          }
      }
      //difference over X axis
      for(int j=0;j<hei;j++){
          for(int i=0;i<=r;i++){
              imDst.at<float>(j,i)=imCum.at<float>(j,i+r);
          }
      }
      for(int j=0;j<hei;j++){
          for(int i=r+1;i<=wid-r-1;i++){
              imDst.at<float>(j,i)=imCum.at<float>(j,i+r)-imCum.at<float>(j,i-r-1);
          }
      }
      for(int j=0;j<hei;j++){
          for(int i=wid-r;i<wid;i++){
              imDst.at<float>(j,i)=imCum.at<float>(j,wid-1)-imCum.at<float>(j,i-r-1);
          }
      }

    return imDst;
}
cv::Mat Dehazor::guildedfilter_color(const cv::Mat &Img, cv::Mat &p, int r, float &epsi)
{

    int hei=p.rows;
    int wid=p.cols;

    cv::Mat matOne(hei,wid,CV_32F,cv::Scalar(1));
    cv::Mat N;

    N=boxfilter(matOne,r);



    cv::Mat mean_I_b(hei,wid,CV_32F);
    cv::Mat mean_I_g(hei,wid,CV_32F);
    cv::Mat mean_I_r(hei,wid,CV_32F);
    cv::Mat mean_p(hei,wid,CV_32F);


    cv::Mat Ip_b(hei,wid,CV_32F);
    cv::Mat Ip_g(hei,wid,CV_32F);
    cv::Mat Ip_r(hei,wid,CV_32F);
    cv::Mat mean_Ip_b(hei,wid,CV_32F);
    cv::Mat mean_Ip_g(hei,wid,CV_32F);
    cv::Mat mean_Ip_r(hei,wid,CV_32F);
    cv::Mat cov_Ip_b(hei,wid,CV_32F);
    cv::Mat cov_Ip_g(hei,wid,CV_32F);
    cv::Mat cov_Ip_r(hei,wid,CV_32F);

    cv::Mat II_bb(hei,wid,CV_32F);
    cv::Mat II_gg(hei,wid,CV_32F);
    cv::Mat II_rr(hei,wid,CV_32F);
    cv::Mat II_bg(hei,wid,CV_32F);
    cv::Mat II_br(hei,wid,CV_32F);
    cv::Mat II_gr(hei,wid,CV_32F);

    cv::Mat var_I_bb(hei,wid,CV_32F);
    cv::Mat var_I_gg(hei,wid,CV_32F);
    cv::Mat var_I_rr(hei,wid,CV_32F);
    cv::Mat var_I_bg(hei,wid,CV_32F);
    cv::Mat var_I_br(hei,wid,CV_32F);
    cv::Mat var_I_gr(hei,wid,CV_32F);

    cv::Mat layb;
    cv::Mat layg;
    cv::Mat layr;
    cv::Mat P_32;

    // create vector of 3 images
    std::vector<cv::Mat> planes;
    // split 1 3-channel image into 3 1-channel images
    cv::split(Img,planes);

    layb=planes[0];
    layg=planes[1];
    layr=planes[2];

    layb.convertTo(layb, CV_32F);
    layg.convertTo(layg, CV_32F);
    layr.convertTo(layr, CV_32F);

    p.convertTo(P_32,CV_32F);
    cv::Mat nom_255(hei,wid,CV_32F,cv::Scalar(255));



    layb=layb.mul(1/nom_255,1);
    layg=layg.mul(1/nom_255,1);
    layr=layr.mul(1/nom_255,1);
    P_32=P_32.mul(1/nom_255,1);

    cv::Mat mean_I_b_temp=boxfilter(layb,r);
    cv::Mat mean_I_g_temp=boxfilter(layg,r);
    cv::Mat mean_I_r_temp=boxfilter(layr,r);
    cv::Mat mean_p_temp=boxfilter(P_32,r);

    mean_I_b=mean_I_b_temp.mul(1/N,1);
    mean_I_g=mean_I_g_temp.mul(1/N,1);
    mean_I_r=mean_I_r_temp.mul(1/N,1);
    mean_p=mean_p_temp.mul(1/N,1);

    Ip_b=layb.mul(P_32,1);
    Ip_g=layg.mul(P_32,1);
    Ip_r=layr.mul(P_32,1);

    cv::Mat mean_Ip_b_temp=boxfilter(Ip_b,r);
    cv::Mat mean_Ip_g_temp=boxfilter(Ip_g,r);
    cv::Mat mean_Ip_r_temp=boxfilter(Ip_r,r);

    mean_Ip_b=mean_Ip_b_temp.mul(1/N,1);
    mean_Ip_g=mean_Ip_g_temp.mul(1/N,1);
    mean_Ip_r=mean_Ip_r_temp.mul(1/N,1);

    cov_Ip_b=mean_Ip_b-mean_I_b.mul(mean_p,1);
    cov_Ip_g=mean_Ip_g-mean_I_g.mul(mean_p,1);
    cov_Ip_r=mean_Ip_r-mean_I_r.mul(mean_p,1);


//     variance of I in each local patch: the matrix Sigma in Eqn (14).
//     Note the variance in each local patch is a 3x3 symmetric matrix:
//                   bb, bg, br
//           Sigma = bg, gg, gr
//                   br, gr, rr
    II_bb=layb.mul(layb,1);
    II_gg=layg.mul(layg,1);
    II_rr=layr.mul(layr,1);
    II_bg=layb.mul(layg,1);
    II_br=layb.mul(layr,1);
    II_gr=layg.mul(layr,1);

    cv::Mat bb_box=boxfilter(II_bb,r);
    cv::Mat gg_box=boxfilter(II_gg,r);
    cv::Mat rr_box=boxfilter(II_rr,r);
    cv::Mat bg_box=boxfilter(II_bg,r);
    cv::Mat br_box=boxfilter(II_br,r);
    cv::Mat gr_box=boxfilter(II_gr,r);

    var_I_bb=bb_box.mul(1/N,1)-mean_I_b.mul(mean_I_b);
    var_I_gg=gg_box.mul(1/N,1)-mean_I_g.mul(mean_I_g);
    var_I_rr=rr_box.mul(1/N,1)-mean_I_r.mul(mean_I_r);
    var_I_bg=bg_box.mul(1/N,1)-mean_I_b.mul(mean_I_g);
    var_I_br=br_box.mul(1/N,1)-mean_I_b.mul(mean_I_r);
    var_I_gr=gr_box.mul(1/N,1)-mean_I_g.mul(mean_I_r);

    cv::Mat a_b(hei,wid,CV_32F);
    cv::Mat a_g(hei,wid,CV_32F);
    cv::Mat a_r(hei,wid,CV_32F);

    cv::Mat b(hei,wid,CV_32F);
    cv::Mat sigma(3,3,CV_32F,cv::Scalar(0));
    cv::Mat inv_sigma(3,3,CV_32F);

    for(int j=0;j<hei;j++){
        for(int i=0;i<wid;i++){
            sigma.at<float>(0,0)=var_I_rr.at<float>(j,i)+epsi;
            sigma.at<float>(0,1)=var_I_gr.at<float>(j,i);
            sigma.at<float>(0,2)=var_I_br.at<float>(j,i);
            sigma.at<float>(1,0)=var_I_gr.at<float>(j,i);
            sigma.at<float>(2,0)=var_I_br.at<float>(j,i);
            sigma.at<float>(1,1)=var_I_gg.at<float>(j,i)+epsi;
            sigma.at<float>(2,2)=var_I_bb.at<float>(j,i)+epsi;
            sigma.at<float>(1,2)=var_I_bg.at<float>(j,i);
            sigma.at<float>(2,1)=var_I_bg.at<float>(j,i);
            inv_sigma=sigma.inv(cv::DECOMP_LU);

            a_r.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,0)+
                               cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,0)+
                               cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,0);
            a_g.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,1)+
                               cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,1)+
                               cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,1);
            a_b.at<float>(j,i)=cov_Ip_r.at<float>(j,i)*inv_sigma.at<float>(0,2)+
                               cov_Ip_g.at<float>(j,i)*inv_sigma.at<float>(1,2)+
                               cov_Ip_b.at<float>(j,i)*inv_sigma.at<float>(2,2);


        }
    }
    b=mean_p-a_b.mul(mean_I_b,1)-a_g.mul(mean_I_g,1)-a_r.mul(mean_I_r,1);

    cv::Mat box_ab=boxfilter(a_b,r);
    cv::Mat box_ag=boxfilter(a_g,r);
    cv::Mat box_ar=boxfilter(a_r,r);
    cv::Mat box_b=boxfilter(b,r);
    cv::Mat q(hei,wid,CV_32F);


    q=box_ab.mul(layb,1)+box_ag.mul(layg,1)+box_ar.mul(layr,1)+box_b;
    q=q.mul(1/N,1);


    return q;

}

 

       截图：