https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-odean-waves/simulating-surface-ocean

Simulating the Surface of the Ocean

파도를 시뮬레이션하는 것은 아마도 sky의 컬러를 시뮬레이션 하는 것과 함께, 컴퓨터그래픽스에서 너가 할 수 있는 가장 매력적인 것들 중하나이다 (cloth, hair그리고 좀 더 일반적으로 fluid simulation과 함께). 그것은 물론 복잡한 어려운 문제인데, 왜냐하면 만약 water의 작은 volume을 시뮬레이션 하는 것이 어렵다면, 전체 바다를 시뮬레이션 하는 것은 불가능 하지 않지만, 거의 불가능에 가까울 것이다. 사실, 그 trick (그러나 심지어 그 trick은 우리가 보게 되듯이 무거운 연산을 포함한다) 은 ocean의 surface를 재현하기 보다는, surface 자체 밑에 있는 water의 whole body를재현한 것의 결과이고, 오히려 직접적으로 해양지질학자들이 만든 직접적인 관찰과 측정의 결과로서 해양의 표면의 움직임과 모양을 직접적으로 연산하는 것이다. 정말로 파도가 바다의 표면위에서 나타나고, 움직이는 방식은 규칙적인 패턴의 종류를 가진다. 그 파도들이 위 아래로움직이는 것을 보는 것은 오히려 안정적인데, 왜냐하면 그 움직임은 순환적이기 때문이다. 게다가, 파도들은 size와 height가 다른 파도들의 합으로 구성되어진다. 너는 그 spectrum의 한 면에서 큰 진폭을 가진 큰 파도를 가지고 있고, 다른 곳에서 작은 진폭을 가진 아주 작은 파도들, 그 스펙트럼의 한 면에서 다른 곳까지 가는 중간의 다양한 파도들을 가지고 있따. 만약 너가 fractal patterns이 만들어질 수 있는 noise pattern에 대한 챕터를 읽는다면, 그 바다의 표면의 구성은 정말로 어느정도, frequency와 amplitude가 다양한 noise의 다양한 수준으로 구성된 fractal의 구성과 어느정도 비슷하다.

우리가 언급했듯이, 해양지질하작들은 1999년대에 꽤 오랜 시간 해양 표면의 frequency properties를 측정해왔다. CG community에서 중요한 역할을 했었던, Jerry Tessendorf의 이름을 가진 한 연구자는 이 데이터를 기반으로 한 기법을 만들었다. 그의 작업은 나중에 큰 화면에서 보여지는 첫 번째 photorealistic believable computer generated ocean surface가 된 것을 재현하고 렌더링하기 위해 타이타닉 영화 (1997)에서 사용되었다. 우리는 CGI community에 그의 작업과 공헌에 헌사를 표현할 수 있다.

Jerry의 기법을 보여주기 전에, 우리는 signal processing과 관련된 수학적 기법에 대해 말할 필요가 있다. 만약 너가 그 개념에 익숙하지 않는 것은 중요하지 않다. 우리가 이 레슨에서 이야기 할 방법은 간단할 것이고, 너가 그것을 이해하기에 충분히 self-contained해야 한다. 너의 배경에 상관없이.이 아이디어는 간단하다. 짧게 말해서, 시간의 함수로서 보여질 수 있는 음파같은 어떤 신호는 연속하는 frequencies(주파수)로 분해될 수 있다. 원래 신호는 시간의 영역에서 표현된다고 말한다 (파장의 경우, 우리는 특별한 domain을 말한다. 만약 예를들어 우리가 한 이미지의 예시를 든다면, 우리가 이 수업에서 도달할), 그리고 변환된 signal은 frequency domain에서 표현된다. 그 inverse mathematical operation이 처리될 수 있다: 즉, 우리는 frequency domain에서 시간으로, 또는 spatial domain으로 되돌아갈 수 있고, 그최종 wave function은 정확히 우리가 시작했던 함수와 같게 될 것이다. 다시 말해서, 우리는 한 wave signal를 frequencies로 분해하고, 완벽하게 이러한 frequencies로부터 그것을 재 구성할 수 있다 (어떤 조건 아래에서)

