Sign in Sign up
Create Art arrow_drop_down
settings

SDF for raymarching (距離関数のスキル)

2020.12.16
gaz

前置き

最近は距離関数の事をSDFと表現することが多く見られるようになりました。以後距離関数をSDFと呼んでいきます。SDFはsigned distance fanctionの略です。
何故、符号付にこだわるかについて軽く書きます。第一に負数が存在しないと法線計算で正しく法線が取れない。もう一点はブーリアン演算のSubtractionが正しく機能しない。Subtractionはmax(-d1,d2)とマイナス距離を必要とするからです。法線計算がちょっとおかしいは許容できてもブーリアン演算の方は許容できません。
SDFはlength()とdot()しか使わないで説明していきます。多少の例外はありますが、その度説明します。length()とdot()だけでSDFが出来るのは不思議でしょうがカラクリがあります。それがfold(折りたたみ)です。これを駆使してプリミティブを作っていきます。
今回説明するのはSDFだけなのでプリミティブが見れる最小限度のshaderを用意しました。以後このmap関数(shadertoyではここにSDFが格納されるが主流なのでそれに倣いました)を書き換えることによってサンプルshaderとしていきます。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
    return length(p)-1.;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

2Dから３Dにプリミティブを拡張

SDFはlength()とdot()しか使わないので、foldの説明が主になります。foldは、ほとんどの場合２D平面でしか使わないので２Dのshaderで良いのですが、それでは芸が無いので２Dを使い３Dプリミティブを作っていきます。作る方法も説明します。それと同時に便利なclamp()のトリックも紹介します。まずはclamp()のトリックから。

p.x-=clamp(p.x,-h,h);

これは説明するよりも見てもらった方が速いですね。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
    p.xy*=rot(time*.5);
    p.xz*=rot(time*.3);
    float h=sin(time)*.5+.5;
    p.x-=clamp(p.x,-h,h);
    return length(vec2(length(p.xy)-.3,p.z))-.1;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}
bvcs1vc3p9f7gigeevmg.gif

トーラスが引き延ばされます。これは面白くて

float map(vec3 p){
  float r=.1;
  vec3 v=vec3(1,2,3);
  p-=clamp(p,-v,v);
  return length(p)-r;
}

こうするとroundBoxのSDFになります。ただし、r=0.;にするとSDFでは無くなります。注意。

bvcs3ts3p9f7gigeevng.gif

単純にp.x-=clamp(p.x,-h,h);とかしてますがhを任意の位置でも出来るので便利に使えたりします。

つぎに2Dから３Dにプリミティブを拡張。これはトーラスのSDFに秘密があります。length(vec2(length(p.xy)-.3,p.z))-.1;これがトーラスのSDF。このvec2(length(p.xy)-.3,p.z))のｘ要素はp.xy平面の２DのSDFになっています。正確に言えば2次元を1次元に投影した座標？とでも言うべきですかね。これと使っていない軸p.zでvec2を作り、このvec2でSDFを作ればトーラスになるという仕組みです。違うモノでやってみます。例外として
https://www.iquilezles.org/www/articles/distfunctions2d/distfunctions2d.htm
ここより２DのSDFを借りてきます。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float sdHexagon( in vec2 p, in float r )
{
    const vec3 k = vec3(-0.866025404,0.5,0.577350269);
    p = abs(p);
    p -= 2.0*min(dot(k.xy,p),0.0)*k.xy;
    p -= vec2(clamp(p.x, -k.z*r, k.z*r), r);
    return length(p)*sign(p.y);
}

float map(vec3 p){
    p.xy*=rot(time*.5);
    p.xz*=rot(time*.3);
    float sdf2d=sdHexagon(p.xy,.8);
    return sdHexagon(vec2(sdf2d,p.z),.3);
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}
bvcs57s3p9f7gigeevo0.gif

これはsdHexagonを借りなくても後で説明するpmodで作れます。

もう一つは押し出しです。まず２DのSDFを作ります。これを単にraymarchingで使うと２Dの形状が断面の無限柱ができます。これをスライスすれば押し出しです。ブーリアン演算の積であるmax()を使います。２Dの形状をxy平面で使い、厚みの半分の値をhとした場合相手のSDFはabs(p.z)-hになります。この２つのSDFのmax()を取れば押し出しです。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
    p.xy*=rot(time*1.);
    p.xz*=rot(time*1.);
    float sdf2d=abs(length(p.xy)-1.)-.2;
    float d=abs(p.z)-.3;
    return max(sdf2d,d);
}


void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}
bvcs64c3p9f7gigeevog.gif

2Dで曲線を描くと歪が出ます。それを補正するスキルがあります。
https://www.iquilezles.org/www/articles/distance/distance.htm
これを使い押出しをしてみます。同時にbevel(角面を落とす)もやってみます。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float bevel(float x)
{
    return max(.15,abs(x));
}

float bevelMax(float a, float b)
{
    return (a+b+bevel(a-b))*.5;
}

float func(float x)
{
    return sin(x*2.)*.3;
}

float map(vec3 p){
    p.xy*=rot(time*1.);
    p.yz*=rot(time*1.);
    p.y -= func(p.x);
    float e = 1e-2;
    // 微分する。
    float g = (func(p.x+e)-func(p.x-e))/(2.*e);
    // 補正 https://www.iquilezles.org/www/articles/distance/distance.htm
    p.y *= 1./sqrt(1.+g*g);
    // clamp()を使い半直線を書く。
    p.x -= clamp(p.x, -5., 5.);
    float sdf2d=length(p.xy)-.6;
    // 押出し(bevel)
    // アーティファクト対策で重みを付ける。
    return bevelMax(sdf2d,abs(p.z)-.05)*.8; 
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution.xy)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-8);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=6./i;
}
bvcs7nk3p9f7gigeevp0.gif

この2つの２Dから３Dにするスキルを使うだけでもかなりのバリエーションが出来ます。

foldについて

