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

1 PIL

PIL（Python Imagine Library Python，图像处理类库）提供了通用的图像处理功能以及大量有用的基本图像操作，比如图像缩放、裁剪、旋转、颜色转换等。

如下介绍一个简单的例子来介绍PIL函数：

from PIL import Image

pil_im = Image.open('1.jpg')
pil_im.show()

pil_im=pil_im.convert('L')
pil_im.show()

在这段代码中，利用PIL中的函数从图像格式的文件中读取数据然后写入最常见的图像格式文件中。PIL中最重要的模块就是Image，可以使用open()函数读取一幅图像，返回值pil_im是一个PIL图像对象，使用convert()来实现颜色的转换。运行结果如下图所示。 

1.1 转换图像格式

通过save()方法，PIL可以将图像保存成多种格式的文件，即保存图像到具有指定文件名的文件。

1.2 创建缩略图

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

from PIL import Image

pil_im = Image.open('2.png')
pil_im.thumbnail((128,128))
pil_im.show()

 

1.3 复制和粘贴图像区域

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

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

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

from PIL import Image
from pylab import *

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

box = (pil_im.size[0]//4,pil_im.size[1]//4,pil_im.size[0]*3//4,pil_im.size[1]*3//4)

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

pil_im.show()

1.4 调整尺寸和旋转

要调整一幅图像的尺寸，可以调用resize()方法，其类似于创建一幅缩略图。该方法的参数是一个元组，用来指定新图像的大小

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

import matplotlib.pyplot as plt
from PIL import Image
from pylab import *

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

img1 = pil_im.resize((128, 128))
img2 = pil_im.rotate(180)

plt.figure("Image")
ax1=subplot(1,2,1)
ax2=subplot(1,2,2)
plt.sca(ax1)
plt.imshow(img1)
plt.sca(ax2)
plt.imshow(img2)
plt.show()

2 Matplotlib

在Python中，Matplotlib是一个很好地类库，在处理数学运算、绘制图表，或者在图像上绘制点、直线和曲线时，具有比PIL更强大的功能。

调用方式为：

plot(x, y, [fmt], data=None, **kwargs)

点和线的坐标由参数x，y提供。可选参数fmt是一个快捷字符串，用于定义颜色、标记符合和线条形状，例如：

基本的颜色格式命令

b	蓝色
g	绿色
r	红色
c	青色
m	品红
y	黄色
k	黑色
w	白色

基本线型格式命令

-	实线
--	虚线
：	点线

 基本绘制标记格式命令

-	点
o	圆圈
s	正方形
*	星形
+	加号
*	叉号

2.1 绘制图像、点和线

from PIL import Image
import matplotlib.pyplot as plt
from pylab import *

im = array(Image.open('1.jpg'))

imshow(im)

x = [100, 100, 400, 400]
y = [200, 300, 200, 300]

plot(x, y, 'r*')
plot(x[:2], y[:2])
axis('off')
title('Plotting: "1.jpg"')

show()

在上述代码中，在x和y列表中给定点的x坐标和y坐标上绘制出红色星状点标记，然后在两个列表表示的前两个点之间绘制一条线段（默认为蓝色）。在PyLab中，约定图像的左上角为坐标原点，图像的坐标轴是一个很有用的调试工具，若不想显示坐标轴可以使用（axis（‘off’））命令 

2.2 图像轮廓和直方图

下边来看两个特别的绘图实例：图像的轮廓和直方图。

绘制轮廓需要对每个坐标[x, y]的像素值施加一个阈值，所以首先需要将图像灰度化。直方图用来表征该图像像素的分布情况。用一定数目的小区间来指定表征像素值的范围，每个小区间会得到落入该小区间表示的像素数目。该图像的直方图可以使用hist()函数绘制。hist()函数的第二个参数指定小区间的数目。需要注意的是，hist()只接受一维数组作为输入，所以在绘制直方图之前，必须先对图像进行压平处理。flatten()方法将任意数组按照行优先准则转换成一维数组。

from PIL import Image
import matplotlib.pyplot as plt
from pylab import *

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

figure()

gray()

contour(im, origin = 'image')

axis('equal')
axis('off')
figure()

hist(im.flatten(), 128)

show()

所绘制出的轮廓图和直方图如下： 

 

2.3 交互式标注

有时用户和某些应用交互，例如在一幅图像中标记一些点，或者标记一些训练数据。PyLab中的ginput()函数就可以实现交互式标注。

from PIL import Image
import matplotlib.pyplot as plt
from pylab import *

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

输出图像后等待用户在绘图窗口的图像区域点击三次，程序将这些点击的坐标[x, y]自动保存在x列表里

3 NumPy

NumPy是非常有名的Python科学计算工具包，其中包含了大量有用的思想，比如数组对象以及线性代数函数。数组对象可以完成比如矩阵乘积、转置、解方程系统、向量乘积和归一化，这为图像变形、对变化进行建模、图像分类、图像聚类等提供了基础。

3.1 图像数组表示

NumPy中的数组对象是多维的，可以用来表示向量、矩阵和图像。一个数组对象很像一个列表，但是数组中的所有元素必须具有相同的数据类型。除非创建数组对象时指定数据类型，否则数据类型会按照数据的类型自动确定。

例如：

from PIL import Image
import matplotlib.pyplot as plt
from pylab import *

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

控制台输出结果为：

每行的第一个元组表示图像数组的大小（行、列、颜色通道），紧接着的字符串表示数组元素的数据类型。因为图像通常被编码成无符号八位整数，所以在第一种情况下载入图像并将其转换到数组中，数组的数据类型为“uint8”。在第二种情况，对象进行灰度化处理，并且在创建数组时使用额外的参数“f”，该参数将数据类型转换为浮点型。

 数组中的元素可以使用下标访问。位于坐标i，j，以及颜色通道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])	# 倒数第2行

