**레이 트레이싱: 남은 인생**
[Peter Shirley][], [Trevor David Black][], [Steve Hollasch][]
버전 4.0.2, 2025-04-25
Copyright 2018-2024 Peter Shirley. 모든 권리를 보유합니다.
개요
====================================================================================================
_주말 동안 배우는 레이 트레이싱_ 과 _레이 트레이싱: 다음 주_ 에서 여러분은 “진짜” 레이 트레이서를 만들었습니다.
의욕이 있다면 이 책들에 담긴 소스와 정보를 활용하여 원하는 어떤 시각 효과든 구현할 수 있습니다. 이 소스는 소규모 취미 프로젝트용 레이 트레이서를 구축하는 데 의미 있고 견고한 기반을 제공합니다. 상용 레이 트레이서에서 볼 수 있는 시각 효과 대부분은 앞의 두 책에서 설명한 기법에 기반합니다. 그러나 수학적 기반이 없으면 표면 하 산란이나 중첩된 유전체처럼 점점 더 복잡해지는 시각 효과를 추가하는 능력은 크게 제한될 것입니다. 이 책에서는 여러분이 이 분야에 관심이 매우 많은 학생이거나 레이 트레이싱과 관련된 진로를 추구하는 사람이라고 가정합니다. 우리는 매우 본격적인 레이 트레이서를 만드는 데 필요한 수학을 깊이 파고들 것입니다. 이 책을 마치면 영화, 텔레비전, 제품 디자인, 건축 산업 등 여러 인기 있는 분야에서 사용되는 다양한 상용 레이 트레이서를 사용하고 수정할 수 있는 역량을 충분히 갖추게 될 것입니다.
이 짧은 책에서 다루지 않는 것은 정말 많습니다. 예를 들어, 몬테카를로 렌더링 프로그램을 작성하는 방법은 매우 많지만, 저는 그중 하나만 깊이 다룹니다. 그림자 광선은 다루지 않으며(대신 광선이 광원 쪽으로 향할 가능성을 높이는 방식을 택했습니다), 양방향 방법, 메트로폴리스 방법, 광자 매핑도 다루지 않습니다. 이러한 기법 중 상당수는 이른바 "본격적인 레이 트레이서"에서 찾아볼 수 있지만, 이 분야의 개념, 수학, 용어를 다루는 것이 더 중요하기 때문에 여기서는 다루지 않습니다. 저는 이 책을 앞으로 여러 차례 접하게 될 깊이 있는 학습 중 첫 번째가 되어야 할 책으로 생각하며, 이 책은 이러한 기법과 그 밖의 흥미로운 기법을 공부하는 데 필요한 개념, 수학, 용어 중 일부를 여러분에게 갖춰 줄 것입니다.
여러분도 저처럼 이 수학을 매력적으로 느끼시기를 바랍니다.
이 프로젝트, GitHub 저장소, 디렉터리 구조, 빌드 및 실행, 그리고 수정 사항과 기여를 작성하거나 참조하는 방법에 대한 정보는 [프로젝트 README][readme] 파일을 참조하십시오.
이전과 마찬가지로, 프로젝트 관련 추가 자료는 [추가 읽을거리 위키 페이지][wiki-further]를 참조하십시오.
이 책들은 브라우저에서 바로 인쇄해도 잘 출력되도록 서식이 지정되어 있습니다. 또한 [각 릴리스마다][releases] "Assets" 섹션에 각 책의 PDF를 포함하고 있습니다.
이 프로젝트에 도움을 주신 모든 분께 감사드립니다. 이분들은 이 책 끝의 감사의 글 절에서 확인하실 수 있습니다.
간단한 몬테카를로 프로그램
====================================================================================================
가장 간단한 몬테카를로 프로그램 중 하나부터 시작하겠습니다. 몬테카를로 프로그램에 익숙하지 않으시다면, 잠시 멈추고 기초를 짚고 넘어가는 것이 좋겠습니다. 무작위 알고리즘에는 몬테카를로와 라스베이거스, 두 가지 종류가 있습니다. 무작위 알고리즘은 컴퓨터 그래픽스 곳곳에서 찾아볼 수 있으므로, 탄탄한 기초를 다져 두는 것도 나쁘지 않습니다. 무작위 알고리즘은 계산 과정에서 어느 정도의 무작위성을 사용합니다. 라스베이거스 무작위 알고리즘은 항상 올바른 결과를 내는 반면, 몬테카를로 알고리즘은 올바른 결과를 낼 _수도_ 있지만--자주 틀리기도 합니다! 하지만 레이 트레이싱처럼 특히 복잡한 문제에서는, 완벽하게 정확한 답을 얻는 것보다 합리적인 시간 안에 답을 얻는 것을 더 중요하게 여길 수 있습니다. 라스베이거스 알고리즘은 결국 올바른 결과에 도달하지만, 거기까지 얼마나 오래 걸릴지는 크게 보장할 수 없습니다. 라스베이거스 알고리즘의 고전적인 예는 _퀵소트_ 정렬 알고리즘입니다. 퀵소트 알고리즘은 항상 완전히 정렬된 리스트를 만들고 끝나지만, 완료하는 데 걸리는 시간은 무작위입니다. 라스베이거스 알고리즘의 또 다른 좋은 예는 단위 원판 안의 무작위 점을 고르는 데 사용하는 코드입니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
inline vec3 random_in_unit_disk() {
while (true) {
auto p = vec3(random_double(-1,1), random_double(-1,1), 0);
if (p.length_squared() < 1)
return p;
}
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [las-vegas-algo]: [vec3.h] 라스베이거스 알고리즘]
이 코드는 결국에는 항상 단위 원판 안의 무작위 점에 도달하지만, 얼마나 걸릴지는 미리 말할 수 없습니다. 단 1번의 반복으로 끝날 수도 있고, 2번, 3번, 4번, 또는 그보다 더 오래 걸릴 수도 있습니다. 반면 몬테카를로 프로그램은 답에 대한 통계적 추정값을 제공하며, 이 추정값은 오래 실행할수록 점점 더 정확해집니다. 즉, 어느 시점에서는 답이 _충분히_ 정확하다고 판단하고 그만둘 수 있다는 뜻입니다. 단순한 프로그램이 노이즈는 있지만 점점 더 나아지는 답을 내놓는다는 이 기본적인 특성이 바로 몬테카를로의 핵심이며, 높은 정확도가 필요하지 않은 그래픽스 같은 응용 분야에 특히 적합합니다.
파이 추정하기
--------------
몬테카를로 알고리즘의 대표적인 예는 $\pi$를 추정하는 것이므로, 그것을 해 보겠습니다. $\pi$를 추정하는 방법은 많으며, 그중 _뷔퐁의 바늘 문제_ 는 고전적인 사례 연구입니다. 뷔퐁의 바늘 문제에서는 폭이 모두 같은 평행한 마룻널 조각들로 이루어진 바닥이 주어집니다. 바늘을 바닥에 무작위로 떨어뜨렸을 때, 바늘이 두 널빤지에 걸쳐 놓일 확률은 얼마입니까? (이 문제에 대한 더 많은 정보는 간단한 인터넷 검색으로 찾을 수 있습니다.)
이 방법에서 영감을 받은 변형을 해 보겠습니다. 정사각형 안에 내접한 원이 있다고 가정해 보십시오:
![Figure [circ-square]: 정사각형 안의 원으로 $\pi$ 추정하기
](../images/fig-3.01-circ-square.jpg)
이제 정사각형 안에서 무작위 점들을 고른다고 가정해 보십시오. 그 무작위 점들 중 원 안에 들어가는 점의 비율은 원의 넓이에 비례해야 합니다. 정확한 비율은 사실 원의 넓이와 정사각형 넓이의 비여야 합니다:
$$ \frac{\pi r^2}{(2r)^2} = \frac{\pi}{4} $$
$r$이 상쇄되므로, 계산하기 편한 값을 아무것이나 고를 수 있습니다. 원점을 중심으로 하는 $r=1$을 사용하겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
int main() {
std::cout << std::fixed << std::setprecision(12);
int inside_circle = 0;
int N = 100000;
for (int i = 0; i < N; i++) {
auto x = random_double(-1,1);
auto y = random_double(-1,1);
if (x*x + y*y < 1)
inside_circle++;
}
std::cout << "Estimate of Pi = " << (4.0 * inside_circle) / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [estpi-1]: [pi.cc] $\pi$ 추정하기 -- 첫 번째 버전]
계산되는 $\pi$ 값은 초기 난수 시드에 따라 컴퓨터마다 다릅니다. 제 컴퓨터에서는 `Estimate of Pi = 3.143760000000`이라는 답이 나옵니다.
수렴 과정 보기
--------------------
프로그램이 끝없이 실행되면서 현재까지의 추정값을 계속 출력하도록 바꾸면 다음과 같습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
int main() {
std::cout << std::fixed << std::setprecision(12);
int inside_circle = 0;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
int runs = 0;
while (true) {
runs++;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
auto x = random_double(-1,1);
auto y = random_double(-1,1);
if (x*x + y*y < 1)
inside_circle++;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
if (runs % 100000 == 0)
std::cout << "\rEstimate of Pi = " << (4.0 * inside_circle) / runs;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ delete
std::cout << "Estimate of Pi = " << (4.0 * inside_circle) / N << '\n';
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [estpi-2]: [pi.cc] $\pi$ 추정하기 -- 두 번째 버전]
층화 샘플(지터링)
-------------------------------
값은 매우 빠르게 $\pi$ 근처에 도달하지만, 그 뒤로는 더 천천히 그 값에 가까워집니다. 이는 _수확 체감의 법칙_ 의 한 예로, 각 샘플이 주는 도움은 직전 샘플보다 적습니다. 이것이 몬테카를로 방법의 가장 나쁜 점입니다. 샘플을 _층화_ 하면(흔히 _지터링_ 이라고 부릅니다) 이러한 수확 체감을 완화할 수 있습니다. 이 방법에서는 무작위 샘플을 취하는 대신 격자를 만들고 각 칸 안에서 샘플을 하나씩 취합니다:
![Figure [jitter]: 지터링된 점으로 영역 샘플링하기](../images/fig-3.02-jitter.jpg)
이렇게 하면 샘플 생성 방식이 바뀌지만, 격자를 알아야 하므로 샘플을 몇 개 취할지 미리 알아야 합니다. 백만 개를 취해서 두 가지 방법으로 모두 시도해 보겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
int main() {
std::cout << std::fixed << std::setprecision(12);
int inside_circle = 0;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
int inside_circle_stratified = 0;
int sqrt_N = 1000;
for (int i = 0; i < sqrt_N; i++) {
for (int j = 0; j < sqrt_N; j++) {
auto x = random_double(-1,1);
auto y = random_double(-1,1);
if (x*x + y*y < 1)
inside_circle++;
x = 2*((i + random_double()) / sqrt_N) - 1;
y = 2*((j + random_double()) / sqrt_N) - 1;
if (x*x + y*y < 1)
inside_circle_stratified++;
}
}
std::cout
<< "Regular Estimate of Pi = "
<< (4.0 * inside_circle) / (sqrt_N*sqrt_N) << '\n'
<< "Stratified Estimate of Pi = "
<< (4.0 * inside_circle_stratified) / (sqrt_N*sqrt_N) << '\n';
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [estpi-3]: [pi.cc] $\pi$ 추정하기 -- 세 번째 버전]
제 컴퓨터에서는 다음과 같은 결과가 나옵니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Regular Estimate of Pi = 3.141184000000
Stratified Estimate of Pi = 3.141460000000
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
참고로 파이의 소수점 이하 첫 12자리는 다음과 같습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3.141592653589
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
흥미롭게도 층화 방법은 더 나을 뿐만 아니라 점근적으로도 더 빠른 속도로 수렴합니다! 안타깝게도 이러한 이점은 문제의 차원이 높아질수록 줄어듭니다(예를 들어 3D 구의 부피를 구하는 버전에서는 차이가 더 작을 것입니다). 이를 _차원의 저주_ 라고 합니다. 레이 트레이싱은 차원이 매우 높은 알고리즘으로, 반사가 일어날 때마다 $\phi_o$와 $\theta_o$라는 두 개의 새로운 차원이 추가됩니다. 이 책에서는 출력 반사 각도를 층화하지 않습니다. 단지 다루기에 조금 너무 복잡하기 때문이며, 현재 이 분야에서는 흥미로운 연구가 많이 진행되고 있습니다.
그 중간 단계로, 각 픽셀 위치 주변의 샘플링 위치를 층화하겠습니다. 익숙한 코넬 박스 장면부터 시작하겠습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include "camera.h"
#include "hittable_list.h"
#include "material.h"
#include "quad.h"
#include "sphere.h"
int main() {
hittable_list world;
auto red = make_shared(color(.65, .05, .05));
auto white = make_shared(color(.73, .73, .73));
auto green = make_shared(color(.12, .45, .15));
auto light = make_shared(color(15, 15, 15));
// Cornell box sides
world.add(make_shared(point3(555,0,0), vec3(0,0,555), vec3(0,555,0), green));
world.add(make_shared(point3(0,0,555), vec3(0,0,-555), vec3(0,555,0), red));
world.add(make_shared(point3(0,555,0), vec3(555,0,0), vec3(0,0,555), white));
world.add(make_shared(point3(0,0,555), vec3(555,0,0), vec3(0,0,-555), white));
world.add(make_shared(point3(555,0,555), vec3(-555,0,0), vec3(0,555,0), white));
// Light
world.add(make_shared(point3(213,554,227), vec3(130,0,0), vec3(0,0,105), light));
// Box 1
shared_ptr box1 = box(point3(0,0,0), point3(165,330,165), white);
box1 = make_shared(box1, 15);
box1 = make_shared(box1, vec3(265,0,295));
world.add(box1);
// Box 2
shared_ptr box2 = box(point3(0,0,0), point3(165,165,165), white);
box2 = make_shared(box2, -18);
box2 = make_shared(box2, vec3(130,0,65));
world.add(box2);
camera cam;
cam.aspect_ratio = 1.0;
cam.image_width = 600;
cam.samples_per_pixel = 64;
cam.max_depth = 50;
cam.background = color(0,0,0);
cam.vfov = 40;
cam.lookfrom = point3(278, 278, -800);
cam.lookat = point3(278, 278, 0);
cam.vup = vec3(0, 1, 0);
cam.defocus_angle = 0;
cam.render(world);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [estpi-3]: [main.cc] 코넬 박스 다시 보기]
이 프로그램을 실행해 층화하지 않은 렌더링 결과를 만들고 비교를 위해 저장하십시오.
이제 층화 샘플링 절차를 구현하기 위해 다음과 같이 변경하십시오:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
public:
...
void render(const hittable& world) {
initialize();
std::cout << "P3\n" << image_width << ' ' << image_height << "\n255\n";
for (int j = 0; j < image_height; j++) {
std::clog << "\rScanlines remaining: " << (image_height - j) << ' ' << std::flush;
for (int i = 0; i < image_width; i++) {
color pixel_color(0,0,0);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
for (int s_j = 0; s_j < sqrt_spp; s_j++) {
for (int s_i = 0; s_i < sqrt_spp; s_i++) {
ray r = get_ray(i, j, s_i, s_j);
pixel_color += ray_color(r, max_depth, world);
}
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
write_color(std::cout, pixel_samples_scale * pixel_color);
}
}
std::clog << "\rDone. \n";
}
private:
int image_height; // Rendered image height
double pixel_samples_scale; // Color scale factor for a sum of pixel samples
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
int sqrt_spp; // Square root of number of samples per pixel
double recip_sqrt_spp; // 1 / sqrt_spp
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
point3 center; // Camera center
...
void initialize() {
image_height = int(image_width / aspect_ratio);
image_height = (image_height < 1) ? 1 : image_height;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
sqrt_spp = int(std::sqrt(samples_per_pixel));
pixel_samples_scale = 1.0 / (sqrt_spp * sqrt_spp);
recip_sqrt_spp = 1.0 / sqrt_spp;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
center = lookfrom;
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
ray get_ray(int i, int j, int s_i, int s_j) const {
// Construct a camera ray originating from the defocus disk and directed at a randomly
// sampled point around the pixel location i, j for stratified sample square s_i, s_j.
auto offset = sample_square_stratified(s_i, s_j);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
auto pixel_sample = pixel00_loc
+ ((i + offset.x()) * pixel_delta_u)
+ ((j + offset.y()) * pixel_delta_v);
auto ray_origin = (defocus_angle <= 0) ? center : defocus_disk_sample();
auto ray_direction = pixel_sample - ray_origin;
auto ray_time = random_double();
return ray(ray_origin, ray_direction, ray_time);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
vec3 sample_square_stratified(int s_i, int s_j) const {
// Returns the vector to a random point in the square sub-pixel specified by grid
// indices s_i and s_j, for an idealized unit square pixel [-.5,-.5] to [+.5,+.5].
auto px = ((s_i + random_double()) * recip_sqrt_spp) - 0.5;
auto py = ((s_j + random_double()) * recip_sqrt_spp) - 0.5;
return vec3(px, py, 0);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
vec3 sample_square() const {
...
}
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [render-estpi-3]: [camera.h] 픽셀 내부의 샘플 층화하기]
층화를 적용한 결과를 비교하면 다음과 같습니다:
눈을 가늘게 뜨고 보면 평면의 가장자리와 상자의 가장자리에서 더 선명한 대비를 확인할 수 있을 것입니다. 이 효과는 변화의 빈도가 높은 위치에서 더 두드러지게 나타납니다. 고주파 변화는 높은 정보 밀도로 생각할 수도 있습니다. 우리의 코넬 박스 장면에서는 모든 재질이 무광이고 머리 위에 부드러운 면광원이 있으므로, 정보 밀도가 높은 위치는 물체의 가장자리뿐입니다. 이 효과는 텍스처와 반사 재질을 사용할 때 더 뚜렷해질 것입니다.
단일 반사나 그림자 처리, 또는 엄밀한 2D 문제를 다룬다면 반드시 층화를 사용하는 것이 좋습니다.
1차원 몬테카를로 적분
====================================================================================================
우리가 변형한 뷔퐁의 바늘 문제는 원의 넓이와 외접 정사각형의 넓이의 비를 구함으로써 $\pi$를 계산하는 방법입니다:
$$ \frac{\operatorname{area}(\mathit{circle})}{\operatorname{area}(\mathit{square})}
= \frac{\pi}{4}
$$
우리는 외접 정사각형 안에서 무작위 점을 여러 개 골라, 그중 단위 원 안에도 들어가는 점의 비율을 세었습니다. 이 비율은 점이 추가될수록 $\frac{\pi}{4}$에 가까워지는 추정값이었습니다. 원의 넓이를 모르더라도 위의 비를 이용해 여전히 구할 수 있습니다. 단위 원과 외접 정사각형의 넓이의 비가 $\frac{\pi}{4}$라는 것을 알고, 외접 정사각형의 넓이가 $4r^2$라는 것도 알고 있으므로, 이 두 값을 이용해 원의 넓이를 구할 수 있습니다:
$$ \frac{\operatorname{area}(\mathit{circle})}{\operatorname{area}(\mathit{square})}
= \frac{\pi}{4}
$$
$$ \frac{\operatorname{area}(\mathit{circle})}{(2r)^2} = \frac{\pi}{4} $$
$$ \operatorname{area}(\mathit{circle}) = \frac{\pi}{4} 4r^2 $$
$$ \operatorname{area}(\mathit{circle}) = \pi r^2 $$
반지름이 $r = 1$인 원을 선택하면 다음을 얻습니다:
$$ \operatorname{area}(\mathit{circle}) = \pi $$
위의 작업은 원의 넓이를 구하는 수단으로서와 마찬가지로 $\pi$를 구하는 수단으로서도 똑같이 유효합니다. 따라서 우리 pi 프로그램의 초기 버전 중 하나에서 다음과 같이 대체할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ delete
std::cout << "Estimate of Pi = " << (4.0 * inside_circle) / N << '\n';
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
std::cout << "Estimated area of unit circle = " << (4.0 * inside_circle) / N << '\n';
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [estunitcircle]: [pi.cc] 단위 원의 넓이 추정하기]
기댓값
--------------
한 걸음 물러서서 우리의 몬테카를로 알고리즘을 조금 더 일반적으로 생각해 봅시다.
다음을 모두 가지고 있다고 가정해 보겠습니다:
1. 원소 $x_i$를 포함하는 값의 리스트 $X$:
$$ X = (x_0, x_1, ..., x_{N-1}) $$
2. 리스트의 원소를 입력으로 받는 연속 함수 $f(x)$:
$$ y_i = f(x_i) $$
3. 리스트 $X$를 입력으로 받아 리스트 $Y$를 출력으로 내놓는 함수 $F(X)$:
$$ Y = F(X) $$
4. 여기서 출력 리스트 $Y$는 원소 $y_i$를 가집니다:
$$ Y = (y_0, y_1, ..., y_{N-1}) = (f(x_0), f(x_1), ..., f(x_{N-1})) $$
위의 내용을 모두 가정하면, 다음과 같이 리스트 $Y$의 산술 평균, 즉 평균을 구할 수 있습니다:
$$ \operatorname{average}(Y_i) = E[Y] = \frac{1}{N} \sum_{i=0}^{N-1} y_i $$
$$ = \frac{1}{N} \sum_{i=0}^{N-1} f(x_i) $$
$$ = E[F(X)] $$
여기서 $E[Y]$를 $Y$의 _기댓값_ 이라고 합니다.
여기서 _평균값_ 과 _기댓값_ 의 미묘한 차이에 유의하십시오:
- 하나의 집합에는 그 집합에서 선택한 서로 다른 부분집합이 많이 있을 수 있습니다. 각 부분집합에는 _평균값_ 이 있으며, 이는 선택된 모든 항목의 합을 선택된 항목의 수로 나눈 값입니다. 주어진 항목은 부분집합에서 0번, 1번 또는 여러 번 나타날 수 있다는 점에 유의하십시오.
- 하나의 집합에는 _기댓값_ 이 단 하나만 있습니다: 집합의 _모든_ 원소의 합을 그 집합의 전체 항목 수로 나눈 값입니다. 달리 말하면, _기댓값_ 은 집합의 _모든_ 원소의 _평균값_ 입니다.
집합에서 뽑는 무작위 샘플의 수가 늘어날수록 집합의 평균값이 기댓값으로 수렴한다는 점에 유의하는 것이 중요합니다.
$x_i$의 값이 모든 $i$에 대해 $ a \leq x_i \leq b $를 만족하도록 연속 구간 $[a,b]$에서 무작위로 선택된다면, $E[F(X)]$는 같은 구간 $ a \leq x' \leq b $에서 연속 함수 $f(x')$의 평균을 근사하게 됩니다.
$$ E[f(x') | a \leq x' \leq b] \approx E[F(X) | X =
\{\small x_i | a \leq x_i \leq b \normalsize \} ] $$
$$ \approx E[Y = \{\small y_i = f(x_i) | a \leq x_i \leq b \normalsize \} ] $$
$$ \approx \frac{1}{N} \sum_{i=0}^{N-1} f(x_i) $$
샘플 수 $N$을 취하고 $N$이 $\infty$로 갈 때의 극한을 취하면 다음을 얻습니다:
$$ E[f(x') | a \leq x' \leq b] = \lim_{N \to \infty} \frac{1}{N} \sum_{i=0}^{N-1} f(x_i) $$
연속 구간 $[a,b]$ 내에서 연속 함수 $f(x')$의 기댓값은 구간 내의 무한히 많은 무작위 점을 합산함으로써 완벽하게 표현될 수 있습니다. 이 점의 수가 $\infty$에 가까워질수록 출력값의 평균은 정확한 답에 수렴합니다. 이것이 몬테카를로 알고리즘입니다.
무작위 점을 샘플링하는 것만이 구간에 대한 기댓값을 구하는 유일한 방법은 아닙니다. 샘플링 점을 어디에 배치할지 직접 선택할 수도 있습니다. 구간 $[a,b]$에 대해 $N$개의 샘플이 있다면, 구간 전체에 걸쳐 점을 균등한 간격으로 배치하도록 선택할 수 있습니다:
$$ x_i = a + i \Delta x $$
$$ \Delta x = \frac{b - a}{N} $$
그런 다음 이들의 기댓값을 구합니다:
$$ E[f(x') | a \leq x' \leq b] \approx \frac{1}{N} \sum_{i=0}^{N-1} f(x_i)
\Big|_{x_i = a + i \Delta x} $$
$$ E[f(x') | a \leq x' \leq b] \approx \frac{\Delta x}{b - a} \sum_{i=0}^{N-1} f(x_i)
\Big|_{x_i = a + i \Delta x} $$
$$ E[f(x') | a \leq x' \leq b] \approx \frac{1}{b - a} \sum_{i=0}^{N-1} f(x_i) \Delta x
\Big|_{x_i = a + i \Delta x} $$
$N$이 $\infty$에 가까워질 때의 극한을 취합니다
$$ E[f(x') | a \leq x' \leq b] = \lim_{N \to \infty} \frac{1}{b - a} \sum_{i=0}^{N-1}
f(x_i) \Delta x \Big|_{x_i = a + i \Delta x} $$
물론 이것은 그저 일반적인 적분입니다:
$$ E[f(x') | a \leq x' \leq b] = \frac{1}{b - a} \int_{a}^{b} f(x) dx $$
기초 미적분학 수업을 떠올려 보면, 함수의 적분은 해당 구간에서 곡선 아래의 넓이입니다:
$$ \operatorname{area}(f(x), a, b) = \int_{a}^{b} f(x) dx$$
따라서 구간에 대한 평균은 그 구간에서 곡선 아래의 넓이와 본질적으로 연결되어 있습니다.
$$ E[f(x) | a \leq x \leq b] = \frac{1}{b - a} \cdot \operatorname{area}(f(x), a, b) $$
함수의 적분과 그 함수의 몬테카를로 샘플링은 모두 특정 구간에 대한 평균을 구하는 데 사용될 수 있습니다. 적분은 구간을 무한히 작게 나눈 무한히 많은 조각의 합으로 평균을 구하는 반면, 몬테카를로 알고리즘은 구간 내에서 계속 늘어나는 무작위 샘플 점의 합을 구함으로써 같은 평균을 근사합니다. 물체 내부에 들어가는 점의 수를 세는 것만이 그 평균이나 넓이를 측정하는 유일한 방법은 아닙니다. 적분 또한 이러한 목적을 위한 일반적인 수학적 도구입니다. 문제에 대한 닫힌 형식이 존재한다면, 적분이 문제를 공식화하는 가장 자연스럽고 깔끔한 방법인 경우가 많습니다.
몇 가지 예제가 도움이 될 것이라고 생각합니다.
x² 적분하기
---------------
고전적인 적분 하나를 살펴보겠습니다:
$$ I = \int_{0}^{2} x^2 dx $$
이것은 적분을 사용하여 풀 수 있습니다:
$$ I = \frac{1}{3} x^3 \Big|_{0}^{2} $$
$$ I = \frac{1}{3} (2^3 - 0^3) $$
$$ I = \frac{8}{3} $$
또는 몬테카를로 접근 방식을 사용하여 이 적분을 풀 수도 있습니다. 컴퓨터 과학식 표기법으로는 이것을 다음과 같이 쓸 수 있습니다:
$$ I = \operatorname{area}( x^2, 0, 2 ) $$
다음과 같이 쓸 수도 있습니다:
$$ E[f(x) | a \leq x \leq b] = \frac{1}{b - a} \cdot \operatorname{area}(f(x), a, b) $$
$$ \operatorname{average}(x^2, 0, 2) = \frac{1}{2 - 0} \cdot \operatorname{area}( x^2, 0, 2 ) $$
$$ \operatorname{average}(x^2, 0, 2) = \frac{1}{2 - 0} \cdot I $$
$$ I = 2 \cdot \operatorname{average}(x^2, 0, 2) $$
몬테카를로 접근 방식:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
int main() {
int a = 0;
int b = 2;
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
auto x = random_double(a, b);
sum += x*x;
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "I = " << (b - a) * (sum / N) << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-xsq-1]: [integrate_x_sq.cc] $x^2$ 적분하기]
예상대로 이것은 적분으로 얻는 정확한 답, _즉_ $I = 2.666… = \frac{8}{3}$에 근사한 값을 산출합니다. 이 예제를 보고 실제로는 적분이 몬테카를로보다 훨씬 수고가 적다고 지적하셔도 옳습니다. 함수가 $f(x) = x^2$인 경우에는 그럴 수 있지만, $f(x) = \sin^5(x)$처럼 적분보다 몬테카를로로 푸는 것이 더 간단할 수 있는 함수도 많이 존재합니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
for (int i = 0; i < N; i++) {
auto x = random_double(a, b);
sum += std::pow(std::sin(x), 5.0);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-sin5]: $\sin^5$ 적분하기]
$f(x) = \ln(\sin(x))$처럼 해석적 적분이 존재하지 않는 함수에도 몬테카를로 알고리즘을 사용할 수 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
for (int i = 0; i < N; i++) {
auto x = random_double(a, b);
sum += std::log(std::sin(x));
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-ln-sin]: $\ln(\sin)$ 적분하기]
그래픽스에서는 명시적으로 적을 수는 있지만 해석적으로 적분하기가 복잡한 함수를 자주 다룹니다. 그만큼 자주 평가는 _할 수 있지만_ 명시적으로 적을 수는 _없는_ 함수도 다루며, _오직_ 확률적으로만 평가할 수 있는 함수를 마주하는 경우도 빈번합니다. 앞선 두 권의 책에서 다룬 `ray_color` 함수는 확률적으로만 결정할 수 있는 함수의 한 예입니다. 임의의 위치에서 모든 방향으로 어떤 색이 보이는지는 알 수 없지만, 특정한 한 위치에서 특정한 한 방향으로 어떤 색이 보이는지는 통계적으로 추정할 수 있습니다.
밀도 함수
------------------
앞선 두 권의 책에서 작성한 `ray_color` 함수는 단순하면서도 우아하지만, 상당히 _중대한_ 문제를 안고 있습니다. 작은 광원은 노이즈를 너무 많이 만들어 냅니다. 이는 균일 샘플링이 이러한 광원을 충분히 자주 샘플링하지 못하기 때문입니다. 광원은 광선이 그쪽으로 산란될 때에만 샘플링되는데, 작은 광원이나 멀리 떨어진 광원의 경우에는 그럴 가능성이 낮을 수 있습니다. 배경색이 검은색이라면, 장면에서 실제 빛의 원천은 장면 곳곳에 실제로 배치된 광원뿐입니다. 표면 위의 가까운 두 지점에서 교차하는 두 광선이 있을 수 있는데, 하나는 무작위로 광원 쪽으로 반사되고 다른 하나는 그렇지 않을 수 있습니다. 광원 쪽으로 반사된 광선은 매우 밝은 색으로 나타날 것입니다. 다른 곳으로 반사된 광선은 매우 어두운 색으로 나타날 것입니다. 두 밝기는 사실 그 중간 어딘가여야 합니다. 두 광선을 모두 광원 쪽으로 향하게 하면 이 문제를 줄일 수 있겠지만, 그렇게 하면 장면이 부정확하게 밝아집니다.
어떤 광선이든 보통은 카메라에서 출발하여 장면을 거쳐 광원에서 끝나도록 추적합니다. 하지만 이 광선을 반대로 광원에서 출발하여 장면을 거쳐 카메라에서 끝나도록 추적한다고 상상해 보십시오. 이 광선은 밝은 세기로 출발하여 장면 안에서 튕길 때마다 에너지를 잃게 됩니다. 결국 이 광선은 여러 표면에서 반사되면서 어두워지고 색이 입혀진 채로 카메라에 도달합니다. 이제 이 광선이 가능한 한 빨리 카메라 쪽으로 튕기도록 강제된다고 상상해 보십시오. 연속된 튕김으로 어두워지지 않았기 때문에 이 광선은 부정확하게 밝게 나타날 것입니다. 이는 광원 쪽으로 더 많은 무작위 샘플을 보내는 것과 유사합니다. 이렇게 하면 밝은 픽셀 옆에 어두운 픽셀이 있는 문제를 해결하는 데 크게 도움이 되겠지만, 결국 _모든_ 픽셀을 밝게 만들어 버릴 것입니다.
과다 샘플링을 보정하도록 해당 샘플의 가중치를 낮추면 이 부정확성을 없앨 수 있습니다. 이 보정은 어떻게 할까요? 먼저 _확률 밀도 함수_ 의 개념을 이해해야 합니다. 그런데 _확률 밀도 함수_ 의 개념을 이해하려면, 먼저 _밀도 함수_ 가 무엇인지 알아야 합니다.
_밀도 함수_ 는 히스토그램의 연속 버전일 뿐입니다. 다음은 위키백과의 히스토그램 페이지에 실린 히스토그램의 예입니다:
![Figure [histogram]: 히스토그램 예시](../images/fig-3.03-histogram.jpg)
데이터 소스에 항목이 더 많다면 구간의 수는 그대로이지만, 각 구간에 속하는 항목의 빈도는 더 높아질 것입니다. 데이터를 더 많은 구간으로 나눈다면 구간은 더 많아지지만, 각 구간에 속하는 항목의 빈도는 더 낮아질 것입니다. 구간의 수를 무한대로 늘린다면 빈도가 0인 구간이 무한히 많아질 것입니다. 이 문제를 해결하기 위해 _이산 함수_ 인 히스토그램을 _이산 밀도 함수_ 로 대체하겠습니다. _이산 밀도 함수_ 는 y축을 각 구간의 전체 개수가 아니라 전체에 대한 비율이나 백분율, _즉_ 밀도로 정규화한다는 점에서 _이산 함수_ 와 다릅니다. _이산 함수_ 를 _이산 밀도 함수_ 로 변환하는 것은 간단합니다:
$$ \text{Density of Bin i} = \frac{\text{Number of items in Bin i}}
{\text{Number of items total}} $$
_이산 밀도 함수_ 를 얻고 나면 이산 값을 연속 값으로 바꾸어 이를 _밀도 함수_ 로 변환할 수 있습니다.
$$ \text{Bin Density} = \frac{(\text{Fraction of trees between height }H\text{ and }H’)}
{(H-H’)} $$
즉, _밀도 함수_ 는 모든 값이 전체에 대해 정규화된 연속 히스토그램입니다. 높이를 알고 싶은 특정 나무가 있다면, 그 나무가 특정 구간에 속할 가능성이 얼마나 되는지 알려 주는 _확률 함수_ 를 만들 수 있습니다.
$$ \text{Probability of Bin i} = \frac{\text{Number of items in Bin i}}
{\text{Number of items total}} $$
_확률 함수_ 와 (연속) _밀도 함수_ 를 결합하면, 이를 나무 높이에 대한 통계적 예측 도구로 해석할 수 있습니다:
$$ \text{Probability a random tree is between } H \text{ and } H’ =
\text{Bin Density}\cdot(H-H’)$$
실제로 이 연속 확률 함수를 사용하면, 임의의 나무가 여러 구간에 걸친 임의의 범위 안에 드는 높이를 가질 가능성을 이제 말할 수 있습니다. 이것이 바로 _확률 밀도 함수_ (이하 _PDF_)입니다. 요컨대 PDF는 적분을 통해 어떤 결과가 특정 적분 구간에서 나타날 가능성이 얼마나 되는지 구할 수 있는 연속 함수입니다.
PDF 구성하기
-------------------
직관을 기르기 위해 PDF를 하나 만들어 이리저리 다루어 보겠습니다. 다음 함수를 사용하겠습니다:
![Figure [linear-pdf]: 선형 PDF](../images/fig-3.04-linear-pdf.jpg)
이 함수는 무엇을 할까요? PDF는 임의의 값 범위에 대한 가능성을 정의하는 연속 함수일 뿐이라는 것을 우리는 알고 있습니다. 이 함수 $p(r)$는 $0$과 $2$ 사이로 제한되며 그 구간을 따라 선형적으로 증가합니다. 따라서 이 함수를 PDF로 사용하여 난수를 생성한다면, 0 근처의 수를 얻을 _확률_ 은 2 근처의 수를 얻을 확률보다 낮을 것입니다.
PDF $p(r)$는 $r=0$에서 $0$으로 시작하여 $r=2$일 때 가장 높은 지점인 $p(2)$까지 단조 증가하는 선형 함수입니다. $p(2)$의 값은 무엇입니까? $p(r)$의 값은 무엇입니까? 어쩌면 $p(2)$는 2일까요? PDF는 0에서 2까지 선형적으로 증가하므로 $p(2)$의 값이 2라고 추측하는 것은 합리적으로 보입니다. 적어도 0일 수는 없어 보입니다.
PDF는 확률 함수라는 점을 기억하십시오. 우리는 PDF가 [0,2] 범위 안에 놓이도록 제한하고 있습니다. PDF는 확률적 목록에 대한 연속 밀도 함수를 나타냅니다. 그 목록의 모든 것이 0과 2 사이에 포함된다는 것을 안다면, 0과 2 사이의 값을 얻을 확률이 100%라고 말할 수 있습니다. 따라서 곡선 아래 넓이의 합은 1이어야 합니다:
$$ \operatorname{area}(p(r), 0, 2) = 1 $$
모든 선형 함수는 상수항에 변수를 곱한 형태로 나타낼 수 있습니다.
$$ p(r) = C \cdot r $$
$C$의 값을 구해야 합니다. 적분을 사용하여 거꾸로 풀어 나갈 수 있습니다.
$$ 1 = \operatorname{area}(p(r), 0, 2) $$
$$ = \int_{0}^{2} C \cdot r dr $$
$$ = C \cdot \int_{0}^{2} r dr $$
$$ = C \cdot \frac{r^2}{2} \Big|_{0}^{2} $$
$$ = C ( \frac{2^2}{2} - \frac{0}{2} ) $$
$$ C = \frac{1}{2} $$
이로써 $p(r) = r/2$라는 PDF를 얻습니다. 히스토그램에서와 마찬가지로 영역을 합산(적분)하여 $r$이 어떤 구간 $[x_0,x_1]$ 안에 있을 확률을 구할 수 있습니다:
$$ \operatorname{Probability} (r | x_0 \leq r \leq x_1 )
= \operatorname{area}(p(r), x_0, x_1)
$$
$$ \operatorname{Probability} (r | x_0 \leq r \leq x_1 ) = \int_{x_0}^{x_1} \frac{r}{2} dr $$
이해한 내용을 확인하려면 $r=0$부터 $r=2$까지의 영역에 대해 적분해 보십시오. 확률 1이 나와야 합니다.
PDF를 충분히 오래 다루다 보면 PDF를 변수 $r$이 값 $x$일 확률, _즉_ $p(r=x)$로 부르기 시작할 수도 있습니다. 그렇게 하지 마십시오. 연속 함수에서 변수가 특정한 값을 가질 확률은 항상 0입니다. PDF는 변수가 주어진 구간 안에 들어갈 확률만 알려 줄 수 있습니다. 확인하려는 구간이 단 하나의 값이라면, 그 "구간(bin)"이 무한히 얇기(폭이 0이기) 때문에 PDF는 항상 확률 0을 반환합니다. 이 사실에 대한 간단한 수학적 증명은 다음과 같습니다:
$$ \operatorname{Probability} (r = x) = \int_{x}^{x} p(r) dr $$
$$ = P(r) \Big|_{x}^{x} $$
$$ = P(x) - P(x) $$
$$ = 0 $$
x를 둘러싼 영역의 확률은 0이 아닐 수 있습니다:
$$ \operatorname{Probability} (r | x - \Delta x < r < x + \Delta x ) =
\operatorname{area}(p(r), x - \Delta x, x + \Delta x) $$
$$ = P(x + \Delta x) - P(x - \Delta x) $$
샘플 선택하기
--------------------
관심 있는 함수에 대한 PDF가 있다면, 그 함수가 임의의 구간 안의 값을 반환할 확률을 알 수 있습니다. 이를 사용하여 어디를 샘플링해야 할지 결정할 수 있습니다. 이 모든 것이 매우 밝은 픽셀 옆에 매우 어두운 픽셀이 생기지 않도록 장면을 샘플링하는 최선의 방법을 찾으려는 탐구에서 시작되었다는 점을 기억하십시오. 장면에 대한 PDF가 있다면, 이미지를 부정확하게 밝게 만들지 않으면서 샘플을 확률적으로 빛 쪽으로 유도할 수 있습니다. 샘플을 빛 쪽으로 유도하면 이미지가 부정확하게 밝아지게 _된다고_ 이미 말씀드렸습니다. 이러한 부정확성을 일으키지 않으면서 샘플을 유도하는 방법을 알아내야 하며, 이는 조금 뒤에 설명하겠습니다. 지금은 PDF가 있을 때 샘플을 생성하는 데 집중하겠습니다. PDF를 따르는 난수는 어떻게 생성할까요? 이를 위해서는 도구가 조금 더 필요합니다. 걱정하지 마십시오 -- 이 과정이 끝없이 이어지지는 않습니다!
난수 생성기 `random_double()`은 0과 1 사이의 무작위 double 값을 생성합니다. 이 난수 생성기는 0과 1 사이에서 균등하므로, 0과 1 사이의 모든 수가 나올 가능성이 같습니다. PDF가 어떤 정의역, 예를 들어 $[0,10]$에서 균일하다면, 다음과 같이 이 균일 PDF에 대한 완벽한 샘플을 간단히 만들 수 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
10.0 * random_double()
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
이는 쉬운 경우이지만, 우리가 다루게 될 경우의 대다수는 균등하지 않습니다. 균등 난수 생성기를 PDF로 분포가 정의되는 비균등 난수 생성기로 변환하는 방법을 알아내야 합니다. 균일한 입력을 받아 PDF로 가중된 불균일 분포를 만들어 내는 함수 $f(d)$가 존재한다고 그냥 _가정_ 하겠습니다. 이제 $f(d)$를 구하는 방법만 알아내면 됩니다.
위에서 주어진 PDF, 즉 $p(r) = \frac{r}{2}$에서는 무작위 샘플이 나올 확률이 0 쪽보다 2 쪽으로 갈수록 더 높습니다. 0.0과 0.2 사이의 수를 얻을 확률보다 1.8과 2.0 사이의 수를 얻을 확률이 더 큽니다. 잠시 수학자의 모자를 벗고 컴퓨터 과학자의 모자를 써 본다면, PDF를 분할하는 영리한 방법을 찾아낼 수 있을지도 모릅니다. 0 근처보다 2 근처에서 확률이 더 높다는 것은 알고 있지만, 확률을 절반으로 나누는 값은 무엇입니까? 난수가 그 값보다 클 확률이 50%이고 작을 확률이 50%인 값은 무엇입니까? 다음을 만족하는 $x$는 무엇입니까:
$$ 50\% = \int_{0}^{x} \frac{r}{2} dr = \int_{x}^{2} \frac{r}{2} dr $$
이를 풀면 다음을 얻습니다:
$$ 0.5 = \frac{r^2}{4} \Big|_{0}^{x} $$
$$ 0.5 = \frac{x^2}{4} $$
$$ x^2 = 2$$
$$ x = \sqrt{2}$$
대략적인 근사로서, `double d = random_double()`을 입력으로 받는 함수 `f(d)`를 만들 수 있습니다. `d`가 0.5보다 작거나 같으면 $[0,\sqrt{2}]$ 범위의 균일한 수를 생성하고, `d`가 0.5보다 크면 $[\sqrt{2}, 2]$ 범위의 균일한 수를 생성합니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
double f(double d)
{
if (d <= 0.5)
return std::sqrt(2.0) * random_double();
else
return std::sqrt(2.0) + (2 - std::sqrt(2.0)) * random_double();
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [crude-approx]: 불균일 PDF에 대한 대략적인 1차 근사]
$\frac{r}{2}$에 대한 새로운 대략적 근사는 (아주 약간이지만) 불균일합니다:
![Figure [approx-f]: $r/2$에 대한 불균일 분포
](../images/fig-3.06-nonuniform-dist.jpg)
위의 적분에 대해서는 해석적 해가 있었으므로, 50% 값을 매우 쉽게 구할 수 있었습니다. 하지만 이 50% 값을 실험적으로 구할 수도 있습니다. 적분을 풀 수 없거나 풀고 싶지 않은 함수도 있을 것입니다. 이러한 경우에는 실제 값에 가까운 실험적 결과를 얻을 수 있습니다. 다음 함수를 살펴보겠습니다:
$$ p(x) = e^{\frac{-x}{2 \pi}} \sin^2(x) $$
이 함수는 대략 다음과 같은 모습입니다:
![Figure [exp-sin2]: 해석적으로 풀고 싶지 않은 함수
](../images/fig-3.07-exp-sin2.jpg)
이제는 곡선 아래 넓이를 실험적으로 구하는 방법에 익숙할 것입니다. 기존 코드를 약간 수정하여 50% 값의 추정치를 구하겠습니다. 곡선 아래 전체 넓이의 절반을 주는 $x$ 값을 구하고자 합니다. N개의 샘플에 대한 누적 합을 구해 나가면서, 각 개별 샘플을 그 `p(x)` 값과 함께 저장하겠습니다. 전체 합을 구한 뒤에는 샘플을 정렬하고, 넓이가 전체의 절반이 될 때까지 더해 나가겠습니다. 예를 들어 $0$부터 $2\pi$까지:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
#include
#include
struct sample {
double x;
double p_x;
};
bool compare_by_x(const sample& a, const sample& b) {
return a.x < b.x;
}
int main() {
const unsigned int N = 10000;
sample samples[N];
double sum = 0.0;
// Iterate through all of our samples.
for (unsigned int i = 0; i < N; i++) {
// Get the area under the curve.
auto x = random_double(0, 2*pi);
auto sin_x = std::sin(x);
auto p_x = exp(-x / (2*pi)) * sin_x * sin_x;
sum += p_x;
// Store this sample.
sample this_sample = {x, p_x};
samples[i] = this_sample;
}
// Sort the samples by x.
std::sort(std::begin(samples), std::end(samples), compare_by_x);
// Find out the sample at which we have half of our area.
double half_sum = sum / 2.0;
double halfway_point = 0.0;
double accum = 0.0;
for (unsigned int i = 0; i < N; i++){
accum += samples[i].p_x;
if (accum >= half_sum) {
halfway_point = samples[i].x;
break;
}
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "Average = " << sum / N << '\n';
std::cout << "Area under curve = " << 2 * pi * sum / N << '\n';
std::cout << "Halfway = " << halfway_point << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [est-halfway]: [estimate_halfway.cc] 함수의 50% 지점 추정하기]
이 코드 조각은 이전에 작성한 코드와 크게 다르지 않습니다. 여전히 구간(0에서 $2\pi$까지)에 대한 합을 구하고 있습니다. 다만 이번에는 모든 샘플을 입력과 출력에 따라 저장하고 정렬하기도 합니다. 이를 사용하여 샘플들의 부분합이 전체 구간에 걸친 합의 절반이 되는 지점을 구합니다. 처음 $j$개의 샘플의 합이 전체 합의 절반이 된다는 것을 알면, $j\text{th}$ $x$가 대략 중간 지점에 해당한다는 것을 알 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Average = 0.314686555791
Area under curve = 1.977233943713
Halfway = 2.016002314977
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
$0$에서 $2.016$까지의 적분과 $2.016$에서 $2\pi$까지의 적분을 구하면 두 결과가 거의 정확히 같게 나올 것입니다.
이제 PDF를 절반으로 나누는 중간 지점을 구하는 방법이 생겼습니다. 원한다면 이를 사용하여 PDF의 중첩된 이진 분할을 만들 수 있습니다:
1. PDF의 중간 지점을 구합니다
2. 아래쪽 절반으로 재귀하여 1단계를 반복합니다
3. 위쪽 절반으로 재귀하여 1단계를 반복합니다
적당한 깊이, 예를 들어 6–10 정도에서 멈춥니다. 짐작하시듯이 이는 계산 비용이 상당히 클 수 있습니다. 위 코드의 계산 병목은 아마도 샘플을 정렬하는 부분일 것입니다. 단순한 정렬 알고리즘은 $\mathcal{O}(\mathbf{n^2})$ 시간의 알고리즘 복잡도를 가질 수 있으며, 이는 엄청나게 비용이 큽니다. 다행히 표준 라이브러리에 포함된 정렬 알고리즘은 대개 $\mathcal{O}(\mathbf{n\log{}n})$ 시간에 훨씬 가깝지만, 특히 수백만 또는 수십억 개의 샘플에 대해서는 여전히 상당히 비용이 클 수 있습니다. 하지만 이 방법은 불균일한 수의 괜찮은 불균일 분포를 만들어 냅니다. 불균일 분포를 만드는 이 분할 정복 방법은 실제로 꽤 흔히 사용되지만, 단순한 이진 분할보다 훨씬 더 효율적인 방법들도 있습니다. 분포의 PDF로 사용하고 싶은 임의의 함수가 있다면 _메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)_ 을 조사해 보십시오.
분포 근사하기
---------------------------
몇 가지 개념을 쌓아 올리기 위해 많은 수학과 작업이 필요했습니다. 처음의 PDF로 돌아가 보겠습니다. 명시적인 확률이 없는 구간에서는 PDF가 0이라고 가정합니다. 따라서 이 장의 처음에 나온 예제에서는 $r \notin [0,2]$에 대해 $p(r) = 0$입니다. $p(r)$을 구간별 형태로 다시 쓸 수 있습니다:
$$ p(r)=\begin{cases}
0 & r < 0 \\
\frac{r}{2} & 0 \leq r \leq 2 \\
0 & 2 < r \\
\end{cases}
$$
이전 절에서 하려던 것을 생각해 보면, 많은 수학이 0부터의 _누적_ 넓이(또는 _누적_ 확률)를 중심으로 이루어졌습니다. 다음 함수의 경우
$$ f(x) = e^{\frac{-x}{2 \pi}} \sin^2(x) $$
$0$에서 $2\pi$까지의 누적 확률(100%)과 $0$에서 $2.016$까지의 누적 확률(50%)에 관심이 있었습니다. 이를 중요한 용어로 일반화할 수 있는데, _누적 분포 함수(Cumulative Distribution Function)_ $P(x)$는 다음과 같이 정의됩니다:
$$ P(x) = \int_{-\infty}^{x} p(x') dx' $$
또는,
$$ P(x) = \operatorname{area}(p(x'), -\infty, x) $$
이는 $-\infty$부터의 _누적_ 확률의 양입니다. 미적분 규칙 때문에 적분을 $x$ 대신 $x'$에 대해 다시 썼습니다. 무슨 의미인지 잘 모르겠더라도 걱정하지 마십시오. 같은 것으로 취급해도 됩니다. 위에서 설명한 적분을 수행하면 구간별로 정의된 $P(r)$을 얻습니다:
$$ P(r)=\begin{cases}
0 & r < 0 \\
\frac{r^2}{4} & 0 \leq r \leq 2 \\
1 & 2 < r \\
\end{cases}
$$
_확률 밀도 함수_ (PDF)는 어떤 수의 구간이 선택될 가능성이 얼마나 되는지를 설명하는 확률 함수입니다. _누적 분포 함수_ (CDF)는 입력값보다 작은 모든 수가 선택될 가능성이 얼마나 되는지를 설명하는 분포 함수입니다. PDF에서 CDF로 가려면 $-\infty$부터 $x$까지 적분해야 하지만, CDF에서 PDF로 가려면 도함수를 구하기만 하면 됩니다:
$$ p(x) = \frac{d}{dx}P(x) $$
CDF $P(r)$을 $r = 1.0$에서 계산하면 다음을 얻습니다:
$$ P(1.0) = \frac{1}{4} $$
이는 _우리의 PDF에서 뽑은 확률 변수가 1 이하일 확률이 25%_ 라는 뜻입니다. 우리는 0과 1 사이의 균일 분포(_즉_ `f(random_double())`)를 입력으로 받아, CDF가 $P(x) = \frac{x^2}{4}$인 분포에 따라 무작위 값을 반환하는 함수 $f(d)$를 원합니다. 함수 $f(d)$가 해석적으로 무엇인지는 아직 모르지만, 반환하는 값의 25%는 1.0보다 작아야 하고 75%는 1.0보다 커야 한다는 것은 알고 있습니다. 마찬가지로, 반환하는 값의 50%는 $\sqrt{2}$보다 작아야 하고 50%는 $\sqrt{2}$보다 커야 한다는 것도 알고 있습니다. $f(d)$가 단조 증가한다면 $f(0.25) = 1.0$이고 $f(0.5) = \sqrt{2}$일 것으로 예상할 수 있습니다. 이를 일반화하면 가능한 모든 입력에 대해 $f(d)$를 구할 수 있습니다:
$$ f(P(x)) = x $$
샘플을 몇 개 더 구해 보겠습니다:
$$ P(0.0) = 0 $$
$$ P(0.5) = \frac{1}{16} $$
$$ P(1.0) = \frac{1}{4} $$
$$ P(1.5) = \frac{9}{16} $$
$$ P(2.0) = 1 $$
따라서 함수 $f()$는 다음 값을 가집니다
$$ f(P(0.0)) = f(0) = 0 $$
$$ f(P(0.5)) = f(\frac{1}{16}) = 0.5 $$
$$ f(P(1.0)) = f(\frac{1}{4}) = 1.0 $$
$$ f(P(1.5)) = f(\frac{9}{16}) = 1.5 $$
$$ f(P(2.0)) = f(1) = 2.0 $$
이 중간값들을 사용하고 그 사이를 보간하여 $f(d)$를 근사할 수 있습니다:
![Figure [approx f]: 비균일 f()의 근사](../images/fig-3.08-approx-f.jpg)
PDF를 해석적으로 풀 수 없다면 CDF도 해석적으로 풀 수 없습니다. 결국 CDF는 PDF의 적분일 뿐이기 때문입니다. 하지만 PDF를 근사하는 분포는 여전히 만들 수 있습니다. PDF를 구하고자 하는 무작위 함수에서 샘플을 많이 추출하면, 그 샘플들의 히스토그램을 구한 뒤 PDF로 변환하여 PDF를 근사할 수 있습니다. 또는 위에서 했던 것처럼 모든 샘플을 정렬할 수도 있습니다.
등식을 더 자세히 살펴보겠습니다:
$$ f(P(x)) = x $$
이는 $f()$가 $P()$가 하는 일을 그대로 되돌린다는 뜻일 뿐입니다. 따라서 $f()$는 역함수입니다:
$$ f(d) = P^{-1}(x) $$
우리의 목적상, PDF $p()$와 누적 분포 함수 $P()$가 있다면 이 "역함수"를 난수와 함께 사용하여 원하는 결과를 얻을 수 있습니다:
$$ f(d) = P^{-1} (\operatorname{random\_double}()) $$
우리의 PDF $p(r) = r/2$와 이에 대응하는 $P(r)$에 대해, $P(r)$의 역함수를 계산해야 합니다. 다음과 같다면
$$ y = \frac{r^2}{4} $$
$r$을 $y$에 대해 풀어 역함수를 얻습니다:
$$ r = \sqrt{4y} $$
즉, CDF의 역함수($ICD(x)$라고 부르겠습니다)는 다음과 같이 정의됩니다
$$ P^{-1}(r) = \operatorname{ICD}(r) = \sqrt{4y} $$
따라서 밀도가 $p(r)$인 난수 생성기는 다음과 같이 만들 수 있습니다:
$$ \operatorname{ICD}(d) = \sqrt{4 \cdot \operatorname{random\_double}()} $$
이 값은 기대한 대로 0부터 2까지의 범위를 가진다는 점에 주목하십시오. 검증해 보면 `random_double()`을 $1/4$로 바꾸면 1을 얻고, $1/2$로 바꾸면 $\sqrt{2}$를 얻어 예상과 정확히 일치합니다.
중요도 샘플링
--------------------
이제 여러분은 해석적인 PDF를 가지고 그 분포를 따르는 난수를 생성하는 함수를 만드는 방법을 충분히 이해했을 것입니다. 더 잘 이해하기 위해 원래의 적분으로 돌아가서 몇 가지 다른 PDF로 시도해 보겠습니다:
$$ I = \int_{0}^{2} x^2 dx $$
지난번에 적분을 풀려고 했을 때는 구간 $[0, 2]$에서 균일하게 샘플링하는 몬테카를로 방식을 사용했습니다. 당시에는 몰랐지만, 우리는 암묵적으로 0과 2 사이의 균일 PDF를 사용하고 있었습니다. 이는 범위 $[0,2]$에서 PDF = $1/2$를 사용하고 있다는 뜻이며, 따라서 CDF는 $P(x) = x/2$이고 $\operatorname{ICD}(d) = 2d$입니다. 이를 알았으니, 이 균일 PDF를 명시적으로 만들 수 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double icd(double d) {
return 2.0 * d;
}
double pdf(double x) {
return 0.5;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ delete
int a = 0;
int b = 2;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto x = icd(random_double());
sum += x*x / pdf(x);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
std::cout << std::fixed << std::setprecision(12);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
std::cout << "I = " << (sum / N) << '\n';
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-xsq-2]: [integrate_x_sq.cc] $x^2$에 대한 명시적 균일 PDF]
강조해야 할 중요한 점이 몇 가지 있습니다. $x$의 각 값은 분포 $[0, 2]$ 안에서 함수 $x^2$의 샘플 하나를 나타냅니다. 이 분포 안에서 샘플을 무작위로 선택하기 위해 함수 $\operatorname{ICD}$를 사용합니다. 이전에는 최종 답을 얻기 위해 구간에 대한 평균(`sum / N`)에 구간의 길이(`b - a`)를 곱했습니다. 여기서는 구간 길이를 곱할 필요가 없습니다. 즉, 더 이상 평균에 $2$를 곱할 필요가 없습니다.
$x$의 PDF가 균일하지 않다는 점을 고려해야 합니다. 이 비균일성을 고려하지 않으면 장면에 편향이 생깁니다. 실제로 이 편향이 부정확하게 밝은 이미지의 원인입니다. 비균일성을 고려하면 정확한 결과를 얻을 수 있습니다. PDF는 샘플을 분포의 특정 부분으로 "유도"하며, 이는 더 빠르게 수렴하게 해 주지만 편향이 생긴다는 대가가 따릅니다. 이 편향을 제거하려면 더 자주 샘플링하는 곳에는 가중치를 낮추고, 덜 자주 샘플링하는 곳에는 가중치를 높여야 합니다. 새로운 비균일 난수 생성기에서는 PDF가 특정 부분을 얼마나 많이 또는 적게 샘플링할지를 정의합니다. 따라서 가중치 함수는 $1/\mathit{pdf}$에 비례해야 합니다. 사실 가중치 함수는 _정확히_ $1/\mathit{pdf}$입니다. 이것이 `x*x`를 `pdf(x)`로 나누는 이유입니다.
선형 PDF $p(r) = \frac{r}{2}$를 사용하여 적분을 풀어 볼 수 있으며, 이 PDF에 대해서는 CDF와 그 역함수인 ICD를 구할 수 있었습니다. 그렇게 하려면 함수를 $\operatorname{ICD}(d) = \sqrt{4d}$와 $p(x) = x/2$로 바꾸기만 하면 됩니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double icd(double d) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
return std::sqrt(4.0 * d);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
double pdf(double x) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
return x / 2.0;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
int main() {
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto z = random_double();
if (z == 0.0) // Ignore zero to avoid NaNs
continue;
auto x = icd(z);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
sum += x*x / pdf(x);
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "I = " << sum / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-xsq-3]: [integrate_x_sq.cc] 선형 PDF로 $x^2$ 적분하기]
균일 PDF와 선형 PDF의 실행 결과를 비교해 보았다면, 아마 선형 PDF가 더 빠르게 수렴한다는 것을 발견했을 것입니다. 생각해 보면, 선형 PDF는 균일 PDF보다 이차 함수를 더 잘 근사할 가능성이 높으므로 더 빠르게 수렴할 것으로 예상할 수 있습니다. 그렇다면 PDF를 이차 함수로 바꾸어 PDF가 피적분 함수와 일치하도록 만들어 보아야 합니다:
$$ p(r)=\begin{cases}
0 & r < 0 \\
C \cdot r^2 & 0 \leq r \leq 2 \\
0 & 2 < r \\
\end{cases}
$$
선형 PDF와 마찬가지로, 구간에 대해 적분한 값이 1이 되도록 하여 상수 $C$를 구하겠습니다:
$$ 1 = \int_{0}^{2} C \cdot r^2 dr $$
$$ = C \cdot \int_{0}^{2} r^2 dr $$
$$ = C \cdot \frac{r^3}{3} \Big|_{0}^{2} $$
$$ = C ( \frac{2^3}{3} - \frac{0}{3} ) $$
$$ C = \frac{3}{8} $$
그러면 다음을 얻습니다:
$$ p(r)=\begin{cases}
0 & r < 0 \\
\frac{3}{8} r^2 & 0 \leq r \leq 2 \\
0 & 2 < r \\
\end{cases}
$$
그리고 이에 대응하는 CDF를 얻습니다:
$$ P(r) = \frac{r^3}{8} $$
그리고
$$ P^{-1}(x) = \operatorname{ICD}(d) = 8d^\frac{1}{3} $$
샘플을 하나만 사용하면 다음과 같은 결과를 얻습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double icd(double d) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
return 8.0 * std::pow(d, 1.0/3.0);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
double pdf(double x) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
return (3.0/8.0) * x*x;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
int main() {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
int N = 1;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
auto sum = 0.0;
for (int i = 0; i < N; i++) {
auto z = random_double();
if (z == 0.0) // Ignore zero to avoid NaNs
continue;
auto x = icd(z);
sum += x*x / pdf(x);
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "I = " << sum / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [integ-xsq-5]: [integrate_x_sq.cc] $x^2$ 적분하기, 최종 버전]
이 방법은 항상 정확한 답을 반환합니다.
균일하지 않은 PDF는 PDF가 큰 곳으로는 더 많은 샘플을, PDF가 작은 곳으로는 더 적은 샘플을 "유도"합니다. 이러한 샘플링을 통해 PDF가 큰 곳에서는 노이즈가 적고 PDF가 작은 곳에서는 노이즈가 많을 것으로 예상할 수 있습니다. 장면에서 노이즈가 많은 부분에서는 더 높고 노이즈가 적은 부분에서는 더 낮은 PDF를 선택하면, 더 적은 샘플로 장면 전체의 노이즈를 줄일 수 있습니다. 이는 균일 PDF를 사용할 때보다 _더 빠르게_ 올바른 장면으로 수렴할 수 있다는 뜻입니다. 사실상 우리는 분포에서 더 _중요한_ 부분으로 샘플을 유도하고 있는 것입니다. 이것이 신중하게 선택한 균일하지 않은 PDF를 사용하는 방법을 흔히 _중요도 샘플링_ 이라고 부르는 이유입니다.
앞서 제시한 모든 예제에서 우리는 항상 올바른 답인 $8/3$ 으로 수렴했습니다. 균일 PDF와 "올바른" PDF(즉, $\operatorname{ICD}(d) = 8d^{\frac{1}{3}}$)를 모두 사용했을 때 같은 답을 얻었습니다. 둘 다 같은 답으로 수렴했지만, 균일 PDF는 훨씬 오래 걸렸습니다. 결국 적분과 완벽하게 일치하는 PDF에서는 샘플 하나만 있으면 충분했습니다. 이는 납득할 만한 결과입니다. 우리는 분포의 중요한 부분을 더 자주 샘플링하도록 선택한 반면, 균일 PDF는 중요도를 고려하지 않고 전체 분포를 동일하게 샘플링했기 때문입니다.
실제로 여러분이 만드는 어떤 PDF든 마찬가지입니다--결국에는 모두 수렴합니다. 이것 또한 몬테카를로 알고리즘이 지닌 강력함의 일부입니다. 50% 값을 구해 분포를 $[0, \sqrt{2}]$ 와 $[\sqrt{2}, 2]$ 의 두 부분으로 나눈 단순한 PDF조차도 마찬가지입니다. 그 PDF도 수렴합니다. 그 PDF가 순수한 균일 PDF보다는 빠르게, 선형 PDF(즉, $\operatorname{ICD}(d) = \sqrt{4d}$)보다는 느리게 수렴하는 이유를 직관적으로 이해하셨기를 바랍니다.
완벽한 중요도 샘플링은 답을 이미 알고 있을 때만 가능하지만($p$ 를 해석적으로 적분하여 $P$ 를 얻었습니다), 코드가 제대로 동작하는지 확인하기에 좋은 연습입니다.
몬테카를로 레이 트레이서의 바탕이 되는 주요 개념을 복습해 보겠습니다:
1. 어떤 정의역 $[a,b]$ 에 대한 $f(x)$ 의 적분이 있습니다
2. $[a,b]$ 에서 0이 아니고 음이 아닌 PDF $p$ 를 선택합니다
3. $r$ 이 PDF $p$ 를 따르는 난수일 때, 아주 많은 $\frac{f(r)}{p(r)}$ 값의 평균을 구합니다.
어떤 PDF $p$ 를 선택하든 항상 올바른 답으로 수렴하지만, $p$ 가 $f$ 를 더 가깝게 근사할수록 더 빠르게 수렴합니다.
방향 구면 위에서의 몬테카를로 적분
====================================================================================================
1차원 몬테카를로 적분 장에서 우리는 균등 난수에서 시작하여, 한 장에 걸쳐 점차 더 복잡한 난수 생성 방법을 쌓아 올렸고, 최종적으로 PDF에 대한 직관과 이를 사용하여 임의의 분포를 갖는 난수를 생성하는 방법에 도달했습니다.
그 장에서 다룬 모든 개념은 단일 차원을 넘어 확장하더라도 계속 유효합니다. 앞으로는 2차원, 3차원 또는 그보다 더 높은 차원의 공간에서 점을 선택한 다음, 그 선택에 임의의 PDF로 가중치를 부여해야 할 수도 있습니다. 이것의 중요한 사례 중 하나는--적어도 레이 트레이싱에서는--무작위 방향을 생성하는 것입니다. 처음 두 권의 책에서는 무작위 벡터를 만들고 그것이 단위 구 바깥에 있으면 기각하는 방식으로 무작위 방향을 생성했습니다. 단위 구 안에 있는 무작위 벡터를 찾을 때까지 이 과정을 반복했습니다. 이 벡터를 정규화하면 단위 구 위에 정확히 놓이는 점이 만들어지며, 이 점이 곧 무작위 방향을 나타냅니다. 샘플을 생성하고 원하는 공간 안에 있지 않으면 기각하는 이 과정을 _기각법_ 이라고 하며, 문헌 곳곳에서 찾아볼 수 있습니다. 앞 장에서 다룬 방법은 PDF를 역변환하기 때문에 _역변환법_ 이라고 합니다.
3D 공간의 모든 방향에는 단위 구 위에 대응하는 점이 있으며, 원점에서 그 대응점으로 향하는 벡터를 구하면 해당 방향을 생성할 수 있습니다. 무작위 방향을 선택하는 것은 제한된 2차원 평면, 즉 단위 구를 데카르트 좌표계로 매핑하여 만든 평면에서 무작위 점을 선택하는 것으로 생각할 수 있습니다. 이전과 같은 방법론이 적용되지만, 이제는 2차원에 걸쳐 정의된 PDF를 다루게 될 수 있습니다. 단위 구의 표면에 걸쳐 다음 함수를 적분하려 한다고 가정해 보겠습니다:
$$ f(\theta, \phi) = \cos^2(\theta) $$
몬테카를로 적분을 사용하면 $\cos^2(\theta) / p(r)$를 샘플링하기만 하면 되며, 이때 $p(r)$는 이제 방향에 대한 밀도인 $p(direction)$일 뿐입니다. 그런데 이 맥락에서 _방향_ 이란 무엇입니까? 극좌표를 기반으로 할 수 있으며, 그러면 $p$는 $p(\theta, \phi)$에서처럼 $\theta$와 $\phi$로 표현됩니다. 어떤 좌표계를 사용할지는 중요하지 않습니다. 다만 어떤 방식을 택하든, PDF는 전체 표면에 걸쳐 적분했을 때 1이 되어야 하며 PDF는 해당 방향이 샘플링될 _상대적 확률_ 을 나타낸다는 점을 기억하십시오. 단위 구 $d$ 위에서 균일한 무작위 샘플을 생성하는 `vec3` 함수(`random_unit_vector()`)가 있다는 것을 떠올려 보십시오. 이 균일 샘플의 PDF는 무엇입니까? 단위 구 위의 균일 밀도이므로 구의 $1/\mathit{area}$이며, 이는 $1/(4\pi)$입니다. 피적분 함수가 $\cos^2(\theta)$이고 $\theta$가 $z$ 축과 이루는 각도라면, 스칼라 사영을 사용하여 $\cos^2(\theta)$를 $d_z$로 다시 쓸 수 있습니다:
$$ d_z = \lVert d \rVert \cos \theta = 1 \cdot \cos \theta $$
그런 다음 $1 \cdot \cos \theta$를 $d_z$로 치환하면 다음을 얻습니다:
$$ f(\theta, \phi) = \cos^2 (\theta) = {d_z}^2 $$
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
double f(const vec3& d) {
auto cosine_squared = d.z()*d.z();
return cosine_squared;
}
double pdf(const vec3& d) {
return 1 / (4*pi);
}
int main() {
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
vec3 d = random_unit_vector();
auto f_d = f(d);
sum += f_d / pdf(d);
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "I = " << sum / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [main-sphereimp]: [sphere_importance.cc] 단위 구 위에서 중요도 샘플링된 점 생성하기]
해석적 답은 $\frac{4}{3} \pi = 4.188790204786391$입니다 -- 고급 미적분을 충분히 기억하신다면 직접 확인해 보십시오! 그리고 위의 코드가 바로 그 값을 산출합니다. 여기서 핵심은 모든 적분과 확률, 그리고 그 밖의 모든 것이 단위 구 위에서 정의된다는 점입니다. 3D에서 하나의 방향을 표현하는 방법은 단위 구 위의 대응하는 점을 사용하는 것입니다. 3D에서 방향의 범위를 표현하는 방법은 그 방향들이 통과하는 단위 구 위의 면적을 사용하는 것입니다. 방향이라고 하든, 면적이라고 하든, _입체각_ 이라고 하든 모두 같은 것입니다. 입체각은 문헌에서 흔히 볼 수 있는 용어입니다. 1차원에서는 $\theta$에 대한 라디안(r)이 있고, 2차원에서는 $\theta$ _와_ $\phi$에 대한 _스테라디안_ (sr)이 있습니다(단위 구는 3차원 물체이지만 그 표면은 2차원일 뿐입니다). 입체각은 단지 각도를 2차원으로 확장한 것입니다. 2차원 각도에 익숙하시다면 좋습니다! 그렇지 않다면 제가 하는 것처럼 일련의 방향들이 통과하는 단위 구 위의 면적을 떠올려 보십시오. 입체각 $\omega$와 단위 구 위에 투영된 면적 $A$는 같은 것입니다.
![Figure [solid-angle]: 구의 입체각 / 투영 면적
](../images/fig-3.09-solid-angle.jpg)
이제 우리가 풀고 있는 빛 전달 방정식으로 넘어가겠습니다.
빛 산란
====================================================================================================
이 장에서는 실제로 아무것도 프로그래밍하지 않습니다. 다음 장에서 있을 조명의 큰 변화를 준비하기만 할 것입니다. 앞의 두 책에서 만든 레이 트레이싱 프로그램은 광선이 표면이나 볼륨과 상호작용할 때 광선을 산란시킵니다. 광선 산란은 장면을 통한 빛의 전파를 시뮬레이션하는 데 가장 흔히 사용되는 모델입니다. 이는 자연스럽게 확률적으로 모델링할 수 있습니다. 광선의 확률적 산란을 모델링할 때는 고려해야 할 사항이 많습니다.
알베도
-------
먼저, 빛이 흡수됩니까?
빛이 산란될 확률: $A$
빛이 흡수될 확률: $1-A$
여기서 $A$는 _알베도_ 를 나타냅니다. 첫 번째 책에서 다루었듯이, 알베도는 부분 반사율의 한 형태라는 점을 기억하십시오. 빛의 전파를 시뮬레이션할 때 우리가 하는 일은 공간을 통과하는 광자의 움직임을 시뮬레이션하는 것뿐이라는 점을 잠시 멈추고 떠올려 보면 도움이 될 수 있습니다. 고등학교 물리를 기억하신다면 모든 광자는 플랑크 상수로 연관된 고유한 에너지와 파장을 가진다는 것을 떠올리실 수 있을 것입니다:
$$ E = \frac{hc}{\lambda} $$
개별 광자 하나는 _아주 작은_ 양의 에너지를 가지지만, 충분히 많이 합치면 렌더링의 모든 조명이 됩니다. 표면이나 볼륨(또는 실제로 광자가 상호작용할 수 있는 모든 것)에서 광자가 흡수되는지 산란되는지는 물체의 알베도에 의해 확률적으로 결정됩니다. 어떤 물체는 특정 파장을 더 잘 흡수하기 때문에 알베도는 색상에 따라 달라질 수 있습니다.
대부분의 물리 기반 렌더러에서는 빛의 색상을 표현할 때 RGB 대신 미리 정의된 특정 파장들의 집합을 사용합니다. 예를 들어, _삼자극_ RGB 렌더러를 300nm, 350nm, 400nm, ..., 700nm에서 명시적으로 샘플링하는 방식으로 대체할 수 있습니다. R, G, B를 파장들의 특정한 대수적 혼합으로 생각하면 직관을 확장할 수 있습니다. 여기서 R은 _대부분_ 빨간색 파장이고, G는 _대부분_ 초록색 파장이며, B는 _대부분_ 파란색 파장입니다. 이는 인간 시각 체계를 근사한 것입니다. 인간 시각 체계에는 _원뿔세포_ 라고 불리는 3가지 고유한 색 수용체 집합이 있으며, 각각은 서로 다른 파장의 대수적 혼합에 민감합니다. 이들은 대략 RGB에 해당하지만 장파장, 중파장, 단파장 원뿔세포라고 불립니다(이 이름은 각 원뿔세포가 민감하게 반응하는 파장을 가리키는 것이지, 원뿔세포의 길이를 가리키는 것이 아닙니다). 색을 RGB 색 공간에서의 세기로 표현할 수 있는 것처럼, _LMS 색 공간_ (장파장, 중파장, 단파장)에서 각 원뿔세포 집합이 얼마나 자극되는지로도 색을 표현할 수 있습니다.
산란
-----------
빛이 산란된다면, 입체각에 대한 PDF로 기술할 수 있는 방향 분포를 가지게 됩니다. 이것을 _산란 PDF_ 라고 부르겠습니다: $\operatorname{pScatter}()$. 산란 PDF는 나가는 방향에 따라 달라집니다: $\operatorname{pScatter}(\omega_o)$. 산란 PDF는 _입사 방향_ 에 따라서도 달라질 수 있습니다: $\operatorname{pScatter}(\omega_i, \omega_o)$. 도로에서 반사되는 모습을 보면 입사 방향에 따라 이것이 달라지는 것을 확인할 수 있습니다. 시야각(입사각)이 스치는 각도에 가까워질수록 거울처럼 됩니다. 산란 PDF는 빛의 파장에 따라서도 달라질 수 있습니다: $\operatorname{pScatter}(\omega_i, \omega_o, \lambda)$. 그 좋은 예는 백색광을 굴절시켜 무지개로 만드는 프리즘입니다. 마지막으로, 산란 PDF는 산란 위치에 따라서도 달라질 수 있습니다: $\operatorname{pScatter}(\mathbf{x}, \omega_i, \omega_o, \lambda)$. $\mathbf{x}$는 산란 위치를 나타내는 수학 표기일 뿐입니다: $\mathbf{x} = (x, y, z)$. 물체의 알베도 역시 이러한 양에 따라 달라질 수 있습니다: $A(\mathbf{x}, \omega_i, \omega_o, \lambda)$.
표면의 색은 이 항들을 입사 방향에 대해 단위 반구 위에서 적분하여 구합니다:
$$ \operatorname{Color}_o(\mathbf{x}, \omega_o, \lambda) = \int_{\omega_i}
A(\mathbf{x}, \omega_i, \omega_o, \lambda) \cdot
\operatorname{pScatter}(\mathbf{x}, \omega_i, \omega_o, \lambda) \cdot
\operatorname{Color}_i(\mathbf{x}, \omega_i, \lambda) $$
$\operatorname{Color}_i$ 항을 추가했습니다. 물체 표면에서의 산란 PDF와 알베도는 그 지점을 비추는 빛에 대한 필터 역할을 합니다. 따라서 그 지점을 비추는 빛을 구해야 합니다. 이것은 재귀 알고리즘이며, `ray_color` 함수가 현재 물체의 색에 다음 광선의 색을 곱한 값을 반환하는 이유입니다.
산란 PDF
-------------------
몬테카를로 기본 공식을 적용하면 다음과 같은 통계적 추정값을 얻습니다:
$$ \operatorname{Color}_o(\mathbf{x}, \omega_o, \lambda) \approx \sum
\frac{A(\, \ldots \,) \cdot
\operatorname{pScatter}(\, \ldots \,) \cdot
\operatorname{Color}_i(\, \ldots \,)}
{p(\mathbf{x}, \omega_i, \omega_o, \lambda)} $$
여기서 $p(\mathbf{x}, \omega_i, \omega_o, \lambda)$는 무작위로 생성하는 나가는 방향의 PDF입니다.
램버시안 표면의 경우, $pScatter(\, \ldots \,)$가 코사인 밀도인 특수한 경우에 대해 이미 이 공식을 암묵적으로 구현했습니다. 램버시안 표면의 $\operatorname{pScatter}(\, \ldots \,)$는 $\cos(\theta_o)$에 비례하며, 여기서 $\theta_o$는 표면 법선에 대한 상대 각도입니다($\theta_o \in [0,\pi]$). 각도 $0$은 나가는 방향이 표면 법선과 같은 방향임을 나타내고, 각도 $\pi$는 나가는 방향이 법선 벡터와 정확히 반대임을 나타냅니다.
$C$를 한 번 더 구해 보겠습니다:
$$ \operatorname{pScatter}(\mathbf{x}, \omega_i, \omega_o, \lambda) = C \cdot \cos(\theta_o) $$
모든 2차원 PDF는 전체 표면에 대해 적분했을 때 1이 되어야 합니다($\operatorname{pScatter}$가 PDF라는 점을 기억하십시오). 지평선 아래로 산란되지 않도록 $\operatorname{pScatter}(\frac{\pi}{2} < \theta_o \le \pi) = 0$으로 설정합니다. 이를 고려하면 $\theta \in [0, \frac{\pi}{2}]$ 구간만 적분하면 됩니다.
$$ 1 = \int_{\phi = 0}^{2 \pi} \int_{\theta = 0}^\frac{\pi}{2} C \cdot \cos(\theta) dA $$
반구 위에서 적분하려면 구면 좌표계에서 다음이 성립한다는 것을 기억하십시오:
$$ dA = \sin(\theta) d\theta d\phi $$
따라서:
$$ 1 = C \cdot \int_0^{2 \pi} \int_0^\frac{\pi}{2}
\cos(\theta) \sin(\theta) d\theta d\phi $$
$$ 1 = C \cdot 2 \pi \frac{1}{2} $$
$$ 1 = C \cdot \pi $$
$$ C = \frac{1}{\pi} $$
반구에 대한 $\cos(\theta_o)$의 적분은 $\pi$이므로, $\frac{1}{\pi}$로 정규화해야 합니다. PDF $\operatorname{pScatter}$는 나가는 방향($\omega_o$)에만 의존하므로, 그 표현을 $\operatorname{pScatter}(\omega_o)$로 단순화하겠습니다. 이 모든 것을 종합하면 램버시안 표면에 대한 산란 PDF를 얻습니다:
$$ \operatorname{pScatter}(\omega_o) = \frac{\cos(\theta_o)}{\pi} $$
$p(\mathbf{x}, \omega_i, \omega_o, \lambda)$가 산란 PDF와 같다고 가정하겠습니다:
$$ p(\omega_o) = \operatorname{pScatter}(\omega_o) = \frac{\cos(\theta_o)}{\pi} $$
분자와 분모가 서로 상쇄되어 다음을 얻습니다:
$$ \operatorname{Color}_o(\mathbf{x}, \omega_o, \lambda) \approx \sum
A(\, \ldots \,) \cdot
\operatorname{Color}_i(\, \ldots \,) $$
이것은 원래의 `ray_color()` 함수에 있던 것과 정확히 같습니다!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
return attenuation * ray_color(scattered, depth-1, world);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
위의 설명은 표면과 볼륨에 동일한 수학이 적용되기를 원하기 때문에 다소 비표준적입니다. 문헌을 읽어 보면, 반사가 _양방향 반사율 분포 함수(Bidirectional Reflectance Distribution Function)_ (BRDF)로 정의되는 것을 보게 될 것입니다. 이는 우리의 용어와 꽤 간단하게 연결됩니다:
$$ BRDF(\omega_i, \omega_o, \lambda) = \frac{A(\mathbf{x}, \omega_i, \omega_o, \lambda) \cdot
\operatorname{pScatter}(\mathbf{x}, \omega_i, \omega_o, \lambda)}{\cos(\theta_o)} $$
예를 들어 램버시안 표면의 경우 $BRDF = A / \pi$입니다. 우리의 용어와 BRDF 사이의 변환은 쉽습니다. 참여 매질(볼륨)의 경우, 우리의 알베도는 보통 _산란 알베도_ 라고 불리며, 우리의 산란 PDF는 보통 _위상 함수_ 라고 불립니다.
여기서 한 것은 재질의 램버시안 산란에 대한 PDF의 개요를 설명한 것뿐입니다. 하지만 광원 쪽과 같은 중요한 방향으로 추가 광선을 보낼 수 있도록 일반화해야 합니다.
중요도 샘플링 다루기
====================================================================================================
앞으로 몇 장에 걸친 목표는 이미지의 노이즈가 줄어들도록 광원 쪽으로 많은 추가 광선을 보내게끔 프로그램을 구성하는 것입니다. PDF $\operatorname{pLight}(\omega_o)$를 사용하여 광원 쪽으로 많은 광선을 보낼 수 있다고 가정하겠습니다. 또한 $\operatorname{pScatter}$와 관련된 PDF가 있다고 가정하고, 이를 $\operatorname{pSurface}(\omega_o)$라고 부르겠습니다. PDF의 훌륭한 점은 이들의 선형 혼합을 사용하여 역시 PDF인 혼합 밀도를 만들 수 있다는 것입니다. 예를 들어, 가장 간단한 것은 다음과 같습니다:
$$ p(\omega_o) = \frac{1}{2} \operatorname{pSurface}(\omega_o) + \frac{1}{2}
\operatorname{pLight}(\omega_o)$$
가중치가 양수이고 합이 1이기만 하면, 이러한 PDF의 혼합은 모두 PDF입니다. 기억하십시오, 어떤 PDF든 사용할 수 있습니다: _모든 PDF는 결국 올바른 답으로 수렴합니다_. 따라서 관건은 다음 곱이 큰 곳에서 PDF를 더 크게 만드는 방법을 알아내는 것입니다.
$$ \operatorname{pScatter}(\mathbf{x}, \omega_i, \omega_o) \cdot
\operatorname{Color}_i(\mathbf{x}, \omega_i) $$
가장 큰 곳에서 PDF를 더 크게 만드는 방법을 알아내는 것입니다. 확산 표면의 경우, 이는 주로 $\operatorname{Color}_i(\mathbf{x}, \omega_i)$가 가장 큰 곳을 추측하는 문제입니다. 이는 가장 많은 빛이 어디에서 오는지 추측하는 것과 같습니다.
거울의 경우, $\operatorname{pScatter}()$는 한 방향 근처에서만 매우 크므로, $\operatorname{pScatter}()$가 훨씬 더 중요합니다. 실제로 대부분의 렌더러는 거울을 특수한 경우로 처리하고 $\operatorname{pScatter}()/p()$를 암묵적으로 만듭니다 -- 현재 우리 코드도 그렇게 합니다.
코넬 박스로 돌아가기
-----------------------------
코넬 박스의 몇 가지 매개변수를 조정해 보겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cam.samples_per_pixel = 1000;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [cornell-box]: [main.cc] 리팩터링된 코넬 박스]
600×600 해상도에서 제 코드는 Macbook의 코어 1개로 15분 만에 이 이미지를 생성합니다:
이 노이즈를 줄이는 것이 우리의 목표입니다. 광원 쪽으로 더 많은 광선을 보내는 PDF를 구성하여 이를 달성할 것입니다.
먼저, 코드가 어떤 PDF를 명시적으로 샘플링한 다음 그에 맞게 정규화하도록 코드를 계측해 보겠습니다. 몬테카를로의 기본을 기억하십시오: $\int f(x) \approx \sum f(r)/p(r)$. 램버시안 재질에 대해서는 지금 하는 방식대로 샘플링하겠습니다: $p(\omega_o) = \cos(\theta_o) / \pi$.
이 중요도 샘플링을 가능하게 하도록 기반 클래스 `material`을 수정하고, 램버시안 재질에 대한 산란 PDF 함수를 정의합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
virtual double scattering_pdf(const ray& r_in, const hit_record& rec, const ray& scattered)
const {
return 0;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
};
class lambertian : public material {
public:
lambertian(const color& albedo) : tex(make_shared(albedo)) {}
lambertian(shared_ptr tex) : tex(tex) {}
bool scatter(const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered)
const override {
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double scattering_pdf(const ray& r_in, const hit_record& rec, const ray& scattered)
const override {
auto cos_theta = dot(rec.normal, unit_vector(scattered.direction()));
return cos_theta < 0 ? 0 : cos_theta/pi;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
private:
shared_ptr tex;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-lambertian-impsample]: [material.h] 중요도 샘플링을 위해 수정된 램버시안 재질]
그리고 `camera::ray_color` 함수에는 약간의 수정이 가해집니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world) const {
// If we've exceeded the ray bounce limit, no more light is gathered.
if (depth <= 0)
return color(0,0,0);
hit_record rec;
// If the ray hits nothing, return the background color.
if (!world.hit(r, interval(0.001, infinity), rec))
return background;
ray scattered;
color attenuation;
color color_from_emission = rec.mat->emitted(rec.u, rec.v, rec.p);
if (!rec.mat->scatter(r, rec, attenuation, scattered))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
double pdf_value = scattering_pdf;
color color_from_scatter =
(attenuation * scattering_pdf * ray_color(scattered, depth-1, world)) / pdf_value;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color_from_emission + color_from_scatter;
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-impsample]: [camera.h] 중요도 샘플링을 위해 수정된 ray_color 함수]
정확히 같은 그림을 얻을 것입니다. 이는 _당연한 결과_ 입니다. `ray_color`의 산란 부분에 `scattering_pdf / pdf_value`가 곱해지는데, `pdf_value`가 `scattering_pdf`와 같으므로 1을 곱하는 것과 동일하기 때문입니다.
완벽하게 일치하는 PDF 대신 균일 PDF 사용하기
----------------------------------------------
이제 경험 삼아 다른 샘플링 PDF를 사용해 보겠습니다. 반사된 광선에는 계속해서 램버시안 가중치, 즉 $\cos(\theta_o)$를 적용하고 산란 PDF도 그대로 유지하지만, 분모에는 다른 PDF를 사용하겠습니다. 반구에 대한 균일 PDF를 사용하여 샘플링할 것이므로 분모를 $1/2\pi$로 설정하겠습니다. 우리가 한 일은 PDF를 바꾼 것뿐이므로 여전히 올바른 답으로 수렴하지만, 이제 PDF가 실제 분포와 덜 완벽하게 일치하므로 수렴하는 데 더 오래 걸립니다. 이는 같은 샘플 수에서 더 노이즈가 많은 이미지를 의미합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world) const {
hit_record rec;
// If we've exceeded the ray bounce limit, no more light is gathered.
if (depth <= 0)
return color(0,0,0);
// If the ray hits nothing, return the background color.
if (!world.hit(r, interval(0.001, infinity), rec))
return background;
ray scattered;
color attenuation;
color color_from_emission = rec.mat->emitted(rec.u, rec.v, rec.p);
if (!rec.mat->scatter(r, rec, attenuation, scattered))
return color_from_emission;
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value = 1 / (2*pi);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
color color_from_scatter =
(attenuation * scattering_pdf * ray_color(scattered, depth-1, world)) / pdf_value;
return color_from_emission + color_from_scatter;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-uniform]: [camera.h] 이제 분모에 균일 PDF를 사용하는 ray_color 함수]
이전과 매우 비슷한 결과를 얻을 것이며, 노이즈가 약간 더 많을 뿐이어서 알아보기 어려울 수도 있습니다.
PDF를 반드시 산란 PDF로 되돌려 놓으십시오.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value = scattering_pdf;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-return]: [camera.h] PDF를 산란 PDF와 같도록 되돌리기]
무작위 반구 샘플링
---------------------------
이해한 내용을 확인하기 위해 다른 산란 분포를 시도해 보겠습니다. 이번에는 첫 번째 책에서 다룬 균일 반구 산란을 다시 구현해 보겠습니다. 이 기법 자체에는 아무런 문제가 없지만, 이제 더 이상 물체를 램버시안(Lambertian)으로 취급하지 않게 됩니다. 램버시안은 $\cos(\theta_o)$ 산란 분포를 요구하는 특정한 유형의 확산 재질입니다. 균일 반구 산란은 이와 다른 확산 재질입니다. 지난 절에서 했던 것처럼 재질은 그대로 두고 PDF만 바꾸면 여전히 같은 답으로 수렴하지만, 수렴에 필요한 샘플 수는 더 많아지거나 적어질 수 있습니다. 그러나 재질을 바꾸면 렌더 자체를 근본적으로 바꾼 것이 되며, 알고리즘은 다른 답으로 수렴하게 됩니다. 따라서 램버시안 확산을 균일 반구 확산으로 교체하면 렌더 결과가 _실질적으로_ 달라질 것이라고 예상해야 합니다. 산란 방향과 산란 PDF를 다음과 같이 조정하겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class lambertian : public material {
public:
lambertian(const color& albedo) : tex(make_shared(albedo)) {}
lambertian(shared_ptr tex) : tex(tex) {}
bool scatter(const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered)
const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto scatter_direction = random_on_hemisphere(rec.normal);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
// Catch degenerate scatter direction
if (scatter_direction.near_zero())
scatter_direction = rec.normal;
scattered = ray(rec.p, scatter_direction, r_in.time());
attenuation = tex->value(rec.u, rec.v, rec.p);
return true;
}
double scattering_pdf(const ray& r_in, const hit_record& rec, const ray& scattered)
const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
return 1 / (2*pi);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [scatter-mod]: [material.h] 수정된 PDF와 산란 함수]
이 새로운 확산 재질은 사실 산란 PDF가 $p(\omega_o) = \frac{1}{2\pi}$인 것에 불과합니다. 따라서 램버시안 확산에는 불완전하게 일치했던 균일 PDF가 균일 반구 확산에는 사실 완벽하게 일치합니다. 렌더링하면 약간 다른 이미지를 얻게 될 것입니다.
이전 그림과 상당히 비슷하지만, 단순한 노이즈로 볼 수 없는 차이가 있습니다. 키가 큰 상자의 앞면 색상이 훨씬 더 균일합니다. 재질에 가장 적합한 샘플링 패턴이 무엇인지 확신이 서지 않는다면 그냥 균일 PDF를 가정하는 것도 꽤 합리적인 선택이며, 수렴이 느릴 수는 있어도 렌더를 망치지는 않습니다. 다만 재질에 맞는 올바른 샘플링 패턴이 무엇인지 확신이 서지 않는다면 PDF 선택은 가장 큰 걱정거리가 아닐 것입니다. 산란 함수를 잘못 선택하면 렌더를 _반드시_ 망치게 되기 때문입니다. 적어도 잘못된 결과를 만들어 낼 것입니다. 몬테카를로 프로그램에서 찾기 가장 어려운 종류의 버그, 즉 그럴듯해 보이는 이미지를 만들어 내는 버그와 마주하게 될 수도 있습니다! 버그가 프로그램의 첫 번째 버전에 있는지, 두 번째 버전에 있는지, 아니면 둘 다에 있는지 알 수 없을 것입니다!
이 문제를 해결하기 위한 기반 구조를 만들어 보겠습니다.
무작위 방향 생성하기
====================================================================================================
이번 장과 이어지는 두 장에서는 이해를 다지고 도구를 견고하게 만들겠습니다.
Z축에 대한 무작위 방향
-----------------------------------------
먼저 무작위 방향을 생성하는 방법을 알아보겠습니다. 이미 기각법(rejection method)을 사용해 무작위 방향을 생성하는 방법이 있으므로, 이번에는 역변환법(inversion method)을 사용하는 방법을 만들어 보겠습니다. 간단히 하기 위해 $z$ 축이 표면 법선이고, $\theta$ 는 법선으로부터의 각도라고 가정하겠습니다. 이번 장에서는 모든 것을 $z$ 축을 기준으로 설정하겠습니다. 다음 장에서는 이 방향들이 표면 법선 벡터를 향하도록 맞추겠습니다. $z$ 에 대해 회전 대칭인 분포만 다루겠습니다. 따라서 $p(\omega) = f(\theta)$ 입니다.
구면 위의 방향성 PDF가 주어졌을 때($p(\omega) = f(\theta)$인 경우), $\theta$와 $\phi$에 대한 1차원 PDF는 다음과 같습니다:
$$ a(\phi) = \frac{1}{2\pi} $$
$$ b(\theta) = 2\pi f(\theta)\sin(\theta) $$
균등 난수 $r_1$과 $r_2$에 대해, CDF를 역변환하여 난수 생성기를 유도할 수 있도록 $\theta$와 $\phi$의 CDF를 구합니다.
$$ r_1 = \int_{0}^{\phi} a(\phi') d\phi' $$
$$ = \int_{0}^{\phi} \frac{1}{2\pi} d\phi' $$
$$ = \frac{\phi}{2\pi} $$
역변환하여 $\phi$를 구합니다:
$$ \phi = 2 \pi \cdot r_1 $$
이는 여러분의 직관과 일치할 것입니다. 무작위 $\phi$를 구하려면 구간 [0,1]의 균등 난수를 하나 취하여 $2\pi$를 곱하면 가능한 모든 $\phi$ 값의 전체 범위, 즉 [0,$2\pi$]를 포괄할 수 있습니다. 무작위 $\theta$ 값을 구하는 방법에 대해서는 아직 직관이 완전히 형성되지 않았을 수 있으므로, 이해를 돕기 위해 수학적 과정을 차근차근 살펴보겠습니다. 형식상 이전과 마찬가지로 $\phi$를 $\phi'$로, $\theta$를 $\theta'$로 다시 씁니다. $\theta$에 대해서는 다음과 같습니다:
$$ r_2 = \int_{0}^{\theta} b(\theta') d\theta' $$
$$ = \int_{0}^{\theta} 2 \pi f(\theta') \sin(\theta') d\theta' $$
$f()$에 대해 몇 가지 다른 함수를 시도해 보겠습니다. 먼저 구면 위의 균등 밀도를 시도해 보겠습니다. 단위 구의 넓이는 $4\pi$이므로, 단위 구 위의 균등 밀도는 $p(\omega) = \frac{1}{4\pi}$입니다.
$$ r_2 = \int_{0}^{\theta} 2 \pi \frac{1}{4\pi} \sin(\theta') d\theta' $$
$$ = \int_{0}^{\theta} \frac{1}{2} \sin(\theta') d\theta' $$
$$ = \frac{-\cos(\theta)}{2} - \frac{-\cos(0)}{2} $$
$$ = \frac{1 - \cos(\theta)}{2} $$
$\cos(\theta)$에 대해 풀면 다음과 같습니다:
$$ \cos(\theta) = 1 - 2 r_2 $$
어차피 $\cos(\theta)$만 알면 될 가능성이 높고 불필요한 $\arccos()$ 호출이 여기저기 생기는 것을 원하지 않으므로, 세타에 대해서는 풀지 않습니다.
$(\theta,\phi)$ 방향을 향하는 단위 벡터를 생성하기 위해 데카르트 좌표로 변환합니다:
$$ x = \cos(\phi) \cdot \sin(\theta) $$
$$ y = \sin(\phi) \cdot \sin(\theta) $$
$$ z = \cos(\theta) $$
그리고 항등식 $\cos^2 + \sin^2 = 1$을 사용하면, 무작위 $(r_1,r_2)$로 표현한 다음 식을 얻습니다:
$$ x = \cos(2\pi \cdot r_1)\sqrt{1 - (1-2 r_2)^2} $$
$$ y = \sin(2\pi \cdot r_1)\sqrt{1 - (1-2 r_2)^2} $$
$$ z = 1 - 2 r_2 $$
조금 단순화하면 $(1 - 2 r_2)^2 = 1 - 4r_2 + 4r_2^2$이므로 다음과 같습니다:
$$ x = \cos(2 \pi r_1) \cdot 2 \sqrt{r_2(1 - r_2)} $$
$$ y = \sin(2 \pi r_1) \cdot 2 \sqrt{r_2(1 - r_2)} $$
$$ z = 1 - 2 r_2 $$
이러한 점들을 몇 개 출력해 볼 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
int main() {
for (int i = 0; i < 200; i++) {
auto r1 = random_double();
auto r2 = random_double();
auto x = std::cos(2*pi*r1) * 2 * std::sqrt(r2*(1-r2));
auto y = std::sin(2*pi*r1) * 2 * std::sqrt(r2*(1-r2));
auto z = 1 - 2*r2;
std::cout << x << " " << y << " " << z << '\n';
}
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [rand-unit-sphere-plot]: [sphere_plot.cc] 단위 구 위의 무작위 점]
그리고 plot.ly(3D 산점도를 지원하는 훌륭한 사이트)에서 무료로 플롯할 수 있습니다:
![Figure [rand-pts-sphere]: 단위 구 위의 무작위 점
](../images/fig-3.10-rand-pts-sphere.jpg)
plot.ly 웹사이트에서 이를 회전시켜 보면 균등하게 보인다는 것을 확인할 수 있습니다.
반구의 균등 샘플링
------------------------------
이제 반구 위의 균등 분포를 유도해 보겠습니다. 반구 위에서 밀도가 균등하다는 것은 $p(\omega) = f(\theta) = \frac{1}{2\pi}$를 의미합니다. 세타 방정식의 상수만 바꾸면 다음을 얻습니다:
$$ r_2 = \int_{0}^{\theta} b(\theta') d\theta' $$
$$ = \int_{0}^{\theta} 2 \pi f(\theta') \sin(\theta') d\theta' $$
$$ = \int_{0}^{\theta} 2 \pi \frac{1}{2\pi} \sin(\theta') d\theta' $$
$$ \ldots $$
$$ \cos(\theta) = 1 - r_2 $$
이는 $\cos(\theta)$가 1에서 0까지 변하므로 $\theta$가 0에서 $\pi/2$까지 변한다는 것을 의미하며, 따라서 지평선 아래로 내려가는 것은 없습니다. 이를 플롯하는 대신, 해가 알려진 2D 적분을 풀어 보겠습니다. 반구 위에서 코사인 세제곱을 적분해 보겠습니다(해가 알려진 임의의 함수를 고른 것일 뿐입니다). 먼저 적분을 손으로 풀어 보겠습니다:
$$ \int_\omega \cos^3(\theta) dA $$
$$ = \int_{0}^{2 \pi} \int_{0}^{\pi /2} \cos^3(\theta) \sin(\theta) d\theta d\phi $$
$$ = 2 \pi \int_{0}^{\pi/2} \cos^3(\theta) \sin(\theta) d\theta = \frac{\pi}{2} $$
이제 중요도 샘플링을 사용한 적분입니다. $p(\omega) = \frac{1}{2\pi}$이므로 $f()/p() = \cos^3(\theta) / \frac{1}{2\pi}$의 평균을 구하며, 이를 다음과 같이 테스트할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
double f(double r2) {
// auto x = std::cos(2*pi*r1) * 2 * std::sqrt(r2*(1-r2));
// auto y = std::sin(2*pi*r1) * 2 * std::sqrt(r2*(1-r2));
auto z = 1 - r2;
double cos_theta = z;
return cos_theta*cos_theta*cos_theta;
}
double pdf() {
return 1.0 / (2.0*pi);
}
int main() {
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
auto r2 = random_double();
sum += f(r2) / pdf();
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "PI/2 = " << pi / 2.0 << '\n';
std::cout << "Estimate = " << sum / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [cos-cubed]: [cos_cubed.cc] $\cos^3(x)$를 사용한 적분]
반구의 코사인 샘플링
------------------------------
이제 계속해서 수평선 위에서 코사인 세제곱을 구해 보겠습니다. 다만 이번에는 PDF를 바꾸어 $p(\omega) = f(\theta) = \cos(\theta) / \pi$로 방향을 생성하겠습니다.
$$ r_2 = \int_{0}^{\theta} b(\theta') d\theta' $$
$$ = \int_{0}^{\theta} 2 \pi f(\theta') \sin(\theta') d\theta' $$
$$ = \int_{0}^{\theta} 2 \pi \frac{\cos(\theta')}{\pi} \sin(\theta') d\theta' $$
$$ = 1 - \cos^2(\theta) $$
따라서 다음과 같습니다.
$$ \cos(\theta) = \sqrt{1 - r_2} $$
다음 사실에 주목하면 특정한 경우에 대수 계산을 조금 줄일 수 있습니다.
$$ z = \cos(\theta) = \sqrt{1 - r_2} $$
$$ x = \cos(\phi) \sin(\theta) = \cos(2 \pi r_1) \sqrt{1 - z^2} = \cos(2 \pi r_1) \sqrt{r_2} $$
$$ y = \sin(\phi) \sin(\theta) = \sin(2 \pi r_1) \sqrt{1 - z^2} = \sin(2 \pi r_1) \sqrt{r_2} $$
다음은 이 PDF에 따라 가중치가 부여된 무작위 벡터를 생성하는 함수입니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
inline vec3 random_cosine_direction() {
auto r1 = random_double();
auto r2 = random_double();
auto phi = 2*pi*r1;
auto x = std::cos(phi) * std::sqrt(r2);
auto y = std::sin(phi) * std::sqrt(r2);
auto z = std::sqrt(1-r2);
return vec3(x, y, z);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [random-cosine-direction]: [vec3.h] 무작위 코사인 방향 유틸리티 함수]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "rtweekend.h"
#include
#include
double f(const vec3& d) {
auto cos_theta = d.z();
return cos_theta*cos_theta*cos_theta;
}
double pdf(const vec3& d) {
return d.z() / pi;
}
int main() {
int N = 1000000;
auto sum = 0.0;
for (int i = 0; i < N; i++) {
vec3 d = random_cosine_direction();
sum += f(d) / pdf(d);
}
std::cout << std::fixed << std::setprecision(12);
std::cout << "PI/2 = " << pi / 2.0 << '\n';
std::cout << "Estimate = " << sum / N << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [cos-density]: [cos_density.cc] 코사인 밀도 함수를 이용한 적분]
다른 밀도는 나중에 필요할 때 생성할 수 있습니다. 이 `random_cosine_direction()` 함수는 $\cos(\theta)$에 따라 가중치가 부여된 무작위 방향을 생성합니다. 여기서 $\theta$는 $z$ 축으로부터의 각도입니다.
정규 직교 기저
====================================================================================================
지난 장에서는 $z$ 축을 기준으로 무작위 방향을 생성하는 방법을 개발했습니다. 어떤 표면에서든 반사를 만들어 내려면 이 방법을 더 일반화해야 합니다. 모든 법선이 $z$ 축과 완벽하게 정렬되어 있지는 않기 때문입니다. 따라서 이 장에서는 임의의 표면 법선 벡터를 지원하도록 방법을 일반화하겠습니다.
상대 좌표
---------------------
_정규 직교 기저_ (ONB)는 서로 직교하는 세 개의 단위 벡터로 이루어진 집합입니다. 이는 좌표계의 엄격한 하위 유형입니다. 데카르트 $xyz$ 축은 정규 직교 기저의 한 예입니다. 모든 렌더 결과는 장면 속 물체들의 상대적인 위치와 방향이 카메라의 이미지 평면에 투영된 결과입니다. 이미지 평면으로의 투영이 논리적으로 정의되려면 카메라와 물체가 같은 좌표계에서 기술되어야 합니다. 그렇지 않으면 카메라가 물체를 올바르게 렌더링할 확실한 방법이 없습니다. 카메라를 물체의 좌표계에서 다시 정의하거나, 물체를 카메라의 좌표계에서 다시 정의해야 합니다. 처음부터 둘을 같은 좌표계에 두어 다시 정의할 필요가 없도록 하는 것이 가장 좋습니다. 카메라와 장면이 같은 좌표계에서 기술되기만 하면 아무 문제가 없습니다. 정규 직교 기저는 공간에서 거리와 방향이 어떻게 표현되는지를 정의하지만, 정규 직교 기저만으로는 충분하지 않습니다. 물체와 카메라는 서로 합의된 한 위치로부터의 변위로 기술되어야 합니다. 그 위치가 바로 장면의 원점 $\mathbf{O}$이며, 이는 모든 것이 변위되는 기준인 우주의 중심을 나타냅니다.
원점 $\mathbf{O}$와 데카르트 단위 벡터 $\mathbf{x}$, $\mathbf{y}$, $\mathbf{z}$가 있다고 가정하겠습니다. 어떤 위치가 (3,-2,7)이라고 말할 때, 실제로는 다음을 의미합니다.
$$ \text{Location is } \mathbf{O} + 3\mathbf{x} - 2\mathbf{y} + 7\mathbf{z} $$
원점이 $\mathbf{O}'$이고 기저 벡터가 $\mathbf{u}$, $\mathbf{v}$, $\mathbf{w}$인 다른 좌표계에서 좌표를 측정하려면, 다음을 만족하는 수 $(u,v,w)$를 찾으면 됩니다.
$$ \text{Location is } \mathbf{O}' + u\mathbf{u} + v\mathbf{v} + w\mathbf{w} $$
정규 직교 기저 생성
--------------------------------
그래픽스 입문 강좌를 수강하면 좌표계와 4×4 좌표 변환 행렬에 많은 시간을 쓰게 됩니다. 주의 깊게 공부하십시오. 정말 중요한 내용입니다! 하지만 이 책에서는 그것이 필요하지 않으므로 그것 없이 해 나가겠습니다. 우리에게 필요한 것은 표면 법선 벡터 $\mathbf{n}$을 기준으로 정해진 분포를 따르는 무작위 방향을 생성하는 것입니다. 방향은 상대적이며 특정한 원점을 갖지 않으므로 이 작업에는 원점이 필요하지 않습니다. 우선, 각각 $\mathbf{n}$에 수직이면서 서로에게도 수직인 두 개의 공탄젠트 벡터가 필요합니다.
일부 3D 객체 모델에는 각 정점마다 하나 이상의 공탄젠트 벡터가 함께 제공됩니다. 모델에 공탄젠트 벡터가 하나뿐이라면 ONB를 만드는 과정은 간단하지 않습니다. 길이가 0이 아니고 $\mathbf{n}$과 평행하지 않은 임의의 벡터 $\mathbf{a}$가 있다고 가정해 보겠습니다. $\mathbf{n} \times \mathbf{a}$가 $\mathbf{n}$과 $\mathbf{a}$ 모두에 수직이라는 외적의 성질을 이용하면 $\mathbf{n}$에 수직인 벡터 $\mathbf{s}$와 $\mathbf{t}$를 얻을 수 있습니다:
$$ \mathbf{s} = \operatorname{unit\_vector}(\mathbf{n} \times \mathbf{a}) $$
$$ \mathbf{t} = \mathbf{n} \times \mathbf{s} $$
여기까지는 좋지만, 문제는 모델을 불러올 때 $\mathbf{a}$가 주어지지 않을 수 있고, 현재 프로그램에는 그것을 생성할 방법이 없다는 점입니다. 초기 벡터로 사용할 $\mathbf{a}$를 임의로 고른다면 $\mathbf{n}$과 평행한 $\mathbf{a}$를 얻게 될 수도 있습니다. 따라서 일반적인 방법은 임의의 축을 하나 골라 그것이 (단위 길이라고 가정하는) $\mathbf{n}$과 평행한지 확인하고, 평행하다면 다른 축을 사용하는 것입니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (std::fabs(n.x()) > 0.9)
a = vec3(0, 1, 0)
else
a = vec3(1, 0, 0)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
그런 다음 외적을 취하여 $\mathbf{s}$와 $\mathbf{t}$를 얻습니다
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vec3 s = unit_vector(cross(n, a));
vec3 t = cross(n, s);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
$\mathbf{t}$에 대해서는 단위 벡터를 구할 필요가 없다는 점에 유의하십시오. $\mathbf{n}$과 $\mathbf{s}$가 모두 단위 벡터이므로 그 외적인 $\mathbf{t}$도 단위 벡터가 됩니다. $\mathbf{s}$, $\mathbf{t}$, $\mathbf{n}$으로 이루어진 ONB가 있고 $z$ 축을 기준으로 한 무작위 $(x,y,z)$가 있으면, 다음과 같이 $\mathbf{n}$을 기준으로 한 벡터를 얻을 수 있습니다:
$$ \text{random vector} = x \mathbf{s} + y \mathbf{t} + z \mathbf{n} $$
기억하시겠지만, 카메라에서 광선을 생성할 때도 비슷한 수학을 사용했습니다. 그것을 카메라의 고유 좌표계로의 변환이라고 생각할 수 있습니다.
ONB 클래스
--------------
ONB를 위한 클래스를 만들어야 할까요, 아니면 유틸리티 함수로 충분할까요? 확실하지는 않지만, 유틸리티 함수보다 딱히 더 복잡하지도 않을 것이므로 클래스를 만들겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#ifndef ONB_H
#define ONB_H
class onb {
public:
onb(const vec3& n) {
axis[2] = unit_vector(n);
vec3 a = (std::fabs(axis[2].x()) > 0.9) ? vec3(0,1,0) : vec3(1,0,0);
axis[1] = unit_vector(cross(axis[2], a));
axis[0] = cross(axis[2], axis[1]);
}
const vec3& u() const { return axis[0]; }
const vec3& v() const { return axis[1]; }
const vec3& w() const { return axis[2]; }
vec3 transform(const vec3& v) const {
// Transform from basis coordinates to local space.
return (v[0] * axis[0]) + (v[1] * axis[1]) + (v[2] * axis[2]);
}
private:
vec3 axis[3];
};
#endif
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-onb]: [onb.h] 정규 직교 기저 클래스]
이를 사용하여 램버시안 재질을 다음과 같이 다시 작성할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "hittable.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
#include "onb.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "texture.h"
class material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
virtual bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const {
return false;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
};
class lambertian : public material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const override {
onb uvw(rec.normal);
auto scatter_direction = uvw.transform(random_cosine_direction());
scattered = ray(rec.p, unit_vector(scatter_direction), r_in.time());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
attenuation = tex->value(rec.u, rec.v, rec.p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
pdf = dot(uvw.w(), scattered.direction()) / pi;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return true;
}
...
};
class metal : public material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
};
class dielectric : public material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
};
class diffuse_light : public material {
...
};
class isotropic : public material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [scatter-onb]: [material.h] 정규 직교 기저를 사용한 산란 함수]
그리고 여기서는 카메라 클래스에 이에 따른 변경 사항을 추가합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world) const {
...
ray scattered;
color attenuation;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
color color_from_emission = rec.mat->emitted(rec.u, rec.v, rec.p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
if (!rec.mat->scatter(r, rec, attenuation, scattered, pdf_value))
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color_from_emission;
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
pdf_value = scattering_pdf;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [scatter-ray-color]: [camera.h] 반환된 PDF 값을 사용하도록 수정한 ray_color 함수]
그 결과는 다음과 같습니다:
노이즈를 어느 정도 없애 보겠습니다.
하지만 먼저 `isotropic` 재질을 간단히 수정하겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class isotropic : public material {
public:
isotropic(const color& albedo) : tex(make_shared(albedo)) {}
isotropic(shared_ptr tex) : tex(tex) {}
bool scatter(
const ray& r_in, const hit_record& rec, color& attenuation, ray& scattered, double& pdf
) const override {
scattered = ray(rec.p, random_unit_vector(), r_in.time());
attenuation = tex->value(rec.u, rec.v, rec.p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
pdf = 1 / (4 * pi);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return true;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double scattering_pdf(const ray& r_in, const hit_record& rec, const ray& scattered)
const override {
return 1 / (4 * pi);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
private:
shared_ptr tex;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-isotropic-impsample]: [material.h] 중요도 샘플링을 위해 수정한 등방성 재질]
광원 직접 샘플링하기
====================================================================================================
모든 방향에 대해 균일하게 샘플링할 때의 문제는 광원이 임의의 방향이나 중요하지 않은 방향보다 더 자주 샘플링되지 않는다는 점입니다. 그림자 광선을 사용하여 주어진 지점에서의 직접 조명을 구할 수도 있습니다. 대신, 저는 광원 쪽으로 더 많은 광선을 보내는 PDF를 사용하겠습니다. 그런 다음 그 PDF를 바꾸어 원하는 어떤 방향으로든 더 많은 광선을 보낼 수 있습니다.
광원을 향하는 무작위 방향을 고르는 것은 정말 쉽습니다. 광원 위의 무작위 점을 하나 고르고 그 방향으로 광선을 보내기만 하면 됩니다. 하지만 렌더링 결과가 편향되지 않도록 PDF, 즉 $p(\omega)$를 알아야 합니다. 그렇다면 그것은 무엇입니까?
광원의 PDF 구하기
---------------------------
표면적이 $A$인 광원에서 균일하게 샘플링한다면, 표면 위의 PDF는 단순히 $\frac{1}{A}$입니다. 광원의 전체 표면을 단위 구에 다시 투영하면 얼마만큼의 면적을 차지합니까? 다행히도 다음 도식에 나타낸 것처럼 간단한 대응 관계가 있습니다:
![Figure [shape-onto-pdf]: 광원 형상의 PDF 위로의 투영
](../images/fig-3.11-shape-onto-pdf.jpg)
광원 위의 작은 면적 $dA$를 보면, 이를 샘플링할 확률은 $\operatorname{p_q}(q) \cdot dA$입니다. 구 위에서 구의 작은 면적 $d\omega$를 샘플링할 확률은 $\operatorname{p}(\omega) \cdot d\omega$입니다. $d\omega$와 $dA$ 사이에는 기하학적 관계가 있습니다:
$$ d\omega = \frac{dA \cdot \cos(\theta)}{\operatorname{distance}^2(p,q)} $$
$d\omega$와 $dA$를 샘플링할 확률은 같아야 하므로, 다음이 성립합니다.
$$ \operatorname{p}(\omega) \cdot d\omega = \operatorname{p_q}(q) \cdot dA $$
$$ \operatorname{p}(\omega)
\cdot \frac{dA \cdot \cos(\theta)}{\operatorname{distance}^2(p,q)}
= \operatorname{p_q}(q) \cdot dA $$
광원 위에서 균일하게 샘플링하면 표면 위의 PDF가 $\frac{1}{A}$라는 것을 알고 있습니다:
$$ \operatorname{p_q}(q) = \frac{1}{A} $$
$$ \operatorname{p}(\omega) \cdot \frac{dA \cdot \cos(\theta)}{\operatorname{distance}^2(p,q)}
= \frac{dA}{A} $$
따라서
$$ \operatorname{p}(\omega) = \frac{\operatorname{distance}^2(p,q)}{\cos(\theta) \cdot A} $$
광원 샘플링
---------------
수학과 개념을 올바르게 이해했는지 확인하기 위해, `ray_color()` 함수를 임시로 고쳐 매우 하드코딩된 방식으로 광원을 샘플링하도록 할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world) const {
// If we've exceeded the ray bounce limit, no more light is gathered.
if (depth <= 0)
return color(0,0,0);
hit_record rec;
// If the ray hits nothing, return the background color.
if (!world.hit(r, interval(0.001, infinity), rec))
return background;
ray scattered;
color attenuation;
double pdf_value;
color color_from_emission = rec.mat->emitted(rec.u, rec.v, rec.p);
if (!rec.mat->scatter(r, rec, attenuation, scattered, pdf_value))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto on_light = point3(random_double(213,343), 554, random_double(227,332));
auto to_light = on_light - rec.p;
auto distance_squared = to_light.length_squared();
to_light = unit_vector(to_light);
if (dot(to_light, rec.normal) < 0)
return color_from_emission;
double light_area = (343-213)*(332-227);
auto light_cosine = std::fabs(to_light.y());
if (light_cosine < 0.000001)
return color_from_emission;
pdf_value = distance_squared / (light_cosine * light_area);
scattered = ray(rec.p, to_light, r.time());
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
color color_from_scatter =
(attenuation * scattering_pdf * ray_color(scattered, depth-1, world)) / pdf_value;
return color_from_emission + color_from_scatter;
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-lights]: [camera.h] 광원 샘플링을 적용한 광선 색상]
이 장면을 픽셀당 샘플 10개만으로 테스트하겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
cam.aspect_ratio = 1.0;
cam.image_width = 600;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cam.samples_per_pixel = 10;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
cam.max_depth = 50;
cam.background = color(0,0,0);
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-lights-10spp]: [main.cc] 10spp에서 광원 샘플링을 적용한 광선 색상]
픽셀당 샘플 10개로 렌더링하면 다음 결과를 얻습니다:
이는 광원만 샘플링하는 방식에서 기대할 수 있는 결과와 거의 같으므로, 제대로 동작하는 것으로 보입니다.
단방향 광원으로 전환하기
----------------------------------
천장의 광원 주변에 튀는 노이즈 점들이 생기는 것은 광원이 양면이고 광원과 천장 사이에 작은 틈이 있기 때문입니다. 광원이 아래쪽으로만 빛을 방출하도록 하는 것이 좋을 것입니다. `hittable::emitted()` 함수가 추가 정보를 받도록 하여 이를 구현할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
virtual color emitted(
const ray& r_in, const hit_record& rec, double u, double v, const point3& p
) const {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color(0,0,0);
}
...
};
class diffuse_light : public material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
color emitted(const ray& r_in, const hit_record& rec, double u, double v, const point3& p)
const override {
if (!rec.front_face)
return color(0,0,0);
return tex->value(u, v, p);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [emitted-directional]: [material.h] 재질 방출, 방향성]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
color ray_color(const ray& r, int depth, const hittable& world) const {
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
color color_from_emission = rec.mat->emitted(r, rec, rec.u, rec.v, rec.p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [emitted-ray-color]: [camera.h] 재질 방출, camera::ray_color() 변경 사항]
그러면 다음과 같은 결과를 얻습니다:
혼합 밀도
====================================================================================================
우리는 $\cos(\theta)$와 관련된 PDF와 광원 샘플링과 관련된 PDF를 사용했습니다. 이 둘을 결합한 PDF가 필요합니다.
PDF 클래스
-------------
우리는 이미 꽤 많은 코드에서 PDF를 다루어 왔습니다. 이제 PDF 사용 방식을 어떻게 표준화할지 정하기에 좋은 시점이라고 생각합니다. 표면에 대한 PDF와 광원에 대한 PDF가 필요하다는 것은 이미 알고 있으므로, `pdf` 기본 클래스를 만들어 보겠습니다. 지금까지는 방향을 받아 그 방향에 대한 PDF의 분포 값을 반환하는 `pdf()` 함수를 사용했습니다. 지금까지 이 값은 $1/4\pi$, $1/2\pi$, $\cos(\theta)/\pi$ 중 하나였습니다. 몇몇 예제에서는 PDF의 분포와 다른 분포를 사용하여 무작위 방향을 생성했습니다. 이에 대해서는 중요도 샘플링 다루기 장에서 상당히 자세히 다루었습니다. 일반적으로 무작위 방향의 분포를 알고 있다면, 같은 분포를 가진 PDF를 사용해야 합니다. 그렇게 하면 가장 빠르게 수렴합니다. 이 점을 염두에 두고, 무작위 방향을 생성하고 PDF 값을 결정하는 역할을 담당하는 `pdf` 클래스를 만들겠습니다.
이 모든 것을 종합하면, 모든 `pdf` 클래스는 다음을 담당해야 합니다.
1. 내부 PDF 분포에 따라 가중된 무작위 방향을 반환하는 것, 그리고
2. 해당 방향에 대응하는 PDF 분포 값을 반환하는 것입니다.
내부적으로 이를 수행하는 세부 방식은 $\operatorname{pSurface}$와 $\operatorname{pLight}$에 따라 다르지만, 클래스 계층 구조는 바로 이런 경우를 위해 고안된 것입니다! 추상 클래스에 무엇을 넣어야 할지는 결코 분명하지 않으므로, 저는 욕심을 부려 최소한의 인터페이스로 충분하기를 바라는 방식을 택하며, `pdf`의 경우 이는 다음을 의미합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#ifndef PDF_H
#define PDF_H
#include "onb.h"
class pdf {
public:
virtual ~pdf() {}
virtual double value(const vec3& direction) const = 0;
virtual vec3 generate() const = 0;
};
#endif
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-pdf]: [pdf.h] 추상 PDF 클래스]
하위 클래스를 구체화해 보면서 `pdf`에 다른 것을 추가해야 하는지 살펴보겠습니다. 먼저 단위 구 위의 균일 밀도를 만들겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class sphere_pdf : public pdf {
public:
sphere_pdf() {}
double value(const vec3& direction) const override {
return 1/ (4 * pi);
}
vec3 generate() const override {
return random_unit_vector();
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-uni-pdf]: [pdf.h] sphere_pdf 클래스]
다음으로 코사인 밀도를 시도해 보겠습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class cosine_pdf : public pdf {
public:
cosine_pdf(const vec3& w) : uvw(w) {}
double value(const vec3& direction) const override {
auto cosine_theta = dot(unit_vector(direction), uvw.w());
return std::fmax(0, cosine_theta/pi);
}
vec3 generate() const override {
return uvw.transform(random_cosine_direction());
}
private:
onb uvw;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-cos-pdf]: [pdf.h] cosine_pdf 클래스]
이 코사인 PDF를 `ray_color()` 함수에서 사용해 볼 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "hittable.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
#include "pdf.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "material.h"
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world) const {
// If we've exceeded the ray bounce limit, no more light is gathered.
if (depth <= 0)
return color(0,0,0);
hit_record rec;
// If the ray hits nothing, return the background color.
if (!world.hit(r, interval(0.001, infinity), rec))
return background;
ray scattered;
color attenuation;
double pdf_value;
color color_from_emission = rec.mat->emitted(r, rec, rec.u, rec.v, rec.p);
if (!rec.mat->scatter(r, rec, attenuation, scattered, pdf_value))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cosine_pdf surface_pdf(rec.normal);
scattered = ray(rec.p, surface_pdf.generate(), r.time());
pdf_value = surface_pdf.value(scattered.direction());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
color color_from_scatter =
(attenuation * scattering_pdf * ray_color(scattered, depth-1, world)) / pdf_value;
return color_from_emission + color_from_scatter;
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-cos-pdf]: [camera.h] 코사인 PDF를 사용하는 ray_color 함수]
그리고 렌더링을 다시 픽셀당 1000개의 샘플로 설정합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
cam.aspect_ratio = 1.0;
cam.image_width = 600;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cam.samples_per_pixel = 1000;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
cam.max_depth = 50;
cam.background = color(0,0,0);
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [cosine-density-1000spp]: [main.cc] 샘플링을 다시 1000spp로 재설정]
이렇게 하면 정확히 같은 결과가 나오므로, 지금까지 한 일은 일부 계산을 `cosine_pdf` 클래스로 옮긴 것뿐입니다:
hittable을 향한 방향 샘플링
---------------------------------------
이제 광원처럼 `hittable`을 향하는 방향을 샘플링해 볼 수 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
#include "hittable_list.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "onb.h"
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
class hittable_pdf : public pdf {
public:
hittable_pdf(const hittable& objects, const point3& origin)
: objects(objects), origin(origin)
{}
double value(const vec3& direction) const override {
return objects.pdf_value(origin, direction);
}
vec3 generate() const override {
return objects.random(origin);
}
private:
const hittable& objects;
point3 origin;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-hittable-pdf]: [pdf.h] hittable_pdf 클래스]
광원을 샘플링하려면 `hittable`이 아직 인터페이스를 갖추지 못한 몇 가지 질의에 답할 수 있어야 합니다. 위 코드는 `hittable` 클래스에 아직 구현되지 않은 두 함수 `pdf_value()`와 `random()`이 존재한다고 가정합니다. 프로그램을 컴파일하려면 이 함수들을 추가해야 합니다. 모든 `hittable` 하위 클래스를 하나씩 살펴보며 이 함수들을 추가할 수도 있지만 번거로우므로, `hittable` 기본 클래스에 간단한 함수 두 개만 추가하겠습니다. 이렇게 하면 이전의 순수 추상 구현이 깨지지만, 작업량을 줄일 수 있습니다. 순수 추상 `hittable` 인터페이스 클래스를 원한다면 이 함수들을 하위 클래스까지 자유롭게 구현하십시오.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class hittable {
public:
virtual ~hittable() = default;
virtual bool hit(const ray& r, interval ray_t, hit_record& rec) const = 0;
virtual aabb bounding_box() const = 0;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
virtual double pdf_value(const point3& origin, const vec3& direction) const {
return 0.0;
}
virtual vec3 random(const point3& origin) const {
return vec3(1,0,0);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [hittable-plus2]: [hittable.h] 두 개의 새 메서드가 추가된 hittable 클래스]
그런 다음 `quad`가 이 함수들을 구현하도록 변경합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class quad : public hittable {
public:
quad(const point3& Q, const vec3& u, const vec3& v, shared_ptr mat)
: Q(Q), u(u), v(v), mat(mat)
{
auto n = cross(u, v);
normal = unit_vector(n);
D = dot(normal, Q);
w = n / dot(n,n);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
area = n.length();
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
set_bounding_box();
}
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value(const point3& origin, const vec3& direction) const override {
hit_record rec;
if (!this->hit(ray(origin, direction), interval(0.001, infinity), rec))
return 0;
auto distance_squared = rec.t * rec.t * direction.length_squared();
auto cosine = std::fabs(dot(direction, rec.normal) / direction.length());
return distance_squared / (cosine * area);
}
vec3 random(const point3& origin) const override {
auto p = Q + (random_double() * u) + (random_double() * v);
return p - origin;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
private:
point3 Q;
vec3 u, v;
vec3 w;
shared_ptr mat;
aabb bbox;
vec3 normal;
double D;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double area;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [quad-pdf]: [quad.h] PDF를 갖춘 사각형]
`quad`에는 `pdf_value()`와 `random()`만 추가하면 됩니다. 이를 광원을 중요도 샘플링하는 데 사용하며, 장면에 있는 유일한 광원이 `quad`이기 때문입니다. 다른 광원 기하를 원하거나 다른 객체에 PDF를 사용하려면 해당 클래스에 위 함수들을 구현해야 합니다.
카메라의 `render()` 함수에 `lights` 매개변수를 추가합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
void render(const hittable& world, const hittable& lights) {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
initialize();
std::cout << "P3\n" << image_width << ' ' << image_height << "\n255\n";
for (int j = 0; j < image_height; j++) {
std::clog << "\rScanlines remaining: " << (image_height - j) << ' ' << std::flush;
for (int i = 0; i < image_width; i++) {
color pixel_color(0,0,0);
for (int s_j = 0; s_j < sqrt_spp; s_j++) {
for (int s_i = 0; s_i < sqrt_spp; s_i++) {
ray r = get_ray(i, j, s_i, s_j);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
pixel_color += ray_color(r, max_depth, world, lights);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
}
write_color(std::cout, pixel_samples_scale * pixel_color);
}
}
std::clog << "\rDone. \n";
}
...
private:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
color ray_color(const ray& r, int depth, const hittable& world, const hittable& lights)
const {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
ray scattered;
color attenuation;
double pdf_value;
color color_from_emission = rec.mat->emitted(r, rec, rec.u, rec.v, rec.p);
if (!rec.mat->scatter(r, rec, attenuation, scattered, pdf_value))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
hittable_pdf light_pdf(lights, rec.p);
scattered = ray(rec.p, light_pdf.generate(), r.time());
pdf_value = light_pdf.value(scattered.direction());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
color sample_color = ray_color(scattered, depth-1, world, lights);
color color_from_scatter = (attenuation * scattering_pdf * sample_color) / pdf_value;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color_from_emission + color_from_scatter;
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-lights]: [camera.h] 광원 PDF를 사용하는 ray_color 함수]
천장 가운데에 광원을 만듭니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
// Box 2
shared_ptr box2 = box(point3(0,0,0), point3(165,165,165), white);
box2 = make_shared(box2, -18);
box2 = make_shared(box2, vec3(130,0,65));
world.add(box2);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
// Light Sources
auto empty_material = shared_ptr();
quad lights(point3(343,554,332), vec3(-130,0,0), vec3(0,0,-105), empty_material);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
camera cam;
cam.aspect_ratio = 1.0;
cam.image_width = 600;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cam.samples_per_pixel = 10;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
cam.max_depth = 50;
cam.background = color(0,0,0);
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
cam.render(world, lights);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-hittable-pdf]: [main.cc] 코넬 박스에 광원 추가하기]
픽셀당 샘플 10개로 다음 결과를 얻습니다:
혼합 PDF 클래스
----------------------
중요도 샘플링 다루기 장에서 간단히 언급했듯이, 임의의 PDF들을 선형 혼합하여 그 자체로도 PDF인 혼합 밀도를 만들 수 있습니다. PDF들의 가중 평균은 무엇이든 역시 PDF입니다. 가중치가 양수이고 합이 1이기만 하면 새로운 PDF를 얻습니다.
$$ \operatorname{pMixture}() = w_0 p_0() + w_1 p_1() + w_2 p_2() + \ldots + w_{n-1} p_{n-1}() $$
$$ 1 = w_0 + w_1 + w_2 + \ldots + w_{n-1} $$
예를 들어, 두 밀도를 단순히 평균할 수 있습니다:
$$ \operatorname{pMixture}(\omega_o)
= \frac{1}{2} \operatorname{pSurface}(\omega_o) + \frac{1}{2} \operatorname{pLight}(\omega_o)
$$
이를 위해 코드를 어떻게 구성해야 할까요? 이것이 예상만큼 쉽지 않은 이유가 되는 매우 중요한 세부 사항이 하나 있습니다. 혼합 PDF에 대한 무작위 방향을 생성하는 것은 간단합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (random_double() < 0.5)
pick direction according to pSurface
else
pick direction according to pLight
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
하지만 $\operatorname{pMixture}$의 PDF 값을 구하는 것은 조금 더 까다롭습니다. 단순히
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (direction is from pSurface)
get PDF value of pSurface
else
get PDF value of pLight
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
우선, 무작위 방향이 어느 PDF에서 나왔는지 알아내는 것은 아마 간단하지 않을 것입니다. `generate()`가 원래의 `random_double()` 값이 무엇이었는지 `value()`에 알려 줄 수 있는 연결 장치가 전혀 없으므로, 무작위 방향이 어느 PDF에서 나왔는지 쉽게 말할 수 없습니다. 위의 생각이 옳다고 본다면, 그 방향이 어느 PDF에서 나올 수 있었는지 알아내기 위해 거꾸로 풀어야 할 것입니다. 솔직히 악몽처럼 들리지만, 다행히 그렇게 할 필요는 없습니다. 두 PDF 모두가 생성할 수 있었던 방향들이 있습니다. 예를 들어, 광원을 향하는 방향은 $\operatorname{pLight}$ _또는_ $\operatorname{pSurface}$ 어느 쪽에 의해서도 생성될 수 있었습니다. 무작위 방향에 대해 $\operatorname{pSurface}$의 PDF 값과 $\operatorname{pLight}$의 PDF 값을 구한 다음, PDF 혼합 가중치를 적용하여 그 방향의 전체 PDF 값을 구하는 것으로 충분합니다. 혼합 밀도 클래스는 사실 꽤 간단합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class mixture_pdf : public pdf {
public:
mixture_pdf(shared_ptr p0, shared_ptr p1) {
p[0] = p0;
p[1] = p1;
}
double value(const vec3& direction) const override {
return 0.5 * p[0]->value(direction) + 0.5 * p[1]->value(direction);
}
vec3 generate() const override {
if (random_double() < 0.5)
return p[0]->generate();
else
return p[1]->generate();
}
private:
shared_ptr p[2];
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [class-mixturep-df]: [pdf.h] mixture_pdf 클래스]
이제 코사인 샘플링과 광원 샘플링의 혼합 밀도를 만들고자 합니다. 이를 `ray_color()`에 연결할 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world, const hittable& lights)
const {
...
if (!rec.mat->scatter(r, rec, attenuation, scattered, pdf_value))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto p0 = make_shared(lights, rec.p);
auto p1 = make_shared(rec.normal);
mixture_pdf mixed_pdf(p0, p1);
scattered = ray(rec.p, mixed_pdf.generate(), r.time());
pdf_value = mixed_pdf.value(scattered.direction());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
color sample_color = ray_color(scattered, depth-1, world, lights);
color color_from_scatter = (attenuation * scattering_pdf * sample_color) / pdf_value;
return color_from_emission + color_from_scatter;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-mixture]: [camera.h] 혼합 PDF를 사용하는 ray_color 함수]
`main.cc`를 픽셀당 1000개 샘플로 업데이트하면(코드는 생략) 다음과 같은 결과가 나옵니다:
몇 가지 아키텍처 결정
====================================================================================================
이 장에서는 코드를 작성하지 않습니다. 우리는 갈림길에 서 있으며 몇 가지 아키텍처 결정을 내려야 합니다.
혼합 밀도 방식은 보다 전통적인 그림자 광선을 사용하는 방식의 대안입니다. 그림자 광선은 교차점에서 주어진 광원까지 가로막히지 않은 경로가 있는지 확인하는 광선입니다. 한 점과 주어진 광원 사이에서 물체와 교차하는 광선은 그 교차점이 해당 광원의 그림자 안에 있음을 나타냅니다. 혼합 밀도 방식은 제가 개인적으로 선호하는 방식인데, 광원뿐만 아니라 창문이나 문 아래의 밝은 틈, 또는 그 밖에 밝을 것 같은 -- 혹은 중요할 것 같은 -- 무엇이든 샘플링할 수 있기 때문입니다. 하지만 대부분의 전문적인 경로 추적기에서는 여전히 그림자 광선을 볼 수 있을 것입니다. 일반적으로 이들은 미리 정해진 개수의 그림자 광선(_예를 들어_ 1, 4, 8, 16개)을 두고, 렌더링 과정에서 경로 추적 광선이 교차하는 각 지점마다 이러한 종단 그림자 광선을 장면 내의 무작위 광원으로 보내 해당 교차점이 그 무작위 광원에 의해 비춰지는지 판단합니다. 교차점은 그 광원에 의해 비춰지거나 완전히 그림자 안에 있게 되며, 그림자 광선이 많을수록 더 정확한 조명을 얻을 수 있습니다. 모든 그림자 광선이 (광원 또는 가리는 표면에서) 종료된 후, 초기 경로 추적 광선은 계속 진행하며 다음 교차점에서 더 많은 그림자 광선이 보내집니다. 그림자 광선에게 무엇이 중요한지는 알려 줄 수 없고 무엇이 발광하는지만 알려 줄 수 있으므로, 그림자 광선은 광자 분포가 지나치게 복잡하지 않은 비교적 단순한 장면에서 가장 잘 작동합니다. 그렇기는 하지만, 그림자 광선은 처음 부딪히는 물체에서 종료되고 여기저기 튕기지 않으므로 그림자 광선 하나는 경로 추적 광선 하나보다 비용이 적게 들며, 이것이 일반적으로 경로 추적 광선보다 훨씬 많은 그림자 광선(_예를 들어_ 1, 4, 8, 16개)을 보게 되는 이유입니다. 보다 제한된 장면에서는 혼합 밀도 대신 그림자 광선을 선택할 수도 있으며, 이는 개인적인 설계 취향의 문제입니다. 그림자 광선은 대략적인 결과를 얻는 데 혼합 밀도보다 비용이 적게 드는 경향이 있으며, 실시간 렌더링에서 점점 더 널리 사용되고 있습니다.
코드에는 그 밖에도 몇 가지 이슈가 있습니다.
PDF 구성이 `ray_color()` 함수에 하드 코딩되어 있습니다. 이를 정리해야 합니다.
실수로 정반사 광선(유리와 금속)을 망가뜨려서 이제 더 이상 지원되지 않습니다. 산란 함수를 델타 함수로 만들기만 하면 수학적으로는 계속 성립하겠지만, 온갖 부동 소수점 재앙을 초래할 것입니다. 정반사를 $f()/p()$를 건너뛰는 특수한 경우로 만들 수도 있고, 표면 거칠기를 아주 작지만 -- 0은 아닌 -- 값으로 설정하여 완벽하게 매끄러워 보이면서도 NaN을 생성하지 않는 거의 거울에 가까운 표면을 만들 수도 있습니다. 어느 방법을 택할지에 대해서는 특별한 의견이 없지만(두 방법을 모두 시도해 보았고 둘 다 장점이 있습니다), 어차피 매끄러운 금속과 유리 코드가 있으므로 명시적인 $f()/p()$ 계산을 그냥 건너뛰는 완벽한 정반사 표면을 추가하겠습니다.
또한 환경 맵이나 더 흥미로운 함수형 배경을 추가하고 싶을 경우에 대비한 제대로 된 배경 함수 기반 구조도 없습니다. 일부 환경 맵은 HDR입니다(RGB 성분이 0–255 바이트가 아니라 정규화된 부동소수점 값입니다). 우리의 출력은 처음부터 줄곧 HDR이었으며, 단지 그것을 잘라내고 있었을 뿐입니다.
마지막으로, 우리의 렌더러는 RGB입니다. 좀 더 물리 기반인 렌더러는 -- 자동차 제조사가 사용할 법한 것처럼 -- 아마도 분광 색상과 어쩌면 편광까지 사용해야 할 것입니다. 영화 렌더러의 경우, 대부분의 스튜디오는 여전히 RGB로도 충분히 해내고 있습니다. 두 모드를 모두 갖춘 하이브리드 렌더러를 만들 수도 있지만, 물론 그것은 더 어렵습니다. 지금은 RGB를 고수하겠지만, 이 책의 마지막 부분에서 이 주제를 다루겠습니다.
PDF 관리 정리하기
====================================================================================================
지금까지는 `ray_color()` 함수가 하드코딩된 두 개의 PDF를 생성하도록 했습니다:
1. 광원의 형태와 관련된 `p0()`
2. 법선 벡터 및 표면 유형과 관련된 `p1()`
광원(또는 샘플링하려는 어떤 `hittable`이든)에 관한 정보를 `ray_color()` 함수에 전달할 수 있으며, `material` 함수에 PDF를 요청할 수도 있습니다(그렇게 하려면 계측 코드를 추가해야 합니다). 또한 산란된 광선이 정반사인지 알아야 하는데, 이는 `hit()` 함수에 묻거나 `material` 클래스에 물어서 알아낼 수 있습니다.
확산 대 정반사
------------------------
우리가 허용하고 싶은 것 중 하나는 -- 니스 칠한 나무처럼 -- 부분적으로는 이상적인 정반사(광택)이고 부분적으로는 확산(나무)인 재질입니다. 일부 렌더러는 재질이 두 개의 광선, 즉 정반사 광선 하나와 확산 광선 하나를 생성하도록 합니다. 저는 분기를 좋아하지 않으므로, 재질이 확산인지 정반사인지를 무작위로 결정하도록 하는 편을 택하겠습니다. 이 접근 방식의 함정은 PDF 값을 요청할 때 주의해야 하며, `ray_color()`가 이 광선이 확산인지 정반사인지 알고 있어야 한다는 점입니다. 다행히 확산인 경우에만 `pdf_value()`를 호출하기로 결정했으므로, 이를 암묵적으로 처리할 수 있습니다.
`hittable`에서 했던 것처럼 `material`을 재설계하여 새로운 인수들을 모두 하나의 클래스에 담을 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "hittable.h"
#include "onb.h"
#include "texture.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
class scatter_record {
public:
color attenuation;
shared_ptr pdf_ptr;
bool skip_pdf;
ray skip_pdf_ray;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class material {
public:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
virtual bool scatter(const ray& r_in, const hit_record& rec, scatter_record& srec) const {
return false;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [material-refactor]: [material.h] material 클래스 리팩터링하기]
`lambertian` 재질은 더 간단해집니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "hittable.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ delete
#include "onb.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
#include "pdf.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "texture.h"
...
class lambertian : public material {
public:
lambertian(const color& albedo) : tex(make_shared(albedo)) {}
lambertian(shared_ptr tex) : tex(tex) {}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(const ray& r_in, const hit_record& rec, scatter_record& srec) const override {
srec.attenuation = tex->value(rec.u, rec.v, rec.p);
srec.pdf_ptr = make_shared(rec.normal);
srec.skip_pdf = false;
return true;
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double scattering_pdf(const ray& r_in, const hit_record& rec, const ray& scattered)
const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto cos_theta = dot(rec.normal, unit_vector(scattered.direction()));
return cos_theta < 0 ? 0 : cos_theta/pi;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
private:
shared_ptr tex;
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [lambertian-scatter]: [material.h] 새로운 lambertian scatter() 메서드]
그리고 `ray_color()`의 변경 사항은 적습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world, const hittable& lights)
const {
// If we've exceeded the ray bounce limit, no more light is gathered.
if (depth <= 0)
return color(0,0,0);
hit_record rec;
// If the ray hits nothing, return the background color.
if (!world.hit(r, interval(0.001, infinity), rec))
return background;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
scatter_record srec;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
color color_from_emission = rec.mat->emitted(r, rec, rec.u, rec.v, rec.p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
if (!rec.mat->scatter(r, rec, srec))
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
auto light_ptr = make_shared(lights, rec.p);
mixture_pdf p(light_ptr, srec.pdf_ptr);
ray scattered = ray(rec.p, p.generate(), r.time());
auto pdf_value = p.value(scattered.direction());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
color sample_color = ray_color(scattered, depth-1, world, lights);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
color color_from_scatter =
(srec.attenuation * scattering_pdf * sample_color) / pdf_value;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return color_from_emission + color_from_scatter;
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-mixture]: [camera.h] 혼합 PDF를 사용하는 ray_color 함수]
정반사 처리하기
------------------
아직 정반사 표면도, 표면 법선을 뒤섞는 인스턴스도 다루지 않았습니다. 하지만 이 설계는 전반적으로 깔끔하며, 이런 문제들은 모두 해결할 수 있습니다. 지금은 `specular`만 고치겠습니다. 금속과 유전체 재질은 쉽게 고칠 수 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class metal : public material {
public:
metal(const color& albedo, double fuzz) : albedo(albedo), fuzz(fuzz < 1 ? fuzz : 1) {}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(const ray& r_in, const hit_record& rec, scatter_record& srec) const override {
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
vec3 reflected = reflect(r_in.direction(), rec.normal);
reflected = unit_vector(reflected) + (fuzz * random_unit_vector());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
srec.attenuation = albedo;
srec.pdf_ptr = nullptr;
srec.skip_pdf = true;
srec.skip_pdf_ray = ray(rec.p, reflected, r_in.time());
return true;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
}
private:
color albedo;
double fuzz;
};
class dielectric : public material {
public:
dielectric(double refraction_index) : refraction_index(refraction_index) {}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
bool scatter(const ray& r_in, const hit_record& rec, scatter_record& srec) const override {
srec.attenuation = color(1.0, 1.0, 1.0);
srec.pdf_ptr = nullptr;
srec.skip_pdf = true;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
double ri = rec.front_face ? (1.0/refraction_index) : refraction_index;
vec3 unit_direction = unit_vector(r_in.direction());
double cos_theta = std::fmin(dot(-unit_direction, rec.normal), 1.0);
double sin_theta = std::sqrt(1.0 - cos_theta*cos_theta);
bool cannot_refract = ri * sin_theta > 1.0;
vec3 direction;
if (cannot_refract || reflectance(cos_theta, ri) > random_double())
direction = reflect(unit_direction, rec.normal);
else
direction = refract(unit_direction, rec.normal, ri);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
srec.skip_pdf_ray = ray(rec.p, direction, r_in.time());
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
return true;
}
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [material-scatter]: [material.h] 금속과 유전체의 산란 메서드]
흐릿함(fuzziness)이 0이 아니면 이 표면은 실제로 이상적인 정반사 표면이 아니지만, 암시적 샘플링은 이전과 똑같이 작동한다는 점에 유의하십시오. 사실상 정반사로 취급하는 재질에 대해서는 PDF 작업을 모두 건너뛰는 셈입니다.
`ray_color()`에는 암시적으로 샘플링된 광선을 생성하는 새로운 경우 하나만 추가하면 됩니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class camera {
...
private:
...
color ray_color(const ray& r, int depth, const hittable& world, const hittable& lights)
const {
...
if (!rec.mat->scatter(r, rec, srec))
return color_from_emission;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
if (srec.skip_pdf) {
return srec.attenuation * ray_color(srec.skip_pdf_ray, depth-1, world, lights);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
auto light_ptr = make_shared(lights, rec.p);
mixture_pdf p(light_ptr, srec.pdf_ptr);
...
}
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [ray-color-implicit]: [camera.h] 암시적으로 샘플링된 광선을 사용하는 광선 색상 함수]
블록 하나를 금속으로 바꾸어 작업 결과를 확인해 보겠습니다. 블록 중 하나를 유리 물체로 바꾸고 싶기도 하지만, 그것은 다음 절로 미루겠습니다. 유리 물체는 잘 렌더링하기 어렵기 때문에 유리 물체용 PDF를 만들고 싶지만, 그렇게 하려면 먼저 해야 할 작업이 더 있습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
// Light
world.add(make_shared(point3(213,554,227), vec3(130,0,0), vec3(0,0,105), light));
// Box 1
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
shared_ptr aluminum = make_shared(color(0.8, 0.85, 0.88), 0.0);
shared_ptr box1 = box(point3(0,0,0), point3(165,330,165), aluminum);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
box1 = make_shared(box1, 15);
box1 = make_shared(box1, vec3(265,0,295));
world.add(box1);
// Box 2
shared_ptr box2 = box(point3(0,0,0), point3(165,165,165), white);
box2 = make_shared(box2, -18);
box2 = make_shared(box2, vec3(130,0,65));
world.add(box2);
// Light Sources
auto empty_material = shared_ptr();
quad lights(point3(343,554,332), vec3(-130,0,0), vec3(0,0,-105), empty_material);
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [scene-cornell-al]: [main.cc] 알루미늄 재질을 사용한 코넬 박스 장면]
상자를 향하는 방향이 더 높은 밀도로 샘플링되지 않기 때문에, 결과 이미지에서는 천장에 노이즈가 많은 반사가 나타납니다.
구 물체 샘플링하기
-------------------------
금속 블록에 대한 PDF를 만들면 천장의 노이즈를 줄일 수 있습니다. 블록을 유리로 만든다면 그 블록에 대한 PDF도 필요할 것입니다. 하지만 블록에 대한 PDF를 만드는 것은 상당히 많은 작업이 필요하고 그다지 흥미롭지도 않으므로, 대신 유리 구에 대한 PDF를 만들어 보겠습니다. 이쪽이 더 빠르고 더 흥미로운 렌더링 결과를 만들어 냅니다. 적절한 PDF 분포를 결정하려면 구를 샘플링하는 방법을 알아내야 합니다. 구 바깥의 한 점에서 구를 샘플링하려는 경우, 표면 위의 무작위 점을 하나 고르는 것으로 끝낼 수는 없습니다. 그렇게 하면 구의 반대편에 있는 점을 자주 고르게 되는데, 그 점은 구의 앞면에 가려집니다. 임의의 점에서 보이는 구의 면을 균일하게 샘플링하는 방법이 필요합니다. 구 바깥의 한 점에서 구의 입체각을 균일하게 샘플링하는 것은 사실 원뿔을 균일하게 샘플링하는 것에 불과합니다. 원뿔의 축은 광선의 원점에서 구의 중심을 지나며, 원뿔의 옆면은 구에 접합니다 -- 아래 그림을 참조하십시오. 코드에 `theta_max`가 있다고 합시다. 무작위 방향 생성하기 장에서 $\theta$를 샘플링하려면 다음 식이 성립했음을 떠올려 보십시오:
$$ r_2 = \int_{0}^{\theta} 2 \pi f(\theta') \sin(\theta') d\theta' $$
여기서 $f(\theta')$는 아직 계산하지 않은 상수 $C$이므로:
$$ r_2 = \int_{0}^{\theta} 2 \pi C \sin(\theta') d\theta' $$
미적분으로 풀면:
$$ r_2 = 2\pi \cdot C \cdot (1-\cos(\theta)) $$
따라서
$$ \cos(\theta) = 1 - \frac{r_2}{2 \pi \cdot C} $$
무작위 방향이 $\theta_{max}$보다 작아야 하도록 분포를 제한하고 있습니다. 이는 0부터 $\theta_{max}$까지의 적분이 1이어야 함을 의미하며, 따라서 $r_2 = 1$입니다. 이를 이용하여 $C$를 구할 수 있습니다:
$$ r_2 = 2\pi \cdot C \cdot (1-\cos(\theta)) $$
$$ 1 = 2\pi \cdot C \cdot (1-\cos(\theta_{max})) $$
$$ C = \frac{1}{2\pi \cdot (1-\cos(\theta_{max}))} $$
이로써 $\theta$, $\theta_{max}$, $r_2$ 사이의 등식을 얻습니다:
$$ \cos(\theta) = 1 + r_2 \cdot (\cos(\theta_{max})-1) $$
$\phi$는 이전과 같이 샘플링하므로 다음과 같습니다:
$$ z = \cos(\theta) = 1 + r_2 \cdot (\cos(\theta_{max}) - 1) $$
$$ x = \cos(\phi) \cdot \sin(\theta) = \cos(2\pi \cdot r_1) \cdot \sqrt{1-z^2} $$
$$ y = \sin(\phi) \cdot \sin(\theta) = \sin(2\pi \cdot r_1) \cdot \sqrt{1-z^2} $$
그렇다면 $\theta_{max}$는 무엇입니까?
![Figure [sphere-enclosing-cone]: 구를 감싸는 원뿔
](../images/fig-3.12-sphere-enclosing-cone.jpg)
그림에서 $\sin(\theta_{max}) = R / length(\mathbf{c} - \mathbf{p})$임을 알 수 있습니다. 따라서:
$$ \cos(\theta_{max}) = \sqrt{1 - \frac{R^2}{length^2(\mathbf{c} - \mathbf{p})}} $$
방향에 대한 PDF도 계산해야 합니다. 구를 향한 균일 분포의 경우 PDF는 $1/\text{solid_angle}$입니다. 구의 입체각은 얼마입니까? 이는 위의 $C$와 관련이 있습니다. 입체각은 정의상 단위 구 위의 면적이므로, 적분은 다음과 같습니다
$$ \text{solid_angle} = \int_{0}^{2\pi} \int_{0}^{\theta_{max}} \sin(\theta)
= 2 \pi \cdot (1-\cos(\theta_{max})) $$
이러한 모든 계산은 수식을 검토해 보는 것이 좋습니다. 저는 보통 극단적인 경우를 대입해 봅니다(이 개념을 알려 주신 제 고등학교 물리 선생님 Mr. Horton께 감사드립니다). 반지름이 0인 구의 경우 $\cos(\theta_{max}) = 1$이며, 이는 성립합니다. $\mathbf{p}$에서 접하는 구의 경우 $\cos(\theta_{max}) = 0$이고 $2\pi$는 반구의 면적이므로, 이 역시 성립합니다.
구 코드 업데이트하기
-------------------------
구 클래스에는 PDF 관련 함수 두 개가 필요합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#include "hittable.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
#include "onb.h"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class sphere : public hittable {
public:
...
aabb bounding_box() const override { return bbox; }
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value(const point3& origin, const vec3& direction) const override {
// This method only works for stationary spheres.
hit_record rec;
if (!this->hit(ray(origin, direction), interval(0.001, infinity), rec))
return 0;
auto dist_squared = (center.at(0) - origin).length_squared();
auto cos_theta_max = std::sqrt(1 - radius*radius/dist_squared);
auto solid_angle = 2*pi*(1-cos_theta_max);
return 1 / solid_angle;
}
vec3 random(const point3& origin) const override {
vec3 direction = center.at(0) - origin;
auto distance_squared = direction.length_squared();
onb uvw(direction);
return uvw.transform(random_to_sphere(radius, distance_squared));
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
private:
...
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
static vec3 random_to_sphere(double radius, double distance_squared) {
auto r1 = random_double();
auto r2 = random_double();
auto z = 1 + r2*(std::sqrt(1-radius*radius/distance_squared) - 1);
auto phi = 2*pi*r1;
auto x = std::cos(phi) * std::sqrt(1-z*z);
auto y = std::sin(phi) * std::sqrt(1-z*z);
return vec3(x, y, z);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [sphere-pdf]: [sphere.h] PDF를 갖춘 구]
먼저 광원 대신 구만 샘플링해 볼 수 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
// Light
world.add(make_shared(point3(213,554,227), vec3(130,0,0), vec3(0,0,105), light));
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
// Box
shared_ptr box1 = box(point3(0,0,0), point3(165,330,165), white);
box1 = make_shared(box1, 15);
box1 = make_shared(box1, vec3(265,0,295));
world.add(box1);
// Glass Sphere
auto glass = make_shared(1.5);
world.add(make_shared(point3(190,90,190), 90, glass));
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
// Light Sources
auto empty_material = shared_ptr();
quad lights(point3(343,554,332), vec3(-130,0,0), vec3(0,0,-105), empty_material);
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [sampling-sphere]: [main.cc] 구만 샘플링하기]
이렇게 하면 노이즈가 많은 방이 나오지만, 구 아래의 커스틱은 좋습니다. 제 코드에서는 광원을 샘플링할 때보다 다섯 배의 시간이 걸렸습니다. 이는 아마도 유리에 맞는 광선들의 계산 비용이 크기 때문일 것입니다!
Hittable 리스트에 PDF 함수 추가하기
---------------------------------------
아마도 구와 광원을 모두 샘플링하는 것이 좋을 것입니다. 두 분포의 혼합 밀도를 만들어 이를 수행할 수 있습니다. `ray_color()` 함수에 hittable 리스트를 전달하여 혼합 PDF를 구성하는 방식으로 할 수도 있고, `hittable_list`에 PDF 함수를 추가할 수도 있습니다. 두 방법 모두 잘 작동하리라 생각하지만, 저는 `hittable_list`에 기능을 추가하는 방법을 택하겠습니다.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class hittable_list : public hittable {
public:
...
aabb bounding_box() const override { return bbox; }
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
double pdf_value(const point3& origin, const vec3& direction) const override {
auto weight = 1.0 / objects.size();
auto sum = 0.0;
for (const auto& object : objects)
sum += weight * object->pdf_value(origin, direction);
return sum;
}
vec3 random(const point3& origin) const override {
auto int_size = int(objects.size());
return objects[random_int(0, int_size-1)]->random(origin);
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [density-mixture]: [hittable_list.h] 밀도 혼합 만들기]
`camera::render()`에 전달할 광원 목록을 구성합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
int main() {
...
// Light Sources
auto empty_material = shared_ptr();
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
hittable_list lights;
lights.add(
make_shared(point3(343,554,332), vec3(-130,0,0), vec3(0,0,-105), empty_material));
lights.add(make_shared(point3(190, 90, 190), 90, empty_material));
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
...
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [scene-density-mixture]: [main.cc] 장면 업데이트하기]
그러면 이전과 마찬가지로 1000개의 샘플로 괜찮은 이미지를 얻습니다:
표면 여드름 처리하기
----------------------
한 예리한 독자가 위 이미지에 검은 반점이 몇 개 있다고 지적했습니다. 모든 몬테카를로 레이 트레이서는 다음과 같은 메인 루프를 가지고 있습니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pixel_color = average(many many samples)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
렌더링 결과에 어떤 형태의 여드름이 생기고 이 여드름이 흰색이나 검은색이라면 -- 하나의 "나쁜" 샘플이 픽셀 전체를 망가뜨리는 것처럼 보인다면 -- 그 샘플은 아마도 엄청나게 큰 수이거나 `NaN`(Not A Number)일 것입니다. 이 특정한 여드름은 아마도 `NaN`일 것입니다. 제 경우에는 대략 1,000만~1억 개의 광선마다 한 번 정도 발생하는 것 같습니다.
그러니 큰 결정을 내려야 합니다: 이 버그를 덮어 두고 `NaN`을 검사할 것인지, 아니면 그냥 `NaN`을 없애 버리고 나중에 이것이 다시 문제를 일으키지 않기를 바랄 것인지입니다. 저는 언제나 게으른 전략을 택합니다. 특히 부동 소수점을 다루는 일이 어렵다는 것을 알고 있을 때는 더욱 그렇습니다. 먼저, `NaN`을 어떻게 검사할까요? `NaN`에 대해 제가 항상 기억하는 한 가지는 `NaN`은 자기 자신과 같지 않다는 것입니다. 이 요령을 사용하여, 모든 `NaN` 성분을 0으로 대체하도록 `write_color()` 함수를 업데이트합니다:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
void write_color(std::ostream& out, const color& pixel_color) {
auto r = pixel_color.x();
auto g = pixel_color.y();
auto b = pixel_color.z();
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
// Replace NaN components with zero.
if (r != r) r = 0.0;
if (g != g) g = 0.0;
if (b != b) b = 0.0;
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
// Apply a linear to gamma transform for gamma 2
r = linear_to_gamma(r);
g = linear_to_gamma(g);
b = linear_to_gamma(b);
// Translate the [0,1] component values to the byte range [0,255].
static const interval intensity(0.000, 0.999);
int rbyte = int(256 * intensity.clamp(r));
int gbyte = int(256 * intensity.clamp(g));
int bbyte = int(256 * intensity.clamp(b));
// Write out the pixel color components.
out << rbyte << ' ' << gbyte << ' ' << bbyte << '\n';
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [write-color-nan]: [color.h] NaN을 허용하는 write_color 함수]
다행히 검은 반점이 사라졌습니다:
남은 인생
====================================================================================================
이 책의 목적은 물리 기반 렌더러의 샘플링 방식을 구성할 때 필요한 모든 사소한 세부 사항을 (하나하나 빠짐없이 꼼꼼하게) 짚어 보는 것이었습니다. 이제 이 모든 세부 사항을 바탕으로 다양한 잠재적 경로를 탐구할 수 있을 것입니다.
몬테카를로 방법을 탐구하고 싶다면, 메트로폴리스와 같은 양방향 및 경로 공간 접근법을 살펴보십시오. 여러분의 확률 공간은 입체각에 대한 것이 아니라 경로 공간에 대한 것이 될 것이며, 여기서 경로는 고차원 공간의 다차원 점입니다. 그렇다고 겁먹지 마십시오 -- 어떤 대상을 숫자 배열로 기술할 수 있다면, 수학자들은 그것을 그러한 점들로 이루어진 가능한 모든 배열의 공간 속 한 점이라고 부릅니다. 이는 단지 보여 주기 위한 것이 아닙니다. 일단 그런 깔끔한 추상화를 얻고 나면, 여러분의 코드도 깔끔해질 수 있습니다. 깔끔한 추상화야말로 프로그래밍의 핵심입니다!
영화용 렌더러를 만들고 싶다면, 여러 스튜디오와 Solid Angle에서 나온 논문을 살펴보십시오. 그들은 자신들의 기술에 대해 놀라울 만큼 개방적입니다.
고성능 레이 트레이싱을 하고 싶다면, 먼저 Intel과 NVIDIA의 논문을 살펴보십시오. 그들 역시 놀라울 만큼 개방적입니다.
하드코어한 물리 기반 렌더러를 만들고 싶다면 렌더러를 RGB 방식에서 스펙트럼 방식으로 전환하십시오. 저는 각 광선이 무작위 파장을 갖고, 프로그램에 있는 거의 모든 RGB가 float로 바뀌는 방식을 무척 좋아합니다. 비효율적으로 들리지만 실제로는 그렇지 않습니다!
어떤 방향을 택하든 광택 BRDF 모델을 추가하십시오. 선택할 수 있는 모델이 많으며, 모델마다 장점이 있습니다.
즐겁게 작업하십시오!
[Peter Shirley][] 솔트레이크시티, 2016년 3월
(insert acknowledgments.md.html here)
이 책 인용하기
====================================================================================================
인용 형식을 통일하면 이 저작물의 출처, 위치, 버전을 더 쉽게 확인할 수 있습니다. 이 책을 인용하실 때는 가능하면 다음 형식 중 하나를 사용해 주시기 바랍니다.
기본 정보
-----------
- **제목(시리즈)**: “주말 동안 배우는 레이 트레이싱 시리즈”
- **제목(도서)**: “레이 트레이싱: 남은 인생”
- **저자**: Peter Shirley, Trevor David Black, Steve Hollasch
- **버전/판**: (v4.0.2)
- **날짜**: 2025-04-25
- **URL(시리즈)**:
- **URL(도서)**:
스니펫
---------
### Markdown
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[_Ray Tracing: The Rest of Your Life_](https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### HTML
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Ray Tracing: The Rest of Your Life
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### LaTeX 및 BibTex
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~\cite{Shirley2025}
@misc{Shirley2025,
title = {Ray Tracing: The Rest of Your Life},
author = {Peter Shirley, Trevor David Black, Steve Hollasch},
year = {2025},
month = {April},
note = {\small \texttt{https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html}},
url = {https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html}
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### BibLaTeX
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\usepackage{biblatex}
~\cite{Shirley2025}
@online{Shirley2025,
title = {Ray Tracing: The Rest of Your Life},
author = {Peter Shirley, Trevor David Black, Steve Hollasch},
year = {2025},
month = {April},
url = {https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html}
}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### IEEE
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
“Ray Tracing: The Rest of Your Life.”
raytracing.github.io/books/RayTracingTheRestOfYourLife.html
(accessed MMM. DD, YYYY)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### MLA:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Ray Tracing: The Rest of Your Life. raytracing.github.io/books/RayTracingTheRestOfYourLife.html
Accessed DD MMM. YYYY.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Peter Shirley]: https://github.com/petershirley
[Steve Hollasch]: https://github.com/hollasch
[Trevor David Black]: https://github.com/trevordblack
[readme]: ../README.md
[releases]: https://github.com/RayTracing/raytracing.github.io/releases/
[wiki-further]: https://github.com/RayTracing/raytracing.github.io/wiki/Further-Readings