筆者在半導體沉睡了五分之一世紀，沒想到「哇」這一覺醒來整個世界都變了！最近被CNN、RNN、GAN、Deep Mind、Deep Learning、Deep Dream、Deep Style、Deep Art、Deep Creation、Deep Tutor、Deep Zzzz … 搞得昏頭轉向，過程雖然非常令人振奮（感到有趣）但同時也非常令人憂慮。心想：「完蛋了！我竟然真空了，這些年來我到底在幹嘛？」若您從未缺席就先恭喜你囉，您肯定是未來的希望！（筆者只幫TSMC搞過Deep Sleep，以前商品要潮要取名叫i-xx，現在要改成deep-xx才能順應時代呢）
初次看到「CNN」與「特徵萃取」，筆者瞬間聯想到了二十幾年前舊石器時代用來直接在離散時間信號做處理的「短脈衝響應濾波器」，利如：可以容易達成多解析度信號分解的「小波轉換（Wavelet Transform）」。當然，也就聯想到了「Lena」，而這張廣泛被用來在學術論文間PK的標準間接告訴你為什麼國內在信號處理的基礎研究上一直處於落後狀態！沒錯，點一下鏈結就知道了，是作學問的「格局」與看問題的「角度」！（指導教授應該更嚴格要求你，因為你的高度不夠：正確的高度是2318，不是512）
知識遊牧民族在忙碌於知識追趕的過程，就好像追著「韓劇」或各個「四爺」們不停惡補的過程，沒完沒了。所以，就先暫時放緩腳步來用Python玩一下以小波做Convolution囉！ 用更輕鬆寬廣的心情請機器人畫素描或進行藝術創作。
AQI預測
我們在「機器人(或AI)的知識擷取」文中有提到如何提取開放資料庫，若給定全台詳細的個觀測站AQI表單資料，然後請某人猜測一個不在此表單內某一城鎮的AQI值。我想人類大腦可能會很直覺地設法將有列在表單內的鄰近測站資料拿出來，然後不是很精確地（人腦會盡可能的偷懶以節省能量消耗）估算一下平均（mean）與大約可能的落差（delta）以做為預測參考吧？而少數落差大者所須記憶量將減少或往往被忽略以推算趨勢，有點像是將這68個有明確值的離散資料做回歸（regression）再拓展回全台的整個尺度（scale）。咦，這不就是很典型的「低通濾波器」嗎？
例如，先將已知數據依照測站所在地理位置做區域平均（選k個鄰近點），再跟據該城鎮的所在的位置（落在哪個區塊）猜測其可能的趨勢，這是一種低通濾波器的行為，之後再加上如麥寮、前鎮、竹山等較少數的嚴重汙染原資料（高通濾波器的特徵）等參考做修正。跟據研究，我們人腦會自然做這樣的處理以減少計算量。
看到下面左邊的圖（先把右邊遮起來），它是經過四階層的高低通濾波器與次採樣最後只留下最低頻32x32像素的結果，您能否只憑這32x32像素想像出是甚麼圖案呢？右邊則是直接只取這32x32像素用低通濾波器重建回來的影像（壓縮比256），有比您預測（想像）的結果好嗎？（應該好很多吧？科學報導有時不用把它看的太神）
旋積/卷積傻傻分不清楚?
提到濾波器，自然要提到最廣泛被使用的時域（或空間域）的「有限脈衝響應濾波器」與「Convolution」。筆者可以這樣分嗎: 翻譯成「旋積」的是老骨頭，說「卷積」的是Google時代所教育的年輕人？（因為我們二十幾年前沒有Google時都翻譯成「迴旋積分」）
可以想像若將DWT用神經網路來做，那麼在第一階神經網路，濾波器係數就好比是神經元鍵結的權重而神經元θ值為濾波器係數與影像信號做「旋積（Convolution）」的結果。經過次取樣可以再傳遞給第二階，如此反覆下去。
而CNN好比是「混合多種濾波器」並「巧妙刻意」破壞掉部份「濾波器係數」與「結果」的系統，我們可以預期第一階的神經元輸出相當於時域（或空間域）的濾波器結果（係數），會萃取出最高頻帶的空間特徵。同理，可以預期在第二階的濾波器結果會萃取出次高頻帶的空間特徵。因此，若輸入信號越能符合由第一階濾波器係數重建的結果，將越能激勵（activate）第一階的神經元。可以預期，由第二階濾波器係數重建的結果將包含許多第一階的空間特徵，而包含該特徵的輸入信號自然越能激勵（activate）第二階的神經元。我們可依此類推到更高階的等效效應，這是多重解析度信號分析與重建的基本原理。
離散小波轉換（Discrete Wavelet Transform，簡稱DWT）
自從Daubechies等人將小波轉換與數位信號處理結合，使小波轉換得以廣泛的應用在數學、醫學、或工程等各領域。無奈，賺錢的卻是寫SnapChat的年輕人？（若請機器人模仿筆者寫作應該蠻簡單的，筆者文章裏有許多酸溜溜的模式Pattern）
如下圖，DWT的架構可拆解為分解（Decomposition）: 將輸入的時-空間域信號經由濾波器分解成高/低頻成分，與合成（Reconstruction）: 將分解的高/低頻信號經由濾波器重新組合。由於DWT能將能量有效地集中於低頻，以及多解析度分解的特性尤其適合被運用在影像處理（包括特徵抽取與圖形辨識等）。
用Python實作二維影像的DWT非常簡單，無須安裝其它套件。我們知道「旋積（convolution）」物理上意謂著用一信號對另一信號做採樣（sampling）的動作，但由於影像輸入信號是有邊界的，在做旋積的過程為了邊界採樣順利通常會作mod運算來簡單達成週期性拓展（periodic extension）。
以Daubechies濾波器家族為例，一般稱tap-t的短脈衝響應濾波器其細數長度為n=2t（例如Haar為tap-1的濾波器其長度為2），我們可以透過巧妙安排暫存器索引來將信號分解過程所造成的相位移調整回來。此外我們對暫存器索引做shift運算以避免除法，在達到次取樣目的同時減少乘法運算。下面程式碼是對二維影像作水平方向的DWT轉換的例子，它沒做任何優化（浪費很多記憶體），讀者可以自行改良它。
#
DWT decomposition and sub sampling
def dwt_decomposition_x(inImage,firset):
firL,firH=firset
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=n-1 #offset
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
data[:,:]=inImage[:,:]
dwt=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,rows): # low pass
for col in range(0,hsize):
for i in range(0,n):
dwt[row,col]+=data[row,((col<<1)-t+i)%cols]*firL[i]
for row in range(0,rows): # high pass
for col in range(0,hsize):
for i in range(0,n):
dwt[row,col+hsize]+=data[row,((col<<1)-t+i)%cols]*firH[i]
return dwt
|
因為次取樣（sub sampling）能量強度會減半，所以升頻重建時須將信號振幅放大兩倍（乘2）。而該網站所提供的濾波器係數已經乘上√2並平均分配到分解與重建濾波器係數內，如此可以大量減少所須的乘法運算，但也因為精度不足的關係所以重建時會有些失真（但PSNR都達200左右，也就是均方差約為1e-16水平，肉眼看不出來）。此外，我們必須反轉原分解濾波器序列以得到重建濾波器序列。
#
DWT up sampling and reconstruction
def dwt_reconstruction_x(inImage,firSet):
firDL,firDH=firSet # reverse the filter coefficient sequency
firL=firDL.copy();
firL.reverse()
firH=firDH.copy(); firH.reverse()
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=0
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
dwtL=np.zeros((rows,cols),dtype=np.float64)
data[:,0::2]=inImage[:,0:hsize] # up sampling low pass to even columns
for row in range(0,rows): # low pass
for col in range(0,cols):
for i in range(0,n):
dwtL[row,col]+=data[row,(col-t+i)%cols]*firL[i]
data[:,0::2]=inImage[:,hsize:] # up sampling high pass to even columns
dwtH=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,rows): # high pass
for col in range(0,cols):
for i in range(0,n):
dwtH[row,col]+=data[row,(col-t+i)%cols]*firH[i]
return (dwtL+dwtH)
|
同理我們可以對垂直方向展開，並反覆對低頻信號做分解動做以達成二維影像的多解析度分解。下面的程式碼包含完整讀取影像、執行多解析度分解與重建，讀者可以自行調整所需分解的階層以及更換濾波器係數。
from urllib.request import urlopen
from matplotlib import pyplot as
plt
import numpy as np
url='http://www-2.cs.cmu.edu/~chuck/lennapg/lena_std.tif'
img=plt.imread(urlopen(url),0)
# lean
plt.set_cmap('gray')
plt.title('Original
Image')
plt.imshow(img)
plt.show()
#
Daubechies Scaling/Wavelet (low/high pass) filter
DB={
# Haar
'L1':[ 0.7071067811865476,0.7071067811865476],
'H1':[-0.7071067811865476,0.7071067811865476],
# Daubechies2 Low/High pass
(Scaling/Wavelet filter)
'L2':[-0.1294095226,0.2241438680,
0.8365163037, 0.4829629131], # scaling filter
'H2':[-0.4829629131,0.8365163037,-0.2241438680,-0.1294095226], # wavelet filter
# Daubechies3 Low/High pass filter
'L3':[0.0352262919,-0.0854412739,-0.1350110200,\
0.4598775021,0.8068915093,0.3326705530],
'H3':[-0.3326705530,0.8068915093,-0.4598775021,\
-0.1350110200,0.0854412739,0.0352262919],
# Daubechies4 Low/High pass filter
'L4':[-0.0105974018,0.0328830117,0.0308413818,-0.1870348117,\
-0.0279837694,0.6308807679,0.7148465706, 0.2303778133],
'H4':[-0.2303778133,0.7148465706,-0.6308807679,-0.0279837694,\
0.1870348117,0.0308413818,-0.0328830117,-0.0105974018],
# Daubechies5 Low/High pass filter
'L5':[0.0033357253,-0.0125807520,-0.0062414902,0.0775714938,-0.0322448696,\
-0.2422948871, 0.1384281459,
0.7243085284,0.6038292698, 0.1601023980],
'H5':[-0.1601023980,
0.6038292698,-0.7243085284,0.1384281459,0.2422948871,\
-0.0322448696,-0.0775714938,-0.0062414902,0.0125807520,0.0033357253],
# Daubechies6 Low/High pass filter
'L6':[-0.0010773011, 0.0047772575, 0.0005538422,-0.0315820393,\
0.0275228655,
0.0975016056,-0.1297668676,-0.2262646940,\
0.3152503517, 0.7511339080,
0.4946238904,0.1115407434],
'H6':[-0.1115407434,
0.4946238904,-0.7511339080,0.3152503517,\
0.2262646940,-0.1297668676,-0.0975016056,0.0275228655,\
0.0315820393,0.0005538422,-0.0047772575,-0.0010773011]
}
####################################################
#
DWT decomposition and sub sampling
def
dwt_decomposition_x(inImage,firset):
firL,firH=firset
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=n-1 # offset
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
data[:,:]=inImage[:,:]
dwt=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,rows): # low pass
for col in range(0,hsize):
for i in range(0,n):
dwt[row,col]+=data[row,((col<<1)-t+i)%cols]*firL[i]
for row in range(0,rows): # high pass
for col in range(0,hsize):
for i in range(0,n):
dwt[row,col+hsize]+=data[row,((col<<1)-t+i)%cols]*firH[i]
return dwt
def
dwt_decomposition_y(inImage,firset):
firL,firH=firset
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=n-1 # offset
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
data[:,:]=inImage[:,:]
dwt=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,hsize): # low pass
for col in range(0,cols):
for i in range(0,n):
dwt[row,col]+=data[((row<<1)-t+i)%rows,col]*firL[i]
for row in range(0,hsize): # high pass
for col in range(0,cols):
for i in range(0,n):
dwt[row+hsize,col]+=data[((row<<1)-t+i)%rows,col]*firH[i]
return dwt
####################################################
#
DWT up sampling and reconstruction
def
dwt_reconstruction_x(inImage,firSet):
firDL,firDH=firSet #
reverse the filter coefficient sequency
firL=firDL.copy();
firL.reverse()
firH=firDH.copy(); firH.reverse()
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=0 # offset
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
dwtL=np.zeros((rows,cols),dtype=np.float64)
data[:,0::2]=inImage[:,0:hsize] # up sampling low pass to even columns
for row in range(0,rows): # low pass
for col in range(0,cols):
for i in range(0,n):
dwtL[row,col]+=data[row,(col-t+i)%cols]*firL[i]
data[:,0::2]=inImage[:,hsize:] # up sampling high pass to even columns
dwtH=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,rows): # high pass
for col in range(0,cols):
for i in range(0,n):
dwtH[row,col]+=data[row,(col-t+i)%cols]*firH[i]
return (dwtL+dwtH)
def dwt_reconstruction_y(inImage,firSet):
firDL,firDH=firSet # reverse the filter coefficient
sequency
firL=firDL.copy(); firL.reverse()
firH=firDH.copy(); firH.reverse()
firL=np.array(firL).astype('float64')
firH=np.array(firH).astype('float64')
n=int(len(firL)) #
tap=n/2: vanishing moment
t=0 # offset
rows,cols=inImage.shape
hsize=int(cols/2)
data=np.zeros((rows,cols),dtype=np.float64)
dwtL=np.zeros((rows,cols),dtype=np.float64)
data[0::2,:]=inImage[0:hsize,:] # up sampling low pass to even rows
for row in range(0,rows): # low pass
for col in range(0,cols):
for i in range(0,n):
dwtL[row,col]+=data[(row-t+i)%cols,col]*firL[i]
data[0::2,:]=inImage[hsize:,:] # up sampling high pass to even rows
dwtH=np.zeros((rows,cols),dtype=np.float64)
for row in range(0,rows): # high pass
for col in range(0,cols):
for i in range(0,n):
dwtH[row,col]+=data[(row-t+i)%cols,col]*firH[i]
return (dwtL+dwtH)
####################################################
#
DWT 2D
def
dwt_2d(img,firSet,level):
rows,cols=img.shape
data=np.zeros((rows,cols))
data[:,:]=img[:,:]
while level>0:
buf=data[:rows,:cols]
dwt=dwt_decomposition_x(buf,firDec)
dwt=dwt_decomposition_y(dwt,firDec)
data[:rows,:cols]=dwt[:rows,:cols]
level-=1
rows>>=1
cols>>=1
return data
#
IDWT 2D
def
idwt_2d(img,firSet,level):
rows,cols=img.shape
data=np.zeros((rows,cols))
data[:,:]=img[:rows,:cols]
rows>>=level
cols>>=level
while level>0:
level-=1
rows<<=1
cols<<=1
idwt=np.zeros((rows,cols))
idwt[:,:]=data[:rows,:cols]
idwt=dwt_reconstruction_x(idwt,firDec)
idwt=dwt_reconstruction_y(idwt,firDec)
data[:rows,:cols]=idwt[:rows,:cols]
return idwt
####################################################
#
compute PSNR
import math
def
calculate_psnr(orgImage,idwtImage):
rows,cols=orgImage.shape
s=0
for row in range(0,rows): # calculate MSE
for col in range(0,cols):
s+=pow((idwtImage[row,col]-orgImage[row,col]),2)
mse=s/(rows*cols)
psnr=10*math.log10(pow(255,2)/(mse+1.0e-60))
print("MSE =%.5e" % mse)
print("PSNR=%
return psnr
####################################################
#
convert to gray color
def
convert_to_gray(img):
if len(img.shape)==3:
rows,cols,layers=img.shape
gray=img[:,:,0]*0.299+img[:,:,1]*0.587+img[:,:,2]*0.114
else:
rows,cols=img.shape
gray=img[:,:]
return gray
gray=convert_to_gray(img)
plt.title('Convert
to Gray')
plt.imshow(gray)
plt.show()
####################################################
#
Daubechies coefficients
firDec=(DB['L3'],DB['H3'])
rows,cols=gray.shape
level=3
csize=cols>>level
#
DWT decomposition
dwt=dwt_2d(gray,firDec,level)
#
IDWT reconstruction
idwt=idwt_2d(dwt,firDec,level)
#
PSNR
calculate_psnr(gray,idwt)
####################################################
#
visualization
tmp=np.zeros((rows,cols))
tmp[:,:]=255-dwt[:,:]*10
tmp[:csize,:csize]=0.5/level*dwt[:csize,:csize]
plt.title('DWT
Decomposition (level={})'.format(level))
plt.imshow(tmp,vmin=0,vmax=255)
plt.show()
plt.title('DWT
Reconstruction (level={})'.format(level))
plt.imshow(idwt)
plt.show()
|
如下圖是以DB3係數做完一階DWT影像轉換的結果，影像被分解為最低頻（左上角方塊）與三個最高頻區塊。右上角區塊篩選出垂直的特徵細節，左下角區塊篩選出水平的特徵細節，而右下角區塊則篩選出最後剩餘對角的特徵細節。
若把Lena低頻的資料全部以房子的低頻資料取代，重建回來的影像會發生甚麼事？（想家的Lena或是倩女幽魂嗎？）
若分解時使用DB3係數而還原時使用DB5係數並丟棄最低頻資訊，如此再重建回來的影像會發生什麼事呢？答案是，我們將得到一幅風格迥異像是用木頭雕刻的圖像。這就是GAN會work的原因，透過濾波器間的關聯性分析，可以將一個圖像轉換成另一風格的圖像。所以很多做影像處理的老骨頭們，最近都轉行做AI了！
這不是AI啊!?
放眼世界，許多AI晶片陸續發表，而筆者嚴重脫隊了，但願您沒缺席！