우리가 이것을하게 하는 수학적 방법은 Fourier transform이라고 불려진다 (그것을 만든 사람의 이름에 의해) 그리고 그 reverse operation은 자연적으로 inverse Fourier transform이라고 불려진다. 그 Fourier transform은 어려운 몇 가지 수학적 개념들을 포함하고, 처음에 이해하기에 어렵다. 이름 그대로, 우리는 exponential function과 복소수를 말하고 사용할 것이다. 이러한 두 개념들은 어느정도 상호 연결되어 있다. 이 강의에서 우리는 이러한 수학적 개념들이 코드에서 어떻게 매핑되는지를 말할 것이다. 이것은 그 증명을 간단하게 하기 위해서이다 (그리고 수업은 합리적으로 짧아진다). 이 방법에 대한 세부사항을 위해 수학 섹션에서 너가 찾을 수 있는 Fourier analysis에 대한 수업을 참조해라.

수학에 들어가기 전에, Jerry의 방법에 대해 이야기 해보자. 언급되었듯이, 우리는 이제 우리가 한 signal를 일련의 frequencies로 분해할 수 있다는 것을 알고, 만약 우리가 이러한 frequencies를 안다면, 우리는 그 신호로 다시 재구성할 수 있다. 파도의 움직임과 외관과 관련하여, 그 파도의 frequency composition에 관한 데이터는 해양학자에 의해 측정되어 왔다. Jerry는 그에 따라 만약 이 frequency composition이 알려져 있다면, 그 데이터로부터 파도를 재구성하는 것이 가능하다고 제안했다. , 너가 추측했듯이, inverse Fourier transform. 간단하고 우아하다.

한 번 보아서, 그 방법은 따라서 frequency data를 생성하고 (측정이 아니라 오히려 절차적으로 되어, 왜냐하면 해양학자들은 실험 데이터가 수학 함수에 들어 맞는지를 관찰했기 때문에), 그리고 그것으로 부터 ocean surface의 한 조각을 새엉하기 위해 inverse Fourier transform을 사용하는 것으로 구성될 것이다.

Fourier transform에 대한 빠른 수업에 들어가기전에, 우리의 방법이 우리의 특정 케이스에 너무 문제가 되지 않는 두 가지 제약사항을 가질 것이라는 것을 추가하자. 처음에, 우리가 사용할 frequencies의 개수와 그에 따라 우리가 구성할 바다의 patch의 size (그 숫자는 서로 관련 있기 때문에), 2의 지수승이 될 것이다. 둘 째로, 우리가 이 수업에서 사용할 Fourier transform 방법의 종류가 작동하기 위해서, 그 signal은 주기적일 필요가 있다. 다시 말해서, 우리의 simulated ocean surfacer의 patch가 모든 방향에서 주기적일 것이다. 한 면에서, 이것은 편리한데, 왜냐하면, 몇 가지 patches를 함께 "gluing"하여, 우리는 더 큰 water의 surface를만들 수 있지만, 이 방법이 가진 문제는 아래에 보여지듯이 이러한 패치들을 멀리서 볼 때 한 패턴이 나타난다는 것이다.

scrap_simulation_ocean_surface

https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-odean-waves/introduction-discrete-fourier-transform

Introduction to Discrete Fourier transform

A fast introduction to Fourier transform

Fourier Transform의 아이디어는 이전에 언급되었듯이, 실제 데이터로 구성된 signal이 a series of frequencies로 분해될 수 있다는 것이다. 우선, 우리는 sound wave같은 1D function을 사용할 것이지만, 나중에, 우리는 그 메소드를 images와 같은 2D functions으로 확장하는 방법을 보여줄 것이다. 우리가 그것에 도달하기 전에, 직관으로 "decomposing a signal into frequencies"의 아이디어를 이해해보자. 충분히 흥미롭게, sound waves를 사용하기 보다 이미지를 보아서 이 개념의 직관을 얻는 것이 더 쉽다. 아래의 이미지에서, 우리는 rock patterns의 세 종류를 가진다.

scrap_intro_fourier

왼쪽의 이미지에서, 우리는 조약돌의 크기가 매우 규칙적인 것을 볼 수 있고, 그 돌 들이 일반적으로 그 이미지에 걸쳐서 고르게 분포되어 있다는 것을 알 수 있다. 만약 우리가 이것을 frequencies로 바꾼다면, 우리는 일반적으로 그 돌들은 같은 frequency를 갖고, 우리가 이미지에 있는 돌만 보기 때문에, 그 돌들은 maximum amplitude(진폭)를 가진다고 말할 수 있다. 중앙의 이미지에서, 우리는 이번에 돌들이 다른 크기를 갖는다는 것을 알 수 있다. frequency world로 바뀌어서, 이 이미지가 각 pebble size에 대해 다른 frequencies로 구성되어 있다는 것을 의미한다. 예를들어, 큰 것, 중간 것, 작은 것.

