虛擬實境 建模 SVD 筆記

由於在選讀論文時,和上課都遇到相關問題,但一直對這 SVD 的實際感觸沒很深,了解數學計算,卻不知道特性一直很困惑我。在這裡收集一下曾經學過的內容,當初在線性代數時沒有學到這一塊,至於 $LU$ 分解有沒有上到過,我自己也深感困惑,不過看一下網路資源勉強還能理解。

至於 SVD 則在線性代數非常重要,在許多領域時常會遇到,如巨量資料 (massive data)?機器學習?影像處理?就我所知理解看看吧!

SVD

定義從一般式

$$A_{m, n} = U_{m, m} \Sigma_{m, r} V_{n, r}^T$$

得到縮減過後的相似矩陣

$$A_{m, n} \cong U_{m, r} \Sigma_{r, r} V_{n, r}^T$$

給定 $A$ 矩陣,要變成三個矩陣相乘,可以壓縮儲存空間,一般來說都是拿來得到近似的 $A$ 矩陣,藉由刪除過小的 eigenvalue 來完成,eigenvalue 就是所謂的矩陣特徵值。

在學線性代數時,沒有仔細想過 eigenvalue 的意義所在,每一個 eigenvalue 也會對應一個 eigenvector,這些 eigenvalue 大小對應的分佈的相關性大小,而 eigenvector 則表示新的轉換軸的方向。類似最小平方法得到的結果。

但在巨量資料課程中我們可以知道一件事情,若是稀疏矩陣 $A$,經過 SVD 運作後,得到 $U, V$ 矩陣都是稠密的,因此空間儲存反而划不來,因此會有另外一種解法 CUR,可以上網搜尋 massive data mining coursera 在 這裏 可以知道關於 CUR 的資訊。

從課程中的例子來看,假設在電影評論系統中,則 A 矩陣表示電影對人的喜好度,則 SVD 相當於如下描述 (假想,單純的矩陣概念)

  • $A(i, j)$ 電影編號 i 對用戶 j 的好感度
  • $U(i, j)$ 電影編號 i 對電影屬性 j 的成分程度
  • $\Sigma(i, j)$ 電影屬性 i 跟電影屬性 j 的之間重要性。
  • $V(i, j)$ 表示用戶 i 對電影屬性 j 的偏愛程度。

也就是說有一些電影屬性不是很重要,可以移除掉,這樣的表示法得到一個近似的矩陣 $A$

在巨量資料中一個有趣的運用,可以用在購物推薦系統中,來達到預測相關性,跟某人對其他產品的喜好程度或者是需求性,但是在處理之前會遇到矩陣過大,但又非常稀疏,為了解決紀錄問題,就有奇奇怪怪的 Dimensionality Reduction 來完成,SVD 就是其中一種方法。

其他應用部分,將不重要的特徵值消除,來求得哪些元素是需要的,看到有一篇 文章 在做 PCA 主成份分析。

$$\begin{align} & A_{m, n} V_{n, r} = U_{m, r} \Sigma_{r, r} V_{n, r}^T V_{n, r} \\ & \because V_{n, r}^T V_{n, r} = I \\ & \therefore A_{m, n} V_{n, r} = U_{m, r} \Sigma_{r, r} \\ \end{align}$$

如此一來就可以知道要保留哪幾個主要資訊。

在虛擬實境論文中也是有這樣的處理,做的是 CT 建模,將模型的每個頂點,要計算 QEM (最小二次誤差) 的重要性,可以將其表示成一個矩陣,利用 SVD 方式得到,來決策要保留前 k 個重要性的節點,再將其他節點移除,來讓新得到節點盡可能地可以貼近真實模型。

之所以會這麼做,是要利用內插法來強迫得到邊緣偵測,但這樣會造成很多新頂點,為了要把過多頂點消除,使用 QEM 作為評價方法,利用 SVD 來完成大規模刪除的操作。

Map-Reduce

