この記事は、ChatGPTをつかって壁打ちをしながらGraphics Gems1には、「AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED CURVES」というアルゴリズムを解説する記事を作っていきます。AIが喋る情報は信用しないでください。

せがわ 2022/12/04

Graphics Gems1のAN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED GURVES のアルゴリズムをOpenSiv3Dで実装するとともに解説する記事を書きたい

Graphics Gems1には、「AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED GURVES」という論文が収録されています。この論文では、手書きの曲線を自動的にフィッティングするアルゴリズムが紹介されています。

そこで、このアルゴリズムをOpenSiv3Dで実装する方法を紹介し、実装例を解説します。

アルゴリズムの概要
「AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED GURVES」で紹介されているアルゴリズムは、以下のような流れで曲線をフィッティングします。

手書きの曲線データを取得する。
取得したデータを多項式によって近似する。
近似する多項式の次数を決定する。
多項式を用いて曲線を描画する。
OpenSiv3Dでの実装例
OpenSiv3Dでこのアルゴリズムを実装するには、以下のようなコードを書きます。

// OpenSiv3Dでの手書き曲線フィッティングアルゴリズム

#include <Siv3D.hpp>
#include <Eigen/Dense>

// 手書き曲線のデータを格納する配列
Array<Vec2> curveData;

// 多項式の次数
int polyDegree = 3;

// 多項式の係数
Eigen::VectorXd polyCoeffs(polyDegree + 1);

// 手書き曲線を近似する多項式を計算する
void calcPoly()
{
    // データ数
    const int n = curveData.size();

    // Vandermonde行列を作成
    Eigen::MatrixXd A(n, polyDegree + 1);
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j <= polyDegree; ++j)
        {
            A(i, j) = std::pow(curveData[i].x, j);
        }
    }

    // 右辺ベクトル

せがわ 2022/12/04

ベジェ曲線をいじれるアプリを作成しました。


# include <Siv3D.hpp> // OpenSiv3D v0.6.5

class ControllablePoint
{
public:

	ControllablePoint(const Vec2& point)
		: m_point(point)
	{

	}

	void update()
	{
		if (Circle(m_point, 10.0).leftClicked())
		{
			isActive = true;
			mouseOffset = Cursor::PosF() - m_point;
		}

		if (isActive)
		{
			m_point = Cursor::PosF() + mouseOffset;
		}

		if (MouseL.up())
		{
			isActive = false;
		}
	}

	void draw() const
	{
		if (Circle(m_point, 10.0).mouseOver())
		{
			Circle(m_point, 10.0).draw(Palette::White).drawFrame(4.0, Palette::Black);
		}
		else
		{
			Circle(m_point, 10.0).draw(Palette::White).drawFrame(3.0, Palette::Black);
		}
		
	}

	Vec2 getPoint() const
	{
		return m_point;
	}

private:

	Vec2 m_point = Vec2::Zero();
	Vec2 mouseOffset = Vec2::Zero();
	bool isActive = false;

};

void Main()
{
	// 背景の色を設定 | Set background color
	Scene::SetBackground(ColorF{1.0 });

	ControllablePoint p0(Vec2{ 300, 400 });
	ControllablePoint p1(Vec2{ 400, 400 });
	ControllablePoint p2(Vec2{ 400, 100 });
	ControllablePoint p3(Vec2{ 500, 100 });

	while (System::Update())
	{
		// 3 次ベジェ曲線
		Bezier3{ p0.getPoint(), p1.getPoint(), p2.getPoint(), p3.getPoint() }
		.draw(4, Palette::Black);

		Line(p0.getPoint(), p1.getPoint()).draw(1.0, Palette::Black);
		Line(p1.getPoint(), p2.getPoint()).draw(1.0, Palette::Black);
		Line(p2.getPoint(), p3.getPoint()).draw(1.0, Palette::Black);

		p0.update();
		p1.update();
		p2.update();
		p3.update();

		p0.draw();
		p1.draw();
		p2.draw();
		p3.draw();
	}
}

せがわ 2022/12/09

ベジェ曲線と点

# include <Siv3D.hpp> // OpenSiv3D v0.6.5

class ControllablePoint
{
public:

	ControllablePoint(const Vec2& point)
		: m_point(point)
	{

	}

	virtual void update()
	{
		if (Circle(m_point, 10.0).leftClicked())
		{
			m_isActive = true;
			m_mouseOffset = Cursor::PosF() - m_point;
		}

		if (m_isActive)
		{
			m_point = Cursor::PosF() + m_mouseOffset;
		}

		if (MouseL.up())
		{
			m_isActive = false;
		}
	}

	virtual void draw() const
	{
		if (Circle(m_point, 10.0).mouseOver())
		{
			Circle(m_point, 10.0).draw(Palette::White).drawFrame(4.0, Palette::Black);
		}
		else
		{
			Circle(m_point, 10.0).draw(Palette::White).drawFrame(3.0, Palette::Black);
		}
		
	}

	Vec2 getPoint() const
	{
		return m_point;
	}

protected:

	Vec2 m_point = Vec2::Zero();
	Vec2 m_mouseOffset = Vec2::Zero();
	bool m_isActive = false;

};