若仅使用一个下标，则该下标为行下标，负数切片表示从最后一个元素逆向计数。

3.2 灰度变换

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

from PIL import Image
import matplotlib.pyplot as plt
from pylab import *

im = array(Image.open('1.jpg').convert('L'))
im2 = 255-im
im3 = (100.0/255)*im2+100
im4 = 255.0*(im/255.0)**2
plt.subplot(131)
plt.imshow(im2)
plt.subplot(132)
plt.imshow(im3)
plt.subplot(133)
plt.imshow(im4)
plt.show()

 运行后的结果如图所示：

在上述实验中，第二个例子将灰度图像进行反相处理；第三个例子将图像的像素值变换到100…200区间；第四个例子对图像使用二次函数变换，使较暗的像素变得更小。 

可以使用print(int(im.min()), int(im.max()))命令查看图像中的最小和最大像素值

也可以使用fromarry()函数，

im = Image.fromarray(im)

3.3 图像缩放

NumPy的数组对象使我们处理图像和数据的主要工具。想要对图像进行缩放处理没有现成简单的方法。可以使用之前PIL对图像对象转换的操作来完成。

from PIL import Image
from numpy import *
from pylab import *
import numpy as np

im = array(Image.open('jimei.jpg'))

pil_im = Image.fromarray(im)

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

imshow(img1)

show()

运行结果： 

3.4 直方图均衡化

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

下边是直方图均衡化的具体实现：

def histeq(im,nbr_bins=256):
    imhist,bins=histogram(im.flatten(),nbr_bins,normed=True)
    cdf = imhist.cumsum()
    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 *
from pylab import *
import numpy as np
import scipy
import imtool

im=array(Image.open('1.jpg').convert('L'))
#原图
subplot(234)
imshow(im)
#原图直方图
subplot(231)
hist(im.flatten(),256)
im2,cdf,bins = imtool.histeq(im)
#函数
subplot(232)
plot(bins[:256],cdf)
#变换后
subplot(236)
imshow(im2)
#变换后的直方图
subplot(233)
hist(im2.flatten(),256)
show()

运行结果： 

 如图为直方图均衡化后的结果，上边为直方图均衡化之前和之后的灰度直方图，以及累计概率分布函数映射图像。可以看到，直方图均衡化之后图像的对比度增强了，原先图像灰色区域的细节变得清晰

实验问题：用imshow绘制图像时，不加gray()函数，出来的灰度图是彩色的

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(f"{imname}+...skipped")
        averageim/=len(imlist)
    return array(averageim,'unit8')

