고리가 있는 사실적인 토성 만들기

태양계에는 크게 두 종류의 행성이 있다. 암석형 행성과 가스형 행성이며, 일부는 고리를 가지고 있다. 그 중에서도 토성은 뚜렷한 고리를 가진, 잘 알려진 행성 중 하나이다. 이 글에서는 고리, 표면, 그리고 고리와 행성 사이에 드리워지는 그림자를 상세히 표현하여 사실적인 토성 3D 모델을 만들어본다.

이 글은 태양계 시뮬레이터 시리즈의 일부다:

1. Three.js 좌표계 기초
2. Three.js PBR (물리 기반 렌더링) 재질 기초
3. 사실적인 지구 만들기
4. 사실적인 태양 만들기
5. 고리가 있는 행성 만들기
6. 불규칙한 형태의 위성 만들기
7. 태양을 빛나게 하기
8. 은하수 스카이박스 만들기
9. 타원 궤도 계산하기

토성의 본체

토성의 시각화는 이전 글¹에서 시작한다. 다만 토성은 가스 거성이므로 법선 텍스처, 반사율 텍스처, 구름은 렌더링할 필요가 없다. 따라서 토성의 본체는 아래와 같이 구현할 수 있다.

// vertex shader of Saturn body
varying vec2 vUv;
varying vec3 vNormal;
varying vec3 vPosition;
varying mat3 vTbn;

attribute vec4 tangent; // "geometry.computeTangents()" is needed.

void main() {
  vUv = uv;
  vNormal = normalize(mat3(modelMatrix) * normal);
  vPosition = mat3(modelMatrix) * position;
  
  gl_Position = projectionMatrix * modelViewMatrix * vec4(position, 1.0);

  vec3 t = normalize(tangent.xyz);
  vec3 n = normalize(normal.xyz);
  vec3 b = normalize(cross(t, n));

  t = mat3(modelMatrix) * t;
  b = mat3(modelMatrix) * b;
  n = mat3(modelMatrix) * n;
  vTbn = mat3(t, b, n);
}

uniform sampler2D u_dayTexture;
uniform sampler2D u_normalTexture;
uniform vec3 u_sunRelPosition; // the relative position of light source
uniform vec3 u_position; // the position of this object
uniform float u_sunRadius; // the radius of light source
uniform float u_specRatio;

varying vec2 vUv;
varying vec3 vNormal;
varying vec3 vPosition;
varying mat3 vTbn;

void main( void ) {
  vec3 sunDir = normalize(u_sunRelPosition);

  // 1. Day and night texture
  vec3 dayColor = texture2D( u_dayTexture, vUv ).rgb;
  float cosAngleSunToNormal = dot(vNormal, sunDir); // Compute cosine sun to normal
  float hemisphereLight = 1. / (1. + exp(-20. * (cosAngleSunToNormal+0.05)));
  float mixAmountSurface = hemisphereLight; // Sharpen the edge beween the transition
  vec3 color = mix( vec3(0.), dayColor, mixAmountSurface ); // Select day or night texture based on mix.

  // 2. Specular map texture with reflection
  vec3 surfacePosition = u_position + vPosition;
  vec3 reflectVec = reflect(-sunDir, vNormal); // reflected vector of sunlight
  float specPower = dot(reflectVec, normalize(cameraPosition - surfacePosition)); // dot product between reflected light and camera vector
  color += mixAmountSurface * pow(specPower, 3.0) * u_specRatio * 0.4;

  gl_FragColor = vec4( color, 1.0 );
}

토성의 고리

이제 토성의 나머지 부분인 고리를 구현해 보자. 토성의 고리는 다양한 색상과 밀도를 가진 무수히 많은 작은 입자들로 이루어져 있다. 모든 입자를 렌더링하지 않더라도, 양면 디스크(double-sided disc) 오브젝트를 생성하여 고리의 각 세부 구간의 색상과 투명도를 표현할 것이다.

