用python写图像处理

2018-03-05 21:57:38 jiaoyangwm 阅读数 39470

第1章 基本的图像操作和处理

1.1 PIL:Python图像处理类库

PIL(Python Imaging Library,图像处理库)提供了通用的图像处理功能,以及大量有用的基本图像操作。PIL库已经集成在Anaconda库中,推荐使用Anaconda,简单方便,常用库都已经集成。

PIL简明教程

  • 读入一副图像:
from PIL import Image
from pylab import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
figure()

pil_im = Image.open('E:\python\Python Computer Vision\Image data\empire.jpg')
gray()
subplot(121)
title(u'原图',fontproperties=font)
axis('off')
imshow(pil_im)

pil_im = Image.open('E:\python\Python Computer Vision\Image data\empire.jpg').convert('L')
subplot(122)
title(u'灰度图',fontproperties=font)
axis('off')
imshow(pil_im)

show()

这里写图片描述

1.1.1 转换图像格式——save()函数

from PCV.tools.imtools import get_imlist #导入原书的PCV模块
from PIL import Image
import os
import pickle

filelist = get_imlist('E:/python/Python Computer Vision/test jpg/') #获取convert_images_format_test文件夹下的图片文件名(包括后缀名)
imlist = open('E:/python/Python Computer Vision/test jpg/imlist.txt','wb+')
#将获取的图片文件列表保存到imlist.txt中
pickle.dump(filelist,imlist) #序列化
imlist.close()

for infile in filelist:
    outfile = os.path.splitext(infile)[0] + ".png" #分离文件名与扩展名
    if infile != outfile:
        try:
            Image.open(infile).save(outfile)
        except IOError:
            print ("cannot convert", infile)

其中,test jpg文件夹是作者自己建立的文件夹,存放测试的**.jpg图像,源代码证添加了部分代码以便将获取的图像文件名保存下来,同时将所有的图像转化为.png格式,运行程序后的结果如下:
这里写图片描述

这里写图片描述

PIL中的open()函数用于创建PIL图像对象,sace()方法用于保存如下到指定文件名的文件夹,上述过程将后缀变为.png,但文件名不变

1.1.2 创建缩略图

利用PIL可以很容易的创建缩略图,设置缩略图的大小,并用元组保存起来,调用thumnail()方法即可生成缩略图。创建缩略图的代码见下面。
例如创建最长边为128像素的缩略图,可以使用:

pil_im.thumbnail((128,128))

1.1.3 复制并粘贴图像区域

调用crop()方法即可从一幅图像中进行区域拷贝,拷贝出区域后,可以对区域进行旋转等变换。

box=(100,100,400,400)
region=pil_im.crop(box)

目标区域由四元组来指定,坐标依次为(左,上,右,下),PIL中指定坐标系的左上角坐标为(0,0),可以旋转后利用paste()放回去,具体实现如下:

region=region.transpose(Image.ROTATE_180)
pil_im.paste(region,box)

1.1.4 调整尺寸和旋转

  • 调整尺寸:利用resize()方法,参数是一个元组,用来指定新图像的大小:
out=pil_im.resize((128,128))
  • 旋转:利用rotate()方法,逆时针方式表示角度
out=pil_im.rotate(45)

上述操作的代码如下:

from PIL import Image
from pylab import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
figure()

# 显示原图
pil_im = Image.open('E:/python/Python Computer Vision/Image data/empire.jpg')
print(pil_im.mode, pil_im.size, pil_im.format)
subplot(231)
title(u'原图', fontproperties=font)
axis('off')
imshow(pil_im)

# 显示灰度图
pil_im = Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L')
gray()
subplot(232)
title(u'灰度图', fontproperties=font)
axis('off')
imshow(pil_im)

# 复制并粘贴区域
pil_im = Image.open('E:/python/Python Computer Vision/Image data/empire.jpg')
box = (100, 100, 400, 400)
region = pil_im.crop(box)
region = region.transpose(Image.ROTATE_180)
pil_im.paste(region, box)
subplot(233)
title(u'复制粘贴区域', fontproperties=font)
axis('off')
imshow(pil_im)

# 缩略图
pil_im = Image.open('E:/python/Python Computer Vision/Image data/empire.jpg')
size = 128, 128
pil_im.thumbnail(size)
print(pil_im.size)
subplot(234)
title(u'缩略图', fontproperties=font)
axis('off')
imshow(pil_im)
pil_im.save('E:/python/Python Computer Vision/Image data/empire thumbnail.jpg')# 保存缩略图

#调整图像尺寸
pil_im=Image.open('E:/python/Python Computer Vision/Image data/empire thumbnail.jpg')
pil_im=pil_im.resize(size)
print(pil_im.size)
subplot(235)
title(u'调整尺寸后的图像',fontproperties=font)
axis('off')
imshow(pil_im)

#旋转图像45°
pil_im=Image.open('E:/python/Python Computer Vision/Image data/empire thumbnail.jpg')
pil_im=pil_im.rotate(45)
subplot(236)
title(u'旋转45°后的图像',fontproperties=font)
axis('off')
imshow(pil_im)

show()

运行结果如下:
这里写图片描述

1.2 Matplotlib库

当在处理数学及绘图或在图像上描点、画直线、曲线时,Matplotlib是一个很好的绘图库,它比PIL库提供了更有力的特性。

matplotlib教程

1.2.1 画图、描点和线

from PIL import Image
from pylab import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

# 读取图像到数组中
im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg'))
figure()

# 绘制有坐标轴的
subplot(121)
imshow(im)
x = [100, 100, 400, 400]
y = [200, 500, 200, 500]

# 使用红色星状标记绘制点
plot(x, y, 'r*')


# 绘制连接两个点的线(默认为蓝色)
plot(x[:2], y[:2])
title(u'绘制empire.jpg', fontproperties=font)

# 不显示坐标轴的
subplot(122)
imshow(im)
x = [100, 100, 400, 400]
y = [200, 500, 200, 500]

plot(x, y, 'r*')
plot(x[:2], y[:2])
axis('off')
title(u'绘制empire.jpg', fontproperties=font)

show()
# show()命令首先打开图形用户界面(GUI),然后新建一个窗口,该图形用户界面会循环阻断脚本,然后暂停,
# 直到最后一个图像窗口关闭。每个脚本里,只能调用一次show()命令,通常相似脚本的结尾调用。

这里写图片描述

绘图时还有很多可选的颜色和样式,如表1-1,1-2,1-3所示,应用例程如下:

plot(x,y)          #默认为蓝色实线
plot(x,y,'go-')    #带有圆圈标记的绿线
plot(x,y,'ks:')    #带有正方形标记的黑色虚线

表1-1 用PyLab库绘图的基本颜色格式命令

符号 颜色
‘b’ 蓝色
‘g’ 绿色
‘r’ 红色
‘c’ 青色
‘m’ 品红
‘y’ 黄色
‘k’ 黑色
‘w’ 白色

表1-2 用PyLab库绘图的基本线型格式命令

符号 线型
‘-‘ 实线
‘–’ 虚线
‘:’ 点线

表1-3 用PyLab库绘图的基本绘制标记格式命令

符号 标记
‘.’
‘o’ 圆圈
’s’ 正方形
‘*’ 星型
‘+’ 加号
‘*’ 叉号

1.2.2 图像轮廓和直方图

from PIL import Image
from pylab import *

# 添加中文字体支持
from matplotlib.font_manager import FontProperties

font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)
# 打开图像,并转成灰度图像
im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))

# 新建一个图像
figure()
subplot(121)
# 不使用颜色信息
gray()
# 在原点的左上角显示轮廓图像
contour(im, origin='image')
axis('equal')
axis('off')
title(u'图像轮廓图', fontproperties=font)

subplot(122)
# 利用hist来绘制直方图
# 第一个参数为一个一维数组
# 因为hist只接受一维数组作为输入,所以要用flatten()方法将任意数组按照行优先准则转化成一个一维数组
# 第二个参数指定bin的个数
hist(im.flatten(), 128)
title(u'图像直方图', fontproperties=font)
# plt.xlim([0,250])
# plt.ylim([0,12000])

show()

这里写图片描述

1.2.3 交互式标注

有时候用户需要和应用进行交互,比如在图像中用点做标识,或者在一些训练数据中进行注释,PyLab提供了一个很简介好用的函数gitput()来实现交互式标注。

from PIL import Image
from pylab import *

im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg'))
imshow(im)

print('Please click 3 points')
x = ginput(3)
print('you clicked:', x)
show()

输出:

you clicked: 
[(118.4632306896458, 177.58271393177051), 
(118.4632306896458, 177.58271393177051),
(118.4632306896458, 177.58271393177051)]

上面代码先读取empire.jpg图像,显示读取的图像,然后用ginput()交互注释,这里设置的交互注释数据点设置为3个,用户在注释后,会将注释点的坐标打印出来。

1.3 NumPy库

NumPy在线文档
NumPy是Python一个流行的用于科学计算包。它包含了很多诸如矢量、矩阵、图像等其他非常有用的对象和线性代数函数。

1.3.1 图像数组表示

在前面图像的示例中,我们将图像用array()函数转为NumPy数组对象,但是并没有提到它表示的含义。数组就像列表一样,只不过它规定了数组中的所有元素必须是相同的类型,除非指定以外,否则数据类型灰按照数据类型自动确定。
举例如下:

from PIL import Image
from pylab import *

im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg'))
print (im.shape, im.dtype)
im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'),'f')
print (im.shape, im.dtype)

输出:

(800, 569, 3) uint8
(800, 569) float32

解释:

第一个元组表示图像数组大小(行、列、颜色通道)
第二个字符串表示数组元素的数据类型,因为图像通常被编码为8位无符号整型;
1. uint8:默认类型
2. float32:对图像进行灰度化,并添加了'f'参数,所以变为浮点型
  • 数组元素如何访问——使用下标访问
value=im[i,j,k]
  • 多个数组元素如何发给我——使用数组切片方式访问,返回的是以指定间隔下标访问该数组的元素值
im[i,:] = im[j,:]     #将第j行的数值赋值给第i行
im[:,j] = 100         #将第i列所有数值设为100
im[:100,:50].sum()    #计算前100行、前50列所有数值的和
im[50:100,50:100]     #50~100行,50~100列,不包含第100行和100列
im[i].mean()          #第i行所有数值的平均值
im[:,-1]              #最后一列
im[-2,:]/im[-2]       #倒数第二行

1.3.2 灰度变换

将图像读入NumPy数组对象后,我们可以对它们执行任意数学操作,一个简单的例子就是图像的灰度变换,考虑任意函数f,它将0~255映射到自身,也就是输出区间和输入区间相同。
举例如下:

