C 语言的 Delaunay 三角剖分函数

• 下载演示项目 - 40 KB 
• 下载源代码 - 11 KB 

更新，2014年9月6日

当我创建这个项目时，我刚刚从Borland编译器切换到Visual Studio 2012，我通过修改DirectX SDK中的一个DXUT示例来创建了这个三角剖分项目。

此后我发现微软有一些奇怪的特点，比如在SDK中包含DXUT，却又建议我们不要使用它！

对您来说更重要的是，我最近在一个实际项目中使用了三角剖分函数，并且遇到了内存错误。显然，微软编译器中有太多的设置，导致三角剖分函数在某些项目中有效，但在其他项目中却失败，原因我无法理解。为了使三角剖分函数正常工作，我不得不更改这一行：

        xlm = ((simplex*) ( (char*)xlm + ((-1)) *simplex_size));

改为这样。

        xlm = ((simplex*) ( (char*)xlm - simplex_size ));

Ken Clarkson为1980年代的通用C编译器编写了该函数，而VS 2012并不总是按照他期望的方式解释((-1)) *simplex_size)。

我已将zip文件更新为包含这些更改。

介绍 

此函数接受一个二维或三维点数组，可以是整数值或浮点值，并使用Delaunay三角剖分算法创建一个索引列表，该列表可直接用于需要三角形索引列表的DirectX或OpenGL函数。

它可以非常快速地对数千个点进行三角剖分，以至于您可能可以在实时操作中使用它。它还提供了将三角形方向设置为顺时针或逆时针的选项。

背景

我最近开始学习DirectX，我发现每个形状都必须分解成三角形，并且所有三角形都需要被索引并按特定顺序排列，通常是顺时针。我在网上搜索了C源代码，希望能找到一个函数，我可以给它一个点数组，它就能返回一个三角形索引列表，但我找到的都是独立的程序，它们读取和写入文本文件，或者是用其他语言编写的函数。此外，这些程序似乎并不关心三角形是顺时针还是逆时针，而且教授提供的软件也有限制。

Ken Clarkson免费发布了他的凸包程序的源代码，所以我决定提取他的三角剖分功能并将其制成一个C函数。http://netlib.sandia.gov/voronoi/hull.html
 

使用代码 

• 三角剖分函数在一个文件中：Clarkson-Delaunay.cpp
• 它只有一个头文件：Clarkson-Delaunay.h
• 您只需调用一个函数即可使用它：BuildTriangleIndexList(...)

我将所有变量和函数设置为“static”，以便它们不会干扰您自己的软件中的名称。

为了测试该函数，我编写了test_triangulation.cpp来创建一个随机点数组，调用BuildTriangleIndexList()来处理这些点并生成一个三角形索引列表，然后将三角形显示在屏幕上。然后它重复这个过程。

它可以如此快速地对点进行三角剖分，以至于三角形会模糊不清。您必须按空格键并在每个周期之间进行单步操作才能看到结果。

我使用Visual Studio 2012创建了演示项目。由于我才刚开始学习DirectX，我 lấy了Microsoft DirectX SDK中的Tutorial02，并向其中添加了Clarkson-Delaunay.cpp和test_triangulation.cpp。三角剖分函数只是一个数学函数，但为了满足编译器的要求，我不得不将#include "DXUT.h"包含进去，尽管它并不需要DirectX头文件中的任何内容。我很难理解Video Studio和DirectX的头文件。我甚至无法编译微软的示例，就会收到诸如“C4005: 'DXGI_STATUS_OCCLUDED': macro redefinition”之类的警告。我提到这一点是为了说明，您可能需要更改Clarkson-Delaunay.cpp中的“include”文件才能使其在您的环境中编译。

三角剖分函数是：  WORD *BuildTriangleIndexList(void *, float, int *, int , int , int *);。

它接受六个参数，并返回一个指向三角形索引数组的指针：

WORD *BuildTriangleIndexList (
     void *pointList,          // INPUT, an array of either float or integer "points".
                                  // A "point" is either two values (an X and a Y),
                                  // or three values (XY and Z).
                                  // You must allocate memory and fill this array 
                                  // with XY or XYZ points.
     float factor,             // INPUT, if pointList is a list of integers, set this
                                  // parameter to zero to let the function realize that 
                                  // you are using integers. If you give this a value,
                                  // pointList will be interpreted as an array of 
                                  // floating-point values.
                                  // Clarkson's function works on integers, so if pointList 
                                  // is an array of floating-point values, each value will 
                                  // be multiplied by this factor in order to convert it to 
                                  // an integer. Therefore, provide a factor that is large 
                                  // enough so that when the floating points are converted,
                                  // you don't lose too many digits after the decimal point,
                                  // unless you want to lose some digits.
                                  // Example: if you provide a factor of only 2.0, then the
                                  // floating-point values 3.0001 and 3.499 will become the 
                                  // same integer value: 6
                                  // That would be acceptable if both of points are supposed
                                  // to be the same, but otherwise it could cause trouble.
     int numberOfInputPoints,  // INPUT, the number of points in the list,
                                  // not the number of bytes or values.
     int numDimensions,        // INPUT, 2 for X and Y points, or 3 for XY and Z points
     int clockwise             // There are three options:
                                  //     -1: put triangles in anti-clockwise order
                                  //      0: don't waste time ordering the triangles
                                  //      1: put triangles in clockwise order
     int *numTriangleVertices) // OUTPUT, this does not need to be initialized.
                                  //  BuildTriangleIndexList gives this a value.

