基于MATLAB的图像去雾系统

（一）图像去雾基础原理

1、雾霭的形成机理

雾实际上是由悬浮颗粒在大气中的微小液滴构成的气溶胶，常呈现乳白色，其底部位于地球表面，所以也可以看作是接近地面的云。霭其实跟雾区别不大，它的一种解释是轻雾，多呈现灰白色，与雾的颜色十分接近。广义的雾包括雾、霾、沙尘、烟等一切导致视觉效果受限的物理现象。由于雾的存在，户外图像质量降低，如果不处理，往往满足不了相关研究、应用的要求。在雾的影响下，经过物体表面的光被大气中的颗粒物吸收和反射，导致获取的图像质量差，细节模糊、色彩暗淡。

2、图像去雾算法

图像去雾算法可以分为两大类：一类是图像增强；另一类是图像复原。图1-1介绍了图像去雾算法的分类:

图1-1 去雾算法分类

从图像呈现的低亮度和低对比度的特征考虑，采用增强的方法处理，即图像增强。比较典型的有全局直方图均衡化，同态滤波，Retinex 算法，小波算法等等。

（1）图像增强技术

为了改善视觉效果或者便于人们对图像的判别和分析，根据图像的特征采取简单的改善方法或者加强特征的措施叫做图像增强。图像增强可分为两大类：频率域法和空间域法。空间域处理主要包括：点处理，模块处理即领域处理。频率域处理主要包括：高、低通滤波、同态滤波等等。

图像增强可分为两大类：频率域法和空间域法。空间域处理主要包括：点处理，模块处理即领域处理。频率域处理主要包括：高、低通滤波、同态滤波等等。

（2）图像复原技术

从广义上讲，图像复原是一个求逆问题，逆问题经常存在非唯一解，甚至无解。图像复原的目的是将所观测到的退化图像恢复到退化前的原始图像，这种恢复过程在很多图像处理中的应用十分重要。为了更好的对图像复原的理解，图1-2为图像复原的流程图：

图1-2 图像复原流程图

其中g(x,y)为降质图像函数，f(x,y)为真实图像函数。

图像复原技术可以分为以下几类：

1）在给定退化模型条件下，分为无约束和有约束两大类。

2）根据是否需要外界干预，分为自动和交互两大类。

3）根据处理所在的域，分为频率域和空间域。

（二）从图像增强角度去雾

基于直方图均衡化的算法以概率论为基础，用灰度变换达到图像增强的目的，是图像增强中最常用的算法之一。直方图均衡化处理的“中心思想”是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。

1、图像灰度直方图

定义一：一个灰度级在范围[0,L-1]的数字图像，其直方图是一个离散函数

n是图像的像素总数，是滴k个灰度级，。

定义二：一个灰度级在范围[0,L-1]的数字图像，其直方图是一个离散函数

由于的增量是1，直方图可以表示为：

即图像中不同灰度级像素的出现次数。

2、直方图变换的理论基础

设连续图像的概率分布为：

其中r为灰度

其中A为图像的面积。

均衡化过程分析：

设r和s分别表示原图像灰度级和经直方图均衡化后的图像灰度级，为便于讨论，对r和s进行归一化，使：；对于一幅给定的图像，归一化后灰度级分布在范围内。对[0,1]区间内的任意一个r值进行如下变换：

该变换式应满足条件：

1. 对于，有
2. 在区间内

从s到r的反变换用下式表示

r的概率密度为；s的概率密度可由求出

假定变换函数为

式中：w是积分变量，而就是r的累积分布函数。

给出灰度级在图像中出现的概率密度统计在MATLAB中，imhist函数可以显示一幅图像的直方图。其常见调用方法如下：

imhist(I)

其中I是图像矩阵，该函数返回一幅图像，显示I的直方图。

通过把原图像的直方图通过变换函数修正为分布比较均匀的直方图，从而改变图像整体偏暗或整体偏亮，灰度层次不丰富的情况，这种技术叫直方图均衡化。

在MATLAB中，用于直方图均衡化的函数是histeq，它的常见调用方式如下：

J=histep(I)

其中，I为输入的原图像，J是直方图均衡化后的图像。

3、直方图均衡化的算法步骤

直方图均衡化的算法步骤如下：

