知識圖譜系列---機器學習---scikit-image圖片處理

python skimage圖像處理(一)

 

This blog is from: https://www.jianshu.com/p/f2e88197e81d 

 

基於python腳本語言開發的數字圖片處理包，比如PIL,Pillow, opencv, scikit-image等。
PIL和Pillow只提供最基礎的數字圖像處理，功能有限；opencv實際上是一個c++庫，只是提供了python接口，更新速度非常慢。scikit-image是基於scipy的一款圖像處理包，它將圖片作為numpy數組進行處理，正好與matlab一樣，因此，我們最終選擇scikit-image進行數字圖像處理。

---

Image讀出來的是PIL的類型，而skimage.io讀出來的數據是numpy格式的

import Image as img import os from matplotlib import pyplot as plot from skimage import io,transform #Image和skimage讀圖片 img_file1 = img.open('./CXR_png/MCUCXR_0042_0.png') img_file2 = io.imread('./CXR_png/MCUCXR_0042_0.png') 

輸出可以看出Img讀圖片的大小是圖片的(width, height)；而skimage的是(height,width, channel), [這也是為什么caffe在單獨測試時要要在代碼中設置：transformer.set_transpose('data',(2,0,1))，因為caffe可以處理的圖片的數據格式是(channel,height,width)，所以要轉換數據]

#讀圖片后數據的大小： print "the picture's size: ", img_file1.size print "the picture's shape: ", img_file2.shape 

the picture's size: (4892, 4020) the picture's shape: (4020, 4892) 

#得到像素： print(img_file1.getpixel((500,1000)), img_file2[500][1000]) print(img_file1.getpixel((500,1000)), img_file2[1000][500]) print(img_file1.getpixel((1000,500)), img_file2[500][1000]) 

(0, 139) (0, 0) (139, 139) 

Img讀出來的圖片獲得某點像素用getpixel((w,h))可以直接返回這個點三個通道的像素值
skimage讀出來的圖片可以直接img_file2[0][0]獲得，但是一定記住它的格式，並不是你想的(channel,height,width)

圖片信息

如果我們想知道一些skimage圖片信息

from skimage import io, data img = data.chelsea() io.imshow(img) print(type(img)) #顯示類型 print(img.shape) #顯示尺寸 print(img.shape[0]) #圖片高度 print(img.shape[1]) #圖片寬度 print(img.shape[2]) #圖片通道數 print(img.size) #顯示總像素個數 print(img.max()) #最大像素值 print(img.min()) #最小像素值 print(img.mean()) #像素平均值 print(img[0][0])#圖像的像素值 

PIL image 查看圖片信息，可用如下的方法

print type(img) print img.size #圖片的尺寸 print img.mode #圖片的模式 print img.format #圖片的格式 print(img.getpixel((0,0)))#得到像素： #img讀出來的圖片獲得某點像素用getpixel((w,h))可以直接返回這個點三個通道的像素值 

# 獲取圖像的灰度值范圍 width = img.size[0] height = img.size[1] # 輸出圖片的像素值 count = 0 for i in range(0, width): for j in range(0, height): if img.getpixel((i, j))>=0 and img.getpixel((i, j))<=255: count +=1 print count print(height*width) 

---

skimage提供了io模塊，顧名思義，這個模塊是用來圖片輸入輸出操作的。為了方便練習，也提供一個data模塊，里面嵌套了一些示例圖片，我們可以直接使用。

skimage包的子模塊

skimage包的全稱是scikit-image SciKit (toolkit for SciPy) ，它對scipy.ndimage進行了擴展，提供了更多的圖片處理功能。它是由python語言編寫的，由scipy 社區開發和維護。skimage包由許多的子模塊組成，各個子模塊提供不同的功能。主要子模塊列表如下：

子模塊名稱　 主要實現功能 io 讀取、保存和顯示圖片或視頻 data 提供一些測試圖片和樣本數據 color 顏色空間變換 filters 圖像增強、邊緣檢測、排序濾波器、自動閾值等 draw 操作於numpy數組上的基本圖形繪制，包括線條、矩形、圓和文本等 transform 幾何變換或其它變換，如旋轉、拉伸和拉東變換等 morphology 形態學操作，如開閉運算、骨架提取等 exposure 圖片強度調整，如亮度調整、直方圖均衡等 feature 特征檢測與提取等 measure 圖像屬性的測量，如相似性或等高線等 segmentation 圖像分割 restoration 圖像恢復 util 通用函數 

---

從外部讀取圖片並顯示

讀取單張彩色rgb圖片，使用skimage.io.imread（fname）函數,帶一個參數，表示需要讀取的文件路徑。顯示圖片使用skimage.io.imshow（arr）函數，帶一個參數，表示需要顯示的arr數組（讀取的圖片以numpy數組形式計算）。

from skimage import io img=io.imread('d:/dog.jpg') io.imshow(img) 

讀取單張灰度圖片，使用skimage.io.imread（fname，as_grey=True）函數，第一個參數為圖片路徑，第二個參數為as_grey, bool型值，默認為False

from skimage import io img=io.imread('d:/dog.jpg',as_grey=True) io.imshow(img) 

程序自帶圖片
skimage程序自帶了一些示例圖片，如果我們不想從外部讀取圖片，就可以直接使用這些示例圖片：

astronaut 航員圖片 coffee 一杯咖啡圖片 lena lena美女圖片 camera 拿相機的人圖片 coins 硬幣圖片 moon 月亮圖片 checkerboard 棋盤圖片 horse 馬圖片 page 書頁圖片 chelsea 小貓圖片 hubble_deep_field 星空圖片 text 文字圖片 clock 時鍾圖片 immunohistochemistry 結腸圖片 

顯示這些圖片可用如下代碼，不帶任何參數

from skimage import io, data img=data.lena() io.imshow(img) 

圖片名對應的就是函數名，如camera圖片對應的函數名為camera(). 這些示例圖片存放在skimage的安裝目錄下面，路徑名稱為data_dir,我們可以將這個路徑打印出來看看

from skimage import data_dir print(data_dir) 

保存圖片
使用io模塊的imsave（fname,arr）函數來實現。第一個參數表示保存的路徑和名稱，第二個參數表示需要保存的數組變量。

from skimage import io,data img=data.chelsea() io.imshow(img) io.imsave('d:/cat.jpg',img) 

保存圖片的同時也起到了轉換格式的作用。如果讀取時圖片格式為jpg圖片，保存為png格式，則將圖片從jpg圖片轉換為png圖片並保存。

圖片信息

如果我們想知道一些圖片信息

from skimage import io, data img = data.chelsea() io.imshow(img) print(type(img)) #顯示類型 print(img.shape) #顯示尺寸 print(img.shape[0]) #圖片高度 print(img.shape[1]) #圖片寬度 print(img.shape[2]) #圖片通道數 print(img.size) #顯示總像素個數 print(img.max()) #最大像素值 print(img.min()) #最小像素值 print(img.mean()) #像素平均值 print(img[0][0])#圖像的像素值 

---

圖像像素的訪問與裁剪

圖片讀入程序中后，是以numpy數組存在的。因此對numpy數組的一切功能，對圖片也適用。對數組元素的訪問，實際上就是對圖片像素點的訪問。

彩色圖片訪問方式為：img[i,j,c]
i表示圖片的行數，j表示圖片的列數，c表示圖片的通道數（RGB三通道分別對應0，1，2）。坐標是從左上角開始。

灰度圖片訪問方式為：gray[i,j]

例1：輸出小貓圖片的G通道中的第20行30列的像素值

from skimage import io,data img=data.chelsea() pixel=img[20,30,1] print(pixel) 

例2：顯示紅色單通道圖片

from skimage import io,data img=data.chelsea() R=img[:,:,0] io.imshow(R) 

除了對像素進行讀取，也可以修改像素值。

例3：對小貓圖片隨機添加椒鹽噪聲

from skimage import io,data import numpy as np img=data.chelsea() #隨機生成5000個椒鹽 rows,cols,dims=img.shape for i in range(5000): x=np.random.randint(0,rows) y=np.random.randint(0,cols) img[x,y,:]=255 io.imshow(img) 

這里用到了numpy包里的random來生成隨機數，randint(0,cols)表示隨機生成一個整數，范圍在0到cols之間。
用img[x,y,:]=255這句來對像素值進行修改，將原來的三通道像素值，變為255

通過對數組的裁剪，就可以實現對圖片的裁剪。
例4：對小貓圖片進行裁剪

from skimage import io,data img=data.chelsea() roi=img[80:180,100:200,:] io.imshow(roi) 

對多個像素點進行操作，使用數組切片方式訪問。切片方式返回的是以指定間隔下標訪問 該數組的像素值。下面是有關灰度圖像的一些例子：

img[i,:] = im[j,:] # 將第 j 行的數值賦值給第 i 行 img[:,i] = 100 # 將第 i 列的所有數值設為 100 img[:100,:50].sum() # 計算前 100 行、前 50 列所有數值的和 img[50:100,50:100] # 50~100 行，50~100 列（不包括第 100 行和第 100 列） img[i].mean() # 第 i 行所有數值的平均值 img[:,-1] # 最后一列 img[-2,:] (or im[-2]) # 倒數第二行 

最后我們再看兩個對像素值進行訪問和改變的例子：

例5：將lena圖片進行二值化，像素值大於128的變為1，否則變為0

from skimage import io,data,color img=data.lena() img_gray=color.rgb2gray(img) rows,cols=img_gray.shape for i in range(rows): for j in range(cols): if (img_gray[i,j]<=0.5): img_gray[i,j]=0 else: img_gray[i,j]=1 io.imshow(img_gray) 

使用了color模塊的rgb2gray（）函數，將彩色三通道圖片轉換成灰度圖。轉換結果為float64類型的數組，范圍為[0,1]之間。

將彩色三通道圖片轉換成灰度圖,最后變成unit8, float轉換為unit8是有信息損失的。

img_path = 'data/dpclassifier/newtrain/test/1_0.png' import Image as img import os from matplotlib import pyplot as plot from skimage import io,transform, img_as_ubyte img_file1 = img.open(img_path) img_file1.show() img_file2 = io.imread(img_path) io.imshow(img_file2) print(type(img_file1),img_file1.mode, type(img_file2),img_file2.shape, img_file2.dtype,img_file2.max(),img_file2.min(),img_file2.mean()) img_file22=skimage.color.rgb2gray(img_file2) print(type(img_file22),img_file22.shape,img_file22.dtype,img_file22.max(),img_file22.min(),img_file22.mean() ) dst=img_as_ubyte(img_file22) print(type(dst),dst.shape,dst.dtype, dst.max(), dst.min(), dst.mean()) 

結果

(<class 'PIL.PngImagePlugin.PngImageFile'>, 'RGB', <type 'numpy.ndarray'>, (420, 512, 3), dtype('uint8'), 255, 0, 130.9983863467262) (<type 'numpy.ndarray'>, (420, 512), dtype('float64'), 1.0, 0.0, 0.5137191621440242) (<type 'numpy.ndarray'>, (420, 512), dtype('uint8'), 255, 0, 130.9983863467262) 

例6：

from skimage import io,data img=data.chelsea() reddish = img[:, :, 0] >170 img[reddish] = [0, 255, 0] io.imshow(img) 

這個例子先對R通道的所有像素值進行判斷，如果大於170，則將這個地方的像素值變為[0,255,0], 即G通道值為255，R和B通道值為0。

---

圖像數據類型及顏色空間轉換

在skimage中，一張圖片就是一個簡單的numpy數組，數組的數據類型有很多種，相互之間也可以轉換。這些數據類型及取值范圍如下表所示：

Data type Range uint8 0 to 255 uint16 0 to 65535 uint32 0 to 232 float -1 to 1 or 0 to 1 int8 -128 to 127 int16 -32768 to 32767 int32 -231 to 231 - 1 

一張圖片的像素值范圍是[0,255], 因此默認類型是unit8, 可用如下代碼查看數據類型

from skimage import io,data img=data.chelsea() print(img.dtype.name) 

在上面的表中，特別注意的是float類型，它的范圍是[-1,1]或[0,1]之間。一張彩色圖片轉換為灰度圖后，它的類型就由unit8變成了float
1、unit8轉float

from skimage import data,img_as_float img=data.chelsea() print(img.dtype.name) dst=img_as_float(img) print(dst.dtype.name) 

2、float轉uint8

