[コードリーディング vol.1] レイマーチングによる波表現を読み解く ref: https://qiita.com/edo_m18/items/a575606a60b21f0d2c57

0.0以下の分岐

if (hmid < 0.0)
{
    tx = tmid;
    hx = hmid;
}

file13.txt

float f = hm / (hm - hx);

file16.txt

float f = hm / (hm - hx);

file18.txt

tmid = mix(tm, tx, f);

file3.txt

R = \frac{\biggl(\frac{\sin(\alpha - \beta)}{\sin(\alpha + \beta)}\biggr)^2 + \biggl(\frac{\tan(\alpha - \beta)}{\tan(\alpha + \beta)}\biggr)^2}{2}

file4.txt

F_r(\theta) \approx  F_0 + (1 - F_0)(1 - \cos \theta)^5

GetSeaColor

vec3 getSeaColor(vec3 p, vec3 n, vec3 l, vec3 eye, vec3 dist)
{
    float fresnel = clamp(1.0 - dot(n, -eye), 0.0, 1.0);
    fresnel = pow(fresnel, 3.0) * 0.65;

    vec3 reflected = getSkyColor(reflect(eye, n));    
    vec3 refracted = SEA_BASE + diffuse(n, l, 80.0) * SEA_WATER_COLOR * 0.12; 

    vec3 color = mix(refracted, reflected, fresnel);

    float atten = max(1.0 - dot(dist, dist) * 0.001, 0.0);
    color += SEA_WATER_COLOR * (p.y - SEA_HEIGHT) * 0.18 * atten;

    color += vec3(specular(n, l, eye,60.0));

    return color;
}

GetSkyColor

// Get sky color by 'eye' position.
// So, The color changed smoothly by eye level.
vec3 getSkyColor(vec3 e)
{
    e.y = max(e.y, 0.0);
    float r = pow(1.0 - e.y, 2.0);
    float g = 1.0 - e.y;
    float b = 0.6 + (1.0 - e.y) * 0.4;
    return vec3(r, g, b);
}

heightMapTracing

// Get Height Map
float heightMapTracing(vec3 ori, vec3 dir, out vec3 p)
{
    float tm = 0.0;

    float tx = 1000.0;

    // Calculate 1000m far distance map.
    float hx = map(ori + dir * tx);

    // if hx over 0.0 is that ray on the sky. right?
    if(hx > 0.0)
    {
        p = vec3(0.0);
        return tx;   
    }

    float hm = map(ori + dir * tm);    
    float tmid = 0.0;

    // NUM_STEPS = 8
    for (int i = 0; i < NUM_STEPS; i++)
    {
        // Normalized by max distance (hx is max far position), is it correct?
        float f = hm / (hm - hx);

        tmid = mix(tm, tx, f);                   
        p = ori + dir * tmid;

        float hmid = map(p);

        if (hmid < 0.0)
        {
            tx = tmid;
            hx = hmid;
        }
        else
        {
            tm = tmid;
            hm = hmid;
        }
    }

    return tmid;
}

heightMapTracing関数

// Get Height Map
float heightMapTracing(vec3 ori, vec3 dir, out vec3 p)
{
    float tm = 0.0;

    float tx = 1000.0;

    // Calculate 1000m far distance map.
    float hx = map(ori + dir * tx);

    // if hx over 0.0 is that ray on the sky. right?
    if(hx > 0.0)
    {
        p = vec3(0.0);
        return tx;   
    }

    float hm = map(ori + dir * tm);    
    float tmid = 0.0;

    // NUM_STEPS = 8
    for (int i = 0; i < NUM_STEPS; i++)
    {
        // Normalized by max distance (hx is max far position), is it correct?
        float f = hm / (hm - hx);

        tmid = mix(tm, tx, f);                   
        p = ori + dir * tmid;

        float hmid = map(p);

        if (hmid < 0.0)
        {
            tx = tmid;
            hx = hmid;
        }
        else
        {
            tm = tmid;
            hm = hmid;
        }
    }

    return tmid;
}

map関数

