zoukankan      html  css  js  c++  java
  • Python图像处理库(2)

    1.4 SciPy

    SciPyhttp://scipy.org/) 是建立在 NumPy 基础上,用于数值运算的开源工具包。SciPy 提供很多高效的操作,可以实现数值积分、优化、统计、信号处理,以及对我们来说最重要的图像处理功能。接下来,本节会介绍 SciPy 中大量有用的模块。SciPy 是个开源工具包,可以从http://scipy.org/Download 下载。

    1.4.1 图像模糊

    图像的高斯模糊是非常经典的图像卷积例子。本质上,图像模糊就是将(灰度)图像 I 和一个高斯核进行卷积操作:

    Iσ = I*

    其中 * 表示卷积操作; 是标准差为 σ 的二维高斯核,定义为 :

    G_sigma=frac{1}{2pisigma}e^{-(x^2+y^2)/2sigma^2}

    高斯模糊通常是其他图像处理操作的一部分,比如图像插值操作、兴趣点计算以及很多其他应用。

    SciPy 有用来做滤波操作的 scipy.ndimage.filters 模块。该模块使用快速一维分离的方式来计算卷积。你可以像下面这样来使用它:

    from PIL import Image
    from numpy import *
    from scipy.ndimage import filters
    
    im = array(Image.open('empire.jpg').convert('L'))
    im2 = filters.gaussian_filter(im,5)
    
    

    上面 guassian_filter() 函数的最后一个参数表示标准差。

    图 1-9 显示了随着 σ 的增加,一幅图像被模糊的程度。σ 越大,处理后的图像细节丢失越多。如果打算模糊一幅彩色图像,只需简单地对每一个颜色通道进行高斯模糊:

    im = array(Image.open('empire.jpg'))
    im2 = zeros(im.shape)
    for i in range(3):
      im2[:,:,i] = filters.gaussian_filter(im[:,:,i],5)
    im2 = uint8(im2)
    
    

    在上面的脚本中,最后并不总是需要将图像转换成 uint8 格式,这里只是将像素值用八位来表示。我们也可以使用:

    im2 = array(im2,'uint8')
    
    

    来完成转换。

    关于该模块更多的内容以及不同参数的选择,请查看http://docs.scipy.org/doc/scipy/reference/ndimage.html 上 SciPy 文档中的 scipy.ndimage部分。

     

    图 1-9:使用 scipy.ndimage.filters 模块进行高斯模糊:(a)原始灰度图像;(b)使用 σ=2 的高斯滤波器;(c)使用 σ=5 的高斯滤波器;(d)使用 σ=10 的高斯滤波器

    1.4.2 图像导数

    整本书中可以看到,在很多应用中图像强度的变化情况是非常重要的信息。强度的变化可以用灰度图像 I(对于彩色图像,通常对每个颜色通道分别计算导数)的 x 和 y 方向导数 Ix 和 Iy 进行描述。

    图像的梯度向量为∇I = [IxIy]T。梯度有两个重要的属性,一是梯度的大小

    left|oldsymbol{
abla I}
ight|=sqrt{{oldsymbol{I}_x}^2+{oldsymbol{I}_y}^2}

    它描述了图像强度变化的强弱,一是梯度的角度

    α=arctan2(IyIx)

    描述了图像中在每个点(像素)上强度变化最大的方向。NumPy 中的 arctan2() 函数返回弧度表示的有符号角度,角度的变化区间为 -π...π。

    我们可以用离散近似的方式来计算图像的导数。图像导数大多数可以通过卷积简单地实现:

    Ix=I*Dx 和 Iy=I*Dy

    对于 Dx 和 Dy,通常选择 Prewitt 滤波器:

    D_x=egin{vmatrix}-1&0&1\-1&0&1\ -1&0&1end{vmatrix} 和 D_y=egin{vmatrix}-1&-1&-1\0&0&0\ 1&1&1end{vmatrix}

    或者 Sobel 滤波器:

    D_x=egin{vmatrix}-1&0&1\-2&0&2\ -1&0&1end{vmatrix} 和 D_y=egin{vmatrix}-1&-2&-1\0&0&0\ 1&2&1end{vmatrix}

    这些导数滤波器可以使用 scipy.ndimage.filters 模块的标准卷积操作来简单地实现,例如:

    from PIL import Image
    from numpy import *
    from scipy.ndimage import filters
    
    im = array(Image.open('empire.jpg').convert('L'))
    
    # Sobel 导数滤波器
    imx = zeros(im.shape)
    filters.sobel(im,1,imx)
    
    imy = zeros(im.shape)
    filters.sobel(im,0,imy)
    
    magnitude = sqrt(imx**2+imy**2)
    
    

    上面的脚本使用 Sobel 滤波器来计算 x 和 y 的方向导数,以及梯度大小。sobel() 函数的第二个参数表示选择 x 或者 y 方向导数,第三个参数保存输出的变量。图 1-10 显示了用 Sobel 滤波器计算出的导数图像。在两个导数图像中,正导数显示为亮的像素,负导数显示为暗的像素。灰色区域表示导数的值接近于零。

     

    图 1-10:使用 Sobel 导数滤波器计算导数图像:(a)原始灰度图像;(b)x 导数图像;(c)y导数图像;(d)梯度大小图像

    上述计算图像导数的方法有一些缺陷:在该方法中,滤波器的尺度需要随着图像分辨率的变化而变化。为了在图像噪声方面更稳健,以及在任意尺度上计算导数,我们可以使用高斯导数滤波器:

    Ix=I*Gσx 和 Iy=I*Gσy

    其中,Gσx和 Gσy 表示  在 x 和 y 方向上的导数, 为标准差为 σ 的高斯函数。

    我们之前用于模糊的 filters.gaussian_filter() 函数可以接受额外的参数,用来计算高斯导数。可以简单地按照下面的方式来处理:

    sigma = 5 # 标准差
    
    imx = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    
    imy = zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    
    

    该函数的第三个参数指定对每个方向计算哪种类型的导数,第二个参数为使用的标准差。你可以查看相应文档了解详情。图 1-11 显示了不同尺度下的导数图像和梯度大小。你可以和图 1-9 中做相同尺度模糊的图像做比较。

     

    图 1-11:使用高斯导数计算图像导数:x 导数图像(上),y 导数图像(中),以及梯度大小图像(下);(a)为原始灰度图像,(b)为使用 σ=2 的高斯导数滤波器处理后的图像,(c)为使 用 σ=5 的高斯导数滤波器处理后的图像,(d)为使用 σ=10 的高斯导数滤波器处理后的图像

    1.4.3 形态学:对象计数

    形态学(或数学形态学)是度量和分析基本形状的图像处理方法的基本框架与集合。形态学通常用于处理二值图像,但是也能够用于灰度图像。二值图像是指图像的每个像素只能取两个值,通常是 0 和 1。二值图像通常是,在计算物体的数目,或者度量其大小时,对一幅图像进行阈值化后的结果。你可以从 http://en.wikipedia.org/wiki/Mathematical_morphology 大体了解形态学及其处理图像的方式。

    scipy.ndimage 中的 morphology 模块可以实现形态学操作。你可以使用 scipy.ndimage 中的measurements 模块来实现二值图像的计数和度量功能。下面通过一个简单的例子介绍如何使用它们。

    考虑在图 1-12a3 里的二值图像,计算该图像中的对象个数可以通过下面的脚本实现:

    3这个图像实际上是图像“分割”后的结果。如果你想知道该图像是如何创建的,可以查看 9.3 节。

    from scipy.ndimage import measurements,morphology
    
    # 载入图像,然后使用阈值化操作,以保证处理的图像为二值图像
    im = array(Image.open('houses.png').convert('L'))
    im = 1*(im<</span>128)
    
    labels, nbr_objects = measurements.label(im)
    print "Number of objects:", nbr_objects
    
    
    

    上面的脚本首先载入该图像,通过阈值化方式来确保该图像是二值图像。通过和 1 相乘,脚本将布尔数组转换成二进制表示。然后,我们使用 label() 函数寻找单个的物体,并且按照它们属于哪个对象将整数标签给像素赋值。图 1-12b 是 labels 数组的图像。图像的灰度值表示对象的标签。可以看到,在一些对象之间有一些小的连接。进行二进制开(binary open)操作,我们可以将其移除:

    # 形态学开操作更好地分离各个对象
    im_open = morphology.binary_opening(im,ones((9,5)),iterations=2)
    
    labels_open, nbr_objects_open = measurements.label(im_open)
    print "Number of objects:", nbr_objects_open
    
    

    binary_opening() 函数的第二个参数指定一个数组结构元素。该数组表示以一个像素为中心时,使用哪些相邻像素。在这种情况下,我们在 y 方向上使用 9 个像素(上面 4 个像素、像素本身、下面 4 个像素),在 x 方向上使用 5 个像素。你可以指定任意数组为结构元素,数组中的非零元素决定使用哪些相邻像素。参数 iterations 决定执行该操作的次数。你可以尝试使用不同的迭代次数 iterations 值,看一下对象的数目如何变化。你可以在图 1-12c 与图 1-12d 中查看经过开操作后的图像,以及相应的标签图像。正如你想象的一样,binary_closing() 函数实现相反的操作。我们将该函数和在 morphology 和 measurements 模块中的其他函数的用法留作练习。你可以从 scipy.ndimage 模块文档 http://docs.scipy.org/doc/scipy/reference/ndimage.html 中了解关于这些函数的更多知识。

     

    图 1-12:形态学示例。使用二值开操作将对象分开,然后计算物体的数目:(a)为原始二值图像;(b)为对应原始图像的标签图像,其中灰度值表示物体的标签;(c)为使用开操作后的二值图像;(d)为开操作后图像的标签图像

    1.4.4 一些有用的SciPy模块

    SciPy 中包含一些用于输入和输出的实用模块。下面介绍其中两个模块:io 和 misc

    1. 读写.mat文件

      如果你有一些数据,或者在网上下载到一些有趣的数据集,这些数据以 Matlab 的 .mat 文件格式存储,那么可以使用 scipy.io 模块进行读取。

      data = scipy.io.loadmat('test.mat')
      
      

      上面代码中,data 对象包含一个字典,字典中的键对应于保存在原始 .mat 文件中的变量名。由于这些变量是数组格式的,因此可以很方便地保存到 .mat 文件中。你仅需创建一个字典(其中要包含你想要保存的所有变量),然后使用 savemat() 函数:

      data = {}
      data['x'] = x
      scipy.io.savemat('test.mat',data)
      
      

      因为上面的脚本保存的是数组 x,所以当读入到 Matlab 中时,变量的名字仍为 x。关于scipy.io 模块的更多内容,请参见在线文档http://docs.scipy.org/doc/scipy/reference/io.html

    2. 以图像形式保存数组

      因为我们需要对图像进行操作,并且需要使用数组对象来做运算,所以将数组直接保存为图像文件 4 非常有用。本书中的很多图像都是这样的创建的。

      imsave() 函数可以从 scipy.misc 模块中载入。要将数组 im 保存到文件中,可以使用下面的命令:

      from scipy.misc import imsave
      imsave('test.jpg',im)
      
      

      scipy.misc 模块同样包含了著名的 Lena 测试图像:

      lena = scipy.misc.lena()
      
      

      该脚本返回一个 512×512 的灰度图像数组。

    4所有 Pylab 图均可保存为多种图像格式,方法是点击图像窗口中的“保存”按钮。

    1.5 高级示例:图像去噪

    我们通过一个非常实用的例子——图像的去噪——来结束本章。图像去噪是在去除图像噪声的同时,尽可能地保留图像细节和结构的处理技术。我们这里使用 ROF(Rudin-Osher-Fatemi)去噪模型。该模型最早出现在文献 [28] 中。图像去噪对于很多应用来说都非常重要;这些应用范围很广,小到让你的假期照片看起来更漂亮,大到提高卫星图像的质量。ROF 模型具有很好的性质:使处理后的图像更平滑,同时保持图像边缘和结构信息。

    ROF 模型的数学基础和处理技巧非常高深,不在本书讲述范围之内。在讲述如何基于 Chambolle 提出的算法 [5] 实现 ROF 求解器之前,本书首先简要介绍一下 ROF 模型。

    一幅(灰度)图像 I 的全变差(Total Variation,TV)定义为梯度范数之和。在连续表示的情况下,全变差表示为:

    J(oldsymbol{I})=intleft|