from skimage import img_as_ubyte import numpy as np img = np.array([0, 0.5, 1], dtype=float) print(img.dtype.name) dst=img_as_ubyte(img) print(dst.dtype.name) 

float轉為unit8，有可能會造成數據的損失，因此會有警告提醒。*

除了這兩種最常用的轉換以外，其實有一些其它的類型轉換，如下表：

Function name Description img_as_float Convert to 64-bit floating point. img_as_ubyte Convert to 8-bit uint. img_as_uint Convert to 16-bit uint. img_as_int Convert to 16-bit int. 

如前所述，除了直接轉換可以改變數據類型外，還可以通過圖像的顏色空間轉換來改變數據類型。

常用的顏色空間有灰度空間、rgb空間、hsv空間和cmyk空間。顏色空間轉換以后，圖片類型都變成了float型。

所有的顏色空間轉換函數，都放在skimage的color模塊內。
例：rgb轉灰度圖

from skimage import io,data,color img=data.lena() gray=color.rgb2gray(img) io.imshow(gray) 

其它的轉換，用法都是一樣的，列舉常用的如下：

skimage.color.rgb2grey(rgb) skimage.color.rgb2hsv(rgb) skimage.color.rgb2lab(rgb) skimage.color.gray2rgb(image) skimage.color.hsv2rgb(hsv) skimage.color.lab2rgb(lab) 

實際上，上面的所有轉換函數，都可以用一個函數來代替

skimage.color.convert_colorspace(arr, fromspace, tospace) 

表示將arr從fromspace顏色空間轉換到tospace顏色空間。

例：rgb轉hsv

from skimage import io,data,color img=data.lena() hsv=color.convert_colorspace(img,'RGB','HSV') io.imshow(hsv) 

在color模塊的顏色空間轉換函數中，還有一個比較有用的函數是
skimage.color.label2rgb(arr), 可以根據標簽值對圖片進行着色。以后的圖片分類后着色就可以用這個函數。

例：將lena圖片分成三類，然后用默認顏色對三類進行着色

from skimage import io,data,color import numpy as np img=data.lena() gray=color.rgb2gray(img) rows,cols=gray.shape labels=np.zeros([rows,cols]) for i in range(rows): for j in range(cols): if(gray[i,j]<0.4): labels[i,j]=0 elif(gray[i,j]<0.75): labels[i,j]=1 else: labels[i,j]=2 dst=color.label2rgb(labels) io.imshow(dst) 

---

圖像的繪制

實際上前面我們就已經用到了圖像的繪制，如：

io.imshow(img) 

這一行代碼的實質是利用matplotlib包對圖片進行繪制，繪制成功后，返回一個matplotlib類型的數據。因此，我們也可以這樣寫：

import matplotlib.pyplot as plt plt.imshow(img) 

imshow()函數格式為：

matplotlib.pyplot.imshow(X, cmap=None) 

X: 要繪制的圖像或數組。
cmap: 顏色圖譜（colormap), 默認繪制為RGB(A)顏色空間。
其它可選的顏色圖譜如下列表：

顏色圖譜 描述 autumn 紅-橙-黃 bone 黑-白，x線 cool 青-洋紅 copper 黑-銅 flag 紅-白-藍-黑 gray 黑-白 hot 黑-紅-黃-白 hsv hsv顏色空間， 紅-黃-綠-青-藍-洋紅-紅 inferno 黑-紅-黃 jet 藍-青-黃-紅 magma 黑-紅-白 pink 黑-粉-白 plasma 綠-紅-黃 prism 紅-黃-綠-藍-紫-...-綠模式 spring 洋紅-黃 summer 綠-黃 viridis 藍-綠-黃 winter 藍-綠 

用的比較多的有gray,jet等，如：

plt.imshow(image,plt.cm.gray) plt.imshow(img,cmap=plt.cm.jet) 

在窗口上繪制完圖片后，返回一個AxesImage對象。要在窗口上顯示這個對象，我們可以調用show()函數來進行顯示，但進行練習的時候（ipython環境中），一般我們可以省略show（）函數，也能自動顯示出來。

from skimage import io,data img=data.astronaut() dst=io.imshow(img) print(type(dst)) io.show() 

可以看到，類型是'matplotlib.image.AxesImage'。顯示一張圖片，我們通常更願意這樣寫：

import matplotlib.pyplot as plt from skimage import io,data img=data.astronaut() plt.imshow(img) plt.show() 

matplotlib是一個專業繪圖的庫，相當於matlab中的plot,可以設置多個figure窗口,設置figure的標題，隱藏坐標尺，甚至可以使用subplot在一個figure中顯示多張圖片。一般我們可以這樣導入matplotlib庫：

import matplotlib.pyplot as plt 

也就是說，我們繪圖實際上用的是matplotlib包的pyplot模塊。

用figure函數和subplot函數分別創建主窗口與子圖
分開並同時顯示宇航員圖片的三個通道

from skimage import data import matplotlib.pyplot as plt img=data.astronaut() plt.figure(num='astronaut',figsize=(8,8)) #創建一個名為astronaut的窗口,並設置大小  plt.subplot(2,2,1) #將窗口分為兩行兩列四個子圖，則可顯示四幅圖片 plt.title('origin image') #第一幅圖片標題 plt.imshow(img) #繪制第一幅圖片 plt.subplot(2,2,2) #第二個子圖 plt.title('R channel') #第二幅圖片標題 plt.imshow(img[:,:,0],plt.cm.gray) #繪制第二幅圖片,且為灰度圖 plt.axis('off') #不顯示坐標尺寸 plt.subplot(2,2,3) #第三個子圖 plt.title('G channel') #第三幅圖片標題 plt.imshow(img[:,:,1],plt.cm.gray) #繪制第三幅圖片,且為灰度圖 plt.axis('off') #不顯示坐標尺寸 plt.subplot(2,2,4) #第四個子圖 plt.title('B channel') #第四幅圖片標題 plt.imshow(img[:,:,2],plt.cm.gray) #繪制第四幅圖片,且為灰度圖 plt.axis('off') #不顯示坐標尺寸 plt.show() #顯示窗口 

在圖片繪制過程中，我們用matplotlib.pyplot模塊下的figure（）函數來創建顯示窗口，該函數的格式為：

matplotlib.pyplot.figure(num=None, figsize=None, dpi=None, facecolor=None, edgecolor=None) 

所有參數都是可選的，都有默認值，因此調用該函數時可以不帶任何參數，其中：

num: 整型或字符型都可以。如果設置為整型，則該整型數字表示窗口的序號。如果設置為字符型，則該字符串表示窗口的名稱。用該參數來命名窗口，如果兩個窗口序號或名相同，則后一個窗口會覆蓋前一個窗口。

figsize: 設置窗口大小。是一個tuple型的整數，如figsize=（8，8）
dpi: 整形數字，表示窗口的分辨率。
facecolor: 窗口的背景顏色。
edgecolor: 窗口的邊框顏色。
用figure()函數創建的窗口，只能顯示一幅圖片，如果想要顯示多幅圖片，則需要將這個窗口再划分為幾個子圖，在每個子圖中顯示不同的圖片。我們可以使用subplot（）函數來划分子圖，函數格式為：

matplotlib.pyplot.subplot(nrows, ncols, plot_number) 

nrows: 子圖的行數。
ncols: 子圖的列數。
plot_number: 當前子圖的編號。

如：

plt.subplot(2,2,1) 

則表示將figure窗口划分成了2行2列共4個子圖，當前為第1個子圖。我們有時也可以用這種寫法：

plt.subplot(221) 

兩種寫法效果是一樣的。每個子圖的標題可用title()函數來設置，是否使用坐標尺可用axis()函數來設置，如：

plt.subplot(221) plt.title("first subwindow") plt.axis('off') 

用subplots來創建顯示窗口與划分子圖

除了上面那種方法創建顯示窗口和划分子圖，還有另外一種編寫方法也可以，如下例:

import matplotlib.pyplot as plt from skimage import data,color img = data.immunohistochemistry() hsv = color.rgb2hsv(img) fig, axes = plt.subplots(2, 2, figsize=(7, 6)) ax0, ax1, ax2, ax3 = axes.ravel() ax0.imshow(img) ax0.set_title("Original image") ax1.imshow(hsv[:, :, 0], cmap=plt.cm.gray) ax1.set_title("H") ax2.imshow(hsv[:, :, 1], cmap=plt.cm.gray) ax2.set_title("S") ax3.imshow(hsv[:, :, 2], cmap=plt.cm.gray) ax3.set_title("V") for ax in axes.ravel(): ax.axis('off') fig.tight_layout() #自動調整subplot間的參數 

直接用subplots()函數來創建並划分窗口。注意，比前面的subplot()函數多了一個s，該函數格式為：

matplotlib.pyplot.subplots(nrows=1, ncols=1)

nrows: 所有子圖行數，默認為1。

ncols: 所有子圖列數，默認為1。

返回一個窗口figure, 和一個tuple型的ax對象，該對象包含所有的子圖,可結合ravel()函數列出所有子圖，如：

fig, axes = plt.subplots(2, 2, figsize=(7, 6))
ax0, ax1, ax2, ax3 = axes.ravel()

創建了2行2列4個子圖，分別取名為ax0,ax1,ax2和ax3, 每個子圖的標題用set_title()函數來設置，如：

ax0.imshow(img)
ax0.set_title("Original image")

如果有多個子圖，我們還可以使用tight_layout()函數來調整顯示的布局，該函數格式為：

matplotlib.pyplot.tight_layout(pad=1.08, h_pad=None, w_pad=None, rect=None) 

所有的參數都是可選的，調用該函數時可省略所有的參數。
pad: 主窗口邊緣和子圖邊緣間的間距，默認為1.08
h_pad, w_pad: 子圖邊緣之間的間距，默認為 pad_inches
rect: 一個矩形區域，如果設置這個值，則將所有的子圖調整到這個矩形區域內。
一般調用為：

plt.tight_layout() #自動調整subplot間的參數 

其它方法繪圖並顯示
除了使用matplotlib庫來繪制圖片，skimage還有另一個子模塊viewer，也提供一個函數來顯示圖片。不同的是，它利用Qt工具來創建一塊畫布，從而在畫布上繪制圖像。

例：

from skimage import data from skimage.viewer import ImageViewer img = data.coins() viewer = ImageViewer(img) viewer.show() 

最后總結一下，繪制和顯示圖片常用到的函數有：

