初識SciPy庫

• 1. SciPy
• 2. SciPy任務子模塊
• 3. 文件輸入/輸出：
• 4. 特殊函數：
• 5. 線性代數操作：
• 6. 快速傅立葉變換：
• 7. 優化及擬合：
• • 7.1 尋找標量函數的最小值
  • 7.2 尋找標量函數的根
  • 7.3 曲線擬合
• 8. 統計和隨機數：
• • 8.1 直方圖和概率密度函數
  • 8.2 百分位數
  • 8.3 統計檢驗
• 9. 插值：
• 10. 數值積分：
• 11. 信號處理：
• 12. 圖像處理：scipy.ndimage
• • 12.1 圖像的幾何變換
  • 12.2 圖像濾波器
  • 12.3 數學形態學
  • 12.4 測量圖像
  • 12.5 圖像處理應用：

1. SciPy

---

SciPy包含許多專注于科學計算中的常見問題的工具箱。它的子模塊對應于不同的應用，比如插值、積分、優化、圖像處理、統計和特殊功能等。

SciPy可以與其他標準科學計算包相對比，比如GSL (C和C++的GNU科學計算包), 或者Matlab的工具箱。SciPy是Python中科學程序的核心程序包；這意味著有效的操作NumPy數組，因此，NumPy和SciPy可以一起工作。

在實現一個程序前，有必要確認一下需要的數據處理方式是否已經在SciPy中實現了。作為非專業程序員，科學家通常傾向于重新發明輪子，這產生了小玩具、不優化、很難分享以及不可以維護的代碼。相反，SciPy的程序是優化并且測試過的，因此應該盡可能使用。

警告 這個教程根本不是數值計算的介紹。因為列舉SciPy的不同子模塊和功能將會是非常枯燥的，相反我們將聚焦于列出一些例子，給出如何用SciPy進行科學計算的大概思路。

參考資料

SciPy文檔

2. SciPy任務子模塊

---

關鍵詞參數可選值

scipy.cluster	向量計算
scipy.constants	物理和數學常量
scipy.fftpack	傅里葉變換
scipy.integrate	積分程序
scipy.interpolate	插值
scipy.io	數據輸入和輸出
scipy.linalg	線性代數程序
scipy.ndimage	n-維圖像包
scipy.odr	正交距離回歸
scipy.optimize	優化
scipy.signal	信號處理
scipy.sparse	稀疏矩陣
scipy.spatial	空間數據結構和算法
scipy.special	一些特殊數學函數
scipy.stats	統計

他們全都依賴于NumPy, 但是大多數是彼此獨立的。

導入NumPy和SciPy的標準方式：

import numpy as np from scipy import stats # 其他的子模塊統計類

scipy的主要命名空間通常包含的函數其實是numpy（試一下scipy.cos其實是np.cos) 。這些函數的暴露只是因為歷史原因；通常沒有必要在你的代碼中使用import scipy。

3. 文件輸入/輸出：

scipy.io

載入和保存matlab文件：

from scipy import io as spioa = np.ones((3, 3)) spio.savemat('file.mat', {'a': a}) # savemat expects a dictionary data = spio.loadmat('file.mat', struct_as_record=True) data['a']

更多請見：

• 加載文本文件：numpy.loadtxt()/numpy.savetxt()
• 智能加載文本/csv文件：numpy.genfromtxt()/numpy.recfromcsv()
• 快速有效，但是針對numpy的二進制格式：numpy.save()/numpy.load()

4. 特殊函數：

特殊函數是超驗函數。scipy.special模塊的文檔字符串寫的很詳細，因此我們不會在這里列出所有的函數。常用的一些函數如下：

• 貝塞爾函數，比如scipy.special.jn() (第n個整型順序的貝塞爾函數)
• 橢圓函數 (scipy.special.ellipj() Jacobian橢圓函數, …)
• Gamma 函數: scipy.special.gamma(), 也要注意 scipy.special.gammaln() 將給出更高準確數值的 Gamma的log。
• Erf, 高斯曲線的面積：scipy.special.erf()

5. 線性代數操作：

scipy.linalg 模塊提供了標準的線性代數操作，這依賴于底層的高效實現（BLAS、LAPACK）。