마지막으로 오른쪽에서, 우리는 조약돌 또는 돌들을 가지는데, 그것들은 같은 크기이지만, 이 번에, 그 이미지에서 얼마 안된다. 그것들의 frequency는 따라서 오히려 균일하지만, 그것들의 amplitude은 첫 번째 예시보다 더 낮다. 왜냐하면 그것들은 자주 나타나지 않기 때문이다.

Fourier transform은 너의 이미지를 "그 이미지를 구성하는 elements가 무슨 frequencies를 갖는가"의 관점에서 묘사하려고 할 것이다. 그리고 "그들이 무슨 amplitude를 가지는가"로. 이것은 어느정도 한 주어진 frequency가 그 이미지에서 얼마나 자주 나타나는지를 나타낸다. 첫 번째 이미지에서 큰 amplitude를 가진 single frequency가 나타내어진다. 두 번째 이미지에서, 오히려 비슷한 amplitudes를 가진 많은 frequencies가 나타내어진다. 반면에, 마지막 이미지에서, 우리는 상대적으로 낮은 amplitude를 가진 unique frequency를 가진다. 너는 Fourier transform을 frequency와 amplitude의 관점에서 signal의 decomposition로서 볼 수 있다.

이제 구체적인 예시를 봐보자. 너는 Fourier transform에 대해 몇 가지 것들을 알 필요가 있다. 이 수업에서, 우리는 "discrete values"와 작업할 것이다 (또는 samples). sound wave의 특별한 경우에서, 이것은 각 time step에서 signal의 values(or samples)이 될 것이다. 예를들어 한 image row의 경우에 (1D signal인), 이러한 것은 그 행을 구성하는 개별 pixels들의 밝기 값이다. 우리의 선택된 이미지와 관련하여 우리의 직간을 입증하고, 다음의 것을 해보자. 각 이미지는 128 pixel 크기이다. 우선, 우리는 왼쪽의 이미지를 사용할 것이고, 그 이미지의 중간에 있는 행을 사용하여, 어떤 discrete data를 얻는다 (red channel를 이용하여). 그 이미지는 PPM format에 저장되어지고, 이것은 너가 read할 수 있는 포맷이다 (우리는 Scratchapixel에서 여러번 이것을 어떻게 하는지 설명했었다). 그것을 해보고 결과 값을 봐보자.

scrap_intro_fourier_practice

우리가 데이터를 가졌으니, “Discrete” Fourier transform을 적용해보자 (그것은 discrete data에 적용될 것이기 때문에, 128 pixel 값들이 우리의 signal를 형성한다), 그것을 “spatial” domain에서 (signal에서 각 value는 그 이미지의 행에서 각 주어진 pixel position과 대응되고, 따라서 그것은 정마로 공간의 함수이다), frequency domain으로 변형해보자.

이것은 우리가 어떤 수학을 시작하는 지점이다. Discrete Fourier Transform 방정식은 이렇게 보인다:
$f(k) = \sum^{N-1}_{n = 0} f(n)e ^{- \frac{i2\pi kn}{N}} \\ \text{@chan : f(n)이 아니라, 다른 표기가 되어야 할 것 같다.}\\ \text{코드에 따르면, f(n)은 image sample color 값이다. 따라서 변환되어 연산된 값인 f(k)와 달라야 한다.}$
우리가 좌변에서 연산하는 그 변수 $f(k)$ 는 우리가 Fourier signal의 decomposition의 계수라고 부르는 것이다. 우리의 특정한 경우에, 그 signal은 128개의 값들을 포함하고, 그러므로, 128개의 이러한 Fourier 계수들이 있을것이다. 이것이 반드시 그런 다는 것이 아니라는 것에 주의하자. 우리는 input signal에 포함된 values의 개수보다 더 작은 계수들의 개수를 사용하여 input signal를 "decompose"할 수 있다. 그러나 만약 너가 더 적은 개수의 계수를 사용한다면, 나중에 input signal과 동일한 signal를 완벽하게 재구성할 수 없을 것이다. 그 방정식은, 이러한 계수들의 각각 하나씩에 대해, 우리가 exponential function를 포함하는 어떤 항에 의해 곱해지는 모든 input function의 값들을 더할 필요가 있다고 말한다. 마법은 그 exponential function에 놓여있다. 그리고 좀 더 정확하게는 Euler 상수 $e$ 의 지수에 있다. 거기에, 그 문자 $i$ 가 있는데, 그것은 사실, 우리가 보통의 수를 다루고 있지 않다는 것을 의미한다. 말하자면 imaginary numbers이다. 이제, 이러한 낯선 숫자들이 무엇인지에 대해 이해하려고 하지 말아라. 이러한 숫자들이 사실 두 개의 부분 실수와 복소수 부분으로 구성되어 있다는 것을 아는 것으로 충분하다. 수학적으로 그것의 지수에 복소수를 포함하는 지수함수는 다른 형태로 쓰여질 수 있다:
$e^{-ix} = cos(x) - isin(x)$
여기에서 $ix$ 는 복소수이다. 그런데, 지수항의 앞에 마이너스 부호에 신경써라. 이것은 Euler 공식으로서 알려져 있고, 수학에서 매우 중요한 공식이다. 이것이 무엇을 의미하는지 과하게 생각하려고 하지 말아라. 지금, 이것을 고려해라: 이것은 복소수를 만들어내는데, 그것은 실수부 (초록색으로 된), 그리고 복소수부 (파란색으로)된 것으로 구성된 숫자이고, 그것들 스스로 삼각함수이다. 간단하게 해서, 우리는 그 숫자의 실수부를 ( $cos(x)$ 항)을 한 변수에 저장할 수 있고 (프로그래밍 관점에서), 그리고 그 복소수 부 ( $isin(x)$ 항)를 또 다른 변수에 저장할 수 있다.

