2. SciPy基础

SciPy以NumPy为基础，与NumPy一样，SciPy有着稳定，成熟，且应用广泛的数值运算库。方便、易于使用、专为科学和工程设计的python工具包，它包括了统计、优化、整合以及线性代数模块、傅里叶变换、信号和图像图例，常微分方差的求解等。

SciPy官方教程很完善，可自行参考学习。

2.1 scipy.misc模块

misc模块，提供一些基本的图像相关的读写函数，可以很轻松的读取本地图像文件到Python程序里，也可将数据输出到图像文件。misc模块自带一些灰度图像ascent和彩色的face图，可以scipy.misc.ascent()直接获取爬楼梯ascent图数据到Python程序里，用scipy.misc.face()获取一副 raccoon浣熊face图，这两个函数的返回值都是Numpy的ndarray数组。

face图像是个彩色图像，其数据是个三维数组，是个1024x768的图像，而图像中每个像素的值又是一个数组，分别对应该像素颜色的红、绿、蓝分量。ascent图像是个灰度图像，其数据是个二维数组，分别对应图像中每个像素的灰度值。

使用ascent图

import matplotlib.pyplot as plt

import numpy as np

import scipy.misc as sm

ascent = sm.ascent()

print ascent.shape

plt.figure()

plt.imshow(ascent)

plt.show()

程序执行结果：

(512, 512) # ascent.face

使用face图

import matplotlib.pyplot as plt

import numpy as np

import scipy.misc as sm

face = sm.face()

print face.shape

plt.figure()

plt.imshow(face)

plt.show()

程序执行结果：

(768, 1024, 3) # face.shape

2.2 图像数据属性

scipy.misc.face()返回的是浣熊的脸的彩色图像，是一个Numpy的三维数组。

import matplotlib.pyplot as plt

import numpy as np

import scipy.misc as sm

x = np.linspace(1, 5, 5)

print type(x), x

print x.shape, x.ndim, x.dtype

face = sm.face()

print type(face)

print face.shape, face.ndim, face.dtype

程序执行结果：

[ 1. 2. 3. 4. 5.]

(5,) 1 float64

(768, 1024, 3) 3 uint8

从程序执行结果可以看出，尽管是scipy模块的函数，创建的face是numpy的数组，说明scipy里多数使用的是numpy对象数据。

2.3 NumPy数组操作函数

2.3.1 unique函数

unique函数可以返回NumPy数组对象里的唯一值，数组由那些数据组成。意思有点像Python的set。

#coding:utf-8

import numpy as np

import scipy.misc as sm

#自己实现unique函数

def my_unique(x):

ret = list()

if x.ndim > 1:

x = x.flatten(order='C')

for v in x:

if v not in ret:

ret.append(v)

ret = np.sort(ret)

return ret

y = np.array([1, 3, 2, 1, 4, 2])

print y, "# y"

print np.unique(y), "# unique(y)"

print my_unique(y), "# my_unique(y)"

face = sm.face()

print np.unique(face),"# unique(face)"

print my_unique(face),"# my_unique(face)"

程序执行结果：

[1 3 2 1 4 2] # y

[1 2 3 4] # unique(y)

[1 2 3 4] # my_unique(y)

[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89

90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125

126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143

144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179

180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197

198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215

216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233

234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251

252 253 254 255] # unique(face)

[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89

90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125

126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143

144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179

180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197

198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215

216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233

234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251

252 253 254 255] # my_unique(face)

从np.unique(face)打印的结果可以看出的确灰度图基础数据是由0～255来描述的。

2.3.2 bincount函数

numpy的bincount函数可以统计出数组里的从0到数组最大值n共n+1个自然数出现的次数(数字图像存储的数据都是整数)。

基本原理是先找出数组里的最大值，统计0～最大值间的所有值出现的次数，函数返回共最大值加一个数据组成的数组。

import numpy as np

import scipy.misc as sm

def mybincount(x):

t = np.max(x) +1

i = 0

ret = np.zeros(t, dtype = np.uint8)

for v in x:

ret[v] += 1

return ret

ascent = sm.ascent()

print "max of array", np.max(ascent[0]), ', has', len(ascent[0])

b = np.bincount(ascent[0])

print "return array's length", len(b)

print "result", b

c = mybincount(ascent[0])

if (b == c).all():

print "mybincount realized bincount"

print c

程序执行结果

max of array 118 , has 512

return array's length 119

result [ 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 2 0 0 1 0 1 1 4 1

0 2 2 0 2 1 1 2 3 2 0 8 2 3 2 2 3 2 6 4 5 0 5 2 3

2 20 8 3 2 1 3 5 2 4 11 9 3 1 1 3 3 1 1 2 10 3 3 2 2

0 1 1 1 2 2 1 24 15 3 28 9 27 3 2 27 7 26 2 1 30 2 26 2 1

6 2 0 0 0 0 1 0 0 0 0 0 0 0 0 22 6 50 3]