ablaoldsymbol{I}
ight|	ext{dx}            (1.1)

    在离散表示的情况下,全变差表示为:

    J(oldsymbol{I})=sum_{	ext{x}}left|
ablaoldsymbol{I}
ight|

    其中,上面的式子是在所有图像坐标 x=[x, y] 上取和。

    在 Chambolle 提出的 ROF 模型里,目标函数为寻找降噪后的图像 U,使下式最小:

    min_Uleft|left|oldsymbol{I}-oldsymbol{U}
ight|
ight|^2+2lambda J(oldsymbol{U}),

    其中范数 ||I-U|| 是去噪后图像 U 和原始图像 I 差异的度量。也就是说,本质上该模型使去噪后的图像像素值“平坦”变化,但是在图像区域的边缘上,允许去噪后的图像像素值“跳跃”变化。

    按照论文 [5] 中的算法,我们可以按照下面的代码实现 ROF 模型去噪:

    from numpy import *
    
    def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):
      """ 使用A. Chambolle(2005)在公式(11)中的计算步骤实现Rudin-Osher-Fatemi(ROF)去噪模型
    
        输入:含有噪声的输入图像(灰度图像)、U 的初始值、TV 正则项权值、步长、停业条件
    
        输出:去噪和去除纹理后的图像、纹理残留"""
    
    m,n = im.shape # 噪声图像的大小
    
    # 初始化
    U = U_init
    Px = im # 对偶域的x 分量
    Py = im # 对偶域的y 分量
    error = 1
    
    while (error > tolerance):
      Uold = U
    
      # 原始变量的梯度
      GradUx = roll(U,-1,axis=1)-U # 变量U 梯度的x 分量
      GradUy = roll(U,-1,axis=0)-U # 变量U 梯度的y 分量
    
      # 更新对偶变量
      PxNew = Px + (tau/tv_weight)*GradUx
      PyNew = Py + (tau/tv_weight)*GradUy
      NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))
    
      Px = PxNew/NormNew # 更新x 分量(对偶)
      Py = PyNew/NormNew # 更新y 分量(对偶)
    
      # 更新原始变量
      RxPx = roll(Px,1,axis=1) # 对x 分量进行向右x 轴平移
      RyPy = roll(Py,1,axis=0) # 对y 分量进行向右y 轴平移
    
      DivP = (Px-RxPx)+(Py-RyPy) # 对偶域的散度
      U = im + tv_weight*DivP # 更新原始变量
    
      # 更新误差
      error = linalg.norm(U-Uold)/sqrt(n*m);
    
    return U,im-U # 去噪后的图像和纹理残余
    
    

    在这个例子中,我们使用了 roll() 函数。顾名思义,在一个坐标轴上,它循环“滚动”数组中的元素值。该函数可以非常方便地计算邻域元素的差异,比如这里的导数。我们还使用了linalg.norm() 函数,该函数可以衡量两个数组间(这个例子中是指图像矩阵 U和 Uold)的差异。我们将这个 denoise() 函数保存到 rof.py 文件中。

    下面使用一个合成的噪声图像示例来说明如何使用该函数:

    from numpy import *
    from numpy import random
    from scipy.ndimage import filters
    import rof
    
    # 使用噪声创建合成图像
    im = zeros((500,500))
    im[100:400,100:400] = 128
    im[200:300,200:300] = 255
    im = im + 30*random.standard_normal((500,500))
    
    U,T = rof.denoise(im,im)
    G = filters.gaussian_filter(im,10)
    
    # 保存生成结果
    from scipy.misc import imsave
    imsave('synth_rof.pdf',U)
    imsave('synth_gaussian.pdf',G)
    
    

    原始图像和图像的去噪结果如图 1-13 所示。正如你所看到的,ROF 算法去噪后的图像很好地保留了图像的边缘信息。

     

    图 1-13:使用 ROF 模型对合成图像去噪:(a)为原始噪声图像;(b)为经过高斯模糊的图像(σ=10);(c)为经过 ROF 模型去噪后的图像

    下面看一下在实际图像中使用 ROF 模型去噪的效果:

    from PIL import Image
    from pylab import *
    import rof
    
    im = array(Image.open('empire.jpg').convert('L'))
    U,T = rof.denoise(im,im)
    
    figure()
    gray()
    imshow(U)
    axis('equal')
    axis('off')
    show()
    
    

    经过 ROF 去噪后的图像如图 1-14c 所示。为了方便比较,该图中同样显示了模糊后的图像。可以看到,ROF 去噪后的图像保留了边缘和图像的结构信息,同时模糊了“噪声”。

     

    图 1-14:使用 ROF 模型对灰度图像去噪:(a)为原始噪声图像;(b)为经过高斯模糊的图像(σ=5);(c)为经过 ROF 模型去噪后的图像

  • 相关阅读:
    jQuery的Autocomplete插件的远程url取json数据的问题
    DataTable 转换成 Json的3种方法
    IFRAME练习 各种调用
    <%: Html.DisplayNameFor与 <%: Html.DisplayFor的区别
    MVC中@Html.DisPlayFor(model=>model.newsName)和 @Model.newsName的区别
    编程思想:如何设计一个好的通信网络协议
    编程思想:巧用位运算重构代码
    编程思想:小谈网关项目中的设计模式
    Unity应用架构设计(13)——日志组件的实施
    Unity应用架构设计(12)——AOP思想的实践
  • 原文地址:https://www.cnblogs.com/jankerxp/p/7774026.html
Copyright © 2011-2022 走看看