1）列出原始图像和变换后图像的灰度级：，,其中是灰度级的个数；

2）统计原图像各灰度级的像素个数;

3）计算原始图像直方图：,N为原始图像像素总个数；

4）利用灰度变换函数计算变换后的灰度值，并四舍五入：；

5）确定灰度变换关系,根据此将原图像的灰度值修正为统计变换后各灰度级的像素个数；

6）计算变换后图像的直方图：。

（二）从图像复原角度去雾

说到图像去雾，就不得不提到由何恺明博士等人提出的基于暗通道的图像去雾算法。这个算法因其新颖的思路和理想的效果而广受关注，相关论文也曾于2009年荣获CVPR最佳论文奖，同时也是该奖设立以来，首次由亚洲学者获此殊荣。

随着大气污染的日益严重，设法改善自动获取的图像质量其意义不言而喻。另一方面，随着数码设备的普及，消费类电子产品的市场也催生出许多新的需求，其中人们对所拍照片质量的修正和优化就是一个显而易见的需求。

通过估算参数，反演退化过程，获得退化前的场景清晰图像。场景目标反射光强经过雾区，会受到雾霾颗粒的强散射和吸收作用，到达探测器的光强会受到影响。

1、暗通道概念

在绝大多数非天空的局部区域里，某些像素总会有至少一个颜色通道具有很低的之。换言之，该区域光强度的最小值诗格很小的数。下面给暗通道一个数学定义，对于任意的输入图像J，其暗通道可以用下式表达：

其中表示彩色图像的每个通道，表示以像素为中心的一个窗口。

暗通道先验的理论指出

也就是说以像素点为中心，分别取三个通道内窗口内的最小值，然后再取三个通道的最值作为像素点的暗通道的值

实际生活中造成暗原色中低通道值的因素有很多。例如，汽车、建筑物和城市中玻璃窗户的阴影，或者是树叶、树与岩石等自然景观的投影；色彩鲜艳的物体或表面，在RGB得三个通道中有些通道的值很低（比如绿色的草地、树木等植物，红色或黄色的花朵、果实或者叶子，或者蓝色、绿色的水面）；颜色较暗的物体或者表面，例如灰暗色的树干、石头以及路面。总之，自然景物中到处都是阴影活着彩色，这些经无图像的暗原色总是表现出较为灰暗的状态。

2、暗通道去雾的原理

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

其中，就是现在已经有的待去雾图像，是要恢复的无雾图像，参数是全球大气光成分，为透射率。现在的已知条件就是，要求目标值。根据基本的代数知识可知这是一个有无数解的方程。只有在一些先验信息基础上才能求出定解。

将上式稍作处理，变形为下式：

首先假设在每一个窗口内透射率为常数，将其定义为，并且值已经给定，然后对上式两边求两次最小值运算，得到下式：

根据前述的暗原色先验理论有：

可推导出

把结论带回原式中，得到：

这就是透射率的预估值。

透射率还可表示为,为大气的散射系数，该式表明景物光线是随着景物深度按指数衰减的。

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

上述推论中都是假设全球大气光值是已知的，在实际中，可以借助暗通道图来从有雾图像中获取该值。具体步骤大致为：

1）首先从暗通道图中按照亮度的大小提取最亮的前0.1%像素；

2）在原始有雾图像中寻找对应位置上的具有最高亮度的点的值，并以此作为的值。 当考虑投射图值很小时，会导致的值偏大，从而使图像整体向白场过渡，因此一般可以设置一个阈值，当值小于时，令。因此，最终的图像恢复公式如下：

当透射率图过于粗糙时，对应暗通道图中颜色较深的部分边缘明显不协调，为了获得更为精细的透射率图，何凯明提出了“soft matting”方法，能得到非常细腻的结果，但是该算法的一个致命弱点就是速度慢，所以何在利用导向滤波的方式来获得较好的透射率图过程中使用简单的盒子滤波相应的快速算法。

五、设计步骤

（一）基于直方图均衡化的图像去雾算法

1、转换为灰度图像后对图像进行直方图均衡

由于直方图均衡仅限于灰度图像，所以我们将输入的RGB图像先转变为灰度图像，再进行图像增强。简要过程如图2-1所示。

图2-1 RGB转灰度图像进行直方图均衡化过程

