Study Lab

수면 아래에서 본 바다

snell's window

Makonea
··17분
Step 1 Lab experiment
Loading lab experiment...

스넬의 창 — 아래에서 바라본 바다

Warning

해당 글은 AI로 대부분의 코드가 작성됨을 알림.

사람이 일부 경험하고 덧붙인 부분은 있으나, 수학적 표현 및 코드 대부분 AI가 작성함.

물속에서 위를 올려다보면 하늘 전체 180°가 반각 arcsin(na/nw)\arcsin(n_a/n_w)인 원뿔 안으로 압축된다. 이 국소 법선 기준 원뿔 밖에서는 공기에서 온 빛이 눈까지 도달할 수 없고, 계면은 전반사에 의한 거울이 된다. 물결치는 수면은 그 경계를 휘고 변형시키지만, 국소 임계각 법칙 자체를 바꾸지는 않는다.

물리 모델 — 측정 상수와 명시적인 모델 계수

항목

출처

nwn_w (해수, S=34.998S=34.998‰, 20 °C, 589.3 nm)

1.33938

Quan & Fry 1995; 임계각 48.30°, 창의 지름 96.60°

흡수 a\mathbf a (680/550/440 nm)

(0.4524, 0.0654, 0.0083) m⁻¹

NASA Ocean Optics Protocols에 표로 정리된 Pope & Fry 1997 값

전체 분자 산란 b\mathbf b (680/550/440 nm)

(0.0006, 0.0015, 0.0038) m⁻¹

NASA Ocean Optics Protocols에 표로 정리된 Morel 1974 값

하늘

Nishita 1993 단일 산란

표준 Rayleigh/Mie β

항목

nwn_w (해수, S=34.998S=34.998‰, 20 °C, 589.3 nm)

1.33938

출처

Quan & Fry 1995; 임계각 48.30°, 창의 지름 96.60°

항목

흡수 a\mathbf a (680/550/440 nm)

(0.4524, 0.0654, 0.0083) m⁻¹

출처

NASA Ocean Optics Protocols에 표로 정리된 Pope & Fry 1997 값

항목

전체 분자 산란 b\mathbf b (680/550/440 nm)

(0.0006, 0.0015, 0.0038) m⁻¹

출처

NASA Ocean Optics Protocols에 표로 정리된 Morel 1974 값

항목

하늘

Nishita 1993 단일 산란

출처

표준 Rayleigh/Mie β

이 구현은 의도적으로 맑은 물 하이브리드로 구성했다. 계면 굴절률에는 지정한 해수 조건을 사용하지만, 스펙트럼 감쇠에는 측정된 순수한 물의 흡수와 분자 산란을 사용한다. 특정 바다의 용존 물질, 식물성 플랑크톤, 부유 입자 산란은 포함하지 않는다.

AI 첫 구현은 청색 흡수에 aB=0.015m1a_B=0.015\,\mathrm{m^{-1}}를 사용했고, 스칼라 B_BACK=0.0025도 사용했다. 첫 번째 값은 440 nm보다 긴 청록색 파장에 훨씬 더 가깝고, 두 번째 값은 후방 산란 계수라고 이름 붙였지만 실제로는 전체 소멸값으로 사용됐다. 둘 다 수정했다. 이제 빔 소멸은 명시적인 관계 c=a+b\mathbf c=\mathbf a+\mathbf b를 따른다.

수학 감사: 정확한 검사와 범위가 제한된 근사

반구 전체가 정확히 원뿔 안에 들어간다. 하늘의 천정각 0–90°를 샘플링하고 각각을 물속으로 굴절시키면, 도달 가능한 최대 각도가 기계 정밀도 안에서 임계각과 일치한다.

밝은 가장자리는 그리지 않는다. 수평선 직전 하늘의 마지막 5°는 각반경의 0.51%를 차지한다. 360° 전체의 수평선 파노라마가 얇은 고리 안으로 몰려드는 셈이다. 원근 화면에서의 반지름은 투영 방식에 따라 달라지므로, 0.51%를 보편적인 픽셀 너비 비율로 제시해서는 안 된다.