都明白 SVD 是一個很繁複的做法,有一個 $O(N^3)$ 的算法來求得所有 eigenvalue,若要用一般電腦跑,光是得到一個大矩陣所有 eigenvalue 複雜度相當高,其複雜程度沒有細讀過。有一種迭代的方法,類似馬可夫過程,每一次迭代可以逼近一個 eigenvalue,為了得到 n 個 eigenvalue,想必要迭代非常多次。運用 Map-Reduce 分散到很多台群落電腦加快矩陣乘法,這種方法的複雜度是可以期待的。

詳細操作 請參考

eigenvalue $\lambda$ 滿足

$Av = \lambda v$

現在開始迭代 $v$ 讓其等式趨近於穩定

$A x_{i} = \lambda x_{i+1}$

初始化$x_{i} = (1, 1, ..., 1, 1)$ ,其中$|x_{i}| = 1$ 藉由長度為 1 的向量,來求得$\lambda$,之所以要對 $x$ 進行正規化,其一就是要求出$\lambda$,另一個原因是要看何時穩定,當$|x_{i} - x_{i+1}|$ 差距夠小時,則表示進入穩定狀態。

那就可以得到矩陣 $A$ 的其中一個 eigenvalue,為了求得另外一組解,得到新的$A^* = A - \lambda x x^T$,再對$A^*$ 找到 eigenvalue,不斷地迭代下去。

Read More +

虛擬實境 考古筆記

虛擬實境這一門課可以讓你體驗有如在上另一門課的虛擬實境,前半段講生物結構,來了解軟硬體需要那些特性、效能,後半段在講視覺建模算法,牽涉到計算幾何、線性代數 … 等課程。

在幾何課程中介紹的 Delaunay triangulation 對於建模算法是相當重要,接著在多重解析度部分牽涉到傳輸跟處理速度,在資料結構的儲存上很有趣。為了加速運算效能、漸進式需求,各種轉換壓縮都會跑出來。

