三角网 插值

ndgxy

那太慢了呀，如果三角形上万个，甚至几十万个，那效率太低了呀
有没其他办法呀？

求指点

newer

ndgxy 发表于 2017-11-1 10:41

一点也不慢，你试下 xdrx_points_delaunay创建三角网，几十万的，也不到1秒就创建出来了，这还是创建实体不是搜索，搜索就更快了。
如果你只想通过点就得到，没有别的办法，除非你交互的时候不只给点，还点选下三角网的3D FACE，一点也不慢。遍历几十万实体的数据库，几毫秒。ACAD很多操作都是遍历数据库。如果你想只给一个点，还想改进想快点，那么就把屏幕内看见的三角网选出来，只遍历屏幕内的。

fools

加了注释，帮助你理解

;;                                                                            ;
;; triloc         by ymg  August 2013                                         ;
;;                                                                            ;
;; Locates triangle which encloses point p using Lawson's Walk.               ;
;;                                                                            ;
;; Given p a point, Returns Index in tl of triangle containing the point.     ;
;; If outside the triangulation Return is nil.                                ;
;;                                                                            ;
;; Point list pl, Neigbour list nl and Triangle list tl are defined           ;
;; outside this routine by calling program.                                   ;
;;                                                                            ;
;; Optimized Speed and re-organized code January 2014                         ;
;; Nice but get lost when triangulation is disjointed.                        ;
;; Modified to start on a random edge and operate in a CDT                    ;
;;                                                                            ;
;;pl-端点集，每条数据储存一个三维点
;;tl-三角形集，每条数据储存三角形三个顶点在端点集中的索引号
;;nl-相邻三角形集，每条数据是每个三角形的三条边与哪个三角形相邻
(defun triloc (p / e i notfound p1 p2 p3 x x1 x2 x3 y y1 y2 y3)
  (if (not tn)
    (setq tn (/ (length tl) 2)) ;_对半查找，提高效率
  )
  (setq        x (car p) ;_横坐标
        y (cadr p) ;_纵坐标
        notfound t ;_查找标记
  )
  (while (and notfound tn)
    (setq i  (nth tn tl) ;_定位一个三角形
          p1 (nth (car i) pl) ;_提取该三角形的三个顶点
          p2 (nth (cadr i) pl)
          p3 (nth (caddr i) pl)
          x1 (- (car p1) x) ;_三个端点到插值点的X,Y分量
          y1 (- (cadr p1) y)
          x2 (- (car p2) x)
          y2 (- (cadr p2) y)
          x3 (- (car p3) x)
          y3 (- (cadr p3) y)
          e  (fix (mrand 3)) ;_调用随机数提高查找效率，减少发生最不利查找的概率
    )
    ;;下述语句原理一致，即不断对半查找
    ;;原理：两点确定一条直线，一条无限延伸的直线相当于将平面一分为二
    ;;通过判断在哪条边的外侧后提取相邻三角形，则另一侧的三角形均不需要遍历
    ;;相当于算法复杂度从O(n)降为O(LogN)
    (cond ((= e 0)
           (cond ((minusp (- (* x1 y2) (* y1 x2))) (setq tn (car (nth tn nl))))
                 ((minusp (- (* x2 y3) (* y2 x3))) (setq tn (cadr (nth tn nl))))
                 ((minusp (- (* x3 y1) (* y3 x1))) (setq tn (caddr (nth tn nl))))
                 ((setq notfound nil))
           )
          )
          ((= e 1)
           (cond ((minusp (- (* x2 y3) (* y2 x3))) (setq tn (cadr (nth tn nl))))
                 ((minusp (- (* x3 y1) (* y3 x1))) (setq tn (caddr (nth tn nl))))
                 ((minusp (- (* x1 y2) (* y1 x2))) (setq tn (car (nth tn nl))))
                 ((setq notfound nil))
           )
          )
          ((= e 2)
           (cond ((minusp (- (* x3 y1) (* y3 x1))) (setq tn (caddr (nth tn nl))))
                 ((minusp (- (* x1 y2) (* y1 x2))) (setq tn (car (nth tn nl))))
                 ((minusp (- (* x2 y3) (* y2 x3))) (setq tn (cadr (nth tn nl))))
                 ((setq notfound nil))
           )
          )
    )
  )
  tn
)

