背景说明

《Single Image Haze Removal Using Dark Channel Prior》获得了CVPR 09的Best Paper Award。也是现有算法中去雾效果非常优秀的算法。另外，在作者的主页也可以获取到相关的paper和slide。

灵感来源

这个想法的产生来自于两个偶然的观察。

第一个观察来自一个3D游戏。这个游戏有很多带有雾的场景，但这些场景都是虚构的不实在的东西。计算机生成的3D图像会与自然图像的统计规律有很大区别，但人的视觉系统却仍然能感觉到虚拟图像中存在的雾。这让我相信，人的视觉系统一定有一种有效的机制来感知有雾的图像，而且这种机制一定与现存的去雾方法不一样。前人提出的去雾方法都把重点放在图像的对比度上，但虚拟场景和现实场景在对比度上的统计规律会很不一样。人的视觉系统仍然能够感知虚拟场景中的雾，说明除了对比度以外，人眼一定还在利用别的东西来感知雾。所以我觉得，这个问题里一定有人们未曾发现的更接近本质的东西。

第二个观察来自对前人的去雾方法的研究。之前最有效的去雾方法是Fattal在2008年的Siggraph文章《Single Image Dehazing》中提出来的，这篇文章是我们首要超越的目标。这篇文章里给出的比较结果中，我发现一种叫做Dark Object Subtraction的方法有时候会有更好的效果。这种方法利用了全图最暗的点来去除全局均匀的雾。如果雾的确是均匀的，这种方法就会更有效。其缺点在于它无法处理不均匀的雾，而这却正是去雾问题中的难点。因此自然的想法就是局部利用Dark Object Subtraction处理图像。而恰巧这样做并不需要利用对比度，说明它与之前的方法有了本质的区别。让人吃惊的是，在大量的实验中，我发现这么简单的想法，其效果却非常好。

什么是暗通道

在绝大多数非天空的局部区域里，某一些像素总会有至少一个颜色通道具有很低的值。换言之，该区域光强度的最小值是个很小的数。

暗通道可以由如下公式描述：

$${\bf J}^{dark}(x) = \min_{ y \in \Omega(x) }(\min_{ c \in { r,g,b } }({\bf J}^c (y))),(1) $$

式中${\bf J}^c$表示彩色图像的每个通道，$\Omega(x)$表示以像素X为中心的一个窗口。

什么是暗通道先验

The dark channel prior is based on the following observation on haze-free outdoor images: in most of the non-sky patches, at least one color channel has very low intensity at some pixels. In other words, the minimum intensity in such a patch should has a very low value.

这里要说明的是无雾图像的暗通道才是趋向于0的，有雾图像一般不会趋向于0。

无雾图像的暗通道先验表示如下

$$ \min_{ y \in \Omega(x) }(\min_{ c \in { r,g,b } }({\bf J}^c (y))) \rightarrow 0,(2) $$

暗通道中低亮度像素点产生的原因

shadows. e.g.,the shadows of cars, buildings and the inside of windows in cityscape images.
colorful objects or surfaces. Lacking color in any color channel will result in low values in the dark channel.
dark objects or surfaces. e.g., dark tree trunk and stone.

暗通道先验的证明

暗通道先验的证明是来源于flickr.com上的5000张无雾图片(haze-free landscape and city scape pictures)，手动剪切掉了天空区域，并且调整成500×500的大小。进行暗通道计算后，做亮度直方图(intensity histogram)统计，发现75%的像素的亮度为0,90%的像素的亮度值小于25。

潜在的无效场景

The dark channel prior may be invalid when the scene object is inherently similar to the airlight over a large local region and no shadow is cast on the object.

如果目标场景内在的就和大气光类似，比如雪地、白色背景墙、大海等，则由于前提条件就不正确，因此一般无法获得满意的效果，而对于一般的风景照片这个算法能处理的不错。

去雾算法其他背景知识

在计算机视觉和计算机图形中，下述方程所描述的雾图形成模型被广泛使用：

$${\bf I}(x) = {\bf J}(x)t(x) + {\bf A}(1-t(x)), (3) $$

${\bf I}(x)$ 为雾化图像(observed intensity)
${\bf J}(x)$ 为恢复图像(scene radiance)
$t(x)$ 为透射率(t is the medium transmission describing the portion of the light that is not scattered and reaches the camera)
${\bf A}$ 为全局大气光线(global atmospheric light)

等式右边的第一项${\bf J}(x)t(x)$被称为直接衰减(direct attennuation)，另一项则是大气光线(airlight)，这两者的叠加构成了我们看到的雾化图像。

在下图中对(3)式中各项做了很好的解释

公式推导

将(3)式两边同处以${\bf A}$，可得

$$\frac{ {\bf I}^c }{ {\bf A}^c }
= \frac{ {\bf J}^c }{ {\bf A}^c }
t + 1-t ,(4)$$

