日常生活でもベジェ曲線(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/
目的関数
目的関数は、ユーザーが入力した点と、我々が推測する曲線との誤差です。
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}
では、この2つのパラメータはどのように求めればよいのでしょうか。
いわゆるニュートン法を利用すれば求めることができます。
\frac{\partial }{\partial \boldsymbol{\alpha}}S(\boldsymbol{\alpha}_{ini} + \Delta \boldsymbol{\alpha}) = 0
テイラー展開し
\frac{\partial }{\partial \boldsymbol{\alpha}}S(\boldsymbol{\alpha}_{ini}) + \frac{\partial}{\partial\boldsymbol{\alpha}}\frac{\partial}{\partial \boldsymbol{\alpha}} S(\boldsymbol{\alpha}_{ini}) \Delta \boldsymbol{\alpha} = 0
つまり
\Delta \boldsymbol{\alpha} = - \left( \frac{\partial}{\partial\boldsymbol{\alpha}}\frac{\partial}{\partial \boldsymbol{\alpha}} S(\boldsymbol{\alpha}_{ini}) \right)^{-1} \frac{\partial }{\partial \boldsymbol{\alpha}}S(\boldsymbol{\alpha}_{ini})
実は、このニュートン法を最初から実装したい場合もあるのですが、$u_i$をどのように決定するかが問題になります。実は、$u_i$自体もニュートン法で決めることができるのです。
f(u) = \Big[ d_i - Q(u_i) \Big] \cdot \frac{d}{du}Q(u_i) = 0
ここでも、初期位置$u_{ini}$が$u_i$になるために移動すべき量$\Delta 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$を求めることができます。ここで注意したいのは、この2番目のケースで登場するニュートン法による$u_{ini}$の求め方が厄介なことで、初期位置によっては結果が得られない場合があるます。本論文のアルゴリズムでは、Chord length parameterizationという手法を用いて$u_{ini}$をある程度絞り込み、ベジェ曲線をフィッティングして繰り返し計算し、曲線セグメントの位置を決定しています。