由於在選讀論文時，和上課都遇到相關問題，但一直對這 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，不斷地迭代下去。

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

在幾何課程中介紹的 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 的作法，最後比較他們的優缺點。

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

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 也許能解決這問題。

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

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

單詞概要

• 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

下篇待續