블랙홀 시뮬레이션 - 2. 구현

1편

블랙홀 시뮬레이션 - 1. 일반상대론 설명

그럼 이제 어떻게 구현했는지 알아보자.

 

레이 마칭

https://rito15.github.io/posts/ray-marching/

레이 마칭은 카메라 앞에 스크린이 있다고 가정할 때, 카메라와 스크린 위의 점을 잇는 광선을 만들어서 물체에 부딫힐 때 까지 전진시키는 방법임. 게임에서 접할 수 있는 레이트레이싱과 좀 비슷한 기술인데, 레이 트레이싱은 광선-물체의 부딫힘을 이용하지만, 레이 마칭은 주변 물체와의 거리를 이용해서 광선을 전진시키며 물체와 부딫힐 때 까지 조금씩 전진시키는 방법임.

https://en.wikipedia.org/wiki/Ray_marching(좌), https://michaelwalczyk.com/blog-ray-marching.html(우)

저렇게 써 놓으면 아무도 이해 못하니까 위에 그림을 보자.

 

 

 

 

 

 

 

 

 

 

이건 간단한 구를 그리는 셰이더인데, 한 번 실행해 보면...

구가 잘 나옴.

 

미분 레이 트레이싱

그런데 블랙홀 시뮬레이션에서는 빛이 직선으로 이동하지 않고 미분방정식을 따라 이동함. 

\[
F(r) = \frac{d^2r}{d\tau^2} =  - \frac{3}{2}\frac{h^2}{r^4}
\]

저번에 구한 이 식을 따라 이동함.

 

이 레이 마칭을 위해서는 약간 다른 레이 마칭 기법을 사용해야 함. 빛이 지나가는 경로가 저 미분 방정식을 따라야 함.

vec4 rayMarch(vec3 pos, vec3 dir) {
    vec4 color = vec4(0.0);
    float alpha = 1.0;

    vec3 h = cross(pos, dir);
    float h2 = dot(h, h);

    for (int i = 0; i < iteration; ++i){
        verlet(pos, h2, dir, steps);

        if (dot(pos, pos) < 1.0){
            return color;
        }

        if (disk)
            diskRender(pos, color, alpha, dir);
    }

    vec4 skyColor = texture(cubemap, normalize(dir));
    color += skyColor;

    return color;
}

그럼 일단 레이 마칭 함수부터 보자. 위의 일반적인 레이 마칭 코드와 아직까지 크게 다른 부분은 없다.

vec3 acceleration(float h2, vec3 pos){
    float r2 = dot(pos, pos);
    float r5 = pow(r2, 2.5);
    vec3 acc = -1.5 * h2 * pos / r5;
    return acc;
}

void verlet(inout vec3 pos, float h2, inout vec3 dir, float dt){
    vec3 a = acceleration(h2, pos);
    vec3 pos_new = pos + dir * dt + 0.5 * a * dt * dt;
    vec3 a_new = acceleration(h2, pos_new);

    dir += 0.5 * (a + a_new) * dt;
    pos = pos_new;
}

생각보다 간?단해 보이지?

 

이 코드에 대해 설명하자면 다음과 같음:

1. 현재 카메라의 위치와 보는 방향을 전달받음

2. 각운동량 h를 구함

3. 반복 횟수만큼

  3-1. 빛의 지오데식 방정식을 verlet방법으로 풀고 빛을 천천히 전진시킴

  3-2. 만약 사건의 지평선 구(반지름이 1이고 중심이 0인 원)에 있다면 빛이 튀어나오지 못함.

  3-3. 디스크 렌더링

4. 배경 렌더링

 

매우 간단해 보이고, 실제로 간단하지만, 이건 실제로 쓰이는 코드는 아니고 프로토타입 코드라 중요하지 않은 부분을 생략했음. 그래도 이 정도만 있어도 블랙홀 시뮬레이션은 작동함.

 

강착 원반 렌더링

블랙홀 시뮬레이션은 의외로 쉽고 짧게 만들 수 있지만, 블랙홀을 더욱 사실적으로 만드는 것은 어려움. 흑체 복사라던가, ISCO라던가 등등을 반영해야 함. 대부분은 간단한 값을 넣는 것만으로 해결 가능하지만, 몇 가지 어려운 것들이 있음.

 

