C# で２次元の凸包（グラハムスキャン）を計算する
こんにちは。
仕事の自動化にやりがいと達成感を感じるガッくんです。
この記事の目次
背景・目的
板の上に数本の釘が刺さっていて、その全ての釘の外側から輪ゴムをかけた時に出来る多角形を求めることを凸包と言います。
簡単なものであれば、パターンマッチングや点の範囲を割り出す時に使えます。
範囲と範囲の重なりや離れ具合なんかを計算する時にも使えます。
２次元の平面上で凸包された輪郭線を検出したくなったので C# で作成しました。
動作環境
・Windows 7
・Visual Studio 2017
・OpenCvSharp
プログラム
仕様
以前、３次元点群データを作成する記事を書きました。C# でテスト用３次元点群データを作成する - 解析エンジニアの自動化 blog
せっかくなので、このプログラムで作成した点群データを使って、凸包します。
記事の通りのプログラムを実行すれば、デスクトップに Test-Nodes-Data.txt というファイルが作成されます。
プログラムを使う流れを考えます。
①上の記事のプログラムで作った点の情報が記述されたテキストファイルをＧＵＩのテキストボックスにドラッグ&ドロップする。
②実行ボタンを押す。
以上！！
この仕様に見合ったＧＵＩを考えます。
ＧＵＩ
作成したＧＵＩは図1 の通りです。２種類のコントロールを合計３つ設置しただけのとても簡単なものです。
・テキストボックス
・ボタン（実行）
・ボタン（キャンセル）
ソースコード
描画は OpenCvSharp を使いました。OpenGL（OpenTK）でも良かったのですが、後で OpenCV を使って凸包の結果に対して、さらに何らかの処理をする事を考えて、 OpenCV にしました。
今回はソースコードが長いので、もし、確認するのであれば、 Visual Studio にコピペして確認することをお勧めします。
using System;
using System.Collections;
using System.Windows.Forms;
using System.IO;
using OpenCvSharp.CPlusPlus;
namespace ConvexHuller
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
private void DD(object sender, DragEventArgs e)
{
// ドロップファイルの取得
string[] args = (string[])e.Data.GetData(DataFormats.FileDrop, false);
// ドロップファイル数の取得
int args_cnt = args.Length;
// ドロップファイル数のチェック
if (args_cnt > 1)
{
MessageBox.Show("引数が多すぎます。" + Environment.NewLine + "ファイルを１つだけドロップしてください。");
return;
}
// sender をテキストボックスとして変数化
TextBox txtTgt = sender as TextBox;
// sender がテキストボックスではないなら終了
if (txtTgt == null)
{
return;
}
// テキストボックスにファイル名出力
txtTgt.Text = args[0];
}
private void DE(object sender, DragEventArgs e)
{
if (e.Data.GetDataPresent(DataFormats.FileDrop))
{
e.Effect = DragDropEffects.All;
}
else
{
e.Effect = DragDropEffects.None;
}
}
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
// 凸包
// ConvexHull(points)
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
private double[][] ConvexHull(double[][] points)
{
// 引き渡されたジャグ配列が何行あるかを取得
int Num = points.Length;
// 同じ行数のジャグ配列の作成
double[][] convex = new double[Num][];
// Y の値が小さい、かつ、 X の値が大きいキー座標を探す
int limit = 0; // キー座標の要素番号記憶用変数
for (int i = 1; i < Num; i++)
{
// 小さい Y の値を見つけたら、要素番号を記憶する
if (points[limit][1] > points[i][1]) { limit = i; }
// Y の値が同じであっても、 X の値が大きければ要素番号を記憶する
else if (points[limit][1] == points[i][1] && points[limit][0] < points[i][0])
{ limit = i; }
}
// キー座標からみた角度を求める 以下のような並び順の配列を作る。ただし、キー座標の角度はゼロにする
// 角度 要否ID X座標 Y座標（３次元座標を持つ場合は凸包したい平面における座標値とする）
// なお、要否ID => 要:0 否:-1
double theta = 0; // 角度
for (int i = 0; i < Num; i++)
{
// キー座標なら、角度=0、要否ID=0、X座標、Y座標
if (limit == i)
{
convex[i] = new double[] { 0, 0, points[i][0], points[i][1] };
}
// キー座標以外は角度を算出して配列の作成
else
{
theta = Math.Atan2(points[i][1] - points[limit][1], points[i][0] - points[limit][0]);
if (theta < 0) { theta = 2 * Math.PI + theta; }
convex[i] = new double[] { theta, 0, points[i][0], points[i][1] };
}
}
// ゼロ列目の角度の大きい順にソートする
Array.Sort(convex, StructuralComparisons.StructuralComparer);
// 角HIJ 計算：角度がマイナスのI点はIDを-1として事実上、無きものとする
int del = 0; // 要素番号の増減用変数
int h; //
int j; //
double theta1 = 0; //
double theta2 = 0; //
int cntMinus = 0; // 要否IDが否になった要素数カウンタ
for (int i = 1; i < Num; i++)
{
if (convex[i][1] == -1)
{
while (convex[i][1] == -1) { i += del; }
}
del = 1; // 次のループのために初期化する
h = i - del; // i点の1個前の点の要素番号を設定する
j = i + del; // i点の1個後の点の要素番号を設定する
// h 点の探索
// 要否ID が「否」ならさらに1個前の要素番号を調べて要否ID が「要」の要素を見つける
while (convex[h][1] == -1) { del++; h = i - del; }
del = 1; // 次のループのために初期化する
// i点の1個後のj点の要素番号が配列の外に到達してしまったら
if (j == Num) { j = 0; }
// j 点の探索
// 要否ID が「否」ならさらに1個後の要素番号を調べて要否ID が「要」の要素を見つける
while (convex[j][1] == -1)
{
del++;
j = i + del;
// i点の1個後のj点の要素番号が配列の外に到達してしまったら
if (j == Num) { j = 0; }
}
del = 1; // 次のループのために初期化する
// i-h角度の算出
theta1 = Math.Atan2(convex[i][3] - convex[h][3], convex[i][2] - convex[h][2]);
if (theta1 < 0) { theta1 = 2 * Math.PI + theta1; }
// j-i角度の算出
theta2 = Math.Atan2(convex[j][3] - convex[i][3], convex[j][2] - convex[i][2]);
if (theta2 < 0) { theta2 = 2 * Math.PI + theta2; }
// 角度がマイナスなら要否ID を「否」にして要素数を減少させる
if (theta2 - theta1 < 0)
{
convex[i][1] = -1; cntMinus++;
while (convex[i][1] == -1) { i -= del; }
i--;
}
del = 1; // 次のループのために初期化する
}
// 戻り値用に配列に入れ直す
// ジャグ配列の準備
double[][] ret = new double[Num - cntMinus][];
int cnt = 0; // 要素番号カウンタ
for (int i = 0; i < convex.Length; i++)
{
// 要素ID が「否」なら戻り値から除外する
if (convex[i][1] != -1)
{
ret[cnt] = new double[] { convex[i][2], convex[i][3] };
cnt++;
}
}
// 戻り値
return ret;
}
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
// OpenCV
// 凸包描画関数
// PlotConvexHull(img, points, Ox, Oy)
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
private void PlotConvexHull(Mat img, double[][] points, double Ox, double Oy)
{
for (int i = 0; i < points.Length; i++)
{
// ジャグ配列で最後の要素はジャグ配列の先頭の要素と直線で結ぶため、要素番号を恣意的に操作して直線で結ぶ
if (i + 1 == points.Length)
{
Cv2.Line(img, new OpenCvSharp.CPlusPlus.Point(Ox + points[i][0], Oy - points[i][1]), new OpenCvSharp.CPlusPlus.Point(Ox + points[0][0], Oy - points[0][1]), new Scalar(0, 0, 255), 1);
}
// ジャグ配列の今のi要素とi+1要素を直線で結ぶ
else
{
Cv2.Line(img, new OpenCvSharp.CPlusPlus.Point(Ox + points[i][0], Oy - points[i][1]), new OpenCvSharp.CPlusPlus.Point(Ox + points[i + 1][0], Oy - points[i + 1][1]), new Scalar(0, 0, 255), 1);
}
}
}
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
// OpenCV
// ×点の OpenCV 出力
// Cv2PlotPoint(img, points, Ox, Oy)
//□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
private void Cv2PlotPoint(Mat img, double[][] point, double Ox, double Oy)
{
// ×印で点を描画する際の×の長さ設定
double xLength = 1.5;
double deltaX = xLength / Math.Sqrt(2);
double deltaY = xLength / Math.Sqrt(2);
// ×印で点を描画する
for (int i = 0; i < point.Length; i++)
{
// ×の＼の描画
Cv2.Line(img, new OpenCvSharp.CPlusPlus.Point(Ox + point[i][0] - deltaX, Oy - point[i][1] - deltaY), new OpenCvSharp.CPlusPlus.Point(Ox + point[i][0] + deltaX, Oy - point[i][1] + deltaY), new Scalar(255, 255, 255));
// ×の／の描画
Cv2.Line(img, new OpenCvSharp.CPlusPlus.Point(Ox + point[i][0] + deltaX, Oy - point[i][1] - deltaY), new OpenCvSharp.CPlusPlus.Point(Ox + point[i][0] - deltaX, Oy - point[i][1] + deltaY), new Scalar(255, 255, 255));
}
}
//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
// 実行ボタン
//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
private void button1_Click(object sender, EventArgs e)
{
// このメソッドへのパラメータ
string arg = textBox1.Text;
string sDir;
string sFile;
string[] sReadLines;
int cntNodes = 0;
string[] delimLine;
double sX;
double sZ;
double minX = 0; double maxX = 0; double minZ = 0; double maxZ = 0;
// 引数ファイルの情報取得
sDir = Path.GetDirectoryName(arg);
sFile = Path.GetFileNameWithoutExtension(arg);
// 引数ファイルの読み込み
sReadLines = File.ReadAllLines(arg);
// ノード情報の数をカウント
for (int j = 0; j < sReadLines.Length; j++)
{
if (sReadLines[j].IndexOf('*') == -1)
{
cntNodes++;
}
}
// ノード数分のジャグ配列を用意する
double[][] coorXZ = new double[cntNodes][];
// ノードの X , Z 座標ジャグ配列を作成
cntNodes = 0;
for (int j = 0; j < sReadLines.Length; j++)
{
if (sReadLines[j].IndexOf('*') == -1)
{
delimLine = sReadLines[j].Split(',');
sX = double.Parse(delimLine[1]);
sZ = double.Parse(delimLine[3]);
if (minX > sX) { minX = sX; }
if (maxX < sX) { maxX = sX; }
if (minZ > sZ) { minZ = sZ; }
if (maxZ < sZ) { maxZ = sZ; }
coorXZ[cntNodes] = new double[] { sX, sZ };
cntNodes++;
}
}
// 凸包
double[][] outXZ = ConvexHull(coorXZ);
// 原点とマージン設定
double Ox = 0; double Oz = 0; double margin = 25;
// 画面サイズ
double justHeight = margin + Math.Abs(maxZ - minZ) + margin;
double justWidth = margin + Math.Abs(maxX - minX) + margin;
int screenHeight = (int)Math.Ceiling(justHeight); // 小数点切り上げ
int screenWidth = (int)Math.Ceiling(justWidth); // 小数点切り上げ
if (minX < 0) { Ox = Math.Abs(0 - minX) + margin; }
else { Ox = -minX + margin; }
if (minZ < 0) { Oz = Math.Abs(0 - minZ) + margin; }
else { Oz = screenHeight - minZ + margin; }
// Mat クラスで画像を作成する
Mat img1 = new Mat(screenHeight, screenWidth, MatType.CV_8UC3, new Scalar(0, 0, 0));
Mat imgGray1 = new Mat(screenHeight, screenWidth, MatType.CV_8SC1, new Scalar(0, 0, 0));
Mat imgBinary1 = new Mat(screenHeight, screenWidth, MatType.CV_8SC1, new Scalar(0, 0, 0));
// Mat にノード
Cv2PlotPoint(img1, coorXZ, Ox, Oz);
// Mat に凸包描画
PlotConvexHull(img1, outXZ, Ox, Oz);
Cv2.ImShow("image1", img1);
}
//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
// キャンセルボタン
//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
private void button2_Click(object sender, EventArgs e)
{
this.Close();
}
}
}
結果
３次元点群を X-Z 平面的に見て凸包計算が出来ました。
コメント
３次元の点群データから X 、 Z 座標を選び、平面的に扱いましたが、 X 、 Y 座標でももちろん出来ます。
必要に応じてソースコードを変更してください。
個人的には面白かったので、３次元凸包も試してみたくなりました。
以上