1. 读入图像，将彩色图像进行降维转换成灰度图像；
2. 对灰度图像的直方图进行均衡化处理；
3. 输出均衡化后的灰色图像。

接下来我们考虑想要得到彩色图像，于是从不同空间对图像进行均衡化处理。

2、在RGB空间对图像进行直方图均衡（流程如图2-2）

1. 读入图像，将图像存储于RGB空间；
2. 分别对RGB空间的R、G、B三个分量的灰度直方图进行均衡化处理；
3. 输出图像。

图2-2 在RGB空间对图像进行直方图均衡

3、在HSV空间对图像进行直方图均衡（流程如图2-3）

1. 读入图像，将图像由RGB空间转换到HSV空间；
2. 对HSV空间饱和度和亮度分量（S、V分量）的灰度直方图进行均衡化处理；
3. 将图像由HSV空间转换到RGB空间并输出。

图2-3 在HSI空间对图像进行直方图均衡

4、在YCrCb空间对图像进行直方图均衡（流程如图2-4）

1. 读入图像，将图像由RGB空间转换到YCrCb空间；
2. 对YCrCb空间亮度分量（Y分量）的灰度直方图进行均衡化处理；
3. 将图像由YCrCb空间转换到RGB空间并输出图像。

图2-4 在YCrCb空间对图像进行直方图均衡

（二）基于暗通道先验的图像去雾算法

根据何凯明的暗通道先验算法原理，大致的流程如图2-6所示。

图2-6 暗通道先验过程图

（1）根据原始有雾图像求暗通道

用for循环求出每个像素RGB分量中的最小值，存入一副和原始图像大小相同的灰度图中，然后再对这副灰度图使用ordfilt2进行最小值滤波，滤波的半径由窗口大小决定。

（2）求解全局大气光照

根据原理应该选择暗通道内图像总像素点个数千分之一个最亮的像素点，并记录这些像素点坐标，再根据这些点的坐标分别在原图像I的三个通道内找到这些像素点并加和。

考虑到这个思路在MATLAB实现中在寻找这些像素点时，会找到不止1000点，因为在图像中有很多值相等的像素点；为了方便计算，在求时取最亮值的点集中的某一个点。

（3）求透射率

首先设定原始程序中去雾系数，根据公式先对透射率有一个预估；接下来利用导向滤波（Guided Filter）来对透射率图。

以原始图像的其中一个通道为引导图，使得透射率图的细节更加准确。而为了提高计算速度，在导向滤波中使用盒子滤波器（Box filter），其本质是通过下采样减少像素点，计算后再进行上采样恢复到原有的尺寸大小。

（4）求解无雾图像

根据公式

求解出去雾的图像即可。在初始程序中采用为标准计算。

六、设计结果与分析

（一）基于直方图均衡化的图像去雾算法

直方图均衡化是图像增强中的一种常用方法，该算法以概率论为基础，基于空间不变思想，运用灰度运算来实现直方图的变换，从而达到增强图像的目的，更适合于景物深度变化很小的图像。分为全局直方图均衡化算法和局部直方图均衡化算法（在本次实验中局部均衡方法去雾效果比较差，所以不对该方法进行赘述）。

全局直方图均衡化算法实现简单，对单景深图像的复原效果较好，但对于场景深度多变的图像的复原效果好，但是对于场景深度多变的额图像，这种方法就很难反映图像中局部景物深度的变化。

原始有雾图像如图3-1所示。

(a) 有雾图像1

图3-1 有雾图像

对全局均衡结果图的效果进行分析（程序zft.m）：

对有雾图像RGB通道进行均衡化处理后融合的图像效果最好，清晰度较高，且颜色较符合地物颜色，但天空处有较多的块状噪声，且树叶产生了较大的畸变，结果如图3-1。

图3-1 对RGB通道进行均衡化处理结果

图3-2 对RGB通道进行均衡化结果及其直方图前后对比

将只对HSI空间的亮度分量[I通道]和对HSI空间的亮度、饱和度分量[I、S通道]进行均衡化，通道进行的均衡化处理的图像颜色对比度明显发生的变化，第二种颜色过于鲜艳而发生畸变，第一种较RGB均衡方式的颜色对比度较差，且树叶处畸变较差结果如图3-3。

