nogirlfriend 发表于 2008-4-16 12:09:00

有点失望 居然没有经典的Delaunay三角形

<p>rt </p>

qjchen 发表于 2008-4-23 20:05:00

:)
对,其实计算几何是计算机时代的几何新应用,值得好好学习的,记得Highflybird版主号称编过一套计算几何的程序,不知道会不会共享一下呢:)
以前我在XDCAD发过一个帖子,曾经置顶,不过不知道现在哪里去了,遍找不到,找出记录,如下
Eachy版主说:
找到一个离散点生成三角网的资料,都是英文,看不懂啊,不过里面还有各种语言的实现代码
http://astronomy.swin.edu.au/%7Epbourke/modelling/triangulate/
我说:
多谢eachy版主介绍了这个站点
经过了艰苦的翻译,把Fortran翻译成Lisp语言,其中,深深的感受到数组与表的转换之苦,所以编辑起来有点罗里八索的感觉,以后得好好的来想一个好办
法。这个效果和网站上的basic和dephi范例是相同的。
感受,可能深入了解算法并从一开始就采用表的操作,会更好些,翻译代码是
件很痛苦的事情,调试起来特别麻烦。很想半途而废,坚持下来了也挺开心。
这个算法昨天在《计算几何》中看到了,应该还是通过voronoi图的增量算法
来构建delaunary三角形的,不过由于只是伪代码,算法也有些细节不同,因而没有来得及仔细思考,这也是下一步有空需要研究的问题。
autolisp代码如下,早期写的代码,比较乱,不过没有再重写了
效果如后,对程序简单修改即可变成选一堆点的三角化