除此之外，还可以使用mean()函数计算平均图像。mean()需要将所有图像堆积到一个数组中，会占很大的内存。

3.6 图像的主成分分析（PCA）

PCA是一个非常有用的降维技巧，可以在使用尽可能少维数的前提下，尽量多的保持训练数据的信息，其产生的投影矩阵可以被视为将原始坐标变换成现有的坐标系，坐标系中的各个坐标按照重要性递减排列。

为了对图像数据进行PCA变换，图像需要转换成一维向量表示，可以使用NumPy类库中的flatten()方法进行变换。PCA操作函数如下：

def pca(X):
    num_data,dim=X.shape
    mean_X=X.mean(axis=0)
    X=x-mean_X
    
    if dim>num_data:
        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:
        U,S,V=linalg.svd(X)
        V=V[:num_data]
    return V,S,mean_X

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

from PIL import Image
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')
V,S,immean = pca.pca(immatrix)

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

show()

3.7 使用pickle模块

pickle模块接受几乎所有的Python对象，并且将其转换成字符串表示。该过程叫做封装。从字符串表示中重构该对象，称为拆封。这些字符串表示可以方便地存储和传输。

使用如下代码保存平均图像和主成分：

import pickle

f = open('font_pca_modes.pkl','wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()

pickle中有很多不同的协议可以生成.pkl文件，如果不确定的话最好以二进制文件的形式读取和写入。

打开文件并保存：

import pickle

f = open('font_pca_modes.pkl','rb')
immean = pickle.load(f)
V = pickle.load(f)
f.close()

可以用with open 来简化代码，且不用在最后加上close（）

作为pickle的一种替代方式，NumPy具有读写文本文件的简单函数，即：

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

4 SciPy模块

SciPy是建立在NumPy基础上，用于数值运算的开源工具包，可以实现数值积分、优化、统计、信号处理以及图像处理功能。

4.1 图像模糊

SciPy有用来做滤波操作的scipy.ndimage.filters模块。该模块使用快速一维分离的方式来计算卷积。

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from scipy.ndimage import filters
im=array(Image.open('1.jpg').convert('L'))
im1 = filters.gaussian_filter(im,2)
im2= filters.gaussian_filter(im,5)
im3 = filters.gaussian_filter(im,10)
subplot(141)
gray()
imshow(im)
subplot(142)
gray()
imshow(im1)
subplot(143)
gray()
imshow(im2)
subplot(144)
gray()
imshow(im3)
show()

运行结果：

 随着σ的增加，图像的模糊程度就越大，细节丢失也越多。

若要对一幅彩色图像进行模糊，只需要对每一个颜色通道进行高斯模糊：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from scipy.ndimage import filters
im=array(Image.open('1.jpg'))
im2 = zeros(im.shape)
for i in range(3):
    im2[:,:,i]=filters.gaussian_filter(im[:,:,i],5)
im2=uint8(im2)
subplot(121)
imshow(im)
axis('off')
subplot(122)
imshow(im2)
axis('off')
show()

运行结果：

 

4.2 图像导数

图像的强度变化可以用灰度图像I的x和y方向的导数Ix和Iy进行描述，图像的导数大多通过卷积来实现。通常选用Prewitt滤波器和Sobel滤波器。这些导数滤波器可以使用scipy.ndimage.filters模块的标准卷积操作来简单实现：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from scipy.ndimage import filters
im=array(Image.open('1.jpg').convert('L'))

imx=zeros(im.shape)
filters.sobel(im,1,imx)

imy =zeros(im.shape)
filters.sobel(im,0,imy)
magnitude=sqrt(imx**2+imy**2)

subplot(221)
gray()
axis('off')
imshow(im)
subplot(222)
gray()
axis('off')
imshow(imx)
subplot(223)
gray()
axis('off')
imshow(imy)
subplot(224)
gray()
axis('off')
imshow(magnitude)
show()

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

上述计算图像导数的方法存在一定的缺陷：在该方法中，滤波器的尺度需要随着图像分辨率的变化而变化。为了在处理图像噪声方面更加稳健，以及在任意尺度上计算导数，可以使用高斯导数滤波器。用于模糊的filters.gaussian)filter()函数可以接受额外的参数用来计算高斯导数，可以按照下边的方式处理：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from scipy.ndimage import filters
im=array(Image.open('1.jpg').convert('L'))

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),(0,1),imx)
magnitude=sqrt(imx**2+imy**2)

