Sangil Lee

Sangil Lee

Staff Engineer in S. LSI, Samsung

GLSL에서 간단한 레이 마칭 (Ray Marching) 구현하기

0

레이 마칭 (Ray marching) 은 3D 장면을 2D 화면에 렌더링하는 기법 중 하나다. 이름에서 알 수 있듯이, 광선(ray) 묶음이 카메라에서 점진적으로 전진하며 물체 또는 광원과의 교차점을 탐색하는 기법이다. 이 포스트에서는 GLSL로 레이 마칭을 단계별로 구현하는 방법을 설명한다. 거리 맵(distance map), 법선 맵(normal map)을 계산하는 법과, 이를 이용하여 조명, 셰이딩, 반사 등을 렌더링하는 방법을 예제 코드와 함께 순서대로 다룬다.

이 포스트에서 최종적으로 구현하고자 하는 예제는 아래와 같다. 마우스나 터치를 통해 카메라 앵글을 조절할 수 있다.

Touch or hover your mouse here

부호 있는 거리 함수 (SDF; Signed Distance Function)

GLSL에서 3D 모델은 SDF로 표현할 수 있다. SDF는 3차원 공간 상의 임의의 점에서 3D 모델 표면까지의 최소 거리를 반환한다. 점이 모델 외부에 있으면 양수, 내부에 있으면 음수를 반환하며, 표면 위의 점에서는 0이 된다. 카메라에서 출발한 ray를 따라 SDF 값을 계산하면 3D 모델의 표면 위치를 구할 수 있다.

아래 예시 코드는 s.xyz 점을 중심으로 하는 반지름 s.w의 구(sphere)와 XZ-평면(plane) 두 객체의 SDF를 정의한다. 두 물체를 고려한 최종 SDF 반환값 d는 sphereDist와 planeDist의 최솟값이다. 가까운 물체의 표면에 ray가 도달하면, 더 이상 전진할 필요가 없기 때문이다. 이렇게 하면 두 모델의 SDF를 하나로 합칠 수 있다. 다양한 모델의 SDF 함수는 Inigo Quilez에서 더 찾아볼 수 있다.

float SDF(vec3 p) 
{
  vec4 s = vec4(0,1.2,0,1); //Sphere xyz is position w is radius
  float sphereDist = length(p-s.xyz) - s.w;
  float planeDist = p.y;
  float d = min(sphereDist,planeDist);
  return d;
}

변형 (Deformation)

사막이나 산악 지형처럼 울퉁불퉁한 지형을 만들고 싶다면, 노이즈 함수 (Noise functions 참고) 와 프랙탈 브라운 운동 (fBm) 기법 (Shader patterns 참고) 을 활용하면 쉽게 구현할 수 있다.

float get_plane_dist(in vec3 p) {
  // return p.y;
  return p.y + 0.5*noise(2.*p.xz); // and so on ..
}

레이 마칭

어느 한 점에서 출발한 ray와 어느 물체 간의 SDF를 이용하면, ray의 출발점과 물체 사이의 거리를 계산할 수 있다. 이 때, ray는 원점 ro와 방향 rd로 정의된다. ray의 원점은 카메라의 위치이며, 방향은 카메라 내부 파라미터를 통해 정규화된 이미지 좌표와 연결된다.

각 ray에 대해 ray 위의 점을 샘플링하여, 해당 점이 물체의 표면 외부에 있는지 확인한다. 표면 또는 내부가 아니라면 ray의 진행 방향으로 한 걸음 나아가 다시 샘플링한다. 이 때 ray는 고정된 값으로 전진하는 것이 아니라, SDF 값만큼 전진한다. SDF 값이 가장 가까운 어느 표면까지의 최소 거리를 나타내므로, 현재 점을 중심으로 SDF 값을 반지름으로 하는 구 내부에는 어떠한 물체도 없음이 보장된다.

아래와 같이 for loop 내에서 샘플링된 점이 표면에 임계치 내로 도달하면 ray가 그동안 전진한 거리 d를 반환한다. 렌더링 성능을 위해, 최대 반복 횟수에 도달하거나 깊이가 최대 거리를 초과하면 1./0.을 반환하도록 한다.

float ray_march(vec3 ro, vec3 rd) {
  float d = 0.; // init depth
  for (int i=0; i < MAX_STEPS; ++i) {
    vec3 p = ro + rd * d;
    float ds = SDF(p);
    
    if (ds < SURF_DIST) 
      return d;
      
    d += ds;
  
    if (d > MAX_DIST)
      return 1./0.;
  }
  return 1./0.;
}
평평한 표면 울퉁불퉁한 표면

