LoginSignup
6
2

More than 5 years have passed since last update.

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

Last updated at Posted at 2017-03-09

modeling with distance functionsの距離関数の一覧に沿って記事を書いています.

image01.PNG

cupsuleの距離関数は,

vec3 pa = p - a, ba = b - a;
float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0);
return length( pa - ba*h ) - r;

数式に直すと

\left\{
\begin{array}{ll}
\vec{p}=(p_x, p_y, p_z), \vec{a}=(a_x, a_y, a_z), \vec{b}=(b_x, b_y, b_z) \\
h = min(max(\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}}, 0.0), 1.0) \\
\| (\vec{p} - \vec{a}) - h(\vec{b} - \vec{a}) \| - r = 0
\end{array}
\right.

記号の説明

\| \alpha \|:距離のこと \\
\vec{p}:ベクトル

とりあえず、なんでこれで、cupsuleになるのや?

と思いの方!

とりあえず、

\| (\vec{p} - \vec{a}) - h(\vec{b} - \vec{a}) \| - r

を考えます.

一旦, h は、面倒なので 1.0 とか 2.0 とかの定数として考えます..

すると

\| \vec{p} + (h-1)\vec{a} - h\vec{b} \| - r

で、また面倒なので、

\vec{p} + (h-1)\vec{a} - h\vec{b}  を  \vec{v}=(v_x, v_y, v_z)  とか適当に置きます.

で, 書き直すと

\sqrt{(p_x-v_x)^2 + (p_y-v_y)^2 + (p_z-v_z)^2} - r

となります.

なんとこれ hx,y,z に依存しない定数とみなすと球になるんですね.

ということは、予測ですが、cupsuleの両端の丸みは、球で描かれているのではないか?

と予測できます!

で、本題

cupsuleのあの胴長の円柱の鍵は、h にありです.

h = \min(\max(\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}}, 0.0), 1.0)

minmax を使わずきっちり書いてみましょう.

h = \left\{
\begin{array}{ll}
0.0 & (\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}} < 0) \\
\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}} & (0 < \frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}} < 1) \\
1.0 & (\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}} > 1)
\end{array}
\right.

となります.

ということは、ここでglslの距離関数に立ち返ると,

vec3 pa = p - a, ba = b - a;
float d = dot(pa,ba)/dot(ba,ba);
if(d < 0.0){
    return length(pa) - r;
} else if(d > 1.0){
    return length( pa - ba ) - r;
}
return length( pa - ba * d ) - r;

ということは,
それぞれ

return length(pa) - r;
return length( pa - ba ) - r;
return length( pa - ba * d ) - r;

を描画すると

return length(pa) - r;

image02.PNG

return length( pa - ba ) - r;

image03.PNG

return length( pa - ba * d ) - r;

image04.PNG

とま~、三つが組み合わさっているのがわかります.

で、最後に、

return length( pa - ba * d ) - r;

の形は、なんで円柱型になるの?を説明にして終わります.

length( pa - ba * d ) - r

数式で書くと

\| (\vec{p} - \vec{a}) - d(\vec{b} - \vec{a}) \| - r 

ここで、

\| (\vec{p} - \vec{a}) - (\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}})(\vec{b} - \vec{a}) \| - r   

ここからの計算は、頭が痛くなるので,(だいぶ緩い計算をします.)

\begin{align}
&\| (\vec{p} - \vec{a}) - (\frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}})(\vec{b} - \vec{a}) \| \\
&= \sqrt{\alpha_1 x^2 + \alpha_2 y^2 + \alpha_3 z^2 + \alpha_4 xy + \alpha_5 yz + \alpha_6 zx} \; (\alpha_n (n=1,2,3,4,5,6) は、適当にまとめた定数)
\end{align}

の形になり描画するのと円柱になります(ここら辺の計算書きたくないので、省略します…)

d の導出方法

(\vec{p} - \vec{a}) - d(\vec{b} - \vec{a}) と (\vec{b} - \vec{a}) が直行したと考えると \\
((\vec{p} - \vec{a}) - d(\vec{b} - \vec{a})) \cdot (\vec{b} - \vec{a}) = 0 なので、\\
hについてまとめると \\
d = \frac{(\vec{p} - \vec{a}) \cdot (\vec{b} - \vec{a})}{ \| \vec{b} - \vec{a} \|^{2}}

ちなみに d についての補足

d は書き方を変えると

