ch19-smallest-enclosing-circle: 平面上の点群に対し最小包含円を求めるアルゴリズムの実装 (#1272)

Siv3D · Nov 18, 2024 · 458bec9 · 458bec9
1 parent f28b27f
commit 458bec9
Show file tree

Hide file tree

Showing 3 changed files with 225 additions and 3 deletions.
diff --git a/Siv3D/include/Siv3D/Geometry2D.hpp b/Siv3D/include/Siv3D/Geometry2D.hpp
@@ -2762,6 +2762,56 @@ namespace s3d
 		[[nodiscard]]
 		Polygon ConcaveHull(const Vec2* points, size_t size, double concavity = 2.0, double lengthThreshold = 0.0);
 
+		//////////////////////////////////////////////////
+		//
+		//	SmallestEnclosingCircle
+		//
+		//////////////////////////////////////////////////
+
+		/// @brief 点 p0, p1, p2 の最小包含円を返します。
+		/// @param p0 点 0
+		/// @param p1 点 1
+		/// @param p2 点 2
+		/// @return 点 p0, p1, p2 の最小包含円
+		[[nodiscard]]
+		Circle SmallestEnclosingCircle(const Vec2& p0, const Vec2& p1, const Vec2& p2);
+
+		/// @brief 点 p0, p1, p2, p3 の最小包含円を返します。
+		/// @param p0 点 0
+		/// @param p1 点 1
+		/// @param p2 点 2
+		/// @param p3 点 3
+		/// @param tolerance 点が円に含まれているかの判定時の許容誤差。相対誤差または絶対誤差がこの値以下であれば、点が円に含まれているとみなします。
+		/// @return 点 p0, p1, p2, p3 の最小包含円
+		[[nodiscard]]
+		Circle SmallestEnclosingCircle(const Vec2& p0, const Vec2& p1, const Vec2& p2, const Vec2& p3, double tolerance = 1e-8);
+
+		/// @brief 点群 points の最小包含円を返します。
+		/// @param points 点群
+		/// @param tolerance 点が円に含まれているかの判定時の許容誤差。相対誤差または絶対誤差がこの値以下であれば、点が円に含まれているとみなします。
+		/// @return 点群 points の最小包含円
+		[[nodiscard]]
+		Circle SmallestEnclosingCircle(Array<Vec2> points, double tolerance = 1e-8);
+
+		/// @brief 点群 points の最小包含円を返します。
+		/// @tparam URBG 乱数生成器の型
+		/// @param points 点群
+		/// @param tolerance 点が円に含まれているかの判定時の許容誤差。相対誤差または絶対誤差がこの値以下であれば、点が円に含まれているとみなします。
+		/// @param urbg 乱数生成器。このアルゴリズムには点群の順序をシャッフルする処理が含まれていて、乱数生成器を使用します。
+		/// @return 点群 points の最小包含円
+		SIV3D_CONCEPT_URBG
+		[[nodiscard]]
+		Circle SmallestEnclosingCircle(Array<Vec2> points, double tolerance, URBG&& urbg);
+
+		/// @brief 点群 points の最小包含円を返します。
+		/// @tparam URBG 乱数生成器の型
+		/// @param points 点群
+		/// @param urbg 乱数生成器。このアルゴリズムには点群の順序をシャッフルする処理が含まれていて、乱数生成器を使用します。
+		/// @param tolerance 点が円に含まれているかの判定時の許容誤差。相対誤差または絶対誤差がこの値以下であれば、点が円に含まれているとみなします。
+		SIV3D_CONCEPT_URBG
+		[[nodiscard]]
+		Circle SmallestEnclosingCircle(Array<Vec2> points, URBG&& urbg, double tolerance = 1e-8);
+
 		//////////////////////////////////////////////////
 		//
 		//	Subtract
@@ -2869,4 +2919,4 @@ namespace s3d
 	}
 }
 
-# include "detail/Geometry2D.ipp"
+# include "detail/Geometry2D.ipp"
diff --git a/Siv3D/include/Siv3D/detail/Geometry2D.ipp b/Siv3D/include/Siv3D/detail/Geometry2D.ipp
@@ -41,6 +41,26 @@ namespace s3d
 				&& (point.y < bottom);
 		}
 
+		/// @brief 点 p が円 c に含まれているかを判定します。
+		/// @param c 円
+		/// @param p 点
+		/// @param tolerance 許容誤差（相対誤差または絶対誤差のいずれかが許容誤差以下であれば許容）
+		/// @return 点 p が円 c に含まれている場合 true, それ以外の場合は false
+		[[nodiscard]]
+		inline bool Contains(const Circle& c, const Vec2& p, const double tolerance = 1e-8)
+		{
+			const double dSquared = (c.center - p).lengthSq();
+			const double rSquared = (c.r * c.r);
+			const double err = Max(0.0, (dSquared - rSquared));
+
+			if (rSquared == 0)
+			{
+				return (err <= tolerance);
+			}
+
+			return (((err / rSquared) <= tolerance) || (err <= tolerance));
+		}
+
 		//
 		//	http://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment
 		//
@@ -1752,5 +1772,93 @@ namespace s3d
 			sum += ((p0.x - p3.x) * (p0.y + p3.y));
 			return (sum < 0);
 		}
