概要

今回はちょっと趣向を変えてコードリーディングをしてみたいと思います。

まずは以下の動画を見てください。

https://www.shadertoy.com/view/Ms2SD1

とてもきれいな海の映像ですよね。
これ、レイマーチングと呼ばれる方法でプログラムでリアルタイムに描いている映像なのです。
Shadertoyの投稿作品のページに飛べば、実際にレンダリングされている様子を見ることができます。

投稿者でなくても一時的なコードの変更をすることができるので（当然保存はできない）、パラメータをいじったりするとリアルタイムに内容が反映されるのが分かるかと思います。

今回はこの海のコードを読み解きつつ、内容を理解していくのが目的となります。
なお、この海のレンダリングにはいくつかの技術が使われており、それぞれについて個別に記事をまとめたのでそちらも合わせて見てもらうとより理解が深まるかなと思います。

コードを見てみる

さて、冒頭の作品のURLはShadertoyでも人気の作品です。
ただ解読するにはややノイズ（ノイズ関数ではなくｗ）があり、演出のために設けたパラメータなどもあります。
そこで、理解しやすさを優先して演出部分は削除しつつ、必要なコードのみにしたものを以下に投稿しました。

https://www.shadertoy.com/view/MttfW8

投稿したこちらのコードをベースに解読を行っていこうと思います。
これを実行したものは以下のようになります。

まずは全体感を見るために、コード全文を載せます。

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