const geometry_ring = new THREE.RingGeometry(1.3, 2.2, 128).rotateX(-Math.PI / 2);
const material_ring = new THREE.MeshBasicMaterial({color: 0xffffaa});

const ring_upper = new THREE.Mesh(geometry_ring, material_ring);
const ring_lower = new THREE.Mesh(geometry_ring, material_ring).rotateX(Math.PI);;

saturn.add(ring_upper);
saturn.add(ring_lower);

더 복잡하고 사실적인 고리를 위해 셰이더 재질을 만들어 보자.

// vertex shader of Saturn ring
varying vec2 vUv;
varying vec3 vNormal;
varying vec3 vPosition;
varying mat3 vTbn;

void main() {
  vUv = uv;
  vNormal = normalize(mat3(modelMatrix) * normal);
  vPosition = mat3(modelMatrix) * position;
  
  gl_Position = projectionMatrix * modelViewMatrix * vec4(position, 1.0);

}

// fragment shader of Saturn ring
uniform sampler2D u_colorTexture;
uniform sampler2D u_alphaTexture;
uniform vec3 u_sunRelPosition; // the relative position of light source
uniform vec3 u_moonPosition; // the position of obstacle
uniform vec3 u_position; // the position of this object
uniform float u_sunRadius; // the radius of light source
uniform float u_moonRadius; // the radius of obstacle

varying vec2 vUv;
varying vec3 vNormal;
varying vec3 vPosition;

#define PI (3.141592)

float eclipse(float angleBtw, float angleLight, float angleObs) {
  float angleRatio2 = pow(angleObs / angleLight, 2.);
  float value;
  if (angleBtw > angleLight - angleObs && angleBtw < angleLight + angleObs) {
    if (angleBtw < angleObs - angleLight) {
      value = 0.;
    }else {
      float x = 0.5/angleBtw * (angleBtw*angleBtw + angleLight*angleLight - angleObs*angleObs);
      float ths = acos(x/angleLight);
      float thm = acos((angleBtw-x)/angleObs);
      value = 1./PI * (PI - ths + 0.5 * sin(2. * ths) - thm * angleRatio2 + 0.5 * angleRatio2 * sin(2. * thm));
    }
  } else if (angleBtw > angleLight + angleObs)
    value = 1.;
  else { // angleBtw < angleLight - angleObs
    value = 1. - angleRatio2;
  }

  return clamp(value, 0., 1.);
}

void main( void ) {
  vec3 sunDir = normalize(u_sunRelPosition);

  // 1. Day and night texture
  vec3 ringColor = texture2D( u_colorTexture, vUv ).rgb;
  float ringAlpha = texture2D( u_alphaTexture, vUv ).r;

  float cosAngleSunToNormal = dot(vNormal, sunDir); // Compute cosine sun to normal
  float mixAmountSurface = 0.7 / (1. + exp(-10. * cosAngleSunToNormal)) + 0.3; // Sharpen the edge beween the transition

  // // 2. Eclipse
  vec3 surfacePosition = u_position + vPosition;
  float distSurfaceToSun = length(u_sunRelPosition);
  float cosAngleBtwSunMoon = dot(sunDir, normalize(u_moonPosition - surfacePosition));
  float angleBtwSunMoon = acos(cosAngleBtwSunMoon);
  float distSurfaceToMoon = length(u_moonPosition - surfacePosition);
  mixAmountSurface *= eclipse(angleBtwSunMoon, asin(u_sunRadius / distSurfaceToSun), asin(u_moonRadius / distSurfaceToMoon));

  vec3 color = mix( vec3(0.), ringColor, mixAmountSurface ); // Select day or night texture based on mix.

  // 3. Specular map texture with reflection
  vec3 reflectVec = reflect(-sunDir, vNormal); // reflected vector of sunlight
  float specPower = dot(reflectVec, normalize(cameraPosition - surfacePosition)); // dot product between reflected light and camera vector
  color += mixAmountSurface * pow(specPower, 3.0) * 0.3;

  gl_FragColor = vec4( color, ringAlpha );
}