• scipy.linalg.det() 函數計算方陣的行列式：

from scipy import linalg arr = np.array([[1, 2],[3, 4]]) linalg.det(arr)#求取行列式的值det arr = np.array([[3, 2],[6, 4]]) linalg.det(arr) try:linalg.det(np.ones((3, 4))) except ValueError:print("expected square matrix")

• scipy.linalg.inv() 函數計算逆方陣：

arr = np.array([[1, 2],[3, 4]]) iarr = linalg.inv(arr) iarr

判斷矩陣與矩陣的逆乘積是否為單位矩陣

np.allclose(np.dot(arr, iarr), np.eye(2))

最后計算逆奇異矩陣（行列式為0）將拋出LinAlgError :

arr = np.array([[3, 2],[6, 4]])try:linalg.inv(arr) except np.linalg.LinAlgError as err:print("singular matrix")

• 還有更多高級的操作，奇異值分解（SVD）：

arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1]) uarr, spec, vharr = linalg.svd(arr)

原始矩陣可以用svd和np.dot矩陣相乘的結果重新獲得：

sarr = np.diag(spec) svd_mat = uarr.dot(sarr).dot(vharr) np.allclose(svd_mat, arr)#判斷兩個矩陣是否相同

SVD常被用于統計和信號處理。其他標準分解 (QR, LU, Cholesky, Schur), 以及線性系統的求解器，也可以在scipy.linalg中找到。

6. 快速傅立葉變換：

---

scipy.fftpack 模塊允許計算快速傅立葉變換。例子，一個（有噪音）的信號輸入是這樣：

time_step = 0.02 period = 5. time_vec = np.arange(0, 20, time_step) sig = np.sin(2 * np.pi / period * time_vec) + \0.5 * np.random.randn(time_vec.size)import matplotlib.pyplot as plt plt.figure(figsize=(6, 5)) plt.plot(time_vec, sig, label='Original signal')

觀察者并不知道信號的頻率，只知道抽樣時間步驟的信號sig。假設信號來自真實的函數，因此傅立葉變換將是對稱的。scipy.fftpack.fftfreq() 函數將生成樣本序列，而scipy.fftpack.fft()將計算快速傅立葉變換：

from scipy import fftpack sample_freq = fftpack.fftfreq(sig.size, d=time_step) sig_fft = fftpack.fft(sig)

因為生成的冪是對稱的，尋找頻率只需要使用頻譜為正的部分：

pidxs = np.where(sample_freq > 0) freqs = sample_freq[pidxs] power = np.abs(sig_fft)[pidxs]

尋找信號頻率：

freq = freqs[power.argmax()] np.allclose(freq, 1./period) # 檢查是否找到了正確的頻率

現在高頻噪音將從傅立葉轉換過的信號移除：

sig_fft[np.abs(sample_freq) > freq] = 0

生成的過濾過的信號可以用scipy.fftpack.ifft()函數：

main_sig = fftpack.ifft(sig_fft)

查看結果：

plt.figure(figsize=(6, 5)) plt.plot(time_vec, sig) plt.plot(time_vec, main_sig.real, linewidth=3) plt.xlabel('Time [s]') plt.ylabel('Amplitude') plt.show()

7. 優化及擬合：

優化是尋找最小化或等式的數值解的問題。

scipy.optimize 模塊提供了函數最小化（標量或多維度）、曲線擬合和求根的有用算法。

from scipy import optimize

7.1 尋找標量函數的最小值

讓我們定義下面的函數：

def f(x):return x**2 + 10*np.sin(x)

繪制它：

x = np.arange(-10, 10, 0.1) plt.plot(x, f(x)) plt.show()

這個函數在-1.3附近有一個全局最小并且在3.8有一個局部最小。

找到這個函數的最小值的常用有效方式是從給定的初始點開始進行一個梯度下降。BFGS算法是這樣做的較好方式：

optimize.fmin_bfgs(f, 0)

這個方法的一個可能問題是，如果這個函數有一些局部最低點，算法可能找到這些局部最低點而不是全局最低點，這取決于初始點：

optimize.fmin_bfgs(f, 3, disp=0)