複習

  1. (15%)[Definition]
    詳細說明一個虛擬實境 (virtual reality) 系統需要具備那些特性,並說明一個虛擬實境系統要用到哪些軟體技術。

    提供人們能感知的動態環境,視覺、嗅覺、觸覺,並且能與環境互動中得到反饋,得到逼真感覺。需要以下幾個軟體

    1. 影像處理與繪製
    2. 物理引擎
    3. 擷取資料的整合
  2. (15%)[Human perception]
    說明人類 1. 單眼立體視覺的原理,2. 雙眼立體視覺的原理;並3. 比較其優缺點。

    1. 單眼立體視覺靠單一眼球調整焦距,將數張影像交給大腦來整合,這必須藉由對物體認知來補足。
    2. 雙眼立體視覺的原理主要靠左右兩個眼球不同的相位差,來計算立體影像的距離差距。
    3. 視野範圍以及反應時間,單眼必須藉由調整焦距來補足,雙眼可以同時獲取不同相位差的影像。
  3. (15%)[Hardware]
    說明感應手套 (sensing glove) 與力回饋手套 (force-feedback glove) 的原理,並比較他們在應用上的差別(可舉例說明)。

    感應手套是由操作者的反應傳輸給感測器,藉由手套上感應器之間的變化數據得到操作者給予的狀態。力回饋手套則是將物理特性反饋給使用者,藉由手套上的機械,施予使用者力反饋。

    感應手套為使用者操作,力反饋手套則相反由環境操作。

  4. (15%)[Modeling]
    說明五種 3D 模型的表面表示式 (surface representation),並比較這些表示式的優缺點。

    目前共計主要分成五種表示法:

    • regular square grid meshes (RSG)

    • triangulated regular meshes (TRM)

    • triangulated irregular networks (TINs)

    • polygon meshes

    • parametric patches

      RSG 建造簡單,但會有四點不共平面,影響到貼圖繪製上的問題。TRM 解決四點不共平面問題,但與 RSG 問題無法凸顯邊緣。TINs 解決無法凸顯邊緣問題,但需要經過最佳化算法來得到不規則三角網格,儲存方式必須記錄所有的點座標,由於要經過最佳化算法,處理速度會慢上很多。動態更動座標的效果較差。polygon meshes 建造相對困難,多邊形的最小單位是三角形,頂點儲存量比 TINs 少很多,parametric patches 藉由取樣點套入參數,讓表面接近平滑,針對部分訊息的改造,平滑效果造成的反饋比較逼真。

  5. (15%)[Multiresolution modeling]
    說明並比較以小波轉換為基礎的多重解析度模塑 (wavelet-based multiresolution modeling) 技術與 Lindstrom 等人的對稱三角網格 (symmetric triangulated networks) 多重解析度模塑技術的原理及優缺點。

    • wavelet-based multiresolution modeling
      解析度資料量與面積成正比,考慮到與觀察者之間的距離,用另外三張影響轉換到更高解析度,每一次解析度增加兩倍。
    • symmetric triangulated networks
      考慮到地形高地、平坦性對於觀察者的影響,投影到使用者的平面,判斷是否差異過大需要提供更高解析,解析度資料量所需節點數成正比

      相較之下,小波轉換給予的傳輸資料反應會比較慢,同時也會在使用者感受不到差異的地方提供更高解析度,很容易浪費資源,對稱三角網格在資料鏈結會複雜許多,沒有小波轉換的儲存架構來得容易理解,但對稱三角網格反應速度快,支持快速的解析度變換。

  6. (15%)[Multiresolution modeling]
    說明漸進網格 (progressive mesh) 多重解析度模塑與二次誤差 (quadric-error-metric) 多重解析度模塑的相同與差異處,並比較兩者的優缺點。

    漸進網格藉由拔點、縮邊,提供不同解析度,依序拔點,考慮四個因素來決策拔點順序,這四個因素分別為影響到其他點之間的距離、縮邊的長度總和、頂點材質的差異、對邊界影響程度。

    考慮的是從精細模型轉換到粗糙模型的差異變化,挑選變成粗糙的影響最小的先拔除,這樣的計算量會比較大,拔除一點後,鄰近的點的誤差值計算也要跟著重新計算。

    二次誤差類似漸進網格的方式進行拔點縮點,但縮邊完的的結果不會只有三種可能,縮到兩端點的其中一個、或者是中點,提供一個數學計算來找到縮到哪一個地方更好。

    二次誤差考慮的是從端點誤差值得相加,並且產生一個新的頂點,有點霍夫曼編碼的意味在。計算從粗糙模型上的點到精細模型平面的距離平方,與之前的模型 (Progressive mesh - Hoppe) 恰恰相反,並不是考慮精細模型上的點到粗糙模型平面的距離平方。因此對於 QEM 的誤差可以累加。問題是粗糙模型上的頂點是無法得知,屬於未知數,用內插法。QEM 誤差精確度只是約略,速度快但品質沒有 Progressive mesh 來得好。

  1. (15%)[Multiresolution modeling]
    說明 (a) 動態載入 (dynamic loading) 的意義。(b) 動態載入與多重解析度模塑技術整合會產生什麼問題。(c)解決上述問題的三種方法。

    動態載入為在記憶體不足時,將部分資料存放在磁碟中,當程式運行到需要那些資料進行運算時再從磁碟中讀取,這樣的過程稱為動態載入。

    動態載入對於多重解析度塑模來說,主要是記憶體問題以及所需要的解析度跟數據大小關係,解析度度越高,所需要的資料量越多。但是畫面呈現的解析度不一下,要怎麼進行資料讀取,動態載入的數據結構定義則需要特別設計,例如是要 response time 最小、計算複雜度最低,讀取數量最小 … 等。

    假設需要不同解析度可以靠 Discrete layered multiresolution models,每一種解析度彼此獨立儲存,這樣的方法可以提供基礎的不同解析度,但讀取時間相對長,儲存在磁碟的空間需求大,效果呈現是最快的,但是離散結果會有不連續的解析度變化。

    為了解決解析度連續變化問題,可以使用 progressive mesh 來達到,提供以點、線與解析度呈現線性變化。

    若需要不同解析度,可以利用 hierarchical models 來完成畫面中每一處的解析度不同,這是一種資料結構,當需要更高解析度時,則走訪更深的樹節點來印出所需要的畫面。

    評價效果好壞利用 Delaunay triangulation 計算平滑程度。

  2. (15%)[Rendering]
    說明shading的意義,敘述 Gouraud shading 與 Phong shading 的作法,最後比較他們的優缺點。

