使用全局照明方法有效生成電腦動畫之研究

行政院國家科學委員會專題研究計畫 成果報告

使用全局照明方法有效生成電腦動畫之研究 研究成果報告(精簡版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 99-2218-E-011-005-

執 行 期 間 ： 99 年 02 月 01 日至 100 年 07 月 31 日 執 行 單 位 ： 國立臺灣科技大學資訊工程系

計 畫 主 持 人 ： 賴祐吉

計畫參與人員： 碩士班研究生-兼任助理人員：楊馥聰 碩士班研究生-兼任助理人員：RIDHO 碩士班研究生-兼任助理人員：DOBROMIR 大專生-兼任助理人員：何芳琳

大專生-兼任助理人員：林哲皓 大專生-兼任助理人員：張述淮

報 告 附 件 ： 出席國際會議研究心得報告及發表論文

處 理 方 式 ： 本計畫可公開查詢

中 華 民 國 100 年 07 月 21 日

行政院國家科學委員會補助專題研究計畫 行政院國家科學委員會補助專題研究計畫 行政院國家科學委員會補助專題研究計畫

行政院國家科學委員會補助專題研究計畫 ■ ■成果報告 ■ ■ 成果報告 成果報告 成果報告

□

□ □

□ 期中進度報告 期中進度報告 期中進度報告 期中進度報告

（計畫名稱）

使用全局照明方法有效生成電腦動畫之研究 Efficiently Rendering Animation with Global Illumination

計畫類別：■ 個別型計畫 □ 整合型計畫 計畫編號：

執行期間： 99 年 2 月 1 日至 100 年 7 月 31 日

計畫主持人：賴祐吉 共同主持人：

計畫參與人員：

成果報告類型(依經費核定清單規定繳交)：
精簡報告 □完整報告

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、列管 計畫及下列情形者外，得立即公開查詢

□涉及專利或其他智慧財產權，□一年□二年後可公開查詢

執行單位：台灣科技大學

1 研究背景

電腦圖學已經成為現代電影中不可或缺的技術，特別是用在特效及視覺效果上，而照片擬真的生 成渲染 (photorealistic rendering) 技術更是電影視覺效果的核心，因為不好的生成效果將會使觀 眾失去真實感，而讓電影不賣座。除了電影以外，產品設計師、建築師以及照明設計師也都需要照片 擬真的生成渲染 (photorealistic rendering) 演算法來幫他們預測其產品在各種光線下的外觀，其 中包括建築師需要知道室內設計在不同燈光下的感覺以及建築物在不同天氣狀況下的外觀；舞台設 計師需要瞭解在不同燈光下演唱會舞台的效果為何；汽車設計師想要看到汽車在不同照明環境下會 有怎樣的外觀等。這些都說明了基於物理的生成渲染技術的主要優點在於：即使沒有建立實際的場 景，也能產生影像去準確的預測在真實世界的外觀。但是，傳統基於 CPU 的全局照明 (CPU-based global illumination) 演算法需要很大的計算量，其結果是不可能即時生成複雜的場景，亦即不 可能讓使用全局照明互動地編輯動畫場景。雖然，目前也已經有基於 GPU 的全局照明 (GPU-based global illumination) 演算法可以使用，即使 GPU 平行運算核數大增，但是仍無法達到全局照明的 即時編輯效果。另外，在互動圖學及電腦遊戲兩大領域中，基於 GPU 的局部照明 (GPU-based local illumination) 演算法提供使用者可以與即時生成結果互動的能力，此類方法的主要目的只是在得到 近似真實世界以及偽造出真實感覺的效果，但是在設計之初並沒有考慮到所有可能的光影效果以及物 體間的相互關係，而造成生成的影像會缺乏部分真實的燈光效果。雖然隨著電腦及演算法的改進，基 於 GPU 的局部照明 (GPU-based local illumination) 演算法、基於 CPU 的全局照明 (CPU-based global illumination) 演算法以及基於 GPU 的全局照明 (GPU-based global illumination) 演算 法都不斷的在進步，但是三種都是獨立設計執行，無法綜合使用。因此，計畫中整合這三個系統來建 構一個綜合系統，該系統會使用局部照明 (local illumination) 演算法來提供場景轉換時所需的中 間訊息，局部照明演算法可以提供編輯動態場景 (animated scenes) 時所需的互動性及場景最終擬 真生成結果的預測性，用以節省使用全局生成法時，設定參數及場景所需的試驗時間，讓他們可以用 更簡單、更方便的設定所需的場景；同時基於 CPU 的全局照明 (CPU-based global illumination) 演算法以及基於 GPU 的全局照明 (GPU-based global illumination) 演算法則提供了最後物理真實 的模擬結果，讓設計者可以得到真實的感受。在計畫中，使用電腦遊戲引擎 OGRE3D 所提供的大量的 技術來做近似即時渲染的效果，這讓我們能夠有彈性的編輯欲得到的場景，並用更快的版本來檢查結 果，這提供了基於 CPU 以及基於 GPU 的生成渲染演算法所無法達到的互動性。而後端的真實生成技 術則將使用基於 CPU 的生成渲染函式庫－ PBRT，以及基於 GPU 的生成渲染函式庫－ Optix 來生 成擬真動畫。最終，此專案將產生一個整合基於 GPU 的局部照明 (GPU-based local illumination) 演算法、基於 CPU 的全局照明 (CPU-based global illumination) 演算法以及基於 GPU 的全局照 明 (GPU-based global illumination) 演算法來編輯動畫場景、預知動畫生成結果以及用以生成擬 真照片的動畫系統。除了建構一個整合系統外，此計畫亦研究如何提生動畫生成的效率，產生更佳的 擬真生成效果演算法。而擬真演算法在未來互動 3D 電視中可能成為一個重要技術，計畫中，我們對 3D 電視技術做一 SURVEY，來研究生成技術在未來 3D 電視上的可能性。

文獻回顧與探討

在互動圖學及電腦遊戲的領域中，大部份互動式的編輯工具及遊戲都是使用局部照明生成渲染 (local illumination) 演算法，特別是以 OpenGL 和 DirectX 為生成渲染的基本工具，雖然電腦硬 體大幅進步，運用硬體來生成的相關研究也大量被發表，但是這類研究大部分強調的是互動性而非正 確性，因此其生成的影像只能模擬部分燈光效果而非真實全部效果，以下會回顧部分的局部照明相關 研究。首先是反射，反射在全局照明 (global illumination) 是很重要的效果，Diefenbach [9, 10]

以及 Wojdala 等人 [39] 提出使用 GPU 的多通演算法來生成渲染平面鏡及曲面鏡的反射結果；多位 學者提出了使用環境映射來處理曲面反射的生成演算影像 [21, 23, 24, 6]；另有學者提出在空間 中，在所有可能表面上取樣點，在每個樣本點上，根據所有可能的視線方向以及光的方向，計算所有 可能的燈光效果並儲存。接著在生成渲染的階段，預先計算的資料會在動態、局部、以及低頻的光線 下，用來快速生成渲染物體 [28, 2, 19, 30, 3, 1, 5]；學者們也提出利用紋理 (texture) 技術，

來計算照明以及產生陰影 [17, 34, 35, 20, 22]。在點陣化硬體 (rasterization hardware) 上，

通常考慮到場景的漫射用以加速互動式顯示，有學者 [40, 41, 42, 12] 提出使用基於影像的散焦映 射 (image-based caustics map) 演算來加入玻璃以及水等等的聚光效果；這些基於 GPU 的局部生

結案報告 共 10 頁第 1 頁

成演算法必須省略部分燈光效果，用以加速計算，達到互動效果，因此逼近的生成結果並無法模擬真 實照片的所有燈光效果。

基於 CPU 的生成渲染演算法大部分採用逐幀的方式來生成渲染一系列的動畫 [31]；由於連續 畫格間在時序有很大的關聯性，因此利用此關聯性在改善全域演算法的效率是極重要的工具。實際 上，若不考慮時序上的相關性，整部動畫將可能產生一些引人注意的影像間的人工雜訊，像是閃爍 或是閃耀。過去的研究中顯示，在所有的畫格中反覆使用緩衝光子 (caching photons) 和輻射樣本 (irradiance samples) [38, 13, 27, 15, 36, 26, 7] 能夠大幅改善計算效率並減少時間上的失 真。但是，這類動畫生成演算法基本上會遭遇的困難是緩衝樣本 (cached samples) 會使畫格的最後 結果產生偏差及錯誤。若是動畫的場景中，只有攝影機及燈光來源允許移動，藉由重複評估路徑的 貢獻 [4, 18, 33]，或是基於多重重要性架構 (multiple-importance) 重新加權樣本，有數個基本 的 Monte Carlo (MC) 技術可以用來重複使用燈光路徑。但也因為限制了只有攝影機和燈光來源的位 置可變，而減少了這些技術的實用性。Havran 等人 [18] 則提出從雙向路徑追蹤 (bi-directional path tracing) 重複使用靜態物件的路徑以減少時間上的閃爍，對於每個靜止的候選路徑，首先要先 檢查路徑的正確性 (也就是說，所有路徑不可被物件的任何一部份擋住)，然後重複計算攝影機第一 個照射點的 BRDF 值，用以評估每一個正確路徑的貢獻，但重複加權的方案卻會造成偏差及不良的雜 訊，而且也只能重複使用靜止的路徑，其他的路徑就無法使用了。Ghosh 等人 [16] 利用路徑間時間 的關聯性，提出了連續 Monte Carlo 的架構，這個方法是從上一個畫格中計算目前畫格的樣本，來重 新加權樣本，不過該方法只限制在環境映射的燈光，而無法追縱所有的路徑。另外，還有少部份的研 究是使用圖形硬體架構來做全局照明生成渲染 (global illumination) 的演算 [11, 32]，這類論文 主要的研究重點在於如何利用 GPU 的平行運算功能來加速全局照明生成渲染 (global illumination) 的運算，而 Owens 等 [29] 提供了使用 GPU 平行運算的縱覽。這些全局照明的演算法雖然可以模擬 照片中所有燈光效果，但是其所需的計算量十分巨大，而造成在設計好場景後，要等很久才能得到結 果，每次修改場景或參數，皆需很長的時間才能看到變化，無法即時編輯，因此，通常其輸入是一個 檔案，而非可以動態編輯的程式系統。

因此，本計畫想使用遊戲引擎的互動功能設計出一個可以創造、編輯及預測全局動畫場景的系統，

同時系統會整合基於 CPU 或是基於 GPU 的全局生成演算系統，讓使用者可以決定如何生成最後的擬 真動畫，此系統允許使用者在同一平台上發展三類演算法，甚至混合三類演算法建構一個更加有效率 的生成法。

2 研究成果

本計畫主要著重在有效率地生成動畫畫面及達成即時互動生成的方法，因此除了人員訓練外，我 們分別完成四項工作:

2.1 GPU 及 CPU 的混合系統

Fig 2呈現本計畫所欲建構的可能系統圖及所需要的子系統模組，該系統將以 OGRE3D 為主要管理 者的核心子系統，用以連結各個子系統，生成引擎 OGRE3D 可以提供簡單且方便的使用者界面、擁有 很棒的互動性及適當場景管理資料結構，同時該系統已經被用在很多應用中，具有適當多數的文件，

可以讓開發的過程更順利。在系統中，藍色模組區塊表示是由 OGRE3D 所建構出來的互動子系統。在 計畫中，我們設計建構以下幾個模組：

1. 場景圖模組 (scene graph) 是整個系統的管理核心，其功能負責處理安排場景的結構，在 OGRE3D 中，場景圖 (graph) 或樹 (tree) 來管理，每個模型是形成一個節點群，每個節點可 以有很多子節點，但通常都只有一個父節點；而在父節點上的效果也會套用至所有的子節點，

因此在一個群組上做運算會自動把效果傳播到裡面所有的成員。在場景圖中的每個節點或物件 的特性將會利用模型模組 (models block) 來編輯，此模組會加入時間編輯及控制動畫的時間，

而動畫中的攝影機移動或是模型物體的性質變換包括亮度、材料及位置方向等，則是利用動畫 關鍵圖模組 (animation key block) 來編輯。

2. 模型模組 (models) 用來輸入不同格式旳模型，並用以編輯處理場景中的各個模型物體，用以

Figure 1: 此圖顯示本計畫所欲建構的系統圖。在圖中，大致可以由三個不同系統所構成，由 OGRE3D 所構成的互動編輯子系統、由 PBRT 所構成的 CPU 全局生成渲染演算子系統以及由 OPTIX 所建構的 GPU 全局生成渲染演算子系統。

結案報告 共 10 頁第 3 頁

3. 動畫關鍵圖模組 (animation key) 用來定義物體的移動及改變，包含了會用來控制移動及改變 的關鍵影格。

4. 局部生成渲染模組負責使用 OGRE3D 的局部生成渲染器 (local renderer) 來顯示或是預覽在 OGRE3D 中的場景，雖然其生成在物理上可能不正確，但是可以提供即時回饋的功能。

全局照明 (global illumination) 是 3D 電腦圖學的一種演算法，是為了在 3D 場景中增加更真實 的燈光。這樣的演算法不只考慮到從光源來的光，也考慮到同樣光源的光線隨後會被場景中其他的表 面反射，不論是反光的或非反光的。基於 CPU 的全局照明 (CPU-based global illumination) 生成 渲染程序會用 PBRT 來處理。紅色模組區塊則是由 PBRT 所建構的 CPU 全局照明生成渲染演算子系 統，可以使用一個或多個 CPU 來生成渲染。OGRE3D 會傳送此區塊所需的即時可用的場景圖形資訊，

來完成在 CPU 上生成渲染的程序。基於 GPU 的全局照明 (GPU-based global illumination) 生成 渲染程序會用 OptiX 來處理。綠色模組區塊是由 OptiX 所建構的 GPU 全局照明生成渲染演算子系 統，生成渲染程序需要 NVIDIA GPU 所提供的 CUDA 來運行。基於速度的理由，此程序需要使用 GPU 內的 CUDA 來平行處理，OGRE3D 會傳送此區塊所需的即時可用的場景圖形資訊，來完成在 GPU 上生 成渲染的程序。最後黃色模組區塊則是由本實驗室所開發出的場景儲存及其它生成功能的子系統，其 中比較重要的是儲存模組 (storage) 用以負責讀取、分析及儲存場景檔案，以便日後重覆使用，方 便研究。

該程式使用以下的步驟來整合此三種系統：

1. 以 OGRE3D 為基礎來建構場景編輯系統，計畫將設計適當的使用者介面及互動模式，預測生成 及場景管理等子系統。再加入模型 (model) 子系統讓系統可以讀取 PBRT 場景及其它動畫軟體 模型檔案、刪除場景中的模型以及為指定模型挑選紋理 (texture) 或材質等功能，另外，設計 場景檔案輸出功能，用以儲存編輯後的場景，做為未來研究使用。

2. 結合 PBRT 的單核心或多核心生成渲染演算法，讓系統可以使用 CPU 來生成最後演算結果。學 生欲讓 PBRT 輸入檔案格式當做整個系統主要的輸出及輸入格式，因此系統中需要整合 PBRT 的函數用以讀取場景檔案、分析場景檔案內容。另外，讓使用者可以選擇 PBRT 中不同生成演 算法並且設定適當的參數用以生成渲染動畫結果，其結果將再使用 PBRT 的格式輸出，並且使 用 HDR 方法輸出可見的結果。

3. 結合 NVIDIA 的 OptiX 之 GPU 生成渲染演算法，讓系統可以使用 GPU 來生成最終結果。學生 將要讓 OptiX 系統合併進入系統中，讓使用者可以選取 GPU 的生成演算法並且設定適當的參 數，用以生成結果，同時因為 OptiX 只有基本生成功能，為了要配合 PBRT，所以學生將書寫 配合的相對應輸入場景敍述方法及材質等生成相關的對應功能函數，同時亦書寫適當的傳輸函 數將資料傳到 GPU 記憶體中，最後生成的結果亦以 PBRT 的格式從 GPU 中輸出，以上述的 HDR 方法輸出可見結果。

該系統將被用於本實驗室的未來研究中。

2.2 Population Monte Carlo Image Plane Sampler [25]

本演算法是從統計學中的 population Monte Carlo (PMC) 取樣架構 [14] 演進出來的，PMC 演 算法則是依據在迭代過程中所獲得的資訊自動調整取樣的分配函數，並且有理論上的低錯誤保證 並且具低的演算成本。Population Monte Carlo image plane sampler(PMC-IP) 則是以影像生成面 (image plane) 上的取樣位置當做計算群組。來預測產生一個良好的啟動要點函數 (kernel)，本演 算法乃採用電腦視覺的技術，用以探索時間上相鄰影像間的取樣分配函數的相關性。這個預測要點函 數 (kernel) 可以避免浪費一些多餘的取樣子在於探測在影像中光樣平順圓滑的區域。該法是利用一 張目前使用少量均量分佈取樣子所生成的中間結果圖、前一張影像生成成果圖以及前一張影片的取樣 分配密度圖為基礎去預測初始要點函數 (kernel)。另外，除了預估要點函數 (kernel) 外，生成每 張影像的迭代過程中，使用中間生成影像計算每個像素的視覺感受加權變異值，再根據每一個像素的 加權變異值最後用於調整要點函數 (kernel)，目的是將較大量的影像面 (image plane) 取樣子放置 在視覺高變異值的區域，用以減少該區域的高變異值及其所造成的視覺不舒服。本演算法可以達成將 變異值較平均地分配到影像面上，以增加影像生成的品質。Fig. 3跟 4呈現了使用 PMC-IP 來生成動

Figure 2: 左圖乃在 OGRE 系統中，所可以互動控制場景所生成的影像，而在中圖則是相同場景用 PBRT 所生成的影像成果，而在右圖中則是使用 OPTIX 所生成的影像成果。

2.3 3D 電視技術 [37]

自從阿凡達 3D 立體電影上市風迷世界以來，3D 立體顯示己成近來媒體的寵兒，更成為在 CSE2010 展覽會中最引人注目的科技。幾乎所有的電視大廠都己經或將在 2010 中，推出他們的 3D 電視產品，屆時 3D 內容也可能隨著上市，其可能內容包括 2010 下半有多款藍光 3D 電影即將出現 在市場上，外加全球多家著名電視台亦將推出 3D 立體頻道。當然這波熱潮也同時掃向了個人電腦、

個人手機及電腦遊戲市場。但是目前市面上，最新款式的 3D 立體視覺系統即所謂的二眼式立體視覺 系統，仍在初期發展階段，這類系統通常要求使用者配戴特別的眼鐘或是頭載式設備觀看固定角度之 影片。雖然市場上預期下一代的 3D 立體感受系統將於 3 到 5 年內上市，新一代的系統不需要配 戴特別裝備，且可偵測使用者動作用於播放內容互動。但是目前技術仍需克服許多的挑戰包括了光學 系統的設計、顯示的技術、使用者互動介面、3D 訊號的傳輸、3D 資料的表示及編碼以及 3D 的通訊 系統等。此論文 [37] 呈現了 3D 影像通訊系統中主要技術的縱覽概述，特別是 3D 影片數位內容從 製作端到透過通訊系統傳輸到使用者端的生命週期為何？一般使用者對 3D 數位內容的感受為何？有 那些因素是可以影響使用者觀看的感受？如何表現及壓縮此類 3D 數位內容或虛擬實境？如何將資料 有效率地透過網路傳輸到使用者端？從生成者端、網路傳輸端及接收端的觀點上，如何設計新的技術 以提昇 3D 立體數位內容體驗的品質？在同時本論文也另外強調數個研究的可能方向給讀者，包括了 如何調整未來的網路架構才得以满足 3D 影片創新數位內容的傳輸服務要求？在 3D 影片傳輸系統 中，如何針對有效隱藏傳輸誤差以滿足使用者對擬真的 3D 影片內容之期待，而設計適當的 Quality of Experience (QoE)？如何有效地藉由 3D 數位內容呈現及編碼中不同重要性及內容重覆性去使用 unequal error protection (UEP) 以及動態資源配置來改良系統的表現？何種編碼系統可以被使用 來增加編碼後的訊息流產生錯誤的抵抗力？如何設計在 4G 網路上傳輸 3D 影片數位內容的時程規畫 者及資料配置方式？如何有效地使用 4G 網路中 QoS 來支授 3D 影像內容傳輸？如何利用新的計算 概念像是綠色運算、雲端運算以及分散／合作運算等來設計傳輸源／傳輸者／接受者的點對點 3D 影 片傳輸系統？如何應用易擴縮性的概念來處理 3D 立體影片在異質混合網路通訊系統上的傳輸？

2.4 Voxelization [8]

在做實體的建模編輯跟生成前，一般需要演算方法將物體的表面資料 (surface representation) 轉換成以實體資訊 (volume representation)，此一類的演算法即是所謂的 Voxelization。由於 GPU(Graphic Processing Unit) 的快速進步，目前使用 GPU 為基礎的 voxelization 演算法，幾乎 己經可達成即時的運算得到結果。在論文 [8] 中提出一種新型以 GPU 演算為基礎的 voxelization 演算法。本演算法是利用設定 clipping plane 的方式將表面資料切片，再從中算取實體資料的演算 方法 [?]，在原方法中需要很多次 GPU 的轉移操作或是使用特別的編碼方式，因而限制其有效演算 速度及變換的複雜度。論文提出的演算法則是運用 GPU 新增之硬體功能包括 geometry shader 跟可

結案報告 共 10 頁第 5 頁

Figure 3: 這乃是一系列使用 PMC-IP 演算法所生成之 Cornell Box 動畫中的影像。最上一排是使 用 PMC-IP 法用 16 次迭代後所生成結果動畫中編號 1、31、61 和 91 的影像，第二排則是在演算 完成後的 α_(i,j) 值，第三排則是使用在電腦視覺預估法 [25] 所估算出的 α_(i,j) 值。該論文分析四個 影像的估計 α(i,j) 累積根號誤差分別為 9.28× 10⁻⁷, 1.26× 10⁻⁶, 8.81× 10⁻⁷, and 1.16× 10⁻⁶.

將 3D texture 設定成生成目標等來克服先前 slicing 演算法的限制。本演算法是利用 geometry shader 將繪製的基本圖形根據所可能交會的切片數去複製，再利用複製的三角形去執行 slice 演 算法，來獲取 voxelization 之結果，該法同時記入複雜表面材料資訊，且將資訊存放在可調整 3D texture 中。因為使用了 geometry shader 克服了原本需要多次轉移操作生成的限制，使用了 3D texture 讓應用程式可以根據其需要及硬體限制而自由調整 voxel 的解析度，並且在調整過程中，

無需要負擔如改變 grid-encoding 演算法中複雜編碼手續及結構的負荷。3D texture 的使用亦允許 將大量的附加資料直接存入 voxel 群體中，可存入的資料包括物體表面的 normal 跟 transmission coefficient 等，除此之外，3D texutre 乃是在 GPU 記憶體空間中，讓本法可以將 voxelization 的結果直接應用於更多不同的生成計算中，其所附加的表面細節資料更可用於產生複雜且豔麗的影 像結果。該論文提供了三個不同的應用來驗證本演法之可行性及提昇演算速度。Fig. 5展示了使用 transparency 跟 refraction 的計算結果。

本結果將於今年九月發表於高雄所舉行的 Pacific Graphics。

參考文獻

[1] Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments, New York, NY, USA, 2002. ACM.

[2] Clustered principal components for precomputed radiance transfer, New York, NY, USA, 2003. ACM.

[3] Precomputed local radiance transfer for real-time lighting design, volume 24, New York, NY, USA, 2005. ACM.

Figure 4: 這乃是使用 PMC-IP 演算法所生成的房間動畫。第一排是使用了 PMC-IP 演算法迭代演 算 64 次後所生成的動畫中編號 1、31、61 和 91 的影像。第二排則是使用每個像素 1024 個抽樣 子的傳統 path tracing 演算法所生成的動畫中編號 1、31、61 和 91 的影像。觀察生成的動畫，

PMC-IP 所生成的結果可以大量減少傳統 path tracing 中突出現雜訊，增加結果的可看度。

[4] P. Bekaert, M. Sbert, and J. Halton. Accelerating path tracing by re-using paths.

In Rendering Techniques ’02 (Proceedings of the Thirteenth Eurographics Workshop on Rendering), pages 125–134, 2002.

[5] A. Ben-Artzi, R. Overbeck, and R. Ramamoorthi. Real-time brdf editing in complex lighting. ACM Trans. Graph., 25(3):945–954, 2006.

[6] B. Cabral, M. Olano, and P. Nemec. Reflection space image based rendering. In In SIGGRAPH 99, pages 165–170, 1999.

[7] M. Cammarano and H. W. Jensen. Time dependent photon mapping. In In Proceedings of the 13th Eurographics Workshop on Rendering, pages 135–144, 2002.

[8] S. Chang, Y. Lai, Y. Niu, F. Liu, and K. Hua. Real-time realistic voxel-based rendering. In Pacific Graphics 2011, 2011.

[9] C. Dachsbacher and M. Stamminger. Reflective shadow maps. In I3D ’05: Proceedings of the 2005 symposium on Interactive 3D graphics and games, pages 203–231, New York, NY, USA, 2005. ACM.

[10] C. Dachsbacher and M. Stamminger. Splatting indirect illumination. In I3D ’06:

Proceedings of the 2006 symposium on Interactive 3D graphics and games, pages 93–100, New York, NY, USA, 2006. ACM.

[11] C. Damez, K. Dmitriev, and K. Myszkowski. State of the art in global illumination for interactive applications and high-quality animations. Computer Graphics Forum, 22(1):

55–77, march 2003.

[12] S. T. Davis and C. Wyman. Interactive refractions with total internal reflection. In Proceedings of Graphics Interface 2007, GI ’07, pages 185–190, New York, NY, USA, 2007.

ACM.

[13] K. Dmitriev, S. Brabec, K. Myszkowski, and H.-P. Seidel. Interactive global illumi- nation using selective photon tracing. In Rendering Techniques ’02 (Proceedings of the 13th Eurographics Workshop on Rendering), pages 25–36, 2002.

結案報告 共 10 頁第 7 頁

Figure 5: 左圖是一個透明的笑面佛像在擁有火爐的房間環境光源下，使用 voxelization 結果的生 成圖，右圖則是一隻龍的模型在房間環境光源上，使用 voxelization 結果的生成圖。在二圖中我們 將反射及折射的材質顏色設為不同，產生一個視覺上的特殊效果，同時也顯示本法計算之準確性。

[14] R. Douc, A. Guillin, J. M. Marin, and C. P. Robert. Minimum variance importance sampling via population Monte Carlo. Technical report, University Paris Dauphine, 2005.

[15] P. Gautron, K. Bouatouch, and S. Pattanaik. Temporal radiance caching. Technical Report 1796, IRISA, Rennes, France, 2006.

[16] A. Ghosh, A. Doucet, and W. Heidrich. Sequential sampling for dynamic environment map illumination. In Proc. Eurographics Symposium on Rendering, pages 115–126, 2006.

[17] P. Haeberli and K. Akeley. The accumulation buffer: hardware support for high-quality rendering. In SIGGRAPH ’90: Proceedings of the 17th annual conference on Computer graphics and interactive techniques, pages 309–318, New York, NY, USA, 1990. ACM.

[18] V. Havran, C. Damez, K. Myszkowski, and H.-P. Seidel. An efficient spatio-temporal architecture for animation rendering. In EGRW ’03: Proceedings of the 14th Eurographics workshop on Rendering, pages 106–117. Eurographics Association, 2003.

[19] M. Haˇsan, F. Pellacini, and K. Bala. Direct-to-indirect transfer for cinematic relighting. ACM Trans. Graph., 25(3):1089–1097, 2006.

[20] T. Heidmann. Real Shadows - Real Time. Iris Universe, 1991.

[21] W. Heidrich. Interactive display of global illumination solutions for non-diffuse environments - a survey. Computer Graphics Forum, 20:1067–7055, 2000.

[22] W. Heidrich and H. peter Seidel. Realistic, hardware-accelerated shading and lighting.

In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pages 171–178, 1999.

[23] J. Kautz and M. D. Mccool. Interactive rendering with arbitrary brdfs using separable approximations. In In Eurographics Rendering Workshop, pages 281–292, 1999.

[24] J. Kautz, P. pau Vazquez, W. Heidrich, H. peter Seidel, and A. Udg. A unified approach to prefiltered environment maps. pages 185–196. Springer, 2000.

[25] Y. Lai, S. Chenney, F. Liu, and S. Niu, Y.and Fan. Animation rendering with population monte carlo image-plane sampler. The Visual Computer, pages 543–553, 2010.

[26] W. Martin, W. Martin, S. Parker, S. Parker, E. Reinhard, E. Reinhard, P. Shirley, P. Shirley, W. Thompson, and W. Thompson. Temporally coherent interactive ray tracing.

Journal of Graphics Tools, 2:41–48, 2002.

[27] M. Meyer and J. Anderson. Statistical acceleration for animated global illumination.

ACM Trans. Graph., 25:1075–1080, July 2006.

[28] R. Ng, R. Ramamoorthi, and P. Hanrahan. All-frequency shadows using non-linear wavelet lighting approximation. ACM Trans. Graph., 22(3):376–381, 2003.

[29] J. D. Owens, D. Luebke, N. Govindaraju, M. Harris, J. Kr � ger, A. E. Lefohn, and T. J. Purcell. A survey of general-purpose computation on graphics hardware. Computer Graphics Forum, 26(1):80–113, 2007.

[30] F. Pellacini, K. Vidimˇce, A. Lefohn, A. Mohr, M. Leone, and J. Warren. Lpics: a hybrid hardware-accelerated relighting engine for computer cinematography. In SIGGRAPH

’05: ACM SIGGRAPH 2005 Papers, pages 464–470, New York, NY, USA, 2005. ACM.

[31] M. Pharr and G. Humphreys. Physically Based Rendering from Theory to Implementation. Morgan Kaufmann, 2004.

[32] T. J. Purcell, C. Donner, M. Cammarano, H. W. Jensen, and P. Hanrahan. Photon mapping on programmable graphics hardware. In Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on Graphics Hardware, pages 41–50. Eurographics Association, 2003.

[33] M. Sbert and F. Castro. Reuse of paths in final gathering step with moving light sources. In International Conference on Computational Science 2004, pages 189–196, 2004.

[34] M. Segal, C. Korobkin, R. van Widenfelt, J. Foran, and P. Haeberli. Fast shadows and lighting effects using texture mapping. SIGGRAPH Comput. Graph., 26(2):249–252, 1992.

[35] J. Shade, D. Lischinski, D. H. Salesin, T. DeRose, and J. Snyder. Hierarchical image caching for accelerated walkthroughs of complex environments. In SIGGRAPH

’96: Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pages 75–82, New York, NY, USA, 1996. ACM.

[36] M. Smky, S.-i. Kinuwaki, R. Durikovic, and K. Myszkowski. Temporally coherent irradiance caching for high quality animation rendering. In The European Association for Computer Graphics 26th Annual Conference EUROGRAPHICS 2005, volume 24 of Computer Graphics Forum, pages xx–xx, Dublin, Ireland, 2005. Blackwell.

[37] G.-M. Su, Y.-C. Lai, A. Kwasinski, and H. Wang. 3d video communications: Challenges and opportunities. International Journal of Communication Systems, pages n/a–n/a, 2010.

[38] T. Tawara, K. Myszkowski, K. Dmitriev, V. Havran, C. Damez, and H.-P. Seidel.

Exploiting temporal coherence in global illumination. In SCCG ’04: Proceedings of the 20th spring conference on Computer graphics, pages 23–33. ACM Press, 2004. Invited Talk.

結案報告 共 10 頁第 9 頁

[39] A. Wojdala, M. Gruszewski, and K. Dudkiewicz. Using harware texture mapping for efficient image synthesis and walkthrough with specular effects. In Machine Graphics and Vision, pages 137–151, 1994.

[40] C. Wyman. An approximate image-space approach for interactive refraction. ACM Trans.

Graph., 24:1050–1053, July 2005.

[41] C. Wyman. Hierarchical caustic maps. In Proceedings of the 2008 symposium on Interactive 3D graphics and games, I3D ’08, pages 163–171, New York, NY, USA, 2008. ACM.

[42] C. Wyman and S. Davis. Interactive image-space techniques for approximating caustics.

In Proceedings of the 2006 symposium on Interactive 3D graphics and games, I3D ’06, pages 153–

160, New York, NY, USA, 2006. ACM.

Noname manuscript No.

(will be inserted by the editor)

Animation Rendering With Population Monte Carlo Image-Plane Sampler

Yu-Chi Lai · Stephen Chenney · Feng Liu · Yuzhen Niu · ShaoHua Fan

Received: date / Accepted: date

Abstract Except for the first frame, population Monte Carlo image plane (PMC-IP) sampler renders with a start-up ker- nel function learned from previous results by using motion analysis techniques in the vision community to explore the temporal coherence existing among kernel functions. Pre- dicted kernel function can shift part of the uniformly dis- tributed samples from regions with low visual variance to regions with high visual variance at the start-up iteration and reduce the temporal noise by considering the temporal rela- tion of sample distributions among frames. In the following iterations, the PMC-IP sampler adapts the kernel function to select pixels for refinement according to a perceptually- weighted variance criterion. Our results improve rendering efficiency by a factor between 2 to 5 over existing techniques in single frame rendering. The rendered animations are per- ceptually more pleasant.

Keywords Ray-tracing · PMC · Monte Carlo · Global Illumination

1 Introduction

Global illumination based on Monte Carlo integration pro- vides the most general solution for photorealistic render- Yu-Chi Lai

National Taiwan University of Science and Technology, R.O.C. E- mail: [email protected]

Stephen Chenney

Emergent Game Technology, U.S.A. E-mail: [email protected] Feng Liu

University of Wisconsin-Madison, U.S.A. E-mail: [email protected] Yuzhen Niu

Shandong University, China E-mail: [email protected] ShaoHua fan

Financial Engineering Research Center at Suzhou University, China and NYU, U.S.A E-mail: [email protected]

ing problems. To reduce image noise (variance) at practi- cal computation times is the core of research in the global illumination community. To render an animation the tempo- ral variance at each pixel among consecutive frames must also be considered because our eyes are good at noticing temporal inconsistency among consecutive frames for sur- vival purpose; our algorithm is derived from the Population Monte Carlo (PMC) sampling framework, which is a tech- nique that adapts sampling distributions over iterations, all with theoretical guarantees on error and little computational overhead. Our algorithm generates a start-up kernel function from previous rendered frames by considering the temporal correlation among kernel functions. In addition, our algo- rithm evenly distributes the variance over the image plane in each frame to remove noisy spikes on the image and in turn, reduce the temporal inconsistency generally existing in a frame-by-frame ray-tracing based algorithm.

PMC algorithms iterate on a population of samples. For our sampler, image-plane sampling (PMC-IP), the popula- tion is a set of image-plane locations. To render each frame, the population is initialized with a sample distribution whose predictions are based on the previous rendering results, and then PMC-IP generates an intermediate image to start the it- eration. Any information available at this stage can then be used to adapt a kernel function that produces a new popu- lation. The initial prediction of the kernel function is based on the result of the current frame rendered with a few sam- ples in each pixel, the result of previous frames, and the sample distribution of previous frames. We can explore the temporal coherence of sample distributions among consec- utive frames by using computer vision techniques to gener- ate a good start-up kernel function. This prediction prevents the redundant probed samples on smooth regions of the im- age. In addition, the prediction takes the temporal variance into account: perceptually high variance regions in previ- ous frames have higher probability to be perceptually high

2

variance regions in current frame, and as a result, the algo- rithm puts more samples at these regions to reduce the noisy spikes. In image-plane sampling, the perceptually-weighted variance in the intermediate images is used to construct the kernel function, resulting in more image plane samples in re- gions of high variance. The procedure is then iterated: sam- ple, adapt, sample,. . . . The result is an unbiased adaptive algorithm. This can achieve an evenly-distributed variance over the image plane.

Photon mapping and path tracing have been the indus- trial rendering algorithms for global illumination. Our sam- pler can be easily incorporated into ray-tracing based global illumination with minimal modifications to improve the ef- ficiency of these algorithms. We demonstrate the ease of in- corporation into the current rendering framework by modi- fying the ray-starting pixel positions.

Our contribution is a specific tool for rendering that uses the Population Monte Carlo framework: An image-plane sampler, PMC-IP, that adapts to guide samples to percep- tually high variance image regions, is cheap to compute, maintains stratification, and is unbiased. In addition, we in- corporate motion analysis techniques from the vision com- munity into global illumination algorithms to explore the temporal coherence of sample distributions among consec- utive frames to improve the samples usage and enhance the temporal consistency among consecutive frames.

We include results comparing each algorithm to existing approaches. We find that PMC-based algorithms improve the efficiency in a factor of 2 to 5 over existing methods when rendering a single frame. In addition, our algorithm can generate a more perceptually pleasant animation than others.

2 Related Work

Here we focus on the adaptive image-plane sampling and algorithms using sequential Monte Carlo algorithm. For an overview of Monte Carlo rendering in general, see Pharr and Humphreys [21].

Typically, adaptive image-plane algorithms perform a first pass with a small number of samples per pixel and use the resulting values to label pixels as adequately sampled or in need of further refinement. The algorithm then iterates on the pixels requiring more samples [10, 20, 22, 19, 2, 23].

A common property of these existing algorithms is that they stop sampling a given pixel when some image-derived metric is satisfied. As Kirk and Arvo [12] point out, the evaluation of the image metric relies on random samples, so there is some non-zero probability that the threshold is incorrectly detected and that sampling stops too soon. This introduces bias in the final image, which is a problem when physically accurate renderings are required. Our algorithm

never uses a threshold to stop sampling a pixel and is statis- tically unbiased.

Many metrics have been proposed for the test to trig- ger additional sampling. Lee et al. [15] used a sample vari- ance based metric. Dipp´e and Wold [4] estimated the change in error as sample counts increase. Painter and Sloan [20]

and Purgathofer [22] used a confidence interval test, which Tamstorf and Jensen [29] extended to account for the tone operator. Mitchell [19] proposed a contrast based criteria because humans are more sensitive to contrast than to ab- solute brightness, and Schlick [26] included stratification into an algorithm that used contrast as its metric. Bolin and Meyer [2], Ramasubramanian et al. [23] and Farrugia and P´eroche [6] used models for human visual perception, of which we use a variant.Rigau et al. [24, 25] introduced entropy- based metrics. [18? ] introduced metrics utilizing a model of the human visual system from the visible distortion based on the detection and classification of visible changes in the image structures.

Our algorithm views the image plane as a single sample space for the purposes of sampling. Dayal et al. [3] used a variance-based metric to control a kD-tree subdivision where samples are drawn uniformly within each adaptively sized cell of the subdivision. Stokes et al. [27] also took a global approach with their perceptual metric.

A Sequential Monte Carlo algorithm, similar in spirit to Population Monte Carlo, has recently been applied by Ghosh, Doucet and Heidrich [9] to the problem of sampling environment maps in animated sequences. Their work ex- ploits another property of iterated importance sampling al- gorithms – the ability to re-use samples from one iteration to the next – and is complementary to our approach.

Lai et al [13, 14] adapted the Population Monte Carlo algorithm into the energy redistribution frame-work to adapt the extent of energy redistribution. Their method focuses on the global rendering algorithm but our algorithm focuses on the first step of sample distributions. Our algorithm can be easily adapted to incorporate the algorithm by adjusting the energy of each path according to the kernel function of pixel position distribution.

3 PMC-IP: Image-Plane Sampling

To render a frame in a sequence of animation is to compute the intensity, I(i, j,t), of each pixel (i, j) for each frame t, by estimating the integrals:

I_{i, j,t}= Z

IW_{i, j,t}(u)L(x, ω,t)du, (1)

whereI is the image plane, Wi, j,t(u) is the measure- ment function for pixel (i, j) of t-th frame – non-zero if u is within the support of the reconstruction filter at (i, j) of t- th frame – and L(x, ω,t) is the radiance leaving the point, x,

3

seen through u in the direction −ω at t-th frame, determined by the projection function of the camera. We are ignoring, for discussion purposes, depth of field effects, which would necessitate integration over directions out of the pixel, and motion blur, which would require integration over time. In the following discussion, we will neglect the denotation of t for the simplification of description. And all the adaptation of the sample distribution happens in the process of render- ing a single frame.

An image-plane sampler selects the image-plane loca- tions, x in Equation 1 for a specific frame. For simplicity, assume we are working with a ray-tracing style algorithm that shoots from the eye out into the scene. Adaptive sam- pling aims to send more rays through image locations that have high noise, while avoiding bias in the final result.

Taking an importance sampling view, given a set of sam- ples, {X₁, . . . , X_N} from an importance function p(x) for a single frame, each pixel is estimated using

Iˆi, j=1 n

N

∑

k=1

W_{i, j}(X_k)L(X_k, ω)

p(X_k) (2)

The source of bias in most existing adaptive image-plane samplers is revealed here. Adaptive sampling without bias must avoid decisions to terminate sampling at an individual pixel, and instead look at the entire image plane to decide where a certain number of new samples will be cast. Every pixel with non-zero brightness must have non-zero probabil- ity of being chosen for a sample, regardless of its estimated error. This guarantee that all pixels will have chance to re- ceive samples to achieve unbiasedness.

We also note the Equation 2 can be broken into many integrals, one for the support of each pixel. Provided p(x) is known in each sub-domain, the global nature of p(x) is not important.

Figure 1 summarizes the final PMC-IP algorithm for a single frame:

1 Estimate the start-up kernel function α_k⁽⁰⁾from previous frames 2 for s = 0, · · · , S

4 Use DMS to allocate samples according to α_k^(s) 5 Generate samples from K_IP^(s)(x) and accumulate to image 6 Compute the perceptually-weighted variance image 7 Compute α_k^(s+1)for each pixel k

Fig. 1 The PMC-IP Algorithm for rendering a single frame.

3.1 Integrating the Sampler into a Global Rendering System

The PMC-IP is easily incorporated into a single rendering pipeline and allows us to improve the rendering efficiency

for ray-tracing based algorithms such as path tracing, pho- tonmapping. Figure 2 shows a modern plug-in style Monte Carlo rendering framework. The only core framework modi- fication required to support adaptive sampling is the addition of a feedback path from the output image generator back to the samplers, required to pass information from one sam- pling iteration back to the samplers for the next iteration.

The PMC-IP sampler also contains a tracking algorithm, de- scribed in Section 3.2, to guess the sample distribution from the previous rendering frames and sample distributions. k

The kernel function is the starting point in creating a PMC algorithm for adaptive image-plane sampling. We need a function that has adaptable parameters, is cheap to sample from, and supports stratification. This can be achieved with a mixture model of component distributions, h_{IP,(i, j)}(x), one for each pixel:

K_IP^(s)(x) =

∑

(i, j)∈P

α_{(i, j)}^(s) h_{IP,(i, j)}(x),

∑

(i, j)∈P

α_{(i, j)}^(s) = 1.

Where (i, j) is the pixel coordinate andP is the set of all pixels in this image. Each component is uniform over the domain of a single pixel integral. The parameters to the dis- tribution are all the α_{(i, j)}^(s) values, and these change for each iteration, s. We achieve an unbiased result if every α_{(i, j)}^(s) ≥ ε, where ε is a small positive constant (we use 0.01). We en- force this through the adaptive process, and the use of ε, rather than 0, provides some assurance that we will not over- look important contributions (referred to as defensive sam- pling [11]).

The use of a mixture as the kernel results in a D-kernel PMC [5] algorithm. Sampling from such a distribution is achieved by choosing a pixel, (i, j) according to the α_{(i, j)}^(s) , and then sampling from h_{IP,(i, j)}(x). The latter can be done with a low-discrepancy sampler within each pixel, giving sub-pixel stratification. Stratification across the entire image plane can be achieved through deterministic mixture sam- pling, which we describe shortly. The importance function p(x) in Equation 2 for a given pixel is p(x) = h_{IP,(i, j)}(x).

This can be derived by considering each pixel as an individ- ual integral and observing that only one mixture component has non-zero probability of contributing to each integral.

Notice that this kernel function is not conditional: K_IP(x^(s)|x^(s−1)) = KIP(x^(s)). Hence, for image-plane sampling we do not in-

clude a resampling step in the PMC algorithm because no samples are re-used. The knowledge gained from prior sam- ples is instead used to adapt the kernel function.

3.2 Predict a Good Initial Start-up Kernel

When observing the kernel function for each frame by us- ing PMC-IP algorithm, we realize that the high-probability regions should be temporally correlated among consecutive

4

Fig. 2 A block diagram of a plug-in style Monte Carlo rendering system, following Pharr and Humphreys [21]. The PMC-IP sampler replaces a uniform sample generator with the addition of a feedback path from the sample accumulator in order to calculate the perceptual variance. The PMC-IP sampler also contains a predictor to estimate the start-up sample distribution for rendering a single frame.

frames. Thus, if we can use motion analysis techniques to predict the movement of these regions, we can save the ex- tra cost of probing the entire image-plane to estimate the α_{(i, j)}. In addition, since the high-probability regions should also be regions with high variance in general Monte Carlo methods, to use the prediction can also reduce the tempo- ral inconsistency between noise by putting more samples in these regions.

For the current frame t, we first obtain a coarse rendering result ˆI^t using a small amount of initial samples per pixel such as 4 in our implementation. Then we predict the kernel function of the current frame α_{(i, j)}^t from that of the previous frame α_{(i, j)}^t−1by exploring the correspondence between ˆI^tand I^t−1as follows:

α_{(i, j)}^t = α_M^t−1t

t−1(i, j) (3)

where M_t−1^t describes the correspondence between α^t and α^t−1.

We approximate M_t−1^t by estimating the correspondence between images ˆI^tand I^t−1. Optical flow [17, 1] algorithms in computer vision community provide natural solutions to estimate this correspondence. However, since the coarse ren- dering result ˆI^t is noisy, optical flow estimation is not reli- able. We regularize the optical flow using a homography be- tween ˆI^tand I^t−1. A homography is a 2D perspective matrix described by 8 parameters. It describes the correspondence between ˆI^tand I^t−1as follows:



 si^t s j^t s



=





h11 h12h13

h₂₁ h₂₂h₂₃ h₃₁ h₃₂ 1







 i^t−1 j^t−1 1



 (4)

Vision research has provided rich literatures for estimating the homography [28]. We adopt a feature-based method. SIFT

features [16] are used due to their robustness to noise. Ide- ally, 4 pairs of correctly matched feature points are enough to estimate the homography. In practice, to get a robust es- timation, we extract from each image a dense set of SIFT features. ( ˆI^t is pre-processed using a median filter to re- duce the noise before estimating the homography.) We use a RANSAC [8] algorithm to obtain a robust estimation of the homography.

When the scene is a plane, or the camera only rotates along its optical center, the homography perfectly describes the correspondence between two images. In practice, when the difference between two consecutive frames are small, the homography can be a good approximation. However, when the camera motion or object motion is significant, the ho- mography is not accurate enough. To relieve this problem, we first detect the high density regions of the kernel α^t−1; then we estimate the homography within these regions be- tween two images, and use this new homography to describe the correspondence between the high-density regions. The insight of this strategy is two-fold: first, a homography is likely to be successful when modeling the correspondence between small regions; second, the high density region is more important than the other region, thus deserving accu- rate prediction. An example of the kernel predication is il- lustrated in Figure 3.

3.3 Adapting the PMC-IP Kernel

The adaption method is responsible for determining the value of each α_{(i, j)}^(s) given the populations from previous iterations and any information available from them, such as the im- age computed so far. We need to to define a α_{(i, j)}^(s) for every

5

Fig. 3 This is the prediction of the kernel function at frame 4. The images from left are result at frame 3, rough result at frame 4, the final kernel function at frame 3, the final predicted kernel function at frame 4, and the final kernel function at frame 4. We analyze the deviation by using the root mean squqre error which is 9.02 × 10⁻⁰⁷

pixel, with pixels that require more samples having higher high α_{(i, j)}^(s) for the component that covers the pixel.

An appropriate criteria assigns α_{(i, j)}^(s) proportional to the perceptually-weighted variance at each pixel. The algorithm tracks the sample variance in power seen among samples that contribute to each pixel. To account for perception, the result is divided by the threshold-versus-intensity function tvi(L) introduced by Ferweda et al. [7]. Normalization also accounts for ε.

α_{i, j}⁰ = σ²_{(i, j)} tvi(L_{(i, j)})

α_{i, j}^(s)= ε + (1 − ε)α_{(i, j)}⁰

∑(i⁰, j⁰)∈Pα_(i⁰0, j⁰)

The first iteration for the first frame of the algorithm samples uniformly over the image plane, so this criteria can always be computed. The first iteration for the following frames predict the kernel function according to the result of previous frame and the current rendering algorithm. The left images in Figure 4 show an example of an α_{(i, j)}⁽⁰⁾ map for a given initial image. The perceptual term in the error image prevents very high errors in both bright regions (a problem with unweighted variance) and dark areas (a problem with luminance-weighted variance).

Note that α_{(i, j)}^(s) ≥ ε, so there is a non-zero probability of generating a sample at any given image plane location.

This meets the requirement for importance sampling that the importance function is non-zero everywhere where the inte- grand is non-zero. Furthermore, as the total sample count ap- proaches infinity, the count at any pixel also approaches in- finity. Hence, with the correctly computed importance weights (Equation 2), the algorithm is unbiased.

3.4 Deterministic Mixture Sampling

Randomly sampling from the discrete distribution defined by the α_{(i, j)}^(s) produces excess noise — some pixels get far more or fewer samples than their expected value. This prob- lem is solved with deterministic mixture sampling, DMS, which is designed to give each component (pixel) a num- ber of samples roughly proportional to its α_{(i, j)}^(s) . Determin-

Image Method # SPP T(s) Err P-Eff

Buddha Uniform 10 58.1 0.625 0.027

PMC-IP 2+4+4 62.4 0.116 0.138

Box Uniform 16 163 0.545 0.011

Uniform 32 328 0.255 0.012

PMC-IP 4+6+6 169 0.182 0.033

Table 1 Measurements comparing PMC-IP and uniform image-plane sampling, for equal total sample counts. The P-Eff is perceptual effi- ciency used in [13]. The Buddha image computed direct lighting with the MIS method, with a total of 8 lighting samples for each pixel sam- ple. PMC-IP sampling improves the perceptual-based RMS error by a factor of 5.4 over uniform sampling with only 7.5% more computa- tion time. It corresponds to an improvement in efficiency of 5.01. The Cornell Box images use path tracing to compute global illumination including caustics. Comparing with images of 16 SPPs, PMC-IP im- proves the efficiency by a factor of 2.65.

istic mixture sampling is unbiased and always gives lower variance when compared to random mixture sampling, as proven by Hesterberg [11].

The number of samples per iteration, N, (the population size) is fixed at a small multiple of the number of pixels. We typically use 4 samples per pixel, which balances between spending too much effort on any one iteration and the over- head of computing a new set of kernel parameters. For each pixel, the deterministic sampler computes n⁰_{(i, j)}= Nα_{(i, j)}, the target number of samples for that pixel. It takes bn⁰_{(i, j)}c samples from each pixel (i, j)’s component. The remaining un-allocated samples are sampled from the residual distri- butionwith probability n⁰_{(i, j)}− bn⁰_{(i, j)}c at each pixel (suitably normalized).

4 Results

This section presents the rendering results when we apply our PMC-IP algorithm to render a single frame of scenes and animation scenes by plugging in our sampler into the modern global illumination framework to demonstrate the usefulness of our algorithm

6

Fig. 4 A comparison between adaptive and uniform image-plane sampling on a direct lighting example. Leftmost is the initial image for PMC-IP sampling, and the α_k⁽⁰⁾image. The initial image used 2 samples per pixel. The next image is the result of PMC-IP sampling with two iterations at 4 SPPs on average. Center is a 10 SPPs image uniformly distributed. The zooms show the shadow near the Buddha’s base (PMC-IP top, uniform bottom). To the right are the corresponding variance images. Note that the variance image for the PMC-IP sampler has few high variance regions, and has a lower contrast in general, representing a more even distribution of error.

Fig. 5 A Cornell Box image computed using the PMC-IP algorithm. The left image is the final result using PMC-IP algorithm with 64 iterations, with each iteration averaging 4 SPPs. It is easier to find converged values in diffuse regions than in the caustic region. Thus, we compare the results by focusing on this region. The images in the second row on the top, from top to down, are the cropped images of the caustic region computed using non-adaptive path tracing with 16, 32, 64 and 128 SPPs. The images in the third row on the top, from top to down, are intermediate results from the adaptive algorithm at 4, 8, 16 and 32 iterations when computing the top image. The right bottom demonstrates that our adaptive sampler produces better visual results at lower sample counts: on the top is the result from 256 SPPs, un-adapted, and on the bottom image is the result of 16 adapting iterations at an average of 4 SPPs per iteration.

4.1 Static Image Rendering

Adaptive image-plane sampling can be used in many situa- tions where pixel samples are required and an iterative algo- rithm can be employed. We have implemented it in the con- texts of direct lighting using a Multiple Importance Sampler

(MIS) and global illumination with path tracing, and as part of a complete photon mapping system.

Figure 4 shows the Buddha direct lighting example. The surface is diffuse with an area light source. Each pixel sam- ple used 8 illumination samples, and the images were ren- dered at 256×512, with statistics presented in Table 1.We