用vlisp生成三维地形的方法

在建筑、水利、道路、桥梁等工程设计中，经常会遇到需要由甲方提供的等高线图纸建立出三维数字模型的情况。虽然一些商业软件，如SURFER可以做出很精细的地形模型，但其使用十分复杂，往往还都需要掌握相当专业的相关知识。而在AutoCAD中，虽然我们可以使用实心体的拉伸建立的三维地形模型，但是效果象梯田一样，在多数情况下都不能满足使用要求；而使用直纹曲面命令生成的地形效果较好，但是经常出现很难控制的地形扭曲，而且坡面过渡比较生硬。（图1）可以说在AutoCAD中一直没有有效的方法能够方便的完成这一常见任务。
针对这个问题，经过多次尝试，我们最终用VisualLISP完成一个比较简便的可以由等高线生成三维地形的程序，生成的三维地形过渡自然，而且模型通用性强，能够满足初步的工程应用。
下面是本程序由等高线生成地形的一些基本依据。（图2）
由于等高线之间的座标点的高程不能直接得到，只能依据该点与邻近等高线的位置关系大致确定。由此，我们设定对于等高线图中的任意一点(p)，如果它(p)正好位于某条等高线上，那么其高程就是该等高线的高程。否则，其高程则应大于包含该点的并与其最近一条等高线(a)的高程(ea)，同时小于不包含该点的并与其最近一条等高线(b)的高程(eb)。假设等高线之间的坡度是均匀的话，那么点p的高程(ep)则大体可以根据它与等高线a的距离da、与等高线b的距离db及等高线高程差，由以下公式计算得出。