from PIL import Image
from numpy import *
from pylab import *

im=array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))
print(int(im.min()),int(im.max()))

im2=255-im               #对图像进行反向处理
print(int(im2.min()),int(im2.max())) #查看最大/最小元素

im3=(100.0/255)*im+100   #将图像像素值变换到100...200区间
print(int(im3.min()),int(im3.max()))

im4=255.0*(im/255.0)**2  #对像素值求平方后得到的图像
print(int(im4.min()),int(im4.max()))

figure()
gray()
subplot(131)
imshow(im2)
axis('off')
title(r'$f(x)=255-x$')

subplot(132)
imshow(im3)
axis('off')
title(r'$f(x)=\frac{100}{255}x+100$')

subplot(133)
imshow(im4)
axis('off')
title(r'$f(x)=255(\frac{x}{255})^2$')

show()

输出:

3 255
0 252
101 200
0 255

这里写图片描述

  • array变换的相反操作可以利用PIL的fromarray()函数来完成
pil_im=Image.fromarray(im)
  • 如果之前的操作将”uint8”数据类型转化为其他类型,则在创建PIL图像之前,需要将数据类型转换回来:
pil_im=Image.fromarray(uint8(im))

1.3.3 图像缩放

NumPy数组将成为我们对图像及数据进行处理的最主要工具,但是调整矩阵大小并没有一种简单的方法。我们可以用PIL图像对象转换写一个简单的图像尺寸调整函数:

def imresize(im,sz):
    """    Resize an image array using PIL. """
    pil_im = Image.fromarray(uint8(im))

    return array(pil_im.resize(sz))

上面定义的调整函数,在imtools.py中你可以找到它。

1.3.4 直方图均衡化

直方图均衡化指将一幅图像的灰度直方图变平,使得变换后的图像中每个灰度值的分布概率都相同,该方法是对灰度值归一化的很好的方法,并且可以增强图像的对比度。

  • 变换函数:图像中像素值的累积分布函数(cdf),将像素值的范围映射到目标范围的归一化操作

下面的函数是直方图均衡化的具体实现:

def histeq(im,nbr_bins=256):
  """ 对一幅灰度图像进行直方图均衡化"""

  # 计算图像的直方图
  imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
  cdf = imhist.cumsum()      # 累积分布函数
  cdf = 255 * cdf / cdf[-1]  # 归一化
  # 此处使用到累积分布函数cdf的最后一个元素(下标为-1),其目的是将其归一化到0~1范围

  # 使用累积分布函数的线性插值,计算新的像素值
  im2 = interp(im.flatten(),bins[:-1],cdf)

  return im2.reshape(im.shape), cdf

解释:

  1. 该函数有两个参数

    • 灰度图像
    • 直方图中使用的bin的数目
  2. 函数返回值

    • 均衡化后的图像
    • 用来做像素值映射的累积分布函数

程序实现:

from PIL import Image
from pylab import *
from PCV.tools import imtools

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))
# 打开图像,并转成灰度图像
#im = array(Image.open('../data/AquaTermi_lowcontrast.JPG').convert('L'))
im2, cdf = imtools.histeq(im)

figure()
subplot(2, 2, 1)
axis('off')
gray()
title(u'原始图像', fontproperties=font)
imshow(im)

subplot(2, 2, 2)
axis('off')
title(u'直方图均衡化后的图像', fontproperties=font)
imshow(im2)

subplot(2, 2, 3)
axis('off')
title(u'原始直方图', fontproperties=font)
#hist(im.flatten(), 128, cumulative=True, normed=True)
hist(im.flatten(), 128, normed=True)

subplot(2, 2, 4)
axis('off')
title(u'均衡化后的直方图', fontproperties=font)
#hist(im2.flatten(), 128, cumulative=True, normed=True)
hist(im2.flatten(), 128, normed=True)

show()

结果:
这里写图片描述

这里写图片描述

1.3.5 图像平均

对图像取平均是一种图像降噪的简单方法,经常用于产生艺术效果。假设所有的图像具有相同的尺寸,我们可以对图像相同位置的像素相加取平均,下面是一个演示对图像取平均的例子:

def compute_average(imlist):
  """ 计算图像列表的平均图像"""

  # 打开第一幅图像,将其存储在浮点型数组中
  averageim = array(Image.open(imlist[0]), 'f')

  for imname in imlist[1:]:
    try:
      averageim += array(Image.open(imname))
    except:
      print imname + '...skipped'
  averageim /= len(imlist)

  # 返回uint8 类型的平均图像
  return array(averageim, 'uint8')

这里写图片描述
注意:有可能因为某些图像打不开而导致平均的结果只是某一幅自身或某两幅图像的平均

1.3.5 对图像进行主成分分析

PCA(Principal Component Analysis,主成分分析)是一个非常有用的降维技巧。它可以在使用尽可能少维数的前提下,尽量多地保持训练数据的信息,在此意义上是一个最佳技巧。即使是一幅 100×100 像素的小灰度图像,也有 10 000 维,可以看成 10 000 维空间中的一个点。一兆像素的图像具有百万维。由于图像具有很高的维数,在许多计算机视觉应用中,我们经常使用降维操作。PCA 产生的投影矩阵可以被视为将原始坐标变换到现有的坐标系,坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行 PCA 变换,图像需要转换成一维向量表示。我们可以使用 NumPy 类库中的flatten() 方法进行变换。

将变平的图像堆积起来,我们可以得到一个矩阵,矩阵的一行表示一幅图像。在计算主方向之前,所有的行图像按照平均图像进行了中心化。我们通常使用 SVD(Singular Value Decomposition,奇异值分解)方法来计算主成分;但当矩阵的维数很大时,SVD 的计算非常慢,所以此时通常不使用 SVD 分解。

下面就是 PCA 操作的代码:

from PIL import Image
from numpy import *

def pca(X):
  """ 主成分分析:
    输入:矩阵X ,其中该矩阵中存储训练数据,每一行为一条训练数据
    返回:投影矩阵(按照维度的重要性排序)、方差和均值"""

  # 获取维数
  num_data,dim = X.shape

  # 数据中心化
  mean_X = X.mean(axis=0)
  X = X - mean_X

if dim>num_data:
  # PCA- 使用紧致技巧
  M = dot(X,X.T) # 协方差矩阵
  e,EV = linalg.eigh(M) # 特征值和特征向量
  tmp = dot(X.T,EV).T # 这就是紧致技巧
  V = tmp[::-1] # 由于最后的特征向量是我们所需要的,所以需要将其逆转
  S = sqrt(e)[::-1] # 由于特征值是按照递增顺序排列的,所以需要将其逆转
  for i in range(V.shape[1]):
    V[:,i] /= S
else:
  # PCA- 使用SVD 方法
  U,S,V = linalg.svd(X)
  V = V[:num_data] # 仅仅返回前nun_data 维的数据才合理

# 返回投影矩阵、方差和均值
return V,S,mean_X

该函数首先通过减去每一维的均值将数据中心化,然后计算协方差矩阵对应最大特征值的特征向量,此时可以使用简明的技巧或者 SVD 分解。这里我们使用了 range() 函数,该函数的输入参数为一个整数 n,函数返回整数 0…(n-1) 的一个列表。你也可以使用 arange() 函数来返回一个数组,或者使用 xrange() 函数返回一个产生器(可能会提升速度)。我们在本书中贯穿使用range() 函数。

如果数据个数小于向量的维数,我们不用 SVD 分解,而是计算维数更小的协方差矩阵 XXT 的特征向量。通过仅计算对应前 k(k 是降维后的维数)最大特征值的特征向量,可以使上面的 PCA 操作更快。由于篇幅所限,有兴趣的读者可以自行探索。矩阵 V 的每行向量都是正交的,并且包含了训练数据方差依次减少的坐标方向。

我们接下来对字体图像进行 PCA 变换。fontimages.zip 文件包含采用不同字体的字符 a 的缩略图。所有的 2359 种字体可以免费下载 2。假定这些图像的名称保存在列表 imlist 中,跟之前的代码一起保存传在 pca.py 文件中,我们可以使用下面的脚本计算图像的主成分:

import pickle
from PIL import Image
from numpy import *
from pylab import *
from PCV.tools import imtools,pca

# Uses sparse pca codepath

# 获取图像列表和尺寸
imlist=imtools.get_imlist('E:/python/Python Computer Vision/Image data/fontimages/a_thumbs')
# open ont image to get the size
im=array(Image.open(imlist[0]))
# get the size of the images
m,n=im.shape[:2]
# get the number of images
imnbr=len(imlist)
print("The number of images is %d" % imnbr)

# create matrix to store all flattened images
immatrix = array([array(Image.open(imname)).flatten() for imname in imlist],'f')

# PCA降维
V,S,immean=pca.pca(immatrix)

# 保存均值和主成分
#f = open('../ch01/font_pca_modes.pkl', 'wb')
#pickle.dump(immean,f)
#pickle.dump(V,f)
#f.close()

# Show the images (mean and 7 first modes)
# This gives figure 1-8 (p15) in the book.

figure()
gray()
subplot(241)
axis('off')
imshow(immean.reshape(m,n))
for i in range(7):
    subplot(2,4,i+2)
    imshow(V[i].reshape(m,n))
    axis('off')

show()

注意,这些图像在拉成一维表示后,必须用reshape()函数将它重新转换回来。运行上面代码,可得原书P15 Figure1-8中的结果,即:
这里写图片描述

1.3.6 Pickle模块

如果想要保存一些结果或者数据以方便后续使用,Python 中的 pickle 模块非常有用。pickle模块可以接受几乎所有的 Python 对象,并且将其转换成字符串表示,该过程叫做封装(pickling)。从字符串表示中重构该对象,称为拆封(unpickling)。这些字符串表示可以方便地存储和传输。

我们来看一个例子。假设想要保存上一节字体图像的平均图像和主成分,可以这样来完成:

# 保存均值和主成分数据
f = open('font_pca_modes.pkl','wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()

在上述例子中,许多对象可以保存到同一个文件中。pickle 模块中有很多不同的协议可以生成 .pkl 文件;如果不确定的话,最好以二进制文件的形式读取和写入。在其他 Python 会话中载入数据,只需要如下使用 load() 方法:

# 载入均值和主成分数据
f = open('font_pca_modes.pkl','rb')
immean = pickle.load(f)
V = pickle.load(f)
f.close()

注意,载入对象的顺序必须和先前保存的一样。Python 中有个用 C 语言写的优化版本,叫做cpickle 模块,该模块和标准 pickle 模块完全兼容。关于 pickle 模块的更多内容,参见pickle 模块文档页 http://docs.python.org/library/pickle.html

在本书接下来的章节中,我们将使用 with 语句处理文件的读写操作。这是 Python 2.5 引入的思想,可以自动打开和关闭文件(即使在文件打开时发生错误)。下面的例子使用 with() 来实现保存和载入操作:

# 打开文件并保存
with open('font_pca_modes.pkl', 'wb') as f:
  pickle.dump(immean,f)
  pickle.dump(V,f)

# 打开文件并载入
with open('font_pca_modes.pkl', 'rb') as f:
  immean = pickle.load(f)
  V = pickle.load(f)

上面的例子乍看起来可能很奇怪,但 with() 确实是个很有用的思想。如果你不喜欢它,可以使用之前的 open 和 close 函数。

作为 pickle 的一种替代方式,NumPy 具有读写文本文件的简单函数。如果数据中不包含复杂的数据结构,比如在一幅图像上点击的点列表,NumPy 的读写函数会很有用。保存一个数组 x 到文件中,可以使用:

savetxt('test.txt',x,'%i')

最后一个参数表示应该使用整数格式。类似地,读取可以使用:

x = loadtxt('test.txt')

可以从在线文档了解更多

最后,NumPy 有专门用于保存和载入数组的函数,在线文档中可以查看关于 save()load() 的更多内容。

1.4 SciPy

SciPy(http://scipy.org/) 是建立在 NumPy 基础上,用于数值运算的开源工具包。SciPy 提供很多高效的操作,可以实现数值积分、优化、统计、信号处理,以及对我们来说最重要的图像处理功能。

1.4.1 图像模糊

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

Iδ=IGδ

其中,*表示卷积,Gδ表示标准差为δ的卷积核

  • 滤波操作模块——scipy.ndimage.filters

该模块可以使用快速一维分离的方式来计算卷积,使用方式如下:

from PIL import Image
from numpy import *
from pylab import *
from scipy.ndimage import filters

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font=FontProperties(fname=r"c:\windows\fonts\SimSun.ttc",size=14)

im=array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))

figure()
gray()
axis('off')
subplot(141)
axis('off')
title(u'原图',fontproperties=font)
imshow(im)

for bi,blur in enumerate([2,5,10]):
    im2=zeros(im.shape)
    im2=filters.gaussian_filter(im,blur)
    im2=np.uint8(im2)
    imNum=str(blur)
    subplot(1,4,2+bi)
    axis('off')
    title(u'标准差为'+imNum,fontproperties=font)
    imshow(im2)

#如果是彩色图像,则分别对三个通道进行模糊
#for bi, blur in enumerate([2, 5, 10]):
#  im2 = zeros(im.shape)
#  for i in range(3):
#    im2[:, :, i] = filters.gaussian_filter(im[:, :, i], blur)
#  im2 = np.uint8(im2)
#  subplot(1, 4,  2 + bi)
#  axis('off')
#  imshow(im2)

show()

这里写图片描述

上面第一幅图为待模糊图像,第二幅用高斯标准差为2进行模糊,第三幅用高斯标准差为5进行模糊,最后一幅用高斯标准差为10进行模糊。关于该模块的使用以及参数选择的更多细节,可以参阅SciPy scipy.ndimage文档

1.4.2 图像导数

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

  • 图像的梯度向量为I=[Ix,Iy]T,描述图像在每个像素点上强度变化最大的方向。
  • 梯度有两个重要的属性:
    1. 梯度的大小:
      |I|=Ix2+Iy2
    2. 梯度的方向:
      α=arctan2(Ix,Iy)

NumPy中的arctan2()函数返回弧度表示的有符号角度,角度的变化区间为[π,π]

我们可以用离散近似的方式来计算图像的导数。图像导数大多数可以通过卷积简单地实现:
Ix=IDxIy=IDy
对于,通常选择prewitt滤波器或sobel滤波器
这些导数滤波器可以使用scipy.ndimage.filters模块的标准卷积操作来简单实现

from PIL import Image
from pylab import *
from scipy.ndimage import  filters
import numpy

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font=FontProperties(fname=r"c:\windows\fonts\SimSun.ttc",size=14)

im=array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))
gray()

subplot(141)
axis('off')
title(u'(a)原图',fontproperties=font)
imshow(im)

# sobel derivative filters
imx=zeros(im.shape)
filters.sobel(im,1,imx)
subplot(142)
axis('off')
title(u'(b)x方向差分',fontproperties=font)
imshow(imx)

imy=zeros(im.shape)
filters.sobel(im,0,imy)
subplot(143)
axis('off')
title(u'(c)y方向差分',fontproperties=font)
imshow(imy)

mag=255-numpy.sqrt(imx**2+imy**2)
subplot(144)
title(u'(d)梯度幅值',fontproperties=font)
axis('off')
imshow(mag)

show()

这里写图片描述

高斯差分:

from PIL import Image
from pylab import *
from scipy.ndimage import filters
import numpy

# 添加中文字体支持
#from matplotlib.font_manager import FontProperties
#font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

def imx(im, sigma):
    imgx = zeros(im.shape)
    filters.gaussian_filter(im, sigma, (0, 1), imgx)
    return imgx


def imy(im, sigma):
    imgy = zeros(im.shape)
    filters.gaussian_filter(im, sigma, (1, 0), imgy)
    return imgy


def mag(im, sigma):
    # there's also gaussian_gradient_magnitude()
    #mag = numpy.sqrt(imgx**2 + imgy**2)
    imgmag = 255 - numpy.sqrt(imgx ** 2 + imgy ** 2)
    return imgmag


im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))
figure()
gray()

sigma = [2, 5, 10]

for i in  sigma:
    subplot(3, 4, 4*(sigma.index(i))+1)
    axis('off')
    imshow(im)
    imgx=imx(im, i)
    subplot(3, 4, 4*(sigma.index(i))+2)
    axis('off')
    imshow(imgx)
    imgy=imy(im, i)
    subplot(3, 4, 4*(sigma.index(i))+3)
    axis('off')
    imshow(imgy)
    imgmag=mag(im, i)
    subplot(3, 4, 4*(sigma.index(i))+4)
    axis('off')
    imshow(imgmag)

show()

这里写图片描述

1.4.3 形态学:对象计数

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

scipy.ndimage 中的 morphology 模块可以实现形态学操作
scipy.ndimage 中的measurements 模块来实现二值图像的计数和度量功能

下面通过一个简单的例子介绍如何使用它们:

from scipy.ndimage import measurements,morphology

# 载入图像,然后使用阈值化操作,以保证处理的图像为二值图像
im = array(Image.open('houses.png').convert('L'))
im = 1*(im<128)

labels, nbr_objects = measurements.label(im)
print "Number of objects:", nbr_objects
  1. 上面的脚本首先载入该图像,通过阈值化方式来确保该图像是二值图像。通过和 1 相乘,脚本将布尔数组转换成二进制表示。
  2. 然后,我们使用 label() 函数寻找单个的物体,并且按照它们属于哪个对象将整数标签给像素赋值。
  3. 图 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 模块文档 中了解关于这些函数的更多知识。
from PIL import Image
from numpy import *
from scipy.ndimage import measurements, morphology
from pylab import *

"""   This is the morphology counting objects example in Section 1.4.  """

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

# load image and threshold to make sure it is binary
figure()
gray()
im = array(Image.open('E:/python/Python Computer Vision/Image data/houses.png').convert('L'))
subplot(221)
imshow(im)
axis('off')
title(u'原图', fontproperties=font)
im = (im < 128)

labels, nbr_objects = measurements.label(im)
print ("Number of objects:", nbr_objects)
subplot(222)
imshow(labels)
axis('off')
title(u'标记后的图', fontproperties=font)

# morphology - opening to separate objects better
im_open = morphology.binary_opening(im, ones((9, 5)), iterations=2)
subplot(223)
imshow(im_open)
axis('off')
title(u'开运算后的图像', fontproperties=font)

labels_open, nbr_objects_open = measurements.label(im_open)
print ("Number of objects:", nbr_objects_open)
subplot(224)
imshow(labels_open)
axis('off')
title(u'开运算后进行标记后的图像', fontproperties=font)

show()

输出:

Number of objects: 45
Number of objects: 48

这里写图片描述

1.4.4 有用的SciPy模块

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

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 模块的更多内容,请参见在线文档

2.以图像形式保存数组

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

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

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

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

lena = scipy.misc.lena()

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

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

1.5 高级示例:图像去噪

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

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

降噪综合示例:

from pylab import *
from numpy import *
from numpy import random
from scipy.ndimage import filters
from scipy.misc import imsave
from PCV.tools import rof

""" This is the de-noising example using ROF in Section 1.5. """

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

# create synthetic image with noise
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)


# save the result
#imsave('synth_original.pdf',im)
#imsave('synth_rof.pdf',U)
#imsave('synth_gaussian.pdf',G)


# plot
figure()
gray()

subplot(1,3,1)
imshow(im)
#axis('equal')
axis('off')
title(u'原噪声图像', fontproperties=font)

subplot(1,3,2)
imshow(G)
#axis('equal')
axis('off')
title(u'高斯模糊后的图像', fontproperties=font)

subplot(1,3,3)
imshow(U)
#axis('equal')
axis('off')
title(u'ROF降噪后的图像', fontproperties=font)

show()

这里写图片描述

其中第一幅图示原噪声图像,中间一幅图示用标准差为10进行高斯模糊后的结果,最右边一幅图是用ROF降噪后的图像。上面原噪声图像是模拟出来的图像,现在我们在真实的图像上进行测试:

from PIL import Image
from pylab import *
from numpy import *
from numpy import random
from scipy.ndimage import filters
from scipy.misc import imsave
from PCV.tools import rof

""" This is the de-noising example using ROF in Section 1.5. """

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

im = array(Image.open('E:/python/Python Computer Vision/Image data/empire.jpg').convert('L'))

U,T = rof.denoise(im,im)
G = filters.gaussian_filter(im,10)


# save the result
#imsave('synth_original.pdf',im)
#imsave('synth_rof.pdf',U)
#imsave('synth_gaussian.pdf',G)


# plot
figure()
gray()

subplot(1,3,1)
imshow(im)
#axis('equal')
axis('off')
title(u'原噪声图像', fontproperties=font)

subplot(1,3,2)
imshow(G)
#axis('equal')
axis('off')
title(u'高斯模糊后的图像', fontproperties=font)

subplot(1,3,3)
imshow(U)
#axis('equal')
axis('off')
title(u'ROF降噪后的图像', fontproperties=font)

show()

这里写图片描述
ROF降噪能够保持边缘和图像结构

Number of objects: 45
Number of objects: 48

这里写图片描述

