意外に難しいGLSL, SDFの関節ワザ

###GLSL, SDFの関節ワザ

***意外にに難しいSDFの関節ワザ***

 SDFを使ってモデリングし、特にSDFプリミティブを繋いで関節の動きのアニメーションを作成する場合、手順を誤ると***深い沼にはまってしまう***。
　SDFの種類によっても手順が異なるので注意が必要である。
ツリー構造やコードを見てもどこが問題なのかわかりにくいので、SDFプリミティブごとに実験する。

SDFのエディティングにはSDF　Editorを使用した https://joetech.itch.io/sdf-editor
がもちろん手でタイプしても実験できる。

先ず、一番わかりやすい"Capsule"から試してみる。(この方法はRoundConeでも有効)
下図のように、"Capsule"三つのツリー構造を作る。

　![スクリーンショット 2021-06-29 111108.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/42f3a06e-849c-9af3-4429-3843e966a84b.jpeg)

次に、それぞれのradiusに0.2を設定し、下二つのTranslationのy軸に1.5を設定する。
この***y軸方向に上につないでいくことがポイント***である。

これだけで、下の二つのRotationに任意の角度を設定すると、関節として問題なく表示できる。

![スクリーンショット 2021-06-29 111144.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/c4f59b7b-c17f-9df9-b525-2a7b8bc72df9.jpeg)


これをリアルタイムアニメーションで表示できるGLSL Sandboxのリンク↓

http://glslsandbox.com/e#73657.0

###楕円体 ellipsoidの場合

楕円体の場合、先ず以下のようなツリー構造を作り関節のルートとする。
（楕円体Icenuyまで、尚この名称はSDF Editorが自動的に設定する。）
楕円体のスケールに、
この場合最初の関節の位置の回転を制御するのは、一番上のRotationである。
また、上から二番目のTransltionのx軸に-1をインプットする。x = 1, y = 0.2  z = 0.2
とする。





![scrn2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/1c121d3d-5de7-0d23-bd21-d36e0cc30a09.jpeg)



次に、楕円体Icenuyの上のRotationの下にTranslationとRotationを追加し、さらにその下に楕円体(Etequs)を追加する。

そして、Etequsの上の二つのTranslationのx軸に-1を設定する。
同様に三つ目の楕円体を繋いだ場合の階層構造は以下のようになる。ここで、x軸に-1をインプットするTranslationは、下から一つ目のと二つ目である。つまり、一倍上のTranslation以外のすべてのTranslationのx軸に-1が入る。

![inu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/e3555b6d-1d4d-7694-052c-3aedf1ccd52a.jpeg)