mybincount realized bincount

[ 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 2 0 0 1 0 1 1 4 1

0 2 2 0 2 1 1 2 3 2 0 8 2 3 2 2 3 2 6 4 5 0 5 2 3

2 20 8 3 2 1 3 5 2 4 11 9 3 1 1 3 3 1 1 2 10 3 3 2 2

0 1 1 1 2 2 1 24 15 3 28 9 27 3 2 27 7 26 2 1 30 2 26 2 1

6 2 0 0 0 0 1 0 0 0 0 0 0 0 0 22 6 50 3]

2.3.3 fromfunction函数

函数的作用是依据第一个形参(是个函数或lambda表达式)来产生一个数组。

import numpy as np

def fn(x, y):

print x, "# x"

print y, "# y"

return (x + y) * 3

t = np.fromfunction(fn, (3, 4), dtype=np.uint8)

print "t = ", t

程序的执行结果：

[[0 0 0 0]

[1 1 1 1]

[2 2 2 2]]# x

[[0 1 2 3]

[0 1 2 3]

[0 1 2 3]]# y

t = [[ 0 3 6 9]

[ 3 6 9 12]

[ 6 9 12 15]]

从结果可以看出，函数先产生shape为(3, 4)的数组x和y，x数组里的元素列变化从0～2,y数组的行元素值变化从0～3和shape指定值有关。结果是return语句的表达式的值，本质上是x和y数组的表达式值。

2.3.4 put函数

put函数可以将数组的某些值替换成其他值，函数原型如下：

numpy.put(a, ind, v, mode='raise')

a是要修改的函数，ind是a里要修改的数值的在数组a里的位置信息的列表，而v是要被替换成的数据列表，原理是先将a“展开”成一维a' = a.flatten()，然后在对a'的ind位置处的数据修改成v里的某个值。

import numpy as np

a = np.arange(12).reshape([3, 4])

def myput(x, ind, v):

s = x.shape

if x.ndim > 1:

x = x.flatten()

i = 0

for j in ind:

x[j] = v[i]

i += 1

x = x.reshape(s)

return x

print a, "# a"

ret = myput(a, [1, 3], [100, 300])

print ret, "# a of put"

ret = np.put(a, [1, 3], [100, 300])

print a, "# a of myput"

程序执行结果：

[[ 0 1 2 3]

[ 4 5 6 7]

[ 8 9 10 11]] # a

[[ 0 100 2 300]

[ 4 5 6 7]

[ 8 9 10 11]] # a of put

[[ 0 100 2 300]

[ 4 5 6 7]

[ 8 9 10 11]] # a of myput

numpy.put函数要求ind和v同等长度，如果两者长度不同可以使用mode参数设置为'clip'。

2.3.5 putmask函数

这个函数很有用，其定义原型：

numpy.putmask(a, mask, values)

a是数组，mask是一个关系表达式，values是修改a用到的数据值。

1). putmask基本使用

先看一个简单的程序。

import numpy as np

a = np.arange(1, 13).reshape([3, 4])

print a

np.putmask(a, a % 2 == 0, [0])

print a

程序执行结果：

[[ 1 2 3 4]

[ 5 6 7 8]

[ 9 10 11 12]]

[[ 1 0 3 0]

[ 5 0 7 0]

[ 9 0 11 0]]

程序的意思是将a数组里元素值为偶数的元素设置为0。

2). 替换多个值

上边的例子是满足条件mask的位置上的数据均替换成0,其实values可以是列表，列出可以替换成的数据信息，这时数组a和values的长度不一样长，那么如何替换呢？规则有点像numpy.put函数，但是values没有规定长度，可以小于数组a的程度，只需把values列表扩展成大于等于数组大的长度然后对应位置替换成对应values扩展后的数值。

import numpy as np

a = np.arange(12).reshape([3, 4])

def myputmask(x, mask, v):

s = x.shape

if x.ndim > 1:

x = x.flatten()

if len(x) > len(v):

n = len(x) / len(v) + 1

v = v * n

mask = mask.flatten()

print "x ", x.flatten()

print "mask", mask

print "v ", v

i = 0

j = 0

for tf in mask:

if tf == True:

x[j] = v[i]

i += 1

j += 1

x = x.reshape(s)

return x

print a, "# a"

ret = myputmask(a, a > 3, [100, 300])

print ret, "# a of myputmask"

ret = np.putmask(a, a > 3, [100, 300])

print a, "# a of putmask"

程序执行结果：

[[ 0 1 2 3]

[ 4 5 6 7]

[ 8 9 10 11]] # a