2018-11-18 15:02:40 qq_37482202 阅读数 20701

质量、速度、廉价,选择其中两个

提到图像处理第一个想到的库就是PIL,全称Python Imaging Library Python,图像处理类库,它提供了大量的图像操作,比如图像缩放,裁剪,贴图,模糊等等,很多时候它需要配合numpy库一起使用

1.open()

你可以使用Image.open打开一个图像文件,它会返回PIL图像对象

image = Image.open(image_address)

2.covert()

你可以 covert() 方法转换图像格式,covert() 有三种传参方式

im.convert(mode) ⇒ image

im.convert(“P”, **options) ⇒ image

im.convert(mode, matrix) ⇒ image

最常用的还是第一种,通过该方法你可以将PIL图像转换成九种不同的格式,分别1,L,P,RGB,RGBA,CMYK,YCbCr,I,F。

1.模式“1”

模式“1”为二值图像,非黑即白。但是它每个像素用8个bit表示,0表示黑,255表示白。

2.模式“L”

模式”L”为灰色图像,它的每个像素用8个bit表示,0表示黑,255表示白,其他数字表示不同的灰度。在PIL中,从模式“RGB”转换为“L”模式是按照下面的公式转换的:

L = R * 299/1000 + G * 587/1000+ B * 114/1000

3.模式“p”

模式“P”为8位彩色图像,它的每个像素用8个bit表示,其对应的彩色值是按照调色板查询出来的。

4.模式“RGBA”

模式“RGBA”为32位彩色图像,它的每个像素用32个bit表示,其中24bit表示红色、绿色和蓝色三个通道,另外8bit表示alpha通道,即透明通道。

5.模式“CMYK”

模式“CMYK”为32位彩色图像,它的每个像素用32个bit表示。模式“CMYK”就是印刷四分色模式,它是彩色印刷时采用的一种套色模式,利用色料的三原色混色原理,加上黑色油墨,共计四种颜色混合叠加,形成所谓“全彩印刷”。

四种标准颜色是:C:Cyan = 青色,又称为‘天蓝色’或是‘湛蓝’M:Magenta = 品红色,又称为‘洋红色’;Y:Yellow = 黄色;K:Key Plate(blacK) = 定位套版色(黑色)。

6.模式“YCbCr”

模式“YCbCr”为24位彩色图像,它的每个像素用24个bit表示。YCbCr其中Y是指亮度分量,Cb指蓝色色度分量,而Cr指红色色度分量。人的肉眼对视频的Y分量更敏感,因此在通过对色度分量进行子采样来减少色度分量后,肉眼将察觉不到的图像质量的变化。

模式“RGB”转换为“YCbCr”的公式如下:

Y= 0.257*R+0.504*G+0.098*B+16
Cb = -0.148*R-0.291*G+0.439*B+128
Cr = 0.439*R-0.368*G-0.071*B+128

7.模式“I”

模式“I”为32位整型灰色图像,它的每个像素用32个bit表示,0表示黑,255表示白,(0,255)之间的数字表示不同的灰度。在PIL中,从模式“RGB”转换为“I”模式是按照下面的公式转换的:

I = R * 299/1000 + G * 587/1000 + B * 114/1000

8.模式“F”

模式“F”为32位浮点灰色图像,它的每个像素用32个bit表示,0表示黑,255表示白,(0,255)之间的数字表示不同的灰度。在PIL中,从模式“RGB”转换为“F”模式是按照下面的公式转换的:

F = R * 299/1000+ G * 587/1000 + B * 114/1000

3.调整尺寸、创建缩略图、裁剪、贴图、旋转

PIL库给我们提供了丰富基本图像操作,如果你想调整一张图片的尺寸,你可以使用resize()方法,该方法需要传入你指定新图像宽高的元组

img = img.resize((128,128))

如果你想创建一张图片的缩略图,你可以使用thumbnail()方法,该方法需要传入缩略图的宽高元组

img=img.thumbnail((128,128))

如果你想对一张图片的一部分进行裁剪,你可以使用crop()方法,该方法需要你传入一个元组,该元组指定裁剪区域的左上角坐标和右下角坐标

box = (100,100,400,400)
img = img.crop(box)

如果你想把一张图片覆盖在另一个图片的上面,你可以使用paste()方法,该方法需要传入要贴的图片和位置(左上角坐标和右下角坐标)

img2=img2.paste(img1,(100,100,200,200))

如果你想要旋转一张图片,你可以使用transpose()方法,该方法传入旋转角度

img = img.transpose(Image.ROTATE_180)

不过这些角度很受限制,只可以传下面之中的一个

  • PIL.Image.FLIP_LEFT_RIGHT 
  • PIL.Image.FLIP_TOP_BOTTOM
  • PIL.Image.ROTATE_90
  • PIL.Image.ROTATE_180
  • PIL.Image.ROTATE_270
  • PIL.Image.TRANSPOSE
  • PIL.Image.TRANSVERSE

你也可以使用rotate()方法,该方法更为简单方便,只需要传入一个旋转角度即可

image = image.rotate(45)

4.Numpy

对图像进行变换其实就是对矩阵进行变换,我们需要把一张图片转换成矩阵再进行操作,使用array()方法

image = Image.open(image_address)
imageArray = array(image)

1.反向处理与二值化

图像一般都是三通道的,也就是红绿蓝,他们的值从0-255,所谓反相处理呢,就是把颜色反过来

imageArray = 255 - imageArray

 图像的二值化也很简单,0-255以128为分界,小于128置为0否则置为1

imageArray = 1 * (imageArray < 128)

2.像素值限制范围

如果你想把一个图像的像素值都限制到一个范围内,比如说你想把像素值限制到100-200这个区间上,你可以这么干

imageArray = (100.0 / 255) * imageArray + 100

3.像素值求平方

imageArray = 255.0 * (imageArray / 255.0) ** 2

4.直方图均衡化

图像灰度变换中一个非常有用的例子就是直方图均衡化。直方图均衡化是指将一幅图像的灰度直方图变平,使变换后的图像中每个灰度值的分布概率都相同。在对图像做进一步处理之前,直方图均衡化通常是对图像灰度值进行归一化的一个非常好的方法,并且可以增强图像的对比度。

在这种情况下,直方图均衡化的变换函数是图像中像素值的累积分布函数(cumulative distribution function,简写为 cdf,将像素值的范围映射到目标范围的归一化操作)

def histeq(im,nbr_bins=256):
    """ 对一幅灰度图像进行直方图均衡化"""
    # 计算图像的直方图
    imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
    cdf = imhist.cumsum()
    # cumulative distribution function
    cdf = 255 * cdf / cdf[-1]
    #  归一化
    #  使用累积分布函数的线性插值,计算新的像素值
    im2 = interp(im.flatten(),bins[:-1],cdf)
    return im2.reshape(im.shape), cdf

该函数有两个输入参数,一个是灰度图像,一个是直方图中使用小区间的数目。函数返回直方图均衡化后的图像,以及用来做像素值映射的累积分布函数。注意,函数中使用到累积分布函数的最后一个元素(下标为 -1),目的是将其归一化到 0...1 范围。

直方图均衡化后图像可以使对比度增强,使原先图像灰色区域的细节变得更清晰

5.多种滤波

gaussian滤波是多维的滤波器,是一种平滑滤波,可以消除高斯噪声

通过调节sigma的值来调整滤波效果

imageArray = filters.gaussian_filter(imageArray, 5)

sobel算子可用来检测边缘

edges = filters.sobel(img)

roberts算子、scharr算子、prewitt算子和sobel算子一样,用于检测边缘

edges = filters.roberts(img)
edges = filters.scharr(img)
edges = filters.prewitt(img)

canny算子也是用于提取边缘特征,但它不是放在filters模块,而是放在feature模块

edges1 = feature.canny(img)   #sigma=1
edges2 = feature.canny(img,sigma=3)   #sigma=3

gabor滤波可用来进行边缘检测和纹理特征提取。

通过修改frequency值来调整滤波效果,返回一对边缘结果,一个是用真实滤波核的滤波结果,一个是想象的滤波核的滤波结果。

filt_real, filt_imag = filters.gabor_filter(img,frequency=0.6)   

6.PCA

PCA(Principal Component Analysis,主成分分析)是一个非常有用的降维技巧。它可以在使用尽可能少维数的前提下,尽量多地保持训练数据的信息,在此意义上是一个最佳技巧。即使是一幅 100×100 像素的小灰度图像,也有 10 000 维,可以看成 10 000 维空间中的一个点。一兆像素的图像具有百万维。由于图像具有很高的维数,在许多计算机视觉应用中,我们经常使用降维操作。PCA 产生的投影矩阵可以被视为将原始坐标变换到现有的坐标系,坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行 PCA 变换,图像需要转换成一维向量表示。我们可以使用 NumPy 类库中的 flatten() 方法进行变换。

将变平的图像堆积起来,我们可以得到一个矩阵,矩阵的一行表示一幅图像。在计算主方向之前,所有的行图像按照平均图像进行了中心化。我们通常使用 SVD(Singular Value Decomposition,奇异值分解)方法来计算主成分;但当矩阵的维数很大时,SVD 的计算非常慢,所以此时通常不使用 SVD 分解。下面就是 PCA 操作的代码:

def pca(X):
    """ 主成分分析:    输入:矩阵X ,其中该矩阵中存储训练数据,每一行为一条训练数据
       返回:投影矩阵(按照维度的重要性排序)、方差和均值"""
    # 获取维数
    num_data,dim = X.shape
    # 数据中心化
    mean_X = X.mean(axis=0)
    X = X - mean_X
    if dim<num_data:
        # PCA- 使用紧致技巧
        M = dot(X,X.T)
        # 协方差矩阵
        e,EV = linalg.eigh(M)
        # 特征值和特征向量
        tmp = dot(X.T,EV).T
        # 这就是紧致技巧
        V = tmp[::-1]
        # 由于最后的特征向量是我们所需要的,所以需要将其逆转
        S = sqrt(e)[::-1]
        # 由于特征值是按照递增顺序排列的,所以需要将其逆转
        for i in range(V.shape[1]):
            V[:,i] /= S
    else:
        # PCA- 使用SVD 方法
        U,S,V = linalg.svd(X)
        V = V[:num_data]
        # 仅仅返回前nun_data 维的数据才合理
        #  返回投影矩阵、方差和均值
    return V,S,mean_X

7.图像添加噪声和降噪

添加噪声比降噪简单得多,只需要把图像矩阵上面随机加一些值就好了

imageArray = imageArray + 30 * random.standard_normal(imageArray.shape)