如果我們不知道全局最低點，并且使用其臨近點來作為初始點，那么我們需要付出昂貴的代價來獲得全局最優。要找到全局最優點，最簡單的算法是暴力算法，算法中會評估給定網格內的每一個點：

grid = (-10, 10, 0.1) xmin_global = optimize.brute(f, (grid,)) xmin_global

對于更大的網格，scipy.optimize.brute() 變得非常慢。scipy.optimize.anneal() 提供了一個替代的算法，使用模擬退火。對于不同類型的全局優化問題存在更多的高效算法，但是這超出了scipy的范疇。OpenOpt、IPOPT、PyGMO和PyEvolve是關于全局優化的一些有用的包。

要找出局部最低點，讓我們用scipy.optimize.fminbound將變量限制在(0,10)區間：

xmin_local = optimize.fminbound(f, 0, 10) xmin_local

7.2 尋找標量函數的根

要尋找上面函數f的根，比如f(x)=0的一個點，我們可以用比如scipy.optimize.fsolve()：

7.3 曲線擬合

假設我們有來自f的樣例數據，帶有一些噪音：

xdata = np.linspace(-10, 10, num=20) ydata = f(xdata) + np.random.randn(xdata.size)

現在，如果我們知道這些sample數據來自的函數（這個案例中是$x^2+sin(x)$）的函數形式，而不知道每個數據項的系數，那么我們可以用最小二乘曲線擬合在找到這些系數。首先，我們需要定義函數來擬合：

def f2(x, a, b):return a*x**2 + b*np.sin(x)

然后我們可以使用scipy.optimize.curve_fit()來找到a和b：

guess = [2, 2] params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess) params

現在我們找到了f的最優解和根，并且用曲線去擬合它，我們將這些結果整合在一個圖中：

plt.figure(figsize=(6, 4)) plt.scatter(xdata, ydata, label='Data') plt.plot(xdata, f2(xdata, params[0], params[1]),label='Fitted function')

注：在Scipy中包含所有最小值和尋找根的算法的統一接口：scipy.optimize.minimize()、 scipy.optimize.minimize_scalar()和 scipy.optimize.root()。他們允許通過method關鍵詞容易的比較多種算法。

你可以在scipy.optimize中找到對于多維度問題有相同功能的算法。

8. 統計和隨機數：

---

scipy.stats模塊包含統計工具和隨機過程的概率描述。在numpy.random中可以找到多個隨機數生成器。

8.1 直方圖和概率密度函數

給定隨機過程的觀察值，它們的直方圖是隨機過程的PDF（概率密度函數）的估計值：

a = np.random.normal(size=1000) bins = np.arange(-4, 5) bins histogram = np.histogram(a, bins=bins)[0] bins = 0.5*(bins[1:] + bins[:-1]) bins from scipy import statsb = stats.norm.pdf(bins) # norm 是一種分布 plt.plot(bins, histogram) plt.plot(bins, b) plt.show()

8.2 百分位數

中位數是有一半值在其上一半值在其下的值：

np.median(a)

中數也被稱為百分位數50，因為50%的觀察值在它之下：

stats.scoreatpercentile(a, 50)

同樣，我們也能計算百分位數90：

stats.scoreatpercentile(a, 90)

百分位數是CDF的估計值：累積分布函數。

8.3 統計檢驗

統計檢驗是一個決策指示器。例如，如果我們有兩組觀察值，我們假設他們來自于高斯過程，我們可以用T檢驗來決定這兩組觀察值是不是顯著不同：

a = np.random.normal(0, 1, size=100) b = np.random.normal(1, 1, size=10) stats.ttest_ind(a, b)

生成的結果由以下內容組成：

• T 統計值：一個值，符號與兩個隨機過程的差異成比例，大小與差異的程度有關。
• p 值：兩個過程相同的概率。如果它接近1，那么這兩個過程幾乎肯定是相同的。越接近于0，越可能這兩個過程有不同的平均數。

9. 插值：

---

scipy.interpolate對從實驗數據中擬合函數是非常有用的，因此，評估沒有測量過的點。這個模塊是基于netlib項目的Fortran子程序 FITPACK

假想一個接近sine函數的實驗數據：

measured_time = np.linspace(0, 1, 10) noise = (np.random.random(10)*2 - 1) * 1e-1 measures = np.sin(2 * np.pi * measured_time) + noise