// pはレイの位置
float map(vec3 p)
{
    // 分かりやすさのために、指定されている値をコメントしました
    float freq = SEA_FREQ; // => 0.16
    float amp = SEA_HEIGHT; // => 0.6
    float choppy = SEA_CHOPPY; // => 4.0

    // XZ平面による計算
    vec2 uv = p.xz;

    float d, h = 0.0;    

    // ITER_GEOMETRY = 3
    // ここで「フラクタルブラウン運動」によるノイズの計算を行っている
    for (int i = 0; i < ITER_GEOMETRY; i++)
    {
        // #define SEA_TIME (1.0 + iTime * SEA_SPEED)
        // SEA_SPEED = 0.8
        // 単純に、iTime（時間経過）によってUV値を微妙にずらしてアニメーションさせている
        // iTimeのみを指定してもほぼ同じアニメーションになる
        //
        // SEA_TIMEのプラス分とマイナス分を足して振幅を大きくしている・・・？
        d = sea_octave((uv + SEA_TIME) * freq, choppy);
        d += sea_octave((uv - SEA_TIME) * freq, choppy);

        h += d * amp;

        // octave_m = mat2(1.6, 1.2, -1.2, 1.6);
        // これは回転行列・・・？
        // uv値を回転させている風。これをなくすと平坦な海になる
        uv *= octave_m;

        // fbm関数として、振幅を0.2倍、周波数を2.0倍して次の計算を行う
        freq *= 2.0;
        amp *= 0.2;

        // choppyを翻訳すると「波瀾」という意味
        // これを小さくすると海が「おとなしく」なる
        choppy = mix(choppy, 1.0, 0.2);
    }

    // 最後に、求まった高さ`h`を、現在のレイの高さから引いたものを「波の高さ」としている
    return p.y - h;
}

sea_octave関数

float sea_octave(vec2 uv, float choppy)
{
    uv += noise(uv);        
    vec2 wv = 1.0 - abs(sin(uv));
    vec2 swv = abs(cos(uv));    
    wv = mix(wv, swv, wv);
    return pow(1.0 - pow(wv.x * wv.y, 0.65), choppy);
}

tmidの計算

tmid = mix(tm, tx, f);

グリッド

// Get each grid vertices.
// +---+---+
// | a | b |
// +---+---+
// | c | d |
// +---+---+
vec2 i = floor(p);
float a = hash(i + vec2(0.0,0.0));
float b = hash(i + vec2(1.0,0.0));
float c = hash(i + vec2(0.0,1.0));
float d = hash(i + vec2(1.0,1.0));

ノイズ関数

// Get random value.
float hash(vec2 p)
{
    float h = dot(p, vec2(127.1, 311.7));   
    return fract(sin(h) * 43758.5453123);
}

// Get Noise.
float noise(in vec2 p)
{
    vec2 i = floor(p);
    vec2 f = fract(p);

    // u = -2.0f^3 + 3.0f^2
    vec2 u = f * f * (3.0 - 2.0 * f);

    // Get each grid vertices.
    // +---+---+
    // | a | b |
    // +---+---+
    // | c | d |
    // +---+---+
    float a = hash(i + vec2(0.0,0.0));
    float b = hash(i + vec2(1.0,0.0));
    float c = hash(i + vec2(0.0,1.0));
    float d = hash(i + vec2(1.0,1.0));

    // Interpolate grid parameters with x and y.
    float result = mix(mix(a, b, u.x),
                        mix(c, d, u.x), u.y);
    return (2.0 * result) - 1.0;
}

フレネルの式（コード）

float fresnel = clamp(1.0 - dot(n, -eye), 0.0, 1.0);
fresnel = pow(fresnel, 3.0) * 0.65;

レイマーチングで海を描く

/*
 * "Seascape" by Alexander Alekseev aka TDM - 2014
 * License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
 */

const int NUM_STEPS = 8;
const float PI      = 3.141592;
const float EPSILON = 1e-3;

#define EPSILON_NRM (0.1 / iResolution.x)
#define SEA_TIME (1.0 + iTime * SEA_SPEED)

