CN102156968B - 一种基于颜色立方先验的单一图像能见度复原方法 - Google Patents

Abstract

本发明公开了一种基于颜色立方先验的单一图像能见度复原方法，属于图像信息处理领域。该方法建立在物理模型的基础上，包括以下各步骤：获取原始输入图像I^c(x)；获取原始输入图像I^c(x)每一坐标邻域内的实际分布特征；对图像采集环境中的介质层颜色进行估计，得到介质层颜色A^c；基于颜色立方先验原理，获取表征介质层光学厚度的传输率估计图

对传输率估计图

进行精细化处理，得到传输率图t(x)；提取能见度复原后的图像J^c(x)；本发明方法仅需单一图像作为输入，能够对图像采集环境中介质层的颜色自适应处理，并输出高质量的能见度复原图像，稳定性强、应用范围广。

Description

一种基于颜色立方先验的单一图像能见度复原方法

技术领域

本发明涉及单一图像能见度复原方法，尤其涉及一种基于图像退化的物理模型对单一图像进行复原的方法，属于图像信息处理技术领域。

背景技术

许多应用领域，如机器视觉、消费摄影、海洋工程、远程监控等，需要以高质量的户外图像作为输入，但实际采集到的图像可能受到环境因素，例如烟雾，海水，以及类似传播介质的影响而能见度低下。这些环境因素会严重的降低图像的品质，具体的说，不仅使图像的对比度降低而且还伴随着图像颜色的失真。重要的是，绝大多数受影响的图像将无法顺利进行后续的其它处理，例如：特征抽取与模式识别等。因此，需要对这些低质量的图像进行一定的预处理，以恢复图像的细节特征以及其他有价值的信息。

对于图像能见度的复原问题，目前主要有两类方法：基于图像增强的方法与基于物理模型的方法。图像增强的方法如对比度最大化，直方图均衡化等，能在一定程度上改善图像细节，但图像增强的方法并不直接遵循物理本质，往往造成处理结果的失真。相比之下，基于物理模型的方法，针对图像退化的光学本质，能够取得更为自然的处理结果。

根据光学原理，景物反射的光线在进入采集设备前，不仅被介质层中的悬浮微粒散射了一部分，而且还与介质层的颜色进行了一定比例的混合。直观地说，受介质层影响越大的区域，其颜色与介质层越接近，例如浓雾中的景物通常显现为灰白色，另外，在水下观察到的物体通常偏蓝。图像受到采集环境中介质层影响的过程可以视为一种图像退化，图像处理及相关领域将其建模为式(0)，如下：

I(x)＝J(x)t(x)+A(1-t(x))    (0)

式(0)是定义在图像的RGB颜色通道上的，其中x是图像的二维坐标，I(x)是实际采集到的图像，J(x)是实际景物的图像，A是介质层颜色，或称为大气光颜色(Airlight Color)，t(x)是传输率图(transmission)。归纳如下：采集到的图像I(x)是真实景物的图像J(x)与一个介质层颜色A的线性组合，当t(x)越大时采集到的图像越接近原始景物的图像，反之，当t(x)越小时，采集到的图像越接近介质层颜色。

基于物理模型的单一图像的能见度复原方法就是要根据采集到的输入图像I(x)，分别获取介质层颜色A与传输率图t(x)，并代入物理模型恢复出所需的实际景物图像J(x)，从而达到能见度复原的目的。其中，估计传输率图t(x)是整个方法中的核心步骤，一般根据图像各处的统计特征信息来推断。现有方法中K M He(何恺明)等人提出的方法具有一定的代表性，可以参考论文“Single Image Haze Removal Using Dark Channel Prior[A]”，CVPR[C]，IEEEComputer Society，p1956-1963，June 2009.(2009年6月发表于CVPR会议的《基于暗原色先验原理的单一图像去雾方法》)，该论文的核心方法是以图像的暗通道(dark channel)作为估计传输率图t(x)的主要依据，该方法暗含了一种假设，即介质层颜色A在任意通道上的灰度值不能过低，当这一条件不成立时，例如图像采集的环境存在饱和度较高的彩色烟雾或图像在水下采集等，该方法可能失效。对于水下图像，Nicholas Carlevaris-Bianco等人提出了一种专用的方法，可以针对在海水或类似的介质层中采集的图像进行能见度复原，具体可参考论文“Initial Results in Underwater Single Image Dehazing[A]”，OCEANS[C]，IEEE ComputerSociety，p1-8，September 2010(2010年9月发表于IEEE OCEANS会议的《水下单一图像去雾的初步成果》)，在该论文的方法中输入图像的红色通道具有特殊的地位，实验表明，当介质层颜色与蓝色的色差较大时，该方法不能取得理想的结果。

综上，现有的基于物理模型的单一图像能见度复原方法对于图像采集环境中介质层颜色的适应能力有限，应用场合较为单一。

发明内容

本发明所要解决的技术问题在于提供一种基于颜色立方先验的单一图像能见度复原方法，以便使用本发明方法能够对图像采集环境中介质层颜色自适应处理，从而能够在物理模型的约束下将传播介质层所造成的影响从输入图像中剥离，实现对输入图像的能见度复原。利用对介质层颜色的自适应处理，使本发明方法能获得更强的稳定性，应用范围更加广泛。

为实现所述发明目的，本发明采用下述的技术方案。

本发明基于颜色立方先验的单一图像能见度复原方法，其特点是按如下步骤进行：

步骤1、获取原始输入图像I^c(x)；

其中c表示图像的通道，x为图像的二维坐标；

步骤2、获取原始输入图像I^c(x)每一坐标邻域内的实际分布特征；

所述实际分布特征包括：

