忘れていたわけじゃないんだよ。ただちょっとモチベーションが足りなかっただけ。
今日は、某所（CG 系の研究室）で需要がありそうだったので、以前 Processing で書いていた Delaunay 分割を C++ で実装し直してみる事にしました。
ソースコード上部の「view plain」という文字をクリックするとコピペ用の窓が出ますので、使いたい方はご自由にどうぞ。
【アルゴリズム部分（Delaunay2d.h）】
#ifndef ___TERCEL_DELAUNAY2D
#define ___TERCEL_DELAUNAY2D
#include <cfloat>
#include <cmath>
#include <list>
namespace Tercel
{
struct Vector
{
double x, y;
bool operator==(const Vector& v) const
{
return (x == v.x && y == v.y);
}
};
struct Circle
{
Vector center; // 中心座標
double radius; // 半径
};
class Triangle
{
public:
const Vector* p1, * p2, * p3; // 頂点座標
public:
// ======================================
// 等価性の判定
// ======================================
bool operator==(const Triangle& t) const
{
return(*p1 == *t.p1 && *p2 == *t.p2 && *p3 == *t.p3 ||
*p1 == *t.p2 && *p2 == *t.p3 && *p3 == *t.p1 ||
*p1 == *t.p3 && *p2 == *t.p1 && *p3 == *t.p2 ||
*p1 == *t.p3 && *p2 == *t.p2 && *p3 == *t.p1 ||
*p1 == *t.p2 && *p2 == *t.p1 && *p3 == *t.p3 ||
*p1 == *t.p1 && *p2 == *t.p3 && *p3 == *t.p2 );
}
// ======================================
// 他の三角形と共有点を持つか
// ======================================
bool hasCommonPoints(const Triangle& t)
{
return(*p1 == *t.p1 || *p1 == *t.p2 || *p1 == *t.p3 ||
*p2 == *t.p1 || *p2 == *t.p2 || *p2 == *t.p3 ||
*p3 == *t.p1 || *p3 == *t.p2 || *p3 == *t.p3 );
}
};
class Delaunay2d
{
private:
static void manageDuplicativeTriangles(std::list<Triangle>* newTriangleList,
std::list<Triangle>* duplicativeTriangleList,
const Triangle& t)
{
typedef std::list<Triangle>::iterator triIter;
bool existsInNewTriangleList = false;
for(triIter iter = newTriangleList->begin();
iter != newTriangleList->end(); iter++)
{
if(*iter == t) {
existsInNewTriangleList = true;
bool existsInDuplicativeTriangleList = false;
for(triIter iter2 = duplicativeTriangleList->begin();
iter2 != duplicativeTriangleList->end(); iter2++)
{
if(*iter2 == t)
{
existsInDuplicativeTriangleList = true;
break;
}
}
if(!existsInDuplicativeTriangleList)
{
duplicativeTriangleList->push_back(t);
}
break;
}
}
if(!existsInNewTriangleList) newTriangleList->push_back(t);
}
public:
static void getDelaunayTriangles(const std::list<Vector>& vertexList,
std::list<Triangle>* triangleList)
{
typedef std::list<Vector>::const_iterator cVtxIter;
typedef std::list<Triangle>::iterator triIter;
Triangle hugeTriangle;
{
// ======================================
// 外部三角形を作る
// ======================================
double maxX, maxY; maxX = maxY = DBL_MIN;
double minX, minY; minX = minY = DBL_MAX;
for(cVtxIter it = vertexList.begin(); it != vertexList.end(); ++it)
{
double x = it->x;
double y = it->y;
if(maxX < x) maxX = x; if(minX > x) minX = x;
if(maxY < y) maxY = y; if(minY > y) minY = y;
}
// すべての点を包含する矩形の外接円
double centerX = (maxX - minX) * 0.5; // 中心x座標
double centerY = (maxY - minY) * 0.5; // 中心y座標
double dX = maxX - centerX;
double dY = maxY - centerY;
double radius = sqrt(dX * dX + dY * dY) + 1.0; // 半径
Vector* p1 = new Vector; // メモリ確保（266行目で解放）
p1->x = centerX - sqrt(3.0) * radius;
p1->y = centerY - radius;
Vector* p2 = new Vector; // メモリ確保（267行目で解放）
p2->x = centerX + sqrt(3.0) * radius;
p2->y = centerY - radius;
Vector* p3 = new Vector; // メモリ確保（268行目で解放）
p3->x = centerX;
p3->y = centerY + 2.0 * radius;
hugeTriangle.p1 = p1;
hugeTriangle.p2 = p2;
hugeTriangle.p3 = p3;
}
triangleList->push_back(hugeTriangle);
// --------------------------------------
// 点を逐次添加し、反復的に三角分割を行う
// --------------------------------------
for(cVtxIter vIter = vertexList.begin(); vIter != vertexList.end(); ++vIter)
{
const Vector* p = &(*vIter);
// --------------------------------------
// 追加候補の三角形を保持する一時リスト
// --------------------------------------
std::list<Triangle> newTriangleList; // 新規分割された三角形
std::list<Triangle> duplicativeTriangleList; // 重複リスト
// --------------------------------------
// 現在の三角形セットから要素を一つずつ取り出して、
// 与えられた点が各々の三角形の外接円の中に含まれるかどうか判定
// --------------------------------------
for(triIter tIter = triangleList->begin(); tIter != triangleList->end(); )
{
// 三角形セットから三角形を取りだして…
Triangle t = *tIter;
// その外接円を求める。
Circle c;
{
// 三角形の各頂点座標を (x1, y1), (x2, y2), (x3, y3) とし、
// その外接円の中心座標を (x, y) とすると、
// (x - x1) * (x - x1) + (y - y1) * (y - y1)
// = (x - x2) * (x - x2) + (y - y2) * (y - y2)
// = (x - x3) * (x - x3) + (y - y3) * (y - y3)
// より、以下の式が成り立つ
//
// x = { (y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
// + (y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c
//
// y = { (x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
// + (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c
//
// ただし、
// c = 2 * {(x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)}
double x1 = t.p1->x; double y1 = t.p1->y;
double x2 = t.p2->x; double y2 = t.p2->y;
double x3 = t.p3->x; double y3 = t.p3->y;
double m = 2.0 * ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1));
double x = ((y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
+ (y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)) / m;
double y = ((x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
+ (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)) / m;
c.center.x = x;
c.center.y = y;
// 外接円の半径 r は、半径から三角形の任意の頂点までの距離に等しい
double dx = t.p1->x - x;
double dy = t.p1->y - y;
double radius = sqrt(dx * dx + dy * dy);
c.radius = radius;
}
double dx = c.center.x - p->x;
double dy = c.center.y - p->y;
double dist = sqrt(dx * dx + dy * dy);
// ======================================
// 一時リストを使って重複判定
// ======================================
if(dist < c.radius)
{
// 再分割
Triangle t1;
t1.p1 = p; t1.p2 = t.p1; t1.p3 = t.p2;
manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t1);
Triangle t2;
t2.p1 = p; t2.p2 = t.p2; t2.p3 = t.p3;
manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t2);
Triangle t3;
t3.p1 = p; t3.p2 = t.p3; t3.p3 = t.p1;
manageDuplicativeTriangles(&newTriangleList, &duplicativeTriangleList, t3);
tIter = triangleList->erase(tIter);
}
else ++tIter;
}
// --------------------------------------
// 一時リストのうち、重複のないものを三角形リストに追加
// --------------------------------------
for(triIter iter = newTriangleList.begin();
iter != newTriangleList.end(); ++iter)
{
bool exists = false;
for(triIter iter2 = duplicativeTriangleList.begin();
iter2 != duplicativeTriangleList.end(); ++iter2)
{
if(*iter == *iter2)
{
exists = true;
break;
}
}
if(!exists) triangleList->push_back(*iter);
}
}
// 最後に、外部三角形の頂点を削除
for(triIter tIter = triangleList->begin(); tIter != triangleList->end(); )
{
if(hugeTriangle.hasCommonPoints(*tIter)) tIter = triangleList->erase(tIter);
else ++tIter;
}
// 巨大三角形の頂点を解放
delete hugeTriangle.p1;
delete hugeTriangle.p2;
delete hugeTriangle.p3;
}
};
}
#endif
【テスト用の GUI （Main.cpp）】
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 結果表示用のウィンドウ
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#include <windows.h>
#include "Delaunay2d.h"
// アプリケーションとウィンドウの初期化に必要な定数
// ----------------------------------------
const LPCTSTR APP_NAME = TEXT("APP_TERCEL_TECH");
const LPCTSTR MUTEX_NAME = TEXT("MUTEX_TERCEL_TECH");
// クライアント領域の幅と高さ
// ----------------------------------------
const UINT WIDTH = 400;
const UINT HEIGHT = 300;
// ----------------------------------------
// 頂点集合、およびそれらを基に
// Delaunay 分割された三角形を格納しておくSTLコンテナ
// ----------------------------------------
std::list<Tercel::Vector> vertices;
std::list<Tercel::Triangle> triangles;
// コールバック関数（メッセージ処理）
// ========================================
LRESULT CALLBACK WindowProc(HWND hWnd, UINT uMsg,
WPARAM wParam, LPARAM lParam)
{
switch(uMsg)
{
case WM_DESTROY:
PostQuitMessage(0);
return 0;
case WM_CREATE:
{
// ======================================
// 頂点集合にランダムでデータをセット
// ======================================
for(int i = 0; i < 200; ++i) {
Tercel::Vector v; // ベクトルを宣言
v.x = rand() % WIDTH;
v.y = rand() % HEIGHT;
vertices.push_back(v);
}
// ======================================
// Delaunay分割を行う関数
// ======================================
// 結果は第2引数（std::list<Tercel::Triangle> 型）に格納される。
//
// Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、
// いずれも Tercel::Vertex オブジェクトへのポインタになっている
// これらのポインタは、vertices リストの要素を参照している。
Tercel::Delaunay2d::getDelaunayTriangles(vertices, &triangles);
// ~~~~~~~~ ~~~~~~~~~~
// 頂点集合 三角形集合
// （参照） （ポインタ）
}
return 0;
case WM_PAINT:
{
typedef std::list<Tercel::Triangle> triangleList;
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
HPEN hPen = CreatePen(PS_SOLID, 1, RGB(0, 0, 0));
SelectObject(hdc, hPen);
// ======================================
// Delaunay分割された三角形の集合を一つ一つ描画
// ======================================
for(triangleList::iterator it = triangles.begin();
it != triangles.end(); ++it)
{
Tercel::Triangle t = *it; // 三角形取得
int x, y;
x = int(t.p1->x + 0.5);
y = int(t.p1->y + 0.5);
MoveToEx(hdc, x, y, NULL);
x = int(t.p2->x + 0.5);
y = int(t.p2->y + 0.5);
LineTo(hdc, x, y);
MoveToEx(hdc, x, y, NULL);
x = int(t.p3->x + 0.5);
y = int(t.p3->y + 0.5);
LineTo(hdc, x, y);
MoveToEx(hdc, x, y, NULL);
x = int(t.p1->x + 0.5);
y = int(t.p1->y + 0.5);
LineTo(hdc, x, y);
}
DeleteObject(hPen);
EndPaint(hWnd, &ps);
}
return 0;
}
return DefWindowProc(hWnd, uMsg, wParam, lParam);
}
// エントリポイント
// ========================================
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
LPSTR lpCmdLine, int nCmdShow)
{
// 多重起動防止
HANDLE hMutex = CreateMutex(NULL, TRUE, MUTEX_NAME);
if(GetLastError() == ERROR_ALREADY_EXISTS)
{
// 多重起動を検知した際には、
// 既存のウィンドウを最前面に表示してプログラムを終了する
HWND existingWnd = FindWindow(APP_NAME, NULL);
if(existingWnd != NULL)
{
if(IsIconic(existingWnd))
ShowWindowAsync(existingWnd, SW_RESTORE);
SetForegroundWindow(existingWnd);
}
return -1;
}
// ウィンドウクラス登録
WNDCLASSEX wc;
wc.cbSize = sizeof(WNDCLASSEX);
wc.style = CS_HREDRAW | CS_VREDRAW;
wc.lpfnWndProc = WindowProc;
wc.cbClsExtra = 0;
wc.cbWndExtra = sizeof(LONG_PTR);
wc.hInstance = hInstance;
wc.hIcon = (HICON)LoadImage(NULL,
MAKEINTRESOURCE(IDI_APPLICATION),
IMAGE_ICON,
0,
0,
LR_DEFAULTSIZE | LR_SHARED);
wc.hCursor = (HICON)LoadImage(NULL,
MAKEINTRESOURCE(IDC_ARROW),
IMAGE_CURSOR,
0,
0,
LR_DEFAULTSIZE | LR_SHARED);
wc.hbrBackground = (HBRUSH)COLOR_BACKGROUND + 1;
wc.lpszMenuName = NULL;
wc.lpszClassName = APP_NAME;
wc.hIconSm = (HICON)LoadImage(NULL,
MAKEINTRESOURCE(IDI_APPLICATION),
IMAGE_ICON,
0,
0,
LR_DEFAULTSIZE | LR_SHARED);
if(!RegisterClassEx(&wc)) return -1;
// クライアント領域のサイズから、ウィンドウそのものの幅と高さを計算
UINT windowWidth = WIDTH + GetSystemMetrics(SM_CXFIXEDFRAME) * 2;
UINT windowHeight = HEIGHT + GetSystemMetrics(SM_CYFIXEDFRAME) * 2
+ GetSystemMetrics(SM_CYCAPTION);
// ウィンドウ生成
HWND hWnd = CreateWindow(APP_NAME,
TEXT("はじめてのDelaunay分割"),
WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX,
CW_USEDEFAULT,
CW_USEDEFAULT,
windowWidth,
windowHeight,
NULL,
NULL,
hInstance,
NULL);
if(!hWnd) return -1;
// ウィンドウ表示
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
// メッセージループ
MSG msg;
while(TRUE)
{
if(PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))
{
if(msg.message == WM_QUIT) break;
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
CloseHandle(hMutex);
return (int)msg.wParam;
}
テスト用の GUI プログラムの方は、思いっきり Windows 依存です。
ただ、Delaunay 分割のアルゴリズムは完全に分離してあるので、他の環境で使いたい場合は Delaunay2d.h だけを持っていって、表示部はがんばって自分で作れば OK です。
肝心の使い方は以下の通りです（Main.cpp のコールバック関数内に実際のサンプルがあります）。
- 頂点リスト（std::list
- 先ほどの頂点リストに、適当な頂点オブジェクトを格納します。
- 三角形リスト（std::list
- Tercel::Delaunay2d::getDelaunayTriangles() 関数を呼びます。
一見すると奇妙な仕様ですが、関数内で中身が変更されるオブジェクトに関してはポインタ渡し、そうでないものは参照渡しで構文を使い分けているためです。
特に関数を呼び出す側（ここでは Main.cpp ）から見ると、参照渡しは値渡しのように見えるため、関数の中でオブジェクトが変更されないような印象になります。
Delaunay 分割された結果は三角形リストに格納されます。三角形（Tercel::Triangle 構造体）は、メンバに頂点（Tercel::Vertex オブジェクト）のポインタを保持します。ここでは、最初に作った頂点リストの要素を参照する形になります。
ですので、頂点リストの要素を削除した場合、三角形の頂点は不正なポインタになってしまう事に注意する必要があります。
0 件のコメント:
コメントを投稿
ひとことどうぞφ(・ω・,,)