프레넬. 수치로 계산한 R0R_0는 해석식[(nwna)/(nw+na)]2=0.02105[(n_w-n_a)/(n_w+n_a)]^2=0.02105와 일치하며, R1R\to1은 임계각에서 성립한다.

계면 복사휘도. 프레넬의 전력 비율은 R+T=1R+T=1을 만족하지만, 그것만으로 복사휘도 변환이 모두 설명되지는 았다.. 투과된 하늘빛은 L/n2L/n^2 불변 법칙도 따르므로, 공기 중 복사휘도가 물속으로 전달될 때 nw2n_w^2가 곱해진다.

카우스틱의 초점 거리가 맞는 지점. 진폭 a, 파수 k인 수면 잔물결은 굴절된 햇빛을 f = 1/(0.25339·a·k²)에 모은다. 여기서 사용한 너울 (λ = 0.8–1.8 m, a = 8–30 mm)은 7.7–10.8 m에서 초점을 만들며, 수면과 바닥 사이의 거리와 일치한다. 이 잔물결이 없으면 너울만으로는 100 m 너머에서 초점이 생겨 카우스틱이 전혀 나타나지 않았다.

카우스틱 추정기. 이득은 굴절된 태양광 착지 지도의 유한차분 야코비안을 사용한다. 평평한 수면에서 원시값은 1이고, 곡률이 빛을 모은다. 반환되는 표시값은 여기에 폴드 클램프와 경험적 정규화를 추가하므로 평평한 영역에서 반드시 1일 필요는 없다.FLOOR_NORMBEAM_NORM은 전역 에너지 보존의 증명이 아니라 경험적 보정 계수라는 것을 명시해둔다.

시도했지만 버린 지름길

처음 AI가 구현했을때, 카우스틱 야코비안을 통으로 구현되어있었다. 카우스틱 야코비안은 샘플 하나마다 수면 법선을 세 번 평가하므로 비싸다. 물결장의 헤시안과 수직 입사 계수 11/n1-1/n을 사용하는 해석적 버전은 더 저렴했지만, 기본 굴절 태양광이 비스듬히 들어오기 때문에 유한차분 순방향 지도와 크게 어긋났다.

현재 추정기는 평균 평면에서 역산한 초기 추측으로 시작하고, 하나의 국소 분기만 평가하며, 폴드 특이점을 클램프하고, 오프라인 정규화를 적용한다. 정확한 역 카우스틱 해법이 아니고, 폴드에서 가능한 모든 원상(preimage)을 합산하지도 않는다.

광선 추정기는 의도적으로 넓은 스텐실과 두 옥타브 수면을 사용한다. 그 스텐실보다 작은 옥타브는 명시적인 대역 제한 근사로 제외했다. 재현 가능한 샘플링 산출물이 없는 상태에서 이 선택을 검증된 물리적 등가물이라고 주장하지 않는다.

공식에서 실제 배포 코드까지

아래 방정식은 현재 실행되는 코드와 짝을 이룬다. 시각 모델 코드는 shaders/snells-window.frag가 소유하고, 브라우저 측 감사·로딩·상호작용·오디오 상태는 main.js가 소유한다.

런타임 파일은 각각 하나의 소유자를 가진다.

  • shaders/vertex.vert — 전체 화면 WebGL 버텍스 단계

  • shaders/snells-window.frag — 모든 광학·수면·카우스틱·물 GLSL

  • main.js — 명시적 셰이더 URL, 로딩, 컴파일, 컨트롤, 감사

  • underwatersound.mp3 — 선택 사항인 원본 오디오 녹음

인라인 셰이더 대체 경로는 없다. GLSL 파일이 없으면 오래된 JavaScript 복사본이 두 번째 진실의 원천이 되는 대신, 오류가 눈에 보이게 발생한다.

JavaScript
const SHADER_URLS = Object.freeze({
    vertex: new URL('./shaders/vertex.vert', import.meta.url),
    fragment: new URL('./shaders/snells-window.frag', import.meta.url)
});