函數名 功能 調用格式 figure 創建一個顯示窗口 plt.figure(num=1,figsize=(8,8) imshow 繪制圖片 plt.imshow(image) show 顯示窗口 plt.show() subplot 划分子圖 plt.subplot(2,2,1) title 設置子圖標題(與subplot結合使用） plt.title('origin image') axis 是否顯示坐標尺 plt.axis('off') subplots 創建帶有多個子圖的窗口 fig,axes=plt.subplots(2,2,figsize=(8,8)) ravel 為每個子圖設置變量 ax0,ax1,ax2,ax3=axes.ravel() set_title 設置子圖標題（與axes結合使用） ax0.set_title('first window') tight_layout 自動調整子圖顯示布局 plt.tight_layout() 

---

圖像的批量處理

有些時候，我們不僅要對一張圖片進行處理，可能還會對一批圖片處理。這時候，我們可以通過循環來執行處理，也可以調用程序自帶的圖片集合來處理。
圖片集合函數為：

skimage.io.ImageCollection(load_pattern,load_func=None) 

這個函數是放在io模塊內的，帶兩個參數，第一個參數load_pattern, 表示圖片組的路徑，可以是一個str字符串。第二個參數load_func是一個回調函數，我們對圖片進行批量處理就可以通過這個回調函數實現。回調函數默認為imread(),即默認這個函數是批量讀取圖片。
先看一個例子：

import skimage.io as io from skimage import data_dir str=data_dir + '/*.png' coll = io.ImageCollection(str) print(len(coll)) 

顯示結果為25, 說明系統自帶了25張png的示例圖片，這些圖片都讀取了出來，放在圖片集合coll里。如果我們想顯示其中一張圖片，則可以在后加上一行代碼：

io.imshow(coll[10]) 

顯示為：

 

 

如果一個文件夾里，我們既存放了一些jpg格式的圖片，又存放了一些png格式的圖片，現在想把它們全部讀取出來，該怎么做呢?

import skimage.io as io from skimage import data_dir str='d:/pic/*.jpg:d:/pic/*.png' coll = io.ImageCollection(str) print(len(coll)) 

注意這個地方'd:/pic/.jpg:d:/pic/.png' ，是兩個字符串合在一起的，第一個是'd:/pic/.jpg', 第二個是'd:/pic/.png' ，合在一起后，中間用冒號來隔開，這樣就可以把d:/pic/文件夾下的jpg和png格式的圖片都讀取出來。如果還想讀取存放在其它地方的圖片，也可以一並加進去，只是中間同樣用冒號來隔開。
io.ImageCollection()這個函數省略第二個參數，就是批量讀取。如果我們不是想批量讀取，而是其它批量操作，如批量轉換為灰度圖，那又該怎么做呢？
那就需要先定義一個函數，然后將這個函數作為第二個參數，如：

from skimage import data_dir,io,color def convert_gray(f): rgb=io.imread(f) return color.rgb2gray(rgb) str=data_dir+'/*.png' coll = io.ImageCollection(str,load_func=convert_gray) io.imshow(coll[10]) 

 

 

這種批量操作對視頻處理是極其有用的，因為視頻就是一系列的圖片組合

from skimage import data_dir,io,color class AVILoader: video_file = 'myvideo.avi' def __call__(self, frame): return video_read(self.video_file, frame) avi_load = AVILoader() frames = range(0, 1000, 10) # 0, 10, 20, ...ic =io.ImageCollection(frames, load_func=avi_load) 

這段代碼的意思，就是將myvideo.avi這個視頻中每隔10幀的圖片讀取出來，放在圖片集合中。
得到圖片集合以后，我們還可以將這些圖片連接起來，構成一個維度更高的數組，連接圖片的函數為：

skimage.io.concatenate_images(ic) 

帶一個參數，就是以上的圖片集合，如：

from skimage import data_dir,io,color coll = io.ImageCollection('d:/pic/*.jpg') mat=io.concatenate_images(coll) 

使用concatenate_images(ic)函數的前提是讀取的這些圖片尺寸必須一致，否則會出錯。我們看看圖片連接前后的維度變化：

from skimage import data_dir,io,color coll = io.ImageCollection('d:/pic/*.jpg') print(len(coll)) #連接的圖片數量 print(coll[0].shape) #連接前的圖片尺寸，所有的都一樣 mat=io.concatenate_images(coll) print(mat.shape) #連接后的數組尺寸 

顯示結果：

2 (870, 580, 3) (2, 870, 580, 3) 

可以看到，將2個3維數組，連接成了一個4維數組
如果我們對圖片進行批量操作后，想把操作后的結果保存起來，也是可以辦到的。
例：把系統自帶的所有png示例圖片，全部轉換成256256的jpg格式灰度圖，保存在d:/data/文件夾下*
改變圖片的大小，我們可以使用tranform模塊的resize()函數，后續會講到這個模塊。

from skimage import data_dir,io,transform,color import numpy as np def convert_gray(f): rgb=io.imread(f) #依次讀取rgb圖片  gray=color.rgb2gray(rgb) #將rgb圖片轉換成灰度圖  dst=transform.resize(gray,(256,256)) #將灰度圖片大小轉換為256*256  return dst str=data_dir+'/*.png' coll = io.ImageCollection(str,load_func=convert_gray) for i in range(len(coll)): io.imsave('d:/data/'+np.str(i)+'.jpg',coll[i]) #循環保存圖片 

結果：

 

 

---

圖像的形變與縮放

圖像的形變與縮放，使用的是skimage的transform模塊，函數比較多，功能齊全。
1、改變圖片尺寸resize
函數格式為：

skimage.transform.resize(image, output_shape) 

image: 需要改變尺寸的圖片
output_shape: 新的圖片尺寸

from skimage import transform,data import matplotlib.pyplot as plt img = data.camera() dst=transform.resize(img, (80, 60)) plt.figure('resize') plt.subplot(121) plt.title('before resize') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('before resize') plt.imshow(dst,plt.cm.gray) plt.show() 

將camera圖片由原來的512x512大小，變成了80x60大小。從下圖中的坐標尺，我們能夠看出來：

 

 

2、按比例縮放rescale
函數格式為：

skimage.transform.rescale(image, scale[, ...]) 

scale參數可以是單個float數，表示縮放的倍數，也可以是一個float型的tuple，如[0.2,0.5],表示將行列數分開進行縮放

from skimage import transform,data img = data.camera() print(img.shape) #圖片原始大小  print(transform.rescale(img, 0.1).shape) #縮小為原來圖片大小的0.1 print(transform.rescale(img, [0.5,0.25]).shape) #縮小為原來圖片行數一半，列數四分之一 print(transform.rescale(img, 2).shape) #放大為原來圖片大小的2倍 

結果為：

(512, 512) (51, 51) (256, 128) (1024, 1024) 

3、旋轉 rotate

skimage.transform.rotate(image, angle[, ...],resize=False) 

angle參數是個float類型數，表示旋轉的度數
resize用於控制在旋轉時，是否改變大小 ，默認為False

from skimage import transform,data import matplotlib.pyplot as plt img = data.camera() print(img.shape) #圖片原始大小 img1=transform.rotate(img, 60) #旋轉90度，不改變大小  print(img1.shape) img2=transform.rotate(img, 30,resize=True) #旋轉30度，同時改變大小 print(img2.shape) plt.figure('resize') plt.subplot(121)plt.title('rotate 60') plt.imshow(img1,plt.cm.gray) plt.subplot(122) plt.title('rotate 30') plt.imshow(img2,plt.cm.gray) plt.show() 

顯示結果:

 

 

 

4、圖像金字塔
以多分辨率來解釋圖像的一種有效但概念簡單的結構就是圖像金字塔。圖像金字塔最初用於機器視覺和圖像壓縮，一幅圖像的金字塔是一系列以金字塔形狀排列的分辨率逐步降低的圖像集合。金字塔的底部是待處理圖像的高分辨率表示，而頂部是低分辨率的近似。當向金字塔的上層移動時，尺寸和分辨率就降低。
在此，我們舉一個高斯金字塔的應用實例，函數原型為：

skimage.transform.pyramid_gaussian(image, downscale=2) 

downscale控制着金字塔的縮放比例

import numpy as np import matplotlib.pyplot as plt from skimage import data,transform image = data.astronaut() #載入宇航員圖片 rows, cols, dim = image.shape #獲取圖片的行數，列數和通道數 pyramid = tuple(transform.pyramid_gaussian(image, downscale=2)) #產生高斯金字塔圖像#共生成了log(512)=9幅金字塔圖像，加上原始圖像共10幅，pyramid[0]-pyramid[1] composite_image = np.ones((rows, cols + cols / 2, 3), dtype=np.double) #生成背景composite_image[:rows, :cols, :] = pyramid[0] #融合原始圖像 i_row = 0 for p in pyramid[1:]: n_rows, n_cols = p.shape[:2] composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = p #循環融合9幅金字塔圖像 i_row += n_rows plt.imshow(composite_image) plt.show() 

 

 

 

上右圖，就是10張金字塔圖像，下標為0的表示原始圖像，后面每層的圖像行和列變為上一層的一半，直至變為1
除了高斯金字塔外，還有其它的金字塔，如：

skimage.transform.pyramid_laplacian(image, downscale=2) 

 

---

對比度與亮度調整

圖像亮度與對比度的調整，是放在skimage包的exposure模塊里面
1、gamma調整
原理：I=Ig

對原圖像的像素，進行冪運算，得到新的像素值。公式中的g就是gamma值。
如果gamma>1, 新圖像比原圖像暗
如果gamma<1,新圖像比原圖像亮
函數格式為：

skimage.exposure.adjust_gamma(image, gamma=1) 

gamma參數默認為1，原像不發生變化 。

from skimage import data, exposure, img_as_float import matplotlib.pyplot as plt image = img_as_float(data.moon()) gam1= exposure.adjust_gamma(image, 2) #調暗 gam2= exposure.adjust_gamma(image, 0.5) #調亮plt.figure('adjust_gamma',figsize=(8,8)) plt.subplot(131)plt.title('origin image') plt.imshow(image,plt.cm.gray)plt.axis('off') plt.subplot(132) plt.title('gamma=2') plt.imshow(gam1,plt.cm.gray) plt.axis('off') plt.subplot(133) plt.title('gamma=0.5') plt.imshow(gam2,plt.cm.gray) plt.axis('off') plt.show() 

 

 

2、log對數調整
這個剛好和gamma相反
原理：I=log(I)

from skimage import data, exposure, img_as_float import matplotlib.pyplot as plt image = img_as_float(data.moon()) gam1= exposure.adjust_log(image) #對數調整 plt.figure('adjust_gamma',figsize=(8,8)) plt.subplot(121)plt.title('origin image') plt.imshow(image,plt.cm.gray) plt.axis('off') plt.subplot(122) plt.title('log') plt.imshow(gam1,plt.cm.gray) plt.axis('off') plt.show() 

 

3、判斷圖像對比度是否偏低
函數：is_low_contrast(img)
返回一個bool型值

from skimage import data, exposure image =data.moon() result=exposure.is_low_contrast(image) print(result) 

輸出為False

4、調整強度
函數：

skimage.exposure.rescale_intensity(image, in_range='image', out_range='dtype') 

in_range 表示輸入圖片的強度范圍，默認為'image', 表示用圖像的最大/最小像素值作為范圍
out_range 表示輸出圖片的強度范圍，默認為'dype', 表示用圖像的類型的最大/最小值作為范圍
默認情況下，輸入圖片的[min,max]范圍被拉伸到[dtype.min, dtype.max]，如果dtype=uint8, 那么dtype.min=0, dtype.max=255

import numpy as np from skimage import exposure image = np.array([51, 102, 153], dtype=np.uint8) mat=exposure.rescale_intensity(image) print(mat) 

輸出為[ 0 127 255]
即像素最小值由51變為0，最大值由153變為255，整體進行了拉伸，但是數據類型沒有變，還是uint8
前面我們講過，可以通過img_as_float()函數將unit8類型轉換為float型，實際上還有更簡單的方法，就是乘以1.0

import numpy as np image = np.array([51, 102, 153], dtype=np.uint8) print(image*1.0) 

即由[51,102,153]變成了[ 51. 102. 153.]
而float類型的范圍是[0,1]，因此對float進行rescale_intensity 調整后，范圍變為[0,1],而不是[0,255]

import numpy as np from skimage import exposure image = np.array([51, 102, 153], dtype=np.uint8) tmp=image*1.0 mat=exposure.rescale_intensity(tmp) print(mat) 

結果為[ 0. 0.5 1. ]
如果原始像素值不想被拉伸，只是等比例縮小，就使用in_range參數，如：

import numpy as np from skimage import exposure image = np.array([51, 102, 153], dtype=np.uint8) tmp=image*1.0 mat=exposure.rescale_intensity(tmp,in_range=(0,255)) print(mat) 

輸出為：[ 0.2 0.4 0.6]，即原像素值除以255
如果參數in_range的[main,max]范圍要比原始像素值的范圍[min,max] 大或者小，那就進行裁剪，如：

mat=exposure.rescale_intensity(tmp,in_range=(0,102)) print(mat) 

輸出[ 0.5 1. 1. ]，即原像素值除以102，超出1的變為1
如果一個數組里面有負數，現在想調整到正數，就使用out_range參數。如：

import numpy as np from skimage import exposure image = np.array([-10, 0, 10], dtype=np.int8) mat=exposure.rescale_intensity(image, out_range=(0, 127)) print(mat) 

輸出[ 0 63 127]

---

直方圖與均衡化

在圖像處理中，直方圖是非常重要，也是非常有用的一個處理要素。
在skimage庫中對直方圖的處理，是放在exposure這個模塊中。
1、計算直方圖
函數：

skimage.exposure.histogram(image, nbins=256) 

在numpy包中，也提供了一個計算直方圖的函數histogram(),兩者大同小義。
返回一個tuple（hist, bins_center), 前一個數組是直方圖的統計量，后一個數組是每個bin的中間值

import numpy as np from skimage import exposure,data image =data.camera()*1.0 hist1=np.histogram(image, bins=2) #用numpy包計算直方圖hist2=exposure.histogram(image, nbins=2) #用skimage計算直方圖 print(hist1) print(hist2) 

輸出：

(array([107432, 154712], dtype=int64), array([ 0. , 127.5, 255. ])) (array([107432, 154712], dtype=int64), array([ 63.75, 191.25])) 

分成兩個bin，每個bin的統計量是一樣的，但numpy返回的是每個bin的兩端的范圍值，而skimage返回的是每個bin的中間值

2、繪制直方圖
繪圖都可以調用matplotlib.pyplot庫來進行，其中的hist函數可以直接繪制直方圖。
調用方式：

n, bins, patches = plt.hist(arr, bins=10, normed=0, facecolor='black', edgecolor='black',alpha=1，histtype='bar') 

hist的參數非常多，但常用的就這六個，只有第一個是必須的，后面四個可選

arr: 需要計算直方圖的一維數組 bins: 直方圖的柱數，可選項，默認為10 normed: 是否將得到的直方圖向量歸一化。默認為0 facecolor: 直方圖顏色 edgecolor: 直方圖邊框顏色 alpha: 透明度 histtype: 直方圖類型，‘bar’, ‘barstacked’, ‘step’, ‘stepfilled’ 

返回值 ：

n: 直方圖向量，是否歸一化由參數normed設定 bins: 返回各個bin的區間范圍 patches: 返回每個bin里面包含的數據，是一個list 

from skimage import data import matplotlib.pyplot as plt img=data.camera() plt.figure("hist") arr=img.flatten() n, bins, patches = plt.hist(arr, bins=256, normed=1,edgecolor='None',facecolor='red') plt.show() 

 

 

其中的flatten()函數是numpy包里面的，用於將二維數組序列化成一維數組。
是按行序列，如

mat=[[1 2 3 　　　 4 5 6]] 

經過 mat.flatten()后，就變成了

mat=[1 2 3 4 5 6] 

3、彩色圖片三通道直方圖

一般來說直方圖都是征對灰度圖的，如果要畫rgb圖像的三通道直方圖，實際上就是三個直方圖的疊加。

from skimage import data import matplotlib.pyplot as plt img=data.lena() ar=img[:,:,0].flatten() plt.hist(ar, bins=256, normed=1,facecolor='r',edgecolor='r',hold=1) ag=img[:,:,1].flatten() plt.hist(ag, bins=256, normed=1, facecolor='g',edgecolor='g',hold=1) ab=img[:,:,2].flatten() plt.hist(ab, bins=256, normed=1, facecolor='b',edgecolor='b') plt.show() 

其中，加一個參數hold=1,表示可以疊加

 

 

4、直方圖均衡化

如果一副圖像的像素占有很多的灰度級而且分布均勻，那么這樣的圖像往往有高對比度和多變的灰度色調。直方圖均衡化就是一種能僅靠輸入圖像直方圖信息自動達到這種效果的變換函數。它的基本思想是對圖像中像素個數多的灰度級進行展寬，而對圖像中像素個數少的灰度進行壓縮，從而擴展取值的動態范圍，提高了對比度和灰度色調的變化，使圖像更加清晰。

from skimage import data,exposure import matplotlib.pyplot as plt img=data.moon() plt.figure("hist",figsize=(8,8)) arr=img.flatten() plt.subplot(221) plt.imshow(img,plt.cm.gray) #原始圖像 plt.subplot(222) plt.hist(arr, bins=256, normed=1,edgecolor='None',facecolor='red') #原始圖像直方圖 img1=exposure.equalize_hist(img) arr1=img1.flatten() plt.subplot(223) plt.imshow(img1,plt.cm.gray) #均衡化圖像 plt.subplot(224) plt.hist(arr1, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() 

 

CLAHE

skimage.exposure.``equalize_adapthist(image, kernel_size=None, clip_limit=0.01, nbins=256)

Contrast Limited Adaptive Histogram Equalization (CLAHE).

An algorithm for local contrast enhancement, that uses histograms computed over different tile regions of the image. Local details can therefore be enhanced even in regions that are darker or lighter than most of the image.

image : (M, N[, C]) ndarray

Input image.
kernel_size: integer or list-like, optional
Defines the shape of contextual regions used in the algorithm. If iterable is passed, it must have the same number of elements as image.ndim (without color channel). If integer, it is broadcasted to each image dimension. By default, kernel_size is 1/8 of image height by 1/8 of its width.

clip_limit : float, optional

Clipping limit, normalized between 0 and 1 (higher values give more contrast).
nbins : int, optional
Number of gray bins for histogram (“data range”).

| Returns: |

out : (M, N[, C]) ndarray

Equalized image.

http://scikit-image.org/docs/dev/api/skimage.exposure.html#equalize-adapthist

from skimage import data,exposure import matplotlib.pyplot as plt #%matplotlib notebook clip_limitnumber=0.01 img=data.moon() print(img.shape) plt.figure("hist",figsize=(8,8)) arr=img.flatten() plt.subplot(5,2,1) plt.title('original') plt.imshow(img,plt.cm.gray) #原始圖像 plt.subplot(5,2,2) plt.hist(arr, bins=256, normed=1,edgecolor='None',facecolor='red') #原始圖像直方圖 # #img1=exposure.equalize_hist(img) # img1=exposure.equalize_hist(img) # arr1=img1.flatten() # plt.subplot(6,2,3) # plt.title('equalize_hist') # plt.imshow(img1,plt.cm.gray) #均衡化圖像 # plt.subplot(6,2,4) # plt.hist(arr1, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 # plt.show() img2=exposure.equalize_adapthist(img, kernel_size=256, clip_limit=clip_limitnumber, nbins=256) arr2=img2.flatten() plt.subplot(5,2,3) plt.title('equalize_adapthist-256-'+ str(clip_limitnumber)) plt.imshow(img2,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,4) plt.hist(arr2, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img3=exposure.equalize_adapthist(img, kernel_size=128, clip_limit=clip_limitnumber, nbins=256) arr3=img3.flatten() plt.subplot(5,2,5) plt.title('equalize_adapthist-128-'+ str(clip_limitnumber)) plt.imshow(img3,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,6) plt.hist(arr3, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img4=exposure.equalize_adapthist(img, kernel_size=64, clip_limit=clip_limitnumber, nbins=256) arr4=img4.flatten() plt.subplot(5,2,7) plt.title('equalize_adapthist-64-'+ str(clip_limitnumber)) plt.imshow(img4,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,8) plt.hist(arr4, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img5=exposure.equalize_adapthist(img, kernel_size=32, clip_limit=clip_limitnumber, nbins=256) arr5=img5.flatten() plt.subplot(5,2,9) plt.title('equalize_adapthist-32-'+ str(clip_limitnumber)) plt.imshow(img5,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,10) plt.hist(arr5, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() 

 

1.png

from skimage import data,exposure import matplotlib.pyplot as plt #%matplotlib notebook kernel_sizenuber=64 img=data.moon() print(img.shape) plt.figure("hist",figsize=(8,8)) arr=img.flatten() plt.subplot(5,2,1) plt.title('original') plt.imshow(img,plt.cm.gray) #原始圖像 plt.subplot(5,2,2) plt.hist(arr, bins=256, normed=1,edgecolor='None',facecolor='red') #原始圖像直方圖 # #img1=exposure.equalize_hist(img) # img1=exposure.equalize_hist(img) # arr1=img1.flatten() # plt.subplot(6,2,3) # plt.title('equalize_hist') # plt.imshow(img1,plt.cm.gray) #均衡化圖像 # plt.subplot(6,2,4) # plt.hist(arr1, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 # plt.show() img2=exposure.equalize_adapthist(img, kernel_size=kernel_sizenuber, clip_limit=0.001, nbins=256) arr2=img2.flatten() plt.subplot(5,2,3) plt.title('equalize_adapthist-'+ str(kernel_sizenuber)+'-0.001') plt.imshow(img2,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,4) plt.hist(arr2, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img3=exposure.equalize_adapthist(img, kernel_size=kernel_sizenuber, clip_limit=0.005, nbins=256) arr3=img3.flatten() plt.subplot(5,2,5) plt.title('equalize_adapthist-'+ str(kernel_sizenuber)+'-0.005') plt.imshow(img3,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,6) plt.hist(arr3, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img4=exposure.equalize_adapthist(img, kernel_size=kernel_sizenuber, clip_limit=0.01, nbins=256) arr4=img4.flatten() plt.subplot(5,2,7) plt.title('equalize_adapthist-'+ str(kernel_sizenuber)+'-0.01') plt.imshow(img4,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,8) plt.hist(arr4, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() img5=exposure.equalize_adapthist(img, kernel_size=kernel_sizenuber, clip_limit=0.05, nbins=256) arr5=img5.flatten() plt.subplot(5,2,9) plt.title('equalize_adapthist-'+ str(kernel_sizenuber)+'-0.05') plt.imshow(img5,plt.cm.gray) #均衡化圖像 plt.subplot(5,2,10) plt.hist(arr5, bins=256, normed=1,edgecolor='None',facecolor='red') #均衡化直方圖 plt.show() 

 

2.png

---

參考文獻
python數字圖像處理

python skimage圖像處理(二)

 

This blog is from: https://www.jianshu.com/p/66e6261f0279 

 

圖像簡單濾波

對圖像進行濾波，可以有兩種效果：一種是平滑濾波，用來抑制噪聲；另一種是微分算子，可以用來檢測邊緣和特征提取。
skimage庫中通過filters模塊進行濾波操作。
1、sobel算子
sobel算子可用來檢測邊緣
函數格式為：

skimage.filters.sobel(image, mask=None) 

from skimage import data,filters import matplotlib.pyplot as plt img = data.camera() edges = filters.sobel(img) plt.imshow(edges,plt.cm.gray) 

 

2、roberts算子

roberts算子和sobel算子一樣，用於檢測邊緣
調用格式也是一樣的：

edges = filters.roberts(img) 

3、scharr算子
功能同sobel，調用格式：

edges = filters.scharr(img) 

4、prewitt算子
功能同sobel，調用格式：

edges = filters.prewitt(img) 

5、canny算子
canny算子也是用於提取邊緣特征，但它不是放在filters模塊，而是放在feature模塊
函數格式：

skimage.feature.canny(image，sigma=1.0) 

可以修改sigma的值來調整效果

from skimage import data,filters,feature import matplotlib.pyplot as plt img = data.camera() edges1 = feature.canny(img) #sigma=1 edges2 = feature.canny(img,sigma=3) #sigma=3 plt.figure('canny',figsize=(8,8)) plt.subplot(121)plt.imshow(edges1,plt.cm.gray) plt.subplot(122)plt.imshow(edges2,plt.cm.gray) plt.show() 

 

 

從結果可以看出，sigma越小，邊緣線條越細小。
6、gabor濾波
gabor濾波可用來進行邊緣檢測和紋理特征提取。
函數調用格式：

 

skimage.filters.gabor_filter(image, frequency) 

通過修改frequency值來調整濾波效果，返回一對邊緣結果，一個是用真實濾波核的濾波結果，一個是想象的濾波核的濾波結果。

from skimage import data,filters import matplotlib.pyplot as plt img = data.camera() filt_real, filt_imag = filters.gabor_filter(img,frequency=0.6) plt.figure('gabor',figsize=(8,8)) plt.subplot(121) plt.title('filt_real') plt.imshow(filt_real,plt.cm.gray) plt.subplot(122) plt.title('filt-imag') plt.imshow(filt_imag,plt.cm.gray) plt.show() 

 

 

以上為frequency=0.6的結果圖。

 

 

 

以上為frequency=0.1的結果圖

7、gaussian濾波
多維的濾波器，是一種平滑濾波，可以消除高斯噪聲。
調用函數為：

skimage.filters.gaussian_filter(image, sigma) 

通過調節sigma的值來調整濾波效果

from skimage import data,filters import matplotlib.pyplot as plt img = data.astronaut() edges1 = filters.gaussian_filter(img,sigma=0.4) #sigma=0.4 edges2 = filters.gaussian_filter(img,sigma=5) #sigma=5 plt.figure('gaussian',figsize=(8,8)) plt.subplot(121) plt.imshow(edges1,plt.cm.gray) plt.subplot(122) plt.imshow(edges2,plt.cm.gray) plt.show() 

 

 

可見sigma越大，過濾后的圖像越模糊
8.median
中值濾波，一種平滑濾波，可以消除噪聲。
需要用skimage.morphology模塊來設置濾波器的形狀。

 

from skimage import data,filters import matplotlib.pyplot as plt from skimage.morphology import disk img = data.camera() edges1 = filters.median(img,disk(5)) edges2= filters.median(img,disk(9)) plt.figure('median',figsize=(8,8)) plt.subplot(121) plt.imshow(edges1,plt.cm.gray) plt.subplot(122) plt.imshow(edges2,plt.cm.gray) plt.show() 

 

 

從結果可以看出，濾波器越大，圖像越模糊。

9、水平、垂直邊緣檢測

上邊所舉的例子都是進行全部邊緣檢測，有些時候我們只需要檢測水平邊緣，或垂直邊緣，就可用下面的方法。

水平邊緣檢測：sobel_h, prewitt_h, scharr_h 垂直邊緣檢測： sobel_v, prewitt_v, scharr_v 

from skimage import data,filters import matplotlib.pyplot as plt img = data.camera() edges1 = filters.sobel_h(img) edges2 = filters.sobel_v(img) plt.figure('sobel_v_h',figsize=(8,8)) plt.subplot(121) plt.imshow(edges1,plt.cm.gray) plt.subplot(122) plt.imshow(edges2,plt.cm.gray) plt.show() 

 

 

上邊左圖為檢測出的水平邊緣，右圖為檢測出的垂直邊緣。

10、交叉邊緣檢測

可使用Roberts的十字交叉核來進行過濾，以達到檢測交叉邊緣的目的。這些交叉邊緣實際上是梯度在某個方向上的一個分量。
其中一個核：

0 1 -1 0 

對應的函數：

roberts_neg_diag(image） 

例：

from skimage import data,filters import matplotlib.pyplot as plt img =data.camera() dst =filters.roberts_neg_diag(img) plt.figure('filters',figsize=(8,8)) plt.subplot(121)plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(dst,plt.cm.gray) 

 

 

另外一個核：

1 0 0 -1 

對應函數為：

roberts_pos_diag(image） 

from skimage import data,filters import matplotlib.pyplot as plt img =data.camera() dst =filters.roberts_pos_diag(img) plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(dst,plt.cm.gray) 

 

---

圖像自動閾值分割

圖像閾值分割是一種廣泛應用的分割技術，利用圖像中要提取的目標區域與其背景在灰度特性上的差異，把圖像看作具有不同灰度級的兩類區域(目標區域和背景區域)的組合，選取一個比較合理的閾值，以確定圖像中每個像素點應該屬於目標區域還是背景區域，從而產生相應的二值圖像。
在skimage庫中，閾值分割的功能是放在filters模塊中。
我們可以手動指定一個閾值，從而來實現分割。也可以讓系統自動生成一個閾值，下面幾種方法就是用來自動生成閾值。

1、threshold_otsu
基於Otsu的閾值分割方法，函數調用格式：

skimage.filters.threshold_otsu(image, nbins=256) 

參數image是指灰度圖像，返回一個閾值。

from skimage import data,filters import matplotlib.pyplot as plt image = data.camera() thresh = filters.threshold_otsu(image) #返回一個閾值 dst =(image <= thresh)*1.0 #根據閾值進行分割 plt.figure('thresh',figsize=(8,8)) plt.subplot(121) plt.title('original image') plt.imshow(image,plt.cm.gray) plt.subplot(122) plt.title('binary image') plt.imshow(dst,plt.cm.gray) plt.show() 

返回閾值為87，根據87進行分割得下圖:

 

 

 

2、threshold_yen
使用方法同上：

thresh = filters.threshold_yen(image) 

返回閾值為198，分割如下圖：

 

 

 

3、threshold_li
使用方法同上：

thresh = filters.threshold_li(image) 

返回閾值64.5，分割如下圖：

 

 

 

4、threshold_isodata
閾值計算方法：

threshold = (image[image <= threshold].mean() +image[image > threshold].mean()) / 2.0 

使用方法同上：

thresh = filters.threshold_isodata(image) 

返回閾值為87，因此分割效果和threshold_otsu一樣。
5、threshold_adaptive
調用函數為：

skimage.filters.threshold_adaptive(image, block_size, method='gaussian'） 

block_size: 塊大小，指當前像素的相鄰區域大小，一般是奇數（如3，5，7。。。）
method: 用來確定自適應閾值的方法，有'mean', 'generic', 'gaussian' 和 'median'。省略時默認為gaussian
該函數直接訪問一個閾值后的圖像，而不是閾值。

from skimage import data,filters import matplotlib.pyplot as plt image = data.camera() dst =filters.threshold_adaptive(image, 15) #返回一個閾值圖像 plt.figure('thresh',figsize=(8,8)) plt.subplot(121) plt.title('original image') plt.imshow(image,plt.cm.gray) plt.subplot(122) plt.title('binary image') plt.imshow(dst,plt.cm.gray) plt.show() 

 

 

大家可以修改block_size的大小和method值來查看更多的效果。如：

dst1 =filters.threshold_adaptive(image,31,'mean') dst2 =filters.threshold_adaptive(image,5,'median') 

兩種效果如下：

 

 

---

基本圖形的繪制

圖形包括線條、圓形、橢圓形、多邊形等。
在skimage包中，繪制圖形用的是draw模塊，不要和繪制圖像搞混了。
1、畫線條
函數調用格式為：

skimage.draw.line(r1,c1,r2,c2) 

r1,r2: 開始點的行數和結束點的行數
c1,c2: 開始點的列數和結束點的列數
返回當前繪制圖形上所有點的坐標，如：

rr, cc =draw.line(1, 5, 8, 2) 

表示從（1，5）到（8，2）連一條線，返回線上所有的像素點坐標[rr,cc]

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc =draw.line(1, 150, 470, 450) img[rr, cc] =255 plt.imshow(img,plt.cm.gray) 

 

 

如果想畫其它顏色的線條，則可以使用set_color（）函數，格式為：

skimage.draw.set_color(img, coords, color) 

例：

draw.set_color(img,[rr,cc],[255,0,0]) 

則繪制紅色線條。

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc =draw.line(1, 150, 270, 250) draw.set_color(img,[rr,cc],[0,0,255]) plt.imshow(img,plt.cm.gray) 

 

2、畫圓
函數格式：

skimage.draw.circle(cy, cx, radius） 

cy和cx表示圓心點，radius表示半徑

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc=draw.circle(150,150,50) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

 

3、多邊形
函數格式：

skimage.draw.polygon(Y,X) 

Y為多邊形頂點的行集合，X為各頂點的列值集合。

from skimage import draw,data import matplotlib.pyplot as plt import numpy as np img=data.chelsea() Y=np.array([10,10,60,60]) X=np.array([200,400,400,200]) rr, cc=draw.polygon(Y,X) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

 

我在此處只設置了四個頂點，因此是個四邊形。
4、橢圓
格式：

skimage.draw.ellipse(cy, cx, yradius, xradius） 

cy和cx為中心點坐標，yradius和xradius代表長短軸。

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc=draw.ellipse(150, 150, 30, 80) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

 

5、貝塞兒曲線
格式：

skimage.draw.bezier_curve(y1,x1,y2,x2,y3,x3,weight) 

y1,x1表示第一個控制點坐標
y2,x2表示第二個控制點坐標
y3,x3表示第三個控制點坐標
weight表示中間控制點的權重，用於控制曲線的彎曲度。

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc=draw.bezier_curve(150,50,50,280,260,400,2) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

 

6、畫空心圓
和前面的畫圓是一樣的，只是前面是實心圓，而此處畫空心圓，只有邊框線。
格式：

skimage.draw.circle_perimeter(yx,yc,radius) 

yx,yc是圓心坐標，radius是半徑

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc=draw.circle_perimeter(150,150,50) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

 

7、空心橢圓
格式：

skimage.draw.ellipse_perimeter(cy, cx, yradius, xradius） 

cy,cx表示圓心
yradius,xradius表示長短軸

from skimage import draw,data import matplotlib.pyplot as plt img=data.chelsea() rr, cc=draw.ellipse_perimeter(150, 150, 30, 80) draw.set_color(img,[rr,cc],[255,0,0]) plt.imshow(img,plt.cm.gray) 

 

---

基本形態學濾波

對圖像進行形態學變換。變換對象一般為灰度圖或二值圖，功能函數放在morphology子模塊內。

1、膨脹（dilation)

原理：一般對二值圖像進行操作。找到像素值為1的點，將它的鄰近像素點都設置成這個值。1值表示白，0值表示黑，因此膨脹操作可以擴大白色值范圍，壓縮黑色值范圍。一般用來擴充邊緣或填充小的孔洞。
功能函數：

skimage.morphology.dilation(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。

from skimage import data import skimage.morphology as sm import matplotlib.pyplot as plt img=data.checkerboard() dst1=sm.dilation(img,sm.square(5)) #用邊長為5的正方形濾波器進行膨脹濾波 dst2=sm.dilation(img,sm.square(15)) #用邊長為15的正方形濾波器進行膨脹濾波 plt.figure('morphology',figsize=(8,8)) plt.subplot(131) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(132) plt.title('morphological image') plt.imshow(dst1,plt.cm.gray) plt.subplot(133) plt.title('morphological image') plt.imshow(dst2,plt.cm.gray) 

分別用邊長為5或15的正方形濾波器對棋盤圖片進行膨脹操作，結果如下：

 

 

 

可見濾波器的大小，對操作結果的影響非常大。一般設置為奇數。
除了正方形的濾波器外，濾波器的形狀還有一些，現列舉如下：

morphology.square: 正方形 morphology.disk: 平面圓形 morphology.ball: 球形 morphology.cube: 立方體形 morphology.diamond: 鑽石形 morphology.rectangle: 矩形 morphology.star: 星形 morphology.octagon: 八角形 morphology.octahedron： 八面體 

注意，如果處理圖像為二值圖像（只有0和1兩個值），則可以調用：

skimage.morphology.binary_dilation(image, selem=None） 

用此函數比處理灰度圖像要快。

2、腐蝕（erosion)

函數：

skimage.morphology.erosion(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。
和膨脹相反的操作，將0值擴充到鄰近像素。擴大黑色部分，減小白色部分。可用來提取骨干信息，去掉毛刺，去掉孤立的像素。

from skimage import data import skimage.morphology as sm import matplotlib.pyplot as plt img=data.checkerboard() dst1=sm.erosion(img,sm.square(5)) #用邊長為5的正方形濾波器進行膨脹濾波 dst2=sm.erosion(img,sm.square(25)) #用邊長為25的正方形濾波器進行膨脹濾波 plt.figure('morphology',figsize=(8,8)) plt.subplot(131) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(132) plt.title('morphological image') plt.imshow(dst1,plt.cm.gray) plt.subplot(133) plt.title('morphological image') plt.imshow(dst2,plt.cm.gray) 

 

 

注意，如果處理圖像為二值圖像（只有0和1兩個值），則可以調用：

skimage.morphology.binary_erosion(image, selem=None） 

用此函數比處理灰度圖像要快。
3、開運算（opening)
函數：

skimage.morphology.openning(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。
先腐蝕再膨脹，可以消除小物體或小斑塊。

from skimage import io,color import skimage.morphology as sm import matplotlib.pyplot as plt img=color.rgb2gray(io.imread('d:/pic/mor.png')) dst=sm.opening(img,sm.disk(9)) #用邊長為9的圓形濾波器進行膨脹濾波 plt.figure('morphology',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.axis('off') plt.subplot(122) plt.title('morphological image') plt.imshow(dst,plt.cm.gray) plt.axis('off') 

 

 

注意，如果處理圖像為二值圖像（只有0和1兩個值），則可以調用：

skimage.morphology.binary_opening(image, selem=None） 

用此函數比處理灰度圖像要快。
4、閉運算（closing)
函數：

skimage.morphology.closing(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。
先膨脹再腐蝕，可用來填充孔洞。

from skimage import io,color import skimage.morphology as sm import matplotlib.pyplot as plt img=color.rgb2gray(io.imread('d:/pic/mor.png')) dst=sm.closing(img,sm.disk(9)) #用邊長為5的圓形濾波器進行膨脹濾波 plt.figure('morphology',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.axis('off') plt.subplot(122) plt.title('morphological image') plt.imshow(dst,plt.cm.gray) plt.axis('off') 

 

 

注意，如果處理圖像為二值圖像（只有0和1兩個值），則可以調用：

skimage.morphology.binary_closing(image, selem=None） 

用此函數比處理灰度圖像要快。
5、白帽（white-tophat)
函數：

skimage.morphology.white_tophat(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。
將原圖像減去它的開運算值，返回比結構化元素小的白點

from skimage import io,color import skimage.morphology as sm import matplotlib.pyplot as plt img=color.rgb2gray(io.imread('d:/pic/mor.png')) dst=sm.white_tophat(img,sm.square(21)) plt.figure('morphology',figsize=(8,8)) plt.subplot(121)plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.axis('off') plt.subplot(122) plt.title('morphological image') plt.imshow(dst,plt.cm.gray) plt.axis('off') 

 

 

6、黑帽（black-tophat)
函數：

skimage.morphology.black_tophat(image, selem=None） 

selem表示結構元素，用於設定局部區域的形狀和大小。
將原圖像減去它的閉運算值，返回比結構化元素小的黑點，且將這些黑點反色。

from skimage import io,color import skimage.morphology as sm import matplotlib.pyplot as plt img=color.rgb2gray(io.imread('d:/pic/mor.png')) dst=sm.black_tophat(img,sm.square(21)) plt.figure('morphology',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.axis('off') plt.subplot(122) plt.title('morphological image') plt.imshow(dst,plt.cm.gray) plt.axis('off') 

 

---

高級濾波

本文提供更多更強大的濾波方法，這些方法放在filters.rank子模塊內。
這些方法需要用戶自己設定濾波器的形狀和大小，因此需要導入morphology模塊來設定。

1、autolevel
這個詞在photoshop里面翻譯成自動色階，用局部直方圖來對圖片進行濾波分級。
該濾波器局部地拉伸灰度像素值的直方圖，以覆蓋整個像素值范圍。
格式：

skimage.filters.rank.autolevel(image, selem） 

selem表示結構化元素，用於設定濾波器。

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) auto =sfr.autolevel(img, disk(5)) #半徑為5的圓形濾波器 plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(auto,plt.cm.gray) 

 

 

2、bottomhat 與 tophat
bottomhat: 此濾波器先計算圖像的形態學閉運算，然后用原圖像減去運算的結果值，有點像黑帽操作。

 

bophat: 此濾波器先計算圖像的形態學開運算，然后用原圖像減去運算的結果值，有點像白帽操作。
格式：

skimage.filters.rank.bottomhat(image, selem） skimage.filters.rank.tophat(image, selem） 

selem表示結構化元素，用於設定濾波器。
下面是bottomhat濾波的例子：

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) auto =sfr.bottomhat(img, disk(5)) #半徑為5的圓形濾波器 plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(auto,plt.cm.gray) 

 

3、enhance_contrast

對比度增強。求出局部區域的最大值和最小值，然后看當前點像素值最接近最大值還是最小值，然后替換為最大值或最小值。
函數：

 enhance_contrast(image, selem） 

selem表示結構化元素，用於設定濾波器。

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) auto =sfr.enhance_contrast(img, disk(5)) #半徑為5的圓形濾波器plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(auto,plt.cm.gray) 

 

4、entropy

求局部熵，熵是使用基為2的對數運算出來的。該函數將局部區域的灰度值分布進行二進制編碼，返回編碼的最小值。
函數格式：

entropy(image, selem） 

selem表示結構化元素，用於設定濾波器。

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) dst =sfr.entropy(img, disk(5)) #半徑為5的圓形濾波器 plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(dst,plt.cm.gray) 

 

 

5、equalize
均衡化濾波。利用局部直方圖對圖像進行均衡化濾波。
函數格式：

 

equalize(image, selem） 

selem表示結構化元素，用於設定濾波器。

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) dst =sfr.equalize(img, disk(5)) #半徑為5的圓形濾波器 plt.figure('filters',figsize=(8,8)) plt.subplot(121) plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122) plt.title('filted image') plt.imshow(dst,plt.cm.gray) 

 

 

6、gradient
返回圖像的局部梯度值（如：最大值-最小值），用此梯度值代替區域內所有像素值。
函數格式：

gradient(image, selem） 

selem表示結構化元素，用於設定濾波器。

from skimage import data,color import matplotlib.pyplot as plt from skimage.morphology import disk import skimage.filters.rank as sfr img =color.rgb2gray(data.lena()) dst =sfr.gradient(img, disk(5)) #半徑為5的圓形濾波器 plt.figure('filters',figsize=(8,8)) plt.subplot(121)plt.title('origin image') plt.imshow(img,plt.cm.gray) plt.subplot(122)plt.title('filted image') plt.imshow(dst,plt.cm.gray) 

 

 

7、其它濾波器
濾波方式很多，下面不再一一詳細講解，僅給出核心代碼，所有的函數調用方式都是一樣的。
最大值濾波器（maximum):返回圖像局部區域的最大值，用此最大值代替該區域內所有像素值。

dst =sfr.maximum(img, disk(5)) 

最小值濾波器（minimum)：返回圖像局部區域內的最小值，用此最小值取代該區域內所有像素值。

dst =sfr.minimum(img, disk(5)) 

均值濾波器（mean) : 返回圖像局部區域內的均值，用此均值取代該區域內所有像素值。

dst =sfr.mean(img, disk(5)) 

中值濾波器（median): 返回圖像局部區域內的中值，用此中值取代該區域內所有像素值。

dst =sfr.median(img, disk(5)) 

莫代爾濾波器（modal) : 返回圖像局部區域內的modal值，用此值取代該區域內所有像素值。

dst =sfr.modal(img, disk(5)) 

otsu閾值濾波（otsu): 返回圖像局部區域內的otsu閾值，用此值取代該區域內所有像素值。

dst =sfr.otsu(img, disk(5)) 

閾值濾波（threshhold): 將圖像局部區域中的每個像素值與均值比較，大於則賦值為1，小於賦值為0，得到一個二值圖像。

dst =sfr.threshold(img, disk(5)) 

減均值濾波（subtract_mean): 將局部區域中的每一個像素，減去該區域中的均值。

dst =sfr.subtract_mean(img, disk(5)) 

求和濾波（sum) :求局部區域的像素總和，用此值取代該區域內所有像素值。

dst =sfr.sum(img, disk(5)) 

參考文獻
python數字圖像處理

python skimage圖像處理(三)

 

This blog is from: https://www.jianshu.com/p/7693222523c0 

 

霍夫線變換

在圖片處理中，霍夫變換主要是用來檢測圖片中的幾何形狀，包括直線、圓、橢圓等。
在skimage中，霍夫變換是放在tranform模塊內，本篇主要講解霍夫線變換。
對於平面中的一條直線，在笛卡爾坐標系中，可用y=mx+b來表示，其中m為斜率，b為截距。但是如果直線是一條垂直線，則m為無窮大，所有通常我們在另一坐標系中表示直線，即極坐標系下的r=xcos(theta)+ysin(theta)。即可用（r,theta）來表示一條直線。其中r為該直線到原點的距離，theta為該直線的垂線與x軸的夾角。如下圖所示。

 

 

 

對於一個給定的點（x0,y0), 我們在極坐標下繪出所有通過它的直線（r,theta)，將得到一條正弦曲線。如果將圖片中的所有非0點的正弦曲線都繪制出來，則會存在一些交點。所有經過這個交點的正弦曲線，說明都擁有同樣的(r,theta), 意味着這些點在一條直線上。

 

 

發上圖所示，三個點(對應圖中的三條正弦曲線）在一條直線上，因為這三個曲線交於一點，具有相同的（r, theta)。霍夫線變換就是利用這種方法來尋找圖中的直線。
函數：

skimage.transform.hough_line(img) 

返回三個值：
h: 霍夫變換累積器
theta: 點與x軸的夾角集合，一般為0-179度
distance: 點到原點的距離，即上面的所說的r.
例：

import skimage.transform as st import numpy as np import matplotlib.pyplot as plt #matplotlib inline # 構建測試圖片 image = np.zeros((100, 100)) #背景圖 idx = np.arange(25, 75) #25-74序列 image[idx[::-1], idx] = 255 # 線條\ image[idx, idx] = 255 # 線條/ # hough線變換 h, theta, d = st.hough_line(image) #生成一個一行兩列的窗口（可顯示兩張圖片）. fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 6)) plt.tight_layout() #顯示原始圖片 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示hough變換所得數據 ax1.imshow(np.log(1 + h)) ax1.set_title('Hough transform') ax1.set_xlabel('Angles (degrees)') ax1.set_ylabel('Distance (pixels)') ax1.axis('image') plt.show() 

 

 

從右邊那張圖可以看出，有兩個交點，說明原圖像中有兩條直線。
如果我們要把圖中的兩條直線繪制出來，則需要用到另外一個函數：

skimage.transform.hough_line_peaks(hspace, angles, dists） 

用這個函數可以取出峰值點，即交點，也即原圖中的直線。
返回的參數與輸入的參數一樣。我們修改一下上邊的程序，在原圖中將兩直線繪制出來。

import skimage.transform as st import numpy as np import matplotlib.pyplot as plt # 構建測試圖片 image = np.zeros((100, 100)) #背景圖 idx = np.arange(25, 75) #25-74序列 image[idx[::-1], idx] = 255 # 線條\ image[idx, idx] = 255 # 線條/ # hough線變換 h, theta, d = st.hough_line(image) #生成一個一行三列的窗口（可顯示三張圖片）. fig, (ax0, ax1,ax2) = plt.subplots(1, 3, figsize=(8, 6)) plt.tight_layout() #顯示原始圖片 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示hough變換所得數據 ax1.imshow(np.log(1 + h)) ax1.set_title('Hough transform') ax1.set_xlabel('Angles (degrees)') ax1.set_ylabel('Distance (pixels)') ax1.axis('image') #顯示檢測出的線條 ax2.imshow(image, plt.cm.gray) row1, col1 = image.shape for _, angle, dist in zip(*st.hough_line_peaks(h, theta, d)): y0 = (dist - 0 * np.cos(angle)) / np.sin(angle) y1 = (dist - col1 * np.cos(angle)) / np.sin(angle) ax2.plot((0, col1), (y0, y1), '-r') ax2.axis((0, col1, row1, 0)) ax2.set_title('Detected lines') ax2.set_axis_off() plt.show() 

 

 

注意，繪制線條的時候，要從極坐標轉換為笛卡爾坐標，公式為：

 

 

skimage還提供了另外一個檢測直線的霍夫變換函數，
概率霍夫線變換：

skimage.transform.probabilistic_hough_line(img, threshold=10, line_length=5,line_gap=3) 

參數：

img: 待檢測的圖像。 threshold： 閾值，可先項，默認為10 line_length: 檢測的最短線條長度，默認為50 line_gap: 線條間的最大間隙。增大這個值可以合並破碎的線條。默認為10 

返回：

lines: 線條列表, 格式如((x0, y0), (x1, y0))，標明開始點和結束點。 

下面，我們用canny算子提取邊緣，然后檢測哪些邊緣是直線？

import skimage.transform as st import matplotlib.pyplot as plt from skimage import data,feature #使用Probabilistic Hough Transform. image = data.camera() edges = feature.canny(image, sigma=2, low_threshold=1, high_threshold=25) lines = st.probabilistic_hough_line(edges, threshold=10, line_length=5,line_gap=3) print(len(lines)) # 創建顯示窗口. fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6)) plt.tight_layout() #顯示原圖像 ax0.imshow(image, plt.cm.gray) ax0.set_title('Input image') ax0.set_axis_off() #顯示canny邊緣 ax1.imshow(edges, plt.cm.gray) ax1.set_title('Canny edges') ax1.set_axis_off() #用plot繪制出所有的直線 ax2.imshow(edges * 0) for line in lines: p0, p1 = line ax2.plot((p0[0], p1[0]), (p0[1], p1[1])) row2, col2 = image.shape ax2.axis((0, col2, row2, 0)) ax2.set_title('Probabilistic Hough') ax2.set_axis_off() plt.show() 

 

---

在極坐標中，圓的表示方式為：

x=x0+rcosθ y=y0+rsinθ 

圓心為(x0,y0),r為半徑，θ為旋轉度數，值范圍為0-359
如果給定圓心點和半徑，則其它點是否在圓上，我們就能檢測出來了。在圖像中，我們將每個非0像素點作為圓心點，以一定的半徑進行檢測，如果有一個點在圓上，我們就對這個圓心累加一次。如果檢測到一個圓，那么這個圓心點就累加到最大，成為峰值。因此，在檢測結果中，一個峰值點，就對應一個圓心點。
霍夫圓檢測的函數：

skimage.transform.hough_circle(image, radius) 

radius是一個數組，表示半徑的集合，如[3，4，5，6]
返回一個3維的數組（radius index, M, N), 第一維表示半徑的索引，后面兩維表示圖像的尺寸。
例1：繪制兩個圓形，用霍夫圓變換將它們檢測出來。

import numpy as np import matplotlib.pyplot as plt from skimage import draw,transform,feature img = np.zeros((250, 250,3), dtype=np.uint8) rr, cc = draw.circle_perimeter(60, 60, 50) #以半徑50畫一個圓 rr1, cc1 = draw.circle_perimeter(150, 150, 60) #以半徑60畫一個圓 img[cc, rr,:] =255 img[cc1, rr1,:] =255 fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5)) ax0.imshow(img) #顯示原圖 ax0.set_title('origin image') hough_radii = np.arange(50, 80, 5) #半徑范圍 hough_res =transform.hough_circle(img[:,:,0], hough_radii) #圓變換  centers = [] #保存所有圓心點坐標 accums = [] #累積值 radii = [] #半徑 for radius, h in zip(hough_radii, hough_res): #每一個半徑值，取出其中兩個圓 num_peaks = 2 peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峰值 centers.extend(peaks) accums.extend(h[peaks[:, 0], peaks[:, 1]]) radii.extend([radius] * num_peaks) #畫出最接近的圓 image =np.copy(img) for idx in np.argsort(accums)[::-1][:2]: center_x, center_y = centers[idx] radius = radii[idx] cx, cy =draw.circle_perimeter(center_y, center_x, radius) image[cy, cx] =(255,0,0) ax1.imshow(image) ax1.set_title('detected image') 

結果圖如下：原圖中的圓用白色繪制，檢測出的圓用紅色繪制。

 

 

 

例2，檢測出下圖中存在的硬幣。

 

 

import numpy as np import matplotlib.pyplot as plt from skimage import data, color,draw,transform,feature,util image = util.img_as_ubyte(data.coins()[0:95, 70:370]) #裁剪原圖片 edges =feature.canny(image, sigma=3, low_threshold=10, high_threshold=50) #檢測canny邊緣 fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5)) ax0.imshow(edges, cmap=plt.cm.gray) #顯示canny邊緣 ax0.set_title('original iamge') hough_radii = np.arange(15, 30, 2) #半徑范圍 hough_res =transform.hough_circle(edges, hough_radii) #圓變換  centers = [] #保存中心點坐標 accums = [] #累積值 radii = [] #半徑 for radius, h in zip(hough_radii, hough_res): #每一個半徑值，取出其中兩個圓 num_peaks = 2 peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峰值 centers.extend(peaks) accums.extend(h[peaks[:, 0], peaks[:, 1]]) radii.extend([radius] * num_peaks) #畫出最接近的5個圓 image = color.gray2rgb(image) for idx in np.argsort(accums)[::-1][:5]: center_x, center_y = centers[idx] radius = radii[idx] cx, cy =draw.circle_perimeter(center_y, center_x, radius) image[cy, cx] = (255,0,0) ax1.imshow(image) ax1.set_title('detected image') 

 

橢圓變換是類似的，使用函數為：

skimage.transform.hough_ellipse(img,accuracy, threshold, min_size, max_size) 

輸入參數：

img: 待檢測圖像。 accuracy: 使用在累加器上的短軸二進制尺寸，是一個double型的值，默認為1 thresh: 累加器閾值，默認為4 min_size: 長軸最小長度，默認為4 max_size: 短軸最大長度，默認為None,表示圖片最短邊的一半。 

返回一個 [(accumulator, y0, x0, a, b, orientation)] 數組，accumulator表示累加器，（y0,x0)表示橢圓中心點，（a,b)分別表示長短軸，orientation表示橢圓方向

例：檢測出咖啡圖片中的橢圓杯口

import matplotlib.pyplot as plt from skimage import data,draw,color,transform,feature #加載圖片，轉換成灰度圖並檢測邊緣 image_rgb = data.coffee()[0:220, 160:420] #裁剪原圖像，不然速度非常慢 image_gray = color.rgb2gray(image_rgb) edges = feature.canny(image_gray, sigma=2.0, low_threshold=0.55, high_threshold=0.8) #執行橢圓變換 result =transform.hough_ellipse(edges, accuracy=20, threshold=250,min_size=100, max_size=120) result.sort(order='accumulator') #根據累加器排序 #估計橢圓參數 best = list(result[-1]) #排完序后取最后一個 yc, xc, a, b = [int(round(x)) for x in best[1:5]] orientation = best[5] #在原圖上畫出橢圓 cy, cx =draw.ellipse_perimeter(yc, xc, a, b, orientation) image_rgb[cy, cx] = (0, 0, 255) #在原圖中用藍色表示檢測出的橢圓 #分別用白色表示canny邊緣，用紅色表示檢測出的橢圓，進行對比 edges = color.gray2rgb(edges) edges[cy, cx] = (250, 0, 0) fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4)) ax1.set_title('Original picture') ax1.imshow(image_rgb) ax2.set_title('Edge (white) and result (red)') ax2.imshow(edges) plt.show() 

 

 

霍夫橢圓變換速度非常慢，應避免圖像太大。

---

 

在前面的python數字圖像處理（10）：圖像簡單濾波 中，我們已經講解了很多算子用來檢測邊緣，其中用得最多的canny算子邊緣檢測。
本篇我們講解一些其它方法來檢測輪廓。
1、查找輪廓（find_contours)
measure模塊中的find_contours()函數，可用來檢測二值圖像的邊緣輪廓。
函數原型為：