图3-3 对HSI通道进行均衡化处理结果

图3-4 对HSI通道进行均衡化前后结果对比及其直方图

对YCrCb空间的亮度分量[Y通道]进行均衡处理，所得图像在天空处处理的较差，有大片空白，但是对于树叶的处理较其他几种方式都好，结果如图3-5。

图3-6 对YCrCb通道进行均衡化结果前后对比及其直方图

当考虑将RGB图像转变为灰度图像在进行直方图均衡化的结果如图3-7所示，

图3-7 对灰度进行均衡化结果前后对比及其直方图

图3-8 对灰度进行均衡化结果前后对比及其直方图

如果忽略灰度图直方图均衡化的视觉效果问题，在去雾效果方面灰度变换的结果还是比较好的。

接着再使用有雾图像2的直方图均衡化，结果如图3-9所示，

图3-8 对灰度进行均衡化结果前后对比及其直方图

由于有雾图像2在前景处没有雾的影响，在整体图像增强后虽然也起到了一定去雾的目的，但是图像效果不太好。

总的来说，直方图均衡化的方法对灰度变化不大的树叶和天空处理效果较差，但是对于建筑物等灰度变化较大的区域效果还是不错的。

但是上面的去雾效果仅仅局限于灰度图片，而不能直接对彩色图像进行处理。经过上述直方图处理后，达到均衡化图像的目的，但是图像会在原图基础上失真，即整幅图像的颜色发生改变，有时会严重偏离去雾的初衷。

（二）基于暗通道先验的图像去雾算法

1、实验过程结果图

（1）根据原始有雾图像求暗通道

图4-1 有雾图像1的暗通道图像

图4-2 有雾图像2的暗通道图像

通过两幅图的暗通道图像尤其是有雾图像2中可以看出在近距离的树叶处非常的暗，通过该暗通道图像去估算大气光照和透射率，建立模型就可以很好地在去雾的同时尽可能保留无雾区域的特征细节。

（2）求解全局大气光照

求解有图像1和2的大气光照数A：

2、影响暗通道先验算法去雾效果的参数研究

计算的标准参数值如下表所示：

1	去雾系数ω	0.95
2	滤波窗口大小（形成暗通道的滤波半径）r	9
3	透射率窗口大小（导向半径）R	80
4	透射率下限值	0.1
5	eps

（1）去雾系数ω

粗透射率为

越大，去雾程度越大，图像对比度显得越暗；反之越小，去雾越不彻底，图像越泛白。一般取。

图4-7 图像1的去雾前图像

图4-11 图像1、2去雾前后对比图像

（2）滤波窗口大小x

x越小，选取的细节越多，暗通道构建越细致，使得最终去雾效果越明显；但也由于其包含暗通道的概率越大，暗通道也就越黑，所得去雾图像越暗；

x越大，选取的细节少，暗通道构建越粗糙，导致最终图像越泛白，即去雾的效果越不明显。

图4-12 图像1改变滤波窗口去雾前后对比图像

（三）对比两种去雾算法

	直方图均衡化	暗通道先验
优点	对去雾图像没有要求	不仅在物理上有效，而且能够处理包括在雾浓度很大情况下的远距离物体；可以减小无雾区域的变样
缺点	去雾后的图像或为灰色图像或产生颜色失真，去雾效果不稳定	得到的去雾图像比较暗，无法处理大面积天空或类似区域

直方图均衡化对于灰度图像去雾效果比较好，而暗原色先验对于彩色图像去雾效果比较明显，图像中的特征得到很好的增强，使人们能够很好的观测图像给予的信息，并作出正确的判断。暗原色先验明显的去除雾的干扰,提高图像的清晰度,增强图像色彩和细节,复原得到高质量的图像。