ndgxy

感谢回答
请问p,tl,nl从哪里得到呀？TRIANGULATE吗？没看到呀

求指点，谢谢了

fools

p        c:demoZ -> (setq pnt (grread nil 13 0)) -> (setq p  (cadr pnt))
pl        c:tin -> (setq pl (remduppoint (sortxy pl) fuzz))
tl        triangulate -> (setq  tl (cons (cadddr tr) tl))
nl        (get_neighbour tl)        

ndgxy

谢谢了，我试试看

fools

不客气，triloc选用的算法很好 。我测试过10000点，20000个三角形，1000个插值测试点，找寻三角形的平均次数为124.7次，结果秒出，远远快于遍历20000个三角形。

ndgxy

我加载了“测试里的文件”，用gen命令算出100点测试，用tin命令绘制出三角网和等高线，
并计算出nl (setq nl(getneighbour pl tl))
但是插值好像不对
里边的命令cz是我自己写的
这是为什么呀？求指点
我在triangulate函数里，把tl改为全局变量，估计pl tl nl的参数有问题，不知道问题出在哪里
求指点，多谢了

fools

以下是基于原程序写的简单插值示例，如果想改写原程序尽量少用全局变量，建议用XRECORD保存

1. ;;gen->tin后
   
2. (defun c:tt (/  PTS RANGEX RANGEY TN TR XYZ )
   
3.   (setq        rangex 5000                        ; Extent in X for the points ;
   
4.         rangey (/ rangex 1.6)                ; Extent in Y / Golden Ratio ;
   
5.   )
   
6.   (and (not nl) (setq nl (get_neighbour tl))) ;_生成相邻三角形集
   
7.   (repeat 1000 (setq pts (cons (list (mrand rangex) (mrand rangey)) pts))) ;_产生1000个随机插值点
   
8.   (foreach pt pts
   
9.     (and (setq tn (triloc pt)) ;_找寻所在三角形
   
10.          (setq tr (mapcar '(lambda (a) (nth a pl)) (nth tn tl)) ;_三角形三顶点
    
11.                xyz  (getz pt (car tr) (cadr tr) (caddr tr)) ;_计算插值点高程
    
12.          )
    
13.          (entmakex (list '(0 . "TEXT")
    
14.                          (cons 1 (rtos (caddr xyz) 2 3))
    
15.                          (cons 10 xyz)
    
16.                          (CONS 72 4)
    
17.                          (cons 11 xyz)
    
18.                          (cons 40 10)
    
19.                    )
    
20.          ) ;_文本值为中心点高程
    
21.     )
    
22.   )
    
23. )
    

;;gen->tin后<br /> (defun c:tt (/  PTS RANGEX RANGEY TN TR XYZ )<br />   (setq        rangex 5000                        ; Extent in X for the points ;<br />         rangey (/ rangex 1.6)                ; Extent in Y / Golden Ratio ;<br />   )<br />   (and (not nl) (setq nl (get_neighbour tl))) ;_生成相邻三角形集<br />   (repeat 1000 (setq pts (cons (list (mrand rangex) (mrand rangey)) pts))) ;_产生1000个随机插值点<br />   (foreach pt pts<br />     (and (setq tn (triloc pt)) ;_找寻所在三角形<br />          (setq tr (mapcar '(lambda (a) (nth a pl)) (nth tn tl)) ;_三角形三顶点<br />                xyz  (getz pt (car tr) (cadr tr) (caddr tr)) ;_计算插值点高程<br />          )<br />          (entmakex (list '(0 . "TEXT")<br />                          (cons 1 (rtos (caddr xyz) 2 3))<br />                          (cons 10 xyz)<br />                          (CONS 72 4)<br />                          (cons 11 xyz)<br />                          (cons 40 10)<br />                    )<br />          ) ;_文本值为中心点高程<br />     )<br />   )<br /> )<br />

lrd1861

ndgxy

fools 发表于 2017-11-3 21:12

好的，谢谢了
我试试看。

xyp1964