흑체 복사

흑체 복사는 흑체(=모든 전자기파를 흡수하는 물체)에서 나오는 열 전자기 복사임. 흑체 복사는 온도에 따른 특정한 연속적인 파장을 갖고 있음. 블랙홀 주변의 강착 원반을 이루는 물질을 흑체라고 생각하면 편하므로, 여기서는 흑체 복사를 사용하겠음.

 

위 그래프는 1000K부터 10000K까지의 온도에 대한 색을 그린 그래프임. 이 그래프를 어떻게 그렸는지 알아보자.

 

자세한 유도는 생략하고, 흑체 복사에서 spectral radiance는 다음과 같음.

\[
I(\lambda, T) = \frac{2\pi h c^2}{\lambda^5} \frac{1}{e^{\frac{h c}{\lambda k_B T}} - 1}
\]

이 식을 각 파장에 대해 적분하면 각 파장 별 복사 강도가 나오고, 각 파장에 대한 RGB 색상을 구하면, 파장 별 복사 색상을 얻을 수 있음. 둘을 이용해서 위 그래프처럼 각 복사 파장에 따른 복사 강도로 정규화된 색상을 구할 수 있음. 처음 내가 이 부분에서 실수를 했는데 복사 강도를 고려하지 않고 색상만을 이용해서 흑체 복사를 얻었다가 잘못된 렌더링 이미지를 얻었음.

 

위 그래프는 

 

 

 

다음과 같은 코드로 얻을 수 있음.

 

위  코드에서 get_matching_functinos는 CIE 1931색 공간에서의 각 파장별 XYZ값을 배열로 나타낸 것임. 여기서

X는 빨강-초록 색 자극

Y는 밝기

Z는 파랑-보라 색 자극

정도로 이해하면 될 듯함. 정확히는 사람 눈의 삼자극 값임

각 XYZ값은 다음과 같은 값을 가짐. 이를 sRGB 색 공간으로 변환해야 하는데, 행렬을 통해 변환이 가능함.

$\begin{bmatrix}
R \\
G \\
B
\end{bmatrix}
=
\begin{bmatrix}
  3.2406 & -1.5372 & -0.4986 \\
 -0.9689 &  1.8758 &  0.0415 \\
  0.0557 & -0.2040 &  1.0570
\end{bmatrix}
\begin{bmatrix}
X \\
Y \\
Z
\end{bmatrix}$

이렇게 변환할 수 있음.(DCI-P3같은 다른 색 공간도 테스트 중임)

 

문제는 이렇게 나온 값의 밝기는 10의 6~17제곱 정도로 큰 크기를 가져서 표현하기 어렵다는 거임. 그래서 밝기를 로그 크기로 정규화해 줘서 0~1 값을 갖도록 해 줄것임.(get_blackbody_cY)

 

정리하자면, 인간이 볼 수 있는 가시광선 대역인 380nm~780nm 에 대해 각 파장별로 흑체 복사의 복사 강도를 구하고, 파장에 맞는 CIE XYZ  색 공간 좌표를 얻음. 그 다음 XYZ를 sRGB로 변환하면 화면에 표시할 수 있는 RGB값을 얻음.

 

상대론적 도플러 효과

도플러 효과는 관측자로부터 다가오거나 멀어지는 광원에서 나오는 빛의 파장이 짧아지거나 늘어나는 현상임. 강착 원반에 있는 물질은 상대론적 속도로 회전하기 때문에 이 효과를 사용해야 더 현실적인 블랙홀을 만들 수 있음.

 

이 효과를 적용하면 위의 그림처럼 강착 원반의 한 쪽이 더 짧은 파장으로 보이고, relativistic beaming으로 짧은 파장 쪽이 더 밝게 보임.

 

이 효과도 어렵지 않게 적용 가능한데, 다음의 식을 이용하면 됨:

\[
\frac{\nu}{\nu^\prime}= \frac{g(k, l)}{g(k^\prime, l^\prime)}
\]

이렇게 구해도 되지만, 좀 더 쉽게 하고 싶다면

\[
1 + z = \frac{1}{(1 - \beta^2)^{1/2}}
\]

도플러 효과 식과