scipy.interpolate.interp1d類可以建立一個線性插值函數：

from scipy.interpolate import interp1d linear_interp = interp1d(measured_time, measures)

scipy.interpolate.linear_interp實例需要評估感興趣的時間點：

computed_time = np.linspace(0, 1, 50) linear_results = linear_interp(computed_time)

通過提供可選的參數kind也可以選擇進行立方插值：

cubic_interp = interp1d(measured_time, measures, kind='cubic') cubic_results = cubic_interp(computed_time)

現在結果可以被整合為下面的Matplotlib圖片:

plt.figure(figsize=(6, 4)) plt.plot(measured_time, measures, 'o', ms=6, label='measures') plt.plot(computed_time, linear_results, label='linear interp') plt.plot(computed_time, cubic_results, label='cubic interp')

scipy.interpolate.interp2d 與scipy.interpolate.interp1d類似，但是是用于2-D數組。注意對于interp家族，計算的時間點必須在測量時間段之內。

10. 數值積分：

scipy.integrate

scipy.integrate.quad()是最常見的積分程序：

from scipy.integrate import quad res, err = quad(np.sin, 0, np.pi/2) np.allclose(res, 1)

其他的積分程序可以在fixed_quad、 quadrature、romberg中找到。

scipy.integrate 可提供了常微分公式(ODE)的特色程序。特別的，scipy.integrate.odeint() 是使用LSODA（Livermore Solver for Ordinary Differential equations with Automatic method switching for stiff and non-stiff problems）的通用積分器，更多細節請見ODEPACK Fortran 庫。

odeint解決如下形式的第一順序ODE系統：

$dy/dt=rhs(y1,y2,..,t0,...)$

作為一個介紹，讓我們解一下在初始條件下$y(t=0)=1$，這個常微分公式$dy/dt=?2y$在$t=0..4$時的值。首先，這個函數計算定義位置需要的導數：

def calc_derivative(ypos, time, counter_arr):counter_arr += 1return -2 * ypos

添加了一個額外的參數counter_arr用來說明這個函數可以在一個時間步驟被調用多次，直到收斂。計數器數組定義如下：

counter = np.zeros((1,), dtype=np.uint16)

現在計算軌跡線：

from scipy.integrate import odeint time_vec = np.linspace(0, 4, 40) yvec, info = odeint(calc_derivative, 1, time_vec,args=(counter,), full_output=True)

注意，求解器對于首個時間步驟需要更多的循環。導數答案yvec可以畫出來：

plt.figure(figsize=(4, 3)) plt.plot(time_vec, yvec) plt.xlabel('t: Time') plt.ylabel('y: Position')

11. 信號處理：

scipy.signal

• scipy.signal.detrend(): 從信號中刪除線性趨勢：

from scipy import signalt = np.linspace(0, 5, 100) x = t + np.random.normal(size=100)plt.plot(t, x, linewidth=3) plt.plot(t, signal.detrend(x), linewidth=3) plt.show()

• scipy.signal 有許多窗口函數：scipy.signal.hamming(), scipy.signal.bartlett(), scipy.signal.blackman()…
• scipy.signal 有濾鏡 (中位數濾鏡scipy.signal.medfilt(), Wienerscipy.signal.wiener())

12. 圖像處理：scipy.ndimage

---

scipy中專注于專注于圖像處理的模塊是scipy.ndimage。

from scipy import ndimage

圖像處理程序可以根據他們進行的處理來分類。

12.1 圖像的幾何變換

改變原點，解析度，…

from scipy import miscface = misc.face(gray=True) shifted_face = ndimage.shift(face, (50, 50)) shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest') rotated_face = ndimage.rotate(face, 30) cropped_face = face[50:-50, 50:-50] zoomed_face = ndimage.zoom(face, 2) zoomed_face.shape

