続： C++ で Delaunay 分割

昨日、C++ で Delaunay 分割を実装しました。

で。昨日はデータ構造に線形リスト（std::list）のみを用いており、コーディングしやすい代わりにあまり効率的とはいえない実装になっていました。

今回は、リストの代わりにセット（std::set）を用いる事で、内部処理の効率化を図りました。セットの実体は二分探索木（おそらくは赤黒木）であり、『重複する要素を持たない』『要素の探索が対数時間』というなかなかおいしいデータ構造のため、昨日の段階でネックになっていた冗長な重複判定処理がかなり省かれています。

ではなぜ最初からそれを使わなかったかというと、ずばり std::set の使い方をよく知らなかったからです。……恥ずかしい。

前述の通り、セットの内部構造は二分探索木ですので、セットに格納する要素同士の大小関係を判定できなくてはなりません。単純に比較できない『頂点』や『三角形』の大小関係をどう定義するかで悩みましたが、結局辞書式順序で順序づける事にしました（Tercel::Vector や Tercel::Triangle 構造体の operator< が、比較関数として機能します）。

ソースコードは以下。

---

【Delaunay2d.h】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// Delaunay分割（やや効率向上？版）
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


#ifndef ___TERCEL_DELAUNAY2D
#define ___TERCEL_DELAUNAY2D

#include <cfloat>
#include <cmath>

#include <map>
#include <set>

namespace Tercel 
{
  struct Vector {
    double x;
    double y;

    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Vector& v) const 
    {
      return (x == v.x && y == v.y);
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Vector& v) const 
    {
      return x != v.x ? x < v.x : y < v.y;
    }
  };


  struct Circle 
  {
    Vector center;  // 中心座標
    double radius;  // 半径
  };


  class Triangle {
  public:
    const Vector* p1, * p2, * p3;  // 頂点座標

  private:
    // 最小の頂点を返す   
    // --------------------------------------
    inline const Vector* getMinVertex() const 
    {
      return *p1 < *p2 && *p1 < *p3 ? p1 : *p2 < *p3 ? p2 : p3;
    }

    // まんなかの頂点を返す
    // --------------------------------------
    inline const Vector* getMidVertex() const 
    {
      return *p1 < *p2 && *p2 < *p3 || *p3 < *p2 && *p2 < *p1 ? p2 : 
             *p2 < *p3 && *p3 < *p1 || *p1 < *p3 && *p3 < *p2 ? p3 : p1;
    }

    // 最大の頂点を返す
    // --------------------------------------
    inline const Vector* getMaxVertex() const 
    {
      return *p2 < *p1 && *p3 < *p1 ? p1 : *p3 < *p2 ? 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;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Triangle& t) const 
    {
      return !(*getMinVertex() == *t.getMinVertex()) ? 
                 *getMinVertex() < *t.getMinVertex() :
             !(*getMidVertex() == *t.getMidVertex()) ? 
                 *getMidVertex() < *t.getMidVertex() :
                 *getMaxVertex() < *t.getMaxVertex();
    }

    // ======================================  
    // 他の三角形と共有点を持つか  
    // ====================================== 
    bool hasCommonPoints(const Triangle& t) const
    {
      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 
  {
    // ======================================  
    // 型定義
    // ====================================== 
    typedef const std::set<Vector>         ConstVertexSet;
    typedef ConstVertexSet::const_iterator ConstVertexIter;

    typedef std::set<Triangle>             TriangleSet;
    typedef std::set<Triangle>::iterator   TriangleIter;

    typedef std::map<Triangle, bool>       TriangleMap;

  private:
    // ======================================  
    // 一時マップを使って重複判定  
    // hashMap  
    //  - Key : 三角形  
    //  - Value : 重複していないかどうか  
    //        - 重複していない : true  
    //        - 重複している   : false  
    // ======================================
    static inline void addElementToRedundanciesMap(TriangleMap* pRddcMap,
                                                   Triangle& t) 
    {
      TriangleMap::iterator it = pRddcMap->find(t);
      if(it != pRddcMap->end() && it->second) 
      {
        // 値を (t, true) から (t, false) に置換
        pRddcMap->erase(it);
        pRddcMap->insert(TriangleMap::value_type(t, false));
      } 
      else 
      {
        pRddcMap->insert(TriangleMap::value_type(t, true));
      }
    }

  public:
    static void getDelaunayTriangles(ConstVertexSet& pVertexSet, 
                                     TriangleSet* triangleSet) 
    {

      Triangle hugeTriangle;
      {
        // ======================================  
        // 外部三角形を作る  
        // ======================================  
        double maxX, maxY; maxX = maxY = DBL_MIN;
        double minX, minY; minX = minY = DBL_MAX;
        for(ConstVertexIter it = pVertexSet.begin(); 
            it != pVertexSet.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;  // メモリ確保（314行目で解放）
        p1->x    = centerX - sqrt(3.0) * radius;
        p1->y    = centerY - radius;

        Vector* p2 = new Vector;  // メモリ確保（315行目で解放）
        p2->x    = centerX + sqrt(3.0) * radius;
        p2->y    = centerY - radius;

        Vector* p3 = new Vector;  // メモリ確保（316行目で解放）
        p3->x    = centerX;
        p3->y    = centerY + 2.0 * radius;

        hugeTriangle.p1 = p1;
        hugeTriangle.p2 = p2;
        hugeTriangle.p3 = p3;
      }

      triangleSet->insert(hugeTriangle);

      // --------------------------------------  
      // 点を逐次添加し、反復的に三角分割を行う  
      // --------------------------------------  
      for(ConstVertexIter vIter = pVertexSet.begin(); 
          vIter != pVertexSet.end(); ++vIter) 
      {
        const Vector* p = &*vIter;

        // --------------------------------------  
        // 追加候補の三角形を保持する一時マップ
        // -------------------------------------- 
        TriangleMap rddcMap;

        // --------------------------------------  
        // 現在の三角形セットから要素を一つずつ取り出して、  
        // 与えられた点が各々の三角形の外接円の中に含まれるかどうか判定  
        // --------------------------------------  
        for(TriangleIter tIter = triangleSet->begin(); 
            tIter != triangleSet->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) } / m
            //  
            // y = { (x1-x3) * (x2*x2-x1*x1+y2*y2-y1*y1)  
            //     + (x2-x1) * (x3*x3-x1*x1+y3*y3-y1*y1) } / m  
            //  
            // ただし、  
            // m = 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;
            addElementToRedundanciesMap(&rddcMap, t1);
 
            Triangle t2;
            t2.p1 = p; t2.p2 = t.p2; t2.p3 = t.p3;
            addElementToRedundanciesMap(&rddcMap, t2);
 
            Triangle t3;
            t3.p1 = p; t3.p2 = t.p3; t3.p3 = t.p1;
            addElementToRedundanciesMap(&rddcMap, t3);

            triangleSet->erase(tIter++);
          } 
          else ++tIter;
        }

        // --------------------------------------  
        // 一時マップのうち、重複のないものを三角形リストに追加   
        // --------------------------------------  
        for(TriangleMap::iterator iter = rddcMap.begin();
          iter != rddcMap.end(); ++iter) 
        {
          if(iter->second) triangleSet->insert(iter->first);
        }
      }

      // --------------------------------------  
      // 最後に、外部三角形の頂点を削除
      // --------------------------------------  
      for(TriangleIter tIter = triangleSet->begin(); 
          tIter != triangleSet->end(); ) 
      {
        if(hugeTriangle.hasCommonPoints(*tIter))
          triangleSet->erase(tIter++);
        else ++tIter;
      }

      // 巨大三角形の頂点を解放
      delete hugeTriangle.p1;
      delete hugeTriangle.p2;
      delete hugeTriangle.p3;
    }
  };
}