foldとは何か。日本語だと「折りたたみ」と呼ばれてます。スキル的にはabs()の事です。これは使い方次第ではmod()と似ている為に混同されます。私も長い間、混同していました。混同する人も多いと思うので、一度混同を解いておこうと思います。実際挙動が似ているなら使って問題が無いと思いますが、最近になって解った事の理解の妨げになりそうです。それはマンデルボックスとスムースabs()での理解です。mod()とabs()は別物として扱わなければならなそうです。mod()がダメという事ではありません。mod()はRepetitionというスキルとして住み分けになると思います。トンネルをするならmod()は必要不可欠です。foldはabs()の系列です。そうは言ってもmod系とabs系の違いは説明しておきます。mod系は境界が分断される。abs系は境界が連続している。これは適材適所の使い方になると思います。mod系は無限に分裂できる。abs系は境界が連続してる為にスムージングが出来る。スムージングについては後で説明します。

abs系とは

bvcs9k43p9f7gigeevq0.png

abs()は絶対値の事です。これをグラフで見ればy軸で左右対称になっています。鏡に例えたら良いと思います。これだけを使っていたらプリミティブにバリエーションが生まれてきません。そこでどうするか。対称にする境界に角度を付けていくということをしていくわけです。任意角のfoldです。ｖは正規化必須。

p-=2.0*min(0.0,dot(p,v))*v;

これが良く出回っている任意角のfoldですがmin()を使わないでabs()を使って書き換えます。

float g=dot(p,v);
p-=(g-abs(g))*v;

同じ挙動になります。わざわざ書き換えるのにも意味があります。smoothAbs()を使いたいからです。smoothMin()もあるのでそれを使う方法もありますがsmoothAbs()の方が拡張性の可能性が高そうです。smoothAbs()については後で説明します。
任意角のfoldは非常に扱いが難しい。なので参考になるgifをshaderを作りました。小さい白丸が補助ベクトル。赤はp.xが正の数エリア。青はp.yが正の数のエリア。

bvcsalk3p9f7gigeevr0.gif

コードを置いておくのでダイレクトにベクトル入れるとかマウスで使いたいのであれば改造して使ってください。

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define TAU atan(1.)*8.
float D(float d){
  return smoothstep(0.03,0.0,d);
}
void main(){
  vec2 p=(gl_FragCoord.xy*2.-resolution)/resolution.y;
  p*=1.5;
  float a=time*TAU/6.;
  vec2 v=vec2(cos(a),sin(a));
  fragColor+=D(length(p-vec2(clamp(p.x,-1.,1.),0)));
  fragColor+=D(length(p-vec2(0,clamp(p.y,-1.,1.))));
  fragColor+=D(length(p-clamp(dot(p,v),.0,.5)*v));
  fragColor+=D(length(p-v*.5)-.05);
  float g=dot(p,v);
  p-=(g-abs(g))*v;
  float len=length(p);
  if(len<.5)
  {
    fragColor.xyz+=vec3(1,0,0)*step(0.,p.x);
  }else if(len<1.){
    fragColor.xyz+=vec3(0,0,1)*step(0.,p.y);
  }
}

これはabs系なので折れた所から左右対称になります。なかなか使い道はありません。これを何回か使い角度を分割する使い方はあります。これが鏡面にならずに使えるなら利用価値が出てきます。良いアイデアがありました。
https://twitter.com/7CIT/status/1235459606840610816
これで2度折り3度折りが可能になってきます。これは取り扱いが非常に難しいけど、とても魅力的です。

bvcsc9k3p9f7gigeevrg.jpg

この絵は他にもスキルを使ってますが、この凸凹感がだせます。まだ他の使い道も有ると思います。そのfoldは

    float g=dot(p,v);
    p=(p-(g-abs(g))*v)*vec2(sign(g),1);
bvcrbps3p9f7gigeeve0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float bevel(float x)
{
    return max(.05,abs(x));
}

float bevelMax(float a, float b)
{
    return (a+b+bevel(a-b))*.5;
}


vec2 fold(vec2 p, vec2 v)
{
  float g=dot(p,v);
  return (p-(g-abs(g))*v)*vec2(sign(g),1);
}

float map(vec3 p){
    p.xy*=rot(time*1.);
    p.yz*=rot(time*1.);
    p.xy=fold(p.xy,normalize(vec2(2,1)));
    p.x+=.4;
    p.xy=fold(p.xy,normalize(vec2(2,-1)));
    p.x-=clamp(p.x,-1.,1.);
    float sdf2d=length(p.xy)-.2;
    // 押出し(bevel)
    return bevelMax(sdf2d,abs(p.z)-.15);
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

これに関しては補助ベクトル、オフセット(移動)は使って試すしかありません。補助ベクトルについてはgifを参考に適当に辺りを付けて試すしかなさそうです。オフセットはx軸に沿ってSDFを作るならx軸で移動。ｙ軸に沿ってSDFを作るならｙ軸で移動。これがコツになります。

極座標での分割にfoldを使う。

極座標での分割はpmod()が主流になってます。その方がお手軽です。ですが境界が分断されるので境界にスムージングが出来ません。スムージングを意識した場合abs系が欲しくなってきます。なのでpmod()と同等なのをabs()で書いてみます。せっかくなのでpmod()も載せておきます。nが分割数になります。

// pmod
#define TAU (atan(1.)*8.)
float a=mod(atan(p.y, p.x),TAU/n)-.5 *TAU/n;
p=length(p)*vec2(sin(a),cos(a));

シンプルな仕組みです。これで注意が必要なのは、p=length(p)*vec2(cos(a),sin(a));とvec2が入れ替わっているのが出回っています。挙動が違ってきます。前者は放射状に伸びる軸がｙ軸。後者はｘ軸。同じつもりで作っていると絵が出ないことがあります。
次にabs()を使う極座標分割。

#define TAU (atan(1.)*8.)
float h=floor(log2(n));
float a =TAU*exp2(h)/n;
for(int i=0;i<int(h)+2;i++)
{
    vec2 v = vec2(-cos(a),sin(a));  
    float g=dot(p,v);
    p-=(g-abs(g))*v;
    a*=0.5;
}

少し複雑になってきます。わざわざこれを作る理由は、境界面のスムージングに応用できるからです。スムージングは後で解説します。スムージング無しならpmodの方がお手軽なので、pmodを使ったサンプルを載せます。

bvcseok3p9f7gigeevsg.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))
#define TAU atan(1.)*8.

vec2 pmod(vec2 p, float n)
{
  float a=mod(atan(p.y, p.x),TAU/n)-.5 *TAU/n;
  return length(p)*vec2(sin(a),cos(a));
}