1. 스넬 원뿔과 임계각

공기 굴절률을 nan_a, 물 굴절률을 nwn_w라고 하면,

nasinθa=nwsinθw.n_a\sin\theta_a=n_w\sin\theta_w.

공기 쪽 수평선은 θa=90\theta_a=90^\circ이므로, 물속에서 보이는 이미지는

θc=arcsin(nanw)=48.2979,2θc=96.5959.\theta_c=\arcsin\left(\frac{n_a}{n_w}\right) =48.2979^\circ, \qquad 2\theta_c=96.5959^\circ.

구현 — main.js, runAudit()

JavaScript
const thetaC = Math.asin(1 / N_W);
let maxSeen = 0;
for (let d = 0; d <= 90; d += 0.5) {
    maxSeen = Math.max(
        maxSeen,
        Math.asin(Math.sin(d / DEG) / N_W)
    );
}

셰이더는 역방향 카메라 광선, 즉 물 → 공기 방향을 사용한다. 따라서 GLSL의 refract(I,N,eta)에는 η=nw/na=NW\eta=n_w/n_a=N_W가 들어간다.

GLSL
vec3 rr0 = refract(rd, -n, N_W);

refract()는 제곱근 안의 값이 음수일 때 영벡터를 반환한다. 이는 국소 원뿔 밖에서 전반사가 일어나는 분기다.

2. 정확한 비편광 유전체 프레넬

물 쪽 입사각의 코사인을 cw=cosθwc_w=\cos\theta_w라고 하면,

sinθa=nw1cw2,ca=1sin2θa.\sin\theta_a=n_w\sqrt{1-c_w^2}, \qquad c_a=\sqrt{1-\sin^2\theta_a}.
Rs=(nwcwcanwcw+ca)2,Rp=(cwnwcacw+nwca)2,R=Rs+Rp2.R_s=\left(\frac{n_wc_w-c_a}{n_wc_w+c_a}\right)^2, \qquad R_p=\left(\frac{c_w-n_wc_a}{c_w+n_wc_a}\right)^2, \qquad R=\frac{R_s+R_p}{2}.

수직 입사에서는,

R0=(nwnanw+na)2=0.021046.R_0=\left(\frac{n_w-n_a}{n_w+n_a}\right)^2=0.021046.

구현 — shaders/snells-window.frag, fresnelWA()

GLSL
float fresnelWA(float cosW) {
    float s = sqrt(max(0.0, 1.0 - cosW * cosW)) * N_W;
    if (s >= 1.0) return 1.0;
    float ca = sqrt(1.0 - s * s);
    float rs = (N_W * cosW - ca) / (N_W * cosW + ca);
    float rp = (N_W * ca - cosW) / (N_W * ca + cosW);
    return clamp(0.5 * (rs * rs + rp * rp), 0.0, 1.0);
}

코드의 rp 진폭 부호는 위에 표시한 관례와 반대다. 그러나 전력 항에서는 제곱하므로 RpR_p 자체는 같다.

3. 프레넬 적용 범위와 n2n^2 복사휘도 법칙

손실이 없는 굴절 계면에서는,

Lana2=Lwnw2.\frac{L_a}{n_a^2}=\frac{L_w}{n_w^2}.

프레넬 투과를 적용하면,

Lw=(1R)(nwna)2La.L_w=(1-R)\left(\frac{n_w}{n_a}\right)^2L_a.

셰이더는 세 개의 수면 법선 샘플로 서브픽셀 범위를 근사한다. 투과율과 굴절 방향을 함께 평균한 다음, 가중된 방향에서 하늘을 한 번 평가한다.

구현 — shaders/snells-window.frag, main()

GLSL
float T0 = 1.0 - R0, T1 = 1.0 - R1, T2 = 1.0 - R2;
vec3 rr0 = refract(rd, -n, N_W);
vec3 rr1 = refract(rd, -nA, N_W);
vec3 rr2 = refract(rd, -nB, N_W);
vec3 transmittedDirection = rr0 * T0 + rr1 * T1 + rr2 * T2;
float transmission = (T0 + T1 + T2) / 3.0;
float R = 1.0 - transmission;