subplot(221)
gray()
axis('off')
imshow(im)
subplot(222)
gray()
axis('off')
imshow(imx)
subplot(223)
gray()
axis('off')
imshow(imy)
subplot(224)
gray()
axis('off')
imshow(magnitude)
show()

该函数的第三个参数指定对每个方向计算哪种类型的导数，第二个参数为使用的标准差。下图显示了在标准差为5的高斯导数下的模糊，相较Sobel导数下模糊程度要更大。

4.3 形态学：对象计数

形态学是度量和分析基本形状的图像处理方法的基本框架与集合。形态学通常用于处理二值图像，但是也能够用于灰度图像。

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

在一个二值图像中，计算该图像中的对象个数可以通过下边代码实现：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from scipy.ndimage import filters
from scipy.ndimage import measurements,morphology
im=array(Image.open('1.jpg').convert('L'))
im=1*(im<128)

labels,nbr_object=measurements.label(im)
print("Number of objects:",nbr_object)

运行结果： 

 

首先载入图片，通过阈值化方式来确保该图像是二值图像。通过和1相乘将布尔数组转换成二进制表示，然后使用label()函数寻找单个的物体，并且按照它们属于哪个对象将整数标签给像素赋值。若是一些对象之间有一些小的连接可以使用开操作将其移出：

im_open = morphology.binary_opening(im, ones((9, 5)), iterations = 2)

label_open, nbr_objects_open = measurements.label(im_open)
print("Number of objects:", nbr_objects)

binart_opening()函数的第二个参数指定一个数组结构元素。该数组表示以一个像素为中心时，使用哪些相邻元素。在这种情况下，我们在y方向上使用9个像素，在x方向上使用5个像素。binary_closing()函数实现相反的操作

4.4 一些有用的Scipy模块

Scipy中包含一些用于输入和输出的实用模块，例如io和misc。

4.4.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。

4.4.2 以图像形式保存数组

因为需要对图像进行操作，并且需要使用数组对象来做运算，所以将数组直接保存为图像文件非常有用。

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

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

scipy.misc模块同样包含了Lena测试：

lena = scipy.misc.lena()

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

5 图像去噪

图像去噪是在去除图像噪声的同时，尽可能地保留图像细节和结构的处理技术。在接下来的实验中使用ROF去噪模型。ROF去噪模型具有 很好的性质：是处理后的图像更平滑，同时保持图像边缘和结构信息，其理论在本质上是使去噪后的图像像素值“平坦”变换，但是在图像区域的边缘上，允许去噪后的图像像素值“跳跃”变化。

实现ROF模型的代码如下：

def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):

    m,n=im.shape

    U=U_init
    Px=im
    Py=im
    error = 1
    while(error>tolerance):
        Uold=U

        GradUx=roll(U,-1,axis=1)-U
        GradUy=roll(U,-1,axis=0)-U

        PxNew = Px+(tau/tv_weight)*GradUx
        PyNew = Py+(tau/tv_weight)*GradUy
        NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))

        Px=PxNew/NormNew
        Py=PyNew/NormNew

        RxPx=roll(Px,1,axis=1)
        RyPy=roll(Py,1,axis=0)

        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()函数可以衡量两个数组之间的差异。

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from numpy import random
from scipy.ndimage import filters
import imtool


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=imtool.denoise(im,im)
G=filters.gaussian_filter(im,10)

subplot(131)
gray()
axis('off')
imshow(im)

subplot(132)
gray()
axis('off')
imshow(U)

subplot(133)
gray()
axis('off')
imshow(G)
show()

原始图像和去噪结果如下所示：

对实际图像去噪：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import scipy
from numpy import random
from scipy.ndimage import filters
import imtool


im=array(Image.open('1.jpg').convert('L'))
U,T=imtool.denoise(im,im)

subplot(121)
gray()
axis('off')
imshow(im)

subplot(122)
gray()
axis('off')
imshow(U)

show()

运行结果：

可以明显看出，ROF模型去噪后的图像保留了边缘和图像的结构信息