#endif

ちなみに使い方（というか Windows の GUI による結果表示用プログラム）はこちら。

【Main.cpp】

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// Delaunay分割（やや効率化版）
//           結果表示用のウィンドウ
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#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::set<Tercel::Vector>   vertices;  // ←型に注意！！
std::set<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.insert(v); // ←変更！！
      }

      // ======================================  
      // Delaunay分割を行う関数
      // ======================================  
      
      // 結果は第2引数（std::set<Tercel::Triangle> 型）に格納される。
      // 
      // Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、
      // いずれも Tercel::Vertex オブジェクトへのポインタになっている
      // これらのポインタは、vertices リストの要素を参照している。
      Tercel::Delaunay2d::getDelaunayTriangles(vertices, &triangles);
      //                     ~~~~~~~~  ~~~~~~~~~~
      //                     頂点集合  三角形集合
      //                     （参照） （ポインタ）
    }
    return 0;

  case WM_PAINT:
    {
      typedef std::set<Tercel::Triangle> TriangleSet;

      PAINTSTRUCT ps;
      HDC     hdc  = BeginPaint(hWnd, &ps);
      HPEN    hPen = CreatePen(PS_SOLID, 1, RGB(0, 0, 0));
      
      SelectObject(hdc, hPen);

      // ======================================  
      // Delaunay分割された三角形の集合を一つ一つ描画
      // ======================================  
      for(TriangleSet::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;
}

イテレータのお陰で、データ構造が変わっても Main.cpp 側のソースコードはほとんどそのままです。便利ですね。