この手順で、下図のような楕円体が三つ繋がった関節構造ができる。
ここで、関節の回転情報を入れるのは楕円体の直上のRotationではなく、そのさらに上の階層のRotationである。(直下のRotationは後で述べるように削除してもよい。）

![snap.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/4ebb2f28-97a9-5c32-5060-5eb01d916117.jpeg)

楕円体バージョンのリアルタイムアニメーション↓
GLSL Sandbox   http://glslsandbox.com/e#73662.0
ただし、楕円体の直上のRotationは今回のアニメーションには使っていなので、
取り除いて下図のようにしてもよい。
![new.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/430586/81de3473-6e07-717b-3bc4-693fb4a0aba1.jpeg)




以下はCapsuleバージョンのshadertoy.comで表示できるコード

```

const float epsilon = 0.01;
const float pi = 3.14159265359;
const float halfpi = 1.57079632679;
const float twopi = 6.28318530718;

mat3 rotateMat(vec3 p, float angle, vec3 axis){
    vec3 a = normalize(axis);
    float s = sin(angle);
    float c = cos(angle);
    float r = 1.0 - c;
    mat3 m = mat3(
        a.x * a.x * r + c,
        a.y * a.x * r + a.z * s,
        a.z * a.x * r - a.y * s,
        a.x * a.y * r - a.z * s,
        a.y * a.y * r + c,
        a.z * a.y * r + a.x * s,
        a.x * a.z * r + a.y * s,
        a.y * a.z * r - a.x * s,
        a.z * a.z * r + c
    );
    return m;
}

#define LIGHT normalize(vec3(1.0, 1.0, 0.0))

float displacement(vec3 p){
    return sin(20.141592*p.x)*sin(0.141592*p.y)*sin(20.131592*p.y);
}

//----------------------------------------------------------------------------------
//  from SDF Editor
//----------------------------------------------------------------------------------


float pCapsule(float r, float h, vec3 p)
{
    p.y -= clamp( p.y, 0.0, h );
    return length( p ) - r;
}

float pCylinder(float r, float h, vec3 p)
{
    vec2 d = abs(vec2(length(p.xz), p.y)) - vec2(r, h);
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}



vec3 mTranslation(vec3 inv_translation, vec3 p)
{
    return p + inv_translation;
}

vec3 mRotation(mat3 inv_rotation, vec3 p)
{
    return inv_rotation * p;
}


float oUnion(float d1, float d2)
{
    return min(d1, d2);
}


//-----------------------------------------------------------------
mat2 rotate(float a)
{
    float s=sin(a),c=cos(a);
    return mat2(c,s,-s,c);
}

mat3 lookat(in vec3 eye, in vec3 target)
{
	vec3 w = normalize(target-eye);
	vec3 u = normalize(cross(w,vec3(0.0,1.0,0.0)));
	vec3 v = normalize(cross(u,w));
    return mat3(u,v,w);
}

float lengthN(in vec2 p, in float n)
{
	//p = pow(abs(p), vec2(n));
    p = pow(abs(p), vec2(n));
	return pow(p.x+p.y, 1.0/n);
}
vec3 TrIn_1 = vec3(-0, -0, -0);
mat3 RoIn_2 = mat3(1, 0, 0, 0, 1, 0, 0, 0, 1);
float CaRa_3 = 0.2;
float CaHe_4 = 1.;
vec3 TrIn_5 = vec3(-0, -1.5, -0);
mat3 RoIn_6 = mat3(0.984808, -0.173648, 0, 0.173648, 0.984808, -0, -0, 0, 1);
float CaRa_7 = 0.2;
float CaHe_8 = 1.;
vec3 TrIn_9 = vec3(-0, -1.5, -0);
mat3 RoIn_10 = mat3(0.984808, -0.173648, 0, 0.173648, 0.984808, -0, -0, 0, 1);
float CaRa_11 = 0.2;
float CaHe_12 = 1.;

float sdf(vec3 p0)
{
	float d1;
	float d2;
	float d3;

	{
		vec3 p1 = mTranslation(TrIn_1, p0);
		{
            mat3 mt = rotateMat(p0,sin(iTime),vec3(0.0,0.0,1));
			vec3 p2 = mRotation(mt, p1);
			d1 = pCapsule(CaRa_3, CaHe_4, p2);
			{
				vec3 p3 = mTranslation(TrIn_5, p2);
				{
                     mat3 mt = rotateMat(p2,sin(iTime),vec3(0.0,0.0,1));
					vec3 p4 = mRotation(mt, p3);
					d2 = pCapsule(CaRa_7, CaHe_8, p4);
					{
						vec3 p5 = mTranslation(TrIn_9, p4);
						{
                            mat3 mt = rotateMat(p4,cos(iTime),vec3(0.0,0.0,1));
							vec3 p6 = mRotation(mt, p5);
							d3 = pCapsule(CaRa_11, CaHe_12, p6);
						}
					}
				}
			}
		}
	}
	return oUnion(d3, oUnion(d2, d1));
}
#define ZERO (min(iFrame,0))



vec3 calcNormal( in vec3 pos)
{

    vec2 e = vec2(1.0,-1.0)*0.5773*0.001;
    return normalize( e.xyy*sdf( pos + e.xyy) + 
					  e.yyx*sdf( pos + e.yyx) + 
					  e.yxy*sdf( pos + e.yxy) + 
					  e.xxx*sdf( pos + e.xxx) );

}

float calcOcclusion( in vec3 pos, in vec3 nor)
{
    float res = 0.0;
    vec3 aopos = pos;
    
    for( int i=0; i<4; i++ )
    {   
        aopos = pos + nor*0.2*float(i);
        float d = sdf(aopos);
        res += d;
    }

    return clamp(res, 0.0, 1.0);   
}



float castRay(in vec3 ro, vec3 rd)
{
   float t = 0.0;
    for(int i = 0; i < 100; i++ )
    {
        vec3 pos = ro+t*rd;
        float h = sdf(pos);
       
        if(h < 0.001) break;
            t += h;
        if(t > 20.) break;
        
    }
    if(t > 20.0) t = -1.0;
    return t;

}


vec3 env_color(vec3 dir)
{
	if (dir.y > 0.0)
		return mix(vec3(0.0, 0.5, 1.0), vec3(0.0, 0.1, 0.8), dir.y);
	else
		return mix(vec3(0.0, 0.5, 1.0), vec3(0.8, 0.7, 0.6), pow(-dir.y, 0.5));
}



void mainImage(out vec4 fragColor, in vec2 fragCoord)
{

    vec2 p = (2.0*fragCoord -iResolution.xy)/iResolution.y;
   

    float fi = 10.*iMouse.x/iResolution.x;
    float th = 10.*iMouse.y/iResolution.y;

    vec3 ro = vec3(5.0*cos(fi)*sin(th),1.0*sin(th)*sin(fi),8.0*cos(th));
    
    vec3 ta = vec3(0.0,0.0,0.0);
    vec3 ww = normalize(ta-ro);
    vec3 uu = normalize(cross(ww,vec3(0,1,0)));
    vec3 vv = normalize(cross(uu,ww));
    vec3 rd = normalize(p.x*uu+p.y*vv+1.5*ww);
    //vec3 col = vec3(0.4,0.75,0.9) - 0.5*rd.y;
    vec3 col = vec3(0.25,0.72,0.75) - 0.005*rd.y;
    //col = mix(col, vec3(0.7,0.75,0.8),exp(-10.0*rd.y));
    float t = castRay(ro,rd);
    
    if(t > 0.0){
       
        col = vec3(1.0);
        vec3 pos = ro + t *  rd;
        vec3 nor = calcNormal(pos);
        //float focc = res.w;
        //float focc = 2.5;  //  focc = clamp(4.0*res.z,0.0,1.0); -->candy
        float occ = calcOcclusion(pos, nor);//*focc;
        vec3 mate = vec3(0.18);
        vec3 sun_dir = normalize(vec3(0.8,0.4,0.2));
        float sun_sha = step(castRay(pos+nor*0.001, sun_dir), 0.0);
        float sun_dif = clamp(dot(nor,sun_dir),0.0,1.0);
        float sky_dif = clamp(0.5+0.5*dot(nor,vec3(0.0,1.0,0.0)),0.0,1.0);
        float bou_dif = clamp(0.5+0.5*dot(nor,vec3(0.0,-1.0,0.0)),0.0,1.0);
        col = mate*vec3(7.0,4.5,3.0)*sun_dif*sun_sha;
        col += mate*vec3(0.5,0.8,0.9) *sky_dif*occ;
        col += mate*vec3(0.7,0.3,0.2) *bou_dif*occ;

    } else {
        col = env_color(rd);
        //col = mix(col, vec3(0.7,0.75,0.9),exp(-10.0*rd.y));
    }
    col = pow(col,vec3(0.4545));
    fragColor = vec4(col,1.0);

}



```