skimage.measure.find_contours(array, level) 

array: 一個二值數組圖像
level: 在圖像中查找輪廓的級別值
返回輪廓列表集合，可用for循環取出每一條輪廓。
例1：

import numpy as np import matplotlib.pyplot as plt from skimage import measure,draw #生成二值測試圖像 img=np.zeros([100,100]) img[20:40,60:80]=1 #矩形 rr,cc=draw.circle(60,60,10) #小圓 rr1,cc1=draw.circle(20,30,15) #大圓 img[rr,cc]=1 img[rr1,cc1]=1 #檢測所有圖形的輪廓 contours = measure.find_contours(img, 0.5) #繪制輪廓 fig, (ax0,ax1) = plt.subplots(1,2,figsize=(8,8)) ax0.imshow(img,plt.cm.gray) ax1.imshow(img,plt.cm.gray) for n, contour in enumerate(contours): ax1.plot(contour[:, 1], contour[:, 0], linewidth=2) ax1.axis('image') ax1.set_xticks([]) ax1.set_yticks([]) plt.show() 

結果如下：不同的輪廓用不同的顏色顯示

 

 

 

例2：

import matplotlib.pyplot as plt from skimage import measure,data,color #生成二值測試圖像 img=color.rgb2gray(data.horse()) #檢測所有圖形的輪廓 contours = measure.find_contours(img, 0.5) #繪制輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(img,plt.cm.gray) ax0.set_title('original image') rows,cols=img.shape ax1.axis([0,rows,cols,0]) for n, contour in enumerate(contours): ax1.plot(contour[:, 1], contour[:, 0], linewidth=2) ax1.axis('image') ax1.set_title('contours') plt.show() 

 