Read More +

虛擬實境 論文研讀 Part 2

接續上一篇,接下來會有很多數學,自己也是一知無解狀態。

QEM

QEM 方程如下,當找到許多在等值曲面上的 $p_k$ 時,接著要找到一個最小 $v$ 參數為滿足目標方程

$v^{M} \leftarrow \text{arg } \underset{v}{min} \sum_{k} (n^t_{k} (v - p_{k}))^2$

為了找到這個 $v$

$v^M \leftarrow v^M + Q^+ (q - Q v^M)$

其中矩陣$Q = \sum_{k} n_{k} n^{t}_{k}$

使用 SVD ( singular value decomposition ),奇異值分解去找解,這個計算好像在巨量資料中也存在 SVD,這世界太可怕了! )。 在找到 $v^M$ 之後,進行簡化模型的縮點 (減少節點數)。經由這一步,可以更貼近物體的輪廓線,來強化邊緣的存在。

為了減少退化三角形跟三角網格的品質,可以藉由 Delaunay Tetrahedralization 再刷新一次。修改、刪除的點其實並不多於全局 (大多都是邊緣?),那麼可以利用鬆弛的方式進行更新會比全局刷新還來得快速。

實戰

接下來這篇論文的作者,要進行實戰比較 (一起來吸收別人怎麼去驗證、統計、比較!實驗方法也是相當重要的。),查看每一次迭代的是否可以的過程、速度。

誤差評定

利用模擬的模型 (幾何模型構成,拿真物要找到物體表面的方向量有點麻煩),來看梯度方向的準確性,簡單來說就是拿實際結果與實驗結果中的一個點,於表面上的法向量夾角作為誤差大小。

$\varepsilon_{n}(p) = \text{arccos } (n(p)^t m(p))$

梯度估計

即使是誤差,平滑也是好事。梯度越小越平滑。

接著它拿兩個物體 sphere 和 fandisk mesh (球體和一個 ???) 進行測試,看來在那些特殊的地方的誤差。算法可能會在某些特殊邊緣、平面的運行效果不好 (算法擅長於那些特徵計算、不善長什麼,來釐清這一點)。

視覺化時,利用下列二種方式凸顯梯度 central differences、Sobel operator (索貝爾運算就是在影像處理運用中,檢測邊緣、圖像梯度用的)。論文作者接下來提供一的新穎 CT value 的計算公式,與一般傳統的不同,讓這個誤差梯度下降。經由上面兩種方式的梯度表示,發現論文作者提供的這種會再更好一點點。

討論

  • 提供一個更準確的 CT value 的計算方式
  • 動態運作 ODT, QEM 的方式

比較 Grid-Based

拿其他的建造方案,如 Grid-Based Polygonization,進行邊緣銳利比較。使用 Hausdorff distance 比較好壞,當然地 Grid-Based 沒對齊好,基本上贏不過三角化的版本!Grid-Based 取樣點夠密集才能勝過 vary triangle 的建模吧!在大多共平面的狀況下,三角形建模的方式需要頂點數非常少。

時間消耗

由於精準度上的調校,提出的方法比起常規的建模是消耗更多時間,但 CT 掃描數據的時間通常需要 10min 到 1hr,但是提出的方法比掃描時間還快。

也就是說讀取資料很慢,處理速度不用這麼快也沒關係,算起來有種 $O(n)$ 的效能。

空間消耗

由於過程中需要大量的矩陣、附加值於記憶體中,會消耗上 GB 的使用空間來維護操作。但現代電腦的記憶體可以超過 10GB 以上,因此這個缺點可以忽略不計。

實作細節

為了加速運行,特別將物體表面周圍的取樣多、密集一點,對於其他地方的取樣少點,如此一來計算量會少很多。

插值

$Tiso= (f(gi)+f(gj))/2$ 中提到插值法,但插值法有很多種,如 tricubic interpolation、bicubic interpolation … 等,如果選擇別種的插植法,效果會不會比較好呢?經過它們的梯度估算方法,他們所提到的插值法會稍微好一些。

結論