float map(vec3 p){
    p.xy*=rot(time*1.);
    p.xz*=rot(time*1.);
    p.xy=pmod(p.xy,5.);
    p.y-=.6+sin(time)*.5;
    vec3 v=vec3(.1,.2,.2);
    p-=clamp(p,-v,v);
    return length(p)-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

直角fold

これは90度に折るfoldです。これは任意角foldよりお手軽だし色々と使い道があります。

if(p.x<p.y)p.xy=p.yx;

コードを読めば納得みたいな良いトリックです。これを使えばboxも作れるし、こんなのも出来ます。イメージとしては6面体の中に空間を折り込む感じです。

bvcsfis3p9f7gigeevtg.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
    p.xy*=rot(time*1.);
    p.xz*=rot(time*1.);
    p=abs(p);
    p-=.7;
    if(p.x<p.y)p.xy=p.yx;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    return length(p.xy)-.1;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

この中のfold部分

 p=abs(p)-.7;
 if(p.x<p.y)p.xy=p.yx;
 if(p.x<p.z)p.xz=p.zx;
 if(p.y<p.z)p.yz=p.zy;

はbox foldとか呼ばれてたりします。なので以後box foldと呼んでいきます。

smooth abs()

境界のスムージングについてです。foldの魅力は、これにつきます。mod系では出来ません。少し前までは任意角foldのmin()をスムースminにして使っていたのですがsmooth abs()の妙味に気が付いてから、こればかりです。なのでfold関係を全てabs()を使って書き換えました。そしてabs()をsabs()に置き換えます。

#define sabs(x) sqrt(x*x+1e-2)
bvcsgus3p9f7gigeevu0.png

シンプルな式です。これで境界がスムージングされます。1e-2を変えればスムージングも変わります。

bvcsibc3p9f7gigeevvg.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))
#define sabs(p)sqrt(p*p+1e-2)

// if(p.x<p.y)p.xy=p.yx;
// のスムージング　タイプ
vec2 sfold(vec2 p)
{
  vec2 v=normalize(vec2(1,-1));
  float g=dot(p,v);
  return p-(g-sabs(g))*v;
}

float map(vec3 p){
  p.xy*=rot(time*1.);
  p.xz*=rot(time*1.);
  p=abs(p)-vec3(.7,.7,1.2);
  p.xz=sfold(p.xz);
  p.yz=sfold(p.yz);
  p.xy=sfold(p.xy);
  return p.x;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

こんな感じで使います。今まであげたfoldの関数内のabs()をsabs()と差し替えるだけで済みます。
abs()ですがこれはmin() max()の代替えに使えます。

#define Min(a,b) (a+b-abs(a-b))*.5
#define Max(a,b) (a+b+abs(a-b))*.5

こうなります。やはりabs()をsabs()と差し替える事が出来ます。

#define sabs(p) sqrt(p*p+1e-2)
#define smin(a,b) (a+b-sabs(a-b))*.5
#define smax(a,b) (a+b+sabs(a-b))*.5

これまでは既にsmin() smax()があるので目新しい事ではありません。面白いのはsabs()はカスタムできる事です。名前が無いのもなんなんでeaseAbs()と呼ぶことにしました。スムージング部分にイージングしていきます。このイージングがカスタム可能なのです。foldというよりもブーリアン演算寄りのスキルです。

bvcsjls3p9f7gigef000.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))
#define TAU atan(1.)*8.

float ease(float x)
{
    return exp(-x*x*25.)*.3;
}

float smoothAbs(float x)
{
    return sqrt(x*x+1e-2);
}

float easeAbs(float x)
{
    return smoothAbs(x)+ease(x);
}

float bevel(float x)
{
    return max(.1,abs(x));
}

float easeMin(float a, float b){
    return (a+b-easeAbs(a-b))*.5;
}

float bevelMin(float a, float b){
    return (a+b-bevel(a-b))*.5;
}

float smoothMin(float a, float b){
    return (a+b-smoothAbs(a-b))*.5;
}
vec2 fold90(vec2 p)
{
  vec2 v=normalize(vec2(1,-1));
  float g=dot(p,v);
  return p-(g-smoothAbs(g))*v;
}

float de0(vec3 p)
{
    vec2 q=vec2(length(p.xy)-.8,abs(p.z)-.8);
    q=fold90(q);
    float de1=q.x;
    return q.x;
}

float de1(vec3 p)
{
    p.x-=clamp(p.x,0.,1.5);
    return length(p)-.3;
}