plt.figure(figsize=(15, 3)) plt.subplot(151) plt.imshow(shifted_face, cmap=plt.cm.gray) plt.axis('off')plt.subplot(152) plt.imshow(shifted_face2, cmap=plt.cm.gray) plt.axis('off')plt.subplot(153) plt.imshow(rotated_face, cmap=plt.cm.gray) plt.axis('off')plt.subplot(154) plt.imshow(cropped_face, cmap=plt.cm.gray) plt.axis('off')plt.subplot(155) plt.imshow(zoomed_face, cmap=plt.cm.gray) plt.axis('off')plt.subplots_adjust(wspace=.05, left=.01, bottom=.01, right=.99, top=.99)

12.2 圖像濾波器

from scipy import miscface = misc.face(gray=True) face = face[:512, -512:] # crop out square on right import numpy as np noisy_face = np.copy(face).astype(float) noisy_face += face.std() * 0.5 * np.random.standard_normal(face.shape) blurred_face = ndimage.gaussian_filter(noisy_face, sigma=3) median_face = ndimage.median_filter(noisy_face, size=5) plt.figure(figsize=(9, 3.5)) plt.subplot(131) plt.imshow(noisy_face, cmap=plt.cm.gray) plt.axis('off') plt.title('noisy')plt.subplot(132) plt.imshow(blurred_face, cmap=plt.cm.gray) plt.axis('off') plt.title('Gaussian filter')plt.subplot(133) plt.imshow(median_face, cmap=plt.cm.gray) plt.axis('off') plt.title('median filter')plt.subplots_adjust(wspace=.05, left=.01, bottom=.01, right=.99, top=.99)

在scipy.ndimage.filters 和 scipy.signal 有更多應用于圖像的濾波器。

12.3 數學形態學

數學形態學是集合理論分支出來的一個數學理論。它刻畫并轉換幾何結構。特別是二元的圖像（黑白）可以用這種理論來轉換：被轉換的集合是臨近非零值像素的集合。這個理論也可以被擴展到灰度值圖像。

初級數學形態學操作使用結構化的元素，以便修改其他幾何結構。
首先讓我們生成一個結構化元素。

el = ndimage.generate_binary_structure(2, 1) el el.astype(int)

• 腐蝕(ndimage.binary_erosion)

a = np.zeros((7,7), dtype=int) a[1:6, 2:5] = 1 a ndimage.binary_erosion(a).astype(a.dtype)

#腐蝕移除了比結構小的對象 ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)

• 擴張 (ndimage.binary_dilation)

a = np.zeros((5, 5)) a[2, 2] = 1 a ndimage.binary_dilation(a).astype(a.dtype)

• 開啟(ndimage.binary_opening)

a = np.zeros((5,5), dtype=int) a[1:4, 1:4] = 1; a[4, 4] = 1 a # 開啟移除了小對象 ndimage.binary_opening(a, structure=np.ones((3,3))).astype(int) # 開啟也可以平滑拐角 ndimage.binary_opening(a).astype(int)

• 閉合: ndimage.binary_closing
  開啟操作移除小的結構，而關閉操作填滿了小洞。因此這些用來”清洗“圖像。

a = np.zeros((50, 50)) a[10:-10, 10:-10] = 1 a += 0.25*np.random.standard_normal(a.shape) mask = a>=0.5 opened_mask = ndimage.binary_opening(mask) closed_mask = ndimage.binary_closing(opened_mask) plt.figure(figsize=(12, 3.5)) plt.subplot(141) plt.imshow(a, cmap=plt.cm.gray) plt.axis('off') plt.title('a')plt.subplot(142) plt.imshow(mask, cmap=plt.cm.gray) plt.axis('off') plt.title('mask')plt.subplot(143) plt.imshow(opened_mask, cmap=plt.cm.gray) plt.axis('off') plt.title('opened_mask')plt.subplot(144) plt.imshow(closed_mask, cmap=plt.cm.gray) plt.title('closed_mask') plt.axis('off')plt.subplots_adjust(wspace=.05, left=.01, bottom=.01, right=.99, top=.99)

對于灰度值圖像，腐蝕（區別于擴張）相當于用感興趣的像素周圍的結構元素中的最小（區別于最大）值替換像素。

a = np.zeros((7,7), dtype=int) a[1:6, 1:6] = 3 a[4,4] = 2; a[2,3] = 1 a ndimage.grey_erosion(a, size=(3,3))

12.4 測量圖像

首先讓我們生成一個漂亮的人造二維圖。

