2017/10/23

深度學習引發之深度偏頭痛




若「機器學習」或「AI」這類的課程能在高中或大一當個基礎前導學科來上,相信學生們比較不會有自己所學是否有用的疑慮了!相反的,若您像筆者一樣回過頭來學習「Deep Learning」之類的深層神經網路,有可能會責怪自己為何當初線性代數、微積分、向量梯度最佳化方法等等沒能把他學好呢?而這些原本以為找不到出口的方法與數學工具,在此它們都被充分地運用上了!



秒懂感知器(Perceptron
假設我們發明一個科普兒童玩具,該裝置包含了兩個閘刀式開關「輸入」(ON會點亮LED燈,反之OFF燈滅)與一個「神經元」以旋鈕式可變電阻連接到兩個輸入端與一個偏權值,其中當電阻調整至最大時(標示為0)相當於斷路,而電阻調整至最小時(標示為1)相當於導通。此時給小朋友一個任務:給定一OR閘的輸入/出特性,試著去調整神經元鍵值的大小(旋鈕式可變電阻)並根據四種可能的輸入(亮燈與否)與該特定的輸出將神經網路訓練成OR閘。




相信不用多久,小朋友就能調整出一個設定值使得該神經網路具有OR閘的運算能力,雖然這個系統可能有無限多種解。同理,該神經元可以被訓練成NOR閘、AND閘與NAND閘。沒錯,小朋友在調整這些神經元鍵值最後找到解的過程就好比如下圖中找到一條線,使得該線能成功將輸出分成兩類(亮燈與否),因此任何一條線都是一組解。而其中偏權值使該線不通過座標原點。


若回到60年前「感知器(Perceptron」被提出來的年代,以當時的時空背景揣摩,或許可以感受到為何線性代數大師們會對其嗤之以鼻吧?假設x為二維座標空間所張成之向量,則w∙x=0,表向量w=[1 1]與向量x=[x1 x2]正交(內積為0),w為其法向量。所有符合w∙x=0.5的解,在x座標空間的展開,為法向量為w且通過點(0.5, 0.5)之直線。



對輸出值套用所謂的活化函數(activate function,例如sigmoid),即可以定義出正/負類別。



注意,此時若沒有sigmoid函數(activate function)對y轉換,我們找到的只是一條線,無法判定類別。wx愈像則其內積為正,反之其內積為負,透過sigmoid我們可以將輸出判定為正類或是負類別,因此該系統與線性規劃是等價的但仍屬於窮舉。



多元感知器可以比擬成:由多個前級神經元電訊號的加總來決定是否激勵下一級神經元的處理過程,如下圖。神經網路理論的發明若是如下圖中由下往上類推的過程,你可能會覺得它並沒甚麼了不起。但若我們說它是從上到下的思考推論過程(說個漂亮的故事),媒體或許就會覺得它很聰明很厲害(把它給神話了)。

(Source: Wikipedia)



此外,眼尖的同學應該很快發現單層的感知器根本無法分類像XOR這種樣本分隔兩地的問題。



這類問題須要層疊感知器來處理(增加hidden layer),相當於多個線性規劃的聯集可以解凸邊形(convex)的問題。因此,透過更深的層疊感知器裏所當然可以處理更複雜的問題,例如樣本間彼此分隔多區域(disjoint convex)的分類問題。



但該系統仍屬於窮舉,因此,在20幾年前的時空背景下被具有非常漂亮數學基底的SVMSupport Vector Machine打趴是不覺得意外的。
  


梯度法Gradient Descent
顧名思義,透過微分取得曲面的斜率並往該方向逐步找到區域最佳解,否則神經網路的學習過程與方式將跟隨機的窮舉法沒甚麼兩樣。

相信很多人在撰寫類神經網路程式碼時會跟筆者一樣,剛開始可能有點不太適應,因為許多運算需要跳脫原本「迴圈計算」的思維模式並進入「向量處理」思維模式的腦。以三維曲面f(x0,x1)=x0^2+x1^2為例,下面Python程式碼可以把它畫出來。

# f(x0,x1)=x0**2+x1**2
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import Axes3D

x0=np.arange(-3,3.5,0.5)
x1=np.arange(-3,3.5,0.5)
X,Y=np.meshgrid(x0,x1)
Z=X**2+Y**2

ax=plt.axes(projection='3d')
#ax.plot_surface(X,Y,Z,cmap='coolwarm')
#ax.plot_wireframe(X,Y,Z)
ax.plot_trisurf(X.flatten(),Y.flatten(),Z.flatten())

plt.xlabel('x0')
plt.ylabel('x1')
plt.show()




對某參數取偏微分的程式碼很簡單,對某一參數的微量變化取均值,其餘變數以定值帶入求斜率(向量梯度)。



改用向量的寫法如下:

# f(x0,x1)=x0**2+x1**2
def function(V):
    return (V[0]**2+V[1]**2)

def numerical_gradient(f, V):
    h=1e-4
    grad=np.zeros_like(V)
    grad[0]=(f([V[0]+h,V[1]])-f([V[0]-h,V[1]]))/(2*h) # gradient x0
    grad[1]=(f([V[0],V[1]+h])-f([V[0],V[1]-h]))/(2*h) # gradient x1
    return grad

def gradient_descent(f, init_v, lr=0.01, step_num=20):
    v=init_v
    hist=[]
    for i in range(step_num):
        hist.append(v.copy())
        grad=numerical_gradient(f, v)
        v-= lr*grad
    return v,np.array(hist)


之後,我們對該曲面某個範圍內的所有格子點做偏微分以計算其梯度。

# gradient
x0=np.arange(-3,3.5,0.5)
x1=np.arange(-3,3.5,0.5)
X,Y=np.meshgrid(x0,x1)
X=X.flatten()
Y=Y.flatten()

grad=numerical_gradient(function,[X,Y])
plt.quiver(X,Y,-grad[0],-grad[1],angles="xy",color="#666666")
plt.xlabel('x0')
plt.ylabel('x1')
plt.grid()



下圖顯示以梯度法逼近曲面區域最低點(微分為零)的收斂過程:



# grdient descent
init_v=np.array([-3.0,2.7])
r,hist=gradient_descent(function,init_v,lr=0.2)
plt.plot(hist[:,0],hist[:,1],'--o')
plt.show()




更深度的學習
到此,眼尖的同學應該會發現:雖然梯度法給了尋求最佳化過程一個指標與方向性,但仍然沒有比窮舉高明多少?它相當的沒有效率,直到「倒傳遞(backpropagation」理論被提出來。若讀者有興趣可以更深入的參考「史丹佛CS231n」的網路學習課程,我們將主動重新拾起課本,看看甚麼是「chain-rule」?原來「微積分」這麼實用且有趣!


2017/9/17

機器人(或AI)創作



筆者在半導體沉睡了五分之一世紀,沒想到「哇」這一覺醒來整個世界都變了!最近被CNNRNNGANDeep MindDeep LearningDeep DreamDeep StyleDeep ArtDeep CreationDeep TutorDeep 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(例如Haartap-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=%.2f" % 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係數並丟棄最低頻資訊,如此再重建回來的影像會發生什麼事呢?答案是,我們將得到一幅風格迥異像是用木頭雕刻的圖像。這就是GANwork的原因,透過濾波器間的關聯性分析,可以將一個圖像轉換成另一風格的圖像。所以很多做影像處理的老骨頭們,最近都轉行做AI了!



套用「梵谷」的筆觸(特徵紋理)對濾波器係數做破壞性重建的影像,筆者較偷懶,若有事先做關聯性分析效果會較好。






這不是AI!?
嗯,空拍機電量不足,只能飛二十分鐘。

放眼世界,許多AI晶片陸續發表,而筆者嚴重脫隊了,但願您沒缺席!