2、逼近多邊形曲線
逼近多邊形曲線有兩個函數：subdivide_polygon（)和 approximate_polygon（）
subdivide_polygon（)采用B樣條（B-Splines)來細分多邊形的曲線，該曲線通常在凸包線的內部。
函數格式為：

skimage.measure.subdivide_polygon(coords, degree=2, preserve_ends=False) 

coords: 坐標點序列。
degree: B樣條的度數，默認為2
preserve_ends: 如果曲線為非閉合曲線，是否保存開始和結束點坐標，默認為false
返回細分為的坐標點序列。
approximate_polygon（）是基於Douglas-Peucker算法的一種近似曲線模擬。它根據指定的容忍值來近似一條多邊形曲線鏈，該曲線也在凸包線的內部。
函數格式為:

skimage.measure.approximate_polygon(coords, tolerance) 

coords: 坐標點序列
tolerance: 容忍值
返回近似的多邊形曲線坐標序列。
例：

import numpy as np import matplotlib.pyplot as plt from skimage import measure,data,color #生成二值測試圖像 hand = np.array([[1.64516129, 1.16145833], [1.64516129, 1.59375], [1.35080645, 1.921875], [1.375, 2.18229167], [1.68548387, 1.9375], [1.60887097, 2.55208333], [1.68548387, 2.69791667], [1.76209677, 2.56770833], [1.83064516, 1.97395833], [1.89516129, 2.75], [1.9516129, 2.84895833], [2.01209677, 2.76041667], [1.99193548, 1.99479167], [2.11290323, 2.63020833], [2.2016129, 2.734375], [2.25403226, 2.60416667], [2.14919355, 1.953125], [2.30645161, 2.36979167], [2.39112903, 2.36979167], [2.41532258, 2.1875], [2.1733871, 1.703125], [2.07782258, 1.16666667]]) #檢測所有圖形的輪廓 new_hand = hand.copy() for _ in range(5): new_hand =measure.subdivide_polygon(new_hand, degree=2) # approximate subdivided polygon with Douglas-Peucker algorithm appr_hand =measure.approximate_polygon(new_hand, tolerance=0.02) print("Number of coordinates:", len(hand), len(new_hand), len(appr_hand)) fig, axes= plt.subplots(2,2, figsize=(9, 8)) ax0,ax1,ax2,ax3=axes.ravel() ax0.plot(hand[:, 0], hand[:, 1],'r') ax0.set_title('original hand') ax1.plot(new_hand[:, 0], new_hand[:, 1],'g') ax1.set_title('subdivide_polygon') ax2.plot(appr_hand[:, 0], appr_hand[:, 1],'b') ax2.set_title('approximate_polygon') ax3.plot(hand[:, 0], hand[:, 1],'r') ax3.plot(new_hand[:, 0], new_hand[:, 1],'g') ax3.plot(appr_hand[:, 0], appr_hand[:, 1],'b') ax3.set_title('all') 

 

