EXQT EXQT
cover

Ray Marching - 1

2023-10-10

Ray Marching 이란?

Ray Marching은 Ray Casting과 유사한 렌더링 기법입니다. 두 기법 모두 광선을 쏘아서 광선이 충돌하는 지점의 색상을 계산합니다.

하지만 Ray Casting은 광선이 충돌하는 지점을 정확히 계산하지만 Ray Marching은 SDF 함수를 이용하여 한 스텝에 충돌하지 않을 만큼 안전하게 전진할 수 있는 거리를 계산하고 해당 거리만큼 전진합니다. 이 과정을 반복하여 광선이 충돌하는 지점을 찾습니다.

단순히 생각했을 때 Ray Casting보다 비효율적이지만 한정적인 분야에서서 유용하게 쓰일 수 있습니다. 예를 들어 볼륨감이 있는 현실적인 구름을 표현하거나 3D 프랙탈을 렌더링할 때 유용합니다.

Shadertoy에서 레이마칭으로 구현된 작품들을 많이 찾아볼 수 있습니다.

Signed Distance Function (SDF) 이란?

보통 그래픽스에서는 모델을 Vertex와 Face로 표현하지만 레이마칭에서 모델들은 Signed Distance Function (SDF) 으로 표현됩니다.

SDF는 3차원 공간의 점 p 를 인자로 받아 해당 점과 모델의 표면 사이의 거리를 알려주는 함수입니다. 반환값이 양수이면 점이 외부에 있다는 것이고 음수이면 내부에 있다는 뜻입니다.

예를 들어 원점에 위치한 반지름이 r 인 구의 SDF는 length(p) - r 로 쓸 수 있습니다. length(p) 는 점 p 와 원점 사이의 거리이므로 length(p) - r 은 점 p 와 구의 표면 사이의 거리가 됩니다.

SDF는 꼭 최단 거리를 반환할 필요가 없습니다. 거리의 부호와 안전하게 지나갈 수 있는 거리만 알려준다면 언젠가는 물체에 도달하기 때문입니다. 단, 최단 거리와 가까울수록 더 빠르게 충돌하는 위치를 구할 수 있습니다.

코드 작성 셋팅

이 글에선 레이마칭 코드를 작성하기 위해 Shadertoy 사이트를 이용합니다. 또는 VSCode 에서 Shader Toy Extension 을 설치하여 웹을 이용하지 않는 방법도 있습니다.

코드를 작성하는 언어는 GLSL 이고 Fragment Shader만 작성합니다.

Shadertoy에서 New를 선택하면 아래의 샘플 코드와 함께 플레이어가 보입니다.

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
  // Normalized pixel coordinates (from 0 to 1)
  vec2 uv = fragCoord/iResolution.xy;

  // Time varying pixel color
  vec3 col = 0.5 + 0.5*cos(iTime+uv.xyx+vec3(0,2,4));

  // Output to screen
  fragColor = vec4(col,1.0);
}

mainImage 함수는 픽셀 좌표 fragCoord 을 입력으로 받아 픽셀의 색상 fragColor 을 출력하는 GLSL Fragment Shader입니다.

이외에도 스크린의 크기 iResolution, 지난 시간 iTime, 마우스 위치 iMouse 등등 이 입력으로 주어집니다.

Ray Marching 구현하기

모델 & 씬

간단하게 구와 평면 모델의 SDF 작성해봅시다. 앞서 설명했듯이 구의 SDF는 length(p) - r 이고 y = h 에 위치하고 y축에 수직인 평면은 p.y - h 로 쓸 수 있습니다.

float sdSphere(vec3 p, float s) {
  return length(p) - s;
}

float sdPlane(vec3 p, float h) {
  return p.y - h;
}