x, y = np.indices((100, 100)) sig = np.sin(2*np.pi*x/50.)*np.sin(2*np.pi*y/50.)*(1+x*y/50.**2)**2 mask = sig > 1

現在讓我們看一下圖像中對象的各種信息：

labels, nb = ndimage.label(mask) nb areas = ndimage.sum(mask, labels, range(1, labels.max()+1)) areas maxima = ndimage.maximum(sig, labels, range(1, labels.max()+1)) maxima ndimage.find_objects(labels==4)

sl = ndimage.find_objects(labels==4) plt.imshow(sig[sl[0]]) plt.show()

12.5 圖像處理應用：

計數氣泡和未融化的顆粒

問題描述

打開圖像文件MV_HFV_012.jpg并且瀏覽一下。看一下imshow文檔字符串中的參數，用“右”對齊來顯示圖片（原點在左下角，而不是像標準數組在右上角）。這個掃描元素顯微圖顯示了一個帶有一些氣泡（黑色）和未溶解沙（深灰）的玻璃樣本（輕灰矩陣）。我們想要判斷樣本由三個狀態覆蓋的百分比，并且預測沙粒和氣泡的典型大小和他們的大小等。

修建圖片，刪除帶有測量信息中底部面板。

用中位數過濾稍稍過濾一下圖像以便改進它的直方圖。看一下直方圖的變化。

使用過濾后圖像的直方圖，決定允許定義沙粒像素，玻璃像素和氣泡像素掩蔽的閾限。其他的選項：寫一個函數從直方圖的最小值自動判斷閾限。

將三種不同的相用不同的顏色上色并顯示圖片。

用數學形態學清理不同的相。

為所有氣泡和沙粒做標簽，從沙粒中刪除小于10像素的掩蔽。要這樣做，用ndimage.sum或np.bincount來計算沙粒大小。

計算氣泡的平均大小。

import numpy as np from scipy import ndimage import matplotlib.pyplot as pltdate = plt.imread('Images/MV_HFV_012.jpg')

• 修建圖片，刪除帶有測量信息中底部面板。

date = date[:-60] plt.imshow(date)

• 用中位數過濾稍稍過濾一下圖像以便改進它的直方圖。看一下直方圖的變化。

filtdate = ndimage.median_filter(date, size=(7,7)) hi_date = np.histogram(date, bins=np.arange(256)) hi_filtdate = np.histogram(filtdate, bins=np.arange(256)) plt.plot(hi_date[1][:-1], hi_date[0]) plt.plot(hi_filtdate[1][:-1], hi_filtdate[0])

• 使用過濾后圖像的直方圖，決定允許定義沙粒像素，玻璃像素和氣泡像素掩蔽的閾限。其他的選項：寫一個函數從直方圖的最小值自動判斷閾限。

void = filtdate <= 50 sand = np.logical_and(filtdate > 50, filtdate <= 114) glass = filtdate > 114

• 將三種不同的相用不同的顏色上色并顯示圖片。

phases = void.astype(int) + 2*glass.astype(int) + 3*sand.astype(int) plt.imshow(phases)

單獨顯示砂礫

sand_op = ndimage.binary_opening(sand, iterations=2) plt.imshow(sand_op.astype(int))

• 為所有氣泡和沙粒做標簽，從沙粒中刪除小于10像素的掩蔽。要這樣做，用ndimage.sum或np.bincount來計算沙粒大小。

sand_labels, sand_nb = ndimage.label(sand_op) sand_areas = np.array(ndimage.sum(sand_op, sand_labels, np.arange(sand_labels.max()+1))) mask = sand_areas > 100 remove_small_sand = mask[sand_labels.ravel()].reshape(sand_labels.shape) plt.imshow(remove_small_sand.astype(int))

• 計算氣泡的平均大小。

plt.imshow(void.astype(int))

bubbles_labels, bubbles_nb = ndimage.label(void) bubbles_areas = np.bincount(bubbles_labels.ravel())[1:] mean_bubble_size = bubbles_areas.mean() median_bubble_size = np.median(bubbles_areas) mean_bubble_size, median_bubble_size

參考文獻來自桑鴻乾老師的課件：科學計算和人工智能

總結

以上是生活随笔為你收集整理的【初识SciPy库】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。