图像降噪是在去除图像噪声的同时,尽可能地保留图像细节和结构的处理技术,我们这里使用 ROF去燥模型

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

J(\boldsymbol{I})=\int\left|\nabla\boldsymbol{I}\right|\text{dx} 

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

J(\boldsymbol{I})=\sum_{\text{x}}\left|\nabla\boldsymbol{I}\right|

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

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

\min_U\left|\left|\boldsymbol{I}-\boldsymbol{U}\right|\right|^2+2\lambda J(\boldsymbol{U}),

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

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  # 去噪后的图像和纹理残余

5.Matplotlib

我们队图像进行处理之后往往需要知道处理后变化如何,该库便可以方便地绘制出条形图,饼状图等等呢个图像,还可在上面添加标记等等

尽管 Matplotlib 可以绘制出较好的条形图、饼状图、散点图等,但是对于大多数计算机视觉应用来说,仅仅需要用到几个绘图命令。最重要的是,我们想用点和线来表示一些事物,比如兴趣点、对应点以及检测出的物体。下面是用几个点和一条线绘制图像的例子:

from PIL import Image
from pylab import *

# 读取图像到数组中
im = array(Image.open('empire.jpg'))

# 绘制图像
imshow(im)

# 一些点
x = [100,100,400,400]
y = [200,500,200,500]

# 使用红色星状标记绘制点
plot(x,y,'r*')

# 绘制连接前两个点的线
plot(x[:2],y[:2])

# 添加标题,显示绘制的图像
title('Plotting: "empire.jpg"')
show()

上面的代码首先绘制出原始图像,然后在 x 和 y 列表中给定点的 x 坐标和 y 坐标上绘制出红色星状标记点,最后在两个列表表示的前两个点之间绘制一条线段(默认为蓝色)。该例子的绘制结果如图 1-2 所示。show() 命令首先打开图形用户界面(GUI),然后新建一个图像窗口。该图形用户界面会循环阻断脚本,然后暂停,直到最后一个图像窗口关闭。在每个脚本里,你只能调用一次 show() 命令,而且通常是在脚本的结尾调用。注意,在 PyLab 库中,我们约定图像的左上角为坐标原点。

图像的坐标轴是一个很有用的调试工具;但是,如果你想绘制出较美观的图像,加上下列命令可以使坐标轴不显示:

axis('off')

 下面是我写的一个图像处理的脚本

import PIL.Image as Image
import os
from pylab import *
from numpy import *
from scipy.ndimage import filters
from scipy.ndimage import measurements,morphology


def get_imlist(path):
    # 一级文件夹下有用
    # return [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]
    g = os.walk(path)
    image_list=[]
    for path, d, filelist in g:
        for filename in filelist:
            if filename.endswith('jpg'):
                image_list.append(os.path.join(path, filename))
    return image_list

def histeq(im,nbr_bins=256):
    """ 对一幅灰度图像进行直方图均衡化"""
    # 计算图像的直方图
    imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
    cdf = imhist.cumsum()
    # cumulative distribution function
    cdf = 255 * cdf / cdf[-1]
    #  归一化
    #  使用累积分布函数的线性插值,计算新的像素值
    im2 = interp(im.flatten(),bins[:-1],cdf)
    return im2.reshape(im.shape), cdf

def pca(X):
    """ 主成分分析:    输入:矩阵X ,其中该矩阵中存储训练数据,每一行为一条训练数据
       返回:投影矩阵(按照维度的重要性排序)、方差和均值"""
    # 获取维数
    num_data,dim = X.shape
    # 数据中心化
    mean_X = X.mean(axis=0)
    X = X - mean_X
    if dim<num_data:
        # PCA- 使用紧致技巧
        M = dot(X,X.T)
        # 协方差矩阵
        e,EV = linalg.eigh(M)
        # 特征值和特征向量
        tmp = dot(X.T,EV).T
        # 这就是紧致技巧
        V = tmp[::-1]
        # 由于最后的特征向量是我们所需要的,所以需要将其逆转
        S = sqrt(e)[::-1]
        # 由于特征值是按照递增顺序排列的,所以需要将其逆转
        for i in range(V.shape[1]):
            V[:,i] /= S
    else:
        # PCA- 使用SVD 方法
        U,S,V = linalg.svd(X)
        V = V[:num_data]
        # 仅仅返回前nun_data 维的数据才合理
        #  返回投影矩阵、方差和均值
    return V,S,mean_X

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  # 去噪后的图像和纹理残余




image_list = get_imlist("G:\\最后两种\\")

index=6858
for image_address in image_list:
    index = index + 1
    dealIndex=0
    for x in range(1,17):
        image = Image.open(image_address)
        imageArray = array(image)
        dealIndex+=1
        if x==1:
            # 反相处理
            imageArray = 255 - imageArray
            print("第"+str(index)+"张 反向处理")
        elif x==2:
            # 将图像像素值变换到100...200 区间
            imageArray = (100.0 / 255) * imageArray + 100
            print("第" + str(index) + "张 像素值变换")
        elif x==3:
            # 对图像像素值求平方后得到的图像
            imageArray = 255.0 * (imageArray / 255.0) ** 2
            print("第" + str(index) + "张 像素值求平方")
        elif x==4:
            # 图像旋转
            image = image.rotate(random.randint(0,360))
            imageArray=array(image)
            print("第" + str(index) + "张 图像旋转")
        elif x==5:
            # 直方图均衡化
            imageArray,cdf=histeq(imageArray)
            print("第" + str(index) + "张 直方图均衡化")
        elif x==6:
            # gaussian滤波
            imageArray = filters.gaussian_filter(imageArray, 5)
            print("第" + str(index) + "张 gaussian滤波")
        elif x==7:
            # Sobel 导数滤波器
            imx = zeros(imageArray.shape)
            filters.sobel(imageArray, 1, imx)
            imy = zeros(imageArray.shape)
            filters.sobel(imageArray, 0, imy)
            magnitude = sqrt(imx ** 2 + imy ** 2)
            imageArray=magnitude
            print("第" + str(index) + "张  Sobel导数滤波器")
        elif x==8:
            # 噪声
            imageArray = imageArray + 30 * random.standard_normal(imageArray.shape)
            print("第" + str(index) + "张  噪声")
        elif x==9:
            # 反相处理+像素值变换
            imageArray = 255 - imageArray
            imageArray = (100.0 / 255) * imageArray + 100
            print("第" + str(index) + "张  反相处理+像素值变换")
        elif x==10:
            # 反相处理+像素值求平方
            imageArray = 255 - imageArray
            imageArray = 255.0 * (imageArray / 255.0) ** 2
            print("第" + str(index) + "张  反相处理+像素值求平方")
        elif x==11:
            # 像素值求平方+反相处理
            imageArray = 255.0 * (imageArray / 255.0) ** 2
            imageArray = 255 - imageArray
            print("第" + str(index) + "张  像素值求平方+反相处理")
        elif x==12:
            # 像素值变换+像素值求平方
            imageArray = (100.0 / 255) * imageArray + 100
            imageArray = 255.0 * (imageArray / 255.0) ** 2
            print("第" + str(index) + "张  像素值变换+像素值求平方")
        elif x==13:
            # 图像旋转+反相
            image = image.rotate(random.randint(0, 360))
            imageArray = array(image)
            imageArray = 255 - imageArray
            print("第" + str(index) + "张  图像旋转+反相")
        elif x==14:
            # 图像旋转+噪声
            image = image.rotate(random.randint(0, 360))
            imageArray = array(image)
            imageArray = imageArray + 30 * random.standard_normal(imageArray.shape)
            print("第" + str(index) + "张  图像旋转+噪声")
        elif x==15:
            # 噪声+直方图均衡化
            imageArray = imageArray + 30 * random.standard_normal(imageArray.shape)
            imageArray, cdf = histeq(imageArray)
            print("第" + str(index) + "张  噪声+直方图均衡化")


        imageArray = uint8(imageArray)
        image=Image.fromarray(imageArray)
        image = image.convert('RGB')
        if image_address.rfind("不规则")!= -1:
            image.save("G:\\兔屎图片_二次处理\\不规则\\" + str(index)+"_"+str(dealIndex) + ".jpg")
        elif image_address.rfind("大小不一") != -1:
            image.save("G:\\兔屎图片_二次处理\\大小不一\\" + str(index) +"_"+str(dealIndex)+ ".jpg")
        elif image_address.rfind("拉稀") != -1:
            image.save("G:\\兔屎图片_二次处理\\拉稀\\" + str(index) +"_"+str(dealIndex)+ ".jpg")
        elif image_address.rfind("正常") != -1:
            image.save("G:\\兔屎图片_二次处理\\正常\\" + str(index) +"_"+str(dealIndex)+ ".jpg")

        print("完事一个")

参考文章:https://www.cnblogs.com/xk-bench/p/7825290.html

2007-10-28 23:45:00 lanphaday 阅读数 102269
用Python做图像处理
       最近在做一件比较 evil 的事情——验证码识别,以此来学习一些新的技能。因为我是初学,对图像处理方面就不太了解了,欲要利吾事,必先利吾器,既然只是做一下实验,那用 Python 来作原型开发再好不过了。在 Python 中,比较常用的图像处理库是 PIL(Python Image Library),当前版本是 1.1.6 ,用起来非常方便。大家可以在 http://www.pythonware.com/products/pil/index.htm 下载和学习。
       在这里,我主要是介绍一下做图像识别时可能会用到的一些 PIL 提供的功能,比如图像增强、还有滤波之类的。最后给出使用 Python 做图像处理与识别的优势与劣势。
基本图像处理
       使用 PIL 之前需要 import Image 模块:
import Image
       然后你就可以使用Image.open(‘xx.bmp’) 来打开一个位图文件进行处理了。打开文件你不用担心格式,也不用了解格式,无论什么格式,都只要把文件名丢给 Image.open 就可以了。真所谓 bmp、jpg、png、gif……,一个都不能少。
img = Image.open(‘origin.png’)    # 得到一个图像的实例对象 img
1原图
       图像处理中,最基本的就是色彩空间的转换。一般而言,我们的图像都是 RGB 色彩空间的,但在图像识别当中,我们可能需要转换图像到灰度图、二值图等不同的色彩空间。 PIL 在这方面也提供了极完备的支持,我们可以:
new_img = img.convert(‘L’)
把 img 转换为 256 级灰度图像, convert() 是图像实例对象的一个方法,接受一个 mode 参数,用以指定一种色彩模式,mode 的取值可以是如下几种:
· 1 (1-bit pixels, black and white, stored with one pixel per byte)
· L (8-bit pixels, black and white)
· P (8-bit pixels, mapped to any other mode using a colour palette)
· RGB (3x8-bit pixels, true colour)
· RGBA (4x8-bit pixels, true colour with transparency mask)
· CMYK (4x8-bit pixels, colour separation)
· YCbCr (3x8-bit pixels, colour video format)
· I (32-bit signed integer pixels)
· F (32-bit floating point pixels)
怎么样,够丰富吧?其实如此之处,PIL 还有限制地支持以下几种比较少见的色彩模式:LA (L with alpha), RGBX (true colour with padding) and RGBa (true colour with premultiplied alpha)。
下面看一下 mode 为 ‘1’、’L’、’P’时转换出来的图像:
2 mode = '1'
3 mode = 'L'
4 mode = 'P'
convert() 函数也接受另一个隐含参数 matrix,转换矩阵 matrix 是一个长度为4 或者16 tuple。下例是一个转换 RGB 空间到 CIE XYZ 空间的例子:
    rgb2xyz = (
        0.412453, 0.357580, 0.180423, 0,
        0.212671, 0.715160, 0.072169, 0,
        0.019334, 0.119193, 0.950227, 0 )
    out = im.convert("RGB", rgb2xyz)
       除了完备的色彩空间转换能力外, PIL 还提供了resize()、rotate()等函数以获得改变大小,旋转图片等几何变换能力,在图像识别方面,图像实例提供了一个 histogram() 方法来计算直方图,非常方便实用。
图像增强
       图像增强通常用以图像识别之前的预处理,适当的图像增强能够使得识别过程达到事半功倍的效果。 PIL 在这方面提供了一个名为 ImageEnhance 的模块,提供了几种常见的图像增强方案:
import ImageEnhance
enhancer = ImageEnhance.Sharpness(image)
for i in range(8):
    factor = i / 4.0
    enhancer.enhance(factor).show("Sharpness %f" % factor)
上面的代码即是一个典型的使用 ImageEnhance 模块的例子。 Sharpness 是 ImageEnhance 模块的一个类,用以锐化图片。这一模块主要包含如下几个类:Color、Brightness、Contrast和Sharpness。它们都有一个共同的接口 .enhance(factor) ,接受一个浮点参数 factor,标示增强的比例。下面看看这四个类在不同的 factor 下的效果
5 使用Color 进行色彩增强,factor 取值 [0, 4],步进 0.5
6 用 Birghtness 增强亮度,factor取值[0,4],步进0.5
7用 Contrast 增强对比度, factor 取值 [0,4],步进0.5
8用 Sharpness 锐化图像,factor取值 [0,4],步进0.5
图像 Filter
       PIL 在 Filter 方面的支持是非常完备的,除常见的模糊、浮雕、轮廓、边缘增强和平滑,还有中值滤波、ModeFilter等,简直方便到可以做自己做一个Photoshop。这些 Filter 都放置在 ImageFilter 模块中,ImageFilter主要包括两部分内容,一是内置的 Filter,如 BLUR、DETAIL等,另一部分是 Filter 函数,可以指定不同的参数获得不同的效果。示例如下:
import ImageFilter
im1 = im.filter(ImageFilter.BLUR)
im2 = im.filter(ImageFilter.MinFilter(3))
im3 = im.filter(ImageFilter.MinFilter()) # same as MinFilter(3)
可以看到 ImageFilter 模块的使用非常简单,每一个 Filter 都只需要一行代码就可调用,开发效率非常高。
 
9使用 BLUR
10使用 CONTOUR
11使用 DETAIL
12使用 EMBOSS
13使用 EDGE_ENHANCE
14使用 EDGE_ENHANCE_MORE
15使用 FIND_EDGES
16使用 SHARPEN
17使用 SMOOTH
18使用 SMOOTH_MORE
       以上是几种内置的 Filter 的效果图,除此之外, ImageFilter 还提供了一些 Filter 函数,下面我们来看看这些可以通过参数改变行为的 Filter 的效果:
19使用 Kernel(),参数:size = (3, 3), kernel = (0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)
20使用 MaxFilter,默认参数
21使用 MinFilter,默认参数
22使用 MedianFilter,默认参数
23使用 ModeFilter,参数 size = 3
24使用 RankFilter,参数 size = 3, rank = 3
小结
       到此,对 PIL 的介绍就告一段落了。总的来说,对于图像处理和识别,PIL 内建了强大的支持,从各种增强算法到 Filter ,都让人无法怀疑使用 Python 的可行性。 Python唯一的劣势在于执行时间过慢,特别是当实现一些计算量大的算法时候,需要极强的耐心。我曾用 Hough Transform(霍夫变换)来查找图像中的直线,纯 Python 的实现处理一个 340 * 100 的图片也要花去数秒时间(P4 3.0G + 1G memory)。但使用 PIL 无需关注图像格式、内建的图像增强算法和 Filter 算法,这些优点使 Python 适合用于构造原型和进行实验,在这两方面Python 比 matlab 更加方便。商业的图像识别产品开发,可以考虑已经被 boost accepted的来自 adobe 的开源 C++ 库 gil,可以兼顾执行性能和开发效率。
2018-06-28 19:08:01 weixin_42052460 阅读数 22395

第 1 章 基本的图像操作和处理

本章讲解操作和处理图像的基础知识,将通过大量示例介绍处理图像所需的 Python 工具包,并介绍用于读取图像、图像转换和缩放、计算导数、画图和保存结果等的基本工具。这些工具的使用将贯穿本书的剩余章节。

1.1 PIL:Python图像处理类库

PIL(Python Imaging Library Python,图像处理类库)提供了通用的图像处理功能,以及大量有用的基本图像操作,比如图像缩放、裁剪、旋转、颜色转换等。PIL 是免费的,可以从 http://www.pythonware.com/products/pil/ 下载。

利用 PIL 中的函数,我们可以从大多数图像格式的文件中读取数据,然后写入最常见的图像格式文件中。PIL 中最重要的模块为 Image。要读取一幅图像,可以使用:

from PIL import Image

pil_im = Image.open('empire.jpg')

上述代码的返回值 pil_im 是一个 PIL 图像对象。

图像的颜色转换可以使用 convert() 方法来实现。要读取一幅图像,并将其转换成灰度图像,只需要加上 convert('L'),如下所示:

pil_im = Image.open('empire.jpg').convert('L')

在 PIL 文档中有一些例子,参见 http://www.pythonware.com/library/pil/handbook/index.htm。这些例子的输出结果如图 1-1 所示。

图 1-1:用 PIL 处理图像的例子

1.1.1 转换图像格式

通过 save() 方法,PIL 可以将图像保存成多种格式的文件。下面的例子从文件名列表(filelist)中读取所有的图像文件,并转换成 JPEG 格式:

from PIL import Image
import os

for infile in filelist:
  outfile = os.path.splitext(infile)[0] + ".jpg"
  if infile != outfile:
    try:
      Image.open(infile).save(outfile)
    except IOError:
      print "cannot convert", infile

PIL 的 open() 函数用于创建 PIL 图像对象,save() 方法用于保存图像到具有指定文件名的文件。除了后缀变为“.jpg”,上述代码的新文件名和原文件名相同。PIL 是个足够智能的类库,可以根据文件扩展名来判定图像的格式。PIL 函数会进行简单的检查,如果文件不是 JPEG 格式,会自动将其转换成 JPEG 格式;如果转换失败,它会在控制台输出一条报告失败的消息。

本书会处理大量图像列表。下面将创建一个包含文件夹中所有图像文件的文件名列表。首先新建一个文件,命名为 imtools.py,来存储一些经常使用的图像操作,然后将下面的函数添加进去:

import os
def get_imlist(path):

""" 返回目录中所有JPG 图像的文件名列表"""

return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]

现在,回到 PIL。

1.1.2 创建缩略图

使用 PIL 可以很方便地创建图像的缩略图。thumbnail() 方法接受一个元组参数(该参数指定生成缩略图的大小),然后将图像转换成符合元组参数指定大小的缩略图。例如,创建最长边为 128 像素的缩略图,可以使用下列命令:

pil_im.thumbnail((128,128))

1.1.3 复制和粘贴图像区域

使用 crop() 方法可以从一幅图像中裁剪指定区域:

box = (100,100,400,400)
region = pil_im.crop(box)

该区域使用四元组来指定。四元组的坐标依次是(左,上,右,下)。PIL 中指定坐标系的左上角坐标为(0,0)。我们可以旋转上面代码中获取的区域,然后使用 paste() 方法将该区域放回去,具体实现如下:

region = region.transpose(Image.ROTATE_180)
pil_im.paste(region,box)

1.1.4 调整尺寸和旋转

要调整一幅图像的尺寸,我们可以调用 resize() 方法。该方法的参数是一个元组,用来指定新图像的大小:

out = pil_im.resize((128,128))

要旋转一幅图像,可以使用逆时针方式表示旋转角度,然后调用 rotate() 方法:

out = pil_im.rotate(45)

上述例子的输出结果如图 1-1 所示。最左端是原始图像,然后是灰度图像、粘贴有旋转后裁剪图像的原始图像,最后是缩略图。

1.2 Matplotlib

我们处理数学运算、绘制图表,或者在图像上绘制点、直线和曲线时,Matplotlib 是个很好的类库,具有比 PIL 更强大的绘图功能。Matplotlib 可以绘制出高质量的图表,就像本书中的许多插图一样。Matplotlib 中的 PyLab 接口包含很多方便用户创建图像的函数。Matplotlib 是开源工具,可以从 http://matplotlib.sourceforge.net/ 免费下载。该链接中包含非常详尽的使用说明和教程。下面的例子展示了本书中需要使用的大部分函数。

1.2.1 绘制图像、点和线

尽管 Matplotlib 可以绘制出较好的条形图、饼状图、散点图等,但是对于大多数计算机视觉应用来说,仅仅需要用到几个绘图命令。最重要的是,我们想用点和线来表示一些事物,比如兴趣点、对应点以及检测出的物体。下面是用几个点和一条线绘制图像的例子:

from PIL import Image
from pylab import *

# 读取图像到数组中
im = array(Image.open('empire.jpg'))

# 绘制图像
imshow(im)

# 一些点
x = [100,100,400,400]
y = [200,500,200,500]

# 使用红色星状标记绘制点
plot(x,y,'r*')

# 绘制连接前两个点的线
plot(x[:2],y[:2])

# 添加标题,显示绘制的图像
title('Plotting: "empire.jpg"')
show()

