DigiVFX hw2: Image stitching 自動縫圖

第二個作業真是可怕，從來沒有寫過這麼難寫的程式，課程網頁中也是時間最長的作業 (1個月)。難怪當初有人是因為這個作業寫不出來就退選了…

聽了課才知道自動縫圖技術成熟是最近的事情。以前的方法是整張圖去套用影像轉換然後找誤差最小的轉換。後來的方法是先對圖找特徵點，然後每張圖之前就用特徵點去找最小誤差的轉換，這樣子運算量就小很多了，而已如果特徵點找的好的話會比用整張圖下去算效果來的好。

所以好的找特徵點的演算法就是這次的重點。

上課教了Harris和SIFT兩種，目前表現最好的SIFT演算法是2004年才發明的。找特徵點最基本的原理就是算影像的gradient然後找區域極值，當dx dy都變化很大，那他可能是特徵點的機會就很大。而SIFT多加的概念是對影像做Gaussian Blur，然後同樣算gradient找極值，當dx dy在這麼糊的情況下都還是極值的話，那就表示這個特徵點很大（相當於看的範圍比較大）。而找特徵點演算法除了縫圖也用在很多地方，大至上要辨識圖的內容的相關演算法都可以扯上關系。因為是個很重要的影像相關演算法分類，所以就選表現比較強的SIFT來寫看看。

整個SIFT的流程大概是：

1. 套各種不同sigma大小的Gaussian Blur，然後做到一個地步，就把圖片長寬縮小二分之一以減少運算量，然後繼續做下去。
2. 把第1步做出來的Image兩兩相臨的影像相減，求得Differential of Gaussian Image。
3. 在第3步的DoG圖上求Gradient，找極值做為候選的特徵點。
4. 用一些內插法把不夠強的極值刪掉。
5. 對每個剩下的特徵點求出各別的一個向量來描述這個特徵點的大小和方向。
6. 對每個特徵點求出描述特徵點周遭的128維descriptor向量。

個人覺得SIFT真的是每一步都很難寫

1. 一開始光是要決定要算多少次就縮小影像就花了很多時間，投影片和課程網頁提供的matlab code和paper說的都不太一樣。對圖算Gaussian Blur是個很慢的事情，一開始用自己寫套濾鏡算法真是慢到爆。而且Gaussian Blur用的sigma很難搞懂到底是要怎麼算，而這又是很重要的部份。
2. 在用內插法刪較弱的極值的方法很難看懂 x_x。
3. 最後在算descriptor的真的看不懂怎麼算。投影片講的很簡單，實際做有一堆問題，要怎麼選多大的區塊來算descriptor就難倒我了。
4. 整個演算法裡面有非常多的magic number，什麼數字該怎麼設原作者都幫你tune好了，但是我用他tune好的數字結果很爛…就不知道是哪裡做錯了。
5. 慢…自己寫的程式真是又慢又爛…原作者的程式算不到30秒我算要快10分鐘…

後來在網路上找到siftpp這個有open source的SIFT實作，速度和原作者的SIFT demo相差不遠。不懂的就看source code，看不懂的就copy過來先跳過。最後做出來的結果：

找特徵點，看起來還可以

找每個特徵點的向量，這個就怪怪的…

至於算descriptor的部份無法用看的猜測正確與否，也讓我不知道該如何debug。

當寫完前半段就大概一個多禮拜了，而這作業還沒有完，還有後面的縫圖的部份。

縫圖的大概流程是：

1. 對每張圖算各自的feature point descriptor
2. 兩兩去比對特徵點，找到很像的特徵就當他們match
3. 每兩張圖會有很多組match的特徵點，用這個去算圖片的轉換矩陣
4. 圖片轉換後縫起來，兩兩之間要做漸層轉換，然後還有一些校正。
5. 把一張很長的圖投影到圈柱座標上，然後用一個特竟的全景圖撥放程式去看。

找feature match的部份課程網頁提供了一個C++ kd-tree library可以找。原本的用意是每個特徵點去和kd-tree裡的全部特徵點比對找出最接近的4個點。不過後來用了用才發現，我原本的descriptor是存float [128]，而這個library支援的型態是double，要改型態的話還要改source code…，C++寫的library竟然沒有用template寫…
每個特徵點是影像座標(x,y)加上128維的descriptor，而每個kd-tree的node只有放descriptor去比較距離，後來找不到好的方法可以單純從我的float descriptor找回kd-tree找到的double descriptor原本的(x,y)座標，所以也沒有用。就用最單純的linear search找距離最短的2個點，如果距離ratio < 0.8就算是找到feature match。

話說投影片有提到RANSAC演算法，用途是從一堆sample得到的點回推原本的model是什麼，就算sample到的點有很多不合這個model的也ok。這個演算法我不知道要用在哪裡？我想可能是用在原本的找4個點吧？

接著找到的幾組match的feature座標，有線性代數公式可以算出轉換矩陣。這裡又是用Lapack++這套library，又讓我遇到一個很奇怪的設計。要算一個反矩陣有一個函式LaLUInverseIP ( matrix, pivot vector) 可以用，但是我一用就當機。後來仔細看官網說明才發現，要用 LaLUInverseIP( matrix, pivot vector)之前要先用 LUFactorizeIP(matrix, pivot vector)算過一次才可以算出結果。這就奇怪了，既然參數都一樣，為什麼沒有一個函式LaInverse(matrix)可以直接算？無法理解…

最後是是圖片轉換，一開始不知道怎麼解進圖片轉換超出原本圖的大小的問題。後來想到的方法是，用原本的4個點先丟進轉換矩陣，求出目標地需要的大小，然就就再加個座標offset就好了。不過又是同樣的問題，自己寫的內插真的是有…夠…慢…

到了測試逢圖的時候，先試試只比較兩張圖然後縫起來：

 

結果左邊的圖轉換後變這樣 = =
 

重新debug的過程中發現，是matching的地方錯了，所以求出來的轉換矩陣也錯了，不過求matching是跟據descriptor來的，而descriptor又是SIFT最不懂的地方，也不知道該如何debug。

當然還是有方法可以改的，像是用Harris的演算法找特徵點，Harris的演算法就簡單很多，對原圖的gradient每個點都可以得到一個矩陣，算一下矩陣的trace和det就可以決定這是不是特徵點了。不過這個方法的問題是，每個點沒有descriptor，我也不知道只有點後來要怎麼做matching…

所以這次的hw2就到此為止了！每天真的是起床就在電腦前寫到晚上睡覺，寫了快100小時，結果最後爛掉了 = =，還是設個停損點好了。現在想想，當初老師說沒有課程壓力就會容易半途而廢，好像真的是這樣… >_<

後來我在網路上找到有修這這門課的人，竟然可以用matlab寫…我還以為資工系的課都要求用C/C++寫…害我寫的時候每次要debug都很麻煩，一錯就要全部重算，又是要等個幾分鐘…

還有寫這次作業中間用library遇到一堆問題，回去翻C++ Primer，複習了運算子超載，STL algorithm和iterator用法，繼承，模版。重看的時候發現，以前一堆東西都不知道為什麼是要這樣做，上學期去台大資工聽了一學期的物件導向程式設計，這次再看就清楚多了。下次來看看深入淺出設計模式看會不會比較有感覺一點…