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

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

二次元のグラフの描画は、[Desmos Graphing Calculator](https://www.desmos.com/) を使っています.

# 第一回目の記事の復習
[GLSLレイマーチング研究_距離関数について勉強してみた01(距離関数と球)](http://qiita.com/muripo_life/items/eb7ec78b08dc39beddd0)

早い話、楕円体って、球のx,y,zの軸の長さをそれぞれ変えたものです。  
なので、この記事で、球の風習の兼ねていこうと思います。

# そもそも二次元の楕円の方程式は?
```math
\frac{x^2}{a^2}+\frac{y^2}{b^2} = r^2 \\
min(|a|,|b|):短軸の長さ\\
max(|a|,|b|):長軸の長さ
```

[グラフはこんな感じ](https://www.desmos.com/calculator/ptzlrfjryw)

![iamge01.PNG](https://qiita-image-store.s3.amazonaws.com/0/102634/23c66866-1c2c-2616-0949-55acb5d72553.png)

で、

```math
|a|=|b|
```

の時に、円の方程式になんるんですね～

# 3Dの楕円体に拡張
```math
\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2} = r^2
```

# glslの距離関数にすると

```c++
float sdMathEllipsoid(in vec3 p, in vec3 r)
{
	return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
}
```

## min(min(r.x,r.y),r.z)は？
どうも、遠近感の調節のようです。
`min(min(r.x,r.y),r.z)` と `max(max(r.x,r.y),r.z)` の範囲で描画しているようです。

```c++
float sdStretchMathEllipsoid(in vec3 p, in vec3 r)
{
	return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * ((max(max(r.x,r.y),r.z)-min(min(r.x,r.y),r.z))*abs(sin(time))+min(min(r.x,r.y),r.z));
}
```
![motion01.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/9132616d-6df2-66f0-5cd7-eb8ce191ab41.gif)

# 双曲線
```c++
(sqrt(-p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z)
(sqrt(p.x/r.x*p.x/r.x-p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z)
(sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y-p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z)
```

![iamge02.PNG](https://qiita-image-store.s3.amazonaws.com/0/102634/27430e8c-cebc-6f28-6068-46c558fde9f7.png)

# ちょっとしたアニメーション(球でもやった)
```c++
(sqrt(p.x/r.x*p.x/r.x+sin(p.y/r.y*p.y/r.y+time)+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z)
```

![motion02.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/772f1df1-3369-a5ed-d2dd-c91401de81e9.gif)

# ノイズ?
```c++
(sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * sin(time)
```
![motion03.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/2890e5d8-5e20-a5f7-f7c2-58a8166e19cc.gif)

# 最後に、
縦横の比がある。  
トーラスを考える。
前のトーラス記事を見る。
[GLSLレイマーチング研究_距離関数について勉強してみた05(Torusの関数をいじる(距離関数導出編))](http://qiita.com/muripo_life/items/1736da4175028e3fb2b7)

```c++
float r = 0.3;
float R = 1.0;
return sqrt(p.x*p.x+p.y*p.y+p.z*p.z+R*R - 2.0 * R * sqrt(p.x*p.x+p.y*p.y) ) - r;
```

トーラスって、球に継ぎ足したものって、見ることが出来るんですよね。  
ということは、楕円形のトーラスを作ることが出来るんですよね。

```c++
// Torusの距離関数
float dTorus(vec3 p){
	float r = 0.3;
	float R = 3.0;
	vec3  s = vec3(0.2, 0.45, 0.6);
	return sqrt(p.x/s.x*p.x/s.x+p.y/s.y*p.y/s.y+p.z/s.z*p.z/s.z+R*R - 2.0 * R * sqrt(p.x/s.x*p.x/s.x+p.y/s.y*p.y/s.y))* min(min(s.x,s.y),s.z)- r;
}
```
![motion04.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/17153efa-089b-dc1f-f9f2-cbac884889de.gif)

[ソース](https://github.com/muripoLife/RayMarchingDistanceFunction/blob/master/18/fsp.frag)

```c++
// ============================================================================
// Ellipsoid
// ============================================================================

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 sdEllipsoid(in vec3 p, in vec3 r)
{
	return (length( p/r ) - 1.0) * min(min(r.x,r.y),r.z);
}

float sdMathEllipsoid(in vec3 p, in vec3 r)
{
	return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
}

float sdStretchMathEllipsoid(in vec3 p, in vec3 r)
{
	return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * ((max(max(r.x,r.y),r.z)-min(min(r.x,r.y),r.z))*abs(sin(time))+min(min(r.x,r.y),r.z));
}

// Torusの距離関数
float dTorus(vec3 p){
	float r = 0.3;
	float R = 3.0;
	vec3  s = vec3(0.2, 0.45, 0.6);
	return sqrt(p.x/s.x*p.x/s.x+p.y/s.y*p.y/s.y+p.z/s.z*p.z/s.z+R*R - 2.0 * R * sqrt(p.x/s.x*p.x/s.x+p.y/s.y*p.y/s.y))* min(min(s.x,s.y),s.z)- r;
}


float sdPlay(in vec3 p, in vec3 r)
{
	/*双曲線*/
	// return (sqrt(-p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
	// return (sqrt(p.x/r.x*p.x/r.x-p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
	// return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y-p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
	/*アニメ1*/
	// return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y*sin(time)+p.z/r.z*p.z/r.z) - 1.0) * min(min(r.x,r.y),r.z);
	// return sdEllipsoid(p, vec3(0.1+0.9*cos(time)*cos(time*2.0),0.1+0.9*cos(time)*sin(time*2.0),0.1+0.9*sin(time)));
	return (sqrt(p.x/r.x*p.x/r.x+p.y/r.y*p.y/r.y+p.z/r.z*p.z/r.z) - 1.0) * sin(time);
}


float sdInstanceEllipsoid(in vec3 p)
{
	// p = mat3(1.0,0,0, 0,cos(1.57),-sin(1.57), 0,sin(1.57),cos(1.57) )*p;
	p = mat3(1.0,0,0, 0,cos(time),-sin(time), 0,sin(time),cos(time) )*p;
 	// p = mat3(cos(time),0,-sin(time), 0,1,0, sin(time),0,cos(time))*p;
	p = mat3(cos(time),-sin(time),0, sin(time), cos(time),0 ,0,0,1)*p;
	// return sdEllipsoid(p, vec3(0.3, 0.45, 0.6));
	// return sdPlay(p, vec3(0.3, 0.45, 0.6));
	// return sdEllipsoid(p, vec3(0.1+0.9*cos(time)*cos(time*2.0),0.1+0.9*cos(time)*sin(time*2.0),0.1+0.9*sin(time)));
	return dTorus(p);
}

float distanceHub(vec3 p){
	return sdInstanceEllipsoid(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 (sdInstanceEllipsoid(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);
}
```