独自に定義した曲線関数 MCQBezierのいろいろ

はじめに

・MCQBezier(Midpoint-Controlled Quadratic Bezier Curve)は、指定する三点を通る二次ベジェ曲線です。

・私個人のプロジェクトで重要な役割をするのですが、いちいち紙のノートを探すのが面倒なのでここにまとめようと思います。

・Bezier曲線の入力を変更しただけで、実質的にBezier曲線と違いはありません。

定義

t = 割合, s = 始点, h = 中点, e = 終点 とする
$$
MCQBezier(t,s,h,e) = t^2 (2s + 2e - 4h) - t (3s + e - 4h) + s
$$

または A = 2s + 2e - 4h, B = A + s - eとし

$$
MCQBezier(t,s,h,e) = t^2 A - t B + s
$$

$$
\frac{d}{dt}MCQBezier(t,s,h,e) = 2tA - B
$$

導出

本来の二次ベジェ曲線をQBezierとする。QBezierの指定する定数の三点を、順番にs,b,eとすると、QBezierの定義は次のとおりである。
$$
QBezier(t,s,b,e)=s(1-t)^2 + 2bt(1-t)+et^2
$$
この曲線のt = 0.5 の時の点をhとする。bをhに入れ換えることを考える
$$
h = QBezier(0.5,s,b,e)= \frac{1}{4} s+ \frac{1}{2}b+\frac{1}{4}e
$$$$
2h = b+\frac{1}{2} s+\frac{1}{2}e
$$$$
b = 2h-\frac{1}{2} s-\frac{1}{2}e
$$

このbをQBezier(t,s,b,e)に代入した式が、MCQBezier(t,s,h,e)である。

長さ(近似)

ベジェ曲線の長さは微分の値(Δtに対する変化量)の長さの総和で求められる。
ガウス・ルジャンドル求積法を用いて0 <= t <= 1の範囲の長さを近似する。
$$
二次関数の近似に、必要な2つの定数は±\frac{1}{\sqrt{3}}
$$$$
d(t) = |2At-B|
$$$$
Length = \frac{1}{2}(d(\frac{1+\frac{1}{\sqrt{3}}}{2})+d(\frac{1-\frac{1}{\sqrt{3}}}{2}))
$$

点pに最も近いMCQBezier(t)

距離の微分を立式する

$MCQBezier(t)と点pの距離をL(t)とする$
$$
L(t) = \sqrt{(MCQBezier(t).x-p.x)^2 + (Yの項) + (Zの項)}
$$

$L(t)の微分が0になるとき、最近傍点となる(分子=0を解く)$
$$
\frac{d}{dt} L(t)の分子 = \frac{d}{dt}(MCQBezier(t).x - p.x)^2 + (Yの項) + (Zの項)
$$
$一旦、\frac{d}{dt}(MCQBezier(t) - p)^2を計算する。$

$$
MCQBezier(t) - p = t^2 A - tB+s-p より
$$$$
\frac{d}{dt}(MCQBezier(t) - p)^2 = t^3( 4A^2 ) + t^2( -6AB ) + t( -4Ap + 4As + 2B^2 ) + 2Bp - 2Bs
$$
以上の計算から
$\frac{d}{dt} L(t)の分子=$
$t^3(4A^2.x+ Yの項 + Zの項)$
$+t^2(-6AB.x+ Yの項 + Zの項)$
$+t( (-4Ap + 4As + 2B^2 ).x + Yの項 + Zの項)$
$+(2Bp - 2Bs).x + Yの項 + Zの項)$

これ=0を解の公式を使って解く

$a = 4A^2.x+ yzの項$
$b = -6AB.x+ yzの項$
$c = (-4Ap + 4As + 2B^2 ).x + yzの項$
$d = (2Bp - 2Bs).x + yzの項$

static float[] SolveCubicEquation(float a, float b, float c, float d)
        {
            if (a == 0.0)
            {
                Console.WriteLine("Error:a = 0.0"); 
                Console.WriteLine("This equation is NOT Cubic");
                return null;
            }
            float A, B, C, p, q, D,inv_a, inv_3,inv_27;
            inv_a = 1.0f / a;
            inv_3 = 0.3333f;
            inv_27 = 0.037037f;
            A = b * inv_a;
            B = c * inv_a;
            C = d * inv_a;
            p = B - A * A / 3.0f;
            q = 2.0f * A * A * A * inv_27 - A * B * inv_3 + C;
            D = q * q / 4.0f + p * p * p * inv_27;
            if (D < 0.0)
            {//three real solutions
                float theta = (float)Math.Atan2(Math.Sqrt(-D), -q * 0.5f);
                float[] result = new float[3];
                result[0] = (float)(2f * Math.Sqrt(-p * inv_3) * Math.Cos(theta *inv_3) - A *inv_3);
                result[1] = (float)(2f * Math.Sqrt(-p * inv_3) * Math.Cos((theta + 2f * Math.PI) * inv_3) - A * inv_3);
                result[2] = (float)(2f * Math.Sqrt(-p * inv_3) * Math.Cos((theta + 4f * Math.PI) * inv_3) - A * inv_3);
                foreach(var res in result)
                {
                    Console.WriteLine($"res {res}");
                }
                return result;
            }
            else
            {//single real solution and two imaginary solutions(c.c)
                float u =(float) Cbrt(-q * 0.5 + Math.Sqrt(D)), v = (float)Cbrt(-q * 0.5 - Math.Sqrt(D));
                float[] result = new float[1];
                result[0] = u + v - A *inv_3;
                Console.WriteLine($"res {result[0] }");
                return result;
            }
        }

四次関数の形から、3つの実数解をもつ場合はtが最小と最大の実数解のどちらかが最小のtとなる