[原创]：用vlisp生成三维地形的方法

wkai

在建筑、水利、道路、桥梁等工程设计中，经常会遇到需要由甲方提供的等高线图纸建立出三维数字模型的情况。虽然一些商业软件，如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

源码及示例：

snoopychen

wkai版主的这个地形程序很好
这个帖子应该是wkai版主的投稿文章吧。终于大概知道这个程序的思路了。
我觉得这个程序有一个挺重要的地方就在于求解等高线间点的高程，用过wkai版主的程序
在最里面一圈等高线的高程版主是怎么确定的呢？是靠外插么。

wkai

这个简单程序将所有最内圈内的标高点标高都等同为该等高线的高程。因此是平顶的山峰。
在 0.05版中不会出现这个现象。

snoopychen

现在用的是版主你的aut程序集中那个地形软件
在内圈的等高线内也能体现出非平顶，过度也很自然，这个觉得挺神奇的。
因为最内圈等高线内没有高程，所以所有的算法应该都是一个近似，不同的思路将得到不同的结果，请问wkai版主这个算法是否需要通过下部坡度的梯度大小来外插得到这个过度呢。
假如是这样的话，似乎又算法太复杂。

wkai

的确 0.05版是需要通过下部坡度的梯度大小来确定那内圈标高点的高程。

WHGF

用VLISP做出只有VC、VB等实现的三维效果真的不错。实用！！
下来试用后感觉是对于没有封闭的等高线处理上要由用户自行封闭，使用上稍有不便，不知能否修改，完善，三维曲面的拟合是用相互独立的矩形实现的数据量大，如果用三维线表示会更好，数据量也小，编辑修改都方便。
这是个人意见，谢谢你的好东东。

momokill

我仍不明白最内圈内的标高点标高都等同为该等高线的高程什么意识？谁能再解释一下“将所有最内圈内的标高点标高都等同为该等高线的高程”

wkai

最内圈内的点的标高点的高程默认取最内圈等高线的高程.
最外圈外的点的标高点的高程默认取最外圈等高线的高程.
这是只是这个程序设定而已。

WHGF

建议取消X和Y方向的相邻方形框，数据量会减少近1/2。

wkai

最初由 WHGF 发布
[B]建议取消X和Y方向的相邻方形框，数据量会减少近1/2。 [/B]

这样做相当于将网格单位扩大了一倍，模型精度会减少的、

hounengwei

启动的命令改为3dmap

zml84

等高线还是使用线性插值，如改进为旋转抛物线插值或多项式插值就更好了。
我们期待着。。
另外，在找当前点所在等高区间时，只要在该点周围一定范围就可以了，这样可以节省时间。

alien

这里能用三角网制作吗

fylinwater

算法上处理有点问题， 把整体地形降低了