---

高級形態學處理

形態學處理，除了最基本的膨脹、腐蝕、開/閉運算、黑/白帽處理外，還有一些更高級的運用，如凸包，連通區域標記，刪除小塊區域等。
1、凸包
凸包是指一個凸多邊形，這個凸多邊形將圖片中所有的白色像素點都包含在內。
函數為：

skimage.morphology.convex_hull_image(image) 

輸入為二值圖像，輸出一個邏輯二值圖像。在凸包內的點為True, 否則為False
例：

import matplotlib.pyplot as plt from skimage import data,color,morphology #生成二值測試圖像 img=color.rgb2gray(data.horse()) img=(img<0.5)*1 chull = morphology.convex_hull_image(img) #繪制輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(img,plt.cm.gray) ax0.set_title('original image') ax1.imshow(chull,plt.cm.gray) ax1.set_title('convex_hull image') 

 

 

convex_hull_image()是將圖片中的所有目標看作一個整體，因此計算出來只有一個最小凸多邊形。如果圖中有多個目標物體，每一個物體需要計算一個最小凸多邊形，則需要使用convex_hull_object（）函數。
函數格式：

skimage.morphology.convex_hull_object(image, neighbors=8) 

輸入參數image是一個二值圖像，neighbors表示是采用4連通還是8連通，默認為8連通。
例：