오브젝트 표면 렌더링의 원리가 동일하므로, 고리와 행성 셰이더의 함수도 동일하다는 것을 알 수 있다. 식(蝕)에 대한 자세한 설명은 이전 글¹을 참고하길 바란다. 디스크의 윗면과 아랫면을 구현하기 위해, 회전 각도가 다른 (0도와 180도) 두 개의 디스크를 만들었다는 점에 유의한다. 위 GLSL 셰이더에서 고리의 텍스처는 디스크에 적용되며, 고리 텍스처의 예시는 여기나 여기에서 다운로드할 수 있다. 위 셰이더의 결과로 아래를 얻는다.

그러나 텍스처 이미지와 디스크 오브젝트의 UV 좌표가 서로 다르기 때문에, 행성에서 멀어지는 방향으로 고리 텍스처 이미지가 제대로 적용되지 않는다. 아래 이미지에서 디스크 오브젝트의 UV 좌표가 직교 좌표계 [R, G, B] = [U, V, 0]을 기반으로 함을 볼 수 있다. 따라서 디스크 오브젝트의 UV 좌표를 직교 좌표계에서 극좌표계로 변환해야 한다.

아래는 직교 좌표를 극좌표로 변환한다. 즉, 디스크 바깥쪽 경계의 정점들은 1의 값을, 안쪽 경계의 정점들은 0의 값을 갖는다.

const pos = geometry_ring.attributes.position;
const mid_point = 0.5 * (1.3 + 2.2);
let v3 = new THREE.Vector3();
for (let i = 0; i < pos.count; ++i){
  v3.fromBufferAttribute(pos, i);
  geometry_ring.attributes.uv.setXY(i, v3.length() < mid_point ? 0 : 1, 0);
}

고리의 geometry에 새로운 값을 설정한 후, 디스크 오브젝트의 UV 좌표는 아래와 같이 보정된다.

그러면 토성의 고리는 다음과 같이 제대로 렌더링된다.

고리의 그림자

고리가 토성 표면에 드리우는 그림자를 렌더링하기 위해, 토성 본체의 material 객체에 uniform 변수를 추가한다.

const material_saturn = new THREE.ShaderMaterial({
  uniforms: {
    u_dayTexture: { value: new THREE.TextureLoader().load( './assets/saturnmap.jpg') },
    u_sunRelPosition: { value: new THREE.Vector3(0,0,0)},
    u_position: { value: new THREE.Vector3(0,0,0)},
    u_sunRadius: { value: 0.2},

    // added variables
    u_alphaTexture: { value: new THREE.TextureLoader().load( 'assets/saturnringpattern.png') },
    u_ringNormal: { value: new THREE.Vector3(0,1,0) },
    u_radius1: { value: 0.2 * 1.3 },
    u_radius2: { value: 0.2 * 2.2 },
  },
  vertexShader: vertex,
  fragmentShader: fragment,
})

u_alphaTexture는 고리의 투명도 맵 텍스처, u_ringNormal은 고리가 놓인 평면의 법선 벡터, u_radius1과 u_radius2는 고리 디스크의 내반경 및 외반경이다. 여기서 0.2는 토성의 반지름이고, 1.3과 2.2는 각각 반지름 대비 토성 고리 디스크의 내반경과 외반경이다.

u_ringNormal \(\mathbf{n}\)과 토성 본체의 위치 u_position \(\mathbf{p}_o\)에 대해, 고리의 중심이 본체의 위치와 같다는 가정을 이용하면, 고리의 평면 방정식은

\[\mathbf{n} \cdot (\mathbf{x} - \mathbf{p}_o) = 0\]