법선 맵 (Normal Map)

물체에 조명과 셰이딩 효과를 적용하려면 법선 맵이 필요하다. 3D 물체 표면에서의 법선 벡터를 구하기 위해, 해당 점 주변 국소 영역에서 SDF의 변화량을 계산한다. 정의상 SDF 값은 표면에서 멀어질수록 증가하므로, 표면 근처에서 SDF의 기울기(gradient)가 곧 표면 법선 벡터가 된다. 아래 코드와 같이 충분히 작은 변위를 이용해 수치적으로 기울기를 계산할 수 있다.

vec3 get_normal(vec3 p) { 
  vec2 e = vec2(.01,0); // Epsilon
  vec3 n = vec3(
    SDF(p+e.xyy) - SDF(p-e.xyy),
    SDF(p+e.yxy) - SDF(p-e.yxy),
    SDF(p+e.yyx) - SDF(p-e.yyx)
  );
  
  return normalize(n);
}

조명과 셰이딩

깊이 맵과 법선 맵을 구했으면, Create realistic Earth에서 설명한 것과 동일한 기법을 적용할 수 있다. 해당 포스트의 “Mountain Shadow” 섹션에서는 아래와 같이 조명과 셰이딩, mixAmountTexture, 을 계산하였다.

// Normal map texture
vec3 t_normal = texture2D( u_normalTexture, vUv ).xyz * 2.0 - 1.0;
vec3 normal = normalize(vTbn * t_normal);
float cosAngleSunToSurface = dot(normal, sunDir); // Compute cosine sun to normal
mixAmountTexture *= 1.0 + u_normalPower * (cosAngleSunToSurface - cosAngleSunToNormal);
mixAmountTexture = clamp(mixAmountTexture, 0., 1.);

위와 유사하게, 쉐이딩를 계산하기 위해 광원 방향 벡터와 법선 벡터의 내적을 수행한다. 주변광을 모사하기 위해 clamp의 두 번째 인수를 0.2로 설정한다. 따라서 최소 밝기의 계수는 0.2가 된다.

occlusion에 의한 그림자는 물체 표면에서 출발하여 광원 방향 ld 으로 전진하는 ray marching 을 계산함으로써 감지할 수 있다. ray marching 으로 계산된 거리값이 표면 점과 광원 사이의 실제 거리보다 작다면,이광원과 표면 사이에 빛을 가리는 무언가가 있다는 의미이다.

실제 occlusion을 적용하는 코드에서는 두 가지 기술적인 트릭을 사용했다. 첫째, 물체 표면의 점 p는 이미 SDF(p) < SURF_DIST 조건을 만족하므로, ld 방향으로 ray marching을 실행하면, 함수가 즉시 종료된다. 따라서, n * SURF_DIST 크기의 법선 벡터를 더해 약간의 오프셋을 준다. 둘째, 광원을 등진 구의 표면과 같이, 표면의 점이 셰이딩과 occlusion을 동시에 겪을 수 있는 것을 방지하기 위해 d > 10. * SURF_DIST 조건을 추가한다. 이 조건을 추가함으로써, 표면이 자기가 속한 물체에 의해 occlusion 되는 것을 방지할 수 있다.

float light(vec3 p, vec3 n) {
  // lighting
  vec3 l = vec3(5.*sin(u_time),5.,5.0*cos(u_time)); // Light Position
  vec3 ld = normalize(l-p); // Light Vector
  
  float dif = dot(n,ld); // Diffuse light
  dif = clamp(dif,0.2,1.); // Clamp so it doesnt go below 0.2
  
  // occlusion
  float d = ray_march(p+n*SURF_DIST, ld); 
  if ( d > 10. * SURF_DIST && d < length(l-p)) dif *= 0.3;
  return dif;
}

점 광원 정반사

GLSL에서 정반사를 렌더링하는 방법은 Create realistic Earth의 “Ocean Reflection” 섹션에서도 다루었다. 아래 코드는 해당 섹션에서 발췌하였다.

// Specular map texture with reflection
vec3 surfacePosition = u_position + vPosition; // surface position in world coordinates
float reflectRatio = texture2D(u_specTexture, vUv).r;
reflectRatio = 0.3 * reflectRatio + 0.1;
vec3 reflectVec = reflect(-sunDir, normal); // reflected vector of sunlight