/*
 * "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);
}

だいぶ短いコードですね。
（まぁ読み解く、という視点で見るとやや骨が折れそうな長さではありますが）

しかしながら、たったこれだけのコード量であの映像が描かれているというのは驚きです。

空の色

まずは簡単なところから攻めましょう。
空の色の取得についてです。コードは以下。

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);
}

やっていることはシンプルに、与えられたvec3型の値のうちy成分のみを使って色を決定しています。（おそらくeyeのeだと思う）
つまり、水平線近くは白く、上の方にいくにしたがって青くなっていくような計算になっているわけです。
（なのでよく見てみると、RGBどの成分も1.0 - e.yという式が入っている）

空の色はこれですべてです。簡単ですね。

海の色

続いて海の色です。

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;
}

空の色と比べるとだいぶ複雑です。
ただ、大まかに計算の内容を分類すると以下のようになります。

フレネル（fresnel ）の式を求める
屈折先の海の色を計算する
反射先の空の色を計算する
（フレネルの式の値を元に）それを合成する

フレネルの式（フレネルの式についてはこちら）は、水やガラスなどの「透過」する素材について、どのくらいの角度からだとどのくらいの反射（屈折）をするのか、を求める式です。

Wikipediaから画像を引用させてもらうと以下のイメージです。

なお、フレネルについては以下の記事が画像付きで説明してくれているのでとても分かりやすいと思います。

フレネルの式

少しだけ余談でフレネルの式について。
マルペケさんの記事から引用させてもらうと、

フレネルさんはこう考えました。「透過する平らな物に光が鉛直に近い角度で飛び込んだら、その光の大部分は透過物を貫くだろう。一方で平行に近い角度で飛び込むと、水面で反射してしまい光は透過物内にほとんど入らないだろう」：

この前提を元に、とある式を用いて反射率などを求めるのが「フレネルの式」です。

式は以下の形になります。

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

$R$は反射率、$\alpha$は入射光の入射角度、$\beta$は屈折光の法線との屈折角度です。

なお、コンピュータグラフィックスの世界ではSchlick氏によって導かれた近似式を用いることが多いです。
式は以下。

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

実際、今回のコードでは以下のように計算されており、似た形になっているのが分かります。

フレネルの式（コード）

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

閑話休題。

海は当然「水」なので、フレネルの式でレイの侵入角度から反射と屈折を計算します。

つまり、反射したレイに関してはその先にある空の色を、屈折したレイについてはその先にある海の色を利用するわけです。
そして最後にふたつの色をmix関数を使い、フレネルの式で求めた値を元に合成します。

それ以外にも、diffuse反射の計算やspecular反射の計算が行われていますが、こちらは3Dに関する基礎的なところなので今回は説明を割愛します。

ノイズ関数の意味

メインの処理（レイマーチング）に行く前に、少しだけノイズ関数の話を。
ノイズ関数は以下のように実装されています。

ノイズ関数

// 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;
}

これを実際に実行すると以下のようなノイズが得られます。

https://www.shadertoy.com/view/XlcfW4

やっていることは、入力された座標位置を元にノイズを計算することです。
まず、入力座標の整数部分を使ってグリッドを決定し、さらに少数部分を使ってその内部を補間しています。
コメントにも書いていますが、整数部分を利用してグリッドを決定するイメージは以下のようになります。

グリッド

// 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));

hash関数から、与えられた座標の上下左右の値を取得しているのが分かりますね。
そしてそれを、少数部分によって補間します。

hashと名のつくワケ

ちなみにhashと名前がついているのは、ここで得たいノイズ（ランダム値）は本当のランダムではなく、あくまで疑似乱数、つまり同じ引数なら同じ戻り値を得ることができる関数を求めているためです。
ハッシュ関数も、引数が同じなら同じ値を返してくれる関数ですよね。なので「hash」と名前がついているのだと思います。

閑話休題。
補間している部分のコードを見てみましょう。

補間

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;

まずa,bおよびc,dを横方向（u.x）に補間し、その結果をさらに縦方向（u.y）に補間しています。
最後に少しだけ加工していますが、これは-1.0～1.0に補正する常套手段ですね。

こうすることによって、上で示した図のようなノイズが得られる、というわけです。

なお、このノイズは「バリューノイズ」と呼ばれるシンプルなノイズ関数です。
こちらの記事がとても分かりやすく解説されているので、興味がある方は読んでみるとより理解が進むと思います。

The Book of Shaders: Noise

本丸「heightMapTracing」を読み解く

さて、いよいよ今回のコードの本丸、heightMapTracing関数を読んでいきましょう。
名前から分かるように、計算して求めたノイズから「Height Map」を生成し、そこからレイマーチングを行っているのだと思います。

Height Mapを利用したレイマーチングについては以前、記事に書いたのでHeight Mapを使うことでなぜレイマーチング可能なのかについてはそちらの記事を参照ください。

レイマーチングでHeight Map Distance Field

まずはheightMapTracing関数を見てみましょう。

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関数が、実際にはレイマーチングの処理をしているものと思われます。

続いてmap関数を見てみましょう。

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;
}

コード中でコメントした「フラクタルブラウン運動」についてはQiitaに記事を書いているのでそちらをご覧ください。

map_detailedについて

ちなみに、map_detailedという関数がありますが、こちらはイテレーションの回数が違うのみで実装はまったく同じものです。
イテレーションカウントにITER_FRAGMENTと記載されているように、おそらく処理負荷の問題から、ざっくり計算のmap関数と詳細な値を求めるmap_detailed関数を分けているのだと思います。

わざわざ実装を分けている理由は、CPUでの計算なら関数の引数にイテレーション回数を指定して統一化したいところですが、GLSLではそれが許されていない（for文のイテレーション回数はコンパイル時に決まっている必要がある）ためです。

さて、続いてsea_octave関数を見てみましょう。

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);
}

冒頭では入力のuv値を前述のnoise関数で若干ノイズを加算しています。
しかしながら、微細なノイズを追加するのが目的で、大局的に見たときの「海の波」は続く下の処理で実現しています。

一見、複雑そうな計算に見えますが、sin, cosを使っている計算部分をグラフに落とし込んでみるとなにをしようとしているのかが見えてきます。
こちらのサイトでグラフを可視化してみます。

まずは入力した数式。

数式自体はコードの計算をそのまま変換しただけです。

さぁ、これがどういうグラフになるか見てみましょう。
紫色の、やや尖ったグラフが最終的なグラフです。
（それ以外の色は、途中結果をグラフ化したものです。上の式の左側の色に対応しています）

グラフにしてみると分かりやすいですね。
そもそも海の波もsin, cosも「波」です。
なので、波を表現するのにsin波とcos波を合成してそれっぽい波を表現していた、というわけなんですね。

ちなみに最後のreturn pow(1.0 - pow(wv.x * wv.y, 0.65), choppy);部分ですが、ここは上段のwv.xの値をそのまま返しても似た形のアニメーションになります。
要は、波を複雑に合成したものを返せればそれでいいわけですね。

heightMapTracingをもう一度

heightMapTracing関数内で使われている他の関数を見てきました。
改めてheightMapTracingを見てみましょう。

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;
}

ここからは推測での話になります。
まずは以下の図を見てください。

図の「目」の部分が視点位置です。
そこからレイを飛ばし始めます。

関数の冒頭でhx = map(ori + dir * tx);が計算されています。
図を見てみると一目瞭然ですが、これははるかに遠い位置（tx = 1000.0）の位置のHeight Mapの計算となります。
そしてそれが0.0以上であれば即座にreturnしています。

おそらく、それだけ遠方のレイが0.0以上あるということは、それは海にレイが当たらないと見て即座にリターンしているのだと思います。
仮にもしその値が0.0以下であれば海にレイがヒットするとして計算を行います。

ここから先がまさにレイマーチングによる計算となります。
ちなみにこれは推測ですが、その後の計算でif文による分岐が行われていますが、これは最適化のためのものと思われます。
まだ自分も、しっかりと理解しきれいていないのですが、推測を交えつつ、いったいどんな計算がされているのかを見ていきましょう。

レイの現在位置と到達点、双方からレイの更新位置を決める

・・・などと見出しをつけましたが、あくまで推測です。（なので、もし間違っていたらご指摘ください；）

まず最初にhmを計算します。
これは開始時点でのレイの位置に応じたHeight Mapの距離です。
（そしてレイの位置の係数としてtmが用意されています。スタートは0で、徐々に値が変化していきます）

また前述のようにtx分レイを進めた値をhxとして保持し、ここから開始となります。
つまり、レイの始点と十分に遠い先のHeight Mapの高さを両端点として計算を開始するわけです。
（そしてこのtx先のレイがすでに水平面より上にある場合は、そもそも水面にレイはぶつからないとして早期リターンしているのは前述の通りです）

さて、ではこの求めた2つの値をどう使っていくのか見ていきましょう。

高さの比率からレイの位置を更新しサンプリングする

ループの冒頭でこんな式が出てきます。

float f = hm / (hm - hx);

このfの値を元に計算した値でレイの位置を更新しています。
その結果（hmid）を元にif文で分岐している箇所がありますね。
（なんでもいいけど、なぜhmidもtmidもidがつくんだろう・・・）

この式の意味は、hmとhxの比率に応じてどの位置のレイをサンプルするかを決めるための「重み」を計算することです。

仮に、hxが十分に大きい(※)とその値は限りなく0に近くなります。（分母が無限に大きくなり続けるならば結果はほぼ0になりますね）

※ ... なお、ここではhxを減じていますがこれは、hmはプラス、hxはマイナスであることを想定しているためと思われます。map関数の中を見てもらうと分かりますが、レイの高さ（ray.y）とHeight Mapの高さを減じているので、実質的に、Height Mapの高さからレイの高さへの符号付き距離を求めていることになるためです。そして、十分に遠いレイは水面はるか下にあるはずで、結果的にhxはマイナス距離になる、というわけです。

つまり、hx（tx）が十分に遠い場合にはできるだけ近くをサンプリングし、この比率が落ち着いてきたら徐々にtxに近い位置のサンプリングを行う、というような計算になっています。

なお、以下のコードを見てもらうと分かりますが、hmは正の値、hxは負の値になるよう更新されていきます。

分岐処理抜粋

float hmid = map(p);

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

おそらくここは「最適化」のための処理だと思われます。
レイを真面目に飛ばしてしまうととても負荷が高くなるので、ざっくりとレイを飛ばし、必要なところまでレイをすぐに移動させる、ということかなと。

実際、以前自分が実装したHeightMapによるレイマーチングでの実装は、イテレーション回数が100回を超えていました。
しかし、今回解読したコードはたかだか8しかループしていません。
にも関わらず、これだけきれいな映像が表現できているとうのは、この最適化による賜物なのかなと思っています。
この最適化手法に関してはぜひ習得したいところです。

終わりに

とにもかくにも、これでコードすべてを網羅しました。
最初にコードを見たときは「いったいこれはなにをしているんだ・・・？」と途方に暮れていましたが、ひとつずつ丁寧に読んでいくことでとても学びが多くありました。
レイマーチングのコードリーディングはとても楽しいですね。

次回は、今回学んだことを生かしてレイマーチングによる「雲」の表現をしてみたいと思っています。