if (dot(transmittedDirection, transmittedDirection) > 1e-6) {
    above = skyColor(normalize(transmittedDirection), sd, 1.0)
          * (N_W * N_W);
}
col = above * (1.0 - R) + mirrored * R;

이것은 완전한 슈퍼샘플링이 아니라 커버리지 안티앨리어싱이다. 굴절된 하늘 광선 세 개를 하나의 가중 방향으로 합치며, 반사 분기는 여전히 큰 스케일의 법선 하나만 사용한다.

4. 물의 흡수, 산란, 빔 소멸

축약된 물기둥 모델은 흡수 a\mathbf a, 전체 산란 b\mathbf b, 빔 소멸 c\mathbf c를 나눈다.

c=a+b,T(s)=exp(cs).\mathbf c=\mathbf a+\mathbf b, \qquad \mathbf T(s)=\exp(-\mathbf c s).

8단계 단일 산란 근사는 다음과 같다.

L(s)=L0ecs+kbS(xk)ectkΔs.\mathbf L(s)=\mathbf L_0e^{-\mathbf cs} +\sum_k \mathbf b\,\mathbf S(\mathbf x_k)\, e^{-\mathbf c t_k}\Delta s.

구현 — shaders/snells-window.frag, 상수와 lookDown(), main()

GLSL
const vec3 A_W = vec3(0.4524, 0.0654, 0.0083);
const vec3 B_W = vec3(0.0006, 0.0015, 0.0038);
const vec3 C_W = A_W + B_W;

vec3 sunHere = SUN_I * sunT
             * exp(-C_W * (-p.y) / max(-dMean.y, 0.2));
inSc += B_W * (sunHere * 0.13 * gg + vec3(0.10))
       * exp(-C_W * tt) * ds;
col = col * exp(-C_W * min(pathLen, 60.0)) + inSc;

0.13은 측정된 체적 산란 함수가 아니라 유효 위상/광원 가중치다. 따라서 이것은 일관된 축약 단일 산란 폐쇄 모델이지, 스펙트럼 해양 복사전달방정식(RTE) 해법은 아니다.

5. 평균 중심 수면과 최초 교차 방정식

절차적 옥타브 함수는 양수 값만 반환한다. 평균 보정이 없으면 이전 수면의 평균은 대략 +0.4m+0.4\,\mathrm m였고, 평균 수면 아래 깊이라고 표시한 조절값이 그만큼 어긋났다. 수정된 장(field)은 다음과 같다.

hN(x,t)=j=0N1ajqj(x,t)μN+hchop(x,t),h_N(\mathbf x,t)= \sum_{j=0}^{N-1}a_j q_j(\mathbf x,t)-\mu_N +h_{chop}(\mathbf x,t),

여기서 μN\mu_N은 각 LOD에 대해 측정한 공간 평균이다. 카메라 광선 r(s)=o+sd{\mathbf r}(s)=\mathbf o+s\mathbf d가 수면에 닿는 최초 근은 다음 식의 해다.

F(s)=ry(s)h3(rx(s),rz(s),t)=0.F(s)=r_y(s)-h_3(r_x(s),r_z(s),t)=0.

구현 — shaders/snells-window.frag, seaHeight3()traceSurface()

GLSL
const float SEA_MEAN_3 = 0.4107;

float seaHeight3(vec2 xz, float t) {
    float freq = 0.22, amp = 0.32, choppy = 3.0;
    vec2 uv = vec2(xz.x * 0.75, xz.y);
    float h = 0.0;
    for (int i = 0; i < 3; i++) {
        h += (seaOct((uv + t) * freq, choppy)
            + seaOct((uv - t) * freq, choppy)) * amp;
        uv = OCTM * uv;
        freq *= 1.9;
        amp *= 0.22;
        choppy = mix(choppy, 1.0, 0.2);
    }
    return h - SEA_MEAN_3 + chop(xz, t) * uChop;
}

