Ray Marching - 1
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가지 경우로 나눌 수 있다.
p.x > b.x && p.y > b.y && p.z > b.zp.x < b.x && p.y < b.y && p.z < b.z- 나머지
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);
}