// dot product between reflected light and camera vector
float specPower = clamp(dot(reflectVec, normalize(cameraPosition - surfacePosition)), 0., 1.);
color += mixAmountTexture * pow(specPower, 2.0) * reflectRatio;

위와 유사하게 반사된 광원 벡터 r과 카메라 벡터의 내적을 계산한다. f는 위 코드의 reflectRatio와 동일한 역할로, 반사율을 제어한다.

float reflection(vec3 p, vec3 n, vec3 c, float f) {
  vec3 l = vec3(5.*sin(u_time),5.,5.0*cos(u_time)); // Light Position
  vec3 ld = normalize(l-p); // Light Vector
  vec3 r = reflect(-ld, n); // reflected vector of sunlight

  // dot product between reflected light and camera vector
  float s = clamp(dot(r, normalize(c - p)), 0., 1.); 
  return pow(s, 10.0) * f;
}

색상

여기서는 3D 렌더링의 기본적인 색상 적용 기법을 다룬다. 재질은 단색 RGB나 텍스처로 단순하게 표현하는 방법에 대해 알아본다.

단색 RGB

먼저 각 물체의 색상을 정의한다. 그런 다음 get_color(p - n*SURF_DIST)를 호출해 표면 안쪽의 색상을 추출함으로써 ray의 색상을 계산할 수 있다.

vec3 get_color(vec3 p) {
  vec4 s = get_sphere();
  if (length(p-s.xyz) < s.w) return vec3(0.5, 0.7, 1.0); // sphere color
  if (p.y < 0.0) return vec3(0.57, 0.42, 0.3); // ground color
  return vec3(0.);
}

텍스처

단색 대신 텍스처를 사용하면 물체를 더 사실적으로 표현할 수 있다. 무료 텍스처는 여기에서 찾을 수 있다. 텍스처를 적용하기 전에 UV 맵 \(f: \mathbb{R}^3 \rightarrow [0, 1]^2\)을 정의해야 한다. SDF가 \(SDF(x, y, z) = y\)인 평면의 경우, UV 맵은 fract(vec2(x, z))로 간단히 매핑되어 반복 패턴을 만들어낸다. 구의 경우, 중심 위치와 반지름을 이용해 경도와 위도로 UV 맵을 매핑하여야 한다.

vec2 get_sphere_tex_coord(in vec3 p) {
  vec4 s = vec4(0, 1.2, 0, 1);
  vec3 r = p - s.xyz;
  float phi = atan(sqrt(dot(r.xz, r.xz)), r.y);
  float theta = atan(r.z, r.x);
  return fract(vec2(phi / PI, theta / TWO_PI));
}

이를 추가하여, get_color 함수를 업데이트한다.

vec3 get_color(vec3 p) {
  vec4 s = get_sphere();
  if (length(p-s.xyz) < s.w) return vec3(0.5, 0.7, 1.0); // sphere color
  if (p.y < 0.0) return texture2D(u_texture, fract(0.2*p.xz)).rgb; // ground color
  return vec3(0.);
}

모래 텍스처를 사용해 사막 지형을 표현하면 다음과 같은 결과를 얻을 수 있다.

시점 변환

카메라의 위치와 자세를 결정하기 위해 회전 행렬과 이동 벡터로 구성된 변환 행렬을 사용한다. 단, 주의할 점은, numpy와 달리 원소값이 column major로 저장된다. 그래서 아래와 같이 개행으로 작성한 mat3 데이터의 실제 행렬은 보이는 행렬의 transpose로 생각해야 된다. X, Y, Z 축에 대한 회전 행렬은 다음과 같다.

mat3 Rx = mat3(
  1., 0., 0.,
  0., cos(psi), sin(psi),
  0., -sin(psi), cos(psi)
);
mat3 Ry = mat3(
  cos(theta), 0., -sin(theta),
  0., 1., 0.,
  sin(theta), 0., cos(theta)
);
mat3 Rz = mat3(
  cos(phi), sin(phi), 0.,
  -sin(phi), cos(phi), 0.,
  0., 0., 1.
);

회전 행렬 R과 이동 벡터 t가 주어지면, 변환 행렬 T는 다음과 같이 계산된다.

mat4 T = mat4(
  vec4(R[0], 0.0),
  vec4(R[1], 0.0),
  vec4(R[2], 0.0),
  vec4(t, 1.0)
);