float map(vec3 p){
  p.xy*=rot(time*1.);
  p.yz*=rot(time*1.);
  float d=de0(p);
  d=easeMin(d,de1(p));
  p.xy*=rot(TAU/3.);
  d=bevelMin(d,de1(p));
  p.xy*=rot(TAU/3.);
  d=smoothMin(d,de1(p));
  return d;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

こちらのshadertoyも合わせて見てください。コメント欄にiq氏から'polyAbs()'を教わりました。これも役に立つと思います。
https://www.shadertoy.com/view/3sVBRG

極座標変換

トンネルを作るのに使ったりします。トンネルは人気があるんでSDFとは離れますが書いてみました。

vec2(atan(p.x,p.y)/PI*3., length(p.xy)-1.);

基本xy平面をこれで変換します。これは中心より1.0離れたところにx軸が円状に存在して、x軸の範囲が-3.0~3.0になる。y軸が放射線状に伸びてる。この位が扱いやすい空間だと思います。

bvcsoq43p9f7gigef02g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define PI acos(-1.)
#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float sabs(float x)
{
    return sqrt(x*x+1e-4);
}

void signedSFold(inout vec2 p, vec2 v)
{
    float g=dot(p,v);
    p=(p-(g-sabs(g))*v)*vec2(sign(g),1);
}

void sFold90(inout vec2 p)
{
    vec2 v=normalize(vec2(1,-1)); ;
    float g=dot(p,v);
    p-=(g-sabs(g))*v;
}

float box(vec3 p, vec3 s)
{
    p=abs(p)-s;
    sFold90(p.xz);
    sFold90(p.yz);
    sFold90(p.xy);
    return p.x;
}

float map(in vec3 p)
{
    p.xy*=rot(time*1.);
    p.xz*=rot(time*1.);
    p.xy=vec2(atan(p.x,p.y)/PI*3., length(p.xy)-1.);
    p.x=mod(p.x,1.)-.5;
    p.y=abs(p.y)-.1;
    p.x=abs(p.x);
    p.x-=.2;
    signedSFold(p.xy,normalize(vec2(2,1)));
    p.x+=.05;
    signedSFold(p.xy,normalize(vec2(2,-1)));
    p.z=abs(p.z)-.5;
    return box(p,vec3(.3,.05,.25));
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

IFS(Iterated function system)

これはフラクタルのの一種だけど仕組みはfoldと相性が良いので使われてます。forループの中でfoldしていき最後にSDFで出力する。けっこう適当に書いても、それなりの絵が出きる。派手めなのでVJ,ライブコーディングとかには向いていると思います。forループの中でスケーリングするモノをフラクタル呼びます。フラクタルは別枠で説明します。IFSでフラクタルまで手を出さなくても、かなり多彩なことが出来てしまいます。サンプルはループの中でfoldだけのイージーな奴です。SDFはベクトル方向normalize(vec3(1))のシリンダーだけ。なのにIFSのおかげで、この表現になる。

bvcsm943p9f7gigef010.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

float map(vec3 p)
{
  p.x-=4.;
  p.z+=time*3.;
  p=mod(p,8.)-4.;
  for(int j=0;j<3;j++)
  {
     p.xy=abs(p.xy)-.3;
     p.yz=abs(p.yz)-sin(time*2.)*.3+.1,
     p.xz=abs(p.xz)-.2;
  }
   return length(cross(p,vec3(.5)))-.1;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=10./i;
}

pmod()を使ったIFSをトンネルにした奴。

bvcsn3k3p9f7gigef01g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

#define TAU atan(1.)*8.

vec2 pmod(vec2 p, float n)
{
  float a=mod(atan(p.y, p.x),TAU/n)-.5 *TAU/n;
  return length(p)*vec2(sin(a),cos(a));
}

float map(vec3 p)
{
    p.z-=-time*2.;
    p.z=mod(p.z,2.)-1.0;
    for(int i=0;i<8;i++)
    {
        p.xy=pmod(p.xy,8.);
        p.y-=2.;
    }
    p.yz = pmod(p.yz, 8.);    
    return dot(abs(p),normalize(vec3(7,3,6)))-.7;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=10./i;
}

もう一つpmodのIFS。

bvcsni43p9f7gigef020.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

#define TAU atan(1.)*8.

vec2 pmod(vec2 p, float n)
{
  float a=mod(atan(p.y, p.x),TAU/n)-.5 *TAU/n;
  return length(p)*vec2(sin(a),cos(a));
}

float map(vec3 p)
{
    p.xy*=rot(time*.3);
    p.yz*=rot(time*.2);
    for(int i=0;i<4;i++)
    {
        p.xy = pmod(p.xy,10.);
        p.y-=2.;
        p.yz = pmod(p.yz, 12.);
        p.z-=10.;
    }
    return dot(abs(p),normalize(vec3(13,1,7)))-.7;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=10./i;
}

IFSの理解の為のサンプル

forループの中でfold、移動。forループの中でfold、移動、回転。forループの中でfold、移動、スケーリング。と段々難易度が上がっています。難易度と言っても絵を出す難易度でスキル的には変わりません。この段階のサンプルを書きます。サンプルなので最初のfoldを使い回します。SDFも同じモノを使います。

まず最初にIFSしないでfoldだけの奴。box foldを使い、そして座標を移動させてます。以後このパターンは固定です。

fold、移動のみ

bve308k3p9f7gigef89g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
  p.xy*=rot(time*.5);
  p.xz*=rot(time*.8);
  p=abs(p)-.3;
  if(p.x<p.y)p.xy=p.yx;
  if(p.x<p.z)p.xz=p.zx;
  if(p.y<p.z)p.yz=p.zy;
  p.xy-=.2;
  float h=.5;
  p.x-=clamp(p.x,-h,h);
  // torus SDF
  return length(vec2(length(p.xy)-.5,p.z))-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-4);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

fold、移動のIFS

bve30ek3p9f7gigef8a0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
  p.xy*=rot(time*.5);
  p.xz*=rot(time*.8);
  for(int i=0;i<3;i++){
    p=abs(p)-.3;
    if(p.x<p.y)p.xy=p.yx;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    p.xy-=.2;
  }
  float h=.5;
  p.x-=clamp(p.x,-h,h);
  // torus SDF
  return length(vec2(length(p.xy)-.5,p.z))-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-7);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

fold、移動、回転のIFS

bve30kc3p9f7gigef8ag.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
  p.xy*=rot(time*.5);
  p.xz*=rot(time*.8);
  for(int i=0;i<3;i++){
    p=abs(p)-.3;
    if(p.x<p.y)p.xy=p.yx;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    p.xy-=.2;
    p.xy*=rot(.5);
    p.yz*=rot(.5);
  }
  float h=.5;
  p.x-=clamp(p.x,-h,h);
  // torus SDF
  return length(vec2(length(p.xy)-.5,p.z))-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-7);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

fold、移動、スケーリングのIFS

このスケーリングをするのをフラクタルと呼びます。これはイタレーションの中で単純に倍率を2倍にするだけです。それでループを抜けた後に倍率を戻しSDFを使う。

bve30u43p9f7gigef8b0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
  p.xy*=rot(time*.5);
  p.xz*=rot(time*.8);
  float s=1.;
  for(int i=0;i<3;i++){
    p=abs(p)-.3;
    if(p.x<p.y)p.xy=p.yx;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    p.xy-=.2;
    p*=2.;
    s*=2.;
  }
  p/=s;
  float h=.5;
  p.x-=clamp(p.x,-h,h);
  // torus SDF
  return length(vec2(length(p.xy)-.5,p.z))-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

fold、移動、スケーリング、回転のIFS

