暗通道去霧與自適應直方圖均衡去霧之比較

本文轉載自查看原文 2013-11-05 10:14 3789 圖像處理/ 信號處理

本文主要對比兩種去霧增強算法的優略。

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 性能比較