提出一個可以用 CT 圖去建模的方法,仰賴的不是三維空間 volumetric data,並且可以提供一個更高精度的建模方法,用少許的三角形就能構造出相當銳利的結果。

  1. 以上都是建立在單一素材上,如果由複合的素材構成,可能是一個無法解決的問題。
  2. 對於光束敏感的物體,CT 的建造方法可以支持更好。
  3. 如果一開始的四面體不夠完善,有可能無法捕抓到邊緣。自適應的 meshing 也許能解決這問題。

提出一個速度慢、空間多的算法,但是現代電腦記憶體夠,速度慢沒關係,比輸入還快就行了

Read More +

虛擬實境 論文研讀 Part 1

期末上繳報告啊,不知道要寫多少才對得起自己。理解正確性不敢保證,相當可怕。

單詞概要

  • Computed Tomography - 電腦斷層掃描 (X 光),簡稱 CT
  • Shepp-Logan filter - 影像去噪中,一種運算操作
  • Delaunay - 計算幾何中 Delaunay Triangulation
  • ODT - optimal Delaunay Triangulation
  • Tetrahedron - 四面體
  • Isosufrface - 等價曲面
  • Truncate - 切去頭端
  • Hausdorff distance - 空間中,集合之間的距離描述 wiki
  • Marching cubes - 用於 CT 模型,重建時的建構 3D 時,根據取樣點得到的最小構造單位。
  • QEM - Quadric Error Mactrics

本文

The Sinogram Polygonizer for Reconstructing 3D Shapes 這篇論文主要拿 Grid-base 掃描重建方法以及根據取樣點重建 3D 模型的算法探討。

首先從離散取樣中重建,最基礎的做法是根據 z-index 作切割,得到 xy 平面,在其上作 Grid-base 網格狀的取樣,得到座標 (x, y, z) 是否在該物體內部。

當取樣數量越多時,密集程度越高,能造就的邊緣就越明顯。藉由 Grid-base 取樣的缺點在於邊緣的凸顯仍然不夠明確於三角化,取樣點的數量必須非常多,才能表示非平行三軸的邊緣。

從 CT 圖中,掃描資料為數張正弦圖 (sinogram),將一個物體放置於轉盤上,轉動轉盤將 X 光照射上去,根據物體的阻射程度得到正弦圖,相當於是光線投影到平面的強度 (CT value)。當一個物體只由一種材質構成時,其阻礙程度與阻礙光線的長度成線性關係,如此一來 CT value 就能有效地反應在物體表面的距離。

這裡補充一點,通常物體不會只由一個材質構成,因此在某些重建工作時會特別弱,如人體掃描重建,則會在眼睛材質上難以得到好的效果。阻礙能力會擴散到其他的正弦圖上,在複合材料的重建工作則會相對地困難。

CT 圖存在些許的干擾,受到鄰近材質的阻射影響。看到網路上文章寫到,即便現代的降噪算法相當厲害,對於醫院普遍使用的成像仍然使用幾十年前的算法,其最重要的原因是擔憂降噪算法捨棄掉的部分到底是不是正確的,去掉資訊若是重要的反而會造成判斷失誤,因此學術界使用 CT 圖的降噪算法仍然不普及於醫學界。

CT

目前 Shepp-Logan filter 相當好用,在 CT 去噪效果受到肯定!公式如下:

$$f(p) = \frac{1}{2} \int_{0}^{2 \pi} \alpha(R_\theta p)^{2} S_{\theta}(P(R_{\theta} p)) d \theta \\ S_{\theta}(q) \leftarrow \int_{-\infty }^{\infty} {\beta(\tilde{q}(x)) S_{\theta}(\tilde{q} (x)) h(x - (q)_{x})} dx \\ h(t) = 2 / (\pi^{2} \delta^{2} - 4 t^{2}) \\$$ $S_{\theta}(q)$ 表示在轉盤角度為$\theta$ 時成像於 q 點的數值$\beta(q) = cos(\angle qso)$,點 q 連線到光源 s、再連線到圓盤圓心點 o。並且$\tilde{q}$ 限定為與點 q 具有相同的 y 值,就是說對水平方向上鄰近點數值去噪$\delta$ 為取樣的像素大小於成像畫面。