+
+		//////////////////////////////////////////////////
+		//
+		//	SmallestEnclosingCircle
+		//
+		//////////////////////////////////////////////////
+	        //-----------------------------------------------
+	        //    Authors (OpenSiv3D challenge #19)
+	        //    - あぷりばーど
+	        //    - Nachia
+	        //    - Luke256
+	        //    - ラクラムシ
+	        //    - polyester
+	        //    - sasa
+	        //-----------------------------------------------
+
+		SIV3D_CONCEPT_URBG_
+		Circle SmallestEnclosingCircle(Array<Vec2> points, const double tolerance, URBG&& urbg)
+		{
+			if (points.size() == 0)
+			{
+				return Circle{};
+			}
+
+			if (points.size() == 1)
+			{
+				return Circle{ points[0], 0.0 };
+			}
+
+			if (points.size() == 2)
+			{
+				return Circle{ points[0], points[1] };
+			}
+
+			if (points.size() == 3)
+			{
+				return SmallestEnclosingCircle(points[0], points[1], points[2]);
+			}
+
+			if (points.size() == 4)
+			{
+				return SmallestEnclosingCircle(points[0], points[1], points[2], points[3], tolerance);
+			}
+
+			points.shuffle(std::forward<URBG>(urbg));
+
+			// 適当な 1 点を含む最小包含円から始めて、少しずつ広げていく戦略を取る。
+			// 含まれない点があったら、それが境界上になるように新たに取り直す。
+			Circle circle{ points[0], 0.0 };
+
+			for (size_t i = 1; i < points.size(); ++i)
+			{
+				const Vec2& p0 = points[i];
+
+				if (not detail::Contains(circle, p0, tolerance))
+				{
+					circle = Circle{ p0, 0.0 };
+
+					for (size_t j = 0; j < i; ++j)
+					{
+						const Vec2& p1 = points[j];
+
+						if (not detail::Contains(circle, p1, tolerance))
+						{
+							circle = Circle{ p0, p1 };
+
+							for (size_t k = 0; k < j; ++k)
+							{
+								const Vec2& p2 = points[k];
+
+								if (not detail::Contains(circle, p2, tolerance))
+								{
+									circle = Triangle(p0, p1, p2).getCircumscribedCircle();
+								}
+							}
+						}
+					}
+				}
+			}
+
+			return circle;
+		}
+
+		SIV3D_CONCEPT_URBG_
+		Circle SmallestEnclosingCircle(Array<Vec2> points, URBG&& urbg, double tolerance)
+		{
+			return SmallestEnclosingCircle(std::move(points), tolerance, std::forward<URBG>(urbg));
+		}
 	}
-}
+}
diff --git a/Siv3D/src/Siv3D/Geometry2D/SivGeometry2D.cpp b/Siv3D/src/Siv3D/Geometry2D/SivGeometry2D.cpp
@@ -5278,6 +5278,70 @@ namespace s3d
 			return detail::ConcaveHull(points, size, concavity, lengthThreshold);
 		}
 
+		//////////////////////////////////////////////////
+		//
+		//	SmallestEnclosingCircle
+		//
+		//////////////////////////////////////////////////
+	        //-----------------------------------------------
+	        //    Authors (OpenSiv3D challenge #19)
+	        //    - あぷりばーど
+	        //    - Nachia
+	        //    - Luke256
+	        //    - ラクラムシ
+	        //    - polyester
+	        //    - sasa
+	        //-----------------------------------------------
+
+		Circle SmallestEnclosingCircle(const Vec2& p0, const Vec2& p1, const Vec2& p2)
+		{
+			// 三角形 (p0, p1, p2) に対して鈍角の存在を判定し、もしあればその対辺（最長辺）を弦とする円が最小の円となる。
+			if ((p1 - p0).dot(p2 - p0) <= 0.0)
+			{ 
+				return Circle{ p1, p2 };
+			}
+
+			if ((p0 - p1).dot(p2 - p1) <= 0.0)
+			{
+				return Circle{ p0, p2 };
+			}
+
+			if ((p0 - p2).dot(p1 - p2) <= 0.0)
+			{
+				return Circle{ p0, p1 };
+			}
+
+			// 鋭角三角形の場合は (p0, p1, p2) の外接円が最小となる。
+			return Triangle{ p0, p1, p2 }.getCircumscribedCircle();
+		}
+
+		Circle SmallestEnclosingCircle(const Vec2& p0, const Vec2& p1, const Vec2& p2, const Vec2& p3, const double tolerance)
+		{
+			Circle circle = SmallestEnclosingCircle(p0, p1, p2);
+
+			if (not detail::Contains(circle, p3, tolerance))
+			{
+				circle = SmallestEnclosingCircle(p0, p1, p3);
+
+				if (not detail::Contains(circle, p2, tolerance))
+				{
+					circle = SmallestEnclosingCircle(p0, p2, p3);
+
+					if (not detail::Contains(circle, p1, tolerance))
+					{
+						circle = SmallestEnclosingCircle(p1, p2, p3);
+					}
+				}
+			}
+
+			return circle;
+		}
+
+		Circle SmallestEnclosingCircle(Array<Vec2> points, const double tolerance)
+		{
+			return SmallestEnclosingCircle(std::move(points), tolerance, GetDefaultRNG());
+		}
+
 		//////////////////////////////////////////////////
 		//
 		//	Subtract
@@ -5581,4 +5645,4 @@ namespace s3d
 
 		return { minX, minY, (maxX - minX), (maxY - minY) };
 	}
-}
+}