Author Topic: Triangulation (re-visited)  (Read 312758 times)

0 Members and 1 Guest are viewing this topic.

Didge

  • Bull Frog
  • Posts: 211
Triangulation (re-visited)
« on: March 13, 2006, 09:43:53 AM »
[ Link to the latest version in this thread: http://www.theswamp.org/index.php?topic=9042.msg555891#msg555891   updated 11/17/15]


The attached code creates Triangulated Irregular Networks (TIN), and generally it does a wonderful job of it. Big credits due to the original authors, Messrs  Paul Bourke & Daniele Piazza.

Unfortunately I've found a re-produceable bug, consequently I've spent the past weekend scratching my head with very little sleep.

When selecting the points (see 'Points.dwg' A2k)  all appears well with the TIN, ie no 3Dface edges overlap their neighbouring faces.

Now when I scale those same points down by 0.1 (from any arbitary base point) and re-run the code, I get un-wanted over-lapping of the 3Dfaces.

I suspect much of this coding would benefit from some VL enhancements, but the new Vl-sort was about as far as my vanilla experience would take me.

I would appreciate any suggestions.

Code: [Select]
;******************************************************************;
; TRIANGULATE - Lisp command to create a TIN from 3D points.       ;
; ===========                                                      ;
;                                                                  ;
; Written by Daniele Piazza, ADN member Mechanical Solution s.r.l. ;
; http://pdcode.com/code.htm                                       ;
;                                                                  ;
; Original C coding "Triangulate"  written by PAUL BOURKE          ;
; http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/     ;
;                                                                  ;
; This program triangulates an irregular set of points.            ;
; You can replace some code (sorting, list manipulation,...) with  ;
; VLisp functions to reduce the execution time.                    ;
;                                                                  ;
; This code is not seriously tested, if you find a bug...sorry!!   ;
; Goodbye, Daniele                                                 ;
;*******************************************************************
;
(defun C:TRIANGULATE (/ fuzzy nulllist ss1 ptlst nv supertriangle trianglelst i j k edgelst
                        circle pt flag perc)

 (setq OLDCMD (getvar "CMDECHO"))
 (setvar "CMDECHO" 0)
 (command ".UNDO" "GROUP")
 (setq OLDSNAP (getvar "OSMODE"))
 (setvar "OSMODE" 0)
 (setq fuzzy 1e-8)   ; tolerance in equality test
 (setq nulllist nil)
 
 (princ "\nSelect points...")
 (setq ss1 (ssget '((0 . "POINT"))))
 (setq start (getvar "date") THINK-CNT 0)            ; initiate timer & Progress Spinner Counter
 (setq ptlst (getptlist ss1))                        ; convert selection set to point list
 (setq ptlst (xsort ptlst))                          ; sort point list by X co-ordinate
 (setq nv (length ptlst))              ; number of points
 (setq supertriangle (findsupertriangle ptlst))      ; find super triangle
 (setq ptlst (append ptlst supertriangle))        ; append coordinates to the end of vertex list
 (setq trianglelst (list (list supertriangle nil)))  ; add supertriangle to the triangle list

 (setq i 0)
 (while   (< i nv)
   (THINKING (strcat "Processing TIN - " (itoa (/ (* i 100) nv)) "%    "))  ; update progress spinner
   (setq pt (nth i ptlst))
   (setq edgelst nil)                                ; initialize edge buffer
   (setq j 0)
   (while (and trianglelst (setq triangle (car (nth j trianglelst))))
     (setq flag T)
     (if (not (cadr (nth j trianglelst)))
       (progn
         (setq circle (getcircircumcircle triangle))       ; calculate circumcircle   
    (if (< (+ (caar circle) (cadr circle)) (car pt))  ; test point x and (pt) location
      (setq trianglelst (nth_subst j (list (car (nth j trianglelst)) T) trianglelst))
    )
    (if (isinside pt circle)
      (setq edgelst     (addtriangleedges triangle edgelst)
       trianglelst (nth_del j trianglelst)
       flag        nil
      )
    )
       ) ; end progn
     )   ; end if
     (if flag (setq j (1+ j)) )
   ) ; end while loop
   (setq edgelst (removedoublyedges edgelst fuzzy nulllist))           ; remove all doubly specified edges
   (setq trianglelst (addnewtriangles pt edgelst trianglelst))         ; form new triangles for current point
   (setq i (1+ i))                                                     ; get next vertex
 ) ; end while loop
 (setq trianglelst (purgetrianglelst trianglelst supertriangle fuzzy)) ; remove triangles with supertriangles edges

 (foreach triangle (mapcar 'car trianglelst)                           ; draw triangles
   (drawtriangle triangle)
 )

 (setvar "OSMODE" OLDSNAP)
 (setq OLDSNAP nil)
 (command ".UNDO" "END")
 (setq stop (getvar "date"))
 (princ (strcat "\r TIN Complete - Elapsed time: " (rtos (* 86400.0 (- stop start)) 2 2) " secs."))
 (setvar "CMDECHO" OLDCMD) 
 (princ)
)



;   XSORT - Original Shell Sort function replaced with VLISP sort (much quicker  ;
;                                                                                        ;
(defun XSORT ( PTLST /)
  (vl-sort PTLST (function (lambda (e1 e2) (< (car e1) (car e2)) ) ) )
)

;   NTH_DEL                              ;
;                                 ;
;   delete the n item in the list (by position, not by value!!)         ;
;                                 ;
;   Elimina l'oggetto che si trova nella posizione N della lista LST. L'utilizzo di   ;
;   funzioni ricorsive,oltre a non assicurare maggiore velocità, può creare problemi;
;   di overflow dello stack in caso di liste molto lunghe.            ;
(defun NTH_DEL (N LST / l)
 (repeat n
  (setq l (cons (car lst) l)
   lst (cdr lst)
  )
 )
 (append (reverse l)(cdr lst))
)

;   NTH_SUBST                           ;
;                                 ;
;   Replace the index element in the list with new element. This function is    ;
;   recursive this is not a great solution with a large amount of data.      ;
;                                 ;
(defun NTH_SUBST (index new Alist)
 (cond
  ((minusp index) Alist)
  ((zerop index)(cons new (cdr Alist)))
  (T (cons (car Alist)(nth_subst (1- index) new (cdr Alist))))
 )
)

;   GETPTLIST                           ;
;                                 ;
;   sset -> list (p1 p2 p3 ... pn)                     ;
;                                 ;
(defun GETPTLIST (ss1 / i pt ptlst)
 (if (not (zerop (sslength ss1)))
  (progn
   (setq i 0)
   (while
     (setq pt (ssname ss1 i))
     (setq ptlst (cons (cdr (assoc 10 (entget pt))) ptlst))
     (setq i (1+ i))
   )
  )
 )
 ptlst
)

;   FINDSUPERTRIANGLE                        ;
;                                 ;
;   Search the supertriangle that contain all points in the data set      ;
;                                 ;
(defun FINDSUPERTRIANGLE (ptlst / xmax xmin ymax ymin dx dy dmax xmid ymid
                   trx1 trx2 trx3 try1 try2 try3 trz1 trz2 trz3
          )
 (setq xmax (apply 'max (mapcar 'car ptlst))
       xmin (apply 'min (mapcar 'car ptlst))
       ymax (apply 'max (mapcar 'cadr ptlst))
       ymin (apply 'min (mapcar 'cadr ptlst))
       dx (- xmax xmin)
       dy (- ymax ymin)
       dmax (max dx dy)
       xmid (* (+ xmax xmin) 0.5)
       ymid (* (+ ymax ymin) 0.5)
       trx1 (- xmid (* dmax 2.0))
       try1 (- ymid dmax)
       trz1 0.0
       trx2 xmid
       try2 (+ ymid dmax)
       trz2 0.0
       trx3 (+ xmid (* dmax 2.0))
       try3 (- ymid dmax)
       trz3 0.0       
 )
 (list (list trx1 try1 trz1)
       (list trx2 try2 trz2)
       (list trx3 try3 trz3)
 )
)

;   GETCIRCIRCUMCIRCLE                        ;
;                                 ;
;   Calculate the circumcircle (center, radius) of triangle in input      ;
;                                 ;
(defun GETCIRCIRCUMCIRCLE (triangle / p1 p2 p3 p1x p2x p3x p1y p2y p3y d xc yc rad)
 (setq p1 (car triangle)
       p2 (cadr triangle)
       p3 (caddr triangle)
       p1x (car p1) p1y (cadr p1)
       p2x (car p2) p2y (cadr p2)
       p3x (car p3) p3y (cadr p3)
       d (* 2.0 (+ (* p1y p3x)
                   (* p2y p1x)
                   (- (* p2y p3x))
                   (- (* p1y p2x))
                   (- (* p3y p1x))
                   (* p3y p2x)
                )
         )
       xc (/ (+ (* p2y p1x p1x )
                (- (* p3y p1x p1x))
                (- (* p2y p2y p1y))
                (* p3y p3y p1y)
                (* p2x p2x p3y)
                (* p1y p1y p2y)
                (* p3x p3x p1y)
                (- (* p3y p3y p2y))
                (- (* p3x p3x p2y))
                (- (* p2x p2x p1y))
                (* p2y p2y p3y)
                (- (* p1y p1y p3y))
             )
             d
          )
       yc (/ (+ (* p1x p1x p3x)
                (* p1y p1y p3x)
                (* p2x p2x p1x)
                (- (* p2x p2x p3x))
                (* p2y p2y p1x)
                (- (* p2y p2y p3x))
                (- (* p1x p1x p2x))
                (- (* p1y p1y p2x))
                (- (* p3x p3x p1x))
                (* p3x p3x p2x)
                (- (* p3y p3y p1x))
                (* p3y p3y p2x)
             )
             d
          )
       rad (sqrt (+ (* (- p1x xc)(- p1x xc))
                    (* (- p1y yc)(- p1y yc))
                 )
           )
 )
 (list (list xc yc) rad)
)

;   ISINSIDE                           ;
;                                 ;
;   test if pt is inside a circle                     ;
;                                 ;
(defun ISINSIDE (pt circle)
 (setq ctr (car circle)
       rad (cadr circle)
 )
 (< (distance pt ctr) rad)
)

;   ADDTRIANGLEEDGES                        ;
;                                 ;
;   add triangle edges at the edge queue                  ;
;                                 ;
(defun ADDTRIANGLEEDGES (triangle edgelst)
 (append edgelst (list (list (car triangle)  (cadr triangle))
                       (list (cadr triangle) (caddr triangle))
                       (list (caddr triangle)(car triangle))
                 )
 )
)

;   DRAWTRIANGLE                           ;
;                                 ;
;   the fun side of the algorithm. Draw triangulation.            ;
;                                 ;
(defun DRAWTRIANGLE (triangle)
  (entmake (list (cons 0 "3DFACE") (cons 10 (car triangle))  (cons 11 (caddr triangle))
                                   (cons 12 (cadr triangle)) (cons 13 (cadr triangle))))
)

;   EQUALMEMBER                           ;
;                                 ;
;   Check if "item" is in "lista" or not by equality test. With real number the   ;
;   standard fuction "member" not work correctly.               ;
;                                 ;
(defun EQUALMEMBER (item lista fuzzy /)
 (apply 'or (mapcar '(lambda (x) (equal x item fuzzy)) lista))
)         

