暗通道去雾与自适应直方图均衡去雾之比较

本文主要对比两种去雾增强算法的优略。

1. 暗通道去雾方法

这种方法目前研究的比较多。虽然大家认为去雾已经被研究烂了，但是从工程实践的角度看，仍然存在很多问题。却不说暗通道先验对天空本身的不合理假设，就单说时间复杂度。对于视频图像处理来说，视频编解码就已经占了大量的CPU时间，而如果其它一些图像增强的算法也要占用大量的计算时间，就会导致了很多所谓的the state of art 的算法在工程实践上，依然并不可行。举一个很简单的例子，目前普通PC单核的处理能力只能编码2-4路的720P的H264码流，如果图像增强算法的处理时间大于15ms，则意味着在普通PC上这种算法并不可行。

关于暗通道去雾的原理，下面几篇网络文章可以参考。

http://www.cnblogs.com/Imageshop/p/3281703.html

http://blog.csdn.net/celerychen2009/article/details/8839098

2. 对比度受限自适应直方图增强去雾方法

直方图增强是一种古老的技术。不过给我的印象是，对于大部分图像，简单的直方图增强的效果很不好。不过，前几天看到网络上一位牛人的文章，让我对直方图增强有了新的认识。当然，不是直接用直方图，而是对比度受限自适应直方图增强。

自适应直方图是使用多个直方图，例如对整个图像可以划分为8x8个小块，称为8x8的tiles，则总共会有64个直方图。用这64个直方图来做，tiles的边界处用线性插值的方法实现平滑的过渡。具体原理可以参考下面的博文：

http://www.cnblogs.com/Imageshop/archive/2013/04/07/3006334.html

http://www.cnblogs.com/Imageshop/p/3395968.html

3. 算法实现

3.1 暗通道去雾的快速算法。

暗通道去雾的快速算法改进主要集中在对透射率的refine上。有很多的非线性边缘保持滤波可以用。我起初使用的是非线性各向异性扩散方程，非线性扩散在KAZE特征提取中被采用，其相关的网络资源可以参考：http://blog.csdn.net/chenyusiyuan/article/details/8710462。采用加法分裂算子来离散化偏微分方程，仅需很少的几次迭代就可以收敛，如果采用跳跃格式的迭代算法，则可以不必求解三对角线性方程组，这种方法是快速稳定高效的。

采用双边滤波也是一种不错的选择，然而双边滤波的快速算法实在是快不起来。一种更快更好的方案是采用引导滤波guided filter. guided filter应该是目前最快的边缘保持滤波了，我的实现方案表明，其计算时间比我的快速高斯滤波算法还要快。

3.2 对比度受限自适应直方图的实现

这个算法实现相对比较简单，没有暗通道去雾那么多过程。

例如对于2x2的tiles，如果图像是64x64的，则每个block的大小是32x32，对每一个32x32的小块，计算直方图，然后根据相应的参数重新调整直方图再映射回去。当然，如果图像的尺寸不是tiles的整数倍的时候需要进行对称拓展。最后，在tiles的边缘处进行插值是比较耗时的操作。本文直接给出笔者的实现代码：

static clahe_interpolation( int width, int height, int stride, cv8u* src, cv8u* dst,int channel, 
						   Size2i* tileSize, int tiles_y, int tiles_x, cv8u* ptr_lut_data, int lut_step)
{
	int i,j;
	if (channel == 1){
		for (j = 0; j < height; j++){
			cv8u* pTmpSrc = src + j * stride;
			cv8u* pTmpDst = dst + j * stride;
			int tyf_i =  (j << CLAHE_Q)/tileSize->height - (1 << (CLAHE_Q -1));
			int ty1_i = tyf_i >> CLAHE_Q;
			int ty2_i = ty1_i + 1;

			int ya_i = tyf_i - (ty1_i << CLAHE_Q);
			cv8u* lutPlane1,*lutPlane2;

			ty1_i = (ty1_i < 0) ? 0 : (ty1_i);
			ty2_i = (ty2_i > tiles_y - 1) ? (tiles_y - 1) : (ty2_i);


			lutPlane1 = &ptr_lut_data[ty1_i * tiles_x * lut_step];
			lutPlane2 = &ptr_lut_data[ty2_i * tiles_x * lut_step];

			for (i = 0; i < width; i++){
				int srcVal;
				int ind1,ind2;
				int res_i;
				int txf_i = (i << CLAHE_Q)/tileSize->width - (1 << (CLAHE_Q -1));
				int tx1_i = txf_i >> CLAHE_Q;
				int tx2_i = tx1_i + 1;
				int xa_i = txf_i - (tx1_i << CLAHE_Q);
				int tf_1_i,tf_2_i,tf_3_i,tf_4_i;
				
				srcVal = *pTmpSrc++;

				tx1_i = (tx1_i < 0) ? 0 : (tx1_i);
				tx2_i = (tx2_i > tiles_x - 1) ? (tiles_x - 1) : (tx2_i);

				ind1 = tx1_i * lut_step + srcVal;
				ind2 = tx2_i * lut_step + srcVal;


				tf_1_i = (((1 << CLAHE_Q) - xa_i) * ( (1 << CLAHE_Q) - ya_i)) >> CLAHE_Q;
				tf_2_i = ((xa_i) * ((1 << CLAHE_Q) - ya_i)) >> CLAHE_Q;
				tf_3_i = (((1 << CLAHE_Q) - xa_i) * (ya_i)) >> CLAHE_Q;
				tf_4_i = ((xa_i) * (ya_i)) >> CLAHE_Q;

				res_i  = lutPlane1[ind1] * tf_1_i;
				res_i += lutPlane1[ind2] * tf_2_i;
				res_i += lutPlane2[ind1] * tf_3_i;
				res_i += lutPlane2[ind2] * tf_4_i;
				res_i = (res_i + (1 << (CLAHE_Q - 1))) >> CLAHE_Q;

				*pTmpDst++ = FCV_CLAMP0255_INT(res_i);
			}
		}
    }
}

3.3 性能比较