通过COM调用在ACAD里画图实在是慢，一个三角网中有几千个三角，等它全画出来，我都可以学唐伯虎去洗个澡再换身行头了。

昨天我疯狂了一把：设计了三角网“一笔画”算法，它会找出一个把三角网的所有边都穿过的路径，这样我只要画一根超级多段线就把整个三角网画出来了。

平常说的一笔画是不带回头的，但是大多数三角形网都不会符合一笔画的图论条件。实际上，在我的这个算法中，走回头路是标准做法。

这是基于深度优先搜索策略的递归算法，在递归的每一层中，它以当前节点为中心，每添加一条边（发散开去），就再添加一条反向的边（回到中心），等把以它为起点的边都添加到超级多段线以后，再依次以那些相邻节点为中心，做同样的事。

对了，在这篇日志——总结delaunay三角网程序中提到生成三角网的边列表时，用矩阵记录边的状态很耗内存，对于有几万个顶点的三角网，这个矩阵太大太稀疏了，现在我先全部添加，然后排序，最后扫描一遍列表就可以剔除重复项目，实际中这个做法快得多。

就数据量来说，这个超级多段线显然是臃肿的，但它真的很适合向ACAD里画，虽然比不上OBJECTARX或VBA的速度，但已经比原来快成百上千倍了。

PS. 调试递归程序时，跟踪总是很慢吗。根据昨天的体验，越深越慢，递归到20层时，按下F10，好久（几秒钟）才运行一行简单的代码，这对我的信念和耐心实在是一种考验。

搜索三角形的效率

在三角剖分阶段，我放弃了周期性地更新分区索引，每次搜索都从上次的结果开始，尝试了对点集进行预处理，主要是排序。效果分两种情形：

1、随机坐标点

速度一般，事先排个序的话，速度会有明显提高，在我的两次测试中，10k个点提高到了原来的2.5倍，50k个点提高到了原来的4倍多。

2、从ACAD等高线读取来的坐标点

由于这些点是关联的，也就是说索引相近的点，它们的坐标值一般也是相近的，所以增加一个排序处理后，总时间反而延长了。但是延长比例不太大，6385个点延长了15%，22750个点延长了20%。

那个排序，我是通过到原点的“曼哈顿距离”来比较点坐标的大小的，即X+Y，这显然不是准确的，但：

（1）我只需要在统计上，XY坐标值相邻的点在点集中的索引也相近；

（2）怎么准确呢?坐标是二维的，排序却是线性的，字典比较可以确定地比较两个坐标的大小，但在我这个问题上，它不比曼哈顿距离比较有明显优势。

比较之后，扫描一遍点集（我需要为点对象重新分配索引），可以顺便剔除一部分重复点，剩下的重复点仍然留到三角剖分阶段处理，在那里它们将被一个不漏地识别出来，并且标记为冗余（实际上我不是用布尔型，而是用一个整型的RedundantBy成员标记谁把它冗余了）。

调试

昨天下载了一个NATVIE C++单元测试工具，TUT(C++ Template Unit Test Framework),其实它就是三个头文件，使用及其方便。

至于那个传说中的WinUnit……我靠，MS的代码下载页混乱不堪，我点“打印版本”才看到了清爽的内容，但不能下载。CSDN下载也失效了。在GOOGLE上搜索WinUnit，铺天盖地地出来的是一个同名的用于单位转换的小工具（敢说那小玩意也挺强悍的），把它过滤后就剩4条结果，一个MS，一个CSDN，两个日文的……被这种情形折磨个大半天，一般人是会抓狂的。

今天用TUT给这个程序做了比较充实的单元测试，修改了不少逻辑错误——在纯几何上，它们是没问题的，但一旦计算机世界里有了浮点误差这东西，事情就复杂多了，感觉上，我的代码的80%是为处理容差而存在的。

莫非调试程序是件修身养性的事？

不管怎么说，现在的代码是彪悍多了。上回提到的诡异的bug，貌似已经消除了，但我不确定是哪一个改动生效的，总之这次在处理浮点数容差的问题上，比起前有进步，抽空我会总结一下。

PS 忙啊，昨天发热了，头疼欲裂，今天恢复了，起了个大早，看起来也睡了个大晚 :) good night to myself.

总结一下那个以delaunay三角网为核心的应用程序，程序本身还在完善中。

一、三角网的数据结构