该函数不需要初始化或关闭。但是，调用函数负责使用free()释放返回值。

您很可能需要至少4个变量来使用它。在我的演示程序中，我使用了这些：  

WORD *triangleIndexList;   // OUTPUT, this does not need initialization
int *testPointsXY;         // INPUT, I fill this with random points for testing
int g_numTestPoints;       // INPUT, the number of points.
                           // I set this to 16 at the top of the test_triangulation.cpp, 
                           // and you can adjust its value as the program is running by
                           // pressing i or d on the keyboard
int numTriangleVertices;   // OUTPUT, this does not need initialization

当我传递浮点值给它时，我调用它的方式如下：    

triangleIndexList =  BuildTriangleIndexList(
          (void*)testPointsXY,       // The array of points
          (float)MAX_RAND,          // Multiplication factor. For testing: max random value
          g_numTestPoints,         // The number of points
          2,                      // 2: the list is XY points, not XYZ points
          1,                     // 1, because I want the triangles clockwise
          &numTriangleVertices);

返回值由BuildTriangleIndexList()分配。

在演示项目zip文件中，有test_triangulation_integers.cpp，它与test_triangulation.cpp相同，只是它将一个整数数组传递给三角剖分函数，以防您想比较传递整数数组和浮点数组之间的差异。

BuildTriangleIndexList()的返回值是一个指向点列表的索引数组。您不需要对该数组做任何事情。只需将该变量和顶点数量放入设置三角形索引缓冲区的那段代码中，以及绘制三角形的语句中。例如，在我的演示程序中： 

bd.Usage = D3D11_USAGE_DEFAULT;
bd.ByteWidth = sizeof( WORD ) * numTriangleVertices;     // <---- from the function
bd.BindFlags = D3D11_BIND_INDEX_BUFFER;
<..>
InitData.pSysMem =  triangleIndexList;                       // <---- same here
pd3dDevice->CreateBuffer( &bd, &InitData, &g_pTestVectorIndexBuffer );
<..>

g_pImmediateContext->DrawIndexed( numTriangleVertices, 0, 0);  // <---- same here

另外，释放返回值

free ( triangleIndexList );

如何使用演示程序

窗口的标题栏显示了每个周期创建的点数以及已显示的周期数。

键盘命令

• 空格键：在周期之间单步操作，而不是全速循环
• i：增加每个周期的点数
• d：减少每个周期的点数
• t：在显示实心三角形或线框之间切换。当三角形是实心的时，您可以轻松地确定它们是否都是顺时针的。如果方向错误，它们将不会显示。
• ESC：退出
• 任何其他键：尽可能快地循环。

注释

输出是一个16位值的三角形索引列表，每个三角形需要三个整数，一个用于每个顶点。因此，为了分配输出内存，函数需要将三角形数量乘以3个WORD。

bytes needed for output array = (number_of_triangles * 3 * sizeof (WORD))

我不知道三角形的数量，但在我的测试中，三角形的数量通常是输入点的1.7到2倍之间，所以我假设三角形的数量永远不会超过输入点的3倍。因此，我为输出数组分配内存如下：

bytes needed for output array = (number_of_triangles * 3 * sizeof (WORD)) * 3

这就是为什么您在第976行看到这个

maxOutputEntries = numPointsProcessed * 3*3;

Clarkson提到他的程序由于尽可能使用精确算术而损失了2到3倍的性能。Clarkson将他的程序设计为精确的，这在您需要精确度时非常棒，但在计算仅短暂显示在屏幕上的图像的三角形时，这种精度是否必要？了解他函数的人可能能够修改它，使其接受一个标志来执行精确的任务，或者执行快速而粗糙的任务。或者，如果他的程序本质上是精确的，那么就需要开发另一个程序来执行快速而粗糙的任务。

static void triangleList_out (int v0, int v1, int v2, int v3) {...}

当您给Clarkson的程序一个二维点数组（X和Y值）时，它会返回一个包含三个值的数组，每个三角形一个，每个顶点的索引，这正是您所期望的。然而，当您提供三维点（XYZ值）时，它会返回第四个值。我不知道那个第四个值是什么，所以triangleList_out()忽略了它。