;;; Purpose: To Generate the triangle by incresing point
;;; Version: 0.1
;;; Credit to Paul Bourke (pbourke@swin.edu.au) for the original C Program :))
;;; http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/
;;; The following codes are translate from C to AutoLisp by QJCHEN
;;; South China University of Technology
;;; Thanks : Eachy at xdcad.net introduce the source code pages
;;; and the STDLIB Function of Reini Urban at http://xarch.tu-graz.ac.at/autocad/stdlib/archive/
;;; 2006.06.30
(defun c:test (/ tpoints temp howmany ij p1 p2 p3)
(setq tpoints 1
vertex (givever)
triangle (givetri)
edges (giveedg)
)
(while (setq temp (getpoint))
    (setq vertex (qj-setnmth (nth 0 temp) tpoints 1 vertex))
    (setq vertex (qj-setnmth (nth 1 temp) tpoints 2 vertex))
    (if (> tpoints 2)
      (progn
(setq howmany (Triangulate tpoints))
      )
    )
    (setq tpoints (1+ tpoints))
    (setq ij 0)
    (command "redraw")
    (if (>= tpoints 4)
      (progn
(repeat howmany
   (setq ij (1+ ij))
   (setq p1 (nth (1- (nth 0 (nth (1- ij) triangle))) vertex))
   (setq p2 (nth (1- (nth 1 (nth (1- ij) triangle))) vertex))
   (setq p3 (nth (1- (nth 2 (nth (1- ij) triangle))) vertex))
   (grdraw p2 p1 1)
   (grdraw p1 p3 1)
   (grdraw p2 p3 1)
)
      )
    )                           ; (grdraw p1 p3 1)
                     ; (grdraw p2 p3 1)
                     ; (grdraw p3 p1 1)
)
)
;|The main function|;
(defun Triangulate (nvert / xmin ymin xmax ymax i dx dy xmid ymid complete
   ntri inc nedge i j Triangulate1
   )
(setq xmin (xofv vertex 1))
(setq ymin (yofv vertex 1))
(setq xmax xmin
ymax ymin
)
(setq i 2)
(while (<= i nvert)
    (if (< (xofv vertex i) xmin)
      (setq xmin (xofv vertex i))
    )
    (if (> (xofv vertex i) xmax)
      (setq xmax (xofv vertex i))
    )
    (if (< (yofv vertex i) ymin)
      (setq ymin (yofv vertex i))
    )
    (if (> (yofv vertex i) ymax)
      (setq ymax (yofv vertex i))
    )
    (setq i (1+ i))
)
(setq dx (- xmax xmin))
(setq dy (- ymax ymin))
(if (> dx dy)
    (setq dmax dx)
    (setq dmax dy)
)
(setq xmid (/ (+ xmax xmin) 2))
(setq ymid (/ (+ ymax ymin) 2))
(setq vertex (qj-setnmth (- xmid (* dmax 2)) (1+ nvert) 1 vertex))
(setq vertex (qj-setnmth (- ymid dmax) (1+ nvert) 2 vertex))
(setq vertex (qj-setnmth xmid (+ nvert 2) 1 vertex))
(setq vertex (qj-setnmth (+ ymid (* 2 dmax)) (+ nvert 2) 2 vertex))
(setq vertex (qj-setnmth (+ xmid (* 2 dmax)) (+ nvert 3) 1 vertex))
(setq vertex (qj-setnmth (- ymid dmax) (+ nvert 3) 2 vertex))
(setq triangle (qj-setnmth (+ nvert 1) 1 1 triangle))
(setq triangle (qj-setnmth (+ nvert 2) 1 2 triangle))
(setq triangle (qj-setnmth (+ nvert 3) 1 3 triangle))
(setq complete (append
   complete
   (list nil)
   )
)
(setq ntri 1);;;;;;;;;;;start loop i
(setq i 1)
(while (<= i nvert)
    (setq nedge 0);;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    (setq j 0
   temp (- 1)
    )
    (while (< temp ntri)
      (setq j (1+ j)
   temp j
      )
      (if (/= (nth (1- j) complete) T)
(progn
   (setq inc (InCircle1 (xofv vertex i) (yofv vertex i) (xof vertex triangle
            j 1
            )
          (yof vertex triangle j 1) (xof vertex
             triangle j 2
      ) (yof vertex
      triangle j 2
          ) (xof vertex
          triangle j
          3
            ) (yof vertex triangle
            j 3
            )
      )
   )
)
      )
      (if inc
(progn
   (setq edges (qj-setnmth (nth 0 (nth (1- j) triangle)) 1
      (+ nedge 1) edges
      )
   )
   (setq edges (qj-setnmth (nth 1 (nth (1- j) triangle)) 2
      (+ nedge 1) edges
      )
   )
   (setq edges (qj-setnmth (nth 1 (nth (1- j) triangle)) 1
      (+ nedge 2) edges
      )
   )
   (setq edges (qj-setnmth (nth 2 (nth (1- j) triangle)) 2
      (+ nedge 2) edges
      )
   )
   (setq edges (qj-setnmth (nth 2 (nth (1- j) triangle)) 1
      (+ nedge 3) edges
      )
   )
   (setq edges (qj-setnmth (nth 0 (nth (1- j) triangle)) 2
      (+ nedge 3) edges
      )
   )
   (setq Nedge (+ Nedge 3))
   (setq triangle (qj-setnmth ( triangle ntri 1) j 1 triangle))
   (setq triangle (qj-setnmth ( triangle ntri 2) j 2 triangle))
   (setq triangle (qj-setnmth ( triangle ntri 3) j 3 triangle))
   (setq complete (std-setnth (nth (1- ntri) complete) (1- j)
         complete
    )
   )
   (setq j (1- j)
temp j
   )
   (setq ntri (1- ntri))
)
      )
    );;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    (setq j 1)
    (while (<= j (1- Nedge))
      (if (and
   (/= ( edges 1 j) 0)
   (/= ( edges 2 j) 0)
   )
(progn
   (setq k (1+ j))
   (while (<= k Nedge)
   (if (and
    (/= ( edges 1 k) 0)
    (/= ( edges 2 k) 0)
)
       (if (= ( edges 1 j) ( edges 2 k))
(if (= ( edges 2 j) ( edges 1 k))
    (progn
      (setq edges (qj-setnmth 0 1 j edges))
      (setq edges (qj-setnmth 0 2 j edges))
      (setq edges (qj-setnmth 0 1 k edges))
      (setq edges (qj-setnmth 0 1 k edges))
    )
)
       )
   )
   (setq k (1+ k))
   )
)
      )
      (setq j (1+ j))
    );;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;
    (setq j 1)
    (while (<= j Nedge)
      (if (and
   (/= ( edges 1 j) 0)
   (/= ( edges 2 j) 0)
   )
(progn
   (setq ntri (1+ ntri))
   (setq triangle (qj-setnmth ( edges 1 j) ntri 1 triangle))
   (setq triangle (qj-setnmth ( edges 2 j) ntri 2 triangle))
   (setq triangle (qj-setnmth i ntri 3 triangle))
   (setq complete (std-setnth nil (1- ntri) complete))
)
      )
      (setq j (1+ j))
    );;;;;;;;;;;;;;;;;;;;;;;;;;;
    (setq i (1+ i))
);;;;;end of loop i
(setq i 0
temp (- 1)
)
(while (< temp ntri)
    (setq i (1+ i)
   temp i
    )
    (if (or
   (> ( triangle i 1) nvert)
   (> ( triangle i 2) nvert)
   (> ( triangle i 3) nvert)
)
      (progn
(setq triangle (qj-setnmth ( triangle ntri 1) i 1 triangle))
(setq triangle (qj-setnmth ( triangle ntri 2) i 2 triangle))
(setq triangle (qj-setnmth ( triangle ntri 3) i 3 triangle))
(setq i (1- i)
       temp i
)
(setq ntri (1- ntri))
      )
    )
)
(setq Triangulate1 ntri)
Triangulate1
)
;;; std&acute;ú&Igrave;&aelig;±í&Ouml;&ETH;&micro;&Uacute;n&cedil;&ouml;&Ocirc;&ordf;&Euml;&Oslash;&micro;&Auml;&ordm;&macr;&Ecirc;&yacute;
(defun std-%setnth (new i lst / fst len)
(cond
    ((minusp i)
      lst
    )
    ((> i (setq len (length lst)))
      lst
    )
    ((> i (/ len 2))
      (reverse (std-%setnth new (1- (- len i)) (reverse lst)))
    )
    (t
      (append
(progn
   (setq fst nil)         ; ; possible vl lsa compiler bug
   (repeat (rem i 4)
   (setq fst (cons (car lst) fst)
    lst (cdr lst)
   )
   )
   (repeat (/ i 4)
   (setq fst (cons (cadddr lst) (cons (caddr lst) (cons
         (cadr lst)
         (cons
               (car lst)
               fst
         )
          )
      )
      )
    lst (cddddr lst)
   )
   )
   (reverse fst)
)
(if (listp new)
   new
   (list new)
)                   ; v0.4001
(cdr lst)
      )
    )
)
)
(defun std-setnth (new i lst)
(std-%setnth (list new) i lst)
)
;;; &acute;ú&Igrave;&aelig;&para;&thorn;&Icirc;&not;±í&Ouml;&ETH;&micro;&Uacute;i&ETH;&ETH;&micro;&Uacute;j&Aacute;&ETH;&Ocirc;&ordf;&Euml;&Oslash;&micro;&Auml;&ordm;&macr;&Ecirc;&yacute;&pound;¨i&ordm;&Iacute;j&acute;&Oacute;1&iquest;&ordf;&Ecirc;&frac14;&pound;&copy;
(defun qj-setnmth (new i j lst / listb lista)
(setq listb lst)
(setq i (1- i))
(setq j (1- j))
(setq lista (nth i lst))
(setq lista (std-setnth new j lista))
(setq listb (std-setnth lista i listb))
listb
)
;;; &raquo;&ntilde;&Egrave;&iexcl;&Auml;&sup3;&cedil;&ouml;&Ecirc;&yacute;×é±í&micro;&Uacute;&frac14;&cedil;&Iuml;&icirc;&micro;&Uacute;&frac14;&cedil;&Iuml;&icirc;&micro;&Auml;&Ecirc;&yacute;&Ouml;&micro;
(defun (a n m / i)         ; n&Ecirc;&Ccedil;&ETH;&ETH;&pound;&not;m&Ecirc;&Ccedil;&Aacute;&ETH;
(setq i (nth (1- m) (nth (1- n) a)))
i
)
;;; &raquo;&ntilde;&Egrave;&iexcl;&Auml;&sup3;&cedil;&ouml;&micro;&yen;&Aacute;&ETH;&Ecirc;&yacute;×é&micro;&Uacute;&frac14;&cedil;&Iuml;&icirc;&micro;&Auml;&Ecirc;&yacute;&Ouml;&micro;
(defun (a n / i)               ; n&Ecirc;&Ccedil;&ETH;&ETH;&pound;&not;m&Ecirc;&Ccedil;&Aacute;&ETH;
(setq i (nth (1- n) a))
i
)
;|Vertex has the form '((x1 y1)(x2 y2)(x3 y3)(x4 y4)())
The function xofv is to get the x value of the i element,i start from 1|;
(defun xofv (vertex i / res)
(setq res (nth 0 (nth (- i 1) vertex)))
res
)
;|Vertex has the form '((x1 y1)(x2 y2)(x3 y3)(x4 y4)())
The function yofv is to get the y value of the i element,i start from 1|;
(defun yofv (vertex i / res)
(setq res (nth 1 (nth (- i 1) vertex)))
res
)
;|Lis has the form '(((x11 y11)(x12 y12)(x13 y13))((x21 y21)(x22 y22)(x23
y23))(()()()))
The function xof is to get the x value of the i,j element,i and j start from
1
and j is the outer sequence, and i is the inter sequence, total 3|;
(defun xof (lisa lisb j v123 / res1 res2 res)
(setq res1 (nth (1- j) lisb))
(setq res2 (nth (1- v123) res1))
(setq res3 (nth (1- res2) lisa))
(setq res (nth 0 res3))
res
)
;|Lis has the form '(((x11 y11)(x12 y12)(x13 y13))((x21 y21)(x22 y22)(x23
y23))(()()()))
The function xof is to get the y value of the i,j element,i and j start from
1
and j is the outer sequence, and i is the inter sequence, total 3|;
(defun yof (lisa lisb j v123 / res1 res2 res)
(setq res1 (nth (1- j) lisb))
(setq res2 (nth (1- v123) res1))
(setq res3 (nth (1- res2) lisa))
(setq res (nth 1 res3))
res
)
;(defun append1 (new n lis / res1 res2 res)
;
;(setq res1 (nth (1- n) lis))
;(setq res2 (append
;         res1
;         (list new)
;         )
;)
;(setq res (std-setnth res2 (1- n) lis))
;res
;)
;
;|Return TRUE if the point (xp,yp) lies inside the circumcircle
made up by points (x1,y1) (x2,y2) (x3,y3)
The circumcircle centre is returned in (xc,yc) and the radius r
NOTE: A point on the edge is inside the circumcircle|;
(defun InCircle1 (xp yp x1 y1 x2 y2 x3 y3 / InCircle eps mx2 my2 xc yc m1
       mx1 my1 m2 mx2 my2 dx dy rsqr r drsqr
   )
(setq eps 0.000001)
(setq InCircle nil)
(if (and
(< (abs (- y1 y2)) eps)
(< (abs (- y2 y3)) eps)
      )
    (alert "INCIRCUM - F - Points are coincident !!")
    (progn
      (cond
((< (abs (- y2 y1)) eps)
   (setq m2 (/ (- x2 x3) (- y3 y2)))
   (setq mx2 (/ (+ x2 x3) 2))
   (setq my2 (/ (+ y2 y3) 2))
   (setq xc (/ (+ x2 x1) 2))
   (setq yc (+ my2 (* m2 (- xc mx2))))
)
((< (abs (- y3 y2)) eps)
   (setq m1 (/ (- x1 x2) (- y2 y1)))
   (setq mx1 (/ (+ x1 x2) 2))
   (setq my1 (/ (+ y1 y2) 2))
   (setq xc (/ (+ x3 x2) 2))
   (setq yc (+ my1 (* m1 (- xc mx1))))
)
(T
   (setq m1 (/ (- x1 x2) (- y2 y1)))
   (setq m2 (/ (- x2 x3) (- y3 y2)))
   (setq mx1 (/ (+ x1 x2) 2))
   (setq mx2 (/ (+ x2 x3) 2))
   (setq my1 (/ (+ y1 y2) 2))
   (setq my2 (/ (+ y2 y3) 2))
   (setq xc (/ (- (+ (* m1 mx1) my2) my1 (* m2 mx2)) (- m1 m2)))
   (setq yc (+ my1 (* m1 (- xc mx1))))
)
      )
      (setq dx (- x2 xc))
      (setq dy (- y2 yc))
      (setq rsqr (+ (* dx dx) (* dy dy)))
      (setq r (sqrt rsqr))
      (setq dx (- xp xc))
      (setq dy (- yp yc))
      (setq drsqr (+ (* dx dx) (* dy dy)))
      (if (<= drsqr rsqr)
(setq InCircle T)
      )
    )
)
InCircle
)
;|Determines which side of a line the point (xp,yp) lies.
The line goes from (x1,y1) to (x2,y2)
Returns -1 for a point to the left
         0 for a point on the line
      +1 for a point to the right|;
(defun whichside (xp yp x1 y1 x2 y2 / equation)
(setq equation (- (* (- yp y1) (- x2 x1)) (* (- y2 y1) (- xp x1))))
(cond
    ((> equation 0)
      (setq whichside (- 0 1))
    )
    ((= equation 0)
      (setq whichside 0)
    )
    (T
      (setq whichside 1)
    )
)
whichside
)
(defun givetri (/ lis)
(repeat 200
    (setq lis (append
lis
(list (list nil nil nil))
       )
    )
)
lis
)
(defun givever (/ lis)
(repeat 200
    (setq lis (append
lis
(list (list nil nil))
       )
    )
)
lis
)
(defun giveedg (/ lis lis1 lis2)
(repeat 200
    (setq lis1 (append
   lis1
   (list nil)
      )
    )
)
(setq lis2 lis1)
(setq lis (append
       lis
       (list lis1)
   )
)
(setq lis (append
       lis
       (list lis2)
   )
)
lis
)

页: [1]
查看完整版本: 有点失望 居然没有经典的Delaunay三角形