《熊大胡说|听王老师解读：数据是什么？》 中提到：
“数据的定义有强烈的时代特征……图像是一种重要的数据。但是在一百年前，我们认为图像不是数据。为什么？因为我们肉眼所见的这个美妙的世界，根本记录不下来。没有办法记录图像，怎么谈得上分析呢？但是今天，数码成像技术的成熟让所有的图像都能够以非常高的分辨率被记录下来，并进行分析，从而支撑很多有趣的应用。”
1、熊大图像处理
1.什么是图像？
图像类型多种多样，风景照、自拍照或者二维码等等，不同图像蕴含的信息也不同。简单的，像二维码，整个图像蕴含的就是一个ID信息；复杂的，像风景照、人物照，不同的人看到会有不同的理解。比如，图4-6熊大这张自拍照，在人类的帮助下，机器也许能认出这是一只熊，却距离理解图像仍有很长的道路，例如，这只熊是在真诚地笑还是假笑还是苦笑？
图1 生活中的图像
在对图像信息的提取过程中，图像处理至关重要。图像处理可以帮助我们提取、理解甚至改变图像信息。以修图为例，有以下几种操作：合照中把别人马赛克，留下自己——这是帮助我们突出理解图片信息；把自己拍丑的图片P美 —— 除去噪音；把自己和男神女神“P”张合影 ——这是添加图片信息。正因为大家对图片数据如此看重，因此都希望能够控制自己传达给别人的数据信息，因此图像信息处理手段自然而然成为大家的需求，这也是美颜修图软件层出不穷、自拍图像质量成为手机新卖点的原因。
图2 生活中的图像处理
过去Matlab往往是图像处理的主力，但其实R里也有着丰富的图像处理与分析的函数，下面我们来看看图片的常见处理操作。
2.图像的读入
图像处理的第一步自然是图像读入，R中对不同格式的图像读入都有专门的包，例如jpg、png、tiff、bmp 分别对应着jpeg、png、tiff、bmp等包，里面会有专门的读取函数，例如readJPEG 或readPNG等。一个简单的例子如下：.
图3 读入jpg格式图像
此处，有一个基本知识点：一幅图片可以看做由R（红色）、G（绿色）、B（蓝色）三原色通道叠加而成（具体内容下节会详细解释），需要注意的是，不同的包读出来数据结构也是不同的，例如熊大.jpg格式用jpeg包读出来是301*296*3的大小，只包含红绿蓝像素的信息，而熊大.png由png包读出来就是301*296*4的大小，还包含透明度的信息（即画图常用的alpha参数）。
图片数据往往会需要批量的图片处理，可以利用下面的函数：
除了读函数，每个包里也带有写函数，例如：writeJPEG，writePNG等。用法比较简单，感兴趣的读者可以自行探索。
R中的图像处理包众多，本节以图像处理包EBImage为例。与一般的在CRAN上的包不同，这个包托管在生物统计的bioconductor 上，因此包的安装也有所不同。
EBImage本身是为医学图像分析准备的，因此其本身的图像处理功能就十分强大。EBImage中，图片的读入是通过readImage实现的，目前支持三种图片格式：JPEG，TIFF和PNG；同时也能通过resize函数自由变化图像尺寸；通过display函数，可以将图片对象便利地绘制在R中的图形设备上即时检查，而不需要写出后再看。
3.图像的三基色和灰度图像
前面提到，图片实际上是由RGB三个通道混合而成的，每个通道都反映着各自的颜色信息，每个像素点中三种颜色的混合构成了图片中的各种颜色。每个通道的取值范围都是[0, 255]之间的整数，每个取值代表一个特定的颜色。在R中，会默认把数值转化成[0, 1]之间。如果读者对于这三个通道组合后是什么颜色很好奇，在R中有rgb这个函数帮我们来看不同强度的三通道会是什么颜色。当三个通道都为0时，就是#000000黑色；当三个通道都为1（R中的1相当于255），就是#FFFFFF白色。
也可以提取图片中的各个通道的颜色，只要将其它通道的数值置为0（如图4）。
图4 一张图片的单通道示意图
图5 将RGB图转换为灰度图
图像的灰度图可以由RGB彩色图像计算得到，计算公式有很多，例如取三个通道的均值或者加权均值等，感兴趣的读者可以在相关数字图像处理专著中学习。R中有现成的包和函数实现这项任务，像图5中所示的两种方法就包括RgoogleMaps里的RGB2GRAY函数以及EBImage包中的channel 函数。除此之外，还有imager包中的grayscale函数等。
4.马赛克背后的故事——分辨率
在买手机时，都会关注前置摄像头几千万像素、后置摄像头几千万像素，其实这里所涉及的就是图像分辨率。图像分辨率指图像中存储的信息量，与图像像素点有关，更本质地说，和图片背后的像素矩阵中的数据有关，而这也是传说中“马赛克”的基础。
马赛克可以通过如下简单方式实现：首先将原来301*296大小的图片（大约）划分为如图6所示的10*10块。
图6 划分图像
然后将每个格点中的像素都取平均值，得到模糊后的熊大图片——马赛克熊（见图7）。
图7 马赛克处理操作示意图
同样，也可以将其划分为20*20、50*50等规模，重复上述的代码，图像变得越来越清晰，即清晰度显著地提升，图像信息也越来越丰富。这个例子反映了图像矩阵中的数值和我们所看到的图像之间的关系。目前图像处理的一个非常热门的方向就是图片解码，Google一直在从事这个方向的工作，并且已经有了很大的进展 。通过这个技术，就能把以前像素低的老照片、珍贵的视频的分辨率提升。
图8 不同分辨率的图像
5.图像的加减乘除
因为图像是以矩阵或者数组的形式存储的，对其进行加减乘除，会有什么效果呢？请看下面的例子（见图9），这里使用display函数进行展示。.
图9 图像的加减
从图9可以看出，图像加上0.5，像素点数值增大，图像开始偏白；图像减去0.5，像素点数值减小，图像开始发黑。对图像进行乘除的效果类似（见图10）。
图10 1-图像和乘方处理图像
这两个方法可以有效地提取图像的信息，通过1-图像，实现反色，可以使外面白色背景三通道数值变为0（黑色），而里面深色变亮；而乘方操作，会使数值小的像素点变得更小最后趋于0，而本身为1的白像素点保持不变。由此可以预见，乘方足够大的话，图片中会只剩下白色和黑色。
6.图像间的合并和差分
同样我们能合并两个图像或者比较两个图像之间的差异。这里首先要保证，两幅图片的尺寸一样，即矩阵或数组的维度一致，因此首先用resize函数，将新读入的图片大小进行处理（见图11和图12）。
图11 图像的合并
图12 图像的差分
7.图像变换的魔术
在了解了图像处理的基本知识后，下面介绍一个R包：magick。这个R包里有一些有趣的函数。
首先读入图片（见图13），接下来我们可以使用函数对图片完成剪裁。
图13 图像的读入与裁剪
也可以尝试给图像加边框或者旋转图像（见图14）。
图14 图像的加框与旋转
还能实现图片的反色和加标签（见图15）。
图4 20 图像的反色和添加标签
当然还有很多其他的功能，如给图像添加噪点甚至转换为素描风格（见图16）。
图16 图片的加噪和素描转换
图片数据作为目前越来越火热的数据种类，其应用在未来有着广阔的前景。处理图片的方法也是数不胜数，在下节中，将介绍一个处理雾霾图像的实例。
2、 看图识雾霾
雾霾不断肆虐，走上街头越来越多的人开始带上“防毒面具”。雾霾的主要成分是PM2.5 。PM2.5能较长时间悬浮于空气中，其在空气中含量浓度越高，则空气污染越严重。根据中国环境监测总站发布的《2016年11月空气质量状况》，2016年11月，北京市空气质量优良天数比例仅43.3%；其余天数中，重度污染6天，严重污染1天，其首要污染物就是PM2.5。截至2016年11月底，全国共设1436个国家空气质量监测站点，由国家统一运行维护，监测数据直报国家并对外公开。但有监测站是否就意味着PM2.5的情况被严格管控了呢？事实上并没有，全国大量的监测站点数据质量其实难以保证。能否通过其他方式来判断雾霾的情况呢？比如有没有可能通过图片来判别雾霾情况呢？
前面介绍了图像操作的基本知识，本节将结合雾霾图片的实际案例，介绍如何从图片中提取我们想要的特征。下面将主要使用imager包中的函数来完成。
首先介绍一下数据。数据来自于2016年12月到2017年3月从熊大办公室拍摄的窗外照片，共206张照片，每张图片按拍摄时间命名。
在图17中，左右两边的图片分别是在空气优良和雾霾严重时拍摄的。这两张图片显然是比较容易分辨的，其标准可能是：左边天空蓝，层次清晰，右边天空灰，一团浆糊；左边视觉更饱满，右边让人昏昏欲睡……那么，从图片中，我们能提取出并量化这样的特征吗？让我们试试寻找这些变量。
图17 不同空气质量下的北京远景
1.计算灰度差分的方差
第一个解释性变量代表风景轮廓是否清晰，具体通过计算灰度图像差分的方差得到。通过对图17的观察可以发现，雾霾影响最明显的地方之一就是远处物体的边缘，正如图17中远方的树木的边缘和山的轮廓，天气晴朗时能清晰地分辨出这些景物，而雾霾时图像就变得特别模糊。为了提取这种轮廓特征，就需要先将灰度差分。前文已经介绍了灰度图片——灰度是指一张黑白图像中的颜色深度。对于彩色照片，灰度就是色彩（红绿蓝三色）的线性组合，可以认为，亮部灰度值高，暗部灰度值低。灰度的变化就可以反映出图像中明暗的变化，雾霾越严重，图像越不清晰，灰度变化也就越小。进一步思考，灰度变化最剧烈的地方其实正是这些物体边界的地方（同一个区域的物体（绿树），因为很相似，灰度变化不大；但是绿树和蓝天交界处会有很大的变化）。先对灰度差分下一个定义：某一像素点的灰度差分等于下格像素点灰度值减去上格像素点灰度值。这个差分就可以帮助提取这种边缘的轮廓信息（见图18）。
图18 灰度差分计算示意图
图19中的两幅图展示了差分后像素值的分布情况。
图19 像素差值分布图
非常直观地可以看出，好天气图像中由于细节轮廓多，像素点差分后的值所变动的范围就大（图19中左图框内区域）；而对于雾霾天气，由于其细节会变少，导致大部分像素差分后的值都聚集在0左右。那么如何反映这种区别呢？统计学告诉我们，变量的方差是衡量变动的一个指标，只需要求灰度差分的方差即可。求图像差分以及差分方差的函数实现方式如下所示（注意，这里用了函数包装）。我们将灰度差分的方差作为一个变量。例如图17中晴朗天气的灰度差分约为 ，而雾霾天气的灰度差分仅为 。
将不同的PM2.5等级对应灰度差分方差用箱线图展现出来（见图20）。从图20可以看出，随着PM2.5污染程度的增加，灰度差分的方差呈减小趋势。
图20 不同污染水平下的灰度差分分组箱线图
2.清晰度
第二个解释性变量是清晰度。从直觉上说，就是一张图片到底清不清楚。清晰度是如何计算的呢？这一变量来源于计算机视觉图像去雾领域的一篇论文名作 （下文简称论文），它表示反射光未衰减的比例。
图21 清晰度计算示意图
拍照的原理是物体的反射光被相机所接收，对于像素点x，设反射光为J(x)（不同景物有着不同的反射光），但它在传播过程中会有一定的衰减，我们将最后到达相机的剩余反射光的比例记为t(x)，即清晰度（其本质是透射系数）。最终到达相机的物体的反射光即为J(x)t(x)。除此之外，还有一部分光也会到达相机，这就是大气背景光，设其为A。到达的大气背景光占所有接收到的光的比例为(1-t(x))，即到达相机的大气背景光为A(1-t(x))。两者相加为相机收到的总光强，等于J(x)t(x)+A(1-t(x))。简单来说，由于雾霾天气对t(x)的影响较大，雾霾越大，t(x)越小（即景物反射光的比例越低，大气背景光比例越高）。因此在这里可以估计t(x)，作为其中一个重要解释变量。根据论文中的暗通道先验结果，t(x)可以通过下式估计（经过一定简化）：
这里I是已有的图片，A是大气背景光，Ω是基于中心像素点x的像素块（大小是31×31，但可以根据图片的大小进行调整）。理论上来说，图片在R，G，B（红、绿、蓝）3个不同的通道有3个不同的分量，这里用 （像素点z的色彩通道值）表示（C∈{R, G, B}）。下面的代码详细阐释了计算的过程：对于输入的图片M，维度为a×b，我们可以计算每个像素点在三个通道上的最小值，得到矩阵minMatrix。再对这个矩阵图像的每个像素点进行遍历操作，即在以每个小像素点(m, n)为中心的31×31小方格范围里（即公式中出现的 ）找到最小值，并将其保存下来，得到一个新的维数为a×b的矩阵。因此这个矩阵被称为图像的暗通道矩阵（当然这个矩阵也可以用灰度图像表示出来，如图22所示）。
图22 暗通道矩阵图像
将暗通道的最大值估计为大气光强A，并仿照上述方式得到t(x)的估计。在实际应用过程中，我们没有按照论文中那样选择t(x)的最小值作为整个图片的清晰度，经过调参，我们选取了70%分位数。同时，这里对A的估计做了一定的简化 。针对我们的展示图片，晴朗天气对应值为0.87，而雾霾天对应值为0.75。
同样地，箱线图结果显示该变量可以比较有效地区分不同空气质量状况（见图23）。
图23 不同污染水平下的清晰度分组箱线图
3.彩色图像信息熵
第三个解释性变量是信息熵。众所周知，熵是度量“变化性”的一种指标，当图片颜色变化比较明显时，熵值较大。因此，当天朗气清、惠风和畅时，我们能够看到蓝天白云、红花绿草；而当雾霾严重时，图像中的多数信息被雾霾所掩盖，因此残存的信息较少，我们也就只能看到“一桶浆糊”了。基于此，我们引入了信息熵这一变量，雾霾越严重，图像中的信息越少，则信息熵就越小。直观上，空气污染严重时，雾霾掩盖了图像中的部分信息。图像看上去灰蒙蒙一片，难以辨认具体细节。这里用信息熵来量化图像的平均信息量。
另一方面，我们看到的图片是由R（Red，红色）、G（Green，绿色）以及B（Blue，蓝色）三个色彩通道组成的。虽然一般对于彩色图片，通常的做法是处理成灰度图片；但是，对于“看图识雾霾”而言，有没有哪个颜色通道是更加重要的呢？为了保留更多的彩色图像信息，首先将预处理图片分到R，G，B三个通道下，并分别进行测试。结果符合我们的直觉，发现B（蓝色）空间下效果较好（毕竟，没有雾霾时蓝天白云中的蓝色成分是很重要的），因此这里以B空间下的图像为例进行分析。
信息熵应该如何计算呢？首先需要对图像像素进行过滤。如果一个像素点跟周围像素取值相似，则将这样的像素点灰度设置成0。为达到这一目的，首先需要计算差分因子。设(i,j)是N×M图像中的任意一点，f(i,j)为其蓝色像素值。在以其为中心的5×5窗口中统计像素值在对角线方向的差值的绝对值，称为差分因子D(i), i=1,2,3,4 。
图24 差分因子计算示意图
设阈值为τ，当D(3)< τ 且D(4)< τ 时，重置中心点的灰度f(i,j)=0。否则，不改变中心点灰度。这使得图像中只保留跟周围元素差异较大的像素点。接下来，遍历图像中的所有可行的像素点，利用如下公式求得信息熵：
其中， 是蓝色通道值为i的像素出现的频率。以空气质量状况迥异的两张图片为例，第一张图像的空气质量状况为优，PM2.5 = 7，其信息熵为1.624；第二张图像的空气质量等级达到了较严重污染，PM2.5 = 436，其信息熵为0.820，明显小于第一张照片。下面代码就是这一变量的具体实现。
同时，箱线图也佐证了信息熵对于雾霾具有良好的区分效果（见图25）。
图25 不同污染水平下的信息熵分组箱线图
4.饱和度
第四个解释性变量是饱和度。饱和度指图像中色彩的鲜艳程度，例如图4-31中的这束花的照片，调高它的饱和度，可以看到花朵确实变得娇艳欲滴了。简要而言，饱和度与图像中灰色成分所占的比例有关。雾霾会使图像变灰，色彩鲜艳程度减弱；雾霾越严重，图像饱和度越低。
图26 图片的饱和度
为了获得图片的饱和度，可以使用下述代码，这里主要靠的是imager包中的RGBtoHSV函数，这个函数可以将图像从RGB空间转化为HSV空间。在RGB空间中，每个颜色的度量是RGB通道的数值。HSV空间是对图像的另外一种描述，由色调（H）、饱和度（S）以及明亮程度（V）来表示。
提取其中的S（饱和度，Saturation）空间就能得到图像的饱和度矩阵（对图像的每一个像素，都有一个饱和度的度量）。将这个矩阵转化为向量并取其99.9%分位数和平均值作为描述饱和度的特征。从箱线图27可以看出，天气情况较好时，图片饱和度明显较高。
图27 不同污染水平下的饱和度分组箱线图
5.高频含量
第五个解释性变量是高频含量。如果把图片看成一个二维信号，将这个信号在频域上进行分解，可以得到图像频率分布。
图28 图片的频率分布
这种频率分布其实反映了图像的像素灰度在空间中变化的情况。例如，一面墙壁的图像，由于灰度值分布平坦，其低频成分就较强，而高频成分较弱。换句话说，高频成分决定了图像细节部分。因此，可以很直观地想到，天气晴朗时图像细节非常多，因此高频成分较多；而污染严重时，大多数图像细节被雾霾掩盖，高频成分会大幅减少。此处定义高频含量为频率分布的99.8%分位数减去频率分布的0.2%分位数。
对于所展示的两张图片，傅里叶变换后频率分布分别展示如图29所示。
图29 晴朗与雾霾天气频率分布对比图
如上过程可以通过以上代码实现。先将图像转换为灰度图像，然后通过傅里叶变换函数将图像进行变换，取变换后结果的实数部分，即得到频率分布。这里取频率分布的99.8%分位数和0.2%分位数之差代表高频含量。
图30 不同污染水平下的高频含量分组箱线图
在上文所述的例图中，当PM2.5 = 7时，天气晴朗，细节丰富，可以分辨远山、建筑、天空，此时高频含量为821.20；而PM2.5 = 436时，雾霾很严重，我们只能看清图像中物体的轮廓，难以分辨更多的细节，高频含量为536.97。从箱线图的结果来看，基本满足污染程度越高，则高频含量越小的特点，这表明这个变量具有解释效力。
6.批量处理图像
在上文提到的5种特征的基础上，我们将对图像逐个分析，提取变量，这里会用到前文介绍过的对文件夹中图像做批量处理的知识，用到的文件操作如下所示。
在获得了图像的各个特征后，就能将他们与特定的Y整理到一起，进行下一步的统计分析了。希望进一步了解图像处理内容的读者，可进入狗熊会公众号，输入关键词“雾霾图像分析”阅读精品案例。
联系客服
