GLSLレイマーチング研究_距離関数について勉強してみた17(Capped Coneの関数をいじる)

[modeling with distance functions](http://iquilezles.org/www/articles/distfunctions/distfunctions.htm)の距離関数の一覧に沿って記事を書いています.

# まず、今回初めて使う数学の関数の意味
[`sign`](https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/sign.xhtml)
シグナムと読む（大学時代は、シグナムと呼んでいた）。
* xが0.0より小さい場合は-1.0
* 0.0が0.0の場合は0.0
* xが0.0より大きい場合は+1.0
を返す関数.

たぶん、こんな感じの定義

```c++
float sign(float x){
	if(x>0.0){
		return 1.0;
	} else if(x=0.0){
		return 0.0;
	}
	return -1.0;
}
```

# Capped Coneの距離関数
```c++
float sdCappedCone( in vec3 p, in vec3 c )
{
    vec2 q = vec2( length(p.xz), p.y );
    vec2 v = vec2( c.z*c.y/c.x, -c.z );
    vec2 w = v - q;
    vec2 vv = vec2( dot(v,v), v.x*v.x );
    vec2 qv = vec2( dot(v,w), v.x*w.x );
    vec2 d = max(qv,0.0)*qv/vv;
    return sqrt( dot(w,w) - max(d.x,d.y) ) * sign(max(q.y*v.x-q.x*v.y,w.y));
}
```

ここで、 `vec3 c` は、 `vec3 c = vec3(【頂点の角度】, 【底面の半径】, 【高さ】);` で定義されている。

# 距離関数から元の数式を出してみる



```math
注釈:数式だと計算の時に掛け算と勘違いするので、vvはt、qvはs、と表記\\
c=(\theta, r, h)\\
q=(\sqrt{x^2+z^2}, y)\\
v=(\frac{rh}{\theta}, -h)\\
w=v-q=(\frac{rh}{\theta}-\sqrt{x^2+z^2}, -h-y)\\
t=(|v|^2,|\frac{rh}{\theta}|^2)=(\frac{r^2h^2}{\theta^2}+h^2, \frac{r^2h^2}{\theta^2})\\
s=(v \cdot w, \frac{rh}{\theta}\cdot(-h-y))=(\frac{r^2h^2}{\theta^2}+h^2-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy, -\frac{r}{\theta}h^2-\frac{rhy}{\theta})\\
\begin{align}
d&=(\max(\frac{r^2h^2}{\theta^2}+h^2-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy,0),\max(-\frac{r}{\theta}h^2-\frac{rhy}{\theta},0))
\left(
 \begin{array}{cc}
  1+\frac{-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy}{\frac{r^2h^2}{\theta^2}+h^2} & 0 \\
  0 & -\frac{h+y}{\frac{rh}{\theta}} \\
 \end{array}
\right)\\
&=(\max(\frac{r^2h^2}{\theta^2}+h^2-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy,0)(1+\frac{-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy}{\frac{r^2h^2}{\theta^2}+h^2}),-\max(-\frac{r}{\theta}h^2-\frac{rhy}{\theta},0)\frac{h+y}{\frac{rh}{\theta}})
\end{align}\\
で、距離関数を出すと\\
\sqrt{\frac{r^2h^2}{\theta^2}(x^2+z^2)+(h+y)^2-\max(\max(\frac{r^2h^2}{\theta^2}-\frac{rh}{\theta}\sqrt{x^2+z^2}+h^2+ry,0)(1+\frac{-\sqrt{x^2+z^2}\frac{rh}{\theta}+hy}{\frac{r^2h^2}{\theta^2}+h^2}),-\max(-\frac{r}{\theta}h^2-\frac{rhy}{\theta},0)\frac{h+y}{\frac{rh}{\theta}},0))}\\
+sign(\max(y\frac{rh}{\theta}+\sqrt{x^2+z^2}h, -h-y))\\
となります。
```
なげー

今回は、実装で生で数式かける気がしないので、いじっていくだけにします。

# `sign(max(q.y*v.x-q.x*v.y,w.y))`って何よ?

-1,0,1を返すなら`cos(time)`でも置き換えられんじぇね？

少し面白いアニメーションになりました。

```c++
return sqrt(dot(w,w) - max(d.x,d.y))*cos(time);
```
![motion01.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/2d9924d3-04ce-d549-4a81-ea6d49074ca5.gif)

それぞれの途中式をいじる.

```c++
	vec3 c = vec3(1.0, 1.0, 1.0);
// 	p = mat3(cos(time),0,-sin(time), 0,1,0, sin(time),0,cos(time))*p;
// 	p = mat3(1,0,0, 0,cos(time),-sin(time), 0,sin(time),cos(time))*p;

	vec2 q  = vec2(length(p.xz), p.y);
// 	vec2 q  = vec2(length(p.xz), p.y)*cos(time);
//      vec2 q  = vec2(length(p.xz)*cos(time), p.y);
//      vec2 q  = vec2(length(p.xz), p.y*cos(time));

	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
// 	vec2 v  = vec2(c.z*c.y/c.x, -c.z)*cos(time);
// 	vec2 v  = vec2(c.z*c.y/c.x*cos(time), -c.z);
// 	vec2 v  = vec2(c.z*c.y/c.x, -c.z*cos(time));


	vec2 w  = v - q;
// 	vec2 w  = (v - q)*cos(time);

	vec2 vv = vec2(dot(v,v), v.x*v.x);
// 	vec2 vv = vec2(dot(v,v), v.x*v.x)*cos(time);
// 	vec2 vv = vec2(dot(v,v)*cos(time), v.x*v.x);
// 	vec2 vv = vec2(dot(v,v), v.x*v.x*cos(time));

	vec2 qv = vec2(dot(v,w), v.x*w.x);
// 	vec2 qv = vec2(dot(v,w), v.x*w.x)*cos(time);
// 	vec2 qv = vec2(dot(v,w)*cos(time), v.x*w.x);
// 	vec2 qv = vec2(dot(v,w), v.x*w.x*cos(time));

	vec2 d  = max(qv,0.0)*qv/vv;
// 	vec2 d  = max(qv,0.0)*qv/vv*cos(time);
// 	vec2 d  = max(qv*cos(time),0.0)*qv/vv;
// 	vec2 d  = max(qv,sin(time))*qv/vv;

	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
```

# やっぱり数式に直していじるべき？

そこで、

```c++

/*できるだけ数式に直したもの*/
float sdMathCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y);
	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.z*p.z), -c.z-p.y);
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,0.0)*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
}
```

ま～式を書きやすい所だ直してみる。

```c

float sdPlayMathCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+(p.z*p.z)*cos(time)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x*cos(time)+(p.z*p.z)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+p.y*p.y*cos(time)+(p.z*p.z)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y*p.y*cos(time));	

	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
// 	vec2 v  = vec2(c.z*c.y/c.x, -c.z+(c.z*c.y/c.x)*cos(time));

	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.z*p.z), -c.z-p.y);
// 	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.y*p.y*cos(time)+p.z*p.z), -c.z-p.y);
	
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,0.0)*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
}
```

`q`に変化を加えたとき

```
vec2 q  = vec2(sqrt(p.x*p.x*cos(time)+(p.z*p.z)), p.y);
```

![motion02.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/20815b58-d239-9f16-8968-289757cb6caf.gif)

```c
vec2 v  = vec2(c.z*c.y/c.x, -c.z+(c.z*c.y/c.x)*cos(time));
```

![motion03.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/6842aedc-7f9a-7dd6-4849-013b067b2361.gif)

[コード](https://github.com/muripoLife/RayMarchingDistanceFunction/blob/master/17/fsp.frag)

```c
// ============================================================================
// Capped Cone
// ============================================================================

precision mediump float;
uniform vec2 resolution;     // resolution (512.0, 512.0)
uniform vec2 mouse;          // mouse(-1.0 ~ 1.0)
uniform float time;          // time(1second == 1.0)
uniform sampler2D prevScene; // previous scene texture

float sdCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(length(p.xz), p.y);
	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
	vec2 w  = v - q;
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,0.0)*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
}

/*できるだけ数式に直したもの*/
float sdMathCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y);
	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.z*p.z), -c.z-p.y);
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,0.0)*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
}

/*数式を変形したもの*/
float sdPlayMathCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+(p.z*p.z)*cos(time)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x*cos(time)+(p.z*p.z)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+p.y*p.y*cos(time)+(p.z*p.z)), p.y);
// 	vec2 q  = vec2(sqrt(p.x*p.x+p.z*p.z), p.y*p.y*cos(time));	

	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
// 	vec2 v  = vec2(c.z*c.y/c.x, -c.z+(c.z*c.y/c.x)*cos(time));

	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.z*p.z), -c.z-p.y);
// 	vec2 w  = vec2(c.z*c.y/c.x-sqrt(p.x*p.x+p.y*p.y*cos(time)+p.z*p.z), -c.z-p.y);
	
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,0.0)*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
}


float sdPlayCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	// p = mat3(cos(time),0,-sin(time), 0,1,0, sin(time),0,cos(time))*p;
	// p = mat3(1,0,0, 0,cos(time),-sin(time), 0,sin(time),cos(time))*p;

	vec2 q  = vec2(length(p.xz), p.y);
	// vec2 q  = vec2(length(p.xz), p.y)*cos(time);
	// vec2 q  = vec2(length(p.xz)*cos(time), p.y);
	// vec2 q  = vec2(length(p.xz), p.y*cos(time));

	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
	// vec2 v  = vec2(c.z*c.y/c.x, -c.z)*cos(time);
	// vec2 v  = vec2(c.z*c.y/c.x*cos(time), -c.z);
	// vec2 v  = vec2(c.z*c.y/c.x, -c.z*cos(time));


	vec2 w  = v - q;
	// vec2 w  = (v - q)*cos(time);

	vec2 vv = vec2(dot(v,v), v.x*v.x);
	// vec2 vv = vec2(dot(v,v), v.x*v.x)*cos(time);
	// vec2 vv = vec2(dot(v,v)*cos(time), v.x*v.x);
	// vec2 vv = vec2(dot(v,v), v.x*v.x*cos(time));

	vec2 qv = vec2(dot(v,w), v.x*w.x);
	// vec2 qv = vec2(dot(v,w), v.x*w.x)*cos(time);
	// vec2 qv = vec2(dot(v,w)*cos(time), v.x*w.x);
	// vec2 qv = vec2(dot(v,w), v.x*w.x*cos(time));

	vec2 d  = max(qv,0.0)*qv/vv;
	// vec2 d  = max(qv,0.0)*qv/vv*cos(time);
	// vec2 d  = max(qv*cos(time),0.0)*qv/vv;
	// vec2 d  = max(qv,sin(time))*qv/vv;

	return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y));
	// return -sqrt(dot(w,w) - max(d.x,d.y));
	// return sqrt(dot(w,w) - max(d.x,d.y)*sin(time)) * sign(max(q.y*v.x-q.x*v.y,w.y));
	// return sqrt(dot(w,w) - max(d.x,d.y));
	// return sqrt(dot(w,w) - max(d.x,d.y))*cos(time);
	// return sqrt(dot(w,w) - max(d.x,d.y)) * sign(max(q.y*v.x-q.x*v.y,w.y)*cos(time));
}

// 割と面白いアニメーション
float sdAnimePlayCappedCone(in vec3 p)
{
	vec3 c = vec3(1.0, 1.0, 1.0);
	vec2 q  = vec2(length(p.xz), p.y);
	vec2 v  = vec2(c.z*c.y/c.x, -c.z);
	vec2 w  = v - q;
	vec2 vv = vec2(dot(v,v), v.x*v.x);
	vec2 qv = vec2(dot(v,w), v.x*w.x);
	vec2 d  = max(qv,sin(time))*qv/vv;
	return sqrt(dot(w,w) - max(d.x,d.y))*cos(time);
}

float distanceHub(vec3 p){
    return sdCappedCone(p);
// 	return sdMathCappedCone(p);
// 	return sdPlayMathCappedCone(p);
// 	return sdPlayCappedCone(p);
	// return sdAnimePlayCappedCone(p);
}

// 法線
vec3 genNormal(vec3 p){
	float d = 0.001;
	// 法線を生成
	return normalize(vec3(
		distanceHub(p + vec3(	d, 0.0, 0.0)) - distanceHub(p + vec3( -d, 0.0, 0.0)),
		distanceHub(p + vec3(0.0,	 d, 0.0)) - distanceHub(p + vec3(0.0,	-d, 0.0)),
		distanceHub(p + vec3(0.0, 0.0,	 d)) - distanceHub(p + vec3(0.0, 0.0,	-d))
	));
}

// 図形ごとに色をわける
vec3 doColor(vec3 p){
	float e = 0.001;
	if (sdCappedCone(p)<e){
		vec3 normal	= genNormal(p);
		vec3 light	 = normalize(vec3(1.0, 1.0, 1.0));
		float diff	 = max(dot(normal, light), 0.1);
		return vec3(diff, diff, diff);
	}
// 	if (sdMathCappedCone(p)<e){
// 		vec3 normal	= genNormal(p);
// 		vec3 light	 = normalize(vec3(1.0, 1.0, 1.0));
// 		float diff	 = max(dot(normal, light), 0.1);
// 		return vec3(0.0, diff, diff);
// 	}
//     if (sdPlayMathCappedCone(p)<e){
// 		vec3 normal	= genNormal(p);
// 		vec3 light	 = normalize(vec3(1.0, 1.0, 1.0));
// 		float diff	 = max(dot(normal, light), 0.1);
// 		return vec3(0.0, diff, diff);
// 	}
// 	if (sdPlayCappedCone(p)<e){
// 		vec3 normal	= genNormal(p);
// 		vec3 light	 = normalize(vec3(1.0, 1.0, 1.0));
// 		float diff	 = max(dot(normal, light), 0.1);
// 		return vec3(0.0, diff, diff);
// 	}
// 	if (sdAnimePlayCappedCone(p)<e){
// 		vec3 normal	= genNormal(p);
// 		vec3 light	 = normalize(vec3(1.0, 1.0, 1.0));
// 		float diff	 = max(dot(normal, light), 0.1);
// 		return vec3(diff, diff, diff);
// 	}
	return vec3(0.0);
}

// カメラのワーク
void main(){
	// スクリーンスペースを考慮して座標を正規化
	vec2 p = (gl_FragCoord.xy * 2.0 - resolution) / min(resolution.x, resolution.y);
	// カメラを定義
	vec3 cPos  = vec3(0.0,	0.0,	3.0);   // カメラの位置
	vec3 cDir  = vec3(0.0,	0.0, -1.0);     // カメラの向き(視線)
	vec3 cUp   = vec3(0.0,	1.0,	0.0);   // カメラの上方向
	vec3 cSide = cross(cDir, cUp);          // 外積を使って横方向を算出
	float targetDepth = 1.0;                // フォーカスする深度

	// カメラの情報からレイを定義
	vec3 ray = normalize(cSide * p.x + cUp * p.y + cDir * targetDepth);
	// マーチングループを組む
	float dist = 0.0;	// レイとオブジェクト間の最短距離
	float rLen = 0.0;	// レイに継ぎ足す長さ
	vec3 rPos = cPos; // レイの先端位置(初期位置)

	// レイが進む処理(マーチングループ)
	for(int i = 0; i < 32; ++i){
		dist = distanceHub(rPos);
		rLen += dist;
		rPos = cPos + ray * rLen;
	}

	// レイとオブジェクトの距離を確認
	vec3 color = doColor(rPos);
	gl_FragColor = vec4(color, 1.0);
}
```