将(4)式两边同时计算暗通道，可得

$$ \min_{ \Omega }(\min_{ c }( \frac{ {\bf I}^c }{ {\bf A}^c } )) =
\left(
\min_{ \Omega }(\min_{ c }( \frac{ {\bf J}^c }{ {\bf A}^c } ))
\right)
t + 1 - t ,(5)$$

应用(2)式暗通道先验，可推导得

$$ \left(
\min_{ \Omega }(\min_{ c }( \frac{ {\bf J}^c }{ {\bf A}^c } ))
\right) \rightarrow 0 ,(6) $$

将(6)式代入(5)式,整理后可得透射率估算公式

$$t =1-\min_{ \Omega }(\min_{ c }( \frac{ {\bf I}^c }{ {\bf A}^c } )),(7)$$

获取透射率图像以后经过计算即可获得恢复图像（scene radiance）${\bf J}(x)$

$${\bf J}(x) = \frac{ { {\bf I}(x) } - { \bf A } }{ t(x) } + {\bf A},(8)$$

$\omega$修正

在现实生活中，即使是晴天白云，空气中也存在着一些颗粒，因此，看远处的物体还是能感觉到雾的影响，另外，雾的存在让人类感到景深的存在，因此，有必要在去雾的时候保留一定程度的雾，这可以通过在式(7)中引入一个在[0,1]之间的因子$\omega$，则式(7)修正为：

$$t =1- \omega \min_{ \Omega }(\min_{ c }( \frac{ {\bf I}^c }{ {\bf A}^c } )),(9)$$

这里的$\omega$可以理解为去雾程度，1为完全去雾，0为不去雾。

透射率图的精化

透射率图(transmission)的精化有助于减弱人工光晕(halo artifacts)的现象。具体的方法而言，可以用原文中提到的soft matting方法，还有何恺明先生后来的 guided image filtering 导向滤波方法。

在我的实现中，暂时没有做透射率图精化，所以这两项先保留。

soft matting

guided image filtering

关于全局大气光值的估算

在实际中，我们可以借助于暗通道图来从有雾图像中获取该值。具体步骤如下：

从暗通道图中按照亮度的大小取前0.1%的像素。
在这些位置中，在原始有雾图像I中寻找对应的具有最高亮度的点的值，作为A值。

以下图为例说明一下这么做的理由：如果使用传统的方法，直接选取图像中的亮度值最高的点作为全局大气光值，这样原始有雾图像中的白色物体会对此有影响，使得其值偏高。暗通道的运算可以抹去原始图像中小块的白色物体，所以这样估计的全局大气光值会更准确。

算法实现

算法流程梳理

计算雾化图像的暗通道
估算全局大气光值
按式(7)估算透射率图
利用导向滤波获取更精细的透射率图
按式(8)恢复为去雾图像

暗通道的计算

暗通道的计算主要分成两个步骤，第一是获取BGR三个通道的最小值，第二是以一个窗口做MinFilter。这里窗口大小一般为15（radius为7）。

获取BGR三个通道的最小值就是遍历整个图像，取最小值即可。

# 获取最小值矩阵
def getMinChannel(img):
	# 输入检查
	if len(img.shape)==3 and img.shape[2]==3:
		pass
	else:
		print "bad image shape, input must be color image"
		return None
	imgGray = np.zeros((img.shape[0],img.shape[1]),dtype = np.uint8)
	localMin = 255
	for i in range(0,img.shape[0]):
		for j in range(0,img.shape[1]):
			localMin = 255
			for k in range(0,3):
				if img.item((i,j,k)) < localMin:
					localMin = img.item((i,j,k))
			imgGray[i,j] = localMin
	return imgGray

我对MinFilter的实现是，先将原矩阵按半径扩大，多出来的像素点全部填充255,然后再用四重循环遍历取最小值。如果追求高效的话，可以考虑复杂度为O(1)的MinFilter实现 - STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN THREE COMPARISONS PER ELEMENT

# 获取暗通道
def getDarkChannel(img,blockSize = 3):
	# 输入检查
	if len(img.shape)==2:
		pass
	else:
		print "bad image shape, more than two demensions"
		return None
	# blockSize检查
	if blockSize % 2 == 0 or blockSize < 3:
		print 'blockSize is not odd or too small'
		return None
	# 计算addSize
	addSize = (blockSize-1)/2
	newHeight = img.shape[0] + blockSize - 1
	newWidth = img.shape[1] + blockSize - 1
	# 中间结果
	imgMiddle = np.zeros((newHeight,newWidth))
	imgMiddle[:,:] = 255
	imgMiddle[addSize:newHeight - addSize,addSize:newWidth - addSize] = img
	imgDark = np.zeros((img.shape[0],img.shape[1]),np.uint8)
	localMin = 255
	for i in range(addSize,newHeight - addSize):
		for j in range(addSize,newWidth - addSize):
			localMin = 255
			for k in range(i-addSize,i+addSize+1):
				for l in range(j-addSize,j+addSize+1):
					if imgMiddle.item((k,l)) < localMin:
						localMin = imgMiddle.item((k,l))
			imgDark[i-addSize,j-addSize] = localMin
	return imgDark

