ちょっと待て！彼女を生成するってやつはどうなったんだ？ って人は、前半をご参照ください。
ウェブアプリを作ろうと思って、気づいたらOS作ってた とか、
3Dアプリを開発しようと思って、気づいたらサーバー組んでた とか、
3Dで彼女をプロシージャル生成しようと思って、気づいたら最適化数学で曲線をフィッティングしていた とか、
我々プログラマには良くあることです。
------------------------------
ペンやマウスで描画すると、自動的にきれいなベジェ曲線にしてくれるドローアプリでは、
Blender使いご用達のInkScapeというアプリがあります。
このような、ある頂点群（たとえばペンで描画した描画点や、メッシュ上の選択した頂点）に対して、
ベジェ曲線をフィッティングする問題を考えます。
といってもそんな頭良くないので、InkScapeの実装の元ネタである Graphics Gems 1 の
「An Algorithm for Automatically Fitting Digitized Curves」を理解しつつ実装することにします。
このアルゴリズムは、2つの手法から成り立っています。
1. ニュートン法 ( Newton-Raphson method )
2. 最小二乗法 ( least squares method )
このうち、2のほうがメインのフィッティング処理です。
というわけで2のほうだけやります。
------------------------------
といっても難しいので、まずは2次元で考えます。
こんな感じで2D平面上に頂点があるとき、全頂点にできるだけ近い直線を引くような問題を解くのが最小二乗法です。
具体的には、直線と各点の距離(誤差)の二乗の合計が、最小になるような直線を求めます。
直線を y = ax + b として、x上のある値 xα での、頂点を yα とすると、
yαと、xαでの直線上の点との差分は、
yα-(axα + b) です。
距離の二乗和は、
となり、この式が最小になるような、a、bを求めます。
求め方ですが、上記の∑の式の中身は距離なので、a, bをそれぞれ変化させたときに、それぞれ極小値があります。 よって、上記の∑の式を J と置くと、
となります。これを計算していくと、最後に連立一次方程式となり、解くことができます。
このあたりの詳細は、参考文献(1)が詳しいです。この辺は、結構ウェブでもあります。
------------------------------
次にベジェ曲線を、ベジェ曲線について、少し書きます。
ベジェ曲線は、3次のものしか使わないので、ここでは3次とします。
3次のベジェ曲線は、アンカー２つと、その間のハンドル２つでなります。
ベジェ曲線の数式は、次のように定義されています。
ここで Vi は i番目のベジェ制御点で、実際には3Dなので V(xi, yi, zi) のようなベクトルになります。
i は、0から n まで入ります。 また、Biは、バーンスタイン基底関数とか呼ばれている、0以上1以下の変数 t の関数です。
この変数 t に 0から1までの数を次々に入れると、ベジェ曲線ができます。
n はベジェ曲線の次数で、2次だったら、アンカー2点とハンドル1つ、3次だったらアンカー2つとハンドル2つとなりますが、
今回は3次なので、nに3を入れておきます。
上から、B0, B1, B2, B3 と置いて、pythonの関数で書くと、こうなります。
def B0(t):
return (1.0-t) * (1.0-t) * (1.0-t)
def B1(t):
return 3 * t * (1.0-t) * (1.0-t)
def B2(t):
return 3 * t * t * (1.0-t)
def B3(t):
return t * t * t
また、今回実装予定のアルゴリズムでは、ハンドルはアンカーの接線上に配置します。そうすることで、曲線の連続性が得られ、前後で滑らかに繋がります。
------------------------------
さて、ベジェ曲線は、ご覧のとおり、t による関数なので、直線のときみたいに x とか y とかを使って最小二乗法するわけには行きません。
また、何を最小にするかも一工夫必要です。
任意の頂点群にできるだけ近いベジェ曲線の、ベジェ制御点 V0, V1, V2, V3を求めたいわけですが、
V0とV3については、頂点群がペンでの描画と考えた場合、最初と最後の点で良さそうです。
残るは、V1, V2ですが、
上で、ハンドルはアンカーの接線上に配置するという制約を付けたので、アンカー接線上にあります。
V0、V3の、アンカー接線の単位ベクトルを T0、T3とすると、
V1、V2は、T0、T3の定数倍をV0, V3に足すことで表させます。
つまり、
V1 = V0 + a T0 (aは定数)
V2 = V3 + b T3 (bは定数)
となります。
また、最小にする値は、頂点群のうちの、任意の頂点 P に対応する t 軸上の点を tαとしたとき、
ベジェ曲線上の点 Q(tα) と頂点 P との距離になります。
任意の頂点 P に対応する t 軸上の点を tα は、どうやって出すかというと、
たぶん 0から1の値をとる t を、適当に頂点群の個数で割った値、とかにしてもいいと思うんですが、
この実装では、
tα = 頂点群の最初から頂点Pまで繋げた長さ / 頂点群を最初から最後まで繋げた全体の長さ
としているようです。適当に割るよりは良さそうです。
さて、以上をまとめると、最小二乗法の最小にしたい式が出来ます。
この式が最小になるような、アンカー接線上の点 V1, V2 を求めるんですが、
上に書いた通り、V1, V2 は 単位接線ベクトルの定数倍なので、その定数 a, bを求めれば、求まります。
a, b を変化させることにより、J、つまり距離(誤差)の二乗和の合計が変化し、それぞれ極小値があるので、
が成り立ちます。で、例によってこの先は計算の嵐なので、省略しますが、最終的に２次の場合と同じく、
連立一次方程式が出来て、その解を求めることで、a、b が求まります。
この辺余力があれば書こうと思ったんですが、もう時間無いので、残りの式が気にある方は自分で出すか、
Graphics Gems 1 を入手してください。
------------------------------
このように最小二乗法によるベジェフィッティングを理解しつつ、
Graphics Gems1 のソースコード（Public Domain)の、最小二乗法のところだけ抜き出して移植したのがこちら。
とりあえずGreasePencilで描いたもの一番最後のやつを変換するようにしてますが、今回のアルゴリズムだけでは
カーブは１個しかできないので、あまり使い物にならないです。
そこで次回 1. ニュートン法による近似 を実装して、マシになる予定です…
あと、特殊な場合については無視してるので場合によっては変換上手くいかないです。
GreasePencilのストロークは、最初の点とその次の点で接線を出そうとしたところ、常に(0,1,0)とかになって
最初と最後付近の頂点が変な感じになってることが判明したので、適当に2,3個の点から作ってる程度にやっつけです。
import bpy
from mathutils import Vector
def make_curve(bezier_list):
curvedata = bpy.data.curves.new("curve", type='CURVE')
curvedata.dimensions = '3D'
objectdata = bpy.data.objects.new("bezier", curvedata)
bpy.context.scene.objects.link(objectdata)
polyline = curvedata.splines.new('BEZIER')
polyline.bezier_points.add(len(bezier_list)-1)
for i, (anchor, h1, h2) in enumerate(bezier_list):
#print("bezier",,anchor, h1, h2)
point = polyline.bezier_points[i]
point.co = anchor
point.handle_left = h1
point.handle_right = h2
def get_points():
points = []
pencil = bpy.data.grease_pencil[-1]
for i, stroke in enumerate(pencil.layers[-1].active_frame.strokes):
for p in stroke.points:
#print(p.co.x, p.co.y, 0)
points.append(Vector((p.co.x, p.co.y, p.co.z)))
return points
def B0(t):
return (1.0-t) * (1.0-t) * (1.0-t)
def B1(t):
return 3 * t * (1.0-t) * (1.0-t)
def B2(t):
return 3 * t * t * (1.0-t)
def B3(t):
return t * t * t
def generate_bezier(points, U, tanL, tanR):
# calculate A
A = []
for u in U:
print(u)
A.append([tanL * B1(u), tanR * B2(u)])
# calculate C, X
C0 = 0.0
C1 = 0.0
C2 = 0.0
C3 = 0.0
X0 = 0.0
X1 = 0.0
for i, u in enumerate(U):
C0 = C0 + A[i][0].dot(A[i][0])
C1 = C1 + A[i][0].dot(A[i][1])
C2 = C1
C3 = C3 + A[i][1].dot(A[i][1])
PB0 = points[0] * B0(u)
PB1 = points[0] * B1(u)
PB2 = points[len(points)-1] * B2(u)
PB3 = points[len(points)-1] * B3(u)
print(PB0, PB1, PB2, PB3)
P = points[i] - (PB0 + PB1 + PB2 + PB3)
X0 = X0 + A[i][0].dot(P)
X1 = X1 + A[i][1].dot(P)
det_C0_C1 = C0 * C3 - C2 * C1
det_C0_X = C0 * X1 - C2 * X0
det_X_C1 = X0 * C3 - X1 * C1
alpha_l = det_X_C1 / det_C0_C1
alpha_r = det_C0_X / det_C0_C1;
if det_C0_C1 == 0:
alpha_l = 0
alpha_r = 0
p0 = points[0]
p3 = points[len(points)-1]
return [(p0, p0, p0 + tanL * alpha_l), (p3, p3 + tanR * alpha_r, p3)]
def create_relative_length_list(points):
u = [0]
for i, p in enumerate(points[1:]):
#print(i, u[i], p, points[i])
u.append( u[i] + (p - points[i]).length )
for i, v in enumerate(u):
u[i] = v / u[len(u)-1]
return u
def fit_curve(points, is_gpen):
if len(points) < 5:
return
points[0] = points[1]
U = create_relative_length_list(points)
tanL = (points[1] - points[0]).normalized()
tanR = (points[len(points)-2] - points[len(points)-1]).normalized()
if is_gpen:
tanL = (points[2] - points[1]).normalized()
tanR = (points[len(points)-4] - points[len(points)-3]).normalized()
print("tanL", tanL)
print("tanR", tanR)
bezier_list = generate_bezier(points, U, tanL, tanR)
make_curve(bezier_list)
fit_curve(get_points(), True)
------------------------------結果
------------------------------
参考文献：
(1) 金谷 健一 (2005)「これなら分かる最適化数学―基礎原理から計算手法まで」(共立出版)
(2) Andrew S. Glassner (1990) 「Graphics Gems (Graphics Gems - IBM)」(Morgan Kaufmann)
(3) Philip J. Schneider 「fit_curve.c」 from "Graphics Gems", Academic Press, 1990
------------------------------
Blender Advent Calender 2014 25日は blug.jp さんで、「Blenderドリル・クリスマスVer.(解答付き)」です☆