bve2vrc3p9f7gigef890.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float map(vec3 p){
  p.xy*=rot(time*.5);
  p.xz*=rot(time*.8);
  float s=1.;
  for(int i=0;i<3;i++){
    p=abs(p)-.3;
    if(p.x<p.y)p.xy=p.yx;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    p.xy=abs(p.xy)-.2;
    p.xy*=rot(.3);
    p.yz*=rot(.3);
    p*=2.;
    s*=2.;
  }
  p/=s;
  float h=.5;
  p.x-=clamp(p.x,-h,h);
  // torus SDF
  return length(vec2(length(p.xy)-.5,p.z))-.05;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

こうやって見てみるとIFS自体の仕組みはシンプルです。問題は絵の出し方。パラメータを駆使したり、SDFを変えてみたり、foldを色々と組み合わせたり、順番を入れ替えてみたり可能性はたくさんあります。回転やスケールを使えば全く絵が出てこないゾーンが有ったり、一筋縄ではいかない所がありますが、ハマれば良い絵がとれます。このサンプルはループ数を低めに設定してます。回数を増やせば良い絵が採れることも有りますが、更に難易度が上がっていきます。

本格的なフラクタル

本格的なフラクタルも何種類かあるようなので全部を網羅できませんがIFSで出来るモノの中のアポロニアン、マンデルボックス系について解説します。
まず絶対的に理解しなけれはいけないモノから。

float opScale( in vec3 p, in float s, in sdf3d primitive )
{
    return primitive(p/s)*s;
}

https://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
より抜粋
スケールをかけた空間から得た距離はスケールで割る。上の関数とは逆の事を書いてますが同じです。実際はIFSの中では、そうやってます。GPUの負荷は乗算よりも除算の方が高いので除算を一回で済ませるようにしてると思います。これを忘れると絵は出てきません。dot()を使うSDFに正規化された法線必須なのも同じ理由と思います。空間のスケールを合わせる事だと思います。次にフラクタルのIFSの作法ですね。今から書くコードでは絵は出ませんが模擬コードで手順を書きます。

float map(vec3 p)
{
  p=abs(p); // 無くても良いケースもあるが、大体の場合必要。この部分に違うfoldで細工をするのも面白い。
  float s=1.; // スケールの初期値。アーティファクトが出る場合は数字を上げる。
  vec3 offset=p*1.; // マンデルボックスはこれをオフセットに使ってる。
  for(int i=0;i<8;i++) // ループ回数は15回くらいが標準的なようですが、回数を減らせるトリックもあります。
  {
    p=fold(p); // 使えるfoldは何種類かあって差し替えて使う。でたらめではダメで条件がある。
    float r=transform(p); // ここは負数を使う事を視野に入れる。なのでループを抜けた後`abs()`が必要になることがある。
    s*=r;
    p*=r;
    p+=offset; // ここは任意のベクトルの場合もある。
  }
  s=abs(s); // ｒが負数、イタレーターが奇数の場合は必要。
  // round boxのスキル。これはスカスカの空間を補う役目になる。無くても良い場合がある。
  //float a;
  //p-=clamp(p,-a,a); // p-=min(p,a);も使える場合があり。`clamp()`と違う絵になったりもする。
  return length(p)/s; // 最終出力は点のSDF。スケールでの除算は必須。他のSDFを使うのも有り。
}

これが基本形です。スケールの面倒な処理がなければIFSと同じ手順で単純なのだけど、マンデルボックスのトリックが混迷を招いたと思ってます。GPUが弱い時代だったから座標のvec3とスケールのfloatを一緒にしてvec4にしたのかもしれません。マンデルボックスを座標とスケールを分けて書くことも出来ますが、同じというわけにはなりませんでした。多少は違いました。だけど本物も偽物も無いんだからで単純を取りました。これで第一関門突破でフラクタルがシンプルに見えてきます。
次にfold部分。ここに使われているモノの特徴。

p=2.*clamp(p,-1.,1.)-p;

これはマンデルボックスに使われてるfold。3X3に分割されています。

bvcsqos3p9f7gigef03g.png

他にもフラクタルに使われてるfoldは何点かありますが、それらは

p=-1.+2.*fract(.5-.5*p);
p=mod(p-1.,2.)-1.;
// ..ect
bvcsrlc3p9f7gigef040.png

これと同じ特徴があります。foldする前の原点とfoldした後も中央に位置する空間の原点が一致する。これが条件になるようです。
Apollonianのフラクタルにはp=mod(p-1.,2.)-1.は使えたがマンデルボックスには使えませんでした：マンデルボックスには正負のあった並びが必要みたいです。

bvcss8c3p9f7gigef04g.png

ｘ領域のみの色付けです。青が正。赤が負。この並びのfoldしか受け付けないようです。なので、この並びを模した5X5のfoldを作りました。

p=-sign(p)*(abs(abs(abs(p)-2.)-1.)-1.);

これをマンデルボックスで代替えコードとして使ったら成功しました。foldに対するアプローチはこれで良さそうです。
foldに関しては、ここから先は実験だけです。
foldの正負を確認用shaderです。かなり重要な要素になっています。使った事のあるfoldを載せてあります。ｘ軸限定ですが、赤が正の領域。黒が負の領域。

bvcssvc3p9f7gigef050.png 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

void main(){
  vec2 p=(gl_FragCoord.xy*2.-resolution)/resolution.y;
  p*=10.0;
  fragColor+=step(abs(p.x),.1);
  fragColor+=step(abs(p.y),.1);
  //p=2.*clamp(p,-1.,1.)-p;
  //p=sign(p)*(1.-abs(abs(p)-1.));
  //p-=2.*clamp(round(p/2.),-1.,1.);
  //p = -1.0 + 2.0*fract(0.5 - 0.5*p);
  //p=1.-mod(p-1.,2.);
  //p=-sign(p)*(abs(abs(abs(p)-2.)-1.)-1.);
  //p=1.-abs(abs(p-2.)-1.);
  //p=1.-abs(abs(abs(p-5.)-2.)-2.);
  //p=abs(mod(p-1.,2.)-1.)-1.;
  p=1.-abs(abs(abs(p-7.)-4.)-2.);
  //p=1.-abs(p-1.);
  fragColor+=step(length(p-clamp(p,-.5,.5)),.0);
  fragColor.x+=step(0.,p.x);
}

他にも可能性があると思いますが把握してる限りだとこれくらいです。

次はスケールとSDFについて。

 float s=1.;
  ........
  //float h=1.;
  //p-=clamp(p,-h,h);
  return length(p)/s;

スケールは初期値を1.0に設定しています。アーティストが出た場合は、この数字を上げて対応します。通常はreturn length(p)/s*.5;とかしたくなりますが、sの数値を上げた方がスマートです。この数字を上げるとレイの進行が短くなりステップ数が増えます。なのでmarchinngループの回数を増やさなければならない事が出てきます。GPU負荷を考えれば、なるべく少ない数字にしたいところです。
return length(p)/s;の最後に半径を引き算してないのが気になると思いますが半径は非常に小さい値になります。衝突判定に使われている0.001と同等又はそれ以下なので無視しても良さそうです。基本は return length(p)/s;だけでループ数は15回くらいですが、round boxを使う方法があります。これを使うとループ回数を減らす事が出来ます。ループ回数を極端に落としてもround boxのサイズを上げて形状を変えるのも面白いです。それとスカスカの点だけの形状の時 に使うと有効だったりします。この辺りは色々と試してください。round boxのaですが、ループ回数を増やせば数字が大きくなっていきます。その辺りは慣れですね。単純にlength(p)/sの後に数字を入れるのも有りです。 フラクタルは点の集合で形があるという感じです。基本フラクタルは最後にreturn length(p)/s;で終わりますが、違うものを試すのも面白いです。

return length(p.xy)/s-.01;
return length(cross(p,normalize(vec3(0,.5,1))))/s-.01; // 任意ベクトルのシリンダーのSDF。これは例外になりますが、個人的にはお気に入りです。
return dot(p,normalize(vec3(1,2,3)))/s;

等々。
色付けのスキルについて。
SDFで使っているスケールを使う方法があります。スケールをlog(s)かlog2(s)とかしてパレットを使うという方法が有ります。それと座標を使うのも有りですね。そうはいえケースバイケースだったりします。簡易パレットはshadertoyのテンプレートに使われている奴です。

cos(vec3(0,2,4)+log(s))*.5+.5

こんな感じで使います。本格的に使いたいなら https://www.shadertoy.com/view/ll2GD3
transformとoffsetについて。
transformは1.0/dot(p,p)をclamp()等で拡張するのが主流。ここは試行錯誤の部分なので解説できません。解っている範囲で言えば、アポロニアン系は正の数を使うので拡大して行く。マンデルボックス系は負の数を使う為、イタレーションの中で交互に正負が入れ替わり拡大せずにそこにいる感じ。ちなみにマンハッタン距離にsabs()を使った奴とか、LP normを試した事がありますが、それも面白いです。
カメラワークについて。
foldが画面センターに来る為、ｚ軸に移動すると直ぐにカメラがプリミティブの中に潜ってしまいします。なのでfoldの周期が半分の所に置くのが良いと思います。
box foldした空間の中でフラクタルをするのもマンデルボックスのように6面体みたいになっておもしろいです。マンデルボックスを入れるのも有りです。そこからmod()で増やすのもありだし、foldのスキルで色々と組み合わせられます。
アポロニアン系とマンデルボックス系は同じようにdot(p,p)を使う。違いは座標にかけるスケールが正の数か負の数かの違い。オフセットの取り方がちょっと違うかんじです。この２つのフラクタルしか理解出来てませんが、これだけで充分な量のバリエーションがあります。これに絡めてfoldを使ったり、modを使ったり極座標変換をしてみたり、dot(p,p)の代替えを入れてみたり、使われてるfoldを違うモノに差し替えたり、forループの中で番手によってパラメータを変えたり、出力のSDFを変えたり、ループの中で回転を入れたり・・・結局アイデアだなって感じました。このアイデアを考えるのは楽しいです。これだけヒントがあればフラクタルも難しくは無いと思います。サンプルは多めに用意しました。サンプルの中のほとんどの数値は弄れるので試してみてください。ただfoldの所は正しい理解が必要なので慎重にお願いします。それとfoldを差し替えてみるのも面白いと思います。基本はトライアンドエラーでしかないです。で、感覚でコントール。慣れればライブコーディングもイケるんじゃないかな。なんせコードが短いからね。

bvcsu743p9f7gigef060.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    // トンネル　パターン
    p.z-=-time*3.;
    p.xy=abs(p.xy)-2.;
    if(p.x<p.y)p.xy=p.yx;
    p.z=mod(p.z,4.)-2.;

    p.x-=3.2;
    p=abs(p);
    float s=2.;
    vec3 offset= p*1.5;
    for (float i=0.; i<5.; i++){
        p=1.-abs(p-1.);
        float r=-7.5*clamp(.38*max(1.2/dot(p,p),1.),0.,1.);
        s*=r;
        p*=r;
        p+=offset;
    }
    s=abs(s);
    scale=s;
    float a=100.;
    p-=clamp(p,-a,a);
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(1,9,3)+log2(scale))*.5+.5),.5)*10./i;
}
bvcsup43p9f7gigef070.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float Scale;
float map(vec3 p){
    p.xy*=rot(time*.5);
    p.xz*=rot(time*.5);
    p=abs(p);
    float s=3.;
    vec3  offset = p*.5;
    for (float i=0.; i<5.; i++){
        p=1.-abs(p-1.);
        float r=-3.*clamp(.57*max(3./dot(p,p),.9),0.,1.);
        s*=r;
        p*=r;
        p+=offset;
    }
    s=abs(s);
    Scale=s;
    return length(cross(p,normalize(vec3(1))))/s-.008;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-7);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=(cos(vec3(1,2,3)+log(Scale)*2.)*.5+.5)*15./i;
}
bvcsv843p9f7gigef07g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float Scale;
float map(vec3 p){
    p.yz*=rot(time*0.1+2.);
    p.xz*=rot(time*0.2+2.);
    p=abs(p)-2.;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    if(p.x<p.y)p.xy=p.yx;
    float s=2.5;
    vec3 off=p*2.8;
    for (float i=0.;i<6.;i++)
    {
        p=-sign(p)*(abs(abs(abs(p)-2.)-1.)-1.);
        float r=-11.*clamp(.8*max(2.5/dot(p,p),.2),.3,.6);
        s*=r;
        p*=r;
        p+=off;
        p.yz*=rot(2.1);
    }
    s=abs(s);
    Scale=s;
    float a=30.;
    p-=clamp(p,-a,a);
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-11);
    float d=1.,i;
    for(;++i<120.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=28.*(cos(vec3(9,5,12)+log(Scale)*2.)*.5+.5)/i;
}
bvcsvkc3p9f7gigef080.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    p.yz*=rot(time*.1);
    p.xz*=rot(time*.2);
    p=abs(p)-3.;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    if(p.x<p.y)p.xy=p.yx;
    float s=3.;
    vec3  offset=p*1.2;
    for (float i=0.;i<8.;i++){
        p=1.-abs(p-1.);
        float r=-6.5*clamp(.41*max(1.1/dot(p,p),.8),.0,1.8);
        s*=r;
        p*=r;
        p+=offset;
        p.yz*=rot(-1.2);
    }
    s=abs(s);
    scale=s;
    float a=20.;
    p-=clamp(p,-a,a);
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-10);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=8.*(cos(vec3(9,5,12)+log(scale)*3.)*.5+.5)*5./i;
}
bvct0343p9f7gigef08g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    p.yz*=rot(time*.1);
    p.xz*=rot(time*.2);
    // これは初期値が高い。こういう所にも形が潜んでいることがあるみたい。
    float s=12.;
    p=abs(p);
    vec3 offset=p*3.;
    for (float i=0.; i<5.; i++){
        p=1.-abs(p-1.);
        float r=-5.5*clamp(.3*max(2.5/dot(p,p),.8),0.,1.5);
        p*=r;
        p+=offset;
        s*=r;
    }
    s=abs(s);
    // box fold
    p=abs(p)-3.;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    if(p.x<p.y)p.xy=p.yx;
    scale=s;
    float a=3.;
    p-=clamp(p,-a,a);
    return length(p.xz)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-4);
    float d=1.,i;
    for(;++i<120.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(3,2,5)+log2(scale)*7.)*.5+.5),.5)*20./i;
}
bvct0lc3p9f7gigef090.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    p.yz*=rot(time*.3);
    p.xz*=rot(time*.2);
    float s=3.;
    p=abs(p);
    for (float i=0.; i<9.; i++){
        p-=clamp(p,-1.,1.)*2.;
        float r=6.62*clamp(.12/min(dot(p,p),1.),0.,1.);
        s*=r;
        p*=r;
        p+=1.5;
    }
    s=abs(s);
    scale=s;
    float a=.8;
    p-=clamp(p,-a,a);
    return length(p)/s;
}