그렇다면 여러개의 모델이 있는 하나의 장면의 SDF는 어떻게 될까요? SDF 함수는 안전하게 전진할 수 있는 거리를 알려주어야 하므로 그 중 최소 거리를 택해야 합니다. 추가로 거리 뿐만 아니라 제일 가까운 모델의 id도 반환하도록 하였습니다. 이는 나중에 서로 다른 모델에 다른 색상을 입히기 위함입니다.

vec2 sdScene(vec3 p) {
  float id = 0.;
  float d = sdPlane(p, -1.);

  float dSphere = sdSphere(p, 1.0);
  if (dSphere < d) {
    d = dSphere;
    id = 1.;
  }

  return vec2(d, id);
}

raymarch()

레이마칭의 핵심이 되는 raymarch 함수는 광선의 시작점 p 와 광선의 방향인 dir 를 인자로 받아 광선이 충돌한 지점을 반환하는 함수입니다.

sdScene(p) 으로 매번 안전하게 전진할 수 있는 거리를 구해주고 해당 거리만큼 전진합니다. 이 과정을 반복하여 광선이 충돌하는 지점을 찾습니다.

하지만 광선이 영원히 총돌하지 않는 경우가 있을 수 있습니다. 광선이 무한히 나가는 것을 막기위해 광선이 출발점으로부터 일정 거리 MAX_DIST 를 넘으면 종료하도록 합니다.

또한, 광선이 언젠가 충돌하기는하나 광선이 지나는 공간이 매우 좁아 반복문이 오래돌면 전체적인 퍼포먼스에 지장을 줄 수 있습니다. 최대 반복횟수 MAX_STEP 을 정하여 반복문이 오래 도는 것을 막아주도록 합니다.

#define MAX_STEP 300
#define MAX_DIST 500.0
#define EPS 0.0001
vec4 raymarch(vec3 p, vec3 dir) {
  float id = -1.0;
  for (int st=0; st < MAX_STEP; st++) {
    vec2 d_id = sdScene(p);
    float d = d_id.x;

    if (d < EPS) return vec4(p, d_id.y);
    if (d > MAX_DIST) return vec4(p, -1.);
    p += d * dir;
  }

  return vec4(p, -2.);
}

카메라

각 픽셀에 대해 광선을 쏘기 위해서는 광선의 시작점과 방향을 알아야 합니다. viewMatrix 를 이용하여 카메라 좌표계에서 월드 좌표계로 변환해주고 rayDirection 함수를 이용하여 광선의 방향을 구해줍니다. 이후 raymarch 함수를 이용하여 광선이 충돌하는 지점을 구해줍니다.

mat3 viewMatrix(vec3 eye, vec3 center, vec3 up) {
  vec3 f = normalize(center - eye);
  vec3 s = normalize(cross(f, up));
  vec3 u = cross(s, f);
  return mat3(s, u, -f);
}

vec3 rayDirection(float fieldOfView, vec2 size, vec2 fragCoord) {
  vec2 xy = fragCoord - size / 2.0;
  float z = size.y / tan(radians(fieldOfView) / 2.0);
  return normalize(vec3(xy, -z));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
  vec3 viewDir = rayDirection(45.0, iResolution.xy, fragCoord);
  float eyeTheta = iMouse.x / iResolution.x * 3.14159 * 2.;
  float eyePhi = iMouse.y / iResolution.y + 0.9;
  float eyeD = 15.0;

  vec3 eye = vec3(cos(eyeTheta)*sin(eyePhi)*eyeD, cos(eyePhi)*eyeD, sin(eyePhi)*sin(eyeTheta)*eyeD);
  mat3 viewToWorld = viewMatrix(eye, vec3(0.0, 0.0, 0.0), vec3(0.0, 1.0, 0.0));
  vec3 worldDir = viewToWorld * viewDir;

  vec4 pos_id = raymarch(eye, worldDir);
  ...
}

그림자

