※ 画面は開発中のものです: June, 2008

TKYM's profile※ 画面は開発中のものですPhotosBlogListsMore

Tools

Help

<< First < Previous Next > Last >>

June 29

C# でベジエ曲線（２）

無事、ベジエ曲線とその並行曲線がかけたのが前回まで。

曲線の式が分かることで、列車を滑らかにカーブさせることができそうだが、
まだ曲線の上を一定速度で動かすことができない。
前回のベジエ曲線の式は、パラメータ t で表されていて、
t を一定に増加させても等速度で動いてくれないように見える。

だから、適当な距離 l において、パラメータ t はいくつなのか調べないといけない。
また、曲線全体の距離も知りたい。

そこで、高校の数学教科書の最後のほうのページに、「曲線の長さ」という項目があったことを思い出し、開いてみた。

(b～a の定積分)∫√[ { f’(t) }^2 + { g’(t) }^2 ] dt

次数の高い多項式がルートの中に入ったものを積分するのは非常に難しそう。
つい最近、大学で積分の数値計算を習ったことから、それで計算することにする。

とはいっても、積分は基本的に、地道に面積を計上していくただのループ。
なんとかならないものかという気持ちが残る。

--------------

積分のアルゴリズムには、台形公式とシンプソンの公式と、ほかに各点の重みや積分区間を不等間隔にする何とかという公式があるらしい。

あるグラフを区間ごとに区切って面積を計算する時、その区間のグラフが曲線でないこととして計算すれば、台形公式になる。それでは誤差が大きくなってしまうので、区間上の３点を通る２次曲線の積分として考えれば、少しは誤差は小さくなる。これがシンプソンの公式である。

他にもいい方法があると本には書いてあるが、難しくなると面倒だし、ピクセル単位であっていればいいのでシンプソン法で計算してみることにした。

------------

とりあえず、曲線のパラメータ t での長さ l が何ピクセルなのかを知りたい。

計算のやり方は教科書のまる写しでなんとかする。
グラフを、 (t0, l0), (t1, l1), (t2, l2),,,,,(tN, lN) と区切って、
(h/3)(l0 + 4l1 + l2) + (h/3)(l2 + 4l3 + l4) + … と計算していく。

これを簡単にするとこんな感じ。
(h/3) * {
l0 + lN
+ l偶数番 * 2
+ l奇数番 * 4
}

ちなみに、区間は偶数個に区切らなければならないらしい。

        //パラメータ t を 開始点からの距離 l に変換する h は積分刻み幅
        //シンプソンの公式を用いた
        static public float t2l(float t, PointF[] pt, float h)
        {
            float tcnt, prev =0.0f, tmplen, tmpadd, total = 0.0f;
            int i;
            for (tcnt = 0.0f, i = 0; tcnt <= t; i++, tcnt += h) {
                tmplen = length(bezier.tangental_func(pt, tcnt));
 
                if (tcnt + h > t) 
                { //最後
                    if (i % 2 == 1)
                    { //奇数個で終わる
                        tmpadd = (prev + tmplen * 4 + length(bezier.tangental_func(pt, tcnt + h))) / 2;
                    }
                    else
                    { //偶数個で終わる
                        tmpadd = tmplen;
                    }
                }
                else if (tcnt < h)
                { //最初:そのまま
                    tmpadd = tmplen;
                }
                else if (i % 2 != 0)
                { //奇数
                    tmpadd = tmplen * 4;
                }
                else
                { //偶数
                    tmpadd = tmplen * 2;
                }
                total += tmpadd;
                prev = tmplen;
            }
            return h/3.0f * total;
        }

積分刻み幅を任意に定められることにしたかったので、項数が奇数になってしまった場合、区間を超えて計算し、超えた区間を２で割るようにしてみた。

