日常生活でもベジェ曲線(Bézier curve)
みなさん、ベジェ曲線はご存知ですか？
ベクターグラフィックソフトやペイントツールを普段からお使いの人は、使ったことがあるかもしれません。
このベジェ曲線を使うと何が嬉しいかというと、少ない制御点で美しい曲線が描けるようになることです。例えば、上のペイントツールでは、4つの位置を指定するだけで、これだけ滑らかな曲線が描けるのです。OpenSiv3Dにもベジェ曲線用の機能もあります。
でも、ひとつ疑問がありました。手書きの自由曲線を描いたとして、この曲線をベジェ曲線に変換するにはどうしたらいいのでしょうか？
頑張って点を配置して、それらしい素敵な曲線を手作業で探すしかないのでしょうか？(実は、パワーポイントのスライドを作るときに、手動で曲線に変換していました)
これを調べてみると、この作業を自動化する方法がGraphicsGemsの AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED CURVESで提案されていることがわかりました。
本の論文
https://dl.acm.org/doi/10.5555/90767.90941
Graphics gemsのPDFはこちらのリンクで見れました
https://theswissbay.ch/pdf/Gentoomen%20Library/Game%20Development/Programming/Graphics%20Gems%201.pdf
オリジナルのCのソースコードもあるよ
https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c
http://www.graphicsgems.org/
この論文では、デジタル化された曲線を自動的にフィッティングするためのアルゴリズムが提案されています。今回は、提案したアルゴリズムに基づき、OpenSiv3Dを用いて手描きの点からベジェ曲線を生成するソフトを作成しました。
AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED CURVESの実装
上の動画でもわかるように、ベジェ曲線は手描きの曲線から自動処理で生成されています。この記事で詳しい説明を書こうと思ったのですが、すみません、間に合わなかったので、OpenSiv3Dの開発途中の実装を紹介します。
サンプルコード
Main.cpp
# include <Siv3D.hpp> // OpenSiv3D v0.6.5
double BernsteinPolynomials(int32 i, double u)
{
if (i == 0)
{
double tmp = 1.0 - u;
return (tmp * tmp * tmp);
}
else if (i == 1)
{
double tmp = 1.0 - u;
return (3 * u * (tmp * tmp));
}
else if (i == 2)
{
double tmp = 1.0 - u;
return (3 * u * u * tmp);
}
else if (i == 3)
{
return (u * u * u);
}
else
{
return 0.0;
}
}
class ControllablePoint
{
public:
ControllablePoint(const Vec2& point)
: point(point)
{
}
virtual void update()
{
if (Circle(point, m_radius).leftClicked())
{
m_isActive = true;
m_mouseOffset = Cursor::PosF() - point;
}
if (m_isActive)
{
point = Cursor::PosF() + m_mouseOffset;
}
if (MouseL.up())
{
m_isActive = false;
}
}
virtual void draw() const
{
if (Circle(point, m_radius).mouseOver())
{
Circle(point, m_radius).draw(Palette::White).drawFrame(4.0, Palette::Black);
}
else
{
Circle(point, m_radius).draw(Palette::White).drawFrame(3.0, Palette::Black);
}
}
Vec2 point = Vec2::Zero();
protected:
Vec2 m_mouseOffset = Vec2::Zero();
bool m_isActive = false;
double m_radius = 5.0;
};
class ControllablePointWithAngle : public ControllablePoint
{
public:
ControllablePointWithAngle(const Vec2& point)
: ControllablePoint(point)
{
}
void update() override
{
ControllablePoint::update();
if (Circle(point + Vec2::Right(m_controlRadius).rotated(angle), m_radius).leftClicked())
{
m_angleIsActive = true;
}
if (m_angleIsActive)
{
angle = (Cursor::PosF() - point).getAngle() - Math::HalfPi;
}
if (MouseL.up())
{
m_angleIsActive = false;
}
}
void draw() const override
{
ControllablePoint::draw();
//Circle(point, m_controlRadius).drawFrame(3.0, ColorF(0.0, 0.3));
const auto angleControlPoint = point + Vec2::Right(m_controlRadius).rotated(angle);
//if (Circle(angleControlPoint, m_radius).mouseOver())
//{
// Circle(angleControlPoint, m_radius).drawFrame(4.0, ColorF(0.0, 0.3));
//}
//else
//{
// Circle(angleControlPoint, m_radius).drawFrame(3.0, ColorF(0.0, 0.3));
//}
}
double angle = 0.0;
private:
double m_controlRadius = 50.0;
bool m_angleIsActive = false;
};
class ControllableBezierCurve
{
public:
void update()
{
startControllablePoint.update();
endControllablePoint.update();
}
void draw()
{ // 3 次ベジェ曲線
const auto bezier = getBezier();
bezier.draw(4, Palette::Black);
const auto p0 = startControllablePoint.point;
const auto p1 = startControllablePoint.point + Vec2::Right(t1).rotated(startControllablePoint.angle);
const auto p2 = endControllablePoint.point + Vec2::Right(t2).rotated(endControllablePoint.angle);
const auto p3 = endControllablePoint.point;
Line(p0, p1).draw(1.0, Palette::Black);
Line(p1, p2).draw(1.0, Palette::Black);
Line(p2, p3).draw(1.0, Palette::Black);
startControllablePoint.draw();
//p1.draw();
//p2.draw();
endControllablePoint.draw();
//SimpleGUI::Slider(t1, 0.0, 1000.0, p0 + Vec2(100.0, 100.0));
//SimpleGUI::Slider(t2, 0.0, 1000.0, p3 + Vec2(-100.0, 100.0));
}
Bezier3 getBezier() const
{
return Bezier3{ getV0(), getV1(), getV2(), getV3() };
}
Vec2 getV0() const
{
return startControllablePoint.point;
}
Vec2 getV1() const
{
return startControllablePoint.point + Vec2::Right(t1).rotated(startControllablePoint.angle);
}
Vec2 getV2() const
{
return endControllablePoint.point + Vec2::Right(t2).rotated(endControllablePoint.angle);
}
Vec2 getV3() const
{
return endControllablePoint.point;
}
ControllablePointWithAngle startControllablePoint{ Vec2{ 300, 400 } };
ControllablePointWithAngle endControllablePoint{ Vec2{ 500, 100 } };
double t1 = 100.0;
double t2 = 100.0;
};
double GetSumOfVerticesLength(const Array<Vec2>& vertices)
{
double sumOfLength = 0.0;
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
sumOfLength += (vertices[i + 1] - vertices[i]).length();
}
return sumOfLength;
}
struct Segment
{
Array<Vec2> vertices;
ControllableBezierCurve bezierCurve;
Array<double> parameters;
size_t maxIteration = 4;
void drawParameterPos() const
{
for (size_t i = 0; i < parameters.size(); ++i)
{
Circle(bezierCurve.getBezier().getPos(parameters[i]), 5.0).draw(ColorF(0.8));
}
}
Array<double> CalculateChordLengthParameters() const
{
Array<double> result;
double tLength = 0.0;
const double sumOfLength = GetSumOfVerticesLength();
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
tLength += (vertices[i + 1] - vertices[i]).length();
const double ratio = tLength / sumOfLength;
result << ratio;
}
return result;
}
Array<double> CalculateNewtonParameters() const
{
Array<double> result;
const auto cubicBezier = bezierCurve.getBezier();
const Vec2 bezier2V0 = (bezierCurve.getV1() - bezierCurve.getV0()) * 3.0;
const Vec2 bezier2V1 = (bezierCurve.getV2() - bezierCurve.getV1()) * 3.0;
const Vec2 bezier2V2 = (bezierCurve.getV3() - bezierCurve.getV2()) * 3.0;
const auto bezier2 = Bezier2(bezier2V0, bezier2V1, bezier2V2);
const Vec2 bezier1V0 = (bezier2V1 - bezier2V0) * 2.0;
const Vec2 bezier1V1 = (bezier2V2 - bezier2V1) * 2.0;
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
double parameter = parameters[i];
for (size_t k = 0; k < maxIteration; ++k)
{
const Vec2 cubicBezierPos = cubicBezier.getPos(parameter);
const Vec2 bezier2Pos = bezier2.getPos(parameter);
const Vec2 bezier1Pos = bezier1V0 + (bezier1V1 - bezier1V0) * parameter;
const double numerator = (cubicBezierPos - vertices[i + 1]).dot(bezier2Pos);
const double denominator = bezier2Pos.dot(bezier2Pos) + (cubicBezierPos - vertices[i + 1]).dot(bezier1Pos);
if (denominator == 0.0) break;
parameter += -(numerator / denominator);
//力尽きました
//https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c
//上記のNewtonRaphsonRootFindを参考に実装してください
}
result << parameter;
}
return result;
}
double GetSumOfVerticesLength() const
{
double sumOfLength = 0.0;
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
sumOfLength += (vertices[i + 1] - vertices[i]).length();
}
return sumOfLength;
}
Vec2 FunctionA(const Vec2 direction, double ratio, int32 i)
{
return direction * BernsteinPolynomials(i, ratio);
}
double CalculateC11()
{
double result = 0.0;
const Vec2 startDirection = Vec2::Right(1.0).rotated(bezierCurve.startControllablePoint.angle);
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
result += FunctionA(startDirection, parameters[i], 1).lengthSq();
}
return result;
}
double CalculateC12()
{
double result = 0.0;
const Vec2 startDirection = Vec2::Right(1.0).rotated(bezierCurve.startControllablePoint.angle);
const Vec2 endDirection = Vec2::Right(1.0).rotated(bezierCurve.endControllablePoint.angle);
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
result += FunctionA(startDirection, parameters[i], 1).dot(FunctionA(endDirection, parameters[i], 2));
}
return result;
}
double CalculateC21()
{
return CalculateC12();
}
double CalculateC22()
{
double result = 0.0;
const Vec2 endDirection = Vec2::Right(1.0).rotated(bezierCurve.endControllablePoint.angle);
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
result += FunctionA(endDirection, parameters[i], 2).lengthSq();
}
return result;
}
double CalculateX1()
{
double result = 0.0;
const Vec2 startDirection = Vec2::Right(1.0).rotated(bezierCurve.startControllablePoint.angle);
const Vec2& startPoint = bezierCurve.startControllablePoint.point;
const Vec2& endPoint = bezierCurve.endControllablePoint.point;
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
const Vec2 v = vertices[i + 1] -
(startPoint * BernsteinPolynomials(0, parameters[i])
+ startPoint * BernsteinPolynomials(1, parameters[i])
+ endPoint * BernsteinPolynomials(2, parameters[i])
+ endPoint * BernsteinPolynomials(3, parameters[i]));
result += v.dot(FunctionA(startDirection, parameters[i], 1));
}
return result;
}
double CalculateX2()
{
double result = 0.0;
const Vec2 endDirection = Vec2::Right(1.0).rotated(bezierCurve.endControllablePoint.angle);
const Vec2& startPoint = bezierCurve.startControllablePoint.point;
const Vec2& endPoint = bezierCurve.endControllablePoint.point;
for (size_t i = 0; i < vertices.size() - 1; ++i)
{
const Vec2 v = vertices[i + 1] -
(startPoint * BernsteinPolynomials(0, parameters[i])
+ startPoint * BernsteinPolynomials(1, parameters[i])
+ endPoint * BernsteinPolynomials(2, parameters[i])
+ endPoint * BernsteinPolynomials(3, parameters[i]));
result += v.dot(FunctionA(endDirection, parameters[i], 2));
}
return result;
}
double Det(double v11, double v12, double v21, double v22)
{
return v11 * v22 - v12 * v21;
}
void setRoughApproximation()
{
parameters = CalculateChordLengthParameters();
const double c11 = CalculateC11();
const double c12 = CalculateC12();
const double c21 = CalculateC21();
const double c22 = CalculateC22();
const double x1 = CalculateX1();
const double x2 = CalculateX2();
bezierCurve.t1 = Det(x1, c12, x2, c22) / Det(c11, c12, c21, c22);
bezierCurve.t2 = Det(c11, x1, c21, x2) / Det(c11, c12, c21, c22);
}
void setExactApproximation()
{
parameters = CalculateNewtonParameters();
const double c11 = CalculateC11();
const double c12 = CalculateC12();
const double c21 = CalculateC21();
const double c22 = CalculateC22();
const double x1 = CalculateX1();
const double x2 = CalculateX2();
bezierCurve.t1 = Det(x1, c12, x2, c22) / Det(c11, c12, c21, c22);
bezierCurve.t2 = Det(c11, x1, c21, x2) / Det(c11, c12, c21, c22);
}
};
Array<Array<Vec2>> SplitVertices(const Array<Vec2>& vertices)
{
Array<Array<Vec2>> result;
size_t dividedIndex = 0;
result << Array<Vec2>();
for (size_t i = 0; i < vertices.size(); ++i)
{
const auto& v = vertices[i];
result[dividedIndex] << v;
if (i > 2 && i < vertices.size() - 1 - 2)
{
const Vec2& a1 = (v - vertices[i - 1]).normalized();
const Vec2& a2 = (vertices[i + 1] - v).normalized();
const double angle = a1.getAngle(a2);
if (Abs(angle) > 1.0)
{
++dividedIndex;
result << Array<Vec2>();
result[dividedIndex] << v;
}
}
}
return result;
}
double SimpleCurveFitEvaluationFunction(const Segment& curve)
{
double result = 0.0;
double tLength = 0.0;
if (curve.vertices.size() == 0)
{
return 0.0;
}
const double sumOfLength = curve.GetSumOfVerticesLength();
for (size_t i = 0; i < curve.vertices.size() - 1; ++i)
{
tLength += (curve.vertices[i + 1] - curve.vertices[i]).length();
const double ratio = tLength / sumOfLength;
const Vec2 pos = curve.bezierCurve.getBezier().getPos(ratio);
result += (pos - curve.vertices[i + 1]).lengthSq();
Circle(pos, 10.0).draw(Palette::Red);
Line(pos, curve.vertices[i + 1]).draw(1.0, Palette::Red);
}
return result;
}
Array<Segment> GenerateSegmentsFromDividedVertices(const Array<Array<Vec2>>& devidedVertices)
{
Array<Segment> segments;
for (size_t i = 0; i < devidedVertices.size(); ++i)
{
auto& vertices = devidedVertices[i];
if (vertices.size() < 6)
{
ControllableBezierCurve bezier;
bezier.startControllablePoint.point = vertices[0];
bezier.endControllablePoint.point = vertices[vertices.size() - 1];
bezier.startControllablePoint.angle = (vertices[vertices.size() - 1] - vertices[0]).getAngle() - Math::HalfPi;
bezier.endControllablePoint.angle = (vertices[0] - vertices[vertices.size() - 1]).getAngle() - Math::HalfPi;
bezier.t1 = bezier.startControllablePoint.point.distanceFrom(bezier.endControllablePoint.point) / 4.0;
bezier.t2 = bezier.t1;
Segment segment;
segment.vertices = vertices;
segment.bezierCurve = bezier;
segment.setRoughApproximation();
segments << segment;
continue;
}
ControllableBezierCurve bezier;
bezier.startControllablePoint.point = vertices[0];
bezier.endControllablePoint.point = vertices[vertices.size() - 1];
bezier.startControllablePoint.angle = (vertices[1] + vertices[2] - 2.0 * vertices[0]).getAngle() - Math::HalfPi;
bezier.endControllablePoint.angle = (vertices[vertices.size() - 2] + vertices[vertices.size() - 3] - 2.0 * vertices[vertices.size() - 1]).getAngle() - Math::HalfPi;
Segment segment;
segment.vertices = vertices;
segment.bezierCurve = bezier;
segment.setRoughApproximation();
segments << segment;
}
return segments;
}
Array<Segment> GenerateSegmentsFromG1ContinuousDividedVertices(const Array<Array<Vec2>>& devidedVertices)
{
Array<Segment> segments;
if (devidedVertices.size() <= 1)
{
return segments;
}
for (size_t i = 0; i < devidedVertices.size(); ++i)
{
auto& vertices = devidedVertices[i];
if (vertices.size() < 4)
{
break;
}
ControllableBezierCurve bezier;
bezier.startControllablePoint.point = vertices[0];
bezier.endControllablePoint.point = vertices[vertices.size() - 1];
bezier.startControllablePoint.angle = Math::Fmod((vertices[1] + vertices[2] - 2.0 * vertices[0]).getAngle() - Math::HalfPi + Math::TwoPi, Math::TwoPi);
bezier.endControllablePoint.angle = Math::Fmod((vertices[vertices.size() - 2] + vertices[vertices.size() - 3] - 2.0 * vertices[vertices.size() - 1]).getAngle() - Math::HalfPi + Math::TwoPi, Math::TwoPi);
Segment segment;
segment.vertices = vertices;
segment.bezierCurve = bezier;
segment.setRoughApproximation();
segments << segment;
}
for (size_t i = 1; i < segments.size(); ++i)
{
const double average = (segments[i - 1].bezierCurve.endControllablePoint.angle + Math::Fmod(segments[i].bezierCurve.startControllablePoint.angle + Math::Pi + Math::TwoPi, Math::TwoPi)) / 2.0;
segments[i - 1].bezierCurve.endControllablePoint.angle = average;
segments[i].bezierCurve.startControllablePoint.angle = Math::Fmod(average + Math::Pi, Math::TwoPi);
}
return segments;
}
Array<Segment> SplitCurve(const Segment& segment)
{
Array<double> errors;
Array<Array<Vec2>> devidedVeritces;
size_t dividedIndex = 0;
devidedVeritces << Array<Vec2>();
Array<Segment> segments;
double tLength = 0.0;
if (segment.vertices.size() < 13)
{
segments << segment;
return segments;
}
const double sumOfLength = segment.GetSumOfVerticesLength();
for (size_t i = 0; i < segment.vertices.size() - 1; ++i)
{
tLength += (segment.vertices[i + 1] - segment.vertices[i]).length();
const double ratio = tLength / sumOfLength;
const Vec2 pos = segment.bezierCurve.getBezier().getPos(ratio);
const double error = (i < 3 || i > segment.vertices.size() - 4) ? 0.0 : (pos - segment.vertices[i + 1]).length();
errors << error;
}
//最も大きいエラーのインデックスを返す
auto iter = std::max_element(errors.begin(), errors.end());
size_t index = std::distance(errors.begin(), iter);
if (*iter < 10)
{
segments << segment;
return segments;
}
for (size_t i = 0; i < segment.vertices.size(); ++i)
{
devidedVeritces[dividedIndex] << segment.vertices[i];
if (i == index)
{
++dividedIndex;
devidedVeritces << Array<Vec2>();
devidedVeritces[dividedIndex] << segment.vertices[i];
}
}
segments = GenerateSegmentsFromG1ContinuousDividedVertices(devidedVeritces);
for (auto& segment : segments)
{
segment.setRoughApproximation();
}
return segments;
}
struct Curve
{
Curve(const Array<Vec2>& mouseVertices)
{
const auto devidedVertices = SplitVertices(mouseVertices);
segments = GenerateSegmentsFromDividedVertices(devidedVertices);
}
Array<Segment> segments;
};
void Main()
{
// 背景の色を設定 | Set background color
Scene::SetBackground(ColorF{ 1.0 });
Array<Vec2> mouseVertices;
Array<Curve> curves;
while (System::Update())
{
if (SimpleGUI::Button(U"(1) Split", Vec2{ 10, 10 }, 190))
{
for (auto& curve : curves)
{
Array<Segment> newSegments;
for (auto& segment : curve.segments)
{
Array<Segment> splittedCurves = SplitCurve(segment);
for (auto& c : splittedCurves)
{
newSegments << c;
}
}
curve.segments = newSegments;
}
}
if (SimpleGUI::Button(U"(2) Reparameterize", Vec2{ 210, 10 }, 190))
{
for (auto& curve : curves)
{
for (auto& segment : curve.segments)
{
segment.setExactApproximation();
}
}
}
if (SimpleGUI::Button(U"(3) Clear", Vec2{ 410, 10 }, 190))
{
curves.clear();
mouseVertices.clear();
}
if (Cursor::Pos().y > 60 && MouseL.pressed())
{
if (mouseVertices.size() == 0)
{
mouseVertices << Cursor::PosF();
}
else if (mouseVertices.back().distanceFrom(Cursor::PosF()) > 5.0)
{
mouseVertices << Cursor::PosF();
}
}
if (MouseL.up() && mouseVertices.size() > 5)
{
curves << Curve(mouseVertices);
mouseVertices.clear();
}
for (size_t i = 0; i < mouseVertices.size(); ++i)
{
const auto& v = mouseVertices[i];
Circle(v, 3.0).draw(Palette::Black);
}
for (size_t i = 0; i < curves.size(); ++i)
{
for (size_t j = 0; j < curves[i].segments.size(); ++j)
{
for (size_t k = 0; k < curves[i].segments[j].vertices.size(); ++k)
{
Circle(curves[i].segments[j].vertices[k], 10.0).draw(HSV(i * 200, 0.1, 1.0));
}
}
}
for (auto& curve : curves)
{
for (auto& segment : curve.segments)
{
segment.drawParameterPos();
segment.bezierCurve.draw();
}
}
}
}
補足
このアルゴリズムの重要な点は以下の2点です。
- 複雑な線をどのようにベジェ曲線で表現しやすい領域に分割しているのか？
- その分割した領域の曲線に、どのようにベジェ曲線を当てはめているのか？
これを踏まえて、以下のような前提条件を用意しておくと、よりスムーズに読み進めることができるかもしれません。
ベジェ曲線
3次ベジェ曲線は、2次元平面上に描かれる曲線の一種である。この曲線は4つのデータ点を使って表現されます。
Q(t) = \sum_{i=0}^n V_i B^n_i(t), \quad t \in [0, 1], \\
B^n_i(t) = \frac{n!}{(n-i)!i!} t^i (1-t)^{n-i}, \quad i = 0,...,n.
面白い動画もあるので見てね 傑作動画
https://youtu.be/aVwxzDHniEw
https://youtu.be/jvPPXbo87ds
インタラクティブなデモもあるよ
https://pomax.github.io/bezierinfo/
目的関数
ここでは原論文の導出に沿って説明します。実装上の見方や別解釈は Appendix を参照してください。
目的関数は、ユーザーが入力した点と、我々が推測する曲線との誤差です。
S = \sum_{i=1}^n \Big[ d_i - Q(u_i) \Big]^2
1本のBezier曲線に対して、求めるべき未知数は端点の単位接線方向に沿った内側2制御点までの距離 $\alpha_1, \alpha_2$ だけです（G1 連続）。
\boldsymbol{\alpha} = (\alpha_1, \alpha_2)^\mathrm{T}
これらは目的関数 $S$ を $\alpha_1,\alpha_2$ で偏微分して 0 とおくと得られる 2×2 の線形方程式系を解くことで解析的に決まります。
\frac{\partial S}{\partial \alpha_1}=0,\quad \frac{\partial S}{\partial \alpha_2}=0
端点の単位接線を $\hat{\boldsymbol t}_1,\hat{\boldsymbol t}_2$、バーンスタイン多項式を $B_i^3(u)$ とし、
A_{i,1}=\hat{\boldsymbol t}_1\,B_1^3(u_i),\quad A_{i,2}=\hat{\boldsymbol t}_2\,B_2^3(u_i)
と定義すると、
\begin{bmatrix}
c_{11} & c_{12} \\
c_{21} & c_{22}
\end{bmatrix}
\begin{bmatrix}\alpha_1\\ \alpha_2\end{bmatrix}
=
\begin{bmatrix}X_1\\ X_2\end{bmatrix}
が得られます。ここで
c_{11}=\sum_i A_{i,1}\cdot A_{i,1},\quad
c_{12}=c_{21}=\sum_i A_{i,1}\cdot A_{i,2},\quad
c_{22}=\sum_i A_{i,2}\cdot A_{i,2},
X_1=\sum_i\left[d_i - V_0 B_0^3(u_i) - V_0 B_1^3(u_i) - V_3 B_2^3(u_i) - V_3 B_3^3(u_i)\right]\cdot A_{i,1},
X_2=\sum_i\left[d_i - V_0 B_0^3(u_i) - V_0 B_1^3(u_i) - V_3 B_2^3(u_i) - V_3 B_3^3(u_i)\right]\cdot A_{i,2}.
よって例えばクラメルの公式で
\alpha_1=\frac{\det\begin{bmatrix}X_1&c_{12}\\ X_2&c_{22}\end{bmatrix}}{\det\begin{bmatrix}c_{11}&c_{12}\\ c_{21}&c_{22}\end{bmatrix}},\quad
\alpha_2=\frac{\det\begin{bmatrix}c_{11}&X_1\\ c_{21}&X_2\end{bmatrix}}{\det\begin{bmatrix}c_{11}&c_{12}\\ c_{21}&c_{22}\end{bmatrix}}.
上の実装の CalculateC11〜CalculateC22 と CalculateX1/CalculateX2 はそれぞれ H と X の計算に対応し、t1/t2 が (\alpha_1/\alpha_2) に対応します。
一方で、各点に対応づけるパラメータ $u_i$ は、初期割当（Chord length parameterization）だけでは精度が不十分な場合があるため、ニュートン–ラフソン法で再パラメータ化して改善します。
f(u) = \big[ Q(u) - d_i \big] \cdot \frac{d}{du}Q(u) = 0
初期値 $u_{ini}$ を $u$ に近づける更新は、
f(u_{ini} + \Delta u) = 0
テイラー展開より
f(u_{ini}) + \frac{d}{du}f(u_{ini}) \Delta u = 0
すなわち
\Delta u = - \frac{f(u_{ini})}{\frac{d}{du}f(u_{ini})}
この反復で各 $u_i$ を改善します。重要なのは、ニュートン法は $u_i$ の再推定にのみ用い、$\alpha_1,\alpha_2$ は上記の線形方程式を解いて求める、という役割分担です。初期値には chord length parameterization を用い、フィッティングと再パラメータ化を数回繰り返して収束させます。
Appendix: αをベクトルとして見る
- 固定した $u_i$ のもとでは $Q$ は $\alpha$ に線形、したがって $S(\alpha)$ は二次です。最小二乗の正規方程式 $H\alpha=X$（本文の $c_{ij}$ と $X_i$）で $\alpha$ が決まります。
- 行列表現では $\alpha = H^{-1}X$ と書けます。2×2 なら本文の式（クラメル）でも、任意の安定な線形ソルバ（ガウス/QR/SVD 等）でも同じ解に到達します。
- なお $H$ が一定のため、ニュートン更新でも1回で同じ解になります。実装方針は用途に応じて選べば十分です。