// sea
const int ITER_GEOMETRY = 3;
const int ITER_FRAGMENT = 5;
const float SEA_HEIGHT = 0.6;
const float SEA_CHOPPY = 4.0;
const float SEA_SPEED = 0.8;
const float SEA_FREQ = 0.16;
const vec3 SEA_BASE = vec3(0.1, 0.19, 0.22);
const vec3 SEA_WATER_COLOR = vec3(0.8,0.9,0.6);
const mat2 octave_m = mat2(1.6, 1.2, -1.2, 1.6);

// Get random value.
float hash(vec2 p)
{
    float h = dot(p, vec2(127.1, 311.7));   
    return fract(sin(h) * 43758.5453123);
}

// Get Noise.
float noise(in vec2 p)
{
    vec2 i = floor(p);
    vec2 f = fract(p);

    // u = -2.0f^3 + 3.0f^2
    vec2 u = f * f * (3.0 - 2.0 * f);

    // Get each grid vertices.
    // +---+---+
    // | a | b |
    // +---+---+
    // | c | d |
    // +---+---+
    float a = hash(i + vec2(0.0,0.0));
    float b = hash(i + vec2(1.0,0.0));
    float c = hash(i + vec2(0.0,1.0));
    float d = hash(i + vec2(1.0,1.0));

    // Interpolate grid parameters with x and y.
    float result = mix(mix(a, b, u.x),
                        mix(c, d, u.x), u.y);

    // Normalized to '-1 - 1'.
    return (2.0 * result) - 1.0;
}

// lighting
float diffuse(vec3 n, vec3 l, float p)
{
    return pow(dot(n, l) * 0.4 + 0.6, p);
}

float specular(vec3 n, vec3 l, vec3 e, float s)
{
    float nrm = (s + 8.0) / (PI * 8.0);
    return pow(max(dot(reflect(e, n), l), 0.0), s) * nrm;
}

// Get sky color by 'eye' position.
// So, The color changed smoothly by eye level.
vec3 getSkyColor(vec3 e)
{
    e.y = max(e.y, 0.0);
    float r = pow(1.0 - e.y, 2.0);
    float g = 1.0 - e.y;
    float b = 0.6 + (1.0 - e.y) * 0.4;
    return vec3(r, g, b);
}

// Get sea wave octave.
float sea_octave(vec2 uv, float choppy)
{
    uv += noise(uv);        
    vec2 wv = 1.0 - abs(sin(uv));
    vec2 swv = abs(cos(uv));    
    wv = mix(wv, swv, wv);
    return pow(1.0 - pow(wv.x * wv.y, 0.65), choppy);
}

// p is ray position.
float map(vec3 p)
{
    float freq = SEA_FREQ; // => 0.16
    float amp = SEA_HEIGHT; // => 0.6
    float choppy = SEA_CHOPPY; // => 4.0

    // XZ plane.
    vec2 uv = p.xz;

    float d, h = 0.0;    

    // ITER_GEOMETRY => 3
    for (int i = 0; i < ITER_GEOMETRY; i++)
    {       
        d = sea_octave((uv + SEA_TIME) * freq, choppy);
        d += sea_octave((uv - SEA_TIME) * freq, choppy);
        h += d * amp;
        uv *= octave_m;
        freq *= 2.0;
        amp *= 0.2;
        choppy = mix(choppy, 1.0, 0.2);
    }

    return p.y - h;
}

// p is ray position.
// This function calculate detail map with more iteration count.
float map_detailed(vec3 p)
{
    float freq = SEA_FREQ;
    float amp = SEA_HEIGHT;
    float choppy = SEA_CHOPPY;

    vec2 uv = p.xz;

    float d, h = 0.0;

    // ITER_FRAGMENT = 5
    for (int i = 0; i < ITER_FRAGMENT; i++)
    {        
        d = sea_octave((uv + SEA_TIME) * freq, choppy);
        d += sea_octave((uv - SEA_TIME) * freq, choppy);
        h += d * amp;        
        uv *= octave_m;
        freq *= 2.0;
        amp *= 0.2;
        choppy = mix(choppy, 1.0, 0.2);
    }

    return p.y - h;
}