이를 통해 카메라 원점과 시선 방향에 변환 행렬 T를 곱하면, 월드 좌표계에서의 ray 원점과 방향을 구할 수 있다.

vec3 co = vec3(0,0,0); // camera origin
vec3 cd = normalize(vec3(xy,1)); // camera ray vector

vec3 ro = (T * vec4(co, 1.)).xyz; // ray origin w.r.t world
vec3 rd = (T * vec4(cd, 0.)).xyz; // ray vector w.r.t world

정리하며

더 읽을거리

위에서 구현한 코드는 ray가 물체의 표면에 닿으면서 marching이 중단되기 때문에, 완전히 불투명한 물체에 대해서만 렌더링이 가능하다. 다음 포스트에서는 투명/반투명 재질의 색상 및 밀도를 participating media rendering 기법을 이용해 계산한다. 더불어, 굴절과 정반사 알고리즘을 일반화하고, 거친 표면에서의 난굴절과 난반사를 렌더링한다.

또한 현재까지 점 광원에 대해서만 다루었으므로, 추후에는 체적 광원(volumetric light source)에 대한 조명, 셰이딩, 반사 기법을 다룰 예정이다.

전체 코드

최종 fragment 코드는 아래와 같다.

// Author: Sangil Lee, https://sangillee.com
// Refer to Tim Coster, https://timcoster.com

precision highp float;
uniform vec2 u_resolution; // Width & height of the shader
uniform float u_time; // Time elapsed
uniform sampler2D u_texture;
uniform vec2 u_angle; // camera rotation

// Constants
#define PI 3.1415925359
#define TWO_PI 6.2831852
#define MAX_STEPS 200 // Mar Raymarching steps
#define MAX_DIST 100. // Max Raymarching distance
#define SURF_DIST .01 // Surface distance
#define SAMPLE_DIST 0.05 // Sample distance

float random(in vec2 x) {
  return fract(sin(dot(x, vec2(12.9898,54.233))) * 43758.5453123);
}

float noise(in vec2 x) {
  vec2 i = floor(x);
  vec2 f = fract(x);

  float tl = random(i); // top-left corner
  float tr = random(i + vec2(1.0, 0.0)); // top-right corner
  float bl = random(i + vec2(0.0, 1.0)); // bottom-left corner
  float br = random(i + vec2(1.0, 1.0)); // bottom-right corner

  vec2 u = smoothstep(0., 1., f);

  return
  mix(
    mix(tl, tr, u.x), 
    mix(bl, br, u.x), 
  u.y);
}

vec3 get_light() {
  return vec3(6.*sin(0.3*u_time),5.,6.*cos(0.3*u_time));
}

float get_plane_dist(in vec3 p) {
  return p.y + 0.5*noise(0.5*p.xz) + 0.05*noise(2.*p.xz) + 0.012*noise(4.*p.xz);
}

float get_sphere_dist(in vec3 p) { 
  vec4 s = vec4(0, 1.2, 0, 1);
  return length(p - s.xyz) - s.w;
}

vec2 get_sphere_tex_coord(in vec3 p) {
  vec4 s = vec4(0, 1.2, 0, 1);
  vec3 r = p - s.xyz;
  float phi = atan(sqrt(dot(r.xz, r.xz)), r.y);
  float theta = atan(r.z, r.x);
  return fract(vec2(phi / PI, theta / TWO_PI));
}

float SDF(in vec3 p) {
  float sphereDist = get_sphere_dist(p);
  float planeDist = get_plane_dist(p);
  float d = min(sphereDist, planeDist);
  return d;
}

int get_object_id(in vec3 p) {
  if (get_sphere_dist(p) < 0.0) return 1;
  if (get_plane_dist(p) < 0.0) return 2;
  return 0;
}

vec3 get_color(in vec3 p) {
  if (get_sphere_dist(p) < 0.0) return vec3(0.5, 0.7, 1.0); // sphere color
  // if (get_plane_dist(p) < 0.0) return vec3(0.8, 0.55, 0.3); // ground color
  // if (get_sphere_dist(p) < 0.0) return texture2D(u_texture, get_sphere_tex_coord(p)).rgb; // sphere color
  if (get_plane_dist(p) < 0.0) return texture2D(u_texture, fract(0.2*p.xz)).rgb; // ground color
  return vec3(0.);
}

float get_density(in vec3 p) {
  if (get_sphere_dist(p) < 0.0) return 0.2; // sphere density
  if (get_plane_dist(p) < 0.0) return 10.0; // ground density
  return 0.;
}