從數學公式中看出,受到鄰近點$\tilde{q}$ 影響與光源的偏角差異、與點 q 的距離呈現的關係。

#### 補充 CT 重建 ####

CT 重建方法 中,存在三種常見的方式

迭代法 反投影法 backprojection
濾波反投影法 filtered backprojection

在掃瞄資料中,根據不同角度進行掃描投影,會得到數張 2D 影像,對於實際物體離散取樣,把每一個樣本都當作一個變數 $x_{i}$,根據模擬的方式,可以得到好幾條方程式,解出方程式後,就能知道每個位置上是否有實體 (在物體內部)。但是這樣的效能並不好,會有成千上萬的變數和方程式,甚至有可能發生多組解。


##### 迭代法 #####

迭代法提供一個逼近方法,類似模擬退火法的那種梯度迭代,每一次迭代逐漸將變數數值增加,增加後與實際投影比較,不斷地逼近投影實際值,捨棄掉不合法的數值,每迭代一次讓解更加地接近,當誤差小於某個定值後即可停止。這種作法通常會很慢,但迭代過程中可以加入很多彈性的參數變化,彈性較高。

##### 反投影法 #####

反投影法則是將投影值結果反投影回去,投影的 2D 影像上的某一個點的數值,會反投影到射線上的每一個點,多次反投影後。若該位置存在於物體中,則貢獻於投影像的次數會很多,因此反投影的累加的值也會比較高。接著給定閥值、誤差修正來判斷,效果會更為顯著。反投影法概念上非常簡單,但在邊緣的結果會很模糊。

##### 濾波反投影法 #####

濾波反投影法這部分不是很理解,看起來是先利用 卷積 (Convolution) 操作,再利用傅立葉卷積的反變換,變換過程中需要一個濾波器協助,反變換回來後,再進行反投影法。在傅立葉反變換中需要濾波器,因此選擇濾波器變得相當重要,本篇論文選的是 Shepp-Logan filter 來完成。濾波反投影法的優點是在於輪廓清楚,但會發生 吉布斯現象 (Gibbs phenomenon),逼近原圖時,不連續的函數截面附近會發生震盪。

## 接續 ##

在濾波後,反投影後得到 CT value $f(p)$。當然 $f(p)$ 是一個連續函數的積分,事實上處理時離散的 (因為機器再怎麼精密,齒輪轉角也不會是連續,就看齒數多寡來決定)。

論文中用了 迭代法 濾波反投影法 重建比較,可以得到非常相似的結果。來確認使用的儀器是相當可靠的精準度。

接下來要進行建模,簡單會具有以下幾個步驟。

1. 可以決定要用幾個四面體構成。
2. 三角網格藉由二值化 (CT value 來判定是在物體中還是外部) 抓取,並且加入四面體中。
3. 三角網格到等值曲面的 QEM (Quadric Error Mactrics) 要最小化。
4. 三角網格上的點座標,移動距離購小的時候,算法停止迭代。否則返回第二步,繼續鬆弛新的節點。

第一步進行四面體構成時,通常會採用 * Delaunay Tetrahedralization
(相對於 2D 平面的 Delaunay Trianglution,Delaunay Trianglution 是由 2D Voronoi Diagram,那麼 Delaunay Tetrahedralization 就是由 3D Voronoi Diagram 搭建完成),事實上還有數種建模的方法 (如 Marching cubes 之類的),稍後有機會會提及到。但不管一開始選定的四面體採用哪種策略,經由往後的迭代,將會修整到差不多 QEM。

### 邊緣擷取 ###

對於每一個四面體,找到四面體的重心$g_{i}$ 並且計算 CT value $f(g_{i})$

對於相鄰的四面體,計算 $Tiso= (f(gi)+f(gj))/2$ 接下來要找等值曲面 $Tiso$,找到一個點 p 位於等值曲面上,並且 p 在 gi, gj 的連線上。對於 p 所在的等值曲面上計算法向量 $n_{k}$,等下要用來計算 QEM 用的。

QEM

下篇待續

Read More +