x [ 0 1 2 3 4 5 6 7 8 9 10 11]

mask [False False False False True True True True True True True True]

v [100, 300, 100, 300, 100, 300, 100, 300, 100, 300, 100, 300, 100, 300]

[[ 0 1 2 3]

[100 300 100 300]

[100 300 100 300]] # a of myputmask

[[ 0 1 2 3]

[100 300 100 300]

[100 300 100 300]] # a of putmask

3). 对图像face进行putmask处理，功能是将像素值小于127的替换成0,大于127的设置成255，这样实现图像“锐化”。

import numpy as np

import scipy.misc as sm

import matplotlib.pyplot as plt

ascent = sm.ascent()

a = ascent

b = ascent

c = ascent

plt.gray()

plt.imshow(a)

plt.show()

np.putmask(b, b >= 127, 255)

plt.gray()

plt.imshow(b)

plt.show()

np.putmask(c, c < 127, 0)

plt.gray()

plt.imshow(c)

plt.show()

程序执行结果：

a数组显示结果

b数组显示结果

c数组显示结果

2.3.6 meshgrid函数

输入两个一维数组，输出两个二维数组，XY平面中网格的交点就是meshgrid函数返回数据值，第一个返回值是X轴的坐标集合，第二个Y轴坐标集合，基于这两个返回值计算z，然后用matplotlib的相应函数可以绘制2D或者3D可视化数据图。

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.mlab as mlab

xi= np.linspace(-3, 3, 50)

yi= np.linspace(-2, 2, 50)

print xi

print yi

X, Y = np.meshgrid(xi, yi)

Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)

Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)

Z = 10.0 * (Z2 - Z1)

plt.figure()

CS = plt.contour(X, Y, Z)

plt.clabel(CS, inline=1, fontsize=10)

plt.title('Simplest default with labels')

plt.show()

程序执行结果：

2.3.7 numpy.mgrid获取坐标矩阵

NumPy内建的numpy.lib.index_tricks.nd_grid类的实例对象mgrid可以通过切片技术获取类似meshgrid函数取得的XY坐标平面的坐标集，返回值和meshgrid一样也有两个分布是X轴坐标集和Y轴坐标集。

这段文字比较绕，简而言之: (1)mgrid是NumPy内建实例对象，即可直接使用也可自己创建一个实例对象，具有和mgrid同样的功能。(2) 对mgrid进行切片处理可以得到两个集合：X轴坐标集和Y轴坐标集。

import numpy as np

xgrid = np.lib.index_tricks.nd_grid()

xx, yy = xgrid[1 : 3, 2 : 6]

print "# xgrid\n", xx, "\n", yy

xx, yy = np.mgrid[1 : 3, 2 : 6]

print "# mgrid\n", xx, "\n", yy

程序执行结果：

# xgrid

[[1 1 1 1]

[2 2 2 2]]

[[2 3 4 5]

[2 3 4 5]]

# mgrid

[[1 1 1 1]

[2 2 2 2]]

[[2 3 4 5]

[2 3 4 5]]

2.3.8 numpy.ogrid获取坐标向量

mgrid和ogrid都是NumPy已经创建好，内建的numpy.lib.index_tricks.nd_grid类的实例对象。区别mgrid是nd_grid的形参sparse为False，而ogrid的形参sparse为True。

import numpy as np

xgrid = np.lib.index_tricks.nd_grid(sparse=True)

xx, yy = xgrid[1 : 3, 2 : 6]

print "# xgrid\n", xx, "\n", yy

xx, yy = np.ogrid[1 : 3, 2 : 6]

print "# ogrid\n", xx, "\n", yy

xx, yy = np.mgrid[1 : 3, 2 : 6]

print "# mgrid\n", xx, "\n", yy

程序执行结果：

# xgrid

[[1]

[2]]

[[2 3 4 5]]

# ogrid

[[1]

[2]]

[[2 3 4 5]]

# mgrid

[[1 1 1 1]

[2 2 2 2]]

[[2 3 4 5]

[2 3 4 5]]

ogrid返回的是x和y坐标组成的向量，可以用matplotlib的scatter函数来画出可视化数据。

import numpy as np

import matplotlib.pyplot as plt

np.random.seed(19680801)

N = 50

x, y = np.ogrid[0:50,0:50]

print x

print y

np.random.shuffle(x)

np.random.shuffle(y)

print x

print y

colors = np.random.rand(N)

area = (30 * np.random.rand(N))**2 # 0 to 15 point radii

plt.scatter(x, y, s=area, c=colors, alpha=0.5)

plt.show()

程序执行结果：

2.4 总结

本章学习了很多的NumPy的数组处理函数，是6 NumPy的数学和统计函数衍生剧，呵呵。