\begin{align}
d &= \frac{|\vec{p} - \vec{a}||\vec{b} - \vec{a}| \cos \theta}{|\vec{b} - \vec{a}|^2} \\
&= \frac{|\vec{p} - \vec{a}| \cos \theta}{|\vec{b} - \vec{a}|} \\
\end{align} \\
と|\vec{p} - \vec{a}|の線分baへの射影と|\vec{b} - \vec{a}|の比とも見れますね.

dの変化を見てみる.

vec3 pa   = p - a, ba = b - a;
float ipr = dot(pa,ba)/(dot(ba,ba)*abs(sin(time)));
float h   = clamp(ipr, 0.0, 1.0);
return length( pa - ba*h ) - r;

motion04.gif

dot(ba,ba)*abs(sin(time)) の変化がいかに、円柱部分を作っていたか分かります.

コード

// ============================================================================
// Capsule / Line function_01
// ============================================================================

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

// Capsule / Lineの距離関数
float sdCapsule(vec3 p)
{
    // 回転
    // mat3 m_x = mat3(1,0,0,0,cos(time),-sin(time),0,sin(time),cos(time));
    // p = m_x * p;
    // mat3 m_y = mat3(cos(time),0,-sin(time),0,1,0,sin(time),0,cos(time));
    // p = m_y * p;
    // mat3 m_z = mat3(cos(time),-sin(time),0,sin(time),cos(time),0,0,0,1);
    // p = m_z * p;

    vec3 a  = vec3(0.75, 0.0, 1.0);
    vec3 b  = vec3(-0.75, 0.0, 1.0);
    float r = 1.0;

    /* 元の距離関数 */
    vec3 pa = p - a, ba = b - a;
    float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0);
    return length( pa - ba*h ) - r;

    /* 距離関数1 */
    // vec3 pa = p - a, ba = b - a;
    // float h = min(max( dot(pa,ba)/dot(ba,ba), 0.0), 1.0);
    // return sqrt((p.x-h*ba.x)*(p.x-h*ba.x)+(p.y-h*ba.y)*(p.y-h*ba.y)+(p.z-h*ba.z)*(p.z-h*ba.z))- r;

    /* 距離関数2 */
    // vec3 pa = p - a, ba = b - a;
    // float d = dot(pa,ba)/dot(ba,ba);
    // if(d < 0.0){
    //      return length(pa) - r;
    // } else if(d > 1.0){
    //      return length( pa - ba ) - r;
    // }
    // return length( pa - ba * d ) - r;

    /* dの変化を見る_1 */
    // vec3 pa   = p - a, ba = b - a;
    // float ipr = dot(pa,ba)/(dot(ba,ba)*abs(sin(time)));
    // float h   = clamp(ipr, 0.0, 1.0);
    // return length( pa - ba*h ) - r;

    /* dの変化を見る_2 */
    // vec3 pa = p - a, ba = b - a;
    // float ipr = (dot(pa,ba)/dot(ba,ba))*abs(sin(time));
    // float h   = clamp(ipr, 0.0, 1.0);
    // return length( pa - ba*h ) - r;

    /* dの変化を見る_3 */
    // vec3 pa = p - a, ba = b - a;
    // float ipr = dot(pa,ba)/dot(ba,ba);
    // float h   = clamp(abs(sin(ipr+time)), 0.0, 1.0);
    // return length( pa - ba*h ) - r;
}

// 距離関数を呼び出すハブ関数
float distanceHub(vec3 p){
    return sdCapsule(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))
    ));
}

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 < 32; ++i){
        dist = distanceHub(rPos);
        rLen += dist;
        rPos = cPos + ray * rLen;
    }
    // レイとオブジェクトの距離を確認
    if(abs(dist) < 0.001){
        // 法線を算出
        vec3 normal = genNormal(rPos);
        // ライトベクトルの定義
        vec3 light = normalize(vec3(1.0, 1.0, 1.0));
        // ライトベクトルとの内積を取る
        float diff = max(dot(normal, light), 0.1);
        // gl_FragColor = vec4(vec3(diff, diff, diff), 1.0);
        gl_FragColor = vec4(vec3(diff*177.0/255.0, diff*120.0/255.0, diff*68.0/255.0), 1.0);
    }else{
        // 衝突しなかった場合はそのまま黒
        gl_FragColor = vec4(vec3(0.0, 0.0, 0.0), 1.0);
    }
}

6
2
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
6
2