(a)表征分布最小值的腐蚀图像E^c(x)，

(b)表征分布最大值的膨胀图像D^c(x)，

(c)表征分布加权均值的参考图像R^c(x)；

步骤3、对图像采集环境中的介质层颜色进行估计，得到介质层颜色A^c；

步骤4、基于颜色立方先验原理，获取表征介质层光学厚度的传输率估计图

步骤5、对传输率估计图

进行精细化处理，得到传输率图t(x)；

步骤6、按式(2)提取能见度复原后的图像J^c(x)；

$J^c\left(x\right)=\frac{I^c\left(x\right)-A^c}{t\left(x\right)}+A^c$ c∈{r，g，b}  (2)；

所述步骤2按以下过程进行：

步骤21、对原始输入图像I^c(x)采用形态学腐蚀处理，获得腐蚀图像E^c(x)；

步骤22、对原始输入图像I^c(x)采用形态学膨胀处理，获得膨胀图像D^c(x)；

步骤23、对原始输入图像I^c(x)采用均值滤波法处理，得到均值图像

步骤24、对原始输入图像I^c(x)与均值图像

进行加权平均，得到参考图像R^c(x)；

所述步骤4按如下过程进行：

首先，根据步骤3中所述的介质层颜色A^c与步骤2中获取的参考图像R^c(x)，确定当原始输入图像I^c(x)满足颜色立方先验原理时所应有的理想分布特征；

然后，在图像退化的物理模型框架下，根据步骤2中确定的实际分布特征与所述理想分布特征之间的差异，确定传输率初始图t_ini(x)；

最后，对传输率初始图t_ini(x)采用形态学开处理，得到传输率估计图

本发明方法的特点也在于：

所述步骤24中参考图像R^c(x)按式(6)确定：

$R^c\left(x\right)=w_1\·I^c\left(x\right)+w_2\·I_{\mathrm{mean}}^c\left(x\right)$ c∈{r，g，b}  (6)，

式(6)中w₁的取值范围是0.7-0.9，w₂＝1-w₁。

3、根据权利要求1所述的方法，其特征在于：

所述步骤3按如下步骤31-38进行：

步骤31、按式(7)用膨胀图像D^c(x)减去腐蚀图像E^c(x)获得梯度图像G^c(x)；

G^c(x)＝D^c(x)-E^c(x)c∈{r，g，b}  (7)；

步骤32、在梯度图像G^c(x)的R、G、B通道上取最大值，按式(8)获得最大梯度图像G_max(x)；

$G_{\max }\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{G}^c\left(x\right)---\left(8\right);$

步骤33、对原始输入图像I^c(x)的R、G、B分量图像做加权平均，按式(9)获得强度图像L(x)；

$L\left(x\right)=\frac{I^r\left(x\right)+I^g\left(x\right)+I^b\left(x\right)}{3}---\left(9\right);$

式(7)中I^r(x)、I^g(x)、I^b(x)分别是原始输入图像的R、G、B分量图像；

步骤34、以强度图像L(x)作为引导图像，对最大梯度图像G_max(x)作引导滤波，获得精细梯度图像G_fine(x)；

步骤35、对精细梯度图像G_fine(x)采用形态学膨胀处理，按式(13)获得膨胀梯度图像G_dilate(x)；

$G_{\mathrm{dilate}}\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\max }{G}_{\mathrm{fine}}\left(y\right)$ c∈{r，g，b}    (13)；

步骤36、对精细梯度图像G_fine(x)与膨胀梯度图像G_dilate(x)进行加权平均，按式(14)得到评价图像G_rank(x)；

G_rank(x)＝k₁·G_fine(x)+k₂·G_dilate(x)    (14)，

式(14)中k₁的取值范围为0.3-0.7，k₂＝1-k₁；

步骤37、对所述评价图像G_rank(x)进行二值化处理，按式(15)得到标记图M(x)；