1. Clarkson的程序将MAXDIM定义为8，这为输入点设置了最多7个维度的限制。二维是X和Y，三维是XYZ，但一个点怎么可能有超过3个维度？我将MAXDIM减少到4，因为该变量在调用malloc时被使用，所以这应该会减少其内存需求。如果您使用的是四维或更多维度的点，您将需要更改MAXDIM。
2. 由于Clarkson的程序设计为读写磁盘文件，它并没有费心去计算会有多少个三角形。然而，作为一个函数，最好是知道三角形的数量，以便函数在为三角形索引列表分配内存时只需要调用一次malloc()。
3. 正如您在演示程序中看到的，它可以极快地确定三角形，至少当值是随机且只有几千个点时是如此。当一些点是网格格式时，例如1,2 1,3 1,4 1,5... 2,1, 2,2 2,3等，它的速度会明显变慢。
4. 虽然很少有人会尝试修改Clarkson的逻辑，但修改输入和输出函数很容易。输出函数在第1062行
5. 该函数返回一个指向16位三角形索引数组的指针。我不知道Clarkson的函数在输入点数量方面有什么限制，但我将返回值设置为16位整数，因为我假设没有人需要在一次函数调用中需要超过64,000个三角形。如果您想返回整数，您只需要在Clarkson-Delaunay.cpp中进行搜索和替换，将WORD更改为int。少数对WORD的引用是我添加的，而不是Clarkson添加的。他在这里有一些关于他的凸包程序的能力和限制的信息：http://netlib.sandia.gov/voronoi/hull.html。

关注点

Ken Clarkson的源代码中变量如此晦涩，宏如此复杂，让我想起了Obfuscated C Code Contest中的条目。他的STORAGE宏尤其有趣。这是一个宏：

#define STORAGE(X)						\
								\
size_t	X##_size;						\
X	*X##_list = 0;						\
								\
X *new_block_##X(int make_blocks)				\
{	int i;							\
	static	X *X##_block_table[max_blocks];			\
		X *xlm, *xbt;					\
	static int num_##X##_blocks;				\
/* long before;	*/						\
	if (make_blocks) {					\
/* DEBEXP(-10, num_##X##_blocks) */				\
		assert(num_##X##_blocks<max_blocks);		\
/* before = _memfree(0);*/					\
 DEB(0, before) DEBEXP(0, Nobj * ##X##_size)			\
								\
		xbt = X##_block_table[num_##X##_blocks++] =	\
			(X*)malloc(Nobj * ##X##_size);		\
 			memset(xbt,0,Nobj * ##X##_size);	\
/* DEB(0, after) DEBEXP(0, 8*_howbig((long*)xbt))*/		\
		if (!xbt) {					\
			DEBEXP(-10,num_##X##_blocks)		\
/* memstats(2);	*/						\
		}						\
		assert(xbt);					\
								\
		xlm = INCP(X,xbt,Nobj);				\
		for (i=0;i<Nobj; i++) {				\
			xlm = INCP(X,xlm,(-1));			\
			xlm->next = X##_list;			\
			X##_list = xlm;				\
		}						\
								\
		return X##_list;				\
	};							\
								\
	for (i=0; i<num_##X##_blocks; i++)			\
		free(X##_block_table[i]);			\
	num_##X##_blocks = 0;					\
	X##_list = 0;						\
	return 0;						\
}								\
								\
void free_##X##_storage(void) {new_block_##X(0);}		\

在你能解码那个宏之前，你必须先解码它里面的三个宏：INCP、DEB和DEBEXP

#define INCP(X,p,k) ((##X*) ( (char*)p + (k) * ##X##_size)) /* portability? */
#define DEB(ll,mes)  DEBS(ll) fprintf(DEBOUT,#mes "\n");fflush(DEBOUT); EDEBS
#define DEBEXP(ll,exp) DEBS(ll) fprintf(DEBOUT,#exp "=%G\n", (double) exp); fflush(DEBOUT); EDEBS 

但是等等！在你解码那三个宏之前，你必须解码它们里面的三个宏：DEBS、DEBOUT和EDEBS

#define DEBS(qq)  {if (DEBUG>qq) {
#define EDEBS }}
#define DEBOUT DFILE

最后，要完成这个过程，你需要知道DEBUG和DFILE是什么

#define DEBUG -7
extern FILE *DFILE;

这种宏让我庆幸我是一个人类而不是一个C编译器！

为了将他的程序转换成一个函数，我不得不弄清楚哪些内存没有被释放，哪些变量需要初始化，但我无法理解他的宏，所以我使用了我的Borland 2007编译器来展开他的宏。这就是为什么在我的这个软件变体中几乎没有宏的原因。