图像增强之(三) --- 同态滤波器增强算法

光照反射模型表示一副图像f（x，y）可以由入射光强i（x，y）（照明分量）和物体的反射率r（x，y）（反射分量）表示：

照明反射模型其中,入射光强取决于光源,由于入射光光强分布函数i(x,y)随空间变化较小,在空间频率域,函数i(x,y)的频谱主要集中于低频段,反映了成像的环境条件.反射光r(x,y)的频谱主要集中于相对高频段比较宽的范围,反映了物体本身的特性.为此,只要能把图像中相对于i(x,y)和r(x,y)两部分分开,然后压制较低频段,放大较高频段,就能有效地降低光照不均匀对图像所带来的影响.该模型在灰度校正和同态滤波中起着重要的作用.

上式不能直接对照明分量和反射分量进行操作，于是对等式两边取对数：

为了能够提取高频和低频部分，再次进行傅立叶变换：

记为：

其中I(u,v)表示图像中的低频部分（光照分量），R（u，v）表示图像中的高频部分（反射分量）。
对图像进行高通滤波，可以保留R（u，v）可以有效减少光照的干扰。
设采用的高通滤波器为H(u,v).滤波后的图像为S。

代码来自：http://blog.csdn.net/lilingyu520/article/details/46654265

void my_HomoFilter(Mat srcImg, Mat &dst)  {      srcImg.convertTo(srcImg, CV_64FC1);      dst.convertTo(dst, CV_64FC1);      //1. ln      for (int i = 0; i < srcImg.rows; i++)      {          double* srcdata = srcImg.ptr<double>(i);          double* logdata = srcImg.ptr<double>(i);          for (int j = 0; j < srcImg.cols; j++)          {              logdata[j] = log(srcdata[j]+0.0001);          }      }      //spectrum      //2. dct      Mat mat_dct = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);      dct(srcImg, mat_dct);      imshow("dct", mat_dct);      //3. linear filter      Mat H_u_v;      double gammaH = 1.5;      double gammaL = 0.5;      double C = 1;      double d0 = (srcImg.rows/2)*(srcImg.rows/2) + (srcImg.cols/2)*(srcImg.cols/2);      double d2 = 0;      H_u_v = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);      double totalWeight = 0.0;      for (int i = 0; i < srcImg.rows; i++)      {          double * dataH_u_v = H_u_v.ptr<double>(i);          for (int j = 0; j < srcImg.cols; j++)          {              d2 = pow((i), 2.0) + pow((j), 2.0);              dataH_u_v[j] =  (gammaH - gammaL)*(1 - exp(-C*d2/d0)) + gammaL;              totalWeight += dataH_u_v[j];          }      }      H_u_v.ptr<double>(0)[0] = 1.1;      //H_u_v = Mat::ones(srcImg.rows, srcImg.cols, CV_64FC1);      imshow("H_u_v", H_u_v);      //imshow("before filter", mat_dct);      mat_dct = mat_dct.mul(H_u_v);      //Mat tmp = mat_dct.mul(H_u_v);      //tmp.copyTo(mat_dct);      //4. idct      idct(mat_dct, dst);  #if 0      //spatial high high pass filter      Mat tmp = Mat::zeros(srcImg.rows, srcImg.cols, CV_64FC1);      GaussianBlur(srcImg, tmp, Size(9, 9), 1.5, 1.5);      const double alpha = 0.5;      for (int i = 0; i < srcImg.rows; i++)      {          double* srcdata = srcImg.ptr<double>(i);          double* blurdata = tmp.ptr<double>(i);          double* dstdata = dst.ptr<double>(i);          for (int j = 0; j < srcImg.cols; j++)          {              dstdata[j] = (1+alpha)*srcdata[j] - alpha*blurdata[j];              //dstdata[j] = blurdata[j];          }      }  #endif      //5. exp      for (int i = 0; i < srcImg.rows; i++)      {          double* srcdata = dst.ptr<double>(i);          double* dstdata = dst.ptr<double>(i);          for (int j = 0; j < srcImg.cols; j++)          {              dstdata[j] = exp(srcdata[j]);          }      }      //imshow("dst", dst);      dst.convertTo(dst, CV_8UC1);  }