C++ で 3 次元 Delaunay 分割

某所にて、『2 次元の Delaunay 分割は簡単だから別にいらない。むしろ 3 次元 Delaunay 分割のソースコードが欲しい』という要望を間接的に頂いたので、みんなが大好きな C++ で作ってみましたよー (´-ω-`)

OpenProcessing に投稿した 3D Delaunay Triangulation よりも、若干ですがソースコードがブラッシュアップされているような気がします。

---

今回は 3D という事もあり、結果表示用のプログラムには DirectX 9 を使用しました。最新バージョンを使わなかったのは、Windows XP でも実行できるようにという配慮からです。あ、あくまで環境依存のプログラムは GUI 部分のみであり、アルゴリズム（Delaunay3d.h）は標準 C++ のみで実装しているため、他の環境にもそのまま移植できると思います。

使用方法は昨日の 2 次元 Delaunay 分割とほぼ同じで、以下の通りです（Main.cpp のコールバック関数、および render() 関数に実際のサンプルがあります）。

（※ただし、2 次元と 3 次元のモジュールを共存させようとすると、いたるところで名前が衝突するので注意）

1. 頂点セット（std::set<Tercel::Vector> オブジェクト）を宣言します。
2. 先ほどの頂点セットに、適当な頂点オブジェクトを格納します。
3. 三角形セット（std::set<Tercel::Triangle> オブジェクト）を宣言します。
4. Tercel::Delaunay3d::getDelaunayTriangles() 関数を呼びます。

このとき、getDelaunayTriangles() の第一引数に頂点セットの『参照』を、第2引数には三角形セットの『ポインタ』を渡します。

例によって、三角形（Tercel::Triangle）は、頂点オブジェクトのポインタを保持していますので、頂点セットから要素を削除した場合（あるいはメモリ上の位置が変わるような事が起きた場合）には、三角形の頂点は不正なポインタになってしまいます。

ご注意ください。

---

【Delaunay3d.h】

view plainprint?

1. // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-  
2. //  
3. // 3次元Delaunay分割  
4. //  
5. // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-  
6.   
7.   
8. #ifndef ___TERCEL_DELAUNAY3D  
9. #define ___TERCEL_DELAUNAY3D  
10.   
11. #include <algorithm>  
12. #include <cfloat>  
13. #include <cmath>  
14.   
15. #include <map>  
16. #include <set>  
17.   
18. namespace Tercel   
19. {  
20.   
21.   // ポインタの参照先を比較するためのファンクタ  
22.   class PtrComp  
23.   {  
24.   public:  
25.     template <class Type>   
26.     bool operator()(const Type* a, const Type* b) const  
27.     {  
28.       return *a < *b;  
29.     }  
30.   };  
31.   
32.   
33.   struct Vector {  
34.     double x, y, z;  
35.   
36.     // ======================================    
37.     // 等価性の判定  
38.     // ======================================   
39.     bool operator==(const Vector& v) const   
40.     {  
41.       return x == v.x && y == v.y && z == v.z;  
42.     }  
43.   
44.     // ======================================    
45.     // 大小判定（set / mapを構築する際に使用）  
46.     // ======================================   
47.     bool operator<(const Vector& v) const   
48.     {  
49.       return x != v.x ? x < v.x :  
50.              y != v.y ? y < v.y :  
51.                         z < v.z;  
52.     }  
53.   };  
54.   
55.   
56.   struct Circle   
57.   {  
58.     Vector center;  // 中心座標  
59.     double radius;  // 半径  
60.   };  
61.   
62.   
63.   class Triangle {  
64.   public:  
65.     const static size_t NUM_VERTEX = 3;  
66.     const Vector* p[NUM_VERTEX];  // 頂点座標  
67.       
68.     // ======================================    
69.     // 等価性の判定  
70.     // ======================================   
71.     bool operator==(const Triangle& t) const   
72.     {  
73.       const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];  
74.       for(int i = 0; i < NUM_VERTEX; ++i)  
75.       { // const Vector const * → const Vector *  
76.         p1[i] = const_cast<const Vector*>(p[i]);  
77.         p2[i] = const_cast<const Vector*>(t.p[i]);  
78.       }  
79.       std::sort(p1, p1 + NUM_VERTEX, PtrComp());  
80.       std::sort(p2, p2 + NUM_VERTEX, PtrComp());  
81.   
82.       for(int i = 0; i < NUM_VERTEX; ++i)  
83.         if(!(p1[i] == p2[i])) return false;  
84.       return true;  
85.     }  
86.   
87.     // ======================================    
88.     // 大小判定（set / mapを構築する際に使用）  
89.     // ======================================   
90.     bool operator<(const Triangle& t) const   
91.     {  
92.       const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];  
93.       for(int i = 0; i < NUM_VERTEX; ++i)  
94.       { // const Vector const * → const Vector *  
95.         p1[i] = const_cast<const Vector*>(p[i]);  
96.         p2[i] = const_cast<const Vector*>(t.p[i]);  
97.       }  
98.       std::sort(p1, p1 + NUM_VERTEX, PtrComp());  
99.       std::sort(p2, p2 + NUM_VERTEX, PtrComp());  
100.       for(int i = 0; i < NUM_VERTEX - 1; ++i)  
101.         if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];  
102.       return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];  
103.     }  
104.   };  
105.   
106.   
107.   class Tetrahedron {  
108.   public:  
109.     const static size_t NUM_VERTEX = 4;  
110.     const Vector* p[NUM_VERTEX];  
111.   
112.     // ======================================    
113.     // 等価性の判定  
114.     // ======================================   
115.     bool operator==(const Tetrahedron& t) const   
116.     {  
117.       const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];  
118.       for(int i = 0; i < NUM_VERTEX; ++i)  
119.       { // const Vector const * → const Vector *  
120.         p1[i] = const_cast<const Vector*>(p[i]);  
121.         p2[i] = const_cast<const Vector*>(t.p[i]);  
122.       }  
123.       std::sort(p1, p1 + NUM_VERTEX, PtrComp());  
124.       std::sort(p2, p2 + NUM_VERTEX, PtrComp());  
125.       for(int i = 0; i < NUM_VERTEX; ++i)  
126.         if(!(p1[i] == p2[i])) return false;  
127.       return true;  
128.     }  
129.   
130.     // ======================================    
131.     // 大小判定（set / mapを構築する際に使用）  
132.     // ======================================   
133.     bool operator<(const Tetrahedron& t) const   
134.     {  
135.       const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];  
136.       for(int i = 0; i < NUM_VERTEX; ++i)  
137.       { // const Vector const * → const Vector *  
138.         p1[i] = const_cast<const Vector*>( p[i] );  
139.         p2[i] = const_cast<const Vector*>(t.p[i]);  
140.       }  
141.       std::sort(p1, p1 + NUM_VERTEX, PtrComp());  
142.       std::sort(p2, p2 + NUM_VERTEX, PtrComp());  
143.   
144.       for(int i = 0; i < NUM_VERTEX - 1; ++i)   
145.         if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];  
146.       return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];  
147.     }  
148.   
149.     // ======================================    
150.     // 他の四面体と共有点を持つか    
151.     // ======================================   
152.     bool hasCommonPoints(const Tetrahedron& t) const  
153.     {  
154.       for(int i = 0; i < NUM_VERTEX; ++i)   
155.         for(int j = 0; j < NUM_VERTEX; ++j)  
156.           if(*p[i] == *t.p[j]) return true;  
157.       return false;  
158.     }  
159.   };  
160.     
161.   class Delaunay3d   
162.   {  
163.     // ======================================    
164.     // 型定義  
165.     // ======================================   
166.     typedef const std::set<Vector>         ConstVertexSet;  
167.     typedef ConstVertexSet::const_iterator ConstVertexIter;  
168.   
169.     typedef std::set<Triangle>             TriangleSet;  
170.   
171.     typedef std::set<Tetrahedron>          TetraSet;  
172.     typedef TetraSet::iterator             TetraIter;  
173.   
174.     typedef std::map<Tetrahedron, bool>    TetraMap;  
175.   
176.   
177.   private:  
178.     // ======================================    
179.     // LU分解による三元一次連立方程式の解法  
180.     // ======================================  
181.     static double lu(double a[3][3], int ip[3])   
182.     {  
183.       const int n = 3;  
184.       double weight[n];  
185.      
186.       for(int k = 0; k < n; ++k)  
187.       {  
188.         ip[k] = k;  
189.         double u = 0;  
190.         for(int j = 0; j < n; ++j)   
191.         {  
192.           double t = fabs(a[k][j]);  
193.           if (t > u) u = t;  
194.         }  
195.         if (u == 0) return 0;  
196.         weight[k] = 1/u;  
197.       }  
198.       double det = 1;  
199.       for(int k = 0; k < n; ++k)   
200.       {  
201.         double u = -1;  
202.         int m = 0;  
203.         for(int i = k; i < n; ++i)   
204.         {  
205.           int ii = ip[i];  
206.           double t = fabs(a[ii][k]) * weight[ii];  
207.           if(t>u)   
208.           {   
209.             u = t;   
210.             m = i;  
211.           }  
212.         }  
213.         int ik = ip[m];  
214.         if (m != k)   
215.         {  
216.           ip[m] = ip[k]; ip[k] = ik;  
217.           det = -det;  
218.         }  
219.         u = a[ik][k]; det *= u;  
220.         if (u == 0) return 0;  
221.         for (int i = k+1; i < n; ++i)   
222.         {  
223.           int ii = ip[i];  
224.           double t = (a[ii][k] /= u);  
225.           for(int j = k+1; j < n; ++j)   
226.             a[ii][j] -= t * a[ik][j];  
227.         }  
228.       }  
229.       return det;  
230.     }  
231.   
232.     static void solve(double a[3][3], double b[3], int ip[3], double x[3])   
233.     {  
234.       const int n = 3;  
235.   
236.       for(int i = 0; i < n; ++i)   
237.       {  
238.         int ii = ip[i];  
239.         double t = b[ii];  
240.         for (int j = 0; j < i; ++j)   
241.           t -= a[ii][j] * x[j];  
242.         x[i] = t;  
243.       }  
244.         
245.       for (int i = n-1; i >= 0; --i)   
246.       {  
247.         double t = x[i];  
248.         int ii = ip[i];  
249.         for(int j = i+1; j < n; ++j)   
250.           t -= a[ii][j] * x[j];  
251.         x[i] = t / a[ii][i];  
252.       }  
253.     }  
254.   
255.     static double gauss(double a[3][3], double b[3], double x[3])   
256.     {  
257.       const int n = 3;  
258.       int ip[n];  
259.       double det = lu(a, ip);  
260.      
261.       if(det != 0) solve(a, b, ip, x);  
262.       return det;  
263.     }  
264.   
265.     // ======================================    
266.     // 与えられた四面体の外接球を求める  
267.     // ======================================  
268.     static void getCircumcircle(const Tetrahedron& tetra, Circle *dst)   
269.     {  
270.       const Vector* p0 = tetra.p[0];  
271.       const Vector* p1 = tetra.p[1];  
272.       const Vector* p2 = tetra.p[2];  
273.       const Vector* p3 = tetra.p[3];  
274.   
275.       double A[3][3] = {  
276.         {p1->x - p0->x, p1->y - p0->y, p1->z - p0->z},  
277.         {p2->x - p0->x, p2->y - p0->y, p2->z - p0->z},  
278.         {p3->x - p0->x, p3->y - p0->y, p3->z - p0->z}  
279.       };  
280.   
281.       double b[3] = {  
282.         0.5 * (p1->x*p1->x - p0->x*p0->x + p1->y*p1->y - p0->y*p0->y + p1->z*p1->z - p0->z*p0->z),  
283.         0.5 * (p2->x*p2->x - p0->x*p0->x + p2->y*p2->y - p0->y*p0->y + p2->z*p2->z - p0->z*p0->z),  
284.         0.5 * (p3->x*p3->x - p0->x*p0->x + p3->y*p3->y - p0->y*p0->y + p3->z*p3->z - p0->z*p0->z)  
285.       };  
286.   
287.       double x[3] = {0.0, 0.0, 0.0};  
288.   
289.       if(gauss(A, b, x) == 0)  
290.       {  
291.         dst->center.x = dst->center.y = dst->center.z = 0;  
292.         dst->radius   = -1;  
293.       }   
294.       else   
295.       {  
296.         dst->center.x = x[0];  
297.         dst->center.y = x[1];  
298.         dst->center.z = x[2];  
299.   
300.         double dx = x[0] - p0->x;  
301.         double dy = x[1] - p0->y;  
302.         double dz = x[2] - p0->z;  
303.         dst->radius =  sqrt(dx * dx + dy * dy + dz * dz);  
304.       }  
305.     }  
306.       
307.     // ======================================    
308.     // 四面体の重複管理  
309.     // ======================================  
310.     static inline void addElementToRedundanciesMap(TetraMap* pRddcMap,    
311.         const Tetrahedron& t)     
312.     {    
313.       TetraMap::iterator it = pRddcMap->find(t);    
314.       if(it != pRddcMap->end() && it->second)     
315.       {    
316.         // 値を (t, true) から (t, false) に置換    
317.         pRddcMap->erase(it);    
318.         pRddcMap->insert(TetraMap::value_type(t, false));    
319.       }   
320.       else  
321.         pRddcMap->insert(TetraMap::value_type(t, true));    
322.     }  
323.   
324.   public:  
325.     static void getDelaunayTriangles(ConstVertexSet& pVertexSet,   
326.         TriangleSet* triangleSet)   
327.     {  
328.   
329.       TetraSet tetraSet;  
330.   
331.       // ======================================    
332.       // 巨大な外部四面体を作る    
333.       // ======================================  
334.       Vector p0, p1, p2, p3;  
335.       Tetrahedron hugeTetrahedron = {&p0, &p1, &p2, &p3};  
336.       //               !! 注意 !!  
337.       // --------------------------------------  
338.       // hugeTetrahedronの要素はローカル変数への  
339.       // ポインタなので、これらの頂点を共有して  
340.       // いるTetrahedronオブジェクトは、この関数  
341.       // を抜ける前に必ず除去しておく必要がある。  
342.       {  
343.         double maxX, maxY, maxZ;  
344.         double minX, minY, minZ;  
345.         maxX = maxY = maxZ = DBL_MIN;  
346.         minX = minY = minZ = DBL_MAX;  
347.           
348.         // まず、全ての頂点を包含する球を得る  
349.   
350.         // 中心座標を求める  
351.         for(ConstVertexIter it = pVertexSet.begin();  
352.           it != pVertexSet.end(); ++it)   
353.         {  
354.           if(maxX  < it->x) maxX = it->x;   
355.           if(it->x < minX ) minX = it->x;  
356.             
357.           if(maxY  < it->y) maxY = it->y;   
358.           if(it->y < minY ) minY = it->y;  
359.             
360.           if(maxZ  < it->z) maxZ = it->z;   
361.           if(it->z < minZ ) minZ = it->z;  
362.         }  
363.         double centerX = 0.5 * (maxX - minX);  
364.         double centerY = 0.5 * (maxY - minY);  
365.         double centerZ = 0.5 * (maxZ - minZ);  
366.   
367.         // 半径を求める  
368.         double dx      = centerX - minX;  
369.         double dy      = centerY - minY;  
370.         double dz      = centerZ - minZ;  
371.         double radius  = sqrt(dx * dx + dy * dy + dz * dz) + 0.1;  
372.   
373.         // 4つの頂点  
374.         p0.x = centerX;  
375.         p0.y = centerY + 3.0 * radius;  
376.         p0.z = centerZ;  
377.   
378.         p1.x = centerX - 2.0 * sqrt(2.0) * radius;  
379.         p1.y = centerY - radius;  
380.         p1.z = centerZ;  
381.   
382.         p2.x = centerX + sqrt(2.0) * radius;  
383.         p2.y = centerY - radius;  
384.         p2.z = centerZ + sqrt(6.0) * radius;  
385.   
386.         p3.x = centerX + sqrt(2.0) * radius;  
387.         p3.y = centerY - radius;  
388.         p3.z = centerZ - sqrt(6.0) * radius;  
389.       }  
390.   
391.       tetraSet.insert(hugeTetrahedron);  
392.   
393.       // --------------------------------------    
394.       // 点を逐次添加し、反復的に四面体分割を行う    
395.       // --------------------------------------    
396.       for(ConstVertexIter vIter = pVertexSet.begin();   
397.         vIter != pVertexSet.end(); ++vIter)   
398.       {  
399.         const Vector* p = &*vIter;  
400.   
401.         // --------------------------------------    
402.         // 追加候補の四面体を保持する一時マップ    
403.         // --------------------------------------  
404.         TetraMap rddcMap;  
405.   
406.         // --------------------------------------    
407.         // 現在の四面体セットから要素を一つずつ取り出して、      
408.         // 与えられた点が各々の四面体の外接球の中に含まれるかどうか判定      
409.         // --------------------------------------    
410.         for(TetraIter tIter = tetraSet.begin();     
411.           tIter != tetraSet.end();)   
412.         {  
413.           // 四面体セットから要素を取りだして…    
414.           Tetrahedron t = *tIter;  
415.   
416.           // 外接球を求める  
417.           Circle c;  
418.           getCircumcircle(t, &c);  
419.           double dx  = c.center.x - p->x;  
420.           double dy  = c.center.y - p->y;  
421.           double dz  = c.center.z - p->z;  
422.           double len = sqrt(dx * dx + dy * dy + dz * dz);  
423.   
424.           if(len < c.radius)  
425.           {  
426.             Tetrahedron t0 = {p, t.p[0], t.p[1], t.p[2]};  
427.             addElementToRedundanciesMap(&rddcMap, t0);  
428.   
429.             Tetrahedron t1 = {p, t.p[0], t.p[2], t.p[3]};  
430.             addElementToRedundanciesMap(&rddcMap, t1);  
431.   
432.             Tetrahedron t2 = {p, t.p[0], t.p[3], t.p[1]};  
433.             addElementToRedundanciesMap(&rddcMap, t2);  
434.   
435.             Tetrahedron t3 = {p, t.p[1], t.p[2], t.p[3]};  
436.             addElementToRedundanciesMap(&rddcMap, t3);  
437.   
438.             tetraSet.erase(tIter++);  
439.           } else ++tIter;  
440.         }  
441.   
442.         for(TetraMap::iterator iter = rddcMap.begin();  
443.           iter != rddcMap.end(); ++iter)   
444.         {  
445.           if(iter->second) tetraSet.insert(iter->first);  
446.         }  
447.       }  
448.   
449.       // --------------------------------------    
450.       // 最後に、外部三角形の頂点を削除    
451.       // --------------------------------------    
452.       for(TetraIter tIter = tetraSet.begin();     
453.         tIter != tetraSet.end(); )     
454.       {    
455.         if(hugeTetrahedron.hasCommonPoints(*tIter))    
456.            tetraSet.erase(tIter++);    
457.         else ++tIter;    
458.       }  
459.   
460.       // そして、伝説へ…  
461.       for(TetraIter tIter = tetraSet.begin();     
462.         tIter != tetraSet.end(); ++tIter)   
463.       {  
464.         Tetrahedron tetra = *tIter;  
465.   
466.         Triangle t0 = {tetra.p[0], tetra.p[1], tetra.p[2]};  
467.         triangleSet->insert(t0);  
468.   
469.         Triangle t1 = {tetra.p[0], tetra.p[2], tetra.p[3]};  
470.         triangleSet->insert(t1);  
471.   
472.         Triangle t2 = {tetra.p[0], tetra.p[3], tetra.p[1]};  
473.         triangleSet->insert(t2);  
474.   
475.         Triangle t3 = {tetra.p[1], tetra.p[2], tetra.p[3]};  
476.         triangleSet->insert(t3);  
477.       }  
478.     }  
479.   };  
480. }  
481.   
482. #endif  

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 3次元Delaunay分割
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


#ifndef ___TERCEL_DELAUNAY3D
#define ___TERCEL_DELAUNAY3D

#include <algorithm>
#include <cfloat>
#include <cmath>

#include <map>
#include <set>

namespace Tercel 
{

  // ポインタの参照先を比較するためのファンクタ
  class PtrComp
  {
  public:
    template <class Type> 
    bool operator()(const Type* a, const Type* b) const
    {
      return *a < *b;
    }
  };


  struct Vector {
    double x, y, z;

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

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


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


  class Triangle {
  public:
    const static size_t NUM_VERTEX = 3;
    const Vector* p[NUM_VERTEX];  // 頂点座標
    
    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Triangle& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());

      for(int i = 0; i < NUM_VERTEX; ++i)
        if(!(p1[i] == p2[i])) return false;
      return true;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Triangle& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());
      for(int i = 0; i < NUM_VERTEX - 1; ++i)
        if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];
      return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];
    }
  };


  class Tetrahedron {
  public:
    const static size_t NUM_VERTEX = 4;
    const Vector* p[NUM_VERTEX];

    // ======================================  
    // 等価性の判定
    // ====================================== 
    bool operator==(const Tetrahedron& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>(p[i]);
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());
      for(int i = 0; i < NUM_VERTEX; ++i)
        if(!(p1[i] == p2[i])) return false;
      return true;
    }

    // ======================================  
    // 大小判定（set / mapを構築する際に使用）
    // ====================================== 
    bool operator<(const Tetrahedron& t) const 
    {
      const Vector* p1[NUM_VERTEX], *p2[NUM_VERTEX];
      for(int i = 0; i < NUM_VERTEX; ++i)
      { // const Vector const * → const Vector *
        p1[i] = const_cast<const Vector*>( p[i] );
        p2[i] = const_cast<const Vector*>(t.p[i]);
      }
      std::sort(p1, p1 + NUM_VERTEX, PtrComp());
      std::sort(p2, p2 + NUM_VERTEX, PtrComp());

      for(int i = 0; i < NUM_VERTEX - 1; ++i) 
        if(!(*p1[i] == *p2[i])) return *p1[i] < *p2[i];
      return *p1[NUM_VERTEX - 1] < *p2[NUM_VERTEX - 1];
    }

    // ======================================  
    // 他の四面体と共有点を持つか  
    // ====================================== 
    bool hasCommonPoints(const Tetrahedron& t) const
    {
      for(int i = 0; i < NUM_VERTEX; ++i) 
        for(int j = 0; j < NUM_VERTEX; ++j)
          if(*p[i] == *t.p[j]) return true;
      return false;
    }
  };
  
  class Delaunay3d 
  {
    // ======================================  
    // 型定義
    // ====================================== 
    typedef const std::set<Vector>         ConstVertexSet;
    typedef ConstVertexSet::const_iterator ConstVertexIter;

    typedef std::set<Triangle>             TriangleSet;

    typedef std::set<Tetrahedron>          TetraSet;
    typedef TetraSet::iterator             TetraIter;

    typedef std::map<Tetrahedron, bool>    TetraMap;


  private:
    // ======================================  
    // LU分解による三元一次連立方程式の解法
    // ======================================
    static double lu(double a[3][3], int ip[3]) 
    {
      const int n = 3;
      double weight[n];
   
      for(int k = 0; k < n; ++k)
      {
        ip[k] = k;
        double u = 0;
        for(int j = 0; j < n; ++j) 
        {
          double t = fabs(a[k][j]);
          if (t > u) u = t;
        }
        if (u == 0) return 0;
        weight[k] = 1/u;
      }
      double det = 1;
      for(int k = 0; k < n; ++k) 
      {
        double u = -1;
        int m = 0;
        for(int i = k; i < n; ++i) 
        {
          int ii = ip[i];
          double t = fabs(a[ii][k]) * weight[ii];
          if(t>u) 
          { 
            u = t; 
            m = i;
          }
        }
        int ik = ip[m];
        if (m != k) 
        {
          ip[m] = ip[k]; ip[k] = ik;
          det = -det;
        }
        u = a[ik][k]; det *= u;
        if (u == 0) return 0;
        for (int i = k+1; i < n; ++i) 
        {
          int ii = ip[i];
          double t = (a[ii][k] /= u);
          for(int j = k+1; j < n; ++j) 
            a[ii][j] -= t * a[ik][j];
        }
      }
      return det;
    }

    static void solve(double a[3][3], double b[3], int ip[3], double x[3]) 
    {
      const int n = 3;

      for(int i = 0; i < n; ++i) 
      {
        int ii = ip[i];
        double t = b[ii];
        for (int j = 0; j < i; ++j) 
          t -= a[ii][j] * x[j];
        x[i] = t;
      }
      
      for (int i = n-1; i >= 0; --i) 
      {
        double t = x[i];
        int ii = ip[i];
        for(int j = i+1; j < n; ++j) 
          t -= a[ii][j] * x[j];
        x[i] = t / a[ii][i];
      }
    }

    static double gauss(double a[3][3], double b[3], double x[3]) 
    {
      const int n = 3;
      int ip[n];
      double det = lu(a, ip);
   
      if(det != 0) solve(a, b, ip, x);
      return det;
    }

    // ======================================  
    // 与えられた四面体の外接球を求める
    // ======================================
    static void getCircumcircle(const Tetrahedron& tetra, Circle *dst) 
    {
      const Vector* p0 = tetra.p[0];
      const Vector* p1 = tetra.p[1];
      const Vector* p2 = tetra.p[2];
      const Vector* p3 = tetra.p[3];

      double A[3][3] = {
        {p1->x - p0->x, p1->y - p0->y, p1->z - p0->z},
        {p2->x - p0->x, p2->y - p0->y, p2->z - p0->z},
        {p3->x - p0->x, p3->y - p0->y, p3->z - p0->z}
      };

      double b[3] = {
        0.5 * (p1->x*p1->x - p0->x*p0->x + p1->y*p1->y - p0->y*p0->y + p1->z*p1->z - p0->z*p0->z),
        0.5 * (p2->x*p2->x - p0->x*p0->x + p2->y*p2->y - p0->y*p0->y + p2->z*p2->z - p0->z*p0->z),
        0.5 * (p3->x*p3->x - p0->x*p0->x + p3->y*p3->y - p0->y*p0->y + p3->z*p3->z - p0->z*p0->z)
      };

      double x[3] = {0.0, 0.0, 0.0};

      if(gauss(A, b, x) == 0)
      {
        dst->center.x = dst->center.y = dst->center.z = 0;
        dst->radius   = -1;
      } 
      else 
      {
        dst->center.x = x[0];
        dst->center.y = x[1];
        dst->center.z = x[2];

        double dx = x[0] - p0->x;
        double dy = x[1] - p0->y;
        double dz = x[2] - p0->z;
        dst->radius =  sqrt(dx * dx + dy * dy + dz * dz);
      }
    }
    
    // ======================================  
    // 四面体の重複管理
    // ======================================
    static inline void addElementToRedundanciesMap(TetraMap* pRddcMap,  
        const Tetrahedron& t)   
    {  
      TetraMap::iterator it = pRddcMap->find(t);  
      if(it != pRddcMap->end() && it->second)   
      {  
        // 値を (t, true) から (t, false) に置換  
        pRddcMap->erase(it);  
        pRddcMap->insert(TetraMap::value_type(t, false));  
      } 
      else
        pRddcMap->insert(TetraMap::value_type(t, true));  
    }

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

      TetraSet tetraSet;

      // ======================================  
      // 巨大な外部四面体を作る  
      // ======================================
      Vector p0, p1, p2, p3;
      Tetrahedron hugeTetrahedron = {&p0, &p1, &p2, &p3};
      //               !! 注意 !!
      // --------------------------------------
      // hugeTetrahedronの要素はローカル変数への
      // ポインタなので、これらの頂点を共有して
      // いるTetrahedronオブジェクトは、この関数
      // を抜ける前に必ず除去しておく必要がある。
      {
        double maxX, maxY, maxZ;
        double minX, minY, minZ;
        maxX = maxY = maxZ = DBL_MIN;
        minX = minY = minZ = DBL_MAX;
        
        // まず、全ての頂点を包含する球を得る

        // 中心座標を求める
        for(ConstVertexIter it = pVertexSet.begin();
          it != pVertexSet.end(); ++it) 
        {
          if(maxX  < it->x) maxX = it->x; 
          if(it->x < minX ) minX = it->x;
          
          if(maxY  < it->y) maxY = it->y; 
          if(it->y < minY ) minY = it->y;
          
          if(maxZ  < it->z) maxZ = it->z; 
          if(it->z < minZ ) minZ = it->z;
        }
        double centerX = 0.5 * (maxX - minX);
        double centerY = 0.5 * (maxY - minY);
        double centerZ = 0.5 * (maxZ - minZ);

        // 半径を求める
        double dx      = centerX - minX;
        double dy      = centerY - minY;
        double dz      = centerZ - minZ;
        double radius  = sqrt(dx * dx + dy * dy + dz * dz) + 0.1;

        // 4つの頂点
        p0.x = centerX;
        p0.y = centerY + 3.0 * radius;
        p0.z = centerZ;

        p1.x = centerX - 2.0 * sqrt(2.0) * radius;
        p1.y = centerY - radius;
        p1.z = centerZ;

        p2.x = centerX + sqrt(2.0) * radius;
        p2.y = centerY - radius;
        p2.z = centerZ + sqrt(6.0) * radius;

        p3.x = centerX + sqrt(2.0) * radius;
        p3.y = centerY - radius;
        p3.z = centerZ - sqrt(6.0) * radius;
      }

      tetraSet.insert(hugeTetrahedron);

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

        // --------------------------------------  
        // 追加候補の四面体を保持する一時マップ  
        // --------------------------------------
        TetraMap rddcMap;

        // --------------------------------------  
        // 現在の四面体セットから要素を一つずつ取り出して、    
        // 与えられた点が各々の四面体の外接球の中に含まれるかどうか判定    
        // --------------------------------------  
        for(TetraIter tIter = tetraSet.begin();   
          tIter != tetraSet.end();) 
        {
          // 四面体セットから要素を取りだして…  
          Tetrahedron t = *tIter;

          // 外接球を求める
          Circle c;
          getCircumcircle(t, &c);
          double dx  = c.center.x - p->x;
          double dy  = c.center.y - p->y;
          double dz  = c.center.z - p->z;
          double len = sqrt(dx * dx + dy * dy + dz * dz);

          if(len < c.radius)
          {
            Tetrahedron t0 = {p, t.p[0], t.p[1], t.p[2]};
            addElementToRedundanciesMap(&rddcMap, t0);

            Tetrahedron t1 = {p, t.p[0], t.p[2], t.p[3]};
            addElementToRedundanciesMap(&rddcMap, t1);

            Tetrahedron t2 = {p, t.p[0], t.p[3], t.p[1]};
            addElementToRedundanciesMap(&rddcMap, t2);

            Tetrahedron t3 = {p, t.p[1], t.p[2], t.p[3]};
            addElementToRedundanciesMap(&rddcMap, t3);

            tetraSet.erase(tIter++);
          } else ++tIter;
        }

        for(TetraMap::iterator iter = rddcMap.begin();
          iter != rddcMap.end(); ++iter) 
        {
          if(iter->second) tetraSet.insert(iter->first);
        }
      }

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

      // そして、伝説へ…
      for(TetraIter tIter = tetraSet.begin();   
        tIter != tetraSet.end(); ++tIter) 
      {
        Tetrahedron tetra = *tIter;

        Triangle t0 = {tetra.p[0], tetra.p[1], tetra.p[2]};
        triangleSet->insert(t0);

        Triangle t1 = {tetra.p[0], tetra.p[2], tetra.p[3]};
        triangleSet->insert(t1);

        Triangle t2 = {tetra.p[0], tetra.p[3], tetra.p[1]};
        triangleSet->insert(t2);

        Triangle t3 = {tetra.p[1], tetra.p[2], tetra.p[3]};
        triangleSet->insert(t3);
      }
    }
  };
}

#endif

---

【Main.cpp】

view plainprint?

1. // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-  
2. //  
3. // 3次元Delaunay分割（表示用GUI）  
4. //  
5. // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-  
6. #include <windows.h>  
7. #include <d3d9.h>  
8. #include <d3dx9.h>  
9.   
10. #include "Delaunay3d.h"  
11.   
12. //------ライブラリの設定とか  
13. #pragma comment(lib, "d3dxof.lib")  
14. #pragma comment(lib, "dxguid.lib")  
15. #pragma comment(lib, "d3dx9.lib")  
16. #pragma comment(lib, "d3d9.lib")  
17. #pragma comment(lib, "winmm.lib")  
18.   
19. //------ウィンドウを出すために必要な物あれこれ  
20. const LPCTSTR APP_NAME    = TEXT("DXT2_BY_TERCEL");  
21. const LPCTSTR MUTEX_NAME  = TEXT("MUTEX_DXT2_BY_TERCEL");  
22. const LPCTSTR WINDOW_NAME = TEXT("3D Delaunay Triangulation");  
23.   
24. const UINT SCREEN_WIDTH  = 800;  
25. const UINT SCREEN_HEIGHT = 600;  
26. const UINT WINDOW_WIDTH  = (SCREEN_WIDTH +   
27.               GetSystemMetrics(SM_CXFIXEDFRAME) * 2);  
28. const UINT WINDOW_HEIGHT = (SCREEN_HEIGHT +   
29.               GetSystemMetrics(SM_CYFIXEDFRAME) * 2 +   
30.               GetSystemMetrics(SM_CYCAPTION));  
31.   
32. // ----------------------------------------  
33. // 頂点集合、およびそれらを基に  
34. // Delaunay 分割された三角形を格納しておくSTLコンテナ  
35. // ----------------------------------------  
36. std::set<Tercel::Vector>   vertices;  
37. std::set<Tercel::Triangle> triangles;  
38.   
39.   
40. //------Direct3Dに必要なものあれこれ  
41. LPDIRECT3D9 g_pD3D = NULL;  
42. LPDIRECT3DDEVICE9 g_pD3DDev = NULL;  
43. D3DPRESENT_PARAMETERS g_D3DPresentParam;  
44.   
45.   
46. //------関数プロトタイプ  
47. LRESULT CALLBACK WindowProc(HWND, UINT, WPARAM, LPARAM);  
48. HRESULT InitD3D(HWND);  
49. VOID    Render();  
50.   
51. //------FVF  
52. struct ColorPoint  
53. {  
54.   D3DXVECTOR3 vCoord; // 座標  
55.   DWORD   color;  // 色  
56. };  
57.   
58. //------エントリポイント  
59. int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,   
60.                    LPSTR lpCmdLine, int nCmdShow)  
61. {  
62.   MSG msg;  
63.   
64.   // 多重起動防止  
65.   HANDLE hMutex = CreateMutex(NULL, TRUE, MUTEX_NAME);  
66.   if(GetLastError() == ERROR_ALREADY_EXISTS)  
67.   {  
68.     // 多重起動を検知した際には、  
69.     // 既存のウィンドウを最前面に表示してプログラムを終了する  
70.     HWND existingWnd = FindWindow(APP_NAME, NULL);  
71.     if(existingWnd != NULL)  
72.     {  
73.       if(IsIconic(existingWnd))  
74.         ShowWindowAsync(existingWnd, SW_RESTORE);  
75.   
76.       SetForegroundWindow(existingWnd);  
77.     }  
78.     return -1;  
79.   }  
80.   
81.   // ウィンドウクラス登録   
82.   WNDCLASSEX wc;  
83.   
84.   wc.cbSize = sizeof(WNDCLASSEX);  
85.   wc.style = CS_HREDRAW | CS_VREDRAW;  
86.   wc.lpfnWndProc = WindowProc;  
87.   wc.cbClsExtra = 0;  
88.   wc.cbWndExtra = 0;  
89.   wc.hInstance = hInstance;  
90.   wc.hIcon = (HICON)LoadImage(NULL,  
91.     MAKEINTRESOURCE(IDI_APPLICATION),  
92.     IMAGE_ICON,  
93.     0,  
94.     0,  
95.     LR_DEFAULTSIZE | LR_SHARED);  
96.   wc.hCursor = (HCURSOR)LoadImage(NULL,  
97.     MAKEINTRESOURCE(IDC_ARROW),  
98.     IMAGE_CURSOR,  
99.     0,  
100.     0,  
101.     LR_DEFAULTSIZE | LR_SHARED);  
102.   wc.hbrBackground = NULL;  
103.   wc.lpszMenuName = NULL;  
104.   wc.lpszClassName = APP_NAME;  
105.   wc.hIconSm = (HICON)LoadImage(NULL,  
106.     MAKEINTRESOURCE(IDI_APPLICATION),  
107.     IMAGE_ICON,  
108.     0,  
109.     0,  
110.     LR_DEFAULTSIZE | LR_SHARED);  
111.   if(!RegisterClassEx(&wc)) return -1;  
112.   
113.   // --------------------------------------  
114.   // 固定サイズのウィンドウ作成  
115.   // --------------------------------------  
116.   HWND hWnd = CreateWindow(APP_NAME,  
117.     WINDOW_NAME,  
118.     WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX, // サイズ変更禁止  
119.     CW_USEDEFAULT,  
120.     CW_USEDEFAULT,  
121.     WINDOW_WIDTH,  
122.     WINDOW_HEIGHT,  
123.     NULL,  
124.     NULL,  
125.     hInstance,  
126.     NULL);  
127.   if(!hWnd) return -1;  
128.   ValidateRect(hWnd, NULL);  
129.   
130.   // --------------------------------------  
131.   // Direct3Dを初期化  
132.   // --------------------------------------  
133.   if(FAILED(InitD3D(hWnd))) return -1;  
134.   
135.   // --------------------------------------  
136.   // ウィンドウ表示  
137.   // --------------------------------------  
138.   ShowWindow(hWnd, nCmdShow);  
139.   UpdateWindow(hWnd);  
140.   
141.   while(TRUE)  
142.   {  
143.     if (PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))  
144.     {  
145.       if (msg.message == WM_QUIT) break;  
146.       TranslateMessage(&msg);  
147.       DispatchMessage(&msg);  
148.     }  
149.     else  
150.     {  
151.       Render();  
152.     }  
153.   }  
154.     
155.   if (g_pD3DDev) g_pD3DDev->Release();  
156.   if (g_pD3D)    g_pD3D->Release();  
157.   CloseHandle(hMutex);  
158.   return (int)msg.wParam;  
159. }  
160.   
161. //------Direct3Dを初期化します  
162. HRESULT InitD3D(HWND hWnd)  
163. {  
164.   if((g_pD3D = Direct3DCreate9(D3D_SDK_VERSION)) == NULL)  
165.     return E_FAIL;  
166.   
167.   D3DDISPLAYMODE d3ddm;  
168.   if(FAILED(g_pD3D->GetAdapterDisplayMode(D3DADAPTER_DEFAULT, &d3ddm)))  
169.     return E_FAIL;  
170.   
171.   ZeroMemory(&g_D3DPresentParam, sizeof(g_D3DPresentParam));  
172.   g_D3DPresentParam.BackBufferCount = 1;  
173.   g_D3DPresentParam.Windowed = TRUE;  
174.   g_D3DPresentParam.BackBufferFormat = D3DFMT_UNKNOWN;  
175.   g_D3DPresentParam.SwapEffect = D3DSWAPEFFECT_DISCARD;  
176.   g_D3DPresentParam.EnableAutoDepthStencil = TRUE;  
177.   g_D3DPresentParam.AutoDepthStencilFormat = D3DFMT_D16;  
178.   
179.   if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,  
180.     D3DDEVTYPE_HAL,  
181.     hWnd,  
182.     D3DCREATE_HARDWARE_VERTEXPROCESSING,  
183.     &g_D3DPresentParam,  
184.     &g_pD3DDev)))  
185.   {  
186.   
187.     if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,  
188.       D3DDEVTYPE_HAL,  
189.       hWnd,  
190.       D3DCREATE_SOFTWARE_VERTEXPROCESSING,  
191.       &g_D3DPresentParam,  
192.       &g_pD3DDev)))  
193.     {  
194.       if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,  
195.         D3DDEVTYPE_REF,  
196.         hWnd,  
197.         D3DCREATE_SOFTWARE_VERTEXPROCESSING,  
198.         &g_D3DPresentParam,  
199.         &g_pD3DDev)))  
200.       {  
201.         return E_FAIL;  
202.       }  
203.     }  
204.   }  
205.   
206.   return S_OK;  
207. }  
208.   
209. //------コールバック関数  
210. LRESULT CALLBACK WindowProc(HWND hWnd, UINT uMsg,   
211.               WPARAM wParam, LPARAM lParam)  
212. {  
213.   switch(uMsg)  
214.   {  
215.     case WM_DESTROY:  
216.       PostQuitMessage(0);  
217.       return 0;  
218.   
219.     case WM_CREATE:  
220.       // ======================================    
221.       // 頂点集合にランダムでデータをセット  
222.       // ======================================    
223.       for(int i = 0; i < 200; ++i) {  
224.         Tercel::Vector v; // ベクトルを宣言  
225.         float r   = 0.6f * SCREEN_WIDTH + 0.1f * (rand() % 10);  
226.         float phi = D3DXToRadian(rand() % 180 - 90);  
227.         float theta = D3DXToRadian(rand() % 360);  
228.   
229.         v.x = r*cos(phi)*cos(theta);  
230.         v.y = r*sin(phi);  
231.         v.z = r*cos(phi)*sin(theta);  
232.   
233.         vertices.insert(v);  
234.       }  
235.   
236.       // ======================================    
237.       // Delaunay分割を行う関数  
238.       // ======================================    
239.         
240.       // 結果は第2引数（std::set<Tercel::Triangle> 型）に格納される。  
241.       //   
242.       // Tercel::Triangle のメンバには、3つの頂点p1, p2, p3 があって、  
243.       // いずれも Tercel::Vertex オブジェクトへのポインタになっている  
244.       // これらのポインタは、vertices リストの要素を参照している。  
245.       Tercel::Delaunay3d::getDelaunayTriangles(vertices, &triangles);  
246.       //                                       ~~~~~~~~  ~~~~~~~~~~  
247.       //                                       頂点集合  三角形集合  
248.       //                                       （参照） （ポインタ）  
249.       return 0;  
250.   }  
251.   return DefWindowProc(hWnd, uMsg, wParam, lParam);  
252. }  
253.   
254. //------描画関数  
255. VOID Render() {  
256.   if(FAILED(g_pD3DDev->Clear(0, NULL,  
257.     D3DCLEAR_TARGET | D3DCLEAR_ZBUFFER,  
258.     D3DXCOLOR(0.0f, 0.0f, 0.0f, 1.0f),  
259.     1.0f, 0))) return;  
260.   // ワールド行列  
261.   D3DXMATRIX mWorld;  
262.   D3DXMatrixIdentity(&mWorld);  
263.   D3DXMatrixRotationY(&mWorld,   
264.                       D3DXToRadian(timeGetTime()/50.0f));  
265.   g_pD3DDev->SetTransform(D3DTS_WORLD, &mWorld);  
266.   
267.   // ビュー行列  
268.   const float PI = 3.1415926535f;  
269.   D3DXVECTOR3 vEyePoint(0.0f, 0.0f, 1500.0f);  
270.   D3DXVECTOR3 vLookAtPt(0.0f, 0.0f, 0.0f);  
271.   D3DXVECTOR3 vUp(0.0f, 1.0f, 0.0f);  
272.   
273.   D3DXMATRIX mView;  
274.   D3DXMatrixLookAtLH(&mView, &vEyePoint, &vLookAtPt, &vUp);  
275.   g_pD3DDev->SetTransform(D3DTS_VIEW, &mView);  
276.   
277.   float aspect;  
278.   aspect = (float)SCREEN_WIDTH / (float)SCREEN_HEIGHT;  
279.       
280.   // 射影行列  
281.   D3DXMATRIX mProj;  
282.   D3DXMatrixPerspectiveFovLH(&mProj,  
283.     D3DXToRadian(45.0f),  
284.     aspect,  
285.     1.0f,  
286.     10000.0f);  
287.   g_pD3DDev->SetTransform(D3DTS_PROJECTION, &mProj);  
288.   
289.   if(SUCCEEDED(g_pD3DDev->BeginScene()))  
290.   {  
291.     // ======================================    
292.     // ここから描画処理  
293.     // ======================================  
294.     ColorPoint vPoint[2];  
295.   
296.     vPoint[0].color = D3DCOLOR_ARGB(255, 0, 100, 200);  
297.     vPoint[1].color = D3DCOLOR_ARGB(255, 0, 100, 200);  
298.   
299.     g_pD3DDev->SetFVF(D3DFVF_XYZ | D3DFVF_DIFFUSE);  
300.     g_pD3DDev->SetRenderState(D3DRS_LIGHTING, FALSE);  
301.   
302.     // ======================================    
303.     // Delaunay分割された三角形の集合を一つ一つ取りだして描画  
304.     // ======================================    
305.     for(std::set<Tercel::Triangle>::iterator it = triangles.begin();  
306.       it != triangles.end(); ++it)   
307.     {  
308.       Tercel::Triangle t = *it;  // 三角形取得  
309.       for(int i = 0; i < 3; ++i) {  
310.         float x1 = (float)t.p[i]->x;  
311.         float y1 = (float)t.p[i]->y;  
312.         float z1 = (float)t.p[i]->z;  
313.   
314.         float x2 = (float)t.p[(i+1)%3]->x;  
315.         float y2 = (float)t.p[(i+1)%3]->y;  
316.         float z2 = (float)t.p[(i+1)%3]->z;  
317.         vPoint[0].vCoord = D3DXVECTOR3(x1, y1, z1);  
318.         vPoint[1].vCoord = D3DXVECTOR3(x2, y2, z2);  
319.         g_pD3DDev->DrawPrimitiveUP(D3DPT_LINELIST, 1,   
320.           vPoint, sizeof(ColorPoint));  
321.       }  
322.     }  
323.     g_pD3DDev->EndScene();  
324.   }  
325.   
326.   //----- デバイスロスト対策  
327.   if(g_pD3DDev->Present(NULL, NULL, NULL, NULL) == D3DERR_DEVICELOST)  
328.   {  
329.     if(g_pD3DDev->TestCooperativeLevel() == D3DERR_DEVICELOST)  
330.     {  
331.       return;  
332.     }  
333.     g_pD3DDev->Reset(&g_D3DPresentParam);  
334.   }  
335. }  

// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
//
// 3次元Delaunay分割（表示用GUI）
//
// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
#include <windows.h>
#include <d3d9.h>
#include <d3dx9.h>

#include "Delaunay3d.h"

//------ライブラリの設定とか
#pragma comment(lib, "d3dxof.lib")
#pragma comment(lib, "dxguid.lib")
#pragma comment(lib, "d3dx9.lib")
#pragma comment(lib, "d3d9.lib")
#pragma comment(lib, "winmm.lib")

//------ウィンドウを出すために必要な物あれこれ
const LPCTSTR APP_NAME    = TEXT("DXT2_BY_TERCEL");
const LPCTSTR MUTEX_NAME  = TEXT("MUTEX_DXT2_BY_TERCEL");
const LPCTSTR WINDOW_NAME = TEXT("3D Delaunay Triangulation");

const UINT SCREEN_WIDTH  = 800;
const UINT SCREEN_HEIGHT = 600;
const UINT WINDOW_WIDTH  = (SCREEN_WIDTH + 
              GetSystemMetrics(SM_CXFIXEDFRAME) * 2);
const UINT WINDOW_HEIGHT = (SCREEN_HEIGHT + 
              GetSystemMetrics(SM_CYFIXEDFRAME) * 2 + 
              GetSystemMetrics(SM_CYCAPTION));

// ----------------------------------------
// 頂点集合、およびそれらを基に
// Delaunay 分割された三角形を格納しておくSTLコンテナ
// ----------------------------------------
std::set<Tercel::Vector>   vertices;
std::set<Tercel::Triangle> triangles;


//------Direct3Dに必要なものあれこれ
LPDIRECT3D9 g_pD3D = NULL;
LPDIRECT3DDEVICE9 g_pD3DDev = NULL;
D3DPRESENT_PARAMETERS g_D3DPresentParam;


//------関数プロトタイプ
LRESULT CALLBACK WindowProc(HWND, UINT, WPARAM, LPARAM);
HRESULT InitD3D(HWND);
VOID    Render();

//------FVF
struct ColorPoint
{
  D3DXVECTOR3 vCoord; // 座標
  DWORD   color;  // 色
};

//------エントリポイント
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, 
                   LPSTR lpCmdLine, int nCmdShow)
{
  MSG msg;

  // 多重起動防止
  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 = 0;
  wc.hInstance = hInstance;
  wc.hIcon = (HICON)LoadImage(NULL,
    MAKEINTRESOURCE(IDI_APPLICATION),
    IMAGE_ICON,
    0,
    0,
    LR_DEFAULTSIZE | LR_SHARED);
  wc.hCursor = (HCURSOR)LoadImage(NULL,
    MAKEINTRESOURCE(IDC_ARROW),
    IMAGE_CURSOR,
    0,
    0,
    LR_DEFAULTSIZE | LR_SHARED);
  wc.hbrBackground = NULL;
  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;

  // --------------------------------------
  // 固定サイズのウィンドウ作成
  // --------------------------------------
  HWND hWnd = CreateWindow(APP_NAME,
    WINDOW_NAME,
    WS_OVERLAPPEDWINDOW^WS_THICKFRAME^WS_MAXIMIZEBOX, // サイズ変更禁止
    CW_USEDEFAULT,
    CW_USEDEFAULT,
    WINDOW_WIDTH,
    WINDOW_HEIGHT,
    NULL,
    NULL,
    hInstance,
    NULL);
  if(!hWnd) return -1;
  ValidateRect(hWnd, NULL);

  // --------------------------------------
  // Direct3Dを初期化
  // --------------------------------------
  if(FAILED(InitD3D(hWnd))) return -1;

  // --------------------------------------
  // ウィンドウ表示
  // --------------------------------------
  ShowWindow(hWnd, nCmdShow);
  UpdateWindow(hWnd);

  while(TRUE)
  {
    if (PeekMessage(&msg, NULL, 0U, 0U, PM_REMOVE))
    {
      if (msg.message == WM_QUIT) break;
      TranslateMessage(&msg);
      DispatchMessage(&msg);
    }
    else
    {
      Render();
    }
  }
  
  if (g_pD3DDev) g_pD3DDev->Release();
  if (g_pD3D)    g_pD3D->Release();
  CloseHandle(hMutex);
  return (int)msg.wParam;
}

//------Direct3Dを初期化します
HRESULT InitD3D(HWND hWnd)
{
  if((g_pD3D = Direct3DCreate9(D3D_SDK_VERSION)) == NULL)
    return E_FAIL;

  D3DDISPLAYMODE d3ddm;
  if(FAILED(g_pD3D->GetAdapterDisplayMode(D3DADAPTER_DEFAULT, &d3ddm)))
    return E_FAIL;

  ZeroMemory(&g_D3DPresentParam, sizeof(g_D3DPresentParam));
  g_D3DPresentParam.BackBufferCount = 1;
  g_D3DPresentParam.Windowed = TRUE;
  g_D3DPresentParam.BackBufferFormat = D3DFMT_UNKNOWN;
  g_D3DPresentParam.SwapEffect = D3DSWAPEFFECT_DISCARD;
  g_D3DPresentParam.EnableAutoDepthStencil = TRUE;
  g_D3DPresentParam.AutoDepthStencilFormat = D3DFMT_D16;

  if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
    D3DDEVTYPE_HAL,
    hWnd,
    D3DCREATE_HARDWARE_VERTEXPROCESSING,
    &g_D3DPresentParam,
    &g_pD3DDev)))
  {

    if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
      D3DDEVTYPE_HAL,
      hWnd,
      D3DCREATE_SOFTWARE_VERTEXPROCESSING,
      &g_D3DPresentParam,
      &g_pD3DDev)))
    {
      if(FAILED(g_pD3D->CreateDevice(D3DADAPTER_DEFAULT,
        D3DDEVTYPE_REF,
        hWnd,
        D3DCREATE_SOFTWARE_VERTEXPROCESSING,
        &g_D3DPresentParam,
        &g_pD3DDev)))
      {
        return E_FAIL;
      }
    }
  }

  return S_OK;
}

//------コールバック関数
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; // ベクトルを宣言
        float r   = 0.6f * SCREEN_WIDTH + 0.1f * (rand() % 10);
        float phi = D3DXToRadian(rand() % 180 - 90);
        float theta = D3DXToRadian(rand() % 360);

        v.x = r*cos(phi)*cos(theta);
        v.y = r*sin(phi);
        v.z = r*cos(phi)*sin(theta);

        vertices.insert(v);
      }

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

//------描画関数
VOID Render() {
  if(FAILED(g_pD3DDev->Clear(0, NULL,
    D3DCLEAR_TARGET | D3DCLEAR_ZBUFFER,
    D3DXCOLOR(0.0f, 0.0f, 0.0f, 1.0f),
    1.0f, 0))) return;
  // ワールド行列
  D3DXMATRIX mWorld;
  D3DXMatrixIdentity(&mWorld);
  D3DXMatrixRotationY(&mWorld, 
                      D3DXToRadian(timeGetTime()/50.0f));
  g_pD3DDev->SetTransform(D3DTS_WORLD, &mWorld);

  // ビュー行列
  const float PI = 3.1415926535f;
  D3DXVECTOR3 vEyePoint(0.0f, 0.0f, 1500.0f);
  D3DXVECTOR3 vLookAtPt(0.0f, 0.0f, 0.0f);
  D3DXVECTOR3 vUp(0.0f, 1.0f, 0.0f);

  D3DXMATRIX mView;
  D3DXMatrixLookAtLH(&mView, &vEyePoint, &vLookAtPt, &vUp);
  g_pD3DDev->SetTransform(D3DTS_VIEW, &mView);

  float aspect;
  aspect = (float)SCREEN_WIDTH / (float)SCREEN_HEIGHT;
    
  // 射影行列
  D3DXMATRIX mProj;
  D3DXMatrixPerspectiveFovLH(&mProj,
    D3DXToRadian(45.0f),
    aspect,
    1.0f,
    10000.0f);
  g_pD3DDev->SetTransform(D3DTS_PROJECTION, &mProj);

  if(SUCCEEDED(g_pD3DDev->BeginScene()))
  {
    // ======================================  
    // ここから描画処理
    // ======================================
    ColorPoint vPoint[2];

    vPoint[0].color = D3DCOLOR_ARGB(255, 0, 100, 200);
    vPoint[1].color = D3DCOLOR_ARGB(255, 0, 100, 200);

    g_pD3DDev->SetFVF(D3DFVF_XYZ | D3DFVF_DIFFUSE);
    g_pD3DDev->SetRenderState(D3DRS_LIGHTING, FALSE);

    // ======================================  
    // Delaunay分割された三角形の集合を一つ一つ取りだして描画
    // ======================================  
    for(std::set<Tercel::Triangle>::iterator it = triangles.begin();
      it != triangles.end(); ++it) 
    {
      Tercel::Triangle t = *it;  // 三角形取得
      for(int i = 0; i < 3; ++i) {
        float x1 = (float)t.p[i]->x;
        float y1 = (float)t.p[i]->y;
        float z1 = (float)t.p[i]->z;

        float x2 = (float)t.p[(i+1)%3]->x;
        float y2 = (float)t.p[(i+1)%3]->y;
        float z2 = (float)t.p[(i+1)%3]->z;
        vPoint[0].vCoord = D3DXVECTOR3(x1, y1, z1);
        vPoint[1].vCoord = D3DXVECTOR3(x2, y2, z2);
        g_pD3DDev->DrawPrimitiveUP(D3DPT_LINELIST, 1, 
          vPoint, sizeof(ColorPoint));
      }
    }
    g_pD3DDev->EndScene();
  }

  //----- デバイスロスト対策
  if(g_pD3DDev->Present(NULL, NULL, NULL, NULL) == D3DERR_DEVICELOST)
  {
    if(g_pD3DDev->TestCooperativeLevel() == D3DERR_DEVICELOST)
    {
      return;
    }
    g_pD3DDev->Reset(&g_D3DPresentParam);
  }
}

---

【おまけ】

勢いでリツイートしてしまったがコード付きという神記事だった2011年11月21日 21:10 via TweetDecknekomatu
nekomatu

わわっΣ(･Д･ﾉ)ﾉ 恐縮です（照