import matplotlib.pyplot as plt from skimage import data,color,morphology,feature #生成二值測試圖像 img=color.rgb2gray(data.coins()) #檢測canny邊緣,得到二值圖片 edgs=feature.canny(img, sigma=3, low_threshold=10, high_threshold=50) chull = morphology.convex_hull_object(edgs) #繪制輪廓 fig, axes = plt.subplots(1,2,figsize=(8,8)) ax0, ax1= axes.ravel() ax0.imshow(edgs,plt.cm.gray) ax0.set_title('many objects') ax1.imshow(chull,plt.cm.gray) ax1.set_title('convex_hull image') plt.show() 

 

 

2、連通區域標記
在二值圖像中，如果兩個像素點相鄰且值相同（同為0或同為1），那么就認為這兩個像素點在一個相互連通的區域內。而同一個連通區域的所有像素點，都用同一個數值來進行標記，這個過程就叫連通區域標記。在判斷兩個像素是否相鄰時，我們通常采用4連通或8連通判斷。在圖像中，最小的單位是像素，每個像素周圍有8個鄰接像素，常見的鄰接關系有2種：4鄰接與8鄰接。4鄰接一共4個點，即上下左右，如下左圖所示。8鄰接的點一共有8個，包括了對角線位置的點，如下右圖所示。

 

 

 

在skimage包中，我們采用measure子模塊下的label（）函數來實現連通區域標記。
函數格式：

skimage.measure.label（image,connectivity=None) 

參數中的image表示需要處理的二值圖像，connectivity表示連接的模式，1代表4鄰接，2代表8鄰接。
輸出一個標記數組（labels), 從0開始標記。