\[
1 + z = \frac{1}{(1 - 1 / r)^{1/2}}
\]

중력 적색편이 식을 사용한 후 두 효과를 합치면 됨.

 

슈테판-볼츠만 법칙에 의하면 열복사의 에너지 밀도는 온도의 4제곱에 비례하므로 도플러 효과 식을 구했다면, 해당 지점에서의 색에 도플러 효과항의 4제곱만큼을 원반의 색에서 나누어 주면 됨. 그러면 원반의 한쪽 부분만이 밝게 나오는 relativistic beaming을 구현할 수 있음.

cl = getBlackbodyColor(accretionTempMod, doppler * redshift);// Normalized Blackbody color
cl /= pow(doppler * redshift, 4);//Stefan-Boltzmann Law

 

ISCO

Innermost stable circular orbit(ISCO)는 일반상대론에서 시험 입자가 거대한 물체 주변에서 안정적으로 궤도를 돌 수 있는 가장 작은 반지름의 원궤도임.

 

자세한 유도는 생략하고, 슈바르츠실트 블랙홀에서 ISCO는 슈바르츠실트 반지름의 3배

\[
r_{ms}=6\frac{GM}{c^2}=3R_s
\]

이고 불안정하더라도 궤도를 돌 수 있는 반지름 한계는 

\[
r_{mb} = 4\frac{GM}{c^2}=2R_s
\]

이므로, noise를 사용하여 volumetric cloud로 강착 원반을 만들었을 때, 3r에서 2r까지 밀도를 줄여가는 방식으로 구현하면 됨.

 

photon sphere 같은 경우 블랙홀을 제대로 구현했다면, 슈바르츠실트 반지름의 1.5배 지점에 얇은 원형 구가 형성됨을 확인 가능함. 

 

Simplex Noise와 Volumetric cloud

강착 원반은 입자들의 모음으로 이루어져 있다고 모델링할 수 있는데, 입자들을 전부 궤도역학, 유체역학을 사용하여 시뮬레이션 하기에는 이런 GPU를 쓰지 않는 이상 컴퓨팅 자원이 부족할 것이므로 간단히 궤도 속도를 따라 움직이는 구름 형태로 구현해야 함.

 

구름을 만들기 위해서는 노이즈 함수를 사용해야 하는데, perlin noise나 simplex noise같은 연속적인 노이즈 함수를 사용하면 자연스러운 구름을 만들 수 있음. 이 구름을 레이 마칭을 이용하여 빛이 구름 위를 지나갈 때 색을 더해주면  다음과 같은 성운이나 구름을 얻을 수 있음.

 

먼저 강착 원반의 밀도 프로파일을 만들고