对于以上的文字描述，我们可以用以下程序步骤实施。
1 测定点与各条等高线的位置包含关系。
如果测定点与等高线外侧一点的连线与这条等高线的交点数为奇数，该点则位于等高线的内部；如果交点数为偶数，则位于等高线外部。需要注意的是，这条连线必须与等高线处于相同的水平面中，否则，用于求交点的vla函数VLA-INTERSECTWITH无法求出交点。
下面函数用于获得测试点的与某条等高线的关系，包括嵌套关系、最近距离，并将它们合并成一个信息表返回。
(DEFUN INFO? (PL PT PT_OUT / DIST ELEV INTS PTS TMP_L)
(SETQ ELEV (LIST (VLA-GET-ELEVATION PL));_获取等高线的高程
PT (APPEND PT ELEV);_修改测试点的高程与等高线相同
DIST (DISTANCE (VLAX-CURVE-GETCLOSESTPOINTTO PL PT) PT)
;_测试点与等高线的最近距离
PT (VLAX-3D-POINT PT);_转换为VARIANT
PT_OUT (VLAX-3D-POINT (APPEND PT_OUT ELEV))
;_修改外部点的高程与等高线相同
TMP_L (VLA-ADDLINE *MSPACE* PT PT_OUT);_绘制一条连接两点的线段
INTS (VLA-INTERSECTWITH PL TMP_L 0);_求出交点数目
INTS (VLAX-VARIANT-VALUE INTS))
(VLA-DELETE TMP_L);_删除测试用线段
(IF (> (VLAX-SAFEARRAY-GET-U-BOUND INTS 1) 0)
(PROGN (SETQ INTS (VLAX-SAFEARRAY->LIST INTS)
PTS (REM (/ (LENGTH INTS) 3) 2)))
(SETQ PTS 0))
(IF (= PTS 0);_是否嵌套
(LIST NIL DIST PL)
(LIST T DIST PL)))
2 在包含和不包含测定点的两组等高线中分别找到与测定点最近的一条。下面的函数可以在表中找出第一项(距离)最小的项目。
(DEFUN NEAREST (LST / MIN_INDEX)
(SETQ MIN_INDEX (CAR (VL-SORT-I LST '(LAMBDA (A B) (< (CAR A) (CAR B))))))
(NTH MIN_INDEX LST))
3 通过上面的公式得到测定点的高程。将信息表传递到以下函数中可以计算出高程
(DEFUN GET_ELEV (LSTA LSTB / DA DB EA EB)
(SETQ DA (CAR LSTA) EA (VLA-GET-ELEVATION (CADR LSTA)))
(SETQ DB (CAR LSTB) EB (VLA-GET-ELEVATION (CADR LSTB)))
(+ EA (* (- EB EA) (/ DA (+ DA DB)))))
以上是取得测定点大致高程的关键函数，通过下面函数将它们打包成一个完整的功能核心。以等高线列表，测试点，外部点作为参数，可以得到测试点的高程。
(DEFUN GETZ (PLS PT_TEST PT_OUT / IN OUT DA DB EA EB IN_NEAR OUT_NEAR IN OUT PT_OUT PT_TEST TMP)
(MAPCAR '(LAMBDA (X)
(SETQ TMP (INFO? X PT_TEST PT_OUT))
(IF (CAR TMP);_根据嵌套关系将等高线分为两组IN和OUT
(SETQ IN (CONS (CDR TMP) IN))
(SETQ OUT (CONS (CDR TMP) OUT))
))
PLS)
(COND ((AND IN OUT)(GET_ELEV (NEAREST IN) (NEAREST OUT))) ;_计算测试点Z坐标
(IN (VLA-GET-ELEVATION (CADR (NEAREST IN))));_测试点在所有等高线内侧
(OUT (VLA-GET-ELEVATION (CADR (NEAREST OUT))))));_测试点在所有等高线侧
至此，我们可以结合AutoCAD中的3DMESH命令，通过GETZ函数测定每个节点的高程，从而得到等高线的三维模型。
(DEFUN C:3DMAP (/ SS N SS_LST CM CN DELTAX DELTAY NX NY PT_LD PT_RU TEST_PT X Y Z) (VL-LOAD-COM)
(SETQ *ACADOBJECT* (VLAX-GET-ACAD-OBJECT))
(SETQ *ACADDOCUMENT* (VLA-GET-ACTIVEDOCUMENT *ACADOBJECT*))
(SETQ *MSPACE* (VLA-GET-MODELSPACE *ACADDOCUMENT*))
(COMMAND "UNDO" "BE")
(SETQ CM 50 CN 50);_三维地形的网格密度
(PRINC "\N选择等高线(闭合的多义线)：")
(SETQ SS (SSGET '((0 . "*POLYLINE"))))
(REPEAT (SETQ N (SSLENGTH SS))
(SETQ SS_LST (CONS (VLAX-ENAME->VLA-OBJECT (SSNAME SS (SETQ N (1- N))))
SS_LST)));_将选择集转换为VLA物体列表
(SETQ PT_LD (GETPOINT "生成地形范围的左下角点："))
(SETQ PT_LD (LIST (CAR PT_LD) (CADR PT_LD)));_转换为2维坐标
(SETQ PT_RU (GETCORNER PT_LD "生成地形范围的右上角点："))
(OR PT_OUT
(SETQ PT_OUT (GETPOINT "\N任意选择等高线外部点：")
PT_OUT (LIST (CAR PT_OUT) (CADR PT_OUT))));_用于测试嵌套关系
(SETQ DELTAX (/ (CAR (MAPCAR '- PT_RU PT_LD)) (1- CM)))
(SETQ DELTAY (/ (CADR (MAPCAR '- PT_RU PT_LD)) (1- CN)))
(SETQ X 0 Y 0);_计数器
(COMMAND "3DMESH" CM CN) ;_启动3DMESH命令
(REPEAT CM ;_X方向循环
(SETQ Y 0
NX (* X DELTAX)
X (1+ X))
(REPEAT CN ;_Y方向循环
(SETQ NY (* Y DELTAY))
(SETQ Y (1+ Y))
(SETQ TEST_PT (MAPCAR '+ PT_LD (LIST NX NY)));_当前点的2维坐标
(SETQ Z (GETZ SS_LST TEST_PT PT_OUT));_取得当前点的Z坐标
(COMMAND "NON" (APPEND TEST_PT (LIST Z)))));_输入3DMESH节点3维坐标
(COMMAND "UNDO" "E"))
这段LISP代码虽然可以实现从等高线生成三维地形的功能，但是还需要针对算法进行优化，以加快运算速度和效率，克服现在平顶或者平底的情况。同时也可以增强用户界面易用性方面进行改进，限于篇幅这里就不再详细探讨了。

wkai 附带了这个的图片 (这个图片经过缩小处理，鼠标滚轮缩放图片或点击图片可以查看清晰效果) :