PolygonをCircleとして近似する関数

[sknjpn]

物理エンジンに頼らずに慣性モーメントが扱いたいです。 下は自分が作成したものです。

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

// Polygonが一様な面密度の板と仮定する
double	getInertia(const Polygon& polygon, double mass)
{
	const auto centroid = polygon.centroid();
	const auto numTriangles = polygon.num_triangles();
	const auto area = polygon.area();

	double inertia = 0;
	for (int i = 0; i < numTriangles; ++i)
	{
		const auto triangle = polygon.triangle(i);
		const auto v0 = triangle.p0 - triangle.centroid();
		const auto v1 = triangle.p1 - triangle.centroid();
		const auto v2 = triangle.p2 - triangle.centroid();
		const auto m = mass * triangle.area() / area;
		const auto it = m / 12 * (v0.lengthSq() + v1.lengthSq() + v2.lengthSq());
		const auto id = m * (triangle.centroid() - centroid).lengthSq();

		inertia += it + id;
	}

	return inertia;
}

// Polygonを一様な面密度の板と仮定し、円盤として近似を行う
Circle	approximateAsCircle(const Polygon& polygon)
{
	return Circle(polygon.centroid(), Sqrt(2 * getInertia(polygon, 1.0)));
}

void Main()
{
	Polygon polygon;

	while (System::Update())
	{
		auto c = Circle(Cursor::PosF(), 16.0);

		if (MouseL.pressed())
		{
			polygon.append(c.asPolygon());
		}
		if (MouseR.pressed())
		{
			auto s = Geometry2D::Subtract(polygon, c.asPolygon());

			if (s.empty()) polygon = Polygon();
			else
			{
				int maxIndex = 0;
				for (int i = 1; i < s.size(); ++i)
					if (s[maxIndex].area() < s[i].area()) maxIndex = i;

				polygon = s[maxIndex];
			}
		}

		c.draw(Palette::White);

		polygon.draw(ColorF(Palette::Red, 0.5));
		approximateAsCircle(polygon).draw(ColorF(1.0, 0.25)).drawFrame(2.0, Palette::White);
	}
}