표면의 한 점과 광원 사이에 물체가 없다면 해당 점은 광원에서 비춰지지 않기 때문에 그림자의 영역이 됩니다. 이를 구현하기 위해서는 광원에서 표면의 한 점으로 광선을 쏘아서 다른 물체와 충돌하는지 검사해 주어야 합니다. 충돌한다면 그 점은 광원에서 비춰지지 않는다는 뜻이므로 그림자가 생기게 됩니다.

float getShadow(vec3 ro, vec3 rd, float mint, float maxt) {
  for (float t=mint; t<maxt;) {
    float h = sdScene(ro + rd*t);
    if (h < 0.001) return 0.0;
    t += h;
  }
  return 1.0;
}

광원

광원 모델로는 Phong Reflection Model을 사용하겠습니다. Phong Reflection은 Ambient, Diffuse, Specular 세가지 요소로 이루어져 있습니다.

Ambient는 표면에 직접 빛이 비춰지지 않아도 반사광에 의해서 기본적으로 깔려있는 빛입니다. Diffuse는 표면에 직접 빛이 비춰지는 빛입니다. 빛의 방향과 물체 표면의 방향이 같을수록 더 밝게 비춰집니다. Specular은 표면에 정반사되어 카메라에 비춰지는 빛입니다. 반사되어 들어오는 빛의 방향과 카메라의 방향이 같을수록 더 밝게 비춰집니다.

Diffuse 값을 구할려면 Normal Vector를 꼭 구해야 합니다. 단순한 SDF 라면 식으로 구할 수 있겠지만 복잡한 SDF의 경우 Gradient 을 구하면 되는데 Gradient는 표면에서 Normal Vector와 방향이 같기 때문에 Gradient 를 구한 다음에 normalize 해주면 구할 수 있습니다.

실제론 정확한 Gradient 을 구할 수는 없고 미분 근사치를 구하는 방식과 동일하게 구현해주어야 합니다. 적당한 EPSILON 을 정해주고 p 에서 EPSILON 만큼 떨어진 두 점을 구한 뒤 normalize 해주면 됩니다.

#define EPSILON 0.001
vec3 getNormal(vec3 p) {
  return normalize(vec3(
    sdScene(vec3(p.x + EPSILON, p.y, p.z)) - sdScene(vec3(p.x - EPSILON, p.y, p.z)),
    sdScene(vec3(p.x, p.y + EPSILON, p.z)) - sdScene(vec3(p.x, p.y - EPSILON, p.z)),
    sdScene(vec3(p.x, p.y, p.z + EPSILON)) - sdScene(vec3(p.x, p.y, p.z - EPSILON))
  ));
}

간단하게 광원은 하나로 두겠습니다. 충돌한 곳의 id를 기준으로 색상이나 반사율을 다르게 해줄 수 있습니다. 이후 Phong Reflection Model의 식을 사용하여 색상을 구해줍니다.

vec3 getLight(vec4 p_id, vec3 lightPos, vec3 eye) {
  vec3 p = p_id.xyz;
  float id = p_id.w;

  vec3 col = vec3(1.0);
  if (id == 0.) { // 체크 무늬 바닥
    float t = floor(0.5*p.x) + floor(0.5*p.z);
    if( mod(t, 2.) < 1.0 ) col *= vec3(.5);
  } else if(id == 1.) {
    col *= vec3(1., .8, 1.);
  }

  vec3 lightCol = vec3(1.0);
  vec3 l = normalize(lightPos - p);
  vec3 n = getNormal(p);
  vec3 shad = vec3(getShadow(p + l*0.01, l, 0.001, 100.0));

  vec3 amb = vec3(0.3);
  float dif = max(dot(n, l), 0.0);
  vec3 ref = reflect(-l, n);
  float spe = .0; // pow(max(dot(ref, normalize(eye-p)), 0.0), 16.0);

  col *= amb + lightCol * (dif + spe) * shad;

  return col;
}

Primitive Model과 Operation

구와 평면 이외에도 몇가지 Primitive Model과 Operation 을 알아봅시다.