虽然原始图像受雾影响模糊不清、颜色不真实的图像，但经过暗通道先验算法处理后在细节上更清晰，颜色更真实，视觉效果更令人满意。这种算法适合以上各种场景的图像去雾，无论雾的分布均匀与否，雾的浓度薄或浓，场景中深度变坏如何，都能有效地去除图像中雾，有效地提高了图像的清晰度。从而很大程度上提高了雾天的能见度。但是该算法存在一定的局限性，暗原色先验是一种统计的结果，是对大量户外无雾照片的统计结果，如果目标场景内在的就和大气光类似，比如雪地、白色背景墙、大海等，则由于前提条件就不正确，因此一般无法获得满意的效果，而对于一般的风景照片这个算法能处理的不错。另外，由于景物退化与场景深度呈非线性关系，由此带来的最大问题是很难保证建立的景物退化模型的正确性和宽适性。目前，大多数的图像复原方法都建立在大气散射模型的基础上，并受到了此类模型的限制。

七、附录（源代码）

1、aft.m

%利用直方图均衡化去雾

%===========================================================

close all;

clear all;

wu=imread('E:\Report\课程设计\图片\2.bmp');

%--------------------------------------------------------------------------

%对饱和度与亮度进行直方图均衡处理

%h：色度

%s：饱和度

%v：亮度

hsv=rgb2hsv(wu);

wu1=rgb2gray(wu);

h=hsv(:,:,1);

s=hsv(:,:,2);

v=hsv(:,:,3);

S=histeq(s);

V=histeq(v);

result_hsv=hsv2rgb(h,S,V);

result_hsv1=rgb2gray(result_hsv);

%--------------------------------------------------------------------------

% 对RGB每个通道进行直方图均衡处理

% R：红色

% G：绿色

% B：蓝色

result_rgb=cat(3,R,G,B);%将RGB三个分量放入一个矩阵

result_rgb1=rgb2gray(result_rgb);

%--------------------------------------------------------------------------

% 对YCbCr的亮度进行直方图均衡处理

%y：颜色的亮度，

%Cb：蓝色的浓度偏移成分

%Cr：红色的浓度偏移成分

ycbcr=rgb2ycbcr(wu);

y=ycbcr(:,:,1);

cb=ycbcr(:,:,2);

cr=ycbcr(:,:,3);

Y=histeq(y);

result_ycbcr=ycbcr2rgb(cat(3,Y,cb,cr));%%将RGB三个分量放入一个矩阵

result_ycbcr1=rgb2gray(result_ycbcr);

%将彩色图像转换为灰度图像后，再做直方图均衡处理

whole=rgb2gray(wu);

W=histeq(whole);

%--------------------------------------------------------------------------

% 结果显示

figure(1)

subplot(2,2,1),imshow(wu),title('原始图像')

subplot(2,2,2),imshow(result_hsv),title('HSV')

subplot(2,2,3),imhist(wu1);

subplot(2,2,4),imhist(result_hsv1);

figure(2)

subplot(2,2,1),imshow(wu),title('原始图像')

subplot(2,2,2),imshow(result_rgb),title('RGB')

subplot(2,2,3),imhist(wu1);

subplot(2,2,4),imhist(result_rgb1);

figure(3)

subplot(2,2,1),imshow(wu),title('原始图像')

subplot(2,2,2),imshow(result_ycbcr),title('YCbCr')

subplot(2,2,3),imhist(wu1);

subplot(2,2,4),imhist(result_ycbcr1);

figure(4)

subplot(2,2,1),imshow(wu),title('原始图像')

subplot(2,2,2),imshow(W),title('灰度图像直方图均衡化')

subplot(2,2,3),imhist(wu1);

subplot(2,2,4),imhist(W);

2、antongd.m

%暗通道去雾算法

%===========================================================

%暗原色去雾算法是建立在户外自然场景暗通道优先法则的基础上的去雾方法

%其实就是解一个方程：I(x)=J(x)t(x)+A(1-t(x))

%其中I(x)是受到雾气污染的图像,J(x)是我们需要求的去雾后的图像

%t(x)是天空中云层的透射分布率，A是天空的亮度

%===========================================================

close all;

clear all;

img_name = imread('E:\Report\课程设计\图片\2.bmp'); %原始图像

I = double(img_name)/255;

[h,w,c] = size(I); %获取图像大小

w0 = 0.95; %去雾系数

dehaze = zeros(h,w,c); %初始化结果图像

win_dark = zeros(h,w); %初始化暗影通道图像

%--------------------------------------------------------------------------

%计算暗通道

%将三个通道中最暗的值赋给win_dark(i,j),使得三维图变成了二维图

win_dark = ordfilt2(win_dark,1,ones(9,9),'symmetric'); %9*9最小值滤波

