どうも、脳内で Template Metaprogramming が流行っている melt です。こんばんは。
といいつつも今日は Template Metaprogramming は全く関係なくて、凸包判定についてのお話です。
自分は最近ゲームプログラミングのためのリアルタイム衝突判定という本を読み進めています。
当然内容なんてほとんど分かっていないのですが、ちょっとだけ理解できる部分に面白いことが書いてありました。
多角形が凸かどうかを判定する方法をネット上で調べてみると、各線分ベクトルの外積の符号が全て同じ(内角が180度未満)であれば凸多角形として判断している場合が多いのだけれども、この判断だけだと十分ではない、ということが書かれていました。
五芒星なんか考えると分かりやすいと思います。

確かに内角が全て180度未満ですが、明らかに凸多角形では無いです。
で、これはどうやって解決すればいいのかなぁと思って続きを読んでみたのですが、この本、参考文献だけ示してソースが載ってないのです。
どうしようかなと思ってとりあえず参考文献の情報を調べてみると、どうやらホームページ上でソースコードが公開されているとのこと。
ここから辿れる中の、convex_opt.cというのがそれに当たるようです。
英語読めない人なのでどういう仕組みなのかあまり理解出来ていないのですが、とりあえずこれを C# 用に書き直してテストしてみると、五芒星はちゃんと NotConvex になりました。
他にも点が同じ場所にあるような場合とか複数の点が同一直線上にあるようなケースも試してみましたが、ちゃんと正常に判断されていました。テラスゴス。
以下 C# 用に書き直したコード
public struct Point
{
public Point(double x, double y)
{
X = x;
Y = y;
}
public double X;
public double Y;
}
public enum ConvexType
{
NotConvex = -2, // 凸包ではない
NotConvexDegenerate = -1, // 凸包ではなく、縮退している
ConvexDegenerate = 0, // 凸包であるが、縮退している
ConvexCCW = 1, // 反時計回りの凸包
ConvexCW = 2, // 時計回りの凸包
}
public static ConvexType IsConvex(Point[] points)
{
int nvert = points.Length;
if (nvert < 3)
{
return ConvexType.ConvexDegenerate;
}
int second;
Point dprev;
int iread = 1;
while (true)
{
second = iread++;
dprev = new Point(
points[second].X - points[0].X,
points[second].Y - points[0].Y);
if (dprev.X != 0 || dprev.Y != 0)
{
break;
}
if (iread >= nvert)
{
return ConvexType.ConvexDegenerate;
}
}
int saveSecond = second;
int curDir =
(dprev.X < 0.0) ? -1 :
(dprev.X > 0.0) ? 1 :
(dprev.Y < 0.0) ? -1 :
(dprev.Y > 0.0) ? 1 : 0;
int dirChanges = 0;
int angleSign = 0;
while (true)
{
if (iread >= nvert)
{
iread = 0;
}
int third = iread++;
Point dcur = new Point(
points[third].X - points[second].X,
points[third].Y - points[second].Y);
if (dcur.X == 0.0 && dcur.Y == 0.0)
{
continue;
}
int thisDir =
(dcur.X < 0.0) ? -1 :
(dcur.X > 0.0) ? 1 :
(dcur.Y < 0.0) ? -1 :
(dcur.Y > 0.0) ? 1 : 0;
if (thisDir == -curDir)
{
++dirChanges;
}
curDir = thisDir;
double cross = dprev.Y * dcur.X - dprev.X * dcur.Y;
if (cross > 0)
{
if (angleSign == -1)
{
return ConvexType.NotConvex;
}
angleSign = 1;
}
else if (cross < 0)
{
if (angleSign == 1)
{
return ConvexType.NotConvex;
}
angleSign = -1;
}
second = third;
dprev = dcur;
if (saveSecond == second)
{
break;
}
}
if (dirChanges > 2)
{
if (angleSign == 0)
{
return ConvexType.NotConvexDegenerate;
}
else
{
return ConvexType.NotConvex;
}
}
if (angleSign > 0)
{
return ConvexType.ConvexCCW;
}
if (angleSign < 0)
{
return ConvexType.ConvexCW;
}
return ConvexType.ConvexDegenerate;
}
C 版の classifyPolygon2 と比べて無駄に比較している部分があるので少しだけ遅くなってしまっていますが、許容範囲かなぁという感じです。
あと↑のコードは 0.0 と比較していますが、ほんとはある程度小さい値(Epsilon)より小さければ 0 として考える方がいいと思うんですけどどうなんでしょう。