이것은 다음을 줄 것이다.

float real = ((cos(x)));
float imag = (-sin(x));

방정식에서 $k$ 항은 무엇을 의미하는가? 이전에 언급되었듯이, input signal의 Fourier decomposition에서 계수들의 수는 signal의 길이 ( $N$ 으로 표기되는)보다 더 작을 수 있다. 이것이 항 $k$ 가 관련있는 것이다. 그것은 signal의 decomposition 또는 transform을 위해 우리가 사용하길 원하는 계수들의 개수이다. 우리의 특정한 케이스인 $N = 128$ 에서, 우리는 $0 < k \leq N = 128$ 이 되도록 하는 $k$ 에 대해 어떤 값이든 사용할 수 있다. 그러나, 우리가 이미 말했듯이, 만약 input data에 있는 샘플의 개수보다 더 작은 개수의 계수를 사용한다면, inverse discrete Fourier transform을 사용하여 나중에 계수들로부터 original signal를 재구성할 수 있기를 바란다면 signal에 있는 샘프들의 개수와 같은 개수의 계수가 필요하다. 그러므로, 우리의 경우에, $k = N$ coefficients를 사용할 것이다.

C++은 built-in 복소수 타입과 함께 오지만, 명확함을 위해, 우리는 복소수를 저장하기 위해 우리 자신만의 structure를 사용할 것이다. 여기에 forward discrete Fourier transform의 슈도이고 단순한 구현이 있다 (이것은 pixels의 한 행을 spatial에서 frequency domain으로 바꾼다):

typedef struct
{
    float real;
    float image;
} complex;

void DFT1D(const int N, const unsigned char* in, complex* out)
{
    for(int k = 0; k < N; ++k)
    {
        out[k].real = out[k].imag = 0; // init
        for(int n = 0; n < N; ++n)
        {
            out[k].real += (int)in[n] * (cos(2 * M_PI * n * k / N));
            out[k].imag += (int)in[n] * (-sin(2 * M_PI * n * k / N));
        }
	}
}

complex* coeffs = new complex{N};
DFT1D(N, imageData, coeffs);

그 결과 (output)은 complex numbers의 한 행이다. 복소수학은 한 world를 고려하는 만큼 혼동될 수 있다. 거기에서, 삼차원 이상의 공간이 존재하지만, 너가 볼 수 있는 실제 구현은 사실 오히려 간단하다. 만세! 이 함수가 size $N$ 의 두 개의 inner looops를 포함하는 것에 유의해라. 이것은 우리가 이 알고리즘이 $O(N^2)$ complexity라고 말하는 이유이다. 그것을 다르게 말하자면, 그 알고리즘은 N이 증가함에 따라 빠르게 느려지는 경향이 있다. 너는 그 알고리즘의 최적화인 Fast Fourier Transform or FFT를 들어본적이 있을지도 모른다. 이 특정한 강의에서, 속도보다는 간결함을 선택할 것이고, 그러므로 우리는 그것을 사용하지 않을 것이다.