void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-25);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(7,15,3)+log2(scale))*.5+.5),.8)*30./i;
}
bvct12c3p9f7gigef09g.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p)
{
    // Apollonianの場合はセンターを半周期ずらす。
    p+=vec3(1,1,time);
    float s=3.;
    for(int i = 0; i < 4; i++) {
        p=mod(p-1.,2.)-1.;
        float r=1.2/dot(p,p);
        p*=r;
        s*=r;
    }
    scale=s;
// ここを 0 にすれば Apollonianだってわかります。box foldを後ろに置いた時の効果。
#if 1
    p = abs(p)-0.8;
    if (p.x < p.z) p.xz = p.zx;
    if (p.y < p.z) p.yz = p.zy;
    if (p.x < p.y) p.xy = p.yx;
#endif
    // 任意ベクトルのシリンダーの距離関数。これはお気に入りなので例外。でもlength()は使ってる。
    return length(cross(p,normalize(vec3(0,1,1))))/s-.001;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-6);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(7,15,3)+log2(scale))*.5+.5),.8)*30./i;
}
bvct1d43p9f7gigef0a0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    p.yz*=rot(time*.1);
    p.xz*=rot(time*.2);
    // box foldを前に置くと6面体の辺に沿ってフラクタルが配置される。
    p=abs(p)-1.2;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    if(p.x<p.y)p.xy=p.yx;

    float s=1.;
    for(int i=0;i<6;i++)
    {
      p=abs(p);
      float r=2./clamp(dot(p,p),.1,1.);
      s*=r;
      p*=r;
      p-=vec3(.6,.6,3.5);
    }
    float a=1.5;
    p-=clamp(p,-a,a);
    scale=s;
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(1,2,3)+log2(scale))*.5+.5),.5)*10./i;
}
bvct1mk3p9f7gigef0b0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
#define sabs(p)sqrt((p)*(p)+1e-1)
#define sabs2(p)sqrt((p)*(p)+1e-3)