float get_specular(in vec3 p) {
  if (get_sphere_dist(p) < 0.0) return 0.3; // sphere specular
  if (get_plane_dist(p) < 0.0) return 0.05; // ground specular
  return 0.;
}

float ray_march(in vec3 ro, in vec3 rd) {
  float d = 0.;
  for (int i = 0; i < MAX_STEPS; ++i) {
    vec3 p = ro + rd * d;
    float ds = SDF(p);
    
    if (ds < SURF_DIST)
      return d;
      
    d += ds;
  
    if (d > MAX_DIST)
      return 1./0.;
  }
  return 1./0.;
}

vec3 get_normal(in vec3 p) { 
  vec2 e = vec2(.01,0); // Epsilon
  vec3 n = vec3(
    SDF(p+e.xyy) - SDF(p-e.xyy),
    SDF(p+e.yxy) - SDF(p-e.yxy),
    SDF(p+e.yyx) - SDF(p-e.yyx)
  );
  
  return normalize(n);
}

float compute_contrast(in vec3 p, in vec3 n) {
  // Directional light
  vec3 l = get_light(); // Light Position
  vec3 ld = normalize(l-p); // Light Vector
  
  float dif = dot(n,ld); // Diffuse light
  dif = clamp(dif,0.,1.); // Clamp so it doesnt go below 0
  
  // Shadows
  float d = ray_march(p + 10.*n*SURF_DIST, ld);
  if (d < length(l-p))
    dif *= 0.3;

  return dif;
}

float compute_reflection(in vec3 p, in vec3 n, in vec3 c, in float f) {
  vec3 l = get_light(); // Light Position
  vec3 ld = normalize(l-p); // Light Vector
  vec3 r = reflect(-ld, n); // reflected vector of sunlight
  float s = clamp(dot(r, normalize(c - p)), 0., 1.); // dot product between reflected light and camera vector
  return pow(s, 10.0) * f;
}

void main()
{
  vec3 color = vec3(0.);
    
  // image plane: u ~ [0, u_resolution.x], v ~ [0, u_resolution.y]
  vec2 uv = gl_FragCoord.xy;
    
  // intrinsic parameter
  vec2 f = vec2(600.);
  vec2 c = 0.5*u_resolution.xy;
  
  // normalized image plane
  vec2 xy = (uv - c) / f;
  
  // distortion
  vec2 k_d = vec2(0.3,0.0);
  vec2 p_d = vec2(0.0,0.0);
  float r2 = dot(xy, xy);
  xy = (1. + k_d.x*r2 + k_d.y*r2*r2) * xy + vec2(2.*p_d.x*xy.x*xy.y + p_d.y*(r2 + 2.*xy.x*xy.x), 2.*p_d.y*xy.x*xy.y + p_d.x*(r2 + 2.*xy.y*xy.y));
  
  // camera
  float theta = u_angle.x; // camera yaw
  float phi = u_angle.y; // camera pitch

  mat3 Ry = mat3(
    cos(theta), 0., -sin(theta),
    0., 1., 0.,
    sin(theta), 0., cos(theta)
  );

  mat3 Rx = mat3(
    1., 0., 0.,
    0., cos(phi), sin(phi),
    0., -sin(phi), cos(phi)
  );

  mat3 R = Ry * Rx; // camera rotation
  float r = 6.0; // camera distance
  vec3 t = vec3(-r*cos(phi)*sin(theta), r*sin(phi), -r*cos(phi)*cos(theta)); // camera translation

  mat4 T = mat4(
    vec4(R[0], 0.0),
    vec4(R[1], 0.0),
    vec4(R[2], 0.0),
    vec4(t, 1.0)
  );

  vec3 co = vec3(0,0,0); // camera origin
  vec3 cd = normalize(vec3(xy,1)); // camera ray vector
  
  // ray
  vec3 ro = (T * vec4(co, 1.)).xyz; // ray origin w.r.t world
  vec3 rd = (T * vec4(cd, 0.)).xyz; // ray vector w.r.t world

  // find the surface point
  float d = ray_march(ro, rd);
  vec3 p = ro + rd * d;
  vec3 n = get_normal(p);

  color = get_color(p - n*SURF_DIST); // apply color and texture
  color *= compute_contrast(p, n); // apply lighting and shading
  
  // apply specular lighting
  color += compute_reflection(p, n, ro, get_specular(p-n*SURF_DIST));
  
  gl_FragColor = vec4(color,1.0);
}