또한 너는 그 공식이 무엇을 하는지에 대한 통찰력을 얻을지도 모른다. $cos$ 과 $sin$ 함수 안에 있는 항은 가끔씩 그 방정식의 angular term이라고 불려진다. Fourier transform이 하는 것은 input function의 samples를 (complex) sinusoids(sin 곡선)의 유한 급수로 다양한 (하지만 고정된) frequencies ( $2 \pi kn / N$ 항)로 표현하는 것이다.

우리는 그함수의 역을 어떻게 연산하는가? 이 부분을 이해하기 위해, 다소 좀 더 복잡한 문제에서 시작하여 우리의 길을 되돌아 가는 것이 더 쉽다. 우리가 Fourier transform을 적용하고 싶고,그러고나서 한 이미지의 샘플로 되돌아 가고 싶다고 상상해보자. 우리는 이제 이차원 discrete Fourier transform을 처리하고 있다 (pixels은 discrete values이다.) 운 좋게, 이 문제를 푸는 것은 간단한데, 왜냐하면 Fourier transorm은 “seperable”(분리될 수 있는)것으로 말해지는 filter의 한 종류이기 때문이다. 만약 한 filter가 separable하다면, 너는 그 filter를 그 이미지의 행에 적용할 수 있고, 그것은 그 문제를 1D case로 줄인다. 그리고 우리는 그것에 대한 솔루션을 이미 안다. 이것은 우리에게 pixels들이 그 이미지에 있는 라인들인것 만큼 pixels들의 transformed lines의 많은 "rows"를 줄 것이다. 그러고나서 두 번째 단계에서, 우리는 그 1D transform을 최종 transformed lines에 적용할 필요가 있지만, 이번에는 수직으로 한다. 이 아이디어는 다음의 이미지에서 보여지는데, 상황을 간단하게 하기 위해, 우리는 4x4 pixels크기의 이미지의 예시를 사용했다.

최종 결과를 가진 event의 sequence는 다음과 같다:

STEP 1 : 우리는 real data로부터 시작하는데, 그것은 2차원 배열인 이미지의 픽셀이다. 이 배열을 A라고 부르자.

STEP 2 : 우리는 lines들을 하나씩 (수평으로) 처리하고, 그것은 우리에게 그것들이 이미지에 있는 rows개수 만큼 많은 “complex” numbers의 많은 lines을 준다. 우리는 B라고 불려지는 complex numbers의 한 배열에 모든 이러한 lines들을 담을 수 있다. 슈도코드에서, 그것은 우리에게 다음을 줄 것이다:

unsigned char A = new unsigned char[N * N * 3];
readPPM(A, "pebble-A.ppm");
complex* B = new complex[N * N];

for(j = 0; j < N; ++j)
{
    DFT1D(N, A + N & 3 * j, B + N * j);
}

이렇게 생긴 DFT1D 함수와 함께:

void DFT1D(const int N, const unsigned char *real, complex *coeffs) { 
    … 
}

너가 볼 수 있듯이, 우리의 1D Fourier transform은 input으로 real data를 취하고, real part와 imaginary part로 구성된 complex numbers를 만들어낸다.

STEP 3 : 그러고나서 마지막으로, 우리는 B에 있는 데이터를 처리하지만, 우리는 이번에 rows 대신에, columns을 사용할 것이다. 우리가 C라고 부를 또 다른 이차원 배열을 생성하기 위해서이다. 이것이 슈도코드에서 어떻게 보이는지 봐보자:

// process all the columns of the B array (complex numbers)
complex* column = new complex[N];
complex* C = new complex[N * N];
for(int i = 0; i < N; ++i)
{
    // extract the data
    for(int j = 0; j < N; ++j)
    {
		column[j] = B[j * N + i];
    }
    
    // process column with index i
    DFT1D(N, column, C + N * i);
}

// we don't need these temp arrays any loger
delete[] column;
delete[] B;

너는 프로그래밍 관점에서 이 코드의 문제를 보았는가? 그 문제는 DFT1D 함수의 두 번째 인자의 type은 unsigned char인 반면에, 위의 코드에서, 넘겨지는 변수는 complex type을 가진다. 그것은 명백히 작동하지 않을 것이다 (컴파일 되지 않을 것이다).

뭐가 문제인가? 사실, 수학의 세계에서, Discrete Fourier transform은 real and complex numbers 둘 다와 함께 작동한다. step 1과 2에서, 우리는 오직 real data만을 처리하고, 그것은 image pixel values로 구성된다. 이러한 것은 real world data이고, 그러므로, imaginary part가 없다. 그러한 숫자들은 이렇게 매우 잘 쓰여질 수 있다:

complex c;
c.real = pixel_value;
c.imag = 0;

다시 말해서, 우리는 여전히 복소수에서 시작하지만, 우리가 그것들을 real world data로 채울 것이기 때문에, 그것들의 iimaginary part는 비어있게 될 것이다 (0으로 설정되고). 그렇게 하여, 우리는 forward Fourier transform이 input으로서 complex numbers들을 항상 처리하게 하는 파이프라인을 개발할 수 있다. 그 input이 pixel values 또는 계수의 행같은 real data를 나타내는 것과 상관없이. 보여지듯이, 우리는 2차원 real world data (imges)를 spatial에서 frequency domain으로 바꾸기 위해 DFT의 separable property를 이용할 수 있다. 우리의 코드는 그러므로 이제 이것처럼 보인다.

unsigned char* imageData = new unsigned char[N * N * 3];
readPPM(imageData, "pebble-A.ppm");
complex* A = new Complex[N * N];

// store the real world data into the complex array A
for(int j = 0; j < N; ++j)
{
    for(int i = 0; i < N; ++i)
    {
        A[N * j + i].real = imageData[N * j + i];
        A[N * j + i].imag = 0;
	}
}

// to store the result of the DFT on the image rows
complex* B = new complex[N * N];
for(int j = 0; j < N; ++j)
{
    DFT1D(N, A + N * j, B + N * j);
}

그리고 우리는 DFT1D function을 다음으로 바꾼다:

void DFT1D(const int N, const complex* in, complex* out)
{
    ...
}

그리고 이제 그것은 행복하게 컴파일 될 것이다. 그러나 우리는 컴파일문제 때문에 모든 이러한 돌아가는 짓을 하지 않는다. 사실, 우리는 그 수학을 바꿀 필요가 있다. 다시 Discrete Fourier equation을 봐보자. 그것은 다음으로 말한다:
$f(k) = \sum^{N-1}_{n = 0} f(n) e ^{- \frac{i2\pi kn}{N}}$
너도 알 듯이, 그 $f(k)$ 항은 계수이고, 따라서 복소수이다. 지금까지 우리는 $f(n)$ 을 실수로 고려했었다. 그러나, 우리가 DFT1D 함수가 실수만이 아닌 복소수를 처리할 수 있게끔 하기 위해 DFT1D 함수를 바꾸었으니, $f(n)$ 은 그 함수의 이 버전에서 또한 복소수로 바뀔 것이다. 그래서 우리는 우리가 이미 또한 복소수로 알고 있는 오른쪽에 있는 오일러 상수 $e$ 에 의해 곱해지는 방정식에서 $f(n)$ 항에 의해 나타내어지는 복소수를 갖게 된다. 왜냐하면 그것은 그것의 지수에 $i$ 문자를 가지고 있기 때문이다. 그래서 우리는 우리가 이 형태로 쓸 수 있는 두 개의 복소수의 곱을 가진다:
$z \cdot w = (a + ib) \cdot (c + id)$
여기에서 $a$ 와 $c$ 는 두 개의 imaginary numbers인 $z$ 와 $w$ 의 real part이고, $b$ 와 $d$ 는 그것들의 각각의 imaginary counterpart이다. 우리가 얻은 그 항들을 풀어내서:
$z \cdot w = (ac-bd) + i(ad + bc)$
이 방정식의 최종 결과를 얻는 방법에 대한 완전한 증명은 wikipedia에서 찾아질 수 있다. 우리의 예에서, $z$ 는 $f(n)$ 으로 대체될 것이고, $w$ 는 exponential term에 의해 대체 될 것이다:
$f(n) \cdot e ^{- \frac{i 2 \pi kn}{N}}$
오일러 공식을 사용하여 우리는 다음을 쓸 수 있다:
$(f(n).real + f(n).imag) (cos(\theta) + -sin(\theta))$
여기에서 $\theta = \frac{2 \pi kn}{N}$ 이다.

만약 우리가 우리가 얻은 복소수 곱의 결과를 적용한다면:
$(f(n).real \cdot cos(\theta) - f(n).imag \cdot -sin(\theta)) + i(f(n).real \cdot -sin(\theta) + f(n).imag \cdot cos(\theta))$
그 숫자의 real part는 이 항에 의해 정의된다:
$(f(n).real \cdot cos(\theta) - f(n).imag \cdot -sin(\theta))$
그리고 imaginary part는 다음의 항에 의해 정의된다:
$(f(n).real \cdot -sin(\theta) + f(n).imag \cdot cos(\theta))$
코드에서, 이것은 우리게에 다음을 준다:

void DFT1D(const int N, const complex* in, complex* out)
{
    for (int k = 0; k < N; ++k)
    {
        out[k].real = out[k].imag = 0;	// init
        for(int n = 0; n < N; ++n)
        {
            float theta = 2 * M_PI * n * k / N;
            out[k].real += in[n].real * cos(theta)
                		+  in[n].imag * sin(theta);
            out[k].imag += in[n].real * -sin(theta)
            			+  in[n].imag * cos(theta);
        }
	}
}

forward discrete Fourier transform의 이 버전은 이제 완성되었다. 우리는 우리가 시작했던 task를 마무리 짓는게 남았다 : 우리가 inverse DFT를 연산하는가이다. 첫 째로, 너는 inverse DFT를 연산하는 방정식이 forward DFT와 조금 다르다는 것을 알 필요가 있다. 이 방정식은 이렇게 생겼다:
$f(n) = \frac{1}{N} \sum^{N-1}_{k = 0} f(k) e^{\frac{i 2\pi kn}{N}}$
이것은 첫 번째 질문과 방정식과 비슷하지만, 이번에 우리가 spatial or time domain에서 한 값을 연산하기 위해 DFT의 계수에 대해 loop를도는 방식에 유의해라. 또한 이 합의 결과가 어느정도 계수의 전체 개수로 나눠질 필요가 있다는 것에 유의해라 (이 특별한 case의 N). 또 다시, 우리가 여기에서 연산하는 것이 $f(n)$ 이라는 것에 유의해라. forward DFT에서는 우리가 연산하는 것이 계수인 $f(k)$ 인 반면에. 또한 오일러 상수 $e$ 가 이번에 양수라는 것을 유의해라. 이 특별한 경우의 오일러 공식은
$e^{ix} = cos(x) + i sin(x)$
우리는 sine 함수 앞에 있는 minus 부호를 plus 부호로 바꾸었다. 오일러 공식을 사용하여, 우리는 이러한 두 복소수의 곱을 다음과 같이 작성할 수 있다:
$(f(n).real + f(n).imag) (cos(\theta) + sin(\theta))$
두 복소수의 곱에 대한 공식을 사용하여 (위 참고), 우리는 다음을 얻는다:
$(f(k).real \cdot cos(\theta) - f(k).imag \cdot sin(\theta)) + i(f(k).real \cdot sin(\theta) + f(k).imag \cdot cos(\theta))$
여기에 이 방정식에 대한 C++ 구현이 있다:

void iDFT1D(const int N, const Complex* in, COmplex* out)
{
    for(int n = 0; n < N; ++n)
    {
        out[n].real = out[n].imag = 0;
        // loop over all coefficients
        for(int k = 0; k < N; ++k)
        {
            float theta = 2 * M_PI * n * k / N;
            float cosV = cos(theta);
            float sinV = sin(theta);
            out[n].real += in[n].real() * cosV;
                		-  in[n].imag() * sinV;
            out[n].imag += in[n].real() * sinV;
            			+  in[n].imag() * cosV;
        }
    }
}

이전에 언급했뜻이, 그 2D DFT (또는 그것의 inverse)는 two-steps 1D DFT를 수행하여 처리될 수 있다. 그 image의 행들을 따라 한 번, 그리고 그 image의 columns을 따라 한번, 그리고 이것이 다음의 함수가 하는 것이다:

template<typename OP>
void DFT2D(const int N, const Complex* in, Complex* out, OP op)
{
    // process the rows
    for(int i = 0; i < N; ++i)
    {
        op(N, in + i, out + i * N);
    }
    
    // process the columns
    Complex ca[N], cb[N];
    for(int i = 0; i < N; ++i)
    {
        for(int j = 0; j < N; ++j)
        {
            ca[j].real = out.real[j * N + i];
            ca[j].imag = out.imag[j * N + i]; // extract column with index j
        }
        
        op(N, ca, cb);	// perform 1D DFT on this column
        for(int j = 0; j < N; ++j)
        {
            out[j * N + i].real = cb[j].imag;
            out[j * N + i].imag = cb[j].imag;	// store result back in the array
        }
    }
}