for (int i = 1; i <= 24; i++) {
    float ti = mix(tLo, tHi, float(i) / 24.0);
    vec3 p = ro + rd * ti;
    float h = p.y - seaHeight3(p.xz, t);
    if (h * hPrev <= 0.0) {
        float a = tPrev, b = ti, ha = hPrev;
        for (int j = 0; j < 8; j++) {
            float m = 0.5 * (a + b);
            vec3 pm = ro + rd * m;
            float hm = pm.y - seaHeight3(pm.xz, t);
            if (hm * ha <= 0.0) {
                b = m;
            } else {
                a = m;
                ha = hm;
            }
        }
        return 0.5 * (a + b);
    }
    tPrev = ti;
    hPrev = h;
}

이전의 조용한 대체 경로와 달리, 부호 변화가 없어 교차점을 찾지 못하면 -1.0을 반환한다. 더 이상 상단 탐색 평면에 임의의 교차점을 만들어내지 않는다.

6. 카우스틱 착지 지도와 유한차분 야코비안

수면의 광원 좌표를 s=(sx,sz)\mathbf s=(s_x,s_z), 굴절된 태양 방향을 d(s)\mathbf d(\mathbf s), 수면 높이를 h(s)h(\mathbf s), 바닥 높이를 yfy_f라고 하면, 착지 지도는 다음과 같다.

X(s)=s+dxz(s)h(s)yfdy(s).\mathbf X(\mathbf s)=\mathbf s+ \mathbf d_{xz}(\mathbf s) \frac{h(\mathbf s)-y_f}{-d_y(\mathbf s)}.

기하광학적 원시 국소 이득은 다음과 같이 근사한다.

G(s)1detDX(s).G(\mathbf s)\approx \frac{1}{|\det D\mathbf X(\mathbf s)|}.

실제로 화면에 표시되는 값은 다음과 같다.

Gdisplay=Nmin ⁣(1max(detDX,ϵJ),Gmax),G_{display}=N\,\min\!\left( \frac{1}{\max(|\det D\mathbf X|,\epsilon_J)},G_{max} \right),

여기서 NN, ϵJ\epsilon_J, GmaxG_{max}는 보정·정규화 항이지 광학 상수가 아니다. 따라서 원시 지도는 평평한 수면에서 1을 반환하지만, 보정된 반환값은 반드시 1일 필요가 없다.

구현 — shaders/snells-window.frag, causticGainFloor()

GLSL
float surfaceY = seaHeight4(sp, t);
vec3 d = refract(-sd, seaNormal4(sp, 0.05, t), 1.0 / N_W);
vec2 land = sp + d.xz * ((surfaceY - p.y) / (-d.y));

vec2 jx = (Lx - L0) / eps;
vec2 jz = (Lz - L0) / eps;
return min(
    1.0 / max(abs(jx.x * jz.y - jz.x * jx.y), 0.02),
    FLOOR_GMAX
) * FLOOR_NORM;

이전 코드는 이동 높이로 -p.y를 사용해 모든 광원 지점을 암묵적으로 y=0y=0에 두었다. 이제 계산한 파도 높이를 포함한다. 추정기는 여전히 평균 평면에서 역산한 초기 추측 하나를 사용하고, 모든 카우스틱 폴드의 원상을 찾거나 합산하지 않는다. FLOOR_GMAXFLOOR_NORM은 명시적인 정규화·보정 상수다.

7. 작은 기울기 근사에서의 초점 거리

수직 입사에 가깝고 h(x)=acos(kx)h(x)=a\cos(kx)인 경우, 작은 기울기 굴절 광선의 편향은 다음 초점 거리를 만든다.

f1(11/nw)ak2=10.25339ak2.f\approx\frac{1}{(1-1/n_w)ak^2} =\frac{1}{0.25339\,ak^2}.

구현 — shaders/snells-window.frag, chop()

GLSL
return 0.030 * sin(3.5 * (xz.x * 0.92 + xz.y * 0.39) + t * 1.7)
     + 0.014 * sin(5.5 * (xz.x * 0.31 - xz.y * 0.95) - t * 2.1)
     + 0.008 * sin(8.0 * (xz.x * 0.71 + xz.y * 0.70) + t * 2.7);