1. 点集合；
2. 三角形集合，每个三角形由三个整数表示，代表三个顶点在点集合中的序号；
3. 可以查询每个顶点所关联的三角形；
4. 可以查询三角形的邻居，最多三个，这项功能由前两项派生；

二、按坐标搜索三角形

就是要在2D空间里查找指定点所在的三角形，如果点在三角形的边或顶点上，同时返回边或顶点的序号。

我采用线性搜索方法，搜索路径是一个大体上笔直的三角形带，以三角形为单位衡量，这个三角形带的宽度为1。这个算法的一个可选参数是起始三角形，如果未指定，就从三角形集合里的第一个开始。

这一算法对三角网的数据的“完整性约束”有很高要求，确切地说，查询点的关联三角形，和查询三角形的邻居必须准确无误。我曾经想在三角网数据结构中再添加边集合，因为控制不了数据状态的复杂性而放弃了，现在把边信息实现为行为（函数）而非数据。

昨天作了一项改进，通过建立三角网分区索引，优化对搜索起点的选择。用类似经纬网的网格把三角网划分成多个方块，那么搜索的第一步是确定指定点在哪一个方块中，第二步是以所在方块内的某个三角形为起点，执行线性搜索。

网格是规整的，所以第一步只需要几个简单的算术运算就可完成，再所以，不需要建立多级索引。分区索引的建立可以在对三角形集合的一次遍历中完成。总之，网格的密度对分区索引的建立和访问速度都没有影响。

以下是一些比较毛糙的实验数据，记录的是建立delaunay三角剖分的耗时，以ms为单位。每组有两个点集，大小分别为6385，22750。

---------------------------
无分区
---------------------------
1412/5307
1502/5558
1422/6129
1412/5418
---------------------------
8×8分区
---------------------------
911/3825
941/3915
951/3866
952/3916
932/3886
---------------------------
16×16分区
---------------------------
942/3515
921/3626
932/3655
951/3615
932/3675
---------------------------
32×32分区（偶尔在计算完成后崩溃，原因还不清楚）
---------------------------
942/3395
951/3415
965/3425

可以看出当点集尺寸很大（显然6385还不够大）时，分区密度几乎是越大越好。

这里面有个很重要的因素，就是刷新分区索引的频率。我的delaunay三角剖分是用增量法构建的，目前是每处理1/8的点集时，就刷新一次，没有测试其它频率。

说实话，这项改进措施的实际意义并不明显，因为不管是构建三角网，还是其它后续的操作，只要频繁查询位置的，那些位置点之前的关联都是很有价值的信息，所以，记住上次的搜索结果，作为当前搜索的起点往往会取得更好的效果：

---------------------------
记住上次位置
---------------------------
721/2915
751/2624
772/2624
721/2884
751/2594
731/2554
711/2614
741/2624

三、获取三角网的边集合

目前的实现思路是：遍列三角形->建立顶点矩阵->转化为列表。

边是无方向的，矩阵的作用是防止添加重复的边，不过它极其耗内存，即使用长度为n*(n-1)/2的valarray<bool>实现这个矩阵，对于22750个顶点，也需要60M的内存。在我的384M内存的机器上，大量的时间花在虚拟内存交换上，所以我在考虑省略状态矩阵，用插入法创建有序的边列表。

四、寻找两层三角网的交点

在半个月前第一个实现方案中，我以三角形为线索，以广度优先策略搜索。后来考虑到我只需要离散的交点，并不在乎它们的顺序，所以干脆遍历边集合，计算一个三角网中的每一条边与另一个三角网的交点。这个方法很野蛮，但很可靠，我还在寻找更快的方法。

计算线段与三角网比较麻烦，因为要查找线段与哪些三角形（理论上可以有很多）相交。

无论以三角形还是三角网的边为线索，都涉及一个问题：“线索”来自那个三角网？我想不到有意义的依据来比较两个三角网中，哪一个更适合用来遍历其中的元素，实际上，对调一下双方角色，寻找两遍才能保证一个不漏。

我说的"凸包三角网"指的是, 这个三角网的边缘是凸多边形. 在构建Delanuay三角剖分时, 这是一个可选的要求.

用增量算法构建的三角网的边缘不一定是点集的凸包, 为了填补边缘和凸包之前的空隙, 我试过两种方法:

(1) 用三角形带填充 edge 和 hull. edge 和 hull 是两根逆时针的闭合多边形. 使用两个顶点指针 i 和 j, 分别指示 edge 和 hull 上的当前顶点. 记 edge 的顶点集合为 E, hull 的为 H. 一个事实是, H 中的每个点都属于 E, 反之则不然.

开始时 E[i] 和 H[j] 是重复的, 然后两个指针逐渐递增, 速度一般并不相等, 也不均匀, 能够保证的是:

    (a) E[i] 到 H[j] 的距离比它到 H 中的任何其它点的距离都短, 反之亦然.

    (b) 两个指针同时走完一圈, 那时它们的值分别是 i + E.size() 和 j + H.size().

如果 E[i] == H[j], E[k] == H[j+1], (假设这里的[]运算符能自动对下标取模.) 则点集{ H[j], H[j+1], E[k-1], E[k-2], ... E[i+1] }是一个闭合多变形, 记为S. 要做到事就是用一些三角形填满这个多边形区域.

在我的实现中, 有一个一度让我很抓狂的 BUG: 有时候填充的三角形竟然会与原来的三角网重叠, 这直接导致我的三角网边缘搜索算法失败.这个 BUG 源于我先前在填补三角形没有遵守"最小外接圆"性质, 导致新创建的三角形内部有其它顶点落入. 现在看来, 合理的填充三角形的规则是:

    (a) 如果 E[i] == H[j], 就创建三角形 { H[j], E[i+1], E[i] }, 然后 i++.

    (b) 如果 E[i] != H[j], 就要么创建三角形 { H[j], H[j+1], E[i] } 然后 j++, 要么 { H[j], E[i+1], E[i] } 然后 i++, 取决于线段 <E[i], H[j+1]> 和 <H[j], E[i+1]> 的长度比较.

(2) 迭代三角网边缘, 遇到一个坑就用相应的三角形填上. 这个过程中, 三角网上的坑被逐渐填平. 在我的实现中, 三角网的边缘至少要迭代一遍, 最后一遍只是检查有没有未填平的坑. 在迭代过程中, 每填补一个坑, 边缘上的顶点就减少一个. 这种方法不需要事先寻找点集的凸包, 但是由于集合 E 在迭代过程中是要反复删除其中的项目的, 如果不想每填一个小坑(这个小坑可能是大坑的一部分)就退出当前迭代, 然后从头开始, 就得在迭代器上费点事了. 我使用 vector, 没用迭代器, 而是用下标, 所以没任何问题. 如果用 list, 会麻烦些.

P.S.

开始时, 我在这个项目中用杨师兄的代码做delaunay三角剖分, 很快就发现有缺陷, 手头只有编译了的 ActiveX 控件, 没有源代码. 就到网上找了别人的代码, 用了个把月才发现他的代码也有问题, 而且没有改造的空间. 不得已, 自己写了, 又几个月过去了, 我的代码修正了好几次, 但仍然不敢说它就无敌了, 只能说, 已经碰到的任何情形, 它都能正确快速地处理, 但谁知道以后会不会碰到新的麻烦!

一段正确无误的程序, 需要多少心血啊. 首饰有瑕疵, 可以作为残次品便宜卖. 可代码呢, 你甚至不能说"这段代码95%是正确的", 因为这讲不通, 崩溃就是崩溃, 跟机器没有什么好商量的. 一个由几万个三角形拼成的网, 就因为其中一个三角形没守规矩, 导致所有经过它的搜索线路都死翘翘, 进得去出不来, 导致程序的后续数据处理全部终止. 有时候我会想, 我喜欢敲代码的根本动机, 到底是写出有用的程序, 还是在机器里创建一个完美的秩序?

写好那个DTriNet程序后，一直耿耿于怀它的一个bug，虽然这个bug通过一些补救代码应付过去了，但仍觉得不是个事。

花了两天时间，重构了代码，减化了数据结构——以后实在有必要克制一下设计繁复精妙的数据结构的冲动，吃力不讨好。我爸会说，小鸡吃蚯蚓，自绕自脖子。

调试这种程序，心理考验似乎大于智力考验。F5时，祈祷不要出现亮着红叉的assert fail。多少次痛苦地闭上眼睛，再隐忍地睁开，继续回味bug，最大的动力是：揪出作梗的bug的瞬间是很惬意的！

无论如何，现在的程序挺快，也足够健壮：