估算全局大气光值

在本阶段，先计算有雾图像的暗通道，然后用一个Node的结构记录暗通道图像每个像素的位置和大小，放入list中。

对list进行降序排序。

之后按暗通道亮度前0.1%(用percent参数指定百分比)的位置，在原始有雾图像中查找最大光强值。

[1] 中提到取0.1%的光强的均值，想尝试一下，所以也写了一个均值模式。

# 获取全局大气光强度
def getAtomsphericLight(darkChannel,img,meanMode = False, percent = 0.001):
	size = darkChannel.shape[0]*darkChannel.shape[1]
	height = darkChannel.shape[0]
	width = darkChannel.shape[1]
	nodes = []
	# 用一个链表结构(list)存储数据
	for i in range(0,height):
		for j in range(0,width):
			oneNode = Node(i,j,darkChannel[i,j])
			nodes.append(oneNode)	
	# 排序
	nodes = sorted(nodes, key = lambda node: node.value,reverse = True)
	
	atomsphericLight = 0
	# 原图像像素过少时，只考虑第一个像素点
	if int(percent*size) == 0:
		for i in range(0,3):
			if img[nodes[0].x,nodes[0].y,i] > atomsphericLight:
				atomsphericLight = img[nodes[0].x,nodes[0].y,i]
		return atomsphericLight
	# 开启均值模式
	if meanMode:
		sum = 0
		for i in range(0,int(percent*size)):
			for j in range(0,3):
				sum = sum + img[nodes[i].x,nodes[i].y,j]
		
		atomsphericLight = int(sum/(int(percent*size)*3))
		return atomsphericLight
	# 获取暗通道前0.1%(percent)的位置的像素点在原图像中的最高亮度值
	for i in range(0,int(percent*size)):
		for j in range(0,3):
			if img[nodes[i].x,nodes[i].y,j] > atomsphericLight:
				atomsphericLight = img[nodes[i].x,nodes[i].y,j]
	return atomsphericLight

估算透射率图

根据公式(7)即可

imgGray = getMinChannel(img)
imgDark = getDarkChannel(imgGray, blockSize = blockSize)
atomsphericLight = getAtomsphericLight(imgDark,img,meanMode = meanMode,percent= percent)
imgDark = np.float64(imgDark)
transmission = 1 - omega * imgDark / atomsphericLight

按式(8)恢复为去雾图像

当投射图t 的值很小时，会导致{\bf J}的值偏大，从而使淂图像整体向白场过度，因此一般可设置一阈值$t_0$，当$t$值小于$t_0$时，令$ t = t_0$，本文中所有效果图均以$ t_0 = 0.1 $为标准计算。

因此，最终的恢复公式如下

$${\bf J}(x) = \frac{ { {\bf I}(x) } - { \bf A } }{ \max ( t(x),t_0 ) } + {\bf A},(10)$$

# 恢复原图像
# Omega 去雾比例 参数
# t0 最小透射率值
def getRecoverScene(img,omega = 0.95,t0 = 0.1 , blockSize = 15 , meanMode = False,percent = 0.001):
	imgGray = getMinChannel(img)
	imgDark = getDarkChannel(imgGray, blockSize = blockSize)
	atomsphericLight = getAtomsphericLight(imgDark,img,meanMode = meanMode,percent= percent)
	imgDark = np.float64(imgDark)
	transmission = 1 - omega * imgDark / atomsphericLight
	# 防止出现t小于0的情况
	# 对t限制最小值为0.1
	for i in range(0,transmission.shape[0]):
		for j in range(0,transmission.shape[1]):
			if transmission[i,j] < 0.1:
				transmission[i,j] = 0.1
	sceneRadiance = np.zeros(img.shape)
	for i in range(0,3):
		sceneRadiance[:,:,i] = img[:,:,i]
		sceneRadiance[:,:,i] = (sceneRadiance[:,:,i] - atomsphericLight)/transmission + atomsphericLight
		# 限制透射率 在0～255
		for j in range(0,sceneRadiance.shape[0]):
			for k in range(0,sceneRadiance.shape[1]):
				if sceneRadiance[j,k,i] > 255:
					sceneRadiance[j,k,i] = 255
				if sceneRadiance[j,k,i] < 0:
					sceneRadiance[j,k,i]= 0
	sceneRadiance = np.uint8(sceneRadiance)
	return sceneRadiance