import numpy as np import scipy.ndimage as ndi from skimage import measure,color import matplotlib.pyplot as plt #編寫一個函數來生成原始二值圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] #生成網絡 mask = np.zeros((l, l)) generator = np.random.RandomState(1) #隨機數種子 points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) #高斯濾波 return mask > mask.mean() data = microstructure(l=128)*1 #生成測試圖片 labels=measure.label(data,connectivity=2) #8連通區域標記 dst=color.label2rgb(labels) #根據不同的標記顯示不同的顏色 print('regions number:',labels.max()+1) #顯示連通區域塊數(從0開始標記) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, plt.cm.gray, interpolation='nearest') ax1.axis('off') ax2.imshow(dst,interpolation='nearest') ax2.axis('off') fig.tight_layout() plt.show() 

在代碼中，有些地方乘以1，則可以將bool數組快速地轉換為int數組。
結果如圖：有10個連通的區域，標記為0-9

 

 

 

如果想分別對每一個連通區域進行操作，比如計算面積、外接矩形、凸包面積等，則需要調用measure子模塊的regionprops（）函數。該函數格式為：

skimage.measure.regionprops(label_image) 

返回所有連通區塊的屬性列表，常用的屬性列表如下表：

屬性名稱 類型 描述 area int 區域內像素點總數 bbox tuple 邊界外接框(min_row, min_col, max_row, max_col) centroid array　　 質心坐標 convex_area int 凸包內像素點總數 convex_image ndarray 和邊界外接框同大小的凸包　　 coords ndarray 區域內像素點坐標 Eccentricity float 離心率 equivalent_diameter float 和區域面積相同的圓的直徑 euler_number int　　 區域歐拉數 extent float 區域面積和邊界外接框面積的比率 filled_area int 區域和外接框之間填充的像素點總數 perimeter float 區域周長 label int 區域標記 

3、刪除小塊區域
有些時候，我們只需要一些大塊區域，那些零散的、小塊的區域，我們就需要刪除掉，則可以使用morphology子模塊的remove_small_objects（)函數。
函數格式：

skimage.morphology.remove_small_objects(ar, min_size=64, connectivity=1, in_place=False) 

參數：

ar: 待操作的bool型數組。 min_size: 最小連通區域尺寸，小於該尺寸的都將被刪除。默認為64. connectivity: 鄰接模式，1表示4鄰接，2表示8鄰接 in_place: bool型值，如果為True,表示直接在輸入圖像中刪除小塊區域，否則進行復制后再刪除。默認為False. 返回刪除了小塊區域的二值圖像。 

import numpy as np import scipy.ndimage as ndi from skimage import morphology import matplotlib.pyplot as plt #編寫一個函數來生成原始二值圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] #生成網絡 mask = np.zeros((l, l)) generator = np.random.RandomState(1) #隨機數種子 points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) #高斯濾波 return mask > mask.mean() data = microstructure(l=128) #生成測試圖片 dst=morphology.remove_small_objects(data,min_size=300,connectivity=1) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, plt.cm.gray, interpolation='nearest') ax2.imshow(dst,plt.cm.gray,interpolation='nearest') fig.tight_layout() plt.show() 

在此例中，我們將面積小於300的小塊區域刪除（由1變為0），結果如下圖：

 

 

Keep the labels with 2 largest areas.

 label_image =measure.label(newpr) areas = [r.area for r in regionprops(label_image)] areas.sort() if len(areas) > 2: for region in regionprops(label_image): if region.area < areas[-2]: for coordinates in region.coords: label_image[coordinates[0], coordinates[1]] = 0 binary = label_image > 0 label_image = morphology.remove_small_holes(binary, areas[-2]) 

 areas = [r.area for r in regionprops(label_image)] areas.sort() if len(areas) > 2: for region in regionprops(label_image): if region.area < areas[-2]: for coordinates in region.coords: label_image[coordinates[0], coordinates[1]] = 0 binary = label_image > 0 if plot == True: plots[3].axis('off') plots[3].imshow(binary, cmap=plt.cm.bone) 

https://github.com/vivek14632/LungCancerProject/blob/master/preprocessing2.py

4、綜合示例：閾值分割+閉運算+連通區域標記+刪除小區塊+分色顯示

import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as mpatches from skimage import data,filters,segmentation,measure,morphology,color #加載並裁剪硬幣圖片 image = data.coins()[50:-50, 50:-50] thresh =filters.threshold_otsu(image) #閾值分割 bw =morphology.closing(image > thresh, morphology.square(3)) #閉運算 cleared = bw.copy() #復制 cleared=segmentation.clear_border(cleared) #清除與邊界相連的目標物 cleared = sm.opening(cleared,sm.disk(2)) cleared = sm.closing(cleared,sm.disk(2)) label_image =measure.label(cleared) #連通區域標記 borders = np.logical_xor(bw, cleared) #異或 label_image[borders] = -1 image_label_overlay =color.label2rgb(label_image, image=image) #不同標記用不同顏色顯示 fig,(ax0,ax1)= plt.subplots(1,2, figsize=(8, 6)) ax0.imshow(cleared,plt.cm.gray) ax1.imshow(image_label_overlay) for region in measure.regionprops(label_image): #循環得到每一個連通區域屬性集 #忽略小區域 if region.area < 100: continue #繪制外包矩形 minr, minc, maxr, maxc = region.bbox rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr, fill=False, edgecolor='red', linewidth=2) ax1.add_patch(rect) fig.tight_layout() plt.show() 

 

---

骨架提取與分水嶺算法

骨架提取與分水嶺算法也屬於形態學處理范疇，都放在morphology子模塊內。
1、骨架提取
骨架提取，也叫二值圖像細化。這種算法能將一個連通區域細化成一個像素的寬度，用於特征提取和目標拓撲表示。
morphology子模塊提供了兩個函數用於骨架提取，分別是Skeletonize（）函數和medial_axis（）函數。我們先來看Skeletonize（）函數。
格式為：s

kimage.morphology.skeletonize(image) 

輸入和輸出都是一幅二值圖像。
例1：

from skimage import morphology,draw import numpy as np import matplotlib.pyplot as plt #創建一個二值圖像用於測試 image = np.zeros((400, 400)) #生成目標對象1(白色U型) image[10:-10, 10:100] = 1 image[-100:-10, 10:-10] = 1 image[10:-10, -100:-10] = 1 #生成目標對象2（X型） rs, cs = draw.line(250, 150, 10, 280) for i in range(10): image[rs + i, cs] = 1 rs, cs = draw.line(10, 150, 250, 280) for i in range(20): image[rs + i, cs] = 1 #生成目標對象3（O型） ir, ic = np.indices(image.shape) circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2 circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2 image[circle1] = 1 image[circle2] = 0 #實施骨架算法 skeleton =morphology.skeletonize(image) #顯示結果 fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4)) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') ax1.set_title('original', fontsize=20) ax2.imshow(skeleton, cmap=plt.cm.gray) ax2.axis('off') ax2.set_title('skeleton', fontsize=20) fig.tight_layout() plt.show() 

生成一幅測試圖像，上面有三個目標對象，分別進行骨架提取，結果如下：

 

 

 

例2：利用系統自帶的馬圖片進行骨架提取

from skimage import morphology,data,color import matplotlib.pyplot as plt image=color.rgb2gray(data.horse()) image=1-image #反相 #實施骨架算法 skeleton =morphology.skeletonize(image) #顯示結果 fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(8, 4)) ax1.imshow(image, cmap=plt.cm.gray) ax1.axis('off') ax1.set_title('original', fontsize=20) ax2.imshow(skeleton, cmap=plt.cm.gray) ax2.axis('off') ax2.set_title('skeleton', fontsize=20) fig.tight_layout() plt.show() 

 

 

medial_axis就是中軸的意思，利用中軸變換方法計算前景（1值）目標對象的寬度，格式為：

skimage.morphology.medial_axis(image, mask=None, return_distance=False) 

mask: 掩模。默認為None, 如果給定一個掩模，則在掩模內的像素值才執行骨架算法。
return_distance: bool型值，默認為False. 如果為True, 則除了返回骨架，還將距離變換值也同時返回。這里的距離指的是中軸線上的所有點與背景點的距離。

import numpy as np import scipy.ndimage as ndi from skimage import morphology import matplotlib.pyplot as plt #編寫一個函數，生成測試圖像 def microstructure(l=256): n = 5 x, y = np.ogrid[0:l, 0:l] mask = np.zeros((l, l)) generator = np.random.RandomState(1) points = l * generator.rand(2, n**2) mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 mask = ndi.gaussian_filter(mask, sigma=l/(4.*n)) return mask > mask.mean() data = microstructure(l=64) #生成測試圖像 #計算中軸和距離變換值 skel, distance =morphology.medial_axis(data, return_distance=True) #中軸上的點到背景像素點的距離 dist_on_skel = distance * skel fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) ax1.imshow(data, cmap=plt.cm.gray, interpolation='nearest') #用光譜色顯示中軸 ax2.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest') ax2.contour(data, [0.5], colors='w') #顯示輪廓線 fig.tight_layout() plt.show() 

 

 

2、分水嶺算法
分水嶺在地理學上就是指一個山脊，水通常會沿着山脊的兩邊流向不同的“匯水盆”。分水嶺算法是一種用於圖像分割的經典算法，是基於拓撲理論的數學形態學的分割方法。如果圖像中的目標物體是連在一起的，則分割起來會更困難，分水嶺算法經常用於處理這類問題，通常會取得比較好的效果。
分水嶺算法可以和距離變換結合，尋找“匯水盆地”和“分水嶺界限”，從而對圖像進行分割。二值圖像的距離變換就是每一個像素點到最近非零值像素點的距離，我們可以使用scipy包來計算距離變換。
在下面的例子中，需要將兩個重疊的圓分開。我們先計算圓上的這些白色像素點到黑色背景像素點的距離變換，選出距離變換中的最大值作為初始標記點（如果是反色的話，則是取最小值），從這些標記點開始的兩個匯水盆越集越大，最后相交於分山嶺。從分山嶺處斷開，我們就得到了兩個分離的圓。
例1：基於距離變換的分山嶺圖像分割

import numpy as np import matplotlib.pyplot as plt from scipy import ndimage as ndi from skimage import morphology,feature #創建兩個帶有重疊圓的圖像 x, y = np.indices((80, 80)) x1, y1, x2, y2 = 28, 28, 44, 52 r1, r2 = 16, 20 mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2 mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2 image = np.logical_or(mask_circle1, mask_circle2) #現在我們用分水嶺算法分離兩個圓 distance = ndi.distance_transform_edt(image) #距離變換 local_maxi =feature.peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), labels=image) #尋找峰值 markers = ndi.label(local_maxi)[0] #初始標記點 labels =morphology.watershed(-distance, markers, mask=image) #基於距離變換的分水嶺算法 fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 8)) axes = axes.ravel() ax0, ax1, ax2, ax3 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') ax0.set_title("Original") ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest') ax1.set_title("Distance") ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest') ax2.set_title("Markers") ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest') ax3.set_title("Segmented") for ax in axes: ax.axis('off') fig.tight_layout() plt.show() 

 

 

分水嶺算法也可以和梯度相結合，來實現圖像分割。一般梯度圖像在邊緣處有較高的像素值，而在其它地方則有較低的像素值，理想情況 下，分山嶺恰好在邊緣。因此，我們可以根據梯度來尋找分山嶺。
例2：基於梯度的分水嶺圖像分割

import matplotlib.pyplot as plt from scipy import ndimage as ndi from skimage import morphology,color,data,filter image =color.rgb2gray(data.camera()) denoised = filter.rank.median(image, morphology.disk(2)) #過濾噪聲 #將梯度值低於10的作為開始標記點 markers = filter.rank.gradient(denoised, morphology.disk(5)) <10 markers = ndi.label(markers)[0] gradient = filter.rank.gradient(denoised, morphology.disk(2)) #計算梯度 labels =morphology.watershed(gradient, markers, mask=image) #基於梯度的分水嶺算法 fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(6, 6)) axes = axes.ravel() ax0, ax1, ax2, ax3 = axes ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest') ax0.set_title("Original") ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest') ax1.set_title("Gradient") ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest') ax2.set_title("Markers") ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest') ax3.set_title("Segmented") for ax in axes: ax.axis('off') fig.tight_layout() plt.show() 

 

---

參考文獻
python數字圖像處理