기본 uChop=1에서 이 모드들의 초점은 대략 10.7, 9.3, 7.8 m에 생긴다. 이 관계는 작은 기울기와 수직 입사에 가까운 조건에서의 추정이며, 런타임 카우스틱 해법에 직접 사용되지는 않는다.

AI가 함께하면서 생겼던 문제들

임계각 근처에서는 프레넬 반사율이 약 2° 안에서 0.2에서 1.0까지 급격히 변한다. 따라서 픽셀당 샘플 하나만 사용하면 이 실제 반짝임이 과소 샘플링되어 반점처럼 보인다.

이 문제는 측정하기 전까지 세 번 정도 잘 못나왔다. 처음 세 수정안은 픽셀 단위로 동일한 출력을 만들었다. 중간값을 디버그 이미지로 렌더링하고 나서야 바로 정리됐다. 수면 법선은 완전히 매끄러웠고, RR은 프랙털 경계에서 순수한 이진 0/1 값이었다. 문제는 법선이 아니었다. 그 지점에서 프레넬 항이 사실상 계단 함수처럼 보였던 것이 원인이었다.

이전에 RR만 평균내려던 시도도 실패했다. 그 스텐실은 픽셀 풋프린트에서 유도했는데, ts * pixAng / max(rd.y, 0.12) ≈ 3.6 cm였고 법선은 수 미터 스케일에서 변한다. 실제로 작동한 수정은 법선이 존재하는 스케일까지 스텐실을 넓힌 것이다.

GLSL
float nEps = clamp(0.35 + ts * 0.03, 0.35, 1.6);   // metres, was clamp(ts*0.012, 0.05, 0.9)

해저 카우스틱과 광선도 같은 종류의 오류를 별도로 겪었다. 픽셀 풋프린트 LOD (lodC, lodA)를 적용해 풋프린트가 커질수록 카우스틱 이득이 1로 희미해지도록 수정했다.

창의 가장자리에는 잔여 반점이 남아 있다. 제대로 고치려면 슈퍼샘플링이 필요하지만 아직 적용하지 않았다. 물기둥에서는 여전히 단일 산란만 사용하며 다중 산란은 포함하지 않는다.

현재의 3탭 계면 필터는 투과 커버리지와 굴절 방향을 함께 평균한다.

AI가 작성한 코드는 RR만 평균했다. 중심 탭이 전반사이고 이웃 탭이 투과하는 경우, 중심의 refract()는 영벡터를 반환하지만 평균 뒤에도 (1R)(1-R)은 0이 아니었다. 그 결과 계면 에너지가 조용히 버려졌다.(이부분은 내 개인적인 수정이다)

어디에도 화이트 밸런스는 적용하지 않는다. 적색은 실제로 수 미터 안에서 사라진다. 수정된 계수로 계산한 흡수만 고려한 1/e1/e 길이는 R/G/B에 대해 2.21 / 15.29 / 120.48 m이고, 직접 빔 소멸 길이 1/(a+b)1/(a+b)는 2.21 / 14.95 / 82.64 m다. 그래서 결과가 파랗게 보이며, 잠수부가 적색 필터를 사용하는 이유도 여기에 있다.

개인적인 메모: AI로 작업을 해보면서 느끼는 건데 구현의 레퍼런스가 되는 논문의 값을 추출해서 기록 해두면 나보다 훨씬 빠르게 작업을 한다.

GLSL을 이용한 이러한 장난은 어렸을때부터 프로그래머가 되면 꼭 구현하고 싶었다. 하지만 실제로 재능이 없었기때문에 대부분 레이마칭 기법이나 남들 쓴것의 색깔과 모양 변형정도 밖에 할 수 없었는데, 세상이 참 격변한 느낌이다.

내가 못했던 걸 하게 되니까 즐겁다.

참고문헌

Optical Aspects of Oceanography (1974), 1–24.