float map(vec3 p){
    p.xy*=rot(time*.2);
    p.yz*=rot(time*.1);
    float s=2.;
    p=abs(p);
    for (int i=0; i<4; i++) 
    {
        p=1.-sabs2(p-1.);
        // マンハッタン距離にsabs()で代用。
        float r=-9.*clamp(max(.2/pow(min(min(sabs(p.x),sabs(p.y)),sabs(p.z)),.5), .1), 0., .5);
        s*=r;
        p*=r;
        p+=1.;
    }
    s=abs(s);
    scale=s;
    float a=2.;
    p-=clamp(p,-a,a);
    return length(p)/s-.01;
}


void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(5,2,3)+log2(scale)*3.)*.5+.5),.8)*25./i;
}
bvct21s3p9f7gigef0bg.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float lpNorm(vec3 p, float n)
{
    p = pow(abs(p), vec3(n));
    return pow(p.x+p.y+p.z, 1.0/n);
}

float scale;
float map(vec3 p){
    p.xy *= rot(time*.2);
    p.yz *= rot(time*.1);
    vec3 offset=p*.5;
    float s=2.;
    for (int i=0; i<5; i++)
    {
        p=clamp(p,-1.,1.)*2.-p;
        // LP normで代用。
        float r=-10.*clamp(max(.3/pow(lpNorm(p,5.),2.),.3),.0,.6);
        s*=r;
        p*=r;
        p+=offset;
    }
    s=abs(s);
    scale=s;
    float a=10.;
    p-=clamp(p,-a,a);
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-5);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(7,8,3)+log(scale))*.5+.5),.3)*10./i;
}
bvct28k3p9f7gigef0c0.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

float scale;
float map(vec3 p){
    p.xy *= rot(time*.5);
    p.yz *= rot(time*.3);
    // box fold
    p=abs(p)-15.;
    if(p.x<p.z)p.xz=p.zx;
    if(p.y<p.z)p.yz=p.zy;
    if(p.x<p.y)p.xy=p.yx;
    float s=2.;
    for (int i=0; i<8; i++)
    {
        p=-sign(p)*(abs(abs(abs(p)-2.)-1.)-1.);
        float r=-1.55/max(.41,dot(p,p));
        s*=r;
        p*=r;
        p-=.5;
    }
    s=abs(s);
    scale=s;
    return dot(p,normalize(vec3(1,2,3)))/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-200);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(15,28,3)+log(scale))*.5+.5),.5)*18./i;
}
bvct2gs3p9f7gigef0cg.gif 
#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