이다. 본체 표면의 한 점 \(\mathbf{p}_s\)에 대해, 태양의 위치 \(\mathbf{p}_{\rm sun}\)과 표면 점을 지나는 직선의 방정식은

\[\mathbf{x} = \mathbf{p}_s + (\mathbf{p}_{\rm sun} - \mathbf{p}_o) * t, ~t > 0\]

이다. 그러면 위 평면과 직선의 교점 \(x_p\)는

\[\mathbf{x}_p = \mathbf{p}_s + (\mathbf{p}_{\rm sun} - \mathbf{p}_o) * t_p\]

와

\[\mathbf{n} \cdot (\mathbf{p}_s + (\mathbf{p}_{\rm sun} - \mathbf{p}_o) * t_p - \mathbf{p}_o) = 0\]

를 만족한다. \(t_p\)를 풀면 다음을 얻는다.

\[t_p = - \frac{\mathbf{n}\cdot (\mathbf{p}_s - \mathbf{p}_o)}{\mathbf{n}\cdot (\mathbf{p}_{\rm sun} - \mathbf{p}_o)}.\]

\(t_p > 0\)인 경우 교점은

\[\mathbf{x}_p = \mathbf{p}_s - (\mathbf{p}_{\rm sun} - \mathbf{p}_o) * \frac{\mathbf{n}\cdot (\mathbf{p}_s - \mathbf{p}_o)}{\mathbf{n}\cdot (\mathbf{p}_{\rm sun} - \mathbf{p}_o)}\]

이다. 다음으로, 교점이 고리 위에 있는지 판단하기 위해, 본체 중심으로부터지 교점까지의 거리 \(d\)를 고리의 내반경 및 외반경과 비교한다.

\[d = ||(\mathbf{x}_p - \mathbf{p}_o)||_2\]

\(d\)가 내반경보다 크고 외반경보다 작으면, 고리가 드리우는 그림자의 강도 ringAlpha를 얻을 수 있다.

vec3 surfacePosition = u_position + vPosition;
float s = -dot(u_ringNormal, surfacePosition - u_position) / dot(u_ringNormal, sunDir);
if (s > 0.) {
  vec3 point = surfacePosition + sunDir * s;
  float alphaRatio = (sqrt(dot(point - u_position, point - u_position)) - u_radius1) / (u_radius2 - u_radius1);
  if (alphaRatio > 0. && alphaRatio < 1.) {
    float ringAlpha = texture2D(u_alphaTexture, vec2(alphaRatio, 1.0) ).r;
    mixAmountSurface *= (1. - 0.8 * ringAlpha * hemisphereLight); // apply rings shadow
  }
}

지금까지의 기법들을 적용하면, 결과는 다음과 같다.

고리가 그림자 (고급 버전)

이전 글에서 구체에 의한 그림자 (월식, 일식 등) 을 수학적으로 정밀히 계산했었고, 덕분에 구체 광원이 가까워지거나 커질수록 본체가 고리에 드리우는 그림자는 흐려지도록 정확히 렌더링된다.

먼 광원	가까운 광원

그러나 광원의 겉보기 크기가 달라져도, 고리가 본체 표면에 드리우는 그림자의 경계는 여전히 날카롭고 선명하게 유지된다. 이러한 이유는 광원을 부피가 있는 구체가 아니라 점으로 가정했기 때문이다. 고리 그림자를 더욱 사실적으로 렌더링하기 위해, 고리 평면의 법선 벡터 방향으로 광원을 여러 구간으로 분할하고 각 세부 구간에 대해 셰이딩 계수들을 계산한다. 아래 그림에서 이 방법의 원리를 확인할 수 있다.

분할 방향, \(\mathbf{v}_{\rm sun,perp}\), 은 태양 방향에 수직이며 고리 평면의 법선 벡터를 따라 놓인다. 수학적으로는 다음과 같이 정의할 수 있다.

