制約3:G1連続セグメントの接続点の傾きが連続であることを保証する制約です。これがないと、セグメントとセグメントの間で線が折れてしまいます。任意の$i$について、ある$\lambda_i\in(0,1)$があって、c_{i,2}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}が満たされれば、$c_{i,2}=c_{i+1,0}$における接線は$c_{i,1}$と$c_{i+1,1}$を結ぶ直線に定まります。κ-Curvesにおいてベジェ制御点は入力ではなく出力なので、これで全ての場合が網羅されています。c_{i,2}$を使わずに以下のように表現可能です。思うかもしれませんが、後で使います。制約4:ほぼG2連続セグメントの接続点の曲率が(ほぼ)連続であることを保証する制約です。\kappa_i(t)=\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|(1-t)(c_{i,1}-c_{i,0})+t(c_{i,2}-c_{i,1})\|^3}で、任意の$i$について$\kappa_i(1)=\kappa_{i+1}(0)$であればいいので、\begin{align}\kappa_i(1)&=\kappa_{i+1}(0)\\\\\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|c_{i,2}-c_{i,1}\|^3}&=\frac{\triangle(c_{i+1,0},c_{i+1,1},c_{i+1,2})}{\|c_{i+1,1}-c_{i+1,0}\|^3}\\\\\frac{\triangle(c_{i,0},c_{i,1},c_{i,2})}{\|c_{i,1}-c_{i+1,1}\|^3\lambda_i^3}&=\frac{\triangle(c_{i+1,0},c_{i+1,1},c_{i+1,2})}{\|c_{i,1}-c_{i+1,1}\|^3(1-\lambda_i)^3}\\\\\frac{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}{\lambda_i^2}&=\frac{\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})}{(1-\lambda_i)^2}\end{align}より(最後の式変形は図を書いて面積比に着目するとすぐ分かります)、\lambda_i=\frac{\sqrt{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}}{\sqrt{\triangle(c_{i,0},c_{i,1},c_{i+1,1})}+\sqrt{\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})}}が得られます。ところで、ベジェ曲線の曲率が0になることはないので、隣り合うセグメントの凹凸が逆である場合、曲率は必ず不連続になります。このとき、$\triangle(c_{i,0},c_{i,1},c_{i+1,1}),\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})$の符号が異なるので、$\lambda_i$は実数ではなくなります。一致させられないのであれば、せめて絶対値の差を0にしましょう。$|\kappa_i(1)|=|\kappa_{i+1}(0)|$を解くわけですが、$\kappa_i(1)$と$\kappa_{i+1}(0)$の符号は逆なので、$\kappa_i(1)=-\kappa_{i+1}(0)$を解けばいいことがわかります。すると、曲率が連続の場合とほぼ同様の手順で、以下が求まります。\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}これは先程の式も包含できているので、制約式としてはこちらのみを使えばよさそうです。p_i'=p_i-c_{i,0})\\\\(2)&c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}\\\\(3)&p_i=(1-\lambda_{i-1})(1-t_i)^2c_{i-1,1}+\big(\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\big)c_{i,1}+\lambda_it_i^2c_{i+1,1}\\\\(4)&\lambda_i=\displaystyle\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}\\\end{array}\right.が満たされるようにベジェ制御点$c_{i,j}$を定めればよいことになります。これをこのまま解析的に解くのは困難です。t_i$があれば、(3)で全ての$c_{i,1}$を求められる(方法は後述)全ての$c_{i,0},c_{i,1},c_{i,2}$があれば、(4)で全ての$\lambda_i$を求められるとなります。これらをうまく組み合わせて、全ての$p_i$から全ての$c_{i,j}$を出力したいわけです。適当な初期値から始めて、何度も式を適用することで正解に近づけていく方針を取ります。Step1.各ステップに分けて見ていきましょう。以下のように初期化します。c_{i+1,1})/2$つまり、ユーザ制御点と中央のベジェ制御点が同じ場所にあり、その他のベジェ制御点は隣接制御点の中点にある状態から始まります。これは論文から引用した図で、初期化時の状態の例です。黒い四角がユーザ制御点$p_i$で、$c_{i,1}$と一致しています。緑色の点は各セグメントの曲率極大点$c_i(t_i)$です。今はまだ$p_i$と$c_i(t_i)$が離れていますが、この後のStep1~4のイテレーションを回すことで近づけていきます。左から順に、初期状態、1ループ後、2ループ後、30ループ後(完全収束)です。式(4):\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}を適用します。式(2):c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}を計算します。p_i'=p_i-c_{i,0})を解きます。実解がただ一つ$[0,1]$に存在することが分かっています。カルダノの公式を愚直に組んでもいいですが(論文はそうしてるみたい?)、場合分けや例外処理が面倒なので私は二分探索しました。計測してみたら二分探索が処理時間の8割だったのでさすがにやめました。カルダノの公式を組みました。式(3):p_i=(1-\lambda_{i-1})(1-t_i)^2c_{i-1,1}+\big(\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\big)c_{i,1}+\lambda_it_i^2c_{i+1,1}を、全ての$c_{i,1}$について解きます。めんどくさポイントです。\gamma_i$として見やすく書き直しておきます。p_i=\alpha_ic_{i-1,1}+\beta_ic_{i,1}+\gamma_ic_{i+1,1}これは、以下のような連立方程式にまとめることができます。\left(\begin{matrix}\beta_0&\gamma_0&&&\alpha_0\\\alpha_1&\beta_1&\gamma_1&&\\&&\ddots&&\\&&\alpha_{n-2}&\beta_{n-2}&\gamma_{n-2}\\\gamma_{n-1}&&&\alpha_{n-1}&\beta_{n-1}\end{matrix}\right)\left(\begin{matrix}c_{0,1}\\c_{1,1}\\\vdots\\c_{n-2,1}\\c_{n-1,1}\\\end{matrix}\right)=\left(\begin{matrix}p_0\\p_1\\\vdots\\p_{n-2}\\p_{n-1}\\\end{matrix}\right)行列部分は三重対角行列+角なので、高速にLU分解でき、解を$O(n)$で求めることができます。これを上下に1つずつ拡張して\left(\begin{matrix}1&0&&&\\\alpha_0&\beta_0&\gamma_0&&\\&&\ddots&&\\&&\alpha_{n-1}&\beta_{n-1}&\gamma_{n-1}\\&&&0&1\end{matrix}\right)\left(\begin{matrix}c_{n-1,1}\\c_{0,1}\\\vdots\\c_{n-1,1}\\c_{0,1}\\\end{matrix}\right)=\left(\begin{matrix}c_{n-1,1}\\p_0\\\vdots\\p_{n-1}\\c_{0,1}\\\end{matrix}\right)という形にすれば、行列部分はただの三重対角行列になるので、さらに計算が楽になります。$n\ge5$の場合はメモリ的にも有利になります。(n+2)$要素分で済むのです。LU分解については記事の最後の付録で解説します。実際に実装していきましょう。言語はC#で、座標表現や各種演算にUnityのVector2クラス、Mathfクラスを借りています。Unity特有の何かがあるわけではないので、適宜好きな言語、好きなベクトル表現・数学ライブラリに置き換えてください。ベジェ制御点用構造体ただの配列のラッパーです。計算結果の出力もこのインスタンスで。配列サイズは最低限の$2n+1$ですが、ここまでの解説に合わせて[i,j]でアクセスできるようにしておきます。}}計算空間の確保ユーザ制御点が移動する度に描画を更新するわけなので、計算空間は事前に確保して使い回しましょう。制御点が増減すると確保し直しになりますが、その辺りを考慮するとコードが煩雑になるのでここでは妥協。}}ユーザ制御点・ベジェ制御点の上下拡張Step4の行列計算時にユーザ制御点ベクトルとベジェ制御点ベクトルを上下拡張しますが、その際のメモリ確保をなくしつつちゃんと配列っぽく扱えるようにするためのラッパー構造体です。コンストラクタで上下拡張時の値の初期化もやっています。本質部分ではないしC#の人じゃないとたぶん意味が分からないので適当に読み飛ばしてください。}}メソッドルート処理の根本はこんな感じでよいでしょう。isLoopを受け取り、Step0で初期化し、Step1~4でイテレーションを回し、最適化の結果を返します。各Stepの実装をしましょう。ループしない場合の対応も一緒にやっていきます。初期化内容はこうでした。(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}$非ループの場合、始端と終端はユーザ制御点であってほしいので、初期化時に固定してしまいましょう。$\lambda_0=0$$\lambda_{n-2}=1$$\lambda_{n-1}=\rm{undefined}$非ループの場合、終端→始端のカーブは必要なくなるので、ループの場合よりセグメントが2つ減ることに注意してください。その結果、$\lambda_{n-1}$は参照されなくなります。セグメントが1つではなく2つ減るのは直感的ではありませんが、実際に見れば納得できるでしょう。黄色・水色・ピンクの線が各セグメントのベジェ制御点を結んだものです。非ループの場合、水色とピンクの線が潰れているのが分かるでしょうか。c_{n-1,1}=p_{n-1}$となればよいので、A=\left(\begin{matrix}0&1&0\\0&1&0\\\alpha_1&\beta_1&\gamma_1\\&\vdots&\\\alpha_{n-2}&\beta_{n-2}&\gamma_{n-2}\\0&1&0\\0&1&0\\\end{matrix}\right)とします。λの算出\lambda_i=\frac{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}}{\sqrt{|\triangle(c_{i,0},c_{i,1},c_{i+1,1})|}+\sqrt{|\triangle(c_{i,1},c_{i+1,1},c_{i+1,2})|}}三角形の面積を求める解説は必要ないでしょう。非ループ時は始端と終端のλは更新しません。}}Sqrt計算を減らすために一応有理化して、分母がほぼ0のときは計算させず0.5にしています。ベジェ制御点(両端)の更新c_{i,2}=c_{i+1,0}=(1-\lambda_i)c_{i,1}+\lambda_ic_{i+1,1}やるだけです。$\lambda_i$の方で非ループ対応はしているので、ここでは特に何もしません。1];}※BezierControlsの実体は最後以外の$c_{i,2}$を削った配列なので、最後以外は片方だけに代入しています。p_i'=p_i-c_{i,0})を任意の方法で解きます。$[0,1]$に実解があることが分かっているので、ここでは実装が簡単な二分探索で。処理時間コストがバカ高かったのでカルダノの公式で組み直しました。}}参考:このページとこのページこれを使って、全ての$i$について三次方程式を解きます。//隣接3制御点が重なっている場合、a,b,c,dが全て0になり破綻するので既定値を返す。//今回の場合、a=0ならb=0なので、2次方程式の場合は考えなくてよい。}}座標表現をfloatでしていても、公式はdoubleで実装しないと桁落ちで破綻します。ベジェ制御点(中央)の更新\left(\begin{matrix}1&0&&&\\\alpha_0&\beta_0&\gamma_0&&\\&&\ddots&&\\&&\alpha_{n-1}&\beta_{n-1}&\gamma_{n-1}\\&&&0&1\end{matrix}\right)\left(\begin{matrix}c_{n-1,1}\\c_{0,1}\\\vdots\\c_{n-1,1}\\c_{0,1}\\\end{matrix}\right)=\left(\begin{matrix}c_{n-1,1}\\p_0\\\vdots\\p_{n-1}\\c_{0,1}\\\end{matrix}\right)を解きます。\begin{align}\alpha_i&=(1-\lambda_{i-1})(1-t_i)^2\\\beta_i&=\lambda_{i-1}(1-t_i)^2+(2-(1+\lambda_i)t_i)t_i\\\gamma_i&=\lambda_it_i^2\\\end{align}です。非ループの場合、$i=0,n-1$についてはStep0で初期化したように\begin{align}\alpha_0&=\alpha_{n-1}=0\\\beta_0&=\beta_{n-1}=1\\\gamma_0&=\gamma_{n-1}=0\\\end{align}で固定なので、上書きしないようにします。三重対角行列の連立方程式を解く関数を作っておきます。三重対角行列のLU分解について、詳細は記事の最後に付録として載せてあります。//対角要素の左隣の要素を対応するx(計算済み)にかけて引く。//対角要素の右隣の要素を対応するx(計算済み)にかけて引いて割る。}}あとはAを組み立ててユーザ制御点・ベジェ制御点を上下拡張して実行するだけです。一つ注意点として、Aはフルランクである必要があります。t_{i+1}=0$の場合にランクが落ちるので、微調整をかけます。ループしない場合、$t_{n-2}=1$や$t_1=0$でも同じことが起きます。extendedPs);}プロットこれでκ-Curvesのシステムは完成ですが、まだベジェ曲線の制御点が算出できただけなので、描画する必要があります。実際に画面に映すのは各描画ライブラリにやってもらうとして、そのための点群を用意しなければなりません。言っても、各セグメントについて、↓これにtを順番に突っ込めばいいだけです。c2;}計算スペースのときのように、プロット用のスペースも確保しておきましょう。セグメント数は非ループ時は2つ少なくなることに注意。}}あとは、頂点数が2つ以下のときの例外処理を挟みつつ、プロッティングしていくだけ。以上を全てKCurvesクラスに実装したとして、以下のようにすれば描画用の点群を取得できます。isLoop);お疲れさまでした。結果Unity上で、1セグメントにつき20ステップで描画してみた結果です。ベジェ制御点も表示してみるとこんな感じ。同じ配置でループさせるとこうなります。バグを修正(?)するはい、まだ終わってません。ここまでの実装で実際に描画してみて、ユーザ制御点をめちゃくちゃ動かしまくってみると分かりますが、特定の状況下において、適切な場所に収束しません。なんか、飛び出しています。見ると接続点の傾きも不連続になっています。尖った領域かつユーザ制御点が近接している場合に起こりがちです。この問題は、何故かよく分かりませんが、Step1の回数を抑えると抑制されます。傾きの連続性については、最後にStep2を一度実行することでとりあえず解消はされそうです。space.C;}これで直ります。必要な最適化を削っていることになるので、どれくらいを境にするかはお好みで。