上面的代码首先绘制出原始图像,然后在 x 和 y 列表中给定点的 x 坐标和 y 坐标上绘制出红色星状标记点,最后在两个列表表示的前两个点之间绘制一条线段(默认为蓝色)。该例子的绘制结果如图 1-2 所示。show() 命令首先打开图形用户界面(GUI),然后新建一个图像窗口。该图形用户界面会循环阻断脚本,然后暂停,直到最后一个图像窗口关闭。在每个脚本里,你只能调用一次 show() 命令,而且通常是在脚本的结尾调用。注意,在 PyLab 库中,我们约定图像的左上角为坐标原点。

图像的坐标轴是一个很有用的调试工具;但是,如果你想绘制出较美观的图像,加上下列命令可以使坐标轴不显示:

axis('off')

上面的命令将绘制出如图 1-2 右边所示的图像。

图 1-2:Matplotlib 绘图示例。带有坐标轴和不带坐标轴的包含点和一条线段的图像

在绘图时,有很多选项可以控制图像的颜色和样式。最有用的一些短命令如表 1-1、表 1-2 和表 1-3 所示。使用方法见下面的例子:

plot(x,y)         # 默认为蓝色实线
plot(x,y,'r*')    # 红色星状标记
plot(x,y,'go-')   # 带有圆圈标记的绿线
plot(x,y,'ks:')   # 带有正方形标记的黑色虚线

表1-1:用PyLab库绘图的基本颜色格式命令

颜色

 

'b'

蓝色

'g'

绿色

'r'

红色

'c'

青色

'm'

品红

'y'

黄色

'k'

黑色

'w'

白色

表1-2:用PyLab库绘图的基本线型格式命令

线型

 

'-'

实线

'--'

虚线

':'

点线

表1-3:用PyLab库绘图的基本绘制标记格式命令

标记

 

'.'

'o'

圆圈

's'

正方形

'*'

星形

'+'

加号

'x'

叉号

1.2.2 图像轮廓和直方图

下面来看两个特别的绘图示例:图像的轮廓和直方图。绘制图像的轮廓(或者其他二维函数的等轮廓线)在工作中非常有用。因为绘制轮廓需要对每个坐标 [x, y] 的像素值施加同一个阈值,所以首先需要将图像灰度化:

from PIL import Image
from pylab import *

# 读取图像到数组中
im = array(Image.open('empire.jpg').convert('L'))

# 新建一个图像
figure()
# 不使用颜色信息
gray()
# 在原点的左上角显示轮廓图像
contour(im, origin='image')
axis('equal')
axis('off')

像之前的例子一样,这里用 PIL 的 convert() 方法将图像转换成灰度图像。

图像的直方图用来表征该图像像素值的分布情况。用一定数目的小区间(bin)来指定表征像素值的范围,每个小区间会得到落入该小区间表示范围的像素数目。该(灰度)图像的直方图可以使用 hist() 函数绘制:

figure()
hist(im.flatten(),128)
show()

hist() 函数的第二个参数指定小区间的数目。需要注意的是,因为 hist() 只接受一维数组作为输入,所以我们在绘制图像直方图之前,必须先对图像进行压平处理。flatten() 方法将任意数组按照行优先准则转换成一维数组。图 1-3 为等轮廓线和直方图图像。

图 1-3:用 Matplotlib 绘制图像等轮廓线和直方图

1.2.3 交互式标注

有时用户需要和某些应用交互,例如在一幅图像中标记一些点,或者标注一些训练数据。PyLab 库中的 ginput() 函数就可以实现交互式标注。下面是一个简短的例子:

from PIL import Image
from pylab import *

im = array(Image.open('empire.jpg'))
imshow(im)
print 'Please click 3 points'
x = ginput(3)
print 'you clicked:',x
show()

上面的脚本首先绘制一幅图像,然后等待用户在绘图窗口的图像区域点击三次。程序将这些点击的坐标 [x, y] 自动保存在 x 列表里。

1.3 NumPy

NumPyhttp://www.scipy.org/NumPy/)是非常有名的 Python 科学计算工具包,其中包含了大量有用的思想,比如数组对象(用来表示向量、矩阵、图像等)以及线性代数函数。NumPy 中的数组对象几乎贯穿用于本书的所有例子中 1 数组对象可以帮助你实现数组中重要的操作,比如矩阵乘积、转置、解方程系统、向量乘积和归一化,这为图像变形、对变化进行建模、图像分类、图像聚类等提供了基础。

1PyLab 实际上包含 NumPy 的一些内容,如数组类型。这也是我们能够在 1.2 节使用数组类型的原因。

NumPy 可以从 http://www.scipy.org/Download 免费下载,在线说明文档(http://docs.scipy.org/doc/numpy/)包含了你可能遇到的大多数问题的答案。关于 NumPy 的更多内容,请参考开源书籍 [24]。

1.3.1 图像数组表示

在先前的例子中,当载入图像时,我们通过调用 array() 方法将图像转换成 NumPy 的数组对象,但当时并没有进行详细介绍。NumPy 中的数组对象是多维的,可以用来表示向量、矩阵和图像。一个数组对象很像一个列表(或者是列表的列表),但是数组中所有的元素必须具有相同的数据类型。除非创建数组对象时指定数据类型,否则数据类型会按照数据的类型自动确定。

对于图像数据,下面的例子阐述了这一点:

im = array(Image.open('empire.jpg'))
print im.shape, im.dtype

im = array(Image.open('empire.jpg').convert('L'),'f')
print im.shape, im.dtype

控制台输出结果如下所示:

(800, 569, 3) uint8
(800, 569) float32

每行的第一个元组表示图像数组的大小(行、列、颜色通道),紧接着的字符串表示数组元素的数据类型。因为图像通常被编码成无符号八位整数(uint8),所以在第一种情况下,载入图像并将其转换到数组中,数组的数据类型为“uint8”。在第二种情况下,对图像进行灰度化处理,并且在创建数组时使用额外的参数“f”;该参数将数据类型转换为浮点型。关于更多数据类型选项,可以参考图书 [24]。注意,由于灰度图像没有颜色信息,所以在形状元组中,它只有两个数值。

数组中的元素可以使用下标访问。位于坐标 ij,以及颜色通道 k 的像素值可以像下面这样访问:

value = im[i,j,k]

多个数组元素可以使用数组切片方式访问。切片方式返回的是以指定间隔下标访问该数组的元素值。下面是有关灰度图像的一些例子:

im[i,:] = im[j,:]      # 将第 j 行的数值赋值给第 i 行
im[:,i] = 100          # 将第 i 列的所有数值设为100
im[:100,:50].sum()     # 计算前100 行、前 50 列所有数值的和
im[50:100,50:100]      # 50~100 行,50~100 列(不包括第 100 行和第 100 列)
im[i].mean()           # 第 i 行所有数值的平均值
im[:,-1]               # 最后一列
im[-2,:] (or im[-2])   # 倒数第二行

注意,示例仅仅使用一个下标访问数组。如果仅使用一个下标,则该下标为行下标。注意,在最后几个例子中,负数切片表示从最后一个元素逆向计数。我们将会频繁地使用切片技术访问像素值,这也是一个很重要的思想。

我们有很多操作和方法来处理数组对象。本书将在使用到的地方逐一介绍。你可以查阅在线文档或者开源图书 [24] 获取更多信息。

1.3.2 灰度变换

将图像读入 NumPy 数组对象后,我们可以对它们执行任意数学操作。一个简单的例子就是图像的灰度变换。考虑任意函数 f,它将 0...255 区间(或者 0...1 区间)映射到自身(意思是说,输出区间的范围和输入区间的范围相同)。下面是关于灰度变换的一些例子:

from PIL import Image
from numpy import *

im = array(Image.open('empire.jpg').convert('L'))

im2 = 255 - im # 对图像进行反相处理

im3 = (100.0/255) * im + 100 # 将图像像素值变换到100...200 区间

im4 = 255.0 * (im/255.0)**2 # 对图像像素值求平方后得到的图像

第一个例子将灰度图像进行反相处理;第二个例子将图像的像素值变换到 100...200 区间;第三个例子对图像使用二次函数变换,使较暗的像素值变得更小。图 1-4 为所使用的变换函数图像。图 1-5 是输出的图像结果。你可以使用下面的命令查看图像中的最小和最大像素值:

print int(im.min()), int(im.max())

图 1-4:灰度变换示例。三个例子中所使用函数的图像,其中虚线表示恒等变换

图 1-5:灰度变换。对图像应用图 1-4 中的函数:f(x)=255-x 对图像进行反相处理(左);f(x)=(100/255)x+100 对图像进行变换(中);f(x)=255(x/255)2 对图像做二次变换(右)

如果试着对上面例子查看最小值和最大值,可以得到下面的输出结果:

2 255
0 253
100 200
0 255

array() 变换的相反操作可以使用 PIL 的 fromarray() 函数完成:

pil_im = Image.fromarray(im)

如果你通过一些操作将“uint8”数据类型转换为其他数据类型,比如之前例子中的 im3 或者 im4,那么在创建 PIL 图像之前,需要将数据类型转换回来:

pil_im = Image.fromarray(uint8(im))

如果你并不十分确定输入数据的类型,安全起见,应该先转换回来。注意,NumPy 总是将数组数据类型转换成能够表示数据的“最低”数据类型。对浮点数做乘积或除法操作会使整数类型的数组变成浮点类型。

1.3.3 图像缩放

NumPy 的数组对象是我们处理图像和数据的主要工具。想要对图像进行缩放处理没有现成简单的方法。我们可以使用之前 PIL 对图像对象转换的操作,写一个简单的用于图像缩放的函数。把下面的函数添加到 imtool.py 文件里:

def imresize(im,sz):
  """ 使用PIL 对象重新定义图像数组的大小"""
  pil_im = Image.fromarray(uint8(im))

  return array(pil_im.resize(sz))

我们将会在接下来的内容中使用这个函数。

1.3.4 直方图均衡化

图像灰度变换中一个非常有用的例子就是直方图均衡化。直方图均衡化是指将一幅图像的灰度直方图变平,使变换后的图像中每个灰度值的分布概率都相同。在对图像做进一步处理之前,直方图均衡化通常是对图像灰度值进行归一化的一个非常好的方法,并且可以增强图像的对比度。

在这种情况下,直方图均衡化的变换函数是图像中像素值的累积分布函数(cumulative distribution function,简写为 cdf,将像素值的范围映射到目标范围的归一化操作)。

下面的函数是直方图均衡化的具体实现。将这个函数添加到 imtool.py 里:

def histeq(im,nbr_bins=256):
  """ 对一幅灰度图像进行直方图均衡化"""

  # 计算图像的直方图
  imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
  cdf = imhist.cumsum() # cumulative distribution function
  cdf = 255 * cdf / cdf[-1] # 归一化

  # 使用累积分布函数的线性插值,计算新的像素值
  im2 = interp(im.flatten(),bins[:-1],cdf)

  return im2.reshape(im.shape), cdf