앞서 설명했듯이, SDF는 정확한 거리가 아닌 안전하게 전진할 수 있는 거리를 알려주어도 잘 작동하지만 차후에 여러가지 연산을 적용할 때 잘 작동하지 않을 수 있어서 Primitive Model의 SDF는 가능한 정확한 거리를 알려주는 것이 좋습니다.

서술할 것들 외에도 더 많은 Primitive Model과 Operation이 있습니다.

Box

Box는 원점대칭이므로 sdBox(p) = sdBox(abs(p)) 입니다. 따라서 좌표가 모두 양수인 공간에서의 SDF만 구해주면 됩니다.

Box의 꼭짓점이 vec3 b 에 있다고 하면 이제 3가지 경우로 나눌 수 있습니다.

  1. p.x > b.x && p.y > b.y && p.z > b.z
  2. p.x < b.x && p.y < b.y && p.z < b.z
  3. 나머지

1번과 3번은 박스 외부의 영역이고 2번은 내부의 영역입니다. 쉽게 구하기 위하여 b 를 원점으로 옮겨서 생각해봅시다. p 에서 b 를 뺀 좌표를 q 라고 하면..

1번의 경우 q 와 원점 사이의 거리가 됩니다.

3번의 경우 q 가 원점보다 더 가까운 곳이 있습니다. 음수가 되는 좌표가 있다는 뜻인데 해당 좌표를 0으로 옮겨주면 1번과 동일하게 처리가 됩니다.

2번의 경우 좌표가 모두 음수이고 최단 거리는 축에 평행한 선분이 됩니다. 모든 축 중 최소가 되는 거리 max(q.x, q.y, q.z) 가 됩니다.

float sdBox(vec3 p, vec3 b) {
  vec3 q = abs(p) - b;
  return length(max(q, 0.0)) + // outer
    min(max(q.x, max(q.y, q.z)), 0.0); // inner
}

Cylinder

y축으로 세워진 원통은 y축을 포함한 평면으로 잘랐을 때 단면은 모두 모양이 같습니다. 반지름이 r , 높이가 h/2 라고하면 b = (r, h) 을 꼭짓점으로 하는 2차원 sdBox와 같습니다. b = (length(p.xz), p.y)