#define sabs(x)sqrt((x)*(x)+.2)
void sFold90(inout vec2 p)
{
    vec2 v=normalize(vec2(1,-1));
    float g=dot(p,v);
    p-=(g-sqrt(g*g+1e-1))*v;
}

float scale;
float map(vec3 p){
    p.yz*=rot(time*0.05);
    p.xz*=rot(time*0.1);  
    p=abs(p)-1.8;
    // ここの並びを変えると組み合わせが変わり絵も変わる。
    sFold90(p.zy);
    sFold90(p.xy);
    sFold90(p.zx);
    float s=2.;
    vec3  offset=p*.5;
    for(int i=0;i<8;i++){
        p=1.-abs(p-1.);
        float r=-1.3*max(1.5/dot(p,p),1.5);
        s*=r;
        p*=r;
        p+=offset;
        p.zx*=rot(-1.2);
    }
    s=abs(s);
    scale=s;
    float a=8.5;
    p-=clamp(p,-a,a);
    return length(p)/s;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-10);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor.xyz+=mix(vec3(1),(cos(vec3(8,18,22)+log(scale))*.5+.5),.5)*18./i;
}

アーティファクトについて

ちょっと無理目なSDFをつかうとアーティファクトが出ます。これはレイのオーバーシュート(行き過ぎ)が原因です。この対策ですがSDFに重みをかけます。sdf*0.5とかレイの進む距離に重みをかけて使います。ただ、これを使うとraymarchingのステップ数が増えってGPU負荷が上がります。この辺りは匙加減が必要です。これは特殊ケースですがmod()で分裂させた空間のプリミティブに別々の変化をつけるとアーティファクトの原因になります。この場合はSDFに重みを付けても解決しません。ちゃんと境界までの距離を計算してそれとmin()を取る方法もありますが、簡易的にするトリックがあります。重さで制限をかけるのではなくてレイの進む距離に制限をかける方法です。SDFの数値をd、制限する距離をaとしたら

min(a,d);

とmin()を使うだけです。言われれば納得ですが、これが意外と盲点になっているようです。

あとがき

これだけ理解できればSDFに関しては多彩な表現ができると思います。shaderは他にも色付けとかエフェクト、カメラワーク、イージング等々、色々なジャンルがあります。それらを組み合わせていけば面白いモノが出来ると思います。サンプルの方は自由に改造してかまいません。これを使って投稿するならgazを引用したよって書いておいてください。

追記

直角foldif(p.x<p.y)p.xy=p.yx; のスムージング　タイプで

p=(p.x+p.y+vec2(1,-1)*sabs(p.x-p.y))*.5;

これも使えます。

if(p.x<p.y)p=p.yx;

と

p=vec2(max(p.x,p.y),min(p.x,p.y));

は同等。

#define sabs(p) sqrt((p)*(p)+1e-2)
#define smin(a,b) (a+b-sabs(a-b))*.5
#define smax(a,b) (a+b+sabs(a-b))*.5
p=vec2(smax(p.x,p.y),smin(p.x,p.y));

が 

p=(p.x+p.y+vec2(1,-1)*sabs(p.x-p.y))*.5;

となります。

追記２

今までのスキルを使っての応用で螺旋。トーラスのSDFからの発展。この場合、z方向にmod()で分裂させて、atan()を使いz方向に座標をスライドさせる。この時atan()から出力される値をPIで除算して正規化しておくと感覚で使えるようになる。
bvq7orc3p9f30ks56pgg.gif 

#version 300 es
precision highp float;
uniform vec2 resolution;
uniform float time;
out vec4 fragColor;

#define rot(a) mat2(cos(a),sin(a),-sin(a),cos(a))

#define PI (atan(1.)*4.)

float map(vec3 p){
    p.xy*=rot(time*.5);
    p.xz*=rot(time*.3);
    float c=.2;
    p.z+=atan(p.y,p.x)/PI*c;
    p.z=mod(p.z,c*2.)-c;
    return length(vec2(length(p.xy)-.4,p.z))-.1;
}

void main(){
    vec2 uv=(gl_FragCoord.xy-.5*resolution)/resolution.y;
    vec3 rd=normalize(vec3(uv,1));
    vec3 p=vec3(0,0,-3);
    float d=1.,i;
    for(;++i<99.&&d>.001;)p+=rd*(d=map(p));
    if(d<.001)fragColor+=3./i;
}

追記３

SDFとしてdot()を使う事に対する説明が欠損している事に気付いた。
例えになるがlength()を透視投影とすればdot()は平行投影。
SDFにおいて距離は面、線、点からの物しかない。length()は点、線からの距離、dot()は面から距離になる。
dot()に使う引数は座標と面の法線(正規化必須)。
ここでdot系列のSDFについて書いてみる。
p.yこれはdot(p,vec3(0,1,0))こういう事である。

https://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
ここにある8面体のSDF 

float sdOctahedron( vec3 p, float s)
{
  p = abs(p);
  return (p.x+p.y+p.z-s)*0.57735027;
}

これの(p.x+p.y+p.z-s)*0.57735027をちょっと変えた(p.x+p.y+p.z)*0.57735027をdot()を使って書き換えるとdot(p,normalize(vec3(1)))となる。1/√3≒0.57735027なのでdot(p,vec3(1)*0.57735027)となる。
この辺りの説明は冗長になっちゃうと思うので感覚で感じてください。で、この応用でdot(p,normalize(vec3(1,2,3)))とかをfoldした空間に使うと面白いプリミティブが出来る。ここはロジックというより経験則です。

追記４

SDFの例外として使っているlength(cross(p,normalize(vec3(1))))ですが、これは原点を通るベクトルnormalize(vec3(1))のシリンダーになります。例外と言ってもlength()は使ってます。これはISFで使うと面白い表現が出来るので多用しています。このベクトルをnormalize(vec3(1,sin(time)*7.,2))とかをIFSの中で使うと楽しかったりします。cross(),dot()の用法については別記事で説明してます。以前に書いたモノで理解が足りてない所もありますが参考にはなると思います。
https://qiita.com/gaziya5/items/52ec06b5a7dd3b345d9e

You can support this creator by paying money.

Purchase disabled

Commercial use NG

Artwork of this product

fractal
View other products of this user