[코딩의탑] 3층: 바다 렌더링
바다를 그려봅시다.
수많은 게임들은 바다를 배경으로 하고 있습니다. 굳이 바다를 게임의 메인 배경으로 하지 않더라도, 게임 어느 중간에는 꼭 바다를 볼 수 있는 순간이 나오기 마련이죠.
 |
| <Sea of Thieves> |
이러한 바다를 유체역학적 지식으로 구현하는 것은 쉽지 않습니다. 화면에 보이는 입자가 너무 많아 현실적으로 계산 능력이 부족하기 때문입니다. 따라서 많은 게임들은 바다를 구현하기 위해서 수많은 편법을 활용하고 있습니다. 대표적인 방법은 Displacement Mapping(변위 매핑)을 이용하는 것이죠.
Displacement Mapping
Displacement Mapping이란 Displacement Map을 이용해 물체를 구성하는 일련의 정점들을 재배치하는 것을 의미합니다. 가령, 평평한 격자의 물체가 있고 다음과 같은 두 개의 Displacement Map이 있다면 우리는 파도와 같은 현상을 만들어낼 수 있습니다. 하나는 큰 파도(저주파)를, 다른 하나는 작은 파도(고주파)를 재현합니다. 흰색일수록 위에 있으며, 검은색일수록 아래에 있게 됩니다.
이러한 Displacement Map들을 겹친 후에, 서로 다른 방향으로 움직이게 만들고 겹쳐진 값을 이용해서 Displacement Mapping을 한다면 저렴한 연산으로 파도와 같은 환상을 만들어낼 수 있습니다. 하지만 이 방법은 정적으로 파도가 만들어진다는 한계가 있습니다. 바람의 방향, 파도의 세기 등을 조절하기가 어렵죠. 조금 더 멋진 Displacement Map을 만들 수 있는 방법이 있을까요?
Jerry Tessendorf의 'Simulating Ocean Water'
영화 프로듀서이자 교수인 Tessendorf는 그럴듯한 바다를 시뮬레이션할 수 있는 방법을 논문으로 정리해 제안하였습니다. 제가 다룰 내용을 간단히 요약하면 다음과 같습니다: 파도는 일종의 사인파와 코사인파의 합입니다. 따라서 푸리에 변환의 철학을 이용해서 복잡한 주기적 파도를 생성할 수 있으며, 추가적으로 의미있는 정보들을 활용해 그럴듯한 파도를 생성할 수 있습니다.
푸리에 변환
파도는 바다 위에서 나타나는 파동의 현상입니다. 파동은 사인파와 코사인파를 이용해서 나타낼 수 있습니다. 따라서 파도를 구성하는 각 파동의 진동수와 진폭을 안다면, 시간의 함수를 이용해서 언제든지 파도를 나타낼 수 있을 것입니다. 이 글에서는 푸리에 변환에 대해서 자세히 다루지 않음에 주의하세요. 우리는 다음과 같은 단계를 거칠 것입니다:
- 파도에 대한 함수를 구성합니다. 이 함수는 시간과 위치(주파수 영역이므로 진동수가 될 것입니다)를 변수로 갖습니다.
- 매 프레임에 대해 3~6을 수행합니다:
- 시간과 진동수에 대해 앞서 구성한 함수를 통해 주파수 영역의 결과를 생성합니다.
- 추가로, Normal Mapping을 위한 정보들도 함수를 기반으로 생성합니다.
- 생성된 주파수 영역의 결과들을 역 푸리에 변환합니다.
- 이를 기반으로 Displacement Mapping, Normal Mapping을 수행합니다.
먼저 우리는 파도를 생성하기 이전에, 푸리에 변환을 이용할 수 있는 기반을 갖춰야 합니다. 고속 푸리에 변환에 대해 먼저 알아보도록 하죠.
병렬 고속 푸리에 변환
이산 푸리에 변환은 N^2의 시간복잡도를 가집니다. 우리가 생성할 Displacement Map의 한 변의 크기는 128~512 픽셀이 될 것이며, 이들을 생성하더라도 초당 60프레임 이상으로 렌더링할 수 있어야 합니다. 렌더링을 위한 수많은 연산들을 고려하였을 때, Displacement Map의 생성에 드는 비용은 결코 싸지 않습니다. 따라서 우리는 고속 푸리에 변환(Fast Fourier Transform)을 이용합니다.
고속 푸리에 변환은 푸리에 변환 시 나타나는 재귀적 특성에 기인합니다. 주기가 반절이면서 같은 파동(정확히는 짝수 번째의 샘플과 홀수 번째의 샘플을 갖는 두 파동)을 이용해서, 두 파동이 합쳐진 파동의 푸리에 변환을 알 수 있습니다. 역시 이 글에서는 고속 푸리에 변환과 병렬 컴퓨팅의 세부사항을 다루지 않습니다.
고속 푸리에 변환은 NlogN의 시간 복잡도를 가집니다. 2차원의 경우 N^2logN이며, 이는 단일 CPU가 매 프레임마다 처리하기 어렵습니다. 따라서 우리는 계산 쉐이더(Compute shader)를 통해 이들을 병렬적으로 처리할 것입니다. 고속 푸리에 변환은 다음과 같은 절차를 통해 이루어집니다:
- 주파수 영역의 값을 RWTexture에 매핑합니다.
- Shift를 수행합니다.
- Bit Reversal을 수행합니다.
- 1-D FFT를 수행합니다.
- Transpose를 수행합니다.
- Bit Reversal을 수행합니다.
- 1-D FFT를 수행합니다.
- Transpose를 수행합니다.
그렇다면 각각의 절차 Shift, Bit Reversal, FFT, Transpose가 무엇인지에 대해서 알아보도록 합시다.
Shift
우리가 사용하는 고속 푸리에 변환 알고리즘은 (0, 0)을 원점으로 하고 [0, N)의 범위를 갖습니다. 하지만 Tessendorf가 제안하는 함수의 경우 [-N/2, N/2)의 범위를 갖습니다. 따라서 이들을 이동시켜주어야 합니다. 이들은 주기적 파동 함수이므로, 같은 타일을 늘어 놓은 채 원점만 변경한다고 생각할 수 있습니다.
참고: (N, 1, 1) 규모의 스레드가 (1, N, K)의 블록으로 디스패치됩니다. 우리는 여러 개의 2차원 주파수 정보들을 다룰 것이므로, z축을 이용해서 K개를 병렬적으로 처리하게 됩니다.
HLSL 코드: Shift()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | void Shift(uint3 xyz) { const uint sizeHalf = gSize * 0.5; if (xyz.y < sizeHalf) { if (xyz.x < sizeHalf) { gOutput[xyz] = gInput[uint3(xyz.x + sizeHalf, xyz.y + sizeHalf, xyz.z)]; } else { gOutput[xyz] = gInput[uint3(xyz.x - sizeHalf, xyz.y + sizeHalf, xyz.z)]; } } else { if (xyz.x < sizeHalf) { gOutput[xyz] = gInput[uint3(xyz.x + sizeHalf, xyz.y - sizeHalf, xyz.z)]; } else { gOutput[xyz] = gInput[uint3(xyz.x - sizeHalf, xyz.y - sizeHalf, xyz.z)]; } } } [numthreads(N, 1, 1)] void ShiftCS(int3 dispatchThreadID : SV_DispatchThreadID) { Shift(dispatchThreadID.xyz); } | cs |
Bit Reversal
고속 푸리에 변환은 짝수 번째와 홀수 번째의 샘플을 이용합니다. 비트 표현의 역순을 이용해서 각 샘플에 접근하면, 재귀함수 없이도 각 샘플들을 반복적으로 손쉽게 접근할 수 있습니다.
참고: (N, 1, 1) 규모의 스레드가 (1, N, K)의 블록으로 디스패치됩니다.
HLSL 코드: BitReversal()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | void BitReversal(uint3 xyz) { int rev = 0; for (int j = 1, target = xyz.x; j < gSize; j <<= 1, target >>= 1) { rev = (rev << 1) + (target & 1); } if (xyz.x < rev) { // swap float4 temp = gOutput[xyz]; gOutput[xyz] = gOutput[uint3(rev, xyz.y, xyz.z)]; gOutput[uint3(rev, xyz.y, xyz.z)] = temp; } } [numthreads(N, 1, 1)] void BitReversalCS(int3 dispatchThreadID : SV_DispatchThreadID) { BitReversal(dispatchThreadID.xyz); } | cs |
설명
rev는 target의 비트가 반전된 결과입니다. target의 비트는 우측으로 shift, rev는 좌측으로 shift 하면서 마지막 자리의 비트가 1이면 rev 역시 1을 더해주는 형태로 결과를 계산합니다. 우리는 주어진 행에 대해서 비트의 역순으로 이들을 배치해야 합니다. 따라서 xyz 중 x만을 다룹니다.
추가로, rev가 target보다 작은 경우에 한해서만 swap을 수행하게 됩니다. 예를 들어, 001의 비트 역연산은 100입니다. target이 001일 때, 100일 때 둘 다 swap을 하게 되면 결론적으로 두 번 swap을 하는 것이기 때문에 위치가 변하지 않습니다. 따라서 한쪽의 경우에서만 swap을 수행합니다.
1D-FFT
우리는 2차원을 다루지만, 2차원에 대해서 푸리에 변환을 수행하지 않습니다. 이는 시간 복잡도가 지수적으로 늘어나는 일이며, 선형 시간 안에 해결할 수 있는 방법이 존재하기 때문입니다. 단순히, 각 차원으로 1차원 푸리에 변환을 수행하면 됩니다. 이를 이용해 N^2logN의 시간 복잡도를 얻을 수 있습니다. 따라서 우리는 1차원 고속 푸리에 변환을 구현합니다. 이해를 돕기 위해, 병렬처리를 하지 않는 1차원 고속 푸리에 변환의 코드를 함께 첨부합니다.
참고: (N, 1, 1) 규모의 스레드가 (1, N, K)의 블록으로 디스패치됩니다.
C++ 코드: FastFourierTransform1d()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | // a의 크기는 2^N 이어야 하고, IFFT의 경우 inverse = true void fft(Complex* a, const int size, bool inverse) { constexpr float PI = 3.141592f; int n = size; // Bit Reversal ... // FFT for (int len = 2; len <= n; len <<= 1) { for (int i = 0; i < n; i += len) { for (int j = 0; j < len / 2; j++) { const float theta = 2.0f * PI * j / len * (gIsInverse ? -1 : 1); const Complex wk = { cosf(theta), sinf(theta) }; const Complex even = a[cur]; const Complex odd = a[cur + len / 2]; a[cur] = even + odd * wk; a[cur + len / 2] = even - odd * wk; } } } if (inverse) { for (int i = 0; i < n; i++) { a[i] = a[i] / n; } } } | cs |
HLSL 코드: FastFourierTransform1d()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | void FastFourierTransform1d(uint3 xyz) { static const float PI = 3.141592f; for (int len = 2, lenHalf = 1; len <= gSize; len <<= 1, lenHalf <<= 1) { int j = xyz.x % len; if (j < lenHalf) { float theta = 2.0f * PI * j / len * (gIsInverse ? -1 : 1); float2 wk = { cos(theta), sin(theta) }; float2 even = gOutput[xyz]; float2 odd = gOutput[uint3(xyz.x + lenHalf, xyz.y, xyz.z)]; gOutput[xyz] = float4(even + ComplexMul(wk, odd), 0.0f, 0.0f); gOutput[uint3(xyz.x + lenHalf, xyz.y, xyz.z)] = float4(even - ComplexMul(wk, odd), 0.0f, 0.0f); } DeviceMemoryBarrierWithGroupSync(); } if (gIsInverse) { gOutput[xyz] /= gSize; } } [numthreads(N, 1, 1)] void Fft1dCS(int3 dispatchThreadID : SV_DispatchThreadID) { FastFourierTransform1d(dispatchThreadID.xyz); } | cs |
설명
1차원으로 고속 푸리에 변환이 수행되므로, 각 행들에 대해서는 독립적입니다. 따라서 이들은 병렬적으로 처리가 가능합니다. 이러한 병렬 처리는 블록 단위에서 수행되게 됩니다.
또한, 짝수 번째의 샘플과 홀수 번째의 샘플을 더하는 것 역시 반절의 스레드가 병렬적으로 처리할 수 있습니다. 가령 256개의 샘플이 있다면, 128개의 스레드가 각각 짝수, 홀수를 하나씩 맡아 정해진 위치에 독립적으로 더할 수 있습니다. 따라서 logN의 시간 동안 반복적으로 이들을 더하게 됩니다. 이때 모든 스레드가 전부 계산을 끝마쳐야 다음 반복으로 넘어갈 수 있음에 주의하세요. 따라서 DeviceMemoryBarrierWithGroupSync()를 이용해 스레드들을 동기화해야 합니다.
Transpose
우리는 1차원 고속 푸리에 변환만을 구현하였습니다. 이는 행에 대해서 연산되며, 열에 대해서는 처리하지 않습니다. 하지만 2차원 고속 푸리에 변환을 수행하기 위해, 열에 대해서도 연산을 수행해야 합니다. 편법으로, 행렬을 전치시키고 다시 행에 대해서 같은 연산을 수행하는 방식을 이용합니다.
참고: (N, 1, 1) 규모의 스레드가 (1, N, K)의 블록으로 디스패치됩니다.
HLSL 코드: Transpose()
1 2 3 4 5 6 7 8 9 10 | void Transpose(uint3 xyz) { gInput[xyz] = gOutput[xyz.yxz]; } [numthreads(N, 1, 1)] void TransposeCS(int3 dispatchThreadID : SV_DispatchThreadID) { Transpose(dispatchThreadID.xyz); } | cs |
2-d Fast Fourier Transform
위 HLSL 코드를 이용해서 2차원 푸리에 변환을 수행하는 CPU 코드는 다음과 같습니다. 각각의 PSO들을 매개변수로 받아서, Command List를 통해 GPU에게 명령합니다.
CPU 코드:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | void OceanMap::ComputeOceanDisplacement( ID3D12GraphicsCommandList* cmdList, ID3D12RootSignature* rootSig, ID3D12PipelineState* shiftCsPso, ID3D12PipelineState* bitReversalCsPso, ID3D12PipelineState* fft1dCsPso, ID3D12PipelineState* transposeCsPso ) { FftConstants c = { mWidth, 1 }; cmdList->SetComputeRootSignature(rootSig); { // UAV 매핑에 필요한 복사 및 상태 전이를 수행합니다. } // HLSL 코드의 gInput과 gOutput에 매핑합니다. cmdList->SetComputeRoot32BitConstants(0, 2, &c, 0); cmdList->SetComputeRootDescriptorTable(1, mhGpuUavDisplacementMap0); cmdList->SetComputeRootDescriptorTable(2, mhGpuUavDisplacementMap1); Shift(cmdList, shiftCsPso, c); BitReversal(cmdList, bitReversalCsPso, c); Fft1d(cmdList, fft1dCsPso, c); Transpose(cmdList, transposeCsPso, c); { // gInput에 gOutput을 복사하는 코드 // 이는 ping-pong 방식의 연산과 // 코드 재사용에 이유가 있습니다. } BitReversal(cmdList, bitReversalCsPso, c); Fft1d(cmdList, fft1dCsPso, c); Transpose(cmdList, transposeCsPso, c); { // 이후 Unordered Access 중이던 리소스들을 // Generic Read 상태로 전이합니다. } } void OceanMap::Dispatch(ID3D12GraphicsCommandList* cmdList) { const auto numGroupsX = static_cast<UINT>(ceilf(static_cast<float>(mWidth) / 512.0f)); const auto numGroupsY = mHeight; const auto numGroupsZ = NUM_OCEAN_FREQUENCY; cmdList->Dispatch(numGroupsX, numGroupsY, numGroupsZ); } void OceanMap::Shift(ID3D12GraphicsCommandList* cmdList, ID3D12PipelineState* shiftCsPso, FftConstants& c) { cmdList->SetPipelineState(shiftCsPso); Dispatch(cmdList); } // 이하 생략 | cs |
주파수 생성: Phillips, hTilde, hTilde0, hTilde0*
푸리에 변환을 이용한 방식에서의 높이 함수는 다음과 같이 기술됩니다.
hTilde 함수는 통계적, 실증적인 연구에 의해 다음과 같이 활용할 수 있음을 입증했습니다:
hTilde0 함수에서 그리스 문자 Xi_r, Xi_i는 0을 평균, 1을 표준편차로 하는 정규 분포에서의 난수입니다. Ph(k)는 Phillips 함수입니다. Phillips 함수에서 L은 바람의 속도 V의 제곱을 중력 상수 g로 나눈 것입니다. w hat은 바람의 방향이며, A는 수치적 상수입니다. 조금 더 자세한 내용은 논문을 참고하세요.
우리는 3차원 공간에서의 Displacement Map을 생성해야 합니다. (Height Map이 아닙니다) 따라서 위와 같은 연산을 3개의 텍스처에 수행해주어야 함에 유의하세요.
Phillips Spectrum
|
| Phillips Spectrum |
Phillips Spectrum은 파도를 만들기 위해 사용되는 함수입니다. 이 함수에는 바람의 세기와 방향, 중력 등이 영향을 미칩니다. 추가적으로 Tessendorf는 계산의 용이성을 위해 Amplitude라는 상수를 도입하였습니다.
HLSL 코드: Phillips()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | // 0 <= n, m < res float Phillips(int n, int m, float amp, float2 wind, int res, float len) { float2 k = { 2.0f * PI * (n - 0.5f * res) / len, 2.0f * PI * (m - 0.5f * res) / len }; // k = 2pi * n0 / L, (-res/2 <= n0 <= res/2) float kLen = length(k); if (kLen < EPSILON) return 0; float kLen2 = kLen * kLen; float kLen4 = kLen2 * kLen2; float kDotW = dot(normalize(k), normalize(wind)); float kDotW2 = kDotW * kDotW; float wlen = length(wind); float l = wlen * wlen / G; float l2 = l * l; float damping = 0.001; float L2 = l2 * damping * damping; return amp * exp(-1 / (kLen2 * l2)) / kLen4 * kDotW2 * exp(-kLen2 * L2); } | cs |
위 코드를 통해 만들어진 주파수 영역의 값은 다음과 같습니다. 보시다시피 진동수의 변화가 아주 부드러우며, 따라서 부드러운 파도가 생성되게 됩니다.
 |
| Phillips의 결과 |
hTilde0
|
| hTilde0 |
상기한 Phillips 함수를 이용해서 hTilde0 함수의 결과를 얻을 수 있습니다. hTilde0는 일종의 Noise를 생성하는 역할을 합니다. 따라서 파도가 거칠어지고, 실제의 바다처럼 변하게 됩니다.
HLSL 코드: HTilde0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | float mod(float x, float y) { return x - y * floor(x / y); } float Rand(float2 uv, float randSeed) { float result = sin(mod(12345678.f, dot(uv, float2(12.9898, 78.233) * 2.0))) * (43758.5453 + randSeed); result = frac(result); return result; } float2 RandNegative1ToPositive1(float2 uv, float randSeed) { float2 r; float2 isNegative = Rand(float2(uv.x, uv.y), randSeed); r.x = Rand(float2(uv.x + 31.2452, uv.y + 27.6354), randSeed); r.y = Rand(float2(uv.x + 11.67834, uv.y + 51.3214), randSeed); if (isNegative.x > 0.5f) r.x = -r.x; if (isNegative.y > 0.5f) r.y = -r.y; return float2(r.x, r.y); } float2 HTilde0(int n, int m, float amp, float2 wind, int res, float len, float randSeed) { return ComplexMul( RandNegative1ToPositive1(float2(n, m), randSeed), sqrt(Phillips(n, m, amp, wind, res, len) / 2.0f) ); } | cs |
설명
HLSL에서 의사 난수 함수를 제공하지 않으므로, 직접 만들어서 사용했습니다. 이는 정규 분포를 따르지 않으므로, 한계를 갖습니다. 다음과 같이 Phillips 함수의 결과에서 Noise가 생긴 결과를 얻을 수 있습니다.
hTilde0* 역시 언급되는데, 이는 단순히 k를 뒤집은 것에 불과합니다. 즉 hTilde0*(k, t) = hTilde0(-k, t) = hTilde(n-k, t)로 치환할 수 있습니다.
 |
| HTilde0의 결과 |
hTilde
hTilde의 경우 앞서 나온 hTilde0의 값을 이용해서 구할 수 있습니다. ω(k)의 경우 Dispersion 함수이며, 다음과 같습니다. g는 중력이며, k는 격자의 점입니다.
이산된 격자의 특성 때문에, 주파수를 제대로 배치할 수 없는 문제가 발생한다고 합니다. 따라서 Tessendorf는 다음과 같이 근사값을 이용하는 것을 제안합니다:
[[]] 기호는 정수 부분만을 취하는 기호입니다.
HLSL 코드: HTilde
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | [numthreads(N, 1, 1)] void HTildeCS( int3 dispatchThreadID : SV_DispatchThreadID) { uint3 basisIndex = dispatchThreadID; basisIndex.z %= 3; float2 hTilde0 = gHTilde0[basisIndex].xy; float2 hTilde0Conj = gHTilde0Conj[basisIndex].xy; float omegat = Dispersion(dispatchThreadID.x, dispatchThreadID.y, gResolutionSize, gWaveLength) * gTime * 0.1f; const float cos_ = cos(omegat); const float sin_ = sin(omegat); float2 c0 = { cos_, sin_ }; float2 c1 = { cos_, -sin_ }; float2 res = ComplexMul(hTilde0, c0) + ComplexMul(hTilde0Conj, c1); // 편의를 위해, e^ikx를 푸리에 변환 이전에 미리 곱합니다. float kx = PI * (2.0f * dispatchThreadID.x - gResolutionSize); float kz = PI * (2.0f * dispatchThreadID.y - gResolutionSize); float2 k = { kx, kz }; float2 x = { dispatchThreadID.x / gResolutionSize - 0.5f, dispatchThreadID.y / gResolutionSize - 0.5f }; const float kDotX = dot(k, x); const float cosKDotX = cos(kDotX); const float sinKDotX = sin(kDotX); const float2 kDotXComplex = { cosKDotX, sinKDotX }; res = ComplexMul(res, kDotXComplex); // ... // 이후에 다시 설명할 부분입니다. // ... gHTilde[dispatchThreadID.xyz] = float4(res, 0.0f, 0.0f); } | cs |
설명
본래 hTilde의 값은 e^ikx를 곱하지 않은 값입니다. 하지만 푸리에 변환을 수행하면서 e^ikx를 곱할 것이기 때문에, 편의상 e^ikx를 곱합니다.
Choppy Wave
상술한 수식으로는 고운 파도를 생성할 수 있습니다. 하지만 날씨가 좋지 않을 때에는 조금 더 급격한 변화를 일으키는 파도를 만드는 것이 좋을 것입니다. 또, 높은 파도는 결국 부서지면서 거품 등을 만들게 되죠. 이러한 현상을 표현하기 위해, Tessendorf는 하나의 수학적 변환을 제안합니다.
 |
| Lie Transform |
이 변환은 완만한 변위 곡선을 가파르게 만듭니다. 다음 그래프를 참고해보세요. 점선의 경우 변환을 수행하지 않은 곡선, 실선의 경우 변환을 수행한 곡선입니다. 조금 더 끝이 뾰족해졌음을 확인하세요.
하지만 이러한 변환은 3차원에서 수행되기 어렵습니다. 따라서 우리가 생성에 이용하고 있는 격자의 차원인 x와 z차원으로만 이러한 연산을 수행하도록 합시다. 추가로, 이 변환의 도함수를 이용해 파도가 부서지면서 거품을 내는지의 여부를 확인할 수 있습니다. 해당 사항은 이 글에서 다루지 않습니다.
HLSL 코드: HTilde (revisited)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 | [numthreads(N, 1, 1)] void HTildeCS( int3 dispatchThreadID : SV_DispatchThreadID) { // 전략 float delta = 1.f / gResolutionSize; float2 dx = { delta, 0.0f }; float2 dz = { 0.0f, delta }; float len = length(k); if (basisIndex.z == 0) { if (len < 0.00001f) { // 너무 작은 수로 나누면 문제가 발생할 수 있습니다. } else { const float2 ikk = { 0.0f, -kx / len }; res = ComplexMul(res, ikk); } } if (basisIndex.z == 2) { if (len < 0.00001f) { } else { const float2 ikk = { 0.0f, -kz / len }; res = ComplexMul(res, ikk); } } // ... // 여전히 추가적으로 설명할 코드가 있습니다. // ... gHTilde[dispatchThreadID.xyz] = float4(res, 0.0f, 0.0f); } | cs |
설명
텍스처는 3차원으로 구성되어 있으며 z축을 기준으로 0번 인덱스에는 x축의 Displacement Map, 1번과 2번에는 각각 y축, z축의 Displacement Map이 있습니다. 단순히, x차원과 z차원을 계산할 때에만 추가적인 연산을 수행하도록 합니다.
Slope & Normal
광학적 효과를 극대화하기 위해서는 파도 표면에 수직한 법선 벡터를 파악할 수 있어야 합니다. 파도의 법선 벡터는 미분을 통해서 알아낼 수 있습니다. 미분하여 접선 벡터를 얻은 뒤, 각 차원(x, z)의 두 벡터를 외적하는 것이죠. 2차원 벡터 x에 대해 hTilde 함수를 미분한 결과는 다음과 같습니다:
혹은 이산 샘플링된 결과이므로, 주변의 값들을 이용해서 순간 변화율을 확인하는 방법도 있습니다. 인공지능에 익숙하다면, 일종의 경사 하강법을 활용하는 셈이죠. 이는 메모리를 절약할 수 있지만, 고주파의 파동에 대해서 정확하지 않은 값을 얻게 됩니다.
HLSL 코드: Slope
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | [numthreads(N, 1, 1)] void HTildeCS( int3 dispatchThreadID : SV_DispatchThreadID) { // 전략 // 경사도를 계산합니다. if (dispatchThreadID.z > 5) { const float kDotDz = dot(k, dz); float2 ik = { cos(kDotDz), sin(kDotDz) }; res = ComplexMul(res, float2(0.0f, kz)); } else if (dispatchThreadID.z > 2) { const float kDotDx = dot(k, dx); float2 ik = { cos(kDotDx), sin(kDotDx) }; res = ComplexMul(res, float2(0.0f, kx)); } gHTilde[dispatchThreadID.xyz] = float4(res, 0.0f, 0.0f); } | cs |
설명
경사도를 담은 주파수 영역의 정보들 역시 메모리에 상주해야 합니다. 따라서 z축 인덱스 3~5번에는 x,y,z Displacement Map을 x로 미분한 결과를, 6~8번에는 z로 미분한 결과를 저장합니다. 나중에 쉐이더 단계에서 이들을 이용해 법선 벡터를 계산하게 됩니다.
Shading
바다의 표면을 그리기 위해서는 해당 픽셀의 색상을 결정해야 합니다. 색이라는 것은 해당 색을 지닌 빛이 우리 시신경을 건들여 인지하게 되는 것입니다. 따라서 우리가 어떤 물체의 표면을 바라봤을 때 그 표면의 색을 결정하는 것은, 외부의 광원들이 해당 표면에 관여한 양이 될 것입니다. 이를 해결하기 위해서는 두 가지만 알면 족합니다: 바로 반사와 굴절입니다.
먼저 반사의 경우를 따져보죠. 반사는 표면을 투과하거나 표면에 흡수되지 못하고 튕겨져 나오는 것입니다. 이때 입사각과 반사각은 같다는 특징이 있습니다. 즉, 입사 벡터와 부딪힌 표면의 법선 벡터를 알고 있다면 반사 벡터 역시 얻을 수 있습니다. 이를 수식화해보면 다음과 같습니다. r은 reflection, i는 incident이며 S는 surface입니다. 표면의 경우 파도에 의해서 출렁거리므로, x와 t의 함수로 나타낼 수 있습니다.
이는 HLSL 내장 함수인 reflection을 통해서 수월하게 계산할 수 있습니다.
다음으로 굴절의 경우입니다. 굴절은 외부와 내부의 밀도 차이 등으로 인해 빛의 속도가 달라지면서 일어나는 현상입니다. 굴절각의 경우 스넬의 법칙을 이용해서 계산할 수 있습니다.
θi를 구하는 방법은 단순합니다. 입사 벡터와 법선 벡터를 정규화하고 내적한 뒤, arccos 연산을 수행하면 됩니다. n1/n2는 공기와 물의 밀도 차이를 나타냅니다. 이 경우, 1.34의 상수를 지닙니다.
우리의 눈에 들어오는 빛의 양은 일정합니다. 물의 밖에서 바다를 바라봤을 때 들어온 전체 빛을 1이라고 생각해봅시다. 일부는 물의 아래에서 굴절이 되면서, 나머지는 물의 위에서 반사가 되면서 들어왔을 것입니다. 즉, R을 반사의 정도, T를 굴절의 정도라고 생각한다면 다음과 같습니다.
반사의 정도를 구하는 방법은 다음 수식을 이용합니다. 우리는 사전에 입사각, 반사각, 굴절각을 전부 알 수 있으므로 손쉽게 구할 수 있습니다. 굴절의 정도는 이를 1에서 빼면 됩니다.
이를 전부 활용한 쉐이더 코드는 다음과 같습니다:
HLSL 코드: Pixel Shader
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | float4 PS(DomainOut pin) : SV_Target { // Material 데이터를 얻습니다. MaterialData matData = gMaterialData[gMaterialIndex]; // 의미적으로는 빛의 난반사에 관한 성질이지만, // 바다 아래에서 오는 굴절광이라고 간주합니다. float4 diffuseAlbedo = matData.DiffuseAlbedo; uint diffuseMapIndex = matData.DiffuseMapIndex; pin.NormalW = normalize(pin.NormalW); float3 toEyeW = gEyePosW - pin.PosW; float distToEye = length(toEyeW); toEyeW /= distToEye; float nSnell = 1.34f; float reflectivity; float3 nI = normalize(toEyeW); float3 nN = normalize(pin.NormalW); float costhetai = dot(nI, nN) ; float thetai = acos(costhetai); float sinthetat = sin(thetai) / nSnell; float thetat = asin(sinthetat); if(thetai == 0.0f) { reflectivity = (nSnell - 1) / (nSnell + 1); reflectivity *= reflectivity; } else { float fs = sin(thetat - thetai) / sin(thetat + thetai); float ts = tan(thetat - thetai) / tan(thetat + thetai); reflectivity = 0.5f * (fs * fs + ts * ts); } // Normal의 오차로 인해, reflectivity가 1.0f을 넘는 경우가 존재합니다. // 이를 제거합니다. reflectivity = clamp(reflectivity, 0.0f, 1.0f); // 바다에 적용된 텍스처를 반영합니다. diffuseAlbedo *= gTextureMaps[diffuseMapIndex].Sample(gsamAnisotropicWrap, pin.TexC); float3 r = reflect(-toEyeW, pin.NormalW); // 환경 매핑을 이용해 반사광의 색을 얻습니다. float4 reflectionColor = gCubeMap.Sample(gsamLinearWrap, r); float3 Ci = (reflectivity * reflectionColor + (1 - reflectivity) * diffuseAlbedo); return float4(Ci, 1.0f); } | cs |
테셀레이션
우리가 렌더링할 격자는 일련의 정점의 집합입니다. 결국 정점들을 잇는 삼각형들로 렌더링될 것이므로, 정점의 수가 너무 적다면 파도가 부드럽지 않고 각져 보일 것입니다. 반면 정점의 수가 너무 많다면 이들을 관리하는데 부하가 걸려 CPU 병목을 유발할 것입니다. 예를 들어, 멀리 있는 파도를 세세하게 그릴 필요는 없습니다.
따라서 동적으로 정점의 수를 변화시킬 필요가 있습니다. 이때 정점의 수는 카메라와 파도를 구성하는 삼각형 사이의 거리에 반비례합니다. 거리가 가까우면 정점을 늘리며, 거리가 멀다면 정점의 수를 줄입니다. 테셀레이션에 대한 세부 사항은 다루지 않습니다.
HLSL 코드: Hull Shader
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 | PatchTess ConstantHS(InputPatch<VertexOut, 3> patch, uint patchID : SV_PrimitiveID) { PatchTess pt; float4x4 world = gWorld; float3 centerL = 0.333f * (patch[0].PosL + patch[1].PosL + patch[2].PosL); float3 centerW = mul(float4(centerL, 1.0f), world).xyz; float d = distance(centerW, gEyePosW); const float d0 = 0.0f; const float d1 = 70.0f; float tess = 2.0f + 6.0f * saturate((d1 - d) / (d1 - d0)); pt.EdgeTess[0] = tess; pt.EdgeTess[1] = tess; pt.EdgeTess[2] = tess; pt.InsideTess = tess; return pt; } [domain("tri")] [partitioning("integer")] [outputtopology("triangle_cw")] [outputcontrolpoints(3)] [patchconstantfunc("ConstantHS")] [maxtessfactor(64.0f)] HullOut HS(InputPatch<VertexOut, 3> p, uint i : SV_OutputControlPointID, uint patchId : SV_PrimitiveID) { HullOut hout; hout.PosL = p[i].PosL; hout.NormalL = p[i].NormalL; hout.TangentU = p[i].TangentU; hout.TexC = p[i].TexC; return hout; } | cs |
결과
위 세 줄은 각 축에 대해 hTilde0, hTilde0*, hTilde, Displacement Map입니다. 맨 아래 세 개는 각 축에 대한 x축 미분값입니다.(slope)
Future
제가 구현한 방법은 여러 가지 한계를 갖습니다. 그 한계는 다음과 같습니다:
- Random이 (-1,1)의 uniform random입니다. Jerry Tessendorf의 글에서는 0을 평균, 1을 표준편차로 하는 정규 분포 난수를 활용합니다.
- Displacement, Normal map 생성에 대한 최적화 방안을 강구해야 합니다. (RTX2070super, 512*512 해상도 기준, 70프레임)
- 단 하나의 격자만을 렌더링하였으며, 연속으로 배치하였을 시 변위가 반복되는 현상(타일링)을 없애야 합니다.
- 테셀레이션의 경계에서 변위 매핑의 정밀도가 달라지면서 메쉬 사이에 작은 틈이 생깁니다.
- 거품을 생성하지 않습니다: Tessendorf는 도함수를 이용한 Folding map을 제안합니다. 이를 이용해서 Choppy wave가 뒤집어지는지의 여부를 확인할 수 있으며, 이를 통해서 파도가 부서진다고 판단하고 거품을 렌더링할 수 있습니다.
참고:
https://github.com/iamyoukou/fftWater/tree/master/src
https://github.com/AlphaMistral/Mistral-Water
소스 코드:
https://github.com/SubinHan/graphics-tutorial/tree/main/WindowsProject1/24Ocean
프로젝트 기간: 2주
끝.
댓글
댓글 쓰기