;   REMOVEDOUBLYEDGES                        ;
;                                 ;
;   Test the edge queue to remove duplicates (warning CW & CCW!)         ;
;                                 ;
(defun REMOVEDOUBLYEDGES (edgelst fuzzy nulllist /)
 (setq j 0)
 (while (< j (length edgelst))
  (setq k (1+ j))
  (while (< k (length edgelst))
   (if
    (or (and (equal (car  (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
        )
        (and (equal (car  (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
        )
    )
    (setq edgelst (nth_subst j nulllist edgelst)
          edgelst (nth_subst k nulllist edgelst)
    )
   )
   (setq k (1+ k))
  )
  (setq j (1+ j))
 )
 edgelst
)

;   ADDNEWTRIANGLES                           ;
;                                 ;
;   Add new triangle generated by pt to triangle list.            ;
;                                 ;
(defun ADDNEWTRIANGLES (pt edgelst trianglelst / j triangle )
 (setq j 0)
 (while (< j (length edgelst))
  (if (nth j edgelst)
   (setq triangle    (cons pt (nth j edgelst))
         trianglelst (cons (list triangle nil) trianglelst)
   )
  )
  (setq j (1+ j))
 )
 trianglelst
)

;   PURGETRIANGLELST                        ;
;                                 ;
;   replace all triangles that share a vertex with supertriangle         ;
;                                 ;
(defun PURGETRIANGLELST (trianglelst supertriangle fuzzy /)
 (setq j 0)
 (while (and trianglelst (setq triangle (car (nth j trianglelst))))
  (if (apply 'or
             (mapcar '(lambda (x) (equalmember x supertriangle fuzzy))
                     triangle
             )
      )
   (setq trianglelst (nth_del j trianglelst))
   (setq j (1+ j))
  )
 )
)


;                                       ;
; THINKING - STANDARD PROGRESS SPINNER  ;
;                                       ;
(defun THINKING (prmpt)
  (setq THINK-CNT (1+ THINK-CNT))
  (princ (strcat "\r" (nth (rem THINK-CNT 4) '("\|" "\/" "\-" "\\")) prmpt))
)


; ********************************* END OF CODING *******************************************
(princ "\n'TRIANGULATE' Loaded \n")
(princ)
« Last Edit: November 18, 2015, 09:46:25 PM by CAB »
Think Slow......

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
Re: Triangulation (re-visited)
« Reply #1 on: March 13, 2006, 10:34:40 AM »
Running the lisp with drawing as is. Then undo, scaling by 0.1, running again, I did not see the edge
overlap. But I am not into 3D so perhaps I did not know what to look for. Is the overlap obvious?
As you can see here examining 2 triangles, the first points are identical as are the third & second
points respectively.

Using plain ACAD2000.


Code: [Select]
Command: list
1 found
                  3D FACE   Layer: "GS-ELEV"
                            Space: Model space
                   Handle = AB8
             first point, X=539837.13728  Y=293829.99543  Z=  2.82000
            second point, X=539828.96328  Y=293844.62043  Z=  2.82000
             third point, X=539826.73328  Y=293843.57243  Z=  2.87000
            fourth point, X=539826.73328  Y=293843.57243  Z=  2.87000


Command:
Command: list
1 found
                  3D FACE   Layer: "GS-ELEV"
                            Space: Model space
                   Handle = AB9
             first point, X=539837.13728  Y=293829.99543  Z=  2.82000
            second point, X=539826.73328  Y=293843.57243  Z=  2.87000
             third point, X=539835.71028  Y=293829.00243  Z=  2.87000
            fourth point, X=539835.71028  Y=293829.00243  Z=  2.87000

Oh, it took .99 seconds of the first run & 1.25 sec on the second run.
« Last Edit: March 13, 2006, 10:40:19 AM by CAB »
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

Didge

  • Bull Frog
  • Posts: 211
Re: Triangulation (re-visited)
« Reply #2 on: March 13, 2006, 11:17:57 AM »
Thank-you for taking a look at this one CAB, I've attached a ".jpg" to demonstrate the bug.

The leftmost image is a correct triangulation, the bug appears in the right-hand image as over-lapping magenta edges to the 3Dfaces.

As a temporary fix, I tried scaling-up the co-ordinates in the "getpntlist" function and then scaling them down again in the "drawtriangle" function - sadly this didn't work. The equality "Fuzzy" factor also grabbed my attention but my tweaks failed to correct the issue.



 

Think Slow......

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
Re: Triangulation (re-visited)
« Reply #3 on: March 13, 2006, 11:32:14 AM »
How about posting the scaled drawing.
I assume you posted the drawing before the scaling operation.
I did not get that error.
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

Didge

  • Bull Frog
  • Posts: 211
Re: Triangulation (re-visited)
« Reply #4 on: March 13, 2006, 11:40:15 AM »
As requested, they're the same points but scaled down by 0.1

I've tried these points on 3 different PC's and had the same results on each.

Think Slow......

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
Re: Triangulation (re-visited)
« Reply #5 on: March 13, 2006, 11:52:31 AM »
Hold on, I do see the problem now.
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

Didge

  • Bull Frog
  • Posts: 211
Re: Triangulation (re-visited)
« Reply #6 on: March 13, 2006, 12:07:23 PM »
I'm not sure if thats good or bad news, lol, I was beginning to hope it was something as basic as drawing units.

I noticed various versions of the code on Paul Bourke's site, most languages are represented but alas no lisp version as yet.

When used with the GET-Z function (in previous thread) we have the basis of an open-source 3D ground modelling package in the making.
( "Swamp Desktop" maybe :-)
Think Slow......

LE

  • Guest
Re: Triangulation (re-visited)
« Reply #7 on: March 13, 2006, 12:15:24 PM »
Didge;

I have used this library:

http://www.cs.cmu.edu/~quake/triangle.html

But it is to be use with your C++ solutions.   :-(

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
Re: Triangulation (re-visited)
« Reply #8 on: March 13, 2006, 01:45:14 PM »
Looks like there may be a problem with the GETCIRCIRCUMCIRCLE function.
Not sure though. See this picture, the one on the left has an error where the points in yellow are
the triangle & the circle should pass through those points as it does on the right frame.
I'm out of time so maybe someone else has some help for you.
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
Re: Triangulation (re-visited)
« Reply #9 on: March 13, 2006, 05:05:03 PM »
Test this:

Code: [Select]
;******************************************************************;
; TRIANGULATE - Lisp command to create a TIN from 3D points.       ;
; ===========                                                      ;
;                                                                  ;
; Written by Daniele Piazza, ADN member Mechanical Solution s.r.l. ;
; http://pdcode.com/code.htm                                       ;
;                                                                  ;
; Original C coding "Triangulate"  written by PAUL BOURKE          ;
; http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/     ;
;                                                                  ;
; This program triangulates an irregular set of points.            ;
; You can replace some code (sorting, list manipulation,...) with  ;
; VLisp functions to reduce the execution time.                    ;
;                                                                  ;
; This code is not seriously tested, if you find a bug...sorry!!   ;
; Goodbye, Daniele                                                 ;
;*******************************************************************
;
;;
;;  Changes by CAB 03/13/06
;;  replaced the GETCIRCIRCUMCIRCLE routine
;;
(defun C:TRIANGULATE (/ fuzzy nulllist ss1 ptlst nv supertriangle trianglelst i j k edgelst
                        circle pt flag perc)

(setq OLDCMD (getvar "CMDECHO"))
(setvar "CMDECHO" 0)
(command ".UNDO" "GROUP")
(setq OLDSNAP (getvar "OSMODE"))
(setvar "OSMODE" 0)
(setq fuzzy 1e-8)   ; tolerance in equality test
(setq nulllist nil)
 
(princ "\nSelect points...")
(setq ss1 (ssget '((0 . "POINT"))))
(setq start (getvar "date") THINK-CNT 0)            ; initiate timer & Progress Spinner Counter
(setq ptlst (getptlist ss1))                        ; convert selection set to point list
(setq ptlst (xsort ptlst))                          ; sort point list by X co-ordinate
(setq nv (length ptlst))      ; number of points
(setq supertriangle (findsupertriangle ptlst))      ; find super triangle
(setq ptlst (append ptlst supertriangle))      ; append coordinates to the end of vertex list
(setq trianglelst (list (list supertriangle nil)))  ; add supertriangle to the triangle list

(setq i 0)
(setq cab 0) ; CAB debug
(while (< i nv)
   (THINKING (strcat "Processing TIN - " (itoa (/ (* i 100) nv)) "%    "))  ; update progress spinner
   (setq pt (nth i ptlst))
   (setq edgelst nil)                                ; initialize edge buffer
   (setq j 0)
   (while (and trianglelst (setq triangle (car (nth j trianglelst))))
     (setq flag T)
     (if (not (cadr (nth j trianglelst)))
       (progn
         (setq circle (getcircircumcircle triangle))       ; calculate circumcircle   
(if (< (+ (caar circle) (cadr circle)) (car pt))  ; test point x and (pt) location
   (setq trianglelst (nth_subst j (list (car (nth j trianglelst)) T) trianglelst))
)
(if (isinside pt circle)
   (setq edgelst     (addtriangleedges triangle edgelst)
trianglelst (nth_del j trianglelst)
flag      nil
   )
)
       ) ; end progn
     )   ; end if
     (if flag (setq j (1+ j)) )
   ) ; end while loop
   (setq edgelst (removedoublyedges edgelst fuzzy nulllist))           ; remove all doubly specified edges
   (setq trianglelst (addnewtriangles pt edgelst trianglelst))         ; form new triangles for current point
   (setq i (1+ i))                                                     ; get next vertex
) ; end while loop
(setq trianglelst (purgetrianglelst trianglelst supertriangle fuzzy)) ; remove triangles with supertriangles edges

(foreach triangle (mapcar 'car trianglelst)                           ; draw triangles
   (drawtriangle triangle)
)

(setvar "OSMODE" OLDSNAP)
(setq OLDSNAP nil)
(command ".UNDO" "END")
(setq stop (getvar "date"))
(princ (strcat "\r TIN Complete - Elapsed time: " (rtos (* 86400.0 (- stop start)) 2 2) " secs."))
(setvar "CMDECHO" OLDCMD) 
(princ)
)



; XSORT - Original Shell Sort function replaced with VLISP sort (much quicker :-)  ;
;                                                                                        ;
(defun XSORT ( PTLST /)
  (vl-sort PTLST (function (lambda (e1 e2) (< (car e1) (car e2)) ) ) )
)

; NTH_DEL ;
; ;
; delete the n item in the list (by position, not by value!!) ;
; ;
; Elimina l'oggetto che si trova nella posizione N della lista LST. L'utilizzo di ;
; funzioni ricorsive,oltre a non assicurare maggiore velocità, può creare problemi;
; di overflow dello stack in caso di liste molto lunghe. ;
(defun NTH_DEL (N LST / l)
(repeat n
  (setq l (cons (car lst) l)
lst (cdr lst)
  )
)
(append (reverse l)(cdr lst))
)

; NTH_SUBST ;
; ;
; Replace the index element in the list with new element. This function is ;
; recursive this is not a great solution with a large amount of data. ;
; ;
(defun NTH_SUBST (index new Alist)
(cond
  ((minusp index) Alist)
  ((zerop index)(cons new (cdr Alist)))
  (T (cons (car Alist)(nth_subst (1- index) new (cdr Alist))))
)
)

; GETPTLIST ;
; ;
; sset -> list (p1 p2 p3 ... pn) ;
; ;
(defun GETPTLIST (ss1 / i pt ptlst)
(if (not (zerop (sslength ss1)))
  (progn
   (setq i 0)
   (while
     (setq pt (ssname ss1 i))
     (setq ptlst (cons (cdr (assoc 10 (entget pt))) ptlst))
     (setq i (1+ i))
   )
  )
)
ptlst
)

; FINDSUPERTRIANGLE ;
; ;
; Search the supertriangle that contain all points in the data set ;
; ;
(defun FINDSUPERTRIANGLE (ptlst / xmax xmin ymax ymin dx dy dmax xmid ymid
          trx1 trx2 trx3 try1 try2 try3 trz1 trz2 trz3
)
(setq xmax (apply 'max (mapcar 'car ptlst))
       xmin (apply 'min (mapcar 'car ptlst))
       ymax (apply 'max (mapcar 'cadr ptlst))
       ymin (apply 'min (mapcar 'cadr ptlst))
       dx (- xmax xmin)
       dy (- ymax ymin)
       dmax (max dx dy)
       xmid (* (+ xmax xmin) 0.5)
       ymid (* (+ ymax ymin) 0.5)
       trx1 (- xmid (* dmax 2.0))
       try1 (- ymid dmax)
       trz1 0.0
       trx2 xmid
       try2 (+ ymid dmax)
       trz2 0.0
       trx3 (+ xmid (* dmax 2.0))
       try3 (- ymid dmax)
       trz3 0.0       
)
(list (list trx1 try1 trz1)
       (list trx2 try2 trz2)
       (list trx3 try3 trz3)
)
)




;;=============================================================
;;  Changes by CAB 03/13/06
;;  replaced the GETCIRCIRCUMCIRCLE routine
;;=============================================================

(defun getcircircumcircle (triangle / p1 p2 p3 pr1 pr2 cen rad bisector)
  ;;  return a pt list for a perpendicular bisector 20 units long
  (defun bisector (p1 p2 / perp_ang midpt)
    (setq p1       (list (car p1) (cadr p1)) ; make sure 2d point
          perp_ang (+ (angle p1 p2) (/ pi 2.0))) ; perpendicular angle
    (setq midpt (mapcar '(lambda (pa pb) (+ (/ (- pb pa) 2.0) pa)) p1 p2))
    (list (polar midpt perp_ang 10) (polar midpt (+ pi perp_ang) 10))
  )
  (setq p1  (car triangle)
        p2  (cadr triangle)
        p3  (caddr triangle)
        pr1 (bisector p1 p2)
        pr2 (bisector p1 p3)
        cen (inters (car pr1) (cadr pr1) (car pr2) (cadr pr2) nil)
        rad (distance cen p1)
  )
  (list cen rad)
)
;;=============================================================



; ISINSIDE ;
; ;
; test if pt is inside a circle ;
; ;
(defun ISINSIDE (pt circle)
(setq ctr (car circle)
       rad (cadr circle)
)
(< (distance pt ctr) rad)
)

; ADDTRIANGLEEDGES ;
; ;
; add triangle edges at the edge queue ;
; ;
(defun ADDTRIANGLEEDGES (triangle edgelst)
(append edgelst (list (list (car triangle)  (cadr triangle))
                       (list (cadr triangle) (caddr triangle))
                       (list (caddr triangle)(car triangle))
                 )
)
)

; DRAWTRIANGLE ;
; ;
; the fun side of the algorithm. Draw triangulation. ;
; ;
(defun DRAWTRIANGLE (triangle)
  (entmake (list (cons 0 "3DFACE") (cons 10 (car triangle))  (cons 11 (caddr triangle))
                                   (cons 12 (cadr triangle)) (cons 13 (cadr triangle))))
)

; EQUALMEMBER ;
; ;
; Check if "item" is in "lista" or not by equality test. With real number the ;
; standard fuction "member" not work correctly. ;
; ;
(defun EQUALMEMBER (item lista fuzzy /)
(apply 'or (mapcar '(lambda (x) (equal x item fuzzy)) lista))
)         

; REMOVEDOUBLYEDGES ;
; ;
; Test the edge queue to remove duplicates (warning CW & CCW!) ;
; ;
(defun REMOVEDOUBLYEDGES (edgelst fuzzy nulllist /)
(setq j 0)
(while (< j (length edgelst))
  (setq k (1+ j))
  (while (< k (length edgelst))
   (if
    (or (and (equal (car  (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
        )
        (and (equal (car  (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
        )
    )
    (setq edgelst (nth_subst j nulllist edgelst)
          edgelst (nth_subst k nulllist edgelst)
    )
   )
   (setq k (1+ k))
  )
  (setq j (1+ j))
)
edgelst
)

; ADDNEWTRIANGLES ;
; ;
; Add new triangle generated by pt to triangle list. ;
; ;
(defun ADDNEWTRIANGLES (pt edgelst trianglelst / j triangle )
(setq j 0)
(while (< j (length edgelst))
  (if (nth j edgelst)
   (setq triangle    (cons pt (nth j edgelst))
         trianglelst (cons (list triangle nil) trianglelst)
   )
  )
  (setq j (1+ j))
)
trianglelst
)

; PURGETRIANGLELST ;
; ;
; replace all triangles that share a vertex with supertriangle ;
; ;
(defun PURGETRIANGLELST (trianglelst supertriangle fuzzy /)
(setq j 0)
(while (and trianglelst (setq triangle (car (nth j trianglelst))))
  (if (apply 'or
             (mapcar '(lambda (x) (equalmember x supertriangle fuzzy))
                     triangle
             )
      )
   (setq trianglelst (nth_del j trianglelst))
   (setq j (1+ j))
  )
)
)


;                                       ;
; THINKING - STANDARD PROGRESS SPINNER  ;
;                                       ;
(defun THINKING (prmpt)
  (setq THINK-CNT (1+ THINK-CNT))
  (princ (strcat "\r" (nth (rem THINK-CNT 4) '("\|" "\/" "\-" "\\")) prmpt))
)


; ********************************* END OF CODING *******************************************
(princ "\n'TRIANGULATE' Loaded \n")
(princ)
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

Didge

  • Bull Frog
  • Posts: 211
Re: Triangulation (re-visited)
« Reply #10 on: March 14, 2006, 06:18:47 AM »
Thats seems to have fixed the bug, thank-you very very muchly CAB.
I can't believe how quickly your mind debugs, whatsmore, like a few others on this site, you achieve it with just a few lines of code. I'm now off to study your changes in great detail.

Thanks for the link LE, unfortunately my employer contracts it's  IT management out to an external company, consequently any sign of an exe, com or DLL incurs disciplinary action :-(
We're therefore restricted to plain acad extensions, although it is amazing how much can be achieved through lisp alone.

I use the excellent 'Ezysurf' terrain modelling package on my home PC, but sadly this requires an exe file for it's licence management.

Ironically, I guess nobody this end realises that lisp could do just as much damage to a network as a virus.
Think Slow......

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Re: Triangulation (re-visited)
« Reply #11 on: June 30, 2006, 06:41:17 AM »
CAB, your code is very good:)
and it let me like the solution.
recently, eachy ( a lisp expert) tell me a website
http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/
I try my best to translate it to Lisp
But my codes is too tedious.
I think after some changes, it could done the effect done by CAB.


Code: [Select]
;;; 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 ([n,m] triangle ntri 1) j 1 triangle))
  (setq triangle (qj-setnmth ([n,m] triangle ntri 2) j 2 triangle))
  (setq triangle (qj-setnmth ([n,m] 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
    (/= ([n,m] edges 1 j) 0)
    (/= ([n,m] edges 2 j) 0)
  )
(progn
  (setq k (1+ j))
  (while (<= k Nedge)
    (if (and
  (/= ([n,m] edges 1 k) 0)
  (/= ([n,m] edges 2 k) 0)
)
      (if (= ([n,m] edges 1 j) ([n,m] edges 2 k))
(if (= ([n,m] edges 2 j) ([n,m] 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
    (/= ([n,m] edges 1 j) 0)
    (/= ([n,m] edges 2 j) 0)
  )
(progn
  (setq ntri (1+ ntri))
  (setq triangle (qj-setnmth ([n,m] edges 1 j) ntri 1 triangle))
  (setq triangle (qj-setnmth ([n,m] 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
  (> ([n,m] triangle i 1) nvert)
  (> ([n,m] triangle i 2) nvert)
  (> ([n,m] triangle i 3) nvert)
)
      (progn
(setq triangle (qj-setnmth ([n,m] triangle ntri 1) i 1 triangle))
(setq triangle (qj-setnmth ([n,m] triangle ntri 2) i 2 triangle))
(setq triangle (qj-setnmth ([n,m] triangle ntri 3) i 3 triangle))
(setq i (1- i)
      temp i
)
(setq ntri (1- ntri))
      )
    )
  )
  (setq Triangulate1 ntri)
  Triangulate1
)
;;; stdlib - to substitute the i element of the list to new value
(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)
)
;;; substitute the i row j column of a 2 dimension array
(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
)
;;; get the n row m column of the two dimension array
(defun [n,m] (a n m / i)      
  (setq i (nth (1- m) (nth (1- n) a)))
  i
)
;;; get the n row of the one dimension array
(defun [n] (a n / i)      
  (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
)

effect

http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #12 on: October 17, 2008, 10:13:33 AM »
Recently, I too needed to use the triangulation program...


Code: [Select]
(defun c:test (/ I L S)
 (princ (strcat "\n select points"))
 (if (setq i 0
           s (ssget '((0 . "POINT")))
     ) ;_  setq
  (progn (repeat (sslength s)
          (setq l (cons (cdr (assoc 10 (entget (ssname s i)))) l)
                i (1+ i)
          ) ;_  setq
         ) ;_  repeat
         (eea-delone-triangulate i l)
  ) ;_  progn
 ) ;_  if
) ;_  defun
(defun eea-delone-triangulate (i1  L   /   A   A1  A2  A3
                               I   I2  L1  L2  L3  LP  MA
                               MI  P   S   TI  TR  X1  X2
                               Y1  Y2
                              )
 ;;*********************************************************
 ;;
 ;; Written by  ElpanovEvgeniy
 ;; 17.10.2008
 ;; Program triangulate an irregular set of 3d points.     
 ;;
 ;;*********************************************************
 (if l
  (progn
   (setq ti (car (_VL-TIMES))
         i  1
         i1 (/ i1 100.)
         i2 0
         l  (vl-sort
             (mapcar
              (function (lambda (p)
                         (list (/ (fix (* (car p) 1000)) 1000.)
                               (/ (fix (* (cadr p) 1000)) 1000.)
                               (caddr p)
                         ) ;_  list
                        ) ;_  lambda
              ) ;_  function
              l
             ) ;_  mapcar
             (function (lambda (a b) (>= (car a) (car b))))
            ) ;_  vl-sort
         x2 (caar l)
         y1 (cadar l)
         y2 y1
   ) ;_  setq
   (while l
    (setq a  (fix (caar l))
          a1 (list (car l))
          l  (cdr l)
    ) ;_  setq
    (while (and l (= (fix (caar l)) a))
     (setq a2 (car l))
     (if (<= (cadr a2) y1)
      (setq y1 (cadr a2))
      (if (> (cadr a2) y2)
       (setq y2 (cadr a2))
      ) ;_  if
     ) ;_  if
     (setq a1 (cons (car l) (vl-remove a2 a1))
           l  (cdr l)
     ) ;_  setq
    ) ;_  while
    (foreach a a1 (setq lp (cons a lp)))
   ) ;_  while
   (setq x1 (caar lp)
         a  (list (/ (+ x1 x2) 2) (/ (+ y1 y2) 2))
         a1 (distance a (list x1 y1))
         ma (+ (car a) a1 a1)
         mi (- (car a) a1)
         s  (list (list ma (cadr a) 0)
                  (list mi (+ (cadr a) a1 a1) 0)
                  (list (- (car a) a1) (- (cadr a) a1 a1) 0)
            ) ;_  list
         l  (list (cons x2 (cons a (cons (+ a1 a1) s))))
         ma (1- ma)
         mi (1+ mi)
   ) ;_  setq
   (while lp
    (setq p  (car lp)
          lp (cdr lp)
          l1 nil
    ) ;_  setq
    (while l
     (setq tr (car l)
           l  (cdr l)
     ) ;_  setq
     (cond
      ((< (car tr) (car p)) (setq l2 (cons (cdddr tr) l2)))
      ((< (distance p (cadr tr)) (caddr tr))
       (setq tr (cdddr tr)
             a1 (car tr)
             a2 (cadr tr)
             a3 (caddr tr)
             l1 (cons (list (+ (car a1) (car a2))
                            (+ (cadr a1) (cadr a2))
                            a1
                            a2
                      ) ;_  list
                      (cons (list (+ (car a2) (car a3))
                                  (+ (cadr a2) (cadr a3))
                                  a2
                                  a3
                            ) ;_  list
                            (cons (list (+ (car a3) (car a1))
                                        (+ (cadr a3) (cadr a1))
                                        a3
                                        a1
                                  ) ;_  list
                                  l1
                            ) ;_  cons
                      ) ;_  cons
                ) ;_  cons
       ) ;_  setq
      )
      (t (setq l3 (cons tr l3)))
     ) ;_  cond
    ) ;_  while
    (setq l  l3
          l3 nil
          l1 (vl-sort l1
                      (function (lambda (a b)
                                 (if (= (car a) (car b))
                                  (<= (cadr a) (cadr b))
                                  (< (car a) (car b))
                                 ) ;_  if
                                ) ;_  lambda
                      ) ;_  function
             ) ;_  vl-sort
    ) ;_  setq
    (while l1
     (if (and (= (caar l1) (caadr l1))
              (= (cadar l1) (cadadr l1))
         ) ;_  and
      (setq l1 (cddr l1))
      (setq l  (cons (eea-data-triangle p (cddar l1)) l)
            l1 (cdr l1)
      ) ;_  setq
     ) ;_  if
    ) ;_  while
    (if (and (< (setq i (1- i)) 1) (< i2 100))
     (progn
      (setvar
       "MODEMACRO"
       (strcat
        "     "
        (itoa (setq i2 (1+ i2)))
        " %    "
        (substr
         "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
         1
         i2
        ) ;_  substr
        (substr
         "..."
         1
         (- 100 i2)
        ) ;_  substr
       ) ;_  strcat
      ) ;_  setvar
      (setq i i1)
     ) ;_  progn
    ) ;_  if
   ) ;_  while
   (foreach a l (setq l2 (cons (cdddr a) l2)))
   (setq
    l2 (vl-remove-if-not
        (function
         (lambda (a)
          (and (< mi (caadr a) ma) (< mi (caaddr a) ma))
         ) ;_  lambda
        ) ;_  function
        l2
       ) ;_  vl-remove-if
   ) ;_  setq
   (foreach a l2
    (entmake (list (cons 0 "3DFACE")
                   (cons 10 (car a))
                   (cons 11 (car a))
                   (cons 12 (cadr a))
                   (cons 13 (caddr a))
             ) ;_  list
    ) ;_  entmake
   ) ;_  foreach
  ) ;_  progn
 ) ;_  if
 (setvar "MODEMACRO" "")
 (princ (strcat "\n "
                (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4)
                " secs."
        ) ;_  strcat
 ) ;_  princ
 (princ)
) ;_  defun
(defun eea-data-triangle (P1 l / A A1 P2 P3 P4 S)
 ;;*********************************************************
 ;;
 ;; Written by  ElpanovEvgeniy
 ;; 17.10.2008
 ;; Calculation of the centre of a circle and circle radius
 ;; for program triangulate
 ;;
 ;; (eea-data-triangle (getpoint)(list(getpoint)(getpoint)))
 ;;*********************************************************
 (setq p2 (car l)
       p3 (cadr l)
       p4 (list (car p3) (cadr p3))
 ) ;_  setq
 (if
  (not
   (zerop
    (setq s (sin (setq a (- (angle p2 p4) (angle p2 p1)))))
   ) ;_  zerop
  ) ;_  not
  (progn (setq a  (polar p4
                         (+ -1.570796326794896 (angle p4 p1) a)
                         (setq a1 (/ (distance p1 p4) s 2.))
                  ) ;_  polar
               a1 (abs a1)
         ) ;_  setq
         (list (+ (car a) a1) a a1 p1 p2 p3)
  ) ;_  progn
 ) ;_  if
) ;_  defun

Additional programs on a subject (CHALLENGE: Triangulation)

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #13 on: October 17, 2008, 10:31:30 AM »
The main task of my code - to work quickly.

The report of job of the program on my the worker computer:

AutoCAD 2008
Windows XP sp2
Intel(R) Core(TM)2 Duo CPU
E8400 @ 3000GHz
3.25 Gb  RAM

Code: [Select]
Command:
Command: (LOAD "D:/Work/triangle/test/test.VLX") nil
Command: test
 select points
Select objects: Specify opposite corner: 100000 found
Select objects:
 39.0630 secs.
Command:

« Last Edit: October 17, 2008, 10:59:03 AM by ElpanovEvgeniy »

Arthur Gan

  • Guest
Re: Triangulation (re-visited)
« Reply #14 on: October 17, 2008, 01:16:26 PM »
there's no information about and your email is hidden and i noticed you contributed a lot of codings..
may i know how to address you formally?
Test this:

Code: [Select]
;******************************************************************;
; TRIANGULATE - Lisp command to create a TIN from 3D points.       ;
; ===========                                                      ;
;                                                                  ;
; Written by Daniele Piazza, ADN member Mechanical Solution s.r.l. ;
; http://pdcode.com/code.htm                                       ;
;                                                                  ;
; Original C coding "Triangulate"  written by PAUL BOURKE          ;
; http://astronomy.swin.edu.au/~pbourke/modelling/triangulate/     ;
;                                                                  ;
; This program triangulates an irregular set of points.            ;
; You can replace some code (sorting, list manipulation,...) with  ;
; VLisp functions to reduce the execution time.                    ;
;                                                                  ;
; This code is not seriously tested, if you find a bug...sorry!!   ;
; Goodbye, Daniele                                                 ;
;*******************************************************************
;
;;
;;  Changes by CAB 03/13/06
;;  replaced the GETCIRCIRCUMCIRCLE routine
;;
(defun C:TRIANGULATE (/ fuzzy nulllist ss1 ptlst nv supertriangle trianglelst i j k edgelst
                        circle pt flag perc)

(setq OLDCMD (getvar "CMDECHO"))
(setvar "CMDECHO" 0)
(command ".UNDO" "GROUP")
(setq OLDSNAP (getvar "OSMODE"))
(setvar "OSMODE" 0)
(setq fuzzy 1e-8)   ; tolerance in equality test
(setq nulllist nil)
 
(princ "\nSelect points...")
(setq ss1 (ssget '((0 . "POINT"))))
(setq start (getvar "date") THINK-CNT 0)            ; initiate timer & Progress Spinner Counter
(setq ptlst (getptlist ss1))                        ; convert selection set to point list
(setq ptlst (xsort ptlst))                          ; sort point list by X co-ordinate
(setq nv (length ptlst))      ; number of points
(setq supertriangle (findsupertriangle ptlst))      ; find super triangle
(setq ptlst (append ptlst supertriangle))      ; append coordinates to the end of vertex list
(setq trianglelst (list (list supertriangle nil)))  ; add supertriangle to the triangle list

(setq i 0)
(setq cab 0) ; CAB debug
(while (< i nv)
   (THINKING (strcat "Processing TIN - " (itoa (/ (* i 100) nv)) "%    "))  ; update progress spinner
   (setq pt (nth i ptlst))
   (setq edgelst nil)                                ; initialize edge buffer
   (setq j 0)
   (while (and trianglelst (setq triangle (car (nth j trianglelst))))
     (setq flag T)
     (if (not (cadr (nth j trianglelst)))
       (progn
         (setq circle (getcircircumcircle triangle))       ; calculate circumcircle   
(if (< (+ (caar circle) (cadr circle)) (car pt))  ; test point x and (pt) location
   (setq trianglelst (nth_subst j (list (car (nth j trianglelst)) T) trianglelst))
)
(if (isinside pt circle)
   (setq edgelst     (addtriangleedges triangle edgelst)
trianglelst (nth_del j trianglelst)
flag      nil
   )
)
       ) ; end progn
     )   ; end if
     (if flag (setq j (1+ j)) )
   ) ; end while loop
   (setq edgelst (removedoublyedges edgelst fuzzy nulllist))           ; remove all doubly specified edges
   (setq trianglelst (addnewtriangles pt edgelst trianglelst))         ; form new triangles for current point
   (setq i (1+ i))                                                     ; get next vertex
) ; end while loop
(setq trianglelst (purgetrianglelst trianglelst supertriangle fuzzy)) ; remove triangles with supertriangles edges

(foreach triangle (mapcar 'car trianglelst)                           ; draw triangles
   (drawtriangle triangle)
)

(setvar "OSMODE" OLDSNAP)
(setq OLDSNAP nil)
(command ".UNDO" "END")
(setq stop (getvar "date"))
(princ (strcat "\r TIN Complete - Elapsed time: " (rtos (* 86400.0 (- stop start)) 2 2) " secs."))
(setvar "CMDECHO" OLDCMD) 
(princ)
)



; XSORT - Original Shell Sort function replaced with VLISP sort (much quicker :-)  ;
;                                                                                        ;
(defun XSORT ( PTLST /)
  (vl-sort PTLST (function (lambda (e1 e2) (< (car e1) (car e2)) ) ) )
)

; NTH_DEL ;
; ;
; delete the n item in the list (by position, not by value!!) ;
; ;
; Elimina l'oggetto che si trova nella posizione N della lista LST. L'utilizzo di ;
; funzioni ricorsive,oltre a non assicurare maggiore velocità, può creare problemi;
; di overflow dello stack in caso di liste molto lunghe. ;
(defun NTH_DEL (N LST / l)
(repeat n
  (setq l (cons (car lst) l)
lst (cdr lst)
  )
)
(append (reverse l)(cdr lst))
)

; NTH_SUBST ;
; ;
; Replace the index element in the list with new element. This function is ;
; recursive this is not a great solution with a large amount of data. ;
; ;
(defun NTH_SUBST (index new Alist)
(cond
  ((minusp index) Alist)
  ((zerop index)(cons new (cdr Alist)))
  (T (cons (car Alist)(nth_subst (1- index) new (cdr Alist))))
)
)

; GETPTLIST ;
; ;
; sset -> list (p1 p2 p3 ... pn) ;
; ;
(defun GETPTLIST (ss1 / i pt ptlst)
(if (not (zerop (sslength ss1)))
  (progn
   (setq i 0)
   (while
     (setq pt (ssname ss1 i))
     (setq ptlst (cons (cdr (assoc 10 (entget pt))) ptlst))
     (setq i (1+ i))
   )
  )
)
ptlst
)

; FINDSUPERTRIANGLE ;
; ;
; Search the supertriangle that contain all points in the data set ;
; ;
(defun FINDSUPERTRIANGLE (ptlst / xmax xmin ymax ymin dx dy dmax xmid ymid
          trx1 trx2 trx3 try1 try2 try3 trz1 trz2 trz3
)
(setq xmax (apply 'max (mapcar 'car ptlst))
       xmin (apply 'min (mapcar 'car ptlst))
       ymax (apply 'max (mapcar 'cadr ptlst))
       ymin (apply 'min (mapcar 'cadr ptlst))
       dx (- xmax xmin)
       dy (- ymax ymin)
       dmax (max dx dy)
       xmid (* (+ xmax xmin) 0.5)
       ymid (* (+ ymax ymin) 0.5)
       trx1 (- xmid (* dmax 2.0))
       try1 (- ymid dmax)
       trz1 0.0
       trx2 xmid
       try2 (+ ymid dmax)
       trz2 0.0
       trx3 (+ xmid (* dmax 2.0))
       try3 (- ymid dmax)
       trz3 0.0       
)
(list (list trx1 try1 trz1)
       (list trx2 try2 trz2)
       (list trx3 try3 trz3)
)
)




;;=============================================================
;;  Changes by CAB 03/13/06
;;  replaced the GETCIRCIRCUMCIRCLE routine
;;=============================================================

(defun getcircircumcircle (triangle / p1 p2 p3 pr1 pr2 cen rad bisector)
  ;;  return a pt list for a perpendicular bisector 20 units long
  (defun bisector (p1 p2 / perp_ang midpt)
    (setq p1       (list (car p1) (cadr p1)) ; make sure 2d point
          perp_ang (+ (angle p1 p2) (/ pi 2.0))) ; perpendicular angle
    (setq midpt (mapcar '(lambda (pa pb) (+ (/ (- pb pa) 2.0) pa)) p1 p2))
    (list (polar midpt perp_ang 10) (polar midpt (+ pi perp_ang) 10))
  )
  (setq p1  (car triangle)
        p2  (cadr triangle)
        p3  (caddr triangle)
        pr1 (bisector p1 p2)
        pr2 (bisector p1 p3)
        cen (inters (car pr1) (cadr pr1) (car pr2) (cadr pr2) nil)
        rad (distance cen p1)
  )
  (list cen rad)
)
;;=============================================================



; ISINSIDE ;
; ;
; test if pt is inside a circle ;
; ;
(defun ISINSIDE (pt circle)
(setq ctr (car circle)
       rad (cadr circle)
)
(< (distance pt ctr) rad)
)

; ADDTRIANGLEEDGES ;
; ;
; add triangle edges at the edge queue ;
; ;
(defun ADDTRIANGLEEDGES (triangle edgelst)
(append edgelst (list (list (car triangle)  (cadr triangle))
                       (list (cadr triangle) (caddr triangle))
                       (list (caddr triangle)(car triangle))
                 )
)
)

; DRAWTRIANGLE ;
; ;
; the fun side of the algorithm. Draw triangulation. ;
; ;
(defun DRAWTRIANGLE (triangle)
  (entmake (list (cons 0 "3DFACE") (cons 10 (car triangle))  (cons 11 (caddr triangle))
                                   (cons 12 (cadr triangle)) (cons 13 (cadr triangle))))
)

; EQUALMEMBER ;
; ;
; Check if "item" is in "lista" or not by equality test. With real number the ;
; standard fuction "member" not work correctly. ;
; ;
(defun EQUALMEMBER (item lista fuzzy /)
(apply 'or (mapcar '(lambda (x) (equal x item fuzzy)) lista))
)         

; REMOVEDOUBLYEDGES ;
; ;
; Test the edge queue to remove duplicates (warning CW & CCW!) ;
; ;
(defun REMOVEDOUBLYEDGES (edgelst fuzzy nulllist /)
(setq j 0)
(while (< j (length edgelst))
  (setq k (1+ j))
  (while (< k (length edgelst))
   (if
    (or (and (equal (car  (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
        )
        (and (equal (car  (nth j edgelst)) (cadr (nth k edgelst)) fuzzy)
             (equal (cadr (nth j edgelst)) (car  (nth k edgelst)) fuzzy)
        )
    )
    (setq edgelst (nth_subst j nulllist edgelst)
          edgelst (nth_subst k nulllist edgelst)
    )
   )
   (setq k (1+ k))
  )
  (setq j (1+ j))
)
edgelst
)

; ADDNEWTRIANGLES ;
; ;
; Add new triangle generated by pt to triangle list. ;
; ;
(defun ADDNEWTRIANGLES (pt edgelst trianglelst / j triangle )
(setq j 0)
(while (< j (length edgelst))
  (if (nth j edgelst)
   (setq triangle    (cons pt (nth j edgelst))
         trianglelst (cons (list triangle nil) trianglelst)
   )
  )
  (setq j (1+ j))
)
trianglelst
)

; PURGETRIANGLELST ;
; ;
; replace all triangles that share a vertex with supertriangle ;
; ;
(defun PURGETRIANGLELST (trianglelst supertriangle fuzzy /)
(setq j 0)
(while (and trianglelst (setq triangle (car (nth j trianglelst))))
  (if (apply 'or
             (mapcar '(lambda (x) (equalmember x supertriangle fuzzy))
                     triangle
             )
      )
   (setq trianglelst (nth_del j trianglelst))
   (setq j (1+ j))
  )
)
)


;                                       ;
; THINKING - STANDARD PROGRESS SPINNER  ;
;                                       ;
(defun THINKING (prmpt)
  (setq THINK-CNT (1+ THINK-CNT))
  (princ (strcat "\r" (nth (rem THINK-CNT 4) '("\|" "\/" "\-" "\\")) prmpt))
)


; ********************************* END OF CODING *******************************************
(princ "\n'TRIANGULATE' Loaded \n")
(princ)