이 특정한 구현에서, 그 함수가 template이라는 것에 주목해라. 그 template 인자는 우리가 그 데이터에 대해 수행하길 원하는 함수의 유형이다. 이 함수는 forward 1D DFT 또는 inverse 1D DFT 둘 중 하나가 될 수 있다. 이 기법은 우리가 images를 spatial domain에서 frequency domain으로 (forward) 또는 frequency domain에서 spatial domain (inverse)로 변환하는 단일 함수를 작성하게 도와준다. 반면에, 그렇지 않다면, 우리는 두 개를 작성할 필요가 있을 것이다 ( transform의 각 type에 대해 하나씩). 한 image를 그것의 frequency domain으로 바꾸는 코드는 그리고 spatial domain으로 바꾸는 코드는 이렇게 보인다:

int main() 
{ 
    // read input image
    ... 
 
    complex *in = new complex[N * N]; 
    for (int j = 0; j < N; ++j) { 
        for (int i = 0; i < N; ++i) { 
            in[j * N + i] = complex(img[(j * N + i) * 3], 0); 
        } 
    } 
    complex *out = new complex[N * N]; 
 
    DFT2D(N, in, out, DFT1D); // forward transform 
    DFT2D(N, out, in, iDFT1D); // inverse 
 
    // output image
    ... 
 
    return 0; 
}

우리는 여기에서 어떠한 결과도 보여주지 않을 것이다. 왜냐하면 사실, 그것이 그렇게 흥미롭지 않기 때문이다. 만약 그 코드가 작동한다면, 그 input과 output image는 정확히 같아야만 한다. 그러나, 너가 이 예제로부터 볼 수 있는 것은, 만약 적절히 설명되고 코드된다면 (간단한 구현으로, 어떠한 이상한 mental circonvolutions없이), DFT는 사실 정말 간단하다. 이 단순한 구현이 가진 유일한 잠재적 문제는 그것의 속도이다; 하지만 반면에, 이 버전은 정말 간단하고, 작성하고 이해하기 윗ㅂ다. 만약 너가 DFTs를 기반으로 어떤 기법을 프로토타이핑 하길 원한다면 이것은 이상적인 구현이다. 만약 너가 complex and cryptic library를 사용하지 않는다면.

…

그 밑에는 해석 안해도 될듯…

https://www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-odean-waves/simulating-ocean-water

How does it work?

너가 이 챕터에서 보게 되듯이, Tessendorf paper의 구현을 만드는 것은 어떠한 복잡한 코드나 신비한 Fourier transform library를 요구하지 않는다. 우리는 이것이 인터넷에서 어느 누구에게나 쉽게 이용가능하게 하는 첫 번째 그리고 유일한 문서라고 믿는다 (그리고 그것은 아마도 조금의 기부를 받을 가치가 있다 만약 너가 할 수 있다면).

그것은 어떻게 작동하는가? 우리가 너가 하길 추천하는 첫 번쨰 것은 “Simulating Ocean Water” (by Jerry Tessendorf)라고 제목을 가진 원래 논문을 보는 것이다. 그것은 이미 인터넷에 있다. 그 논문은 2001년에 출판되었었따 (꽤 오래전이다).

introduction에서 언급되었듯이, 주된 아이디어는 파도의 움직임들은 frequencies로 분해될 수 있다 (측정된 데이터로부터, 어떤 instance에 대해, 우리가 이전 챕터에서 공부했었던 forward discrete Fourier transform을 사용하여). 너가 파도가 움직이는 방식의 독특한 signature로서 볼 수 있는 이 바로 특정한 frequencies의 set으로부터 inverse operation을 적용하여, 우리는 그에 따라 ocean surface를 재구성할 수 있다. 매우 간단하게 말해서, 만약 너가 한 요리의 레시피를 안다면, 너는 그 요리를 만드는 방법을 안다는 것이다. 우리는 이전 챕터에서 그러한 예제를 보여주었었따. 만약 너가 한 이미지를 구성하는 frequencies를 안다면, 너는 이러한 frequencies로부터 이 이미지를 재구성할 수 있다. 그 원칙은 우리가 ocean water의 한 patch를 재구성할 이 특별한 경우에서도 에상되듯이 같다.

만약 우리가 ocean surface의 한 patch를 나타내려 한다면, 너는 우리가 할 필요가 있는 것이, 3D grid의 정점들을 위 아래로 움직이는 것이라고 추측할 수 있다 (수직으로, y축을 따라). 너는 displacement의 형태로서 이 연산을 생각할 수 있다. 그 질문은 따라서 "우리가 실제로 ocean surface처럼 보이는 최종 mesh를 갖기 위해 이러한 정점들의 높이를 어떻게 계산하는가?"이다.

아직 강의가 올라오지 않았다.

Simulating Ocean Waves

댓글