%--------------------------------------------------------------------------

%计算大气亮度A

dark_channel = win_dark;

A = max(max(dark_channel));

[i0,j0] = find(dark_channel==A);

i = i0(1);

j = j0(1);

A = mean(I(i,j,:));

%--------------------------------------------------------------------------

%计算透射率t(x)

transmission = 1-w0*win_dark/A; %透射率预估

gray_I = I(:,:,3); %这里gray_I可以是RGB图像中任何一个通道

p = transmission; %透射率图

r = 50;

eps = 10^-3;

transmission_filter = guidedfilter(gray_I,p,r,eps);

%用guided filter对trasmission map做soft matting获得更为精细的透射率图

t0=0.1; %设置阈值

t1 = max(t0,transmission_filter);

%--------------------------------------------------------------------------

for i=1:c

for j=1:h

for l=1:w

dehaze(j,l,i)=(I(j,l,i)-A)/t1(j,l)+A; %无雾图像恢复

end

end

end

%--------------------------------------------------------------------------

%显示图像

figure

subplot(121)

imshow(win_dark);title('暗通道')

subplot(122)

imshow(t1);title('透射率图')

figure

subplot(121)

imshow(I);title('去雾前')

subplot(122)

imshow(dehaze);title('去雾后')

3、guidedfilter.m

%guided filter（导向滤波函数）

%===========================================================

%引导图：I（灰度图/单通道图像）

%输入图像：p（灰度图/单通道图像）

%本地窗口半径：r

%正规化参数：eps

%==========================================================================

function q = guidedfilter(I,p,r,eps)

%N=(2r+1)^2 except for boundary pixels.

mean_I = boxfilter(I,r)./N;

mean_p = boxfilter(p,r)./N;

mean_Ip = boxfilter(I.*p,r)./N;

cov_Ip = mean_Ip-mean_I.*mean_p; %协方差

mean_II = boxfilter(I.*I,r)./N;

var_I = mean_II-mean_I.*mean_I; %方差

a = cov_Ip./(var_I+eps);

b = mean_p-a.*mean_I;

mean_a = boxfilter(a,r)./N;

mean_b = boxfilter(b,r)./N;

q = mean_a.*I+mean_b;

end

4、boxfilter.m

%boxfilter（盒子滤波函数）

%===========================================================

%输入：imSrc

%给定的滑动窗口大小：r

%===========================================================

function imDst = boxfilter(imSrc,r)

[hei,wid] = size(imSrc);

imDst = zeros(size(imSrc)); %初始化

%--------------------------------------------------------------------------

imCum = cumsum(imSrc,1); %y轴累计求和，将每一行的数值一次累加到下一行

%y轴方向差异

imDst(1:r+1,:) = imCum(1+r:2*r+1,:);

%将imCum中的1+r到2*r+1行数据复制到imDst函数中的1到r+1行

imDst(r+2:hei-r,:) = imCum(2*r+2:hei,:)-imCum(1:hei-2*r-1,:);

%将imCum中的2*r+2到hei行数据依次减去1到hei-2*r-1行的值

%把最终结果存放到函数imDst的r+2到hei-r行

imDst(hei-r+1:hei,:) = repmat(imCum(hei,:),[r,1])-imCum(hei-2*r:hei-r-1,:);

%使用repmat函数将imCum变为r*1的格式

%然后减去函数imCum中从hei-2*r到hei-r-1行的数值

%得到的结果放入imDst的hei-r+1到hei行位置

%--------------------------------------------------------------------------

imCum = cumsum(imDst,2); %x轴累计求和，将每一列的数值一次累加到下一列

%y轴方向差异

imDst(:, 1:r+1) = imCum(:,1+r:2*r+1);

%将imCum中的1+r到2*r+1列数据复制到imDst函数中的1到r+1列

%把最终结果存放到函数imDst的r+2到hei-r列

imDst(:, wid-r+1:wid) = repmat(imCum(:,wid),[1,r])-imCum(:,wid-2*r:wid-r-1);

%使用repmat函数将imCum变为r*1的格式

%然后减去函数imCum中从hei-2*r到hei-r-1列的数值

%得到的结果放入imDst的hei-r+1到hei列位置

end