该函数有两个输入参数,一个是灰度图像,一个是直方图中使用小区间的数目。函数返回直方图均衡化后的图像,以及用来做像素值映射的累积分布函数。注意,函数中使用到累积分布函数的最后一个元素(下标为 -1),目的是将其归一化到 0...1 范围。你可以像下面这样使用该函数:

from PIL import Image
from numpy import *

im = array(Image.open('AquaTermi_lowcontrast.jpg').convert('L'))
im2,cdf = imtools.histeq(im)

图 1-6 和图 1-7 为上面直方图均衡化例子的结果。上面一行显示的分别是直方图均衡化之前和之后的灰度直方图,以及累积概率分布函数映射图像。可以看到,直方图均衡化后图像的对比度增强了,原先图像灰色区域的细节变得清晰。

图 1-6:直方图均衡化示例。左侧为原始图像和直方图,中间图为灰度变换函数,右侧为直方图均衡化后的图像和相应直方图

图 1-7:直方图均衡化示例。左侧为原始图像和直方图,中间图为灰度变换函数,右侧为直方图均衡化后的图像和相应直方图

1.3.5 图像平均

图像平均操作是减少图像噪声的一种简单方式,通常用于艺术特效。我们可以简单地从图像列表中计算出一幅平均图像。假设所有的图像具有相同的大小,我们可以将这些图像简单地相加,然后除以图像的数目,来计算平均图像。下面的函数可以用于计算平均图像,将其添加到 imtool.py 文件里:

def compute_average(imlist):
  """ 计算图像列表的平均图像"""

  # 打开第一幅图像,将其存储在浮点型数组中
  averageim = array(Image.open(imlist[0]), 'f')

  for imname in imlist[1:]:
    try:
      averageim += array(Image.open(imname))
    except:
      print imname + '...skipped'
  averageim /= len(imlist)

  # 返回uint8 类型的平均图像
  return array(averageim, 'uint8')

该函数包括一些基本的异常处理技巧,可以自动跳过不能打开的图像。我们还可以使用 mean() 函数计算平均图像。mean() 函数需要将所有的图像堆积到一个数组中;也就是说,如果有很多图像,该处理方式需要占用很多内存。我们将会在下一节中使用该函数。

1.3.6 图像的主成分分析(PCA)

PCA(Principal Component Analysis,主成分分析)是一个非常有用的降维技巧。它可以在使用尽可能少维数的前提下,尽量多地保持训练数据的信息,在此意义上是一个最佳技巧。即使是一幅 100×100 像素的小灰度图像,也有 10 000 维,可以看成 10 000 维空间中的一个点。一兆像素的图像具有百万维。由于图像具有很高的维数,在许多计算机视觉应用中,我们经常使用降维操作。PCA 产生的投影矩阵可以被视为将原始坐标变换到现有的坐标系,坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行 PCA 变换,图像需要转换成一维向量表示。我们可以使用 NumPy 类库中的 flatten() 方法进行变换。

将变平的图像堆积起来,我们可以得到一个矩阵,矩阵的一行表示一幅图像。在计算主方向之前,所有的行图像按照平均图像进行了中心化。我们通常使用 SVD(Singular Value Decomposition,奇异值分解)方法来计算主成分;但当矩阵的维数很大时,SVD 的计算非常慢,所以此时通常不使用 SVD 分解。下面就是 PCA 操作的代码:

from PIL import Image
from numpy import *

def pca(X):
  """ 主成分分析:
    输入:矩阵X ,其中该矩阵中存储训练数据,每一行为一条训练数据
    返回:投影矩阵(按照维度的重要性排序)、方差和均值"""

  # 获取维数
  num_data,dim = X.shape

  # 数据中心化
  mean_X = X.mean(axis=0)
  X = X - mean_X

if dim>num_data:
  # PCA- 使用紧致技巧
  M = dot(X,X.T) # 协方差矩阵
  e,EV = linalg.eigh(M) # 特征值和特征向量
  tmp = dot(X.T,EV).T # 这就是紧致技巧
  V = tmp[::-1] # 由于最后的特征向量是我们所需要的,所以需要将其逆转
  S = sqrt(e)[::-1] # 由于特征值是按照递增顺序排列的,所以需要将其逆转
  for i in range(V.shape[1]):
    V[:,i] /= S
else:
  # PCA- 使用SVD 方法
  U,S,V = linalg.svd(X)
  V = V[:num_data] # 仅仅返回前nun_data 维的数据才合理

# 返回投影矩阵、方差和均值
return V,S,mean_X

该函数首先通过减去每一维的均值将数据中心化,然后计算协方差矩阵对应最大特征值的特征向量,此时可以使用简明的技巧或者 SVD 分解。这里我们使用了 range() 函数,该函数的输入参数为一个整数 n,函数返回整数 0...(n-1) 的一个列表。你也可以使用 arange() 函数来返回一个数组,或者使用 xrange() 函数返回一个产生器(可能会提升速度)。我们在本书中贯穿使用 range() 函数。

如果数据个数小于向量的维数,我们不用 SVD 分解,而是计算维数更小的协方差矩阵 XXT 的特征向量。通过仅计算对应前 kk 是降维后的维数)最大特征值的特征向量,可以使上面的 PCA 操作更快。由于篇幅所限,有兴趣的读者可以自行探索。矩阵 V 的每行向量都是正交的,并且包含了训练数据方差依次减少的坐标方向。

我们接下来对字体图像进行 PCA 变换。fontimages.zip 文件包含采用不同字体的字符 a 的缩略图。所有的 2359 种字体可以免费下载 2。假定这些图像的名称保存在列表 imlist 中,跟之前的代码一起保存传在 pca.py 文件中,我们可以使用下面的脚本计算图像的主成分:

2免费字体图像库由 Martin Solli 收集并上传(http://webstaff.itn.liu.se/~marso/)。

from PIL import Image
from numpy import *
from pylab import *
import pca

im = array(Image.open(imlist[0])) # 打开一幅图像,获取其大小
m,n = im.shape[0:2] # 获取图像的大小
imnbr = len(imlist) # 获取图像的数目

# 创建矩阵,保存所有压平后的图像数据
immatrix = array([array(Image.open(im)).flatten()
               for im in imlist],'f')

# 执行 PCA 操作
V,S,immean = pca.pca(immatrix)

# 显示一些图像(均值图像和前 7 个模式)
figure()
gray()
subplot(2,4,1)
imshow(immean.reshape(m,n))
for i in range(7):
  subplot(2,4,i+2)
  imshow(V[i].reshape(m,n))

show()

注意,图像需要从一维表示重新转换成二维图像;可以使用 reshape() 函数。如图 1-8 所示,运行该例子会在一个绘图窗口中显示 8 个图像。这里我们使用了 PyLab 库的 subplot() 函数在一个窗口中放置多个图像。

图 1-8:平均图像(左上)和前 7 个模式(具有最大方差的方向模式)

1.3.7 使用pickle模块

如果想要保存一些结果或者数据以方便后续使用,Python 中的 pickle 模块非常有用。pickle 模块可以接受几乎所有的 Python 对象,并且将其转换成字符串表示,该过程叫做封装(pickling)。从字符串表示中重构该对象,称为拆封(unpickling)。这些字符串表示可以方便地存储和传输。

我们来看一个例子。假设想要保存上一节字体图像的平均图像和主成分,可以这样来完成:

# 保存均值和主成分数据
f = open('font_pca_modes.pkl', 'wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()

在上述例子中,许多对象可以保存到同一个文件中。pickle 模块中有很多不同的协议可以生成 .pkl 文件;如果不确定的话,最好以二进制文件的形式读取和写入。在其他 Python 会话中载入数据,只需要如下使用 load() 方法:

# 载入均值和主成分数据
f = open('font_pca_modes.pkl', 'rb')
immean = pickle.load(f)
V = pickle.load(f)
f.close()

注意,载入对象的顺序必须和先前保存的一样。Python 中有个用 C 语言写的优化版本,叫做 cpickle 模块,该模块和标准 pickle 模块完全兼容。关于 pickle 模块的更多内容,参见 pickle 模块文档页 http://docs.python.org/library/pickle.html

在本书接下来的章节中,我们将使用 with 语句处理文件的读写操作。这是 Python 2.5 引入的思想,可以自动打开和关闭文件(即使在文件打开时发生错误)。下面的例子使用 with() 来实现保存和载入操作:

# 打开文件并保存
with open('font_pca_modes.pkl', 'wb') as f:
  pickle.dump(immean,f)
  pickle.dump(V,f)

# 打开文件并载入
with open('font_pca_modes.pkl', 'rb') as f:
  immean = pickle.load(f)
  V = pickle.load(f)

上面的例子乍看起来可能很奇怪,但 with() 确实是个很有用的思想。如果你不喜欢它,可以使用之前的 open 和 close函数。

作为 pickle 的一种替代方式,NumPy 具有读写文本文件的简单函数。如果数据中不包含复杂的数据结构,比如在一幅图像上点击的点列表,NumPy 的读写函数会很有用。保存一个数组 x 到文件中,可以使用:

savetxt('test.txt',x,'%i')

最后一个参数表示应该使用整数格式。类似地,读取可以使用:

x = loadtxt('test.txt')

你可以从在线文档 http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html 了解更多内容。

最后,NumPy 有专门用于保存和载入数组的函数。你可以在上面的在线文档里查看关于 save() 和 load() 的更多内容。

1.4 SciPy

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

1.4.1 图像模糊

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

Iσ = I*

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

G_\sigma=\frac{1}{2\pi\sigma}e^{-(x^2+y^2)/2\sigma^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|\boldsymbol{\nabla I}\right|=\sqrt{{\boldsymbol{I}_x}^2+{\boldsymbol{I}_y}^2}

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

α=arctan2(IyIx)

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

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

Ix=I*Dx 和 Iy=I*Dy

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

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

或者 Sobel 滤波器:

D_x=\begin{vmatrix}-1&0&1\\-2&0&2\\ -1&0&1\end{vmatrix} 和 D_y=\begin{vmatrix}-1&-2&-1\\0&0&0\\ 1&2&1\end{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<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(\boldsymbol{I})=\int\left|\nabla\boldsymbol{I}\right|\text{dx}            (1.1)

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

J(\boldsymbol{I})=\sum_{\text{x}}\left|\nabla\boldsymbol{I}\right|

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

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

\min_U\left|\left|\boldsymbol{I}-\boldsymbol{U}\right|\right|^2+2\lambda J(\boldsymbol{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,

Python中的图像处理

阅读数 3145

python 图像处理

阅读数 3160