// p = ray position
// n = surface normal
// l = light
// eye = eye
// dist = ray marching distance
vec3 getSeaColor(vec3 p, vec3 n, vec3 l, vec3 eye, vec3 dist)
{
    float fresnel = clamp(1.0 - dot(n, -eye), 0.0, 1.0);
    fresnel = pow(fresnel, 3.0) * 0.65;

    vec3 reflected = getSkyColor(reflect(eye, n));    
    vec3 refracted = SEA_BASE + diffuse(n, l, 80.0) * SEA_WATER_COLOR * 0.12; 

    vec3 color = mix(refracted, reflected, fresnel);

    float atten = max(1.0 - dot(dist, dist) * 0.001, 0.0);
    color += SEA_WATER_COLOR * (p.y - SEA_HEIGHT) * 0.18 * atten;

    color += vec3(specular(n, l, eye,60.0));

    return color;
}

// tracing
vec3 getNormal(vec3 p, float eps)
{
    vec3 n;
    n.y = map_detailed(p);    
    n.x = map_detailed(vec3(p.x + eps, p.y, p.z)) - n.y;
    n.z = map_detailed(vec3(p.x, p.y, p.z + eps)) - n.y;
    n.y = eps;
    return normalize(n);
}

// Get Height Map
float heightMapTracing(vec3 ori, vec3 dir, out vec3 p)
{
    float tm = 0.0;

    float tx = 1000.0;

    // Calculate 1000m far distance map.
    float hx = map(ori + dir * tx);

    // if hx over 0.0 is that ray on the sky. right?
    if(hx > 0.0)
    {
        p = vec3(0.0);
        return tx;   
    }

    float hm = map(ori + dir * tm);    
    float tmid = 0.0;

    // NUM_STEPS = 8
    for (int i = 0; i < NUM_STEPS; i++)
    {
        // Normalized by max distance (hx is max far position), is it correct?
        float f = hm / (hm - hx);

        tmid = mix(tm, tx, f);                   
        p = ori + dir * tmid;

        float hmid = map(p);

        if (hmid < 0.0)
        {
            tx = tmid;
            hx = hmid;
        }
        else
        {
            tm = tmid;
            hm = hmid;
        }
    }

    return tmid;
}

// main
void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = fragCoord.xy / iResolution.xy;
    uv = uv * 2.0 - 1.0;

    const vec3 light = normalize(vec3(0.0, 1.0, 0.8)); 

    // ray 
    vec3 ori = vec3(0.0, 3.5, 5.0);
    vec3 dir = normalize(vec3(uv.xy, -2.0));

    // tracing
    vec3 p;
    heightMapTracing(ori, dir, p);

    vec3 dist = p - ori;
    vec3 n = getNormal(p, dot(dist, dist) * EPSILON_NRM);

    // color
    vec3 sky = getSkyColor(dir);
    vec3 sea = getSeaColor(p, n, light, dir, dist);

    // This is coefficient for smooth blending sky and sea with 
    float t = pow(smoothstep(0.0, -0.05, dir.y), 0.3);
    vec3 color = mix(sky, sea, t);

    // post
    fragColor = vec4(color, 1.0);
}

分岐処理抜粋

float hmid = map(p);

// サンプリング位置の高さがマイナス距離の場合は`hx`, `tx`を更新する
if (hmid < 0.0)
{
    tx = tmid;
    hx = hmid;
}
// そうでない場合は、`hm`, `tm`を更新する
else
{
    tm = tmid;
    hm = hmid;
}

補間

vec2 f = fract(p);

// u = -2.0f^3 + 3.0f^2
vec2 u = f * f * (3.0 - 2.0 * f);

// Interpolate grid parameters with x and y.
float result = mix(mix(a, b, u.x),
                   mix(c, d, u.x), u.y);
return (2.0 * result) - 1.0;