class ControllablePointWithAngle : public ControllablePoint
{
public:

	ControllablePointWithAngle(const Vec2& point)
		: ControllablePoint(point)
	{

	}

	void update() override
	{
		ControllablePoint::update();


		if (Circle(m_point + Vec2::Right(m_controlRadius).rotated(m_angle), 10.0).leftClicked())
		{
			m_angleIsActive = true;
		}

		if (m_angleIsActive)
		{
			m_angle = (Cursor::PosF() - m_point).getAngle() - Math::HalfPi;
		}

		if (MouseL.up())
		{
			m_angleIsActive = false;
		}
	}

	void draw() const override
	{
		ControllablePoint::draw();

		Circle(m_point, m_controlRadius).drawFrame(3.0, ColorF(0.0,0.3));

		const auto angleControlPoint = m_point + Vec2::Right(m_controlRadius).rotated(m_angle);

		if (Circle(angleControlPoint, 10.0).mouseOver())
		{
			Circle(angleControlPoint, 10.0).drawFrame(4.0, ColorF(0.0, 0.3));
		}
		else
		{
			Circle(angleControlPoint, 10.0).drawFrame(3.0, ColorF(0.0, 0.3));
		}
	}

	double getAngle() const
	{
		return m_angle;
	}

private:

	double m_controlRadius = 50.0;
	double m_angle = 0.0;

	bool m_angleIsActive = false;

};

class ControllableBezierCurve
{
public:

	void update()
	{
		StartPoint.update();
		//p1.update();
		//p2.update();
		EndPoint.update();
	}

	void draw()
	{		// 3 次ベジェ曲線
		const auto bezier = getBezier();

		bezier.draw(4, Palette::Black);

		const auto p0 = StartPoint.getPoint();
		const auto p1 = StartPoint.getPoint() + Vec2::Right(t1).rotated(StartPoint.getAngle());
		const auto p2 = EndPoint.getPoint() + Vec2::Right(t2).rotated(EndPoint.getAngle());
		const auto p3 = EndPoint.getPoint();

		Line(p0, p1).draw(1.0, Palette::Black);
		Line(p1, p2).draw(1.0, Palette::Black);
		Line(p2, p3).draw(1.0, Palette::Black);

		StartPoint.draw();
		//p1.draw();
		//p2.draw();
		EndPoint.draw();

		SimpleGUI::Slider(t1, 0.0, 500.0, p0 + Vec2::Down(100.0));
		SimpleGUI::Slider(t2, 0.0, 500.0, p3 + Vec2::Down(100.0));


	}

	Bezier3 getBezier() const
	{
		return Bezier3{ StartPoint.getPoint()
			,StartPoint.getPoint() + Vec2::Right(t1).rotated(StartPoint.getAngle())
			,EndPoint.getPoint() + Vec2::Right(t2).rotated(EndPoint.getAngle())
			, EndPoint.getPoint() };
	}

private:

	ControllablePointWithAngle StartPoint{ Vec2{ 300, 400 } };
	ControllablePointWithAngle EndPoint{Vec2{ 500, 100 }};

	double t1 = 100.0;
	double t2 = 100.0;

};

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>();
			}
		}
	}

	return result;
}


void Main()
{
	// 背景の色を設定 | Set background color
	Scene::SetBackground(ColorF{1.0 });

	ControllableBezierCurve bezier;


	Array<Vec2> mouseVertices;
	Array<Array<Vec2>> devidedLines;




	while (System::Update())
	{
		if (MouseR.pressed())
		{

			if (mouseVertices.size() == 0)
			{
				mouseVertices << Cursor::PosF();

			}
			else if (mouseVertices.back().distanceFrom(Cursor::PosF()) > 5.0)
			{
				mouseVertices << Cursor::PosF();
			}
		}

		if (MouseR.up() && mouseVertices.size() > 5)
		{
			devidedLines = SplitVertices(mouseVertices);

			mouseVertices.clear();
		}


		if (KeyC.down())
		{
			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 < devidedLines.size(); ++i)
		{
			for (size_t k = 0; k < devidedLines[i].size(); ++k)
			{
				Circle(devidedLines[i][k], 10.0).draw(HSV(i * 200));
			}

		}
		//---------------------------------------------------------

		//ベジェ曲線の操作と干渉させたくない。
		bezier.update();
		bezier.draw();
	}
}

せがわ 2022/12/09

素晴らしい動画
- https://youtu.be/jvPPXbo87ds
- https://youtu.be/aVwxzDHniEw

せがわ 2022/12/09

点をもらい、それをBスプラインで表現する。始点と終点とその方向はある程度点から予測できる。
G1連続性　あとは　パラメータt1,t2をどのように調節するか？の問題に帰着できる。
ニュートン法でやると精度もよくいけるのだが、初期値によって収束しない問題がある。
そのために、ある程度解に近い、良い初期値に単純な計算で誘導したあとに、ニュートン法を実行する。

せがわ 2023/01/23

記事を書きました

せがわ 2023/01/23

補足

このアルゴリズムの重要な点は以下の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}$ をある程度絞り込み、ベジェ曲線をフィッティングして繰り返し計算し、曲線セグメントの位置を決定しています。

このスクラップは2023/01/23にクローズされました