GLSLレイマーチング研究_距離関数について勉強してみた23(sigmoid関数でスライム的なものを作ってみた。)

# 参考URL

* [図のURL](https://www.desmos.com/calculator/xaghf5x2fg)
* [シグモイド関数の意味と簡単な性質](http://mathtrain.jp/sigmoid)
* [シグモイド](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%82%B0%E3%83%A2%E3%82%A4%E3%83%89)
* [シグモイド関数](http://d.hatena.ne.jp/Zellij/20111011/p1)
* [Obey Your MATHEMATICS](http://mathetake.hatenablog.com/entry/2016/11/03/194454)
* [CHAPTER 4ニューラルネットワークが任意の関数を表現できることの視覚的証明](https://nnadl-ja.github.io/nnadl_site_ja/chap4.html)

もともとは、スライム的なものを作る気は、さらさらなく…  
機械学習の関数をレイマーチで出してみようと思ったからでした。

# シグモイド関数って?
生物の神経細胞が持つ性質をモデル化したものとして用いられ、今流行の機械学習で、用いられる関数です。

```math
f(x) = \frac{1}{1+\exp(-ax)}
```

![motion02.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/70391573-5c06-5446-8f7c-fe95b76a2066.gif)

アニメーションをみたらわかるように、係数 `a` の値を大きくしていくと…  
カーブが急になっていることに気づきます。  
ということは、ステップ関数（step function）を滑らかなにしたもの

個人的には、ヘビサイド関数のようなイメージでしたが…  

# Sigmoid関数を定義する.

```glsl
float sigmoid(float x, float k)
{
    return 1./(1.+exp(-(k*x)));
}
```

# とりあえず、描画してみる.

```glsl
float dInstanceSigmoid(vec3 p){
    /*
    ** オーソドックスに書いてみる
    ** f(x) = 1/(1+exp(-x), 0<f(x)<1で描画
    */
    // -------------------------------------------------------------------
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数を変化させてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 10.*sin(time)
    */
    // -------------------------------------------------------------------
    // float a = 10.*sin(time*0.1);
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数の変化を正だけにさせてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 10.*abs(sin(time))
    */
    // -------------------------------------------------------------------
    // float a = 10.*abs(sin(time*0.1));
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数の変化を30まで正だけにさせてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 30.*abs(sin(time))
    */
    // -------------------------------------------------------------------
    // float a = 30.*abs(sin(time*0.1));
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** オーソドックスに書いてみる
    ** f(x) = 1/(1+exp(-x), 0<f(x)<10で描画
    */
    // -------------------------------------------------------------------
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-10.0;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)-1.0;
    // -------------------------------------------------------------------

    /*
    ** f(x) = 1/(1+exp(-x)
    ** 0<f(x)<k*abs(sin(time)), k =10
    */
    // -------------------------------------------------------------------
    // float k = 10.*abs(sin(time*0.1));
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-k;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-k;
    // return sigmoid(p.x, 1.)-k;
    // -------------------------------------------------------------------

    /*
    ** 合わせてみる
    ** f(x) = 1/(1+exp(-x)
    ** 0<f(x)<k*abs(sin(time)), k =10
    */
    // -------------------------------------------------------------------
    float k = 10.*abs(sin(time*0.1));
    float a = 10.*abs(sin(time*0.1));
    return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-k;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-k;
    // return 1./(1.+exp(-a*p.x))-k;
    // -------------------------------------------------------------------
}
```

![motion03.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/457164a2-767f-dc12-cf64-7a0e398f0c77.gif)

なんか微妙だな…  
いろいろな条件で試したものがあるので、コメント外して試してみてください.

# とりあえず、平面にかませてみる

```glsl
float dSequenceSigmoid(vec3 p){
    // sin刻みでSigmoid関数を慣べてみる
    // float s1 = sigmoid(p.x, 1.) + sin(p.z);
    float s = sigmoid(p.x, 32.*sin(time/3.)) + sin(p.z);
    // 高さがわかりやすように強調する
    s = s*1.5;
    // ベースラインを下げてからシグモイド関数で持ち上げてみる.
    return dot(p, vec3(0.0, 1.0, 0.0)) + 4.0 - s;
}
```

![motion04.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/053faeff-5cb0-b77f-3d4e-fbbcdb32b30e.gif)

これは、予想通り

# 球の距離関数にかませてみる

```glsl
float dSigmoidSphere(vec3 p){
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 32.*abs(sin(time/3.))), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x*sin(time), 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y*sin(time), 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z*sin(time), 1.));
    p = vec3(sigmoid(p.x*sin(time), 5.*abs(sin(time/3.)))+cos(time), sigmoid(p.y*sin(time)+sin(time), 5.*abs(sin(time/3.))), sigmoid(p.z, 1.));
    return length(p) - 1.0;
}
```

![motion05.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/cc1dd3e3-724c-c7b3-ae2f-327af1baacf4.gif)

# トーラスの距離関数にかませてみる

```glsl
float dSigmoidTorus(vec3 p){
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x*sin(time), 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y*sin(time), 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z*sin(time), 1.));
    p = vec3(sigmoid(p.x*sin(time), 1.)+cos(time), sigmoid(p.y*sin(time)+sin(time), 1.), sigmoid(p.z, 1.));

    float r = 0.3;
    float R = 1.0;
    vec2 q = vec2(length(p.xy)-R,p.z);
    return length(q)-r;
}
```

![motion06.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/9f641ce7-4482-9525-a703-c64fe2a8b295.gif)

# あとは、適当に混ぜてみたら、スライム的なモノができました。

```glsl
return (dSphere(p)*(0.1+abs(cos(time/3.)*cos(time/2.)))+dSigmoidSphere(p)*(0.1+abs(cos(time/3.)*sin(time/2.)))+dSigmoidTorus(p)*(0.1+abs(sin(time/3.))))-0.3;
```

![motion01.gif](https://qiita-image-store.s3.amazonaws.com/0/102634/d6241621-00ad-1814-1c4a-0d2c54f68793.gif)

コード

```glsl
// ============================================================================
// sigmoid distance function
// f(x) = 1/(1+exp(-ax))
// ============================================================================

precision mediump float;
uniform vec2  resolution;
uniform float time;
uniform sampler2D prevScene;

mat3 rot_x(float arg){
    return mat3(
        1.0,0,0,
        0,cos(arg),-sin(arg),
        0,sin(arg),cos(arg)
    );
}

mat3 rot_y(float arg){
    return mat3(
        cos(arg),0,-sin(arg),
        0,1,0,
        sin(arg),0,cos(arg)
    );
}

mat3 rot_z(float arg){
    return mat3(
        cos(time),-sin(time),0,
        sin(time), cos(time),0
        ,0,0,1
    );
}

float sigmoid(float x, float k)
{
    return 1./(1.+exp(-(k*x)));
}

float sigmoid_atan(float x, float k)
{
    return (1.26*atan(0.315*k*x)+1.)/2.;
}

float dSphere(vec3 p){

    return length(p)-1.;
}

float dSequenceSigmoid(vec3 p){
    // sin刻みでSigmoid関数を慣べてみる
    // float s1 = sigmoid(p.x, 1.) + sin(p.z);
    float s = sigmoid(p.x, 32.*sin(time/3.)) + sin(p.z);
    // 高さがわかりやすように強調する
    s = s*1.5;
    // ベースラインを下げてからシグモイド関数で持ち上げてみる.
    return dot(p, vec3(0.0, 1.0, 0.0)) + 4.0 - s;
}

float dSigmoidSphere(vec3 p){
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 32.*abs(sin(time/3.))), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x*sin(time), 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y*sin(time), 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z*sin(time), 1.));
    p = vec3(sigmoid(p.x*sin(time), 5.*abs(sin(time/3.)))+cos(time), sigmoid(p.y*sin(time)+sin(time), 5.*abs(sin(time/3.))), sigmoid(p.z, 1.));
    return length(p) - 1.0;
}

float dSigmoidTorus(vec3 p){
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x*sin(time), 1.), sigmoid(p.y, 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y*sin(time), 1.), sigmoid(p.z, 1.));
    // p = vec3(sigmoid(p.x, 1.), sigmoid(p.y, 1.), sigmoid(p.z*sin(time), 1.));
    p = vec3(sigmoid(p.x*sin(time), 1.)+cos(time), sigmoid(p.y*sin(time)+sin(time), 1.), sigmoid(p.z, 1.));

    float r = 0.3;
    float R = 1.0;
    vec2 q = vec2(length(p.xy)-R,p.z);
    return length(q)-r;
}


float dInstanceSigmoid(vec3 p){
    /*
    ** オーソドックスに書いてみる
    ** f(x) = 1/(1+exp(-x), 0<f(x)<1で描画
    */
    // -------------------------------------------------------------------
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数を変化させてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 10.*sin(time)
    */
    // -------------------------------------------------------------------
    // float a = 10.*sin(time*0.1);
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数の変化を正だけにさせてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 10.*abs(sin(time))
    */
    // -------------------------------------------------------------------
    // float a = 10.*abs(sin(time*0.1));
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** 係数の変化を30まで正だけにさせてみる
    ** f(x) = 1/(1+exp(-ax)), 0<f(x)<1で描画
    ** a = 30.*abs(sin(time))
    */
    // -------------------------------------------------------------------
    // float a = 30.*abs(sin(time*0.1));
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-1.;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-1.0;
    // return 1./(1.+exp(-a*p.x))-1.0;
    // -------------------------------------------------------------------

    /*
    ** オーソドックスに書いてみる
    ** f(x) = 1/(1+exp(-x), 0<f(x)<10で描画
    */
    // -------------------------------------------------------------------
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-10.0;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-1.0;
    // return sigmoid(p.x, 1.)-1.0;
    // -------------------------------------------------------------------

    /*
    ** f(x) = 1/(1+exp(-x)
    ** 0<f(x)<k*abs(sin(time)), k =10
    */
    // -------------------------------------------------------------------
    // float k = 10.*abs(sin(time*0.1));
    // return sigmoid(p.x, 1.)+sigmoid(p.y, 1.)+sigmoid(p.z, 1.)-k;
    // return sigmoid(p.x, 1.)+sigmoid(p.z, 1.)-k;
    // return sigmoid(p.x, 1.)-k;
    // -------------------------------------------------------------------

    /*
    ** 合わせてみる
    ** f(x) = 1/(1+exp(-x)
    ** 0<f(x)<k*abs(sin(time)), k =10
    */
    // -------------------------------------------------------------------
    float k = 10.*abs(sin(time*0.1));
    float a = 10.*abs(sin(time*0.1));
    return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.y))+1./(1.+exp(-a*p.z))-k;
    // return 1./(1.+exp(-a*p.x))+1./(1.+exp(-a*p.z))-k;
    // return 1./(1.+exp(-a*p.x))-k;
    // -------------------------------------------------------------------
}

float distanceHub(vec3 p){
    // p = rot_x(time) * p;
    // p = rot_y(time) * p;
    // p = rot_z(time) * p;
    // p = vec3(p.x+cos(time)*cos(time/2.), p.y+cos(time)*sin(time/2.), p.z+sin(time));


    // return dInstanceSigmoid(p);
    // return dSequenceSigmoid(p);
    // return dSigmoidSphere(p);
    // return dSigmoidTorus(p);
    // return dSigmoidSphere(p)*abs(sin(time))+dSigmoidTorus(p)*(1.-abs(sin(time)));
    // return (dSphere(p))*(0.1+abs(cos(time/3.)*cos(time/2.)))+dSigmoidSphere(p)*(0.1+abs(cos(time/3.)*sin(time/2.)))+dSigmoidTorus(p)*(0.1+abs(sin(time/3.)));
    return (dSphere(p)*(0.1+abs(cos(time/3.)*cos(time/2.)))+dSigmoidSphere(p)*(0.1+abs(cos(time/3.)*sin(time/2.)))+dSigmoidTorus(p)*(0.1+abs(sin(time/3.))))-0.3;
    
}

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 (distanceHub(p)<e){
        vec3 normal = genNormal(p);
        vec3 light  = normalize(vec3(1.0, 1.0, 1.0));
        float diff  = max(dot(normal, light), 0.1);
        float spec = pow(diff*diff, 15.0);
        return vec3(diff*0.3+spec, diff+spec, diff+spec);
        // return vec3(diff*cos(time)*cos(time/2.)+spec*abs(0.5+0.5*sin(time)), diff*cos(time)*sin(time/2.)+spec*abs(0.5+0.5*sin(time)), diff*sin(time)+spec*abs(0.5+0.5*sin(time)));

    }
    // if (dSigmoidSphere(p)<e){
    //     vec3 normal = genNormal(p);
    //     vec3 light  = normalize(vec3(1.0, 1.0, 1.0));
    //     float diff  = max(dot(normal, light), 0.1);
    //     float spec = pow(diff*diff, 15.0);
    //     return vec3(diff*0.5+spec, diff+spec, diff+spec);
    // }
    // if (dSigmoidTorus(p)<e){
    //     vec3 normal = genNormal(p);
    //     vec3 light  = normalize(vec3(1.0, 1.0, 1.0));
    //     float diff  = max(dot(normal, light), 0.1);
    //     float spec = pow(diff*diff, 15.0);
    //     return vec3(diff+spec, diff*0.5+spec, diff+spec);
    // }
    // if (dSphere(p)<e){
    //     vec3 normal = genNormal(p);
    //     vec3 light  = normalize(vec3(1.0, 1.0, 1.0));
    //     float diff  = max(dot(normal, light), 0.1);
    //     float spec = pow(diff*diff, 15.0);
    //     return vec3(diff*0.5+spec, diff+spec, diff*0.5+spec);
    // }
    return vec3(0.);
}

void main(){
    vec2 p = (gl_FragCoord.xy * 2.0 - resolution) / min(resolution.x, resolution.y);
    vec3 cPos  = vec3(0.0, 0.0,  5.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 < 64; ++i){
    // 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);
}
```