float sdCylinder(vec3 p, vec2 h) {
  vec2 d = abs(vec2(length(p.xz),p.y)) - h;
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

연산 적용 순서

원점이 중점인 큐브를 회전 후 이동 시킨다고 합시다. 월드_좌표계 = 모델_연산(=이동연산 * 회전연산) * 모델_좌표계 이므로 일반적으로 Vertex과 Face로 이뤄진 모델이 있고 해당 식을 사용하여 월드 좌표계로 변환하는 겁니다.

하지만 Ray Marching에서는 상황이 다릅니다. 알고 있는 것은 현재 광선의 위치는 월드 좌표계에 대한 것이고 모델(=SDF)을 계산하고자 한다면 모델 좌표계로 변환해야 합니다. 따라서 월드좌표계를 `모델좌표계 = 모델연산의역연산(=회전연산의 역연산 이동연산의_역연산) 월드_좌표` 로 변환 후 모델 좌표계를 SDF에 넣어 계산해야 합니다.

이동과 회전

+v 만큼 이동하는 것의 역연산은 -v 만큼 이동하는 것이므로 p-v 를 반환해주면 됩니다.

vec3 opTranslate(vec3 p, vec3 v) {
  return p-v;
}

회전변환의 경우도 +theta 만큼 회전 시키는 것의 역연산은 -theta 만큼 회전시키는 것이므로 p 를 회전시킨 뒤에 -r 만큼 회전시켜주면 됩니다.

mat2 rotateMat(float a) {
  float s = sin(a), c = cos(a);
  return mat2(c, -s, s, c);
}

vec3 opRotate(vec3 p, vec3 r) {
  p.xy *= rotateMat(-r.z);
  p.zx *= rotateMat(-r.y);
  p.yz *= rotateMat(-r.x);
  return p;
}

y축을 기준으로 위아래로 움직이며 회전하는 큐브는

// p는 광선의 위치
vec3 boxPos = vec3(0.0, cos(iTime)*0.5 + 1.5, 0.0);
vec3 boxRot = vec3(iTime, iTime*1.5, .0);
vec3 boxP = opTranslate(p, boxPos);
boxP = opRotate(boxP, boxRot);
float dBox = sdBox(boxP, .7);

로 구현할 수 있습니다.

연산의 순서가 거꾸로 적용되는 것을 주의해야 합니다. 행렬을 사용하여 역행렬을 구하는 것이 더 편리할 수 있습니다.

Boolean Operations

수학의 집합과 유사하게 모델들을 합집합, 교집합, 차집합으로 연산할 수 있습니다.

합집합은 둘 중 하나라도 부딧치지 않도록 두 개의 거리 중 최소를 선택합니다.

float opUnion(float da, float db) {
  return min(da, db);
}
float d = opUnion(
  sdBox(p, vec3(1.0)),
  sdSphere(opTranslate(p, vec3(1.0)), 1.0);
);

교집합은 현재 위치가 하나의 모델 안에 있더라도 다른 모델의 표면까지 갈 수 있도록 두 개의 거리 중 최대를 선택합니다.

float opIntersection(float da, float db) {
  return max(da, db);
}

차집합은 SDF는 모델의 내부에서는 표면과의 거리를 음수로 가집니다. 따라서 거리에 - 를 붙여주는 것이 내부와 외부 바꾸는 역할을 합니다. A - B = A ∩ B^c 이므로 B 의 내부와 외부를 바꾼뒤에 교집합을 구하면 차집합이 됩니다.

float opDifference(float da, float db) {
  return max(da, -db);
}

Boolean Operation엔 Smooth 버젼이 있는데 이는 두 모델이 겹치는 부분을 부드럽게 만들어줍니다. 여기 에서 기술적인 내용을 확인할 수 있습니다.

float opSmoothUnion(float d1, float d2, float k) {
  float h = clamp(0.5 + 0.5*(d2-d1)/k, 0.0, 1.0);
  return mix(d2, d1, h) - k*h*(1.0-h);
}

종합

지금까지의 코드들을 종합해보면 아래와 같습니다.

float sdSphere(vec3 p, float radius) {
  return length(p) - radius;
}

float sdPlane(vec3 p, float height) {
   return p.y - height;
}

float sdBox(vec3 p, float b) {
  vec3 d = abs(p) - b;
  return length(max(d,0.0))
         + min(max(d.x,max(d.y,d.z)),0.0);
}

vec3 opTranslate(vec3 p, vec3 v) {
  return p-v;
}

mat2 rotateMat(float a) {
  float s = sin(a), c = cos(a);
  return mat2(c, -s, s, c);
}

vec3 opRotate(vec3 p, vec3 r) {
  p.xy *= rotateMat(-r.z);
  p.zx *= rotateMat(-r.y);
  p.yz *= rotateMat(-r.x);
  return p;
}

float opSmoothUnion(float d1, float d2, float k) {
    float h = clamp(0.5 + 0.5*(d2-d1)/k, 0.0, 1.0);
    return mix(d2, d1, h) - k*h*(1.0-h);
}

vec2 sdScene(vec3 p) {
  float id = 0.;
  float d = sdPlane(p, -1.);

  float dSphere = sdSphere(p, 1.0);
  vec3 boxPos = vec3(0.0, cos(iTime)*0.5 + 1.5, 0.0);
  vec3 boxRot = vec3(iTime, iTime, .0);
  vec3 boxP = opTranslate(p, boxPos);
  boxP = opRotate(boxP, boxRot);
  float dBox = sdBox(boxP, .7);
  float dUnion = opSmoothUnion(dSphere, dBox, 1.0);
  if (dUnion < d) {
    d = dUnion;
    id = 1.;
  }

  return vec2(d, id);
}

#define EPSILON 0.0001
vec3 getNormal(vec3 p) {
  return normalize(vec3(
    sdScene(vec3(p.x + EPSILON, p.y, p.z)).x - sdScene(vec3(p.x - EPSILON, p.y, p.z)).x,
    sdScene(vec3(p.x, p.y + EPSILON, p.z)).x - sdScene(vec3(p.x, p.y - EPSILON, p.z)).x,
    sdScene(vec3(p.x, p.y, p.z  + EPSILON)).x - sdScene(vec3(p.x, p.y, p.z - EPSILON)).x
  ));
}


#define MAX_STEP 300
#define MAX_DIST 500.0
#define EPS 0.0001
vec4 raymarch(vec3 p, vec3 dir) {
  float id = -1.0;
  for (int st=0; st < MAX_STEP; st++) {
    vec2 d_id = sdScene(p);
    float d = d_id.x;

    if (d < EPS) return vec4(p, d_id.y);
    if (d > MAX_DIST) return vec4(p, -1.);
    p += d * dir;
  }

  return vec4(p, -2.);
}

float getShadow(vec3 ro, vec3 rd, float mint, float maxt) {
  for (float t = mint; t < maxt;) {
    float h = sdScene(ro + rd*t).x;
    if (h < 0.001) return 0.0;
    t += h;
  }
  return 1.0;
}

vec3 getLight(vec4 p_id, vec3 lightPos, vec3 eye) {
  vec3 p = p_id.xyz;
  float id = p_id.w;

  vec3 col = vec3(1.0);
  if (id == 0.) {
    float t = floor(0.5*p.x) + floor(0.5*p.z);
    if( mod(t, 2.) < 1.0 ) col *= vec3(.5);
  } else if(id == 1.) {
    col *= vec3(1., .8, 1.);
  }

  vec3 lightCol = vec3(1.0);
  vec3 l = normalize(lightPos - p);
  vec3 n = getNormal(p);
  vec3 shad = vec3(getShadow(p + l*0.01, l, 0.001, 100.0));

  vec3 amb = vec3(0.3);
  float dif = max(dot(n, l), 0.0);
  vec3 ref = reflect(-l, n);
  float spe = .0; // pow(max(dot(ref, normalize(eye-p)), 0.0), 16.0);

  col *= amb + lightCol * (dif + spe) * shad;

  return col;
}

mat3 viewMatrix(vec3 eye, vec3 center, vec3 up) {
  vec3 f = normalize(center - eye);
  vec3 s = normalize(cross(f, up));
  vec3 u = cross(s, f);
  return mat3(s, u, -f);
}

vec3 rayDirection(float fieldOfView, vec2 size, vec2 fragCoord) {
  vec2 xy = fragCoord - size / 2.0;
  float z = size.y / tan(radians(fieldOfView) / 2.0);
  return normalize(vec3(xy, -z));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {
  vec3 viewDir = rayDirection(45.0, iResolution.xy, fragCoord);
  float eyeTheta = iMouse.x / iResolution.x * 3.14159 * 2.;
  float eyePhi = iMouse.y / iResolution.y + 0.9;
  float eyeD = 15.0;

  vec3 eye = vec3(cos(eyeTheta)*sin(eyePhi)*eyeD, cos(eyePhi)*eyeD, sin(eyePhi)*sin(eyeTheta)*eyeD);
  mat3 viewToWorld = viewMatrix(eye, vec3(0.0, 0.0, 0.0), vec3(0.0, 1.0, 0.0));
  vec3 worldDir = viewToWorld * viewDir;

  vec3 col = vec3(0.0);
  vec4 pos_id = raymarch(eye, worldDir);
  vec3 light = getLight(pos_id, vec3(-10., 10.0, 10.), eye);

  col = light;

  fragColor = vec4(col,1.0);
}