\[\mathbf{v}_{\rm sun,perp} = \rm{arg}\max_{\mathbf{v}} (\mathbf{v} \cdot \mathbf{n}_{ring})\]

이고

\[\mathbf{v} \cdot \mathbf{v}_{\rm sun} = 0,\]

이며, 여기서 \(\mathbf{n}_{\rm ring}\)은 고리 평면의 법선 벡터이고 \(\mathbf{v}_{\rm sun}\)은 표면에서 태양의 원점으로 향하는 방향 벡터다. 그러면 \(i\)번째 세부 구간에 대해 \(\mathbf{v}_{\rm sunRay}=\mathbf{v}_{\rm sun}+(i/DIV * R_{\rm sun}) ~\mathbf{v}_{\rm sun,perp}\)를 얻는데, 여기서 \(i\)는 \(-DIV, \ldots, DIV\) 범위의 정수이고, \(DIV\)는 분할 파라미터, 즉 전체 세부 구간의 수는 \(2 * DIV + 1\)이다. 정확한 그림자 효과를 흉내 내려면 \(DIV\)를 큰 값으로 설정하면 된다. 계산 부하와의 트레이드 오프를 고려하여, 값은 7이면 적절하다.

세부 구간의 광원 면적은 서로 다르므로, 각 세부 구간에 의한 그림자도 다른 세기로 더해져야 한다. 여기서는 각 구간의 폭 \(\sqrt{1-(i/DIV)^2}\)을 가중치로 사용했다. 따라서 행성 표면의 한 점에서의 전체 그림자 계수는 고리 알파 값들의 가중합이다. 아래 코드는 프래그먼트 셰이더의 추가 부분이다.

#define DIV 7

#if DIV == 0
  #define INV_DIV 1. // prevent nan
#else
  #define INV_DIV (1./float(DIV))
#endif

float eclipseByRings(vec3 surfacePosition, vec3 sunRadiusPerp) {
  float totalRingAlpha = 0.;
  for (int i = -DIV; i <= DIV; ++i) {
    vec3 sunRay = u_sunRelPosition + INV_DIV * float(i) * sunRadiusPerp;
    float weight = sqrt(1. - pow(float(i)*INV_DIV, 2.));
    float s = -dot(u_ringNormal, surfacePosition - u_position) / dot(u_ringNormal, sunRay);
    if (s > 0.) {
      vec3 point = surfacePosition + sunRay * s;
      float alphaRatio = (sqrt(dot(point - u_position, point - u_position)) - u_radius1) / (u_radius2 - u_radius1);
      if (alphaRatio > 0. && alphaRatio < 1.) {
        float ringAlpha = texture2D( u_alphaTexture, vec2(alphaRatio, 1.0) ).r;
        totalRingAlpha += ringAlpha * weight;
      }
    }
  }
  totalRingAlpha /= (2. * float(DIV) + 1.);
  return totalRingAlpha;
}

따라서 이전 프래그먼트 셰이더의 main 함수는 다음과 같이 된다.

// 2. Shadow of rings
vec3 surfacePosition = u_position + vPosition;
vec3 sunRadiusPerp = u_ringNormal - dot(sunDir, u_ringNormal)/dot(sunDir, sunDir) * sunDir;
sunRadiusPerp = normalize(sunRadiusPerp) * u_sunRadius / length(u_sunRelPosition);
float ringAlpha = eclipseByRings(surfacePosition, sunRadiusPerp);
mixAmountSurface *= (1. - 0.8 * ringAlpha * hemisphereLight);
vec3 color = mix( vec3(0.), dayColor, mixAmountSurface ); // Select day or night texture based on mix.

결과적으로, 광원이 행성에 가까워질수록 고리가 드리우는 그림자가 흐려지는 것을 볼 수 있다.

먼 광원	가까운 광원

1. [Three.js] Create a Realistic Earth with Shaders ↩ ↩²