$M\left(x\right)=\left\{\begin{array}{c}0,G_{\mathrm{rank}}\left(k\right)>\mathrm{Threshole}\\ 1,G_{\mathrm{rank}}\left(x\right)\≤\mathrm{Threshole}\end{array}\right.---\left(15\right),$

式(15)中二值化的阈值Threshold按式(16)确定：

Threshold＝min{z|G_rank(x)中有q％元素的灰度值小于等于Z }(16)，

式(16)中q的取值范围为0.01-2.0；

步骤38、对原始输入图像I^c(x)对应标记图M(x)中非零位置的像素求均值，按式(17)获取介质层颜色A^c；

$A^c=\frac{\underset{x}{\mathrm{\Σ}}(M\left(x\right)\·I^c\left(x\right))}{\underset{x}{\mathrm{\Σ}}M\left(x\right)}$ c∈{r，g，b}     (17)。

本发明方法的特点也在于：

所述步骤4中传输率初始图t_ini(x)，按如下步骤41-42确定：

步骤41、根据分布特征的差异，获取传输率候选图

分为以下三种情况进行讨论：

情况(i)，当A^c＝Imin时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{D^c\left(x\right)-I\min }{I\max -I\min },$ c∈{r，g，b}    (18)；

情况(ii)，当A^c＝Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{I\max -E^c\left(x\right)}{I\max -I\min },$ c∈{r，g，b}    (19)；

情况(iii)，当Imin＜A^c＜Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\left\{\begin{array}{c}\frac{D^c\left(x\right)-A^c}{I\max -A^c},A^c\≤R^c\left(x\right)\\ \frac{A^c-E^c\left(x\right)}{A^c-I\min },A^c>R^c\left(x\right)\end{array}\right.,$ c∈{r，g，b}    (20)；

式(18)至(20)中Imin为图像灰度取值范围的最小值，Imax为图像灰度取值范围的最大值；

步骤42、在传输率候选图

的R、G、B通道上取最大值，按式(21)获得传输率初始图t_ini(x)；

$t_{\mathrm{ini}}\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{t}_{\mathrm{can}}^c\left(x\right)---\left(21\right).$

本发明方法的特点还在于：

所述步骤5中的精细化处理按如下步骤51-53进行：

步骤51、将原始输入图像I^c(x)由RGB图像空间转换为YCrCb图像空间，得到Y分量图像I_Y(x)、Cr分量图像I_Cr(x)、Cb分量图像I_Cb(x)；

步骤52、对传输率估计图

按顺序以所述Cr、Cb、Y分量图像为引导图像进行引导滤波，按如下过程得到第三次精细图像

首先，以所述Cr分量图像I_Cr(x)作为引导图像，对所述传输率估计图

进行引导滤波，得到第一次精细图像

然后，以所述Cb分量图像I_Cb(x)作为引导图像，对所述第一次精细图像

进行引导滤波，得到第二次精细图像

最后，以所述Y分量图像I_Y(x)作为引导图像，对所述第二次精细图像

进行引导滤波，得到第三次精细图像

步骤53、对第三次精细图像

做门限处理，获得传输率图t(x)；

$t\left(x\right)=\max \{{\tilde{t}}_3\left(x\right),\mathrm{LB}\},$ 其中LB的取值范围为0.01-0.1。

与现有技术相比，本发明有益效果体现在：

1、本发明基于颜色立方先验原理获取表征介质层光学厚度的传输率估计图，该方法本质上以图像各处的实际分布特征与理想分布特征之间的差异作为主要依据获取传输率估计图，正是因为所述理想分布特征是根据介质层颜色动态确定的，所以该方法能够对图像采集环境中介质层颜色自适应，具体的说，无论图像的采集环境是在陆地(介质层颜色接近白色)还是在水下(介质层颜色接近蓝色)或者其他更复杂的环境，该方法均能准确的对光学模型中最为关键的传输率图进行估计。通过有关实验，验证了本方法的有效性。在实验中以实际传输率图已知同时介质层颜色不同的合成图像作为实验图像，分别使用本发明的方法与现有技术获取实验图像的传输率估计图，实验结果表明现有技术只能在有限范围的介质层颜色中取得理想的结果，而本发明的方法表现较为稳定，所获取的传输率估计图与实际传输率图之间的平均绝对误差(MAE)在所有实验图像上的均值为6.726％，最大值仅为8.908％，充分说明相较于现有技术本发明的方法更通用，能够对图像采集环境中介质层颜色自适应处理，从而使得本发明整体的稳定性更强，应用范围更广。

2、本发明通过获取梯度图像为依据估计介质层颜色，该梯度图像实际上是原始输入图像的形态学梯度信息，因为提取形态学梯度信息本身并不受图像采集环境中介质层颜色范围的制约，所以基于梯度图像能够相对稳定的估计介质层颜色，与现有的基于暗原色原理或者基于图像亮度筛选的方法相比，本发明的适用面更广。

3、本发明分别利用原始输入图像的Cr、Cb、Y分量图像对传输率估计图进行精细化，能够充分的利用原始输入图像的色度与亮度信息。有关实验验证了该精细化步骤的有效性。对传输率估计图进行精细化后与实际传输率图之间的平均绝对误差(MAE)在所有实验图像上的均值由原先的6.726％下降为5.161％，说明本方法能够通过精细化步骤进一步减小传输率估计图的误差，从而使得本发明能够更为精确的对图像进行能见度复原。

4、本发明的方法基于物理模型，遵循着光在介质层中的传播特性进行，能够将图像采集环境中传播介质层中的悬浮微粒所造成的散射等影响从输入图像中剥离，从而实现对输入图像的能见度复原，由于基于物理模型的方法针对的是图像退化的光学本质，所以相对于另一类基于图像增强的方法能够取得更为自然的能见度复原效果。

附图说明

图1本发明的流程简图；

图2本发明中步骤2的流程简图；

图3本发明中步骤3的流程简图；

图4本发明中步骤4的流程简图；

图5本发明中步骤5的流程简图；

图6本发明方法与其他现有技术对水下图像的处理结果图；

图7本发明方法对于红色烟雾的处理结果图；

图8本发明方法对于普通薄雾下图像的处理结果图；

图9示出了用于合成实验图像的原始数据图像；

图10示出了本发明与现有技术对实验图像获取传输率估计图的平均绝对误差散点图；

图11示出了本发明与现有技术对实验图像获取传输率估计图的直观比较图；

图12示出了本发明所获取的传输率估计图精细化前后的平均绝对误差散点图；

图13示出了图像能见度复原装置的框图；

具体实施方式

本发明所提供的方法主要用于处理户外采集的图像数据。根据光在介质层中的传播特性，景物反射的光线在到达采集设备前将被介质层中的悬浮微粒散射一部分，同时介质层反射的光线也以一定比例掺杂其中一同进入采集设备。受此类影响而退化的图像具有相应的物理模型，可表达为式(1)，如下：

I^c(x)＝J^c(x)t(x)+A^c(1-t(x))c∈{r，g，b}    (1)

式(1)中，x为图像的二维坐标，c为图像的通道，I^c(x)是实际采集到的图像，将用作输入图像，J^c(x)是实际景物的图像，是期望求得的能见度复原图像，A^c是介质层颜色，t(x)是传输率图(transmission)，根据朗伯-比尔(Lambert-Beer)定律，t(x)等于e^-βd(x)，其中d(x)是景深，即景物反射的光至观察者所通过的路径长度，β是介质的衰减系数，可以参考S Narasimhan等人的论文“Vision and the Atmosphere[J]”，IJCV，Kluwer Academic Publishers，vol48(3)，p233-254，Jul 2002(2002年7月发表于IJCV期刊的《视觉与大气研究》)，该论文对模型有更详细的探讨。本发明对单一图像的能见度复原是基于式(1)所表达的物理模型，根据实际采集到的图像I^c(x)，分别获取介质层颜色A^c与传输率图t(x)，然后代回式(1)中解出实际景物的图像J^c(x)作为能见度复原后的图像。

为实现单一图像能见度复原的目的，本实施例按如下过程进行(如图1所示)：

步骤1、获取原始输入图像I^c(x)；

其中c表示图像的通道，x为图像的二维坐标；

步骤2、获取原始输入图像I^c(x)每一坐标邻域内的实际分布特征；

所述实际分布特征包括：

(a)表征分布最小值的腐蚀图像E^c(x)，

(b)表征分布最大值的膨胀图像D^c(x)，

(c)表征分布加权均值的参考图像R^c(x)；

步骤3、对图像采集环境中的介质层颜色进行估计，得到介质层颜色A^c；

步骤4、基于颜色立方先验原理，获取表征介质层光学厚度的传输率估计图

步骤5、对传输率估计图

进行精细化处理，得到传输率图t(x)；

步骤6、根据物理模型，按式(2)提取能见度复原后的图像J^c(x)；

$J^c\left(x\right)=\frac{I^c\left(x\right)-A^c}{t\left(x\right)}+A^c$ c∈{r，g，b}    (2)。

步骤1的实现方式为：

通过摄影或摄像设备获取原始输入图像I^c(x)，其中c表示图像的通道，c∈{r，g，b}，x为图像的二维坐标；图像通常由大量像素组成，每个像素具有灰度值；为了支持后续步骤，不妨设原始输入图像纵向的像素数量为H，横向的像素数量为W，图像灰度取值范围的最小值为Imin，最大值为Imax。

步骤2中的邻域是指以图像的二维坐标x为中心的正方形模板区域Ω(x)，模板宽度取为参照边长S；

其中f的取值范围为40-60，为向下取整操作；不妨设|Ω|为Ω(x)内的像素数量。

步骤2是按步骤21-24进行(如图2所示)：

步骤21、对原始输入图像I^c(x)采用形态学腐蚀处理，获得腐蚀图像E^c(x)；

$E^c\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\min }{I}^c\left(y\right)$ c∈{r，g，b}    (3)；

步骤22、对原始输入图像I^c(x)采用形态学膨胀处理，获得膨胀图像D^c(x)；

$D^c\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\min }{I}^c\left(y\right)$ c∈{r，g，b}    (4)；