そして、その逆変換

        //開始点からの距離 l を パラメータ t に変換する h は積分刻み幅
        //シンプソンの公式
        static public float l2t(float l, PointF[] pt, float h)
        {
            float tcnt, prev =0.0f, tmplen, tmpadd, total = 0.0f, lcnt, endlen = 0.0f;
            int i;
            for (tcnt = 0.0f, i = 0; tcnt <= 1.0f && endlen < l; i++, tcnt += h) {
                tmplen = length(bezier.tangental_func(pt, tcnt));
 
                //ここで終了した場合の長さ
                if (i % 2 == 1)
                { //奇数個で終わる
                    endlen = (total + (prev + tmplen * 4 + length(bezier.tangental_func(pt, tcnt + h))) / 2)
                        * (h / 3.0f);
                }
                else
                { //偶数個で終わる
                    endlen = (total + tmplen)
                         * (h / 3.0f);
                }

                //加算処理
                if (tcnt < h)
                { //最初
                    tmpadd = tmplen;
                }
                else if (i % 2 != 0)
                { //奇数
                    tmpadd = tmplen * 4;
                }
                else
                { //偶数
                    tmpadd = tmplen * 2;
                }

                total += tmpadd;
                prev = tmplen;
            }
            return tcnt - h/2.0f;
        }

距離から t を求める場合、ループするたびにその地点での距離を知りたいので、そのつど終了した場合の距離を計算し、終了判断しなければならない。あと、戻り値は必ず本来の値以上になることから、刻み幅の半分を引いてみた。

t = 0.5 の位置を距離に変換し、戻してみた。
大体あっているが、実用的なのか微妙なところ。

それより、コードがなんか汚い。言い訳としては、ご丁寧に中かっこをエディタが改行してくれるということ。
IntelliSense も C++ に比べるとおせっかいに働いてくれる。

曲線の接線に対して直角に、一定距離進んだところをつないで線にすれば平行になりそう。
接線を求めるには式を微分すればよさそう。

Bi と(n i) を定数として、 { (t^i) * (1-t)^(n-i)} ’ の部分を考える。
単純に考えると
(t^i)’ * (1-t)^(n-i) + t^i * {(1-t)^(n-i)}’ 積の微分
it^(i-1) * (1-t)^(n-i) + t^i * (n-i)(1-t)^(n-i-1)(-1) 合成関数

よって ∑ Bi (n i) { it^(i-1) * (1-t)^(n-i) + t^i * (n-i)(1-t)^(n-i-1)(-1) }

という感じになりそうなのだが、このままコードにすると失敗した。
0乗のときは1なので微分するとそのまま消えないといけないから。

        //ベジエ曲線の、接線を表す関数
        public static SizeF tangental_func(PointF[] pt, float t)
        {
            PointF total = new PointF(0.0f, 0.0f);
            int i;
            int n = pt.Length - 1;

            for (i = 0; i < pt.Length; i++)
            {
                float d1 = (i == 0) ? 0.0f : (float)i * powFI(t, i - 1);
                float d2 = (n - i == 0) ? 0.0f : (float)(n - i) * powFI(1.0f - t, n - i - 1) * (-1);

                total.X += pt[i].X * mCn(n, i) *
                (
                    d1 * powFI(1.0f - t, n - i) + powFI(t, i) * d2
                );


                total.Y += (pt[i].Y * mCn(n, i) *
                (
                    d1 * powFI(1.0f - t, n - i) + powFI(t, i) * d2
                ));
            }
            return new SizeF(total.X, total.Y);
        }

* powFI は float 型をキャストせずにべき乗できるようにした単純な関数。
* d1 と d2 がそれぞれ t^i と (1-t)^(n-i) を微分したものをいれる変数。

この戻り値をそのまま描画したらうまく接線が描画された。
しかし、場所によって長さが変わってしまった。曲線のパラメータ t の 1あたりの長さになってしまっているかららしい。

よって、このx, y に同じ値 l をかけて　長さを１にする必要がある。
そうしないと並行曲線が描けない。

(lx)^2 + (ly)^2 = 1^2
l^2 = 1 / (x^2 + y^2)
l = 1 / √(x^2 + y^2)

となり、(x / √(x^2 + y^2) , y / √(x^2 + y^2) ) が長さ１ピクセル分のベクトルになる。
英語のwikipedia に Parallel curve という項目があって、これを見ると何となく正しいような気がする。
http://en.wikipedia.org/wiki/Parallel_curve

この値に必要な長さをかけて使えばよい。

----------

接線を直角に曲げて、一定距離の地点を線でつないで行けば平行な線が引けそうである。
接線のベクトルが (x, y) なら、　その直角は ( y, –x )。

で、レール（道路）の幅だけ掛け算する。
これで平行線のグラフもわかった。

これらを使って適当に線を引く