(https://jila.colorado.edu/~pja/astr3830/lecture26.pdf)

-> z 방향으로는 $\rho = \rho_{z = 0} e ^ {-z^2 / 2 H^2}$

(https://en.wikipedia.org/wiki/Accretion_disk)

-> r방향으로는 $\rho \propto \rho_{r = 0} r ^{-15/8}$

 

해당 밀도 프로파일을 이용해서 블랙홀 주변에 적절히 만들어주면 됨.

 

HDR 렌더링, 블룸, 톤매핑, 감마 보정 등의 후처리 프로세싱

https://learnopengl.com/Advanced-Lighting/HDR

LearnOpenGL - HDR

https://learnopengl.com/Advanced-Lighting/Bloom

LearnOpenGL - Bloom

렌더링 결과를 더 현실감 있게 만들어줌.

 

렌더링 시 색상 밝기 값은 0.0~1.0사이인데, HDR 렌더링을 적용하면 그 이상의 색상도 렌더링할 수 있음. 따라서 밝고 어두움을 더 정확히 나타낼 수 있게 됨. 하지만, 대부분의 사용자 모니터는 SDR만을 지원하기 때문에, 톤매핑을 해줘야 함. 블랙홀 렌더링 프로젝트는, 매우 밝은 부분과 어두운 부분이 공존하므로, HDR렌더링 후 톤매핑을 해주는 것이 좋음.

 

다양한 톤매핑 함수를 사용 가능하지만, 블랙홀 렌더링 프로젝트에는 Hable이 가장 잘 어울리는 방식임.

 

밝은 부분을 보정하고 번짐 효과를 구현하기 위해 감마 보정과 블룸을 구현했음. 감마 보정의 경우 DCI-P3에서 많이 사용하는 2.4를 사용했고, 블룸은 10단계를 거치도록 했음. 

 

코드 최적화 방법

프래그먼트 셰이더 대신 컴퓨트 셰이더를 사용하여  캐시 적중율을 높이려는 시도를 해 봤음:
https://computergraphics.stackexchange.com/questions/9956/performance-of-compute-shaders-vs-fragment-shaders-for-deferred-rendering

Performance of Compute Shaders vs. Fragment Shaders for Deferred Rendering

실질적인 효과가 있었냐 하면 몰?루 블랙홀 렌더링에만 적용해 보았는데, 평균 1프레임 정도 향상이 있었음. 공유 메모리를 사용해서 속도 향상을 얻을만한 부분이 적어(배경 렌더링 같은 부분 정도?) 성능 향상이 적은 것 같음. 이건 잘 체감하려면 격자 형태의 시뮬레이션이나 행렬 연산이 좋지 레이 트레이싱 시뮬레이션에서는 별로임.

 

if문의 위치를 잘 조절해서 분기예측 성능을 높이기도 적용해 봤음. 블랙홀 주변에 바운딩 볼륨을 쳐서 꼭 필요한 부분만 렌더링하게 하고, if문 분기를 렌더링 부문과 설정 부분 모두 최대한 간소화했음.

 

컴퓨트 셰이더 타일 수는 16 X 16으로 설정함.

https://frguthmann.github.io/posts/compute_shaders/

Exploring ways to optimize compute shaders - Part 1.

nSight graphics로 디버깅 중인데, 좀 더 개선할 부분을 찾는 중임:

Warp Occupancy가 60% 정도로 낮은 원인은 레지스터 부족으로 발생한 건데, OpenGL 특성상 셰이더 디버깅이 좀 많이 까다로워서 어디서 막히는지는 잘 모르겠음. 일단 셰이더 내의 변수를 줄이는 등 최적화를 하고 있기는 한데, 개선이 잘 안됨. 텍스쳐를 빼거나 원반을 끄고 렌더링 해 봐도 프레임 레이트만 증가하고 비슷한 걸 보면 레이트레이싱 문제인데, 디버깅이 안 되니까 어디서 문제인지를 모르겠음.Vulkan을 배울 이유가 또 늘어났다

 

최종 결과물

이렇게 나옴.

 

블랙홀을 정면에서 보거나 위에서 보면 동심원 모양의 무늬가 보이는데, 그건 부동소수점 정밀도 오류로 보여짐. RKF45 등을 사용해 보아도 개선되지 않아서 그냥 놔둠.

 

참고 자료:

https://arxiv.org/pdf/2010.08735

위 사진과 거의 비슷한 결과물을 얻을 수 있음. 그런데 이 자료에 나와있는 코드는 좀 알아먹기 어렵게 만들어져 있더라

 

코드는 계속 업데이트 중이니 앞으로도 조금씩 바뀔 수 있음.

코드: 

https://github.com/hydrogendeuteride/BlackHoleRayTracer

GitHub - hydrogendeuteride/BlackHoleRayTracer: blackhole ray marching

실행 가능 파일:

https://github.com/hydrogendeuteride/BlackHoleRayTracer/releases

Releases · hydrogendeuteride/BlackHoleRayTracer

윈도우와 리눅스는 작동 확인 했는데, 맥은 잘 모르겠음.

되겠지 뭐

이거 만드느라 대학교 마지막 방학 한달 반을 날리긴 했는데(겨울 방학? 한번 더 남음), 이거 말고도 레이 트레이싱 공부, CUDA 공부, 웹 프론트-백엔드 공부도  해놨고 이번 학기가 마지막 학기라 시간이 많으므로 프로젝트 좀 할 수도 있을듯함.취업 준비는 안하고 이상한것만 함.

 

다음에 만들 프로젝트 목록:

이징 모델 GPU 시뮬레이션: CUDA 학습용인데 아마 일주일이면 할듯

GPU 옥트리: 하다가 잠시 멈춤

웹서버, 페이지 만들어서 webGPU로 뭔가 시각화하고 놀기