步骤23、对原始输入图像I^c(x)采用均值滤波法处理，得到均值图像

$I_{\mathrm{mean}}^c\left(x\right)={\left|\mathrm{\Ω}\right|}^{-1}\underset{y\∈\mathrm{\Ω}\left(x\right)}{\mathrm{\Σ}}{I}^c\left(y\right)$ c∈{r，g，b}    (5)；

步骤24、对原始输入图像I^c(x)与均值图像进行加权平均，得到参考图像R^c(x)；

$R^c\left(x\right)=w_1\·I^c\left(x\right)+w_2\·I_{\mathrm{mean}}^c\left(x\right)$ c∈{r，g，b}    (6)，

式(6)中w₁的取值范围是0.7-0.9，w₂＝1-w₁。

步骤3按如下步骤31-38进行(如图3所示)：

步骤31、用膨胀图像D^c(x)减去腐蚀图像E^c(x)，获得梯度图像G^c(x)；

G^c(x)＝D^c(x)-E^c(x)c∈{r，g，b}    (7)；

梯度图像G^c(x)实际是原始输入图像的形态学梯度信息，对其进一步处理可以得到能够稳定的表征图像各处的平滑程度的评价图G_rank(x)，由于受环境介质层影响严重的区域通常较为平滑，因此可以选取原始输入图像相对平滑的区域作为估计介质层颜色的主要依据。因为提取形态学梯度信息本身并不受图像采集环境中介质层颜色范围的制约，所以基于梯度图像能够相对稳定的估计介质层颜色，与现有的基于暗原色原理或者基于图像亮度筛选的方法相比，适用面更广。

步骤32、在梯度图像G^c(x)的R、G、B通道上取最大值，获得最大梯度图像G_max(x)；

$G_{\max }\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{G}^c\left(x\right)---\left(8\right);$

步骤33、对原始输入图像I^c(x)的R、G、B分量图像做加权平均，获得强度图像L(x)；

$L\left(x\right)=\frac{I^r\left(x\right)+I^g\left(x\right)+I^b\left(x\right)}{3}---\left(9\right);$

式(7)中I^r(x)、I^g(x)、I^b(x)分别是原始输入图像的R、G、B分量图像；

步骤34、以强度图像L(x)作为引导图像，对最大梯度图像G_max(x)作引导滤波，获得精细梯度图像G_fine(x)，其中所述引导滤波的滤波半径r取值为4·S，规格化参数ε取值为0.03；

引导滤波可以参考论文“Guided Image Filtering[A]”，11^th ECCV[C]，p1-14，September2010(2010年9月发表于第11届ECCV会议的《图像的引导滤波》)，引导滤波的具体实现方法引述如下：

设：pinput(x)是引导滤波的输入图像，I_guidance(x)是引导滤波的引导图像，滤波半径为r，规格化参数为ε；

首先，按式(10)确定临时矩阵a(x)：

$a\left(x\right)=\frac{{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{I}_{\mathrm{guidance}}\left(y\right)p_{\mathrm{input}}\left(y\right)-{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{I}_{\mathrm{guidance}}\left(y\right)\·{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{p}_{\mathrm{input}}\left(y\right)}{{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{\left(I_{\mathrm{guidance}}\left(y\right)\right)}^2-{\left({\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{I}_{\mathrm{guidance}}\left(y\right)\right)}^2+\mathrm{\ϵ}}---\left(10\right),$

然后，按式(11)确定另一临时矩阵b(x)：

$b\left(x\right)={\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{p}_{\mathrm{input}}\left(y\right)-a\left(x\right)\·{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}{I}_{\mathrm{guidance}}\left(y\right)---\left(11\right),$

最后，按式(12)确定引导滤波的输出图像q_output(x)：

$q_{\mathrm{output}}\left(x\right)={\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}a\left(y\right)\·I_{\mathrm{guidance}}\left(x\right)+{\left|w\right|}^{-1}\underset{y\∈w\left(x\right)}{\mathrm{\Σ}}b\left(y\right)---\left(12\right),$

式(10)至式(12)中w(x)均为以x为中心的正方形模板，模板宽度取值为2.r+1，|w|是w(x)内的像素数量；

步骤35、对精细梯度图像G_fine(x)采用形态学膨胀处理，获得膨胀梯度图像G_dilate(x)；

$G_{\mathrm{dilate}}\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\max }{G}_{\mathrm{fine}}\left(y\right)$ c∈{r，g，b}    (13)；

步骤36、对精细梯度图像G_fine(x)与膨胀梯度图像G_dilate(x)进行加权平均，得到评价图像G_rank(x)；

G_rank(x)＝k₁·G_fine(x)+k₂·G_dilate(x)    (14)，

式(14)中k₁的取值范围为0.3-0.7，k₂＝1-k₁；

步骤37、对所述评价图像G_rank(x)进行二值化处理，得到标记图M(x)；

$M\left(x\right)=\left\{\begin{array}{c}0,G_{\mathrm{rank}}\left(k\right)>\mathrm{Threshole}\\ 1,G_{\mathrm{rank}}\left(x\right)\≤\mathrm{Threshole}\end{array}\right.---\left(15\right),$

式(15)中二值化的阈值Threshold按如下方式确定：

Threshold＝min{Z|G_rank(x)中有q％元素的灰度值小于等于Z}  (16)，

式(16)中q的取值范围为0.01-2.0；

步骤38、对原始输入图像I^c(x)对应标记图M(x)中非零位置的像素求均值，获取介质层颜色A^c；

$A^c=\frac{\underset{x}{\mathrm{\Σ}}(M\left(x\right)\·I^c\left(x\right))}{\underset{x}{\mathrm{\Σ}}M\left(x\right)}$ c∈{r，g，b}    (17)。

步骤4按如下过程进行：

首先，根据步骤3中的介质层颜色A^c与步骤2中获取的参考图像R^c(x)，确定当原始输入图像I^c(x)满足颜色立方先验原理时所应有的理想分布特征；

然后，在图像退化的物理模型框架下，根据步骤2中确定的实际分布特征与所述理想分布特征之间的差异，确定传输率初始图t_ini(x)；

最后，对传输率初始图t_ini(x)采用形态学开处理，得到传输率估计图

颜色立方先验原理为自然图像在未受到环境介质层的影响时任意坐标处的邻域内存在一个像素落在颜色立方的边界上。更具体的说，满足该原理的图像每一坐标邻域内的具有如下理想分布特征：至少在一个颜色通道上，其表征分布最小值的像素灰度值为Imin，或者其表征分布最大值的像素灰度值为Imax。根据图像退化的光学模型，图像在受到环境介质层的影响后，局部的加权均值将由单一的方向靠近介质层颜色，因此在某一颜色通道上，若介质层颜色的灰度值小于等于分布加权均值则所述理想分布特征可以进一步确定为分布最大值的像素灰度值为Imax，反之则为分布最小值的像素灰度值为Imin。可见，上述理想分布特征是根据介质层颜色动态确定的。将实际分布特征与理想分布特征代入光学模型，即可获得在每个通道上的传输率候选图，由于在R、G、B通道上只需存在一个通道满足分布最值特征，因此选取传输率候选图R、G、B通道上的最大值作为传输率初始图。之后再对传输率初始图进行一定后处理，即形态学开处理，得到传输率估计图本实施例将进一步对基于颜色立方先验原理获取表征介质层光学厚度的传输率估计图的有益效果进行量化分析。

图4所示为本发明中步骤4的流程简图。

步骤4中传输率初始图t_ini(x)按如下步骤41-42确定：

步骤41、根据分布特征的差异，获取传输率候选图

分为如下三种情况进行讨论：

情况(i)，当A^c＝Imin时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{D^c\left(x\right)-I\min }{I\max -I\min },$ c∈{r，g，b}    (18)；

情况(ii)，当A^c＝Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{I\max -E^c\left(x\right)}{I\max -I\min },$ c∈{r，g，b}    (19)；

情况(iii)，当Imin＜A^c＜Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\left\{\begin{array}{c}\frac{D^c\left(x\right)-A^c}{I\max -A^c},A^c\≤R^c\left(x\right)\\ \frac{A^c-E^c\left(x\right)}{A^c-I\min },A^c>R^c\left(x\right)\end{array}\right.,$ c∈{r，g，b}    (20)；

式(18)至(20)中Imin为图像灰度取值范围的最小值，Imax为图像灰度取值范围的最大值；

步骤42、在传输率候选图

的R、G、B通道上取最大值，获得传输率初始图t_ini(x)；

$t_{\mathrm{ini}}\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{t}_{\mathrm{can}}^c\left(x\right)---\left(21\right).$

所述步骤4中传输率估计图

按如下步骤43确定：

步骤43、对传输率初始图t_ini(x)采用形态学开处理，得到传输率估计图

首先，对传输率初始图t_ini(x)进行形态学腐蚀处理，得到传输率腐蚀图t_erode(x)；

$t_{\mathrm{erode}}\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\min }{t}_{\mathrm{ini}}\left(y\right)---\left(22\right),$

然后，对传输率腐蚀图t_erode(x)进行形态学膨胀处理，得到传输率估计图

$\tilde{t}\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\max }{t}_{\mathrm{erode}}\left(y\right)---\left(23\right).$

步骤5中的精细化处理按如下步骤51-53进行(如图5所示)：

步骤51、将原始输入图像I^c(x)由RGB图像空间转换为YCrCb图像空间，得到Y分量图像I_Y(x)、Cr分量图像I_Cr(x)、Cb分量图像I_Cb(x)；

I_Y(x)＝0.299·I^r(x)+0.587·I^g(x)+0.114·I^b(x)  (24)，

I_Cr(x)＝(I^r-I_Y)·0.713+Delta                   (25)，

I_Cb(x)＝(I^b-I_Y)·0.564+Delta                   (26)，

其中Delta的值按式(27)确定：

步骤52、对传输率估计图

按顺序以Cr、Cb、Y分量图像为引导图像进行引导滤波，得到第三次精细图像

首先，以Cr分量图像I_Cr(x)作为引导图像，对传输率估计图

进行引导滤波，得到第一次精细图像

其中滤波半径r取值为3·S，规格化参数ε取值为0.01；

然后，以Cb分量图像I_Cb(x)作为引导图像，对第一次精细图像

进行引导滤波，得到第二次精细图像

其中滤波半径r取值为3·S，规格化参数ε取值为0.01；

最后，以Y分量图像I_Y(x)作为引导图像，对第二次精细图像

进行引导滤波，得到第三次精细图像其中滤波半径r取值为3·S，规格化参数ε取值为0.01；

步骤53、对第三次精细图像

做门限处理，获得传输率图t(x)；

$t\left(x\right)=\max \{{\tilde{t}}_3\left(x\right),\mathrm{LB}\},$ 其中LB的取值范围为0.01-0.1。

至此，基于颜色立方先验的单一图像能见度复原方法的流程结束。

以下通过图6、图7和图8所示的内容对本发明技术效果作出更为直观的说明：图6所示为本发明方法与现有技术对水下图像的处理结果。图6中(a1)和(a2)为原始输入图像，(b1)和(b2)为R Fattal提出的方法求得的能见度复原图像，图6中(c1)和(c2)为Nicholas等人提出的方法求取的能见度复原图像，图6中(d1)和(d2)为本发明方法获取的能见度复原图像。通过对比处理结果可见，本方法的最终处理结果更为自然，所产生的失真更少。

图7示出了本发明方法对于红色烟雾的处理结果，图7中(a)为原始输入图像，(b)为本发明方法求取的能见度复原图像。通过处理结果可见船体附近的红色烟雾所造成的能见度下降在处理后得到大幅改善，鉴于彩色烟雾常为军用，可以预期本发明将在军事领域产生一定的应用价值。

图8示出了本发明方法对于普通薄雾下图像的处理结果。图8中(a)为原始输入图像，(b)为本发明方法求取的能见度复原图像。通过实验处理结果可见，本发明方法对于普通薄雾下图像能够起到能见度复原的效果。

以下进一步提供对关键步骤有益效果的量化分析：

为了验证本发明基于颜色立方先验原理获取表征介质层光学厚度的传输率估计图的有益效果，按如下方式进行实验：

首先，获取一组实际传输率图已知并且介质层颜色不同的合成图像作为实验图像；图9所示为用于合成实验图像的原始数据图像。图9中(a)为普通景物图像，(b)为实际传输率图。本实验通过将图9(a)作为实际景物的图像J^c(x)，图9(b)作为传输率图t(x)，再选择一组具有代表性的颜色作为介质层颜色A^c，代入式(1)中，所获得的一组I^c(x)即为实验图像；值得强调的是，这组实验图像的传输率图相同，而介质层颜色不同，并且色差较大，任意两个颜色满足至少在一个通道上相差127；

然后，分别使用本发明方法与现有技术获取实验图像的传输率估计图，这样可以检验一种方法在不同的介质层颜色下获取传输率估计图的准确性，在本实验中通过计算不同方法获取的传输率估计图与实际传输率图的平均绝对误差(MAE)进行度量。

表1

实验结果的数据列于表1，表1中示出了本发明与现有技术对实验图像获取传输率估计图的平均绝对误差。表1中，自左边起的第1列为合成实验图像所使用的介质层颜色(BGR通道的强度值)，自左边起第2列为K M He(何恺明)等人的方法获取的传输率估计图的平均绝对误差，自左边起第3列为Nicholas等人所的方法获取的传输率估计图的平均绝对误差，自左边起第4列为本发明获取的传输率估计图的平均绝对误差。为了进行比较，在图10中示出了本发明与现有技术对实验图像获取传输率估计图的平均绝对误差散点图。图10中X轴为介质层颜色的编号，X＝i处对应了表1第i行的介质层颜色。实验结果表明，K M He(何恺明)等人提出的方法在介质层颜色饱和度较大时无法取得理想的结果，Nicholas等人提出的方法仅在介质层颜色接近蓝色时估计较为准确，而本发明在不同的介质层颜色下均能够稳定而准确的获取传输率估计图，本发明获取的传输率估计图的平均绝对误差在所有实验图像上的均值为6.726％，最大值仅为8.908％，说明相较于现有技术本发明的方法更通用，能够对图像采集环境中介质层颜色自适应处理，从而使得本发明整体的稳定性更强，应用范围更广。

为了更直观地演示实验结果，在图11中示出了本发明与现有技术对实验图像获取传输率估计图的直观比较图。图11中(a1)至(a6)为实验图像，(b1)至(b6)为K M He(何恺明)等人提出的方法所获取的传输率估计图，(c1)至(c6)为Nicholas等人提出的方法求取的传输率估计图，(d1)至(d6)为本发明方法获取的传输率估计图。图11中(a1)至(a6)的介质层颜色分别为：(a1)B：255 G：255 R：255，(a2)B：255 G：0 R：255，(a3)B：0 G：255 R：127，(a4)B：255 G：0 R：127，(a5)B：0 G：127 R：127，(a6)B：255 G：0 R：0；实验图像的实际传输率均为图9(b)。通过实验结果的对比可见，K M He(何恺明)等人提出的方法仅在介质层颜色为白色时，能够取得理想的结果，见图11(b1)，Nicholas等人提出的方法仅在介质层颜色为蓝色时，能够取得理想的结果，见图11(c6)，而本发明方法更通用，能够适用于更广泛的介质层颜色，见图11(d1)至(d6)，介质层颜色的变化对本发明获取取传输率估计图步骤的影响相对较小。

为验证本发明分别利用原始输入图像的Cr、Cb、Y分量图像对传输率估计图进行精细化的有益效果，对同一组实验图像做进一步实验，对本发明的方法所获取的传输率估计图进行精细化，即所述步骤5，获取传输率图，分别计算传输率估计图精细化前后与实际传输率图的平均绝对误差(MAE)。图12示出了本发明所获取的传输率估计图精细化前后的平均绝对误差散点图。图12中X轴为介质层颜色的编号，X＝i处对应了表1第i行的介质层颜色。实验结果显示，对传输率估计图进行精细化后的平均绝对误差在所有实验图像上的均值由原先的6.726％下降为5.161％，说明本方法能够有效地对传输率估计图进行精细化，从而使得本发明能够更为精确的对图像进行能见度复原。

图13所示为一用于实施本发明方法的单一图像能见度复原装置130，该装置包括图像捕获装置131、处理器132、存储器133以及图像输出装置134，其中：

图像捕获装置131可以由公知的摄像装置构成，用于捕获外部的数字图像数据。在本实施例中图像捕获装置131通过USB接口与处理器132相连接。

处理器132用于控制单一图像能见度复原装置130对外部输入的数字图像数据进行处理，所述的处理按图1所示的流程进行。

存储器133可以由公知的易失或者非易失存储芯片构成，能够存储数据信息；本实施例中存储器133存储的数据信息包括以电子信号形式存在的数字图像数据，也包括用于指示处理器132完成图像处理的多个程序模块，其中包括：输入模块，特征提取模块，介质层估计模块，传输率估计模块，精细化模块，复原模块；其中，输入模块通过图像捕获装置131获取原始输入图像I^c(x)；特征提取模块用于获取原始输入图像每一坐标邻域内的实际分布特征；介质层估计模块用于对图像采集环境中的介质层颜色进行估计，得到介质层颜色A^c；传输率估计模块用于基于颜色立方先验原理，获取表征介质层光学厚度的传输率估计图

精细化模块用于对传输率估计图

进行精细化处理，得到传输率图t(x)；复原模块用于提取能见度复原后的图像J^c(x)。

图像输出装置134作为接口与处理器132相连接，可以由公知的显示设备或者通信设备构成，用于将能见度复原后的图像或存储器133所保存的数据进行输出。

Claims (2)

1.一种基于颜色立方先验的单一图像能见度复原方法，其特征在于，按如下步骤进行：

步骤1、获取原始输入图像I^c(x)；

其中c表示图像的通道，x为图像的二维坐标；

步骤2、获取原始输入图像I^c(x)每一坐标邻域内的实际分布特征；

所述实际分布特征包括：

(a)表征分布最小值的腐蚀图像E^c(x)，

(b)表征分布最大值的膨胀图像D^c(x)，

(c)表征分布加权均值的参考图像R^c(x)；

步骤3、对图像采集环境中的介质层颜色进行估计，得到介质层颜色A^c；

步骤4、基于颜色立方先验原理，获取表征介质层光学厚度的传输率估计图

步骤5、对传输率估计图

进行精细化处理，得到传输率图t(x)；

步骤6、按式(2)提取能见度复原后的图像J^c(x)；

$J^c\left(x\right)=\frac{I^c\left(x\right)-A^c}{t\left(x\right)}+A^c$ c∈{r，g，b}    (2)；

所述步骤2按以下过程进行：

步骤21、对原始输入图像I^c(x)采用形态学腐蚀处理，获得腐蚀图像E^c(x)；

步骤22、对原始输入图像I^c(x)采用形态学膨胀处理，获得膨胀图像D^c(x)；

步骤23、对原始输入图像I^c(x)采用均值滤波法处理，得到均值图像

步骤24、对原始输入图像I^c(x)与均值图像

进行加权平均，得到参考图像R^c(x)；

所述步骤4按如下过程进行：

首先，根据步骤3中所述的介质层颜色A^c与步骤2中获取的参考图像R^c(x)，确定当原始输入图像I^c(x)满足颜色立方先验原理时所应有的理想分布特征；

然后，在图像退化的物理模型框架下，根据步骤2中确定的实际分布特征与所述理想分布特征之间的差异，确定传输率初始图t_ini(x)；

最后，对传输率初始图t_ini(x)采用形态学开处理，得到传输率估计图

所述步骤3按如下步骤31-38进行：

步骤31、按式(7)用膨胀图像D^c(x)减去腐蚀图像E^c(x)获得梯度图像G^c(x)；

G^c(x)＝D^c(x)-E^c(x) c∈{r，g，b}        (7)；

步骤32、在梯度图像G^c(x)的R、G、B通道上取最大值，按式(8)获得最大梯度图像G_max(x)；

$G_{\max }\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{G}^c\left(x\right)---\left(8\right);$

步骤33、对原始输入图像I^c(x)的R、G、B分量图像做加权平均，按式(9)获得强度图像L(x)；

$L\left(x\right)=\frac{I^r+I^g+I^b}{3}---\left(9\right);$

式(7)中I^r(x)、I^g(x)、I^b(x)分别是原始输入图像的R、G、B分量图像；

步骤34、以强度图像L(x)作为引导图像，对最大梯度图像G_max(x)作引导滤波，获得精细梯度图像G_fine(x)；

步骤35、对精细梯度图像G_fine(x)采用形态学膨胀处理，按式(13)获得膨胀梯度图像G_dilate(x)；

$G_{\mathrm{dilate}}\left(x\right)=\underset{y\∈\mathrm{\Ω}\left(x\right)}{\max }{G}_{\mathrm{fine}}\left(y\right)c\∈\{r,g,b\}---\left(13\right);$

步骤36、对精细梯度图像G_fine(x)与膨胀梯度图像G_dilate(x)进行加权平均，按式(14)得到评价图像G_rank(x)；

G_rank(x)＝k₁·G_fine(x)+k₂·G_dilate(x)            (14)，

式(14)中k₁的取值范围为0.3-0.7，k₂＝1-k₁；

步骤37、对所述评价图像G_rank(x)进行二值化处理，按式(15)得到标记图M(x)；

$M\left(x\right)=\left\{\begin{array}{c}0,G_{\mathrm{rank}}\left(x\right)>\mathrm{Threshold}\\ 1,G_{\mathrm{rank}}\left(x\right)\≤\mathrm{Threshold}\end{array}\right.---\left(15\right),$

式(15)中二值化的阈值Threshold按式(16)确定：

Threshold＝min{Z|G_rank(x)中有q％元素的灰度值小于等于Z}        (16)，

式(16)中q的取值范围为0.01-2.0；

步骤38、对原始输入图像I^c(x)对应标记图M(x)中非零位置的像素求均值，按式(17)获取介质层颜色A^c；

$A^c=\frac{\underset{x}{\mathrm{\Σ}}(M\left(x\right)\·I^c\left(x\right))}{\underset{x}{\mathrm{\Σ}}M\left(x\right)}$ c∈{r，g，b}                    (17)；

所述步骤4中传输率初始图t_ini(x)，按如下步骤41-42确定：

步骤41、根据分布特征的差异，获取传输率候选图

分为以下三种情况进行讨论：

情况(i)，当A^c＝Imin时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{D^c\left(x\right)-I\min }{I\max -I\min }$ ，c∈{r，g，b}            (18)；

情况(ii)，当A^c＝Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\frac{I\max -E^c\left(x\right)}{I\max -I\min }$ ，c∈{r，g，b}            (19)；

情况(iii)，当Imin＜A^c＜Imax时：

$t_{\mathrm{can}}^c\left(x\right)=\left\{\begin{array}{cc}\frac{D^c\left(x\right)-A^c}{I\max -A^c}& ,A^c{\≤R}^c\left(x\right)\\ \frac{A^c-E^c\left(x\right)}{A^c-I\min }& ,A^c>R^c\left(x\right)\end{array}\right.$ ，c∈{r，g，b}            (20)；

式(18)至(20)中Imin为图像灰度取值范围的最小值，Imax为图像灰度取值范围的最大值；

步骤42、在传输率候选图

的R、G、B通道上取最大值，按式(21)获得传输率初始图t_ini(x)；

$t_{\mathrm{ini}}\left(x\right)=\underset{c\∈\{r,g,b\}}{\max }{t}_{\mathrm{can}}^c\left(x\right)---\left(21\right);$

所述步骤5中的精细化处理按如下步骤51-53进行：

步骤51、将原始输入图像I^c(x)由RGB图像空间转换为YCrCb图像空间，得到Y分量图像I_Y(x)、Cr分量图像I_Cr(x)、Cb分量图像I_Cb(x)；

步骤52、对传输率估计图

按顺序以所述Cr、Cb、Y分量图像为引导图像进行引导滤波，按如下过程得到第三次精细图像

首先，以所述Cr分量图像I_Cr(x)作为引导图像，对所述传输率估计图

进行引导滤波，得到第一次精细图像

然后，以所述Cb分量图像I_Cb(x)作为引导图像，对所述第一次精细图像

进行引导滤波，得到第二次精细图像

最后，以所述Y分量图像I_Y(x)作为引导图像，对所述第二次精细图像进行引导滤波，得到第三次精细图像

步骤53、对第三次精细图像

做门限处理，获得传输率图t(x)；

$t\left(x\right)=\max \{{\tilde{t}}_3\left(x\right),\mathrm{LB}\},$ 其中LB的取值范围为0.01-0.1。

2.根据权利要求1所述的方法，其特征在于：

所述步骤24中参考图像R^c(x)按式(6)确定：

$R^c\left(x\right)=w_1\·I^c\left(x\right)+w_2\·I_{\mathrm{mean}}^c\left(x\right)$ c∈{r，g，b}                (6)，

式(6)中w₁的取值范围是0.7-0.9，w₂＝1-w₁。

Images