TheSwamp

Code Red => AutoLISP (Vanilla / Visual) => Topic started by: Didge on March 13, 2006, 09:43:53 AM

Title: Triangulation (re-visited)
Post by: Didge 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 (http://www.theswamp.org/index.php?topic=9042.msg555712#msg555712)   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)
Title: Re: Triangulation (re-visited)
Post by: CAB 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.
Title: Re: Triangulation (re-visited)
Post by: Didge 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.



 

Title: Re: Triangulation (re-visited)
Post by: CAB 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.
Title: Re: Triangulation (re-visited)
Post by: Didge 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.

Title: Re: Triangulation (re-visited)
Post by: CAB on March 13, 2006, 11:52:31 AM
Hold on, I do see the problem now.
Title: Re: Triangulation (re-visited)
Post by: Didge 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 :-)
Title: Re: Triangulation (re-visited)
Post by: LE 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.   :-(
Title: Re: Triangulation (re-visited)
Post by: CAB 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.
Title: Re: Triangulation (re-visited)
Post by: CAB 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)
Title: Re: Triangulation (re-visited)
Post by: Didge 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.
Title: Re: Triangulation (re-visited)
Post by: qjchen 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

Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy 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) (http://www.theswamp.org/index.php?topic=15784.0)
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy 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:

Title: Re: Triangulation (re-visited)
Post by: Arthur Gan 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)
Title: Re: Triangulation (re-visited)
Post by: CAB on October 17, 2008, 06:15:21 PM
Arthur,
I sent you an email.

Unfortunately I know little about Triangulation but am a very good debugger.  8-)
Title: Re: Triangulation (re-visited)
Post by: XXL66 on January 09, 2010, 08:48:36 AM
@ElpanovEvgeniy : respect...  a triangulation in lisp this fast is amazing ! To bad your site is in russian, seems to be VERY interesting. BTW: you should give it a try in BricsCAD v10, time is reduced 50% compared to ACAD !



Title: Re: Triangulation (re-visited)
Post by: ymg on May 19, 2011, 12:43:13 PM
I've been playing around with the triangulation written by Daniel Piazza and modified by Cab.

Sorry to report I found a few bugs.

Namely in certain case the supertriangle interfere with the triangulation on the convex hull of the points cloud.

Also the PURGETRIANGLELST function sometimes return J instead of trianglelst causing the triangulation to fail.

Find below the revised code:

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
;;
;; Change by ymg 05/19/2011
;;Modified FINDSUPERTRIANGLE and PURGETRIANGLELST
;;Remove recursion in NTH_SUBST
(defun C:TRIANGULATE (/ fuzzy   nulllist  ss1      
      nv ptlst     trianglelst    
      i j   k        circle
      pt flag   perc         edgelst
     )

   (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
   )                               
   (setq ptlst (getptlist ss1))
   (setq ptlst (xsort ptlst))
   (setq nv (length ptlst))
   
   (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)) "%    "))
      (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))
   )
      )
      (setq edgelst (removedoublyedges edgelst fuzzy nulllist))
      (setq trianglelst (addnewtriangles pt edgelst trianglelst))
      (setq i (1+ i))
   )
   
   (setq trianglelst (purgetrianglelst trianglelst supertriangle fuzzy)) ; remove triangles with supertriangles edges

   (foreach triangle (mapcar 'car trianglelst)
     (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 (a b) (< (car a) (car b)))))
)

; 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)))
;      )
;   )
;)
;;;Removed the recursion  ymg --------------------------
(defun NTH_SUBST (index new Alist / temp)
   (cond
      ((minusp index) Alist)
      ((zerop index)(cons new (cdr Alist)))
      ((> index 0)  (while (> index 0)
                       (setq temp  (cons (car alist) temp)
                     alist (cdr alist)
                     index (1- index)
                       )
                    )
                    (append (reverse temp) (list 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 20.0)) ;modified ymg
try1 (- ymid dmax)
trz1 0.0
trx2 xmid
try2 (+ ymid (* dmax 20.0)) ;modified ymg
trz2 0.0
trx3 (+ xmid (* dmax 20.0)) ;modified ymg
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))
      )
   )
   (setq trianglelst trianglelst); modified to return trianglelst was returning J in certain case ymg
   
)






;                                       ;
; 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)

The triangulation still fail on large data set.  I have removed the recursion in the NTH_SUBST as an attempt
to cure this problem.  I'll look at it some more.

Also played with Elpanov triangulation.  It also suffers from problem on the convex hull.
On a small set of 25 points, it missed to triangulate one of the point and 4 triangles on the hull
are missing.

Here is the point set:
     
Code: [Select]
((1652.17 356.759 446.623) (1666.15 431.163 -353.053) (1688.64 379.861 -372.616) (1708.17 888.849 489.959)
 (1763.96 799.643 117.206) (1811.9 678.149 -387.295) (1818.56 140.657 -256.432) (1883.13 226.078 -79.1498)
 (1888.23 124.665 122.761) (1900.26 864.281 -41.0016) (1950.15 730.671 -164.785) (1979.73 671.496 -19.8523)
 (2031.64 260.656 -497.925) (2069.42 69.732 -278.069) (2071.19 123.139 -183.401) (2096.73 383.737 -280.053)
 (2173.55 135.927 283.044) (2241.51 1048.67 47.7767) (2298.4 460.399 211.447) (2304.6 871.301 -156.27)
 (2441.41 517.957 -411.649) (2455.6 695.636 -390.896) (2462.35 249.15 99.3225) (2585.43 387.857 201.498)
 (2591.77 477.032 100.238) (-17456.9 -419.739 0.0) (2121.97 20138.0 0.0) (21700.8 -419.739 0.0))

So far I have made no attempt to find what is wrong in this code.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 19, 2011, 03:32:16 PM
Hello:
With the algorithm of P. Bourke, employs a few seconds.
That may solve your problem?
Best Regards

PD:
;;;Credit to Paul Bourke (pbourke@swin.edu.au) for the original Fortran 77 Program :))
;;;Converted to AutoLISP by Pedro Ferreira (pwp.netcabo.pt/pedro_ferreira)
;;;Check out: http://local.wasp.uwa.edu.au/~pbourke/papers/triangulate/
;;;You can use this code however you like providing the above credits remain in tact
Title: Re: Triangulation (re-visited)
Post by: CAB on May 19, 2011, 03:43:17 PM
Thanks

http://paulbourke.net/papers/triangulate/Triangulator.LSP

Oh, I see there are two.
Title: Re: Triangulation (re-visited)
Post by: ymg on May 19, 2011, 06:35:10 PM
I have tried also those two implementations of the Bourke algorithm and can report that they are OK.

Triangulator the one by Ferreira suffers from the fact that It uses global variables all over the place and could
therefore leave you stranded.

The one by Mihai Popescu is a nicer one.  I have not look in it to analyze the data structure.

Triangulate by Daniele Piazza is also nicely written and works provided that we modify it as per above.

Now the Fast one by Elpanov would certainly be interesting as long as we can fix the problem along the
convex hull.   Also I have not looked at the data structure.

All this to say that we have all kind of triangulation.

However Triangulation are created to ease further treatment of the data like contouring and
creating profiles on the TIN.  In other word we need to create more than just the 3dface.

Look at the following paper http://research.engineering.wustl.edu/~pless/506/l5.html (http://research.engineering.wustl.edu/~pless/506/l5.html)
for the kind of data structure we should be creating, namely a DCEL (Short for Doubly Connected Edge List)

With such a structure it becomes easy to walk on the TIN, go to the neighbor of a 3dface , contour etc.

So far none of the triangulation that we  have do it well.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 20, 2011, 07:36:59 AM
Surprisingly, in my program there is a small error - a typo, but for many years, nobody has corrected the ...
It really is such a complicated program?
 :wink:
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 20, 2011, 08:07:43 AM
Now, especially checked - Playing a random mistake. I accidentally deleted a few lines. The code is similar to the line and the code is understandable, what is missing! If, indeed, no one can understand this program, I promise next week to lay out the full code! But I will be sorry that such a code just copy and do not try to understand...
 :-o
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 20, 2011, 08:15:37 AM
Corrected code faster!  :-D
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 09:20:59 AM
Hi Elpanov,

The code is understandable.  But the fact is that you write very tight code and consequentlly
It is a little bit hard to follow.

However yours is fast, very fast for something written in Autolisp.  I would not have though that
such speed was attainable.  Therefore Hat's off to you !

What I am really interested in is that the triangulation should create a doubly connected edge list
in order to be ready to contour and or profile afterward. Could even be used to create toolpath
on a surface if we contour along the x or y axis then generate gcode.

Also why autolisp? when there is all kind of C++  and arx solutions.  My reason is that autolisp
can be used on any version of autocad without having to worry about recompilation and all.
I have code wich dates back to be before release 12 that still get some use.

Again congratulations on attaining such speed.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 20, 2011, 09:49:48 AM

In other word we need to create more than just the 3dface.

With such a structure it becomes easy to walk on the TIN, go to the neighbor of a 3dface , contour etc.

So far none of the triangulation that we  have do it well.

ymg
Hello :
if you have 3DFACEs you can converted to sliced-3dsolid, to unite all, and you can easily apply Boolean operations or sections. Usually I do and areas of some hectares with resolution of 2 or 5 meters work very well.
Greetings.
PD.: I attached a small example .... final solid (zip) and previous steps (dwg)
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 20, 2011, 09:58:27 AM
oops... :-) a little picture...
(http://img805.imageshack.us/img805/5760/mdtsolido.gif) (http://imageshack.us/photo/my-images/805/mdtsolido.gif/)
Greetings  :-)
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 10:01:30 AM
Here is something I originally posted in general programming on the subject of contouring.
With some corrections...!

If you look at the following link http://paulbourke.net/papers/conrec/ there are some explanation
on what is involved for determining a contour through a single triangle.

Here I have crudely translated the C routine proposed in that page to Autolisp (Tested and No attempt at Optimization)

Code: [Select]
;Create a contour slice through a 3 vertex facet pa, pb, pc
;The contour "level" is a horizontal plane perpendicular to the z axis,
;ie: The equation of the contour plane Ax + By + Cz + D = 0
;has A = 0, B = 0, C = 1, D = -level
;   Return:
;           ((nil nil nil) (nil nil nil)  0)   if the contour plane doesn't cut the facet
;           (  (x y z)       (x y z)      2)   if it does cut the facet
;           ((nil nil nil) (nil nil nil) -1)   for an unexpected occurrence
;If a vertex touches the contour plane nothing need to be drawn

(defun contourfacet (triangle level / pa pb pc sidea sideb sidec p1x p1y p1z p2x p2y p2z)
   (setq pa (car triangle)
         pb (cadr triangle)
         pc (caddr triangle)
   )
   (setq sidea (-(caddr pa)level)
sideb (-(caddr pb)level)
sidec (-(caddr pc)level)
   )
   (cond
      ;---Are all the vertices on one side----------
      ((and (>= sidea 0)  (>= sideb 0) (>= sidec 0))
       (setq ret 0)
      )
      
      ((and (<= sidea 0)  (<= sideb 0) (<= sidec 0))
       (setq ret 0)
      )
      ;---Is pa the only point on a side by itself---
      ((and (/=(sign sidea) (sign sideb))  (/=(sign sidea) (sign sidec)))

          (setq p1x (- (car pa) (* sidea (/ (- (car pc) (car pa)) (- sidec sidea)))))
          (setq p1y (- (cadr pa) (* sidea (/ (- (cadr pc) (cadr pa)) (- sidec sidea)))))
          (setq p1z (- (caddr pa) (* sidea (/ (- (caddr pc) (caddr pa)) (- sidec sidea)))))
          (setq p2x (- (car pa) (* sidea (/ (- (car pb) (car pa)) (- sideb sidea)))))
          (setq p2y (- (cadr pa) (* sidea (/ (- (cadr pb) (cadr pa)) (- sideb sidea)))))
          (setq p2z (- (caddr pa) (* sidea (/ (- (caddr pb) (caddr pa)) (- sideb sidea)))))
          (setq ret 2)
      )
      ;---Is pb the only point on a side by itself---
      ((and (/=(sign sideb) (sign sidea))  (/=(sign sideb) (sign sidec)))

          (setq p1x (- (car pb) (* sideb (/ (- (car pc) (car pb)) (- sidec sideb)))))
          (setq p1y (- (cadr pb) (* sideb (/ (- (cadr pc) (cadr pb)) (- sidec sideb)))))
          (setq p1z (- (caddr pb) (* sideb (/ (- (caddr pc) (caddr pb)) (- sidec sideb)))))
          (setq p2x (- (car pb) (* sideb (/ (- (car pa) (car pb)) (- sidea sideb)))))
          (setq p2y (- (cadr pb) (* sideb (/ (- (cadr pa) (cadr pb)) (- sidea sideb)))))
          (setq p2z (- (caddr pb) (* sideb (/ (- (caddr pa) (caddr pb)) (- sidea sideb)))))
          (setq ret 2)
      )
      ;---Is pc the only point on a side by itself---
      ((and (/=(sign sidec) (sign sidea))  (/=(sign sidec) (sign sideb)))

          (setq p1x (- (car pc) (* sidec (/ (- (car pa) (car pc)) (- sidea sidec)))))
          (setq p1y (- (cadr pc) (* sidec (/ (- (cadr pa) (cadr pc)) (- sidea sidec)))))
          (setq p1z (- (caddr pc) (* sidec (/ (- (caddr pa) (caddr pc)) (- sidea sidec)))))
          (setq p2x (- (car pc) (* sidec (/ (- (car pb) (car pc)) (- sideb sidec)))))
          (setq p2y (- (cadr pc) (* sidec (/ (- (cadr pb) (cadr pc)) (- sideb sidec)))))
          (setq p2z (- (caddr pc) (* sidec (/ (- (caddr pb) (caddr pc)) (- sideb sidec)))))
          (setq ret 2)
      )
      (t (setq ret -1))
   )
   (list (list p1x p1y p1z) (list p2x p2y p2z) ret)
)


;--Define the Signum function-----
(defun sign (x)
   ( cond
      (( minusp x ) -1 )
      (( zerop x)    0 )
      ( t            1 )
   )
)

Now in this next link http://www.originlab.com/www/helponline/Origin/en/UserGuide/Creating_Contour_Graphs.html#Drawing_of_contour_lines
there is a nice explanation on how you should follow a contour on the TIN.

Basically you find a triangle which has a crossing for the level under consideration, calculate the two points p1 and p2 that intersect the triangle
with above routine. The first point is saved as the starting point of a polyline.  This triangle is marked as traversed.

You then move to the neighboring triangle (The one which  has a reversed edge called a twin), calculate the two points again, mark as traversed until you reach the outside of the TIN (You are sitting on  an edge that has no twin) or you reach your starting point.

We  now check if any other triangle still has a crossing for that level, If it is the case we start a new polyline on the same level.


We then repeat the previous steps for the next level.

Now who is up to come up with an efficient way to do this ?
I am quite rusty with my Autolisp and I am just starting with the visual stuff.

Now you know why I am after a DCEL connected structure from the triangulation.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 10:06:34 AM
Sofito_Soft,

That's an interesting way to do it.

How fast would it be if we wanted to generate a contour map from it ?

I will look at it.

Thanks,

ymg
Title: Re: Triangulation (re-visited)
Post by: VovKa on May 20, 2011, 12:14:54 PM
Now, especially checked - Playing a random mistake. I accidentally deleted a few lines. The code is similar to the line and the code is understandable, what is missing! If, indeed, no one can understand this program, I promise next week to lay out the full code! But I will be sorry that such a code just copy and do not try to understand...
 :-o
hmmm, your variables naming style scares us off, Evgeniy :)
once i tested your code and i had to change
this
Code: [Select]
(vl-remove-if-not
    (function (lambda (a) (and (< mi (caadr a) ma) (< mi (caaddr a) ma))))
    l2
)
to this
Code: [Select]
(vl-remove-if
    (function
        (lambda (a)
            (or (vl-some (function (lambda (c) (vl-position c s))) a) (null a))
        )
    )
    l2
)
to make it work

by i don't really remember why i've done that...
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 20, 2011, 02:22:40 PM
Sofito_Soft,

That's an interesting way to do it.

How fast would it be if we wanted to generate a contour map from it ?

I will look at it.

Thanks,

ymg
Hello:
with A2008 in a pc portable cpu 2 duo T5750 2gh and 3 gb ram
for calculating 3DFACEs:
1000 points = 45 seg....
5000 points = 1 hour....
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 20, 2011, 03:16:25 PM
Hello swampy people:
A comparison between the solution Eugeni (cyan ) and Bourke (Magenta )
The original 1000  points (yellow) are a real example of a terrain slightly irregular.
(http://img840.imageshack.us/img840/83/eugenivsbourke.gif) (http://imageshack.us/photo/my-images/840/eugenivsbourke.gif/)
Time is infinitely better for Eugeni (1 sec) Bourke (35 sec)
The surface is substantially the same.
The perimeter is also better for Eugeni ...
On the left overlapping the 2 solutions ( 3dfaces ).
I will stop using the Bourke and I will from now fans of Eugene. :-) :-) :-)
Thank you very much for sharing, Eugeni.
Greetings.
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 05:40:20 PM
Sofito_soft,

That is a lot of time in order to unite al those 3dface.

Whereby if we would preserve the data structure of the triangulation in a double edged list it would take very
little time to contour.  Most of the work is already done by the triangulation.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 20, 2011, 10:38:20 PM
Sofito_soft,

That is a lot of time in order to unite al those 3dface.

Whereby if we would preserve the data structure of the triangulation in a double edged list it would take very
little time to contour.  Most of the work is already done by the triangulation.

ymg
Hello :
That's true. If you just want to contour, no need to tape "3DFACE" as each 3DFACE have what it takes to calculate its intersection with a plane parallel to XY with Z determined. I think what I have done. I'll look and tell you something.
Greetings.
PS: I will paste a nice slope map calculated only with the information of the triangulation.
(http://img841.imageshack.us/img841/9272/mapadependientes.gif) (http://imageshack.us/photo/my-images/841/mapadependientes.gif/)
Best regards.

 
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 10:49:01 PM
I did a preliminary run through the code of Evgenyi and It seems to me that there is a logic error.

He sets up the first two point of the list and proceed to insert point from this point on.

This would explain why we are missing edges on the convex hull.  Unless he uses an algorithm
that is unknown to me.  The initial triangle has to be carefully selected in order to insure a
complete Delaunay triangulation.

Most everybody uses the midpoint of the bounding box and build a triangle (supertriangle) from
this point.  Piazza's had problem because that triangle was not big enough even though
it did contain all the points.

It you look at the example by Bourke in C, he actually offset that midpoint by (* 20 dmax),
dmax being the biggest of (- ymax ymin) or (- xmax xmin).  That first triangle can be as big as we like
but can be too small and interfere with some point on the convex hull, even if the condition that it should
enclosed all the points is met.  Some authors advocate going to infinity.  Obviously we cannot meet
this condition, therefore the (* dmax 20) is reasonable and should be applied to all three points.

Anyway Evgenyi promised us a correction to his code we will see.  If I am right about this conjecture
addind the supertriangle would not slow his code in any significant way since it  only add one more triangle.

Preserving the data structure Vertex, Edges and Faces might be another story.  As far as I can tell we are
left with only L2 which contains the list of faces.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on May 20, 2011, 11:05:07 PM
Sofito_Soft,

Quote
Hello :
That's true. If you just want to contour, no need to tape "3DFACE" as each 3DFACE have what it takes to calculate its intersection with a plane parallel to XY with Z determined. I think what I have done. I'll look and tell you something.

Indeed the 3dface has all the information necessary to find if any contour goes through it.
Look at the little routine that I have attached in previous post called Contourfacet.

What is more difficult is to trace that contour from face to face.  We do not want to go by brute force
because It's gonna take way too long.

Look  at this link that I posted before http://www.originlab.com/www/helponline/Origin/en/UserGuide/Creating_Contour_Graphs.html#Drawing_of_contour_lines (http://www.originlab.com/www/helponline/Origin/en/UserGuide/Creating_Contour_Graphs.html#Drawing_of_contour_lines)

Our triangulation could preserve all that is necessary to go from face to face rendering the process efficient.

We don't need to re-invent the wheel, very smart people have come up with way to accomplish it efficiently.
All we need is to do it in Autolisp.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 21, 2011, 01:34:53 AM
Sofito_Soft,

Quote
Hello :
That's true. If you just want to contour, no need to tape "3DFACE" as each 3DFACE have what it takes to calculate its intersection with a plane parallel to XY with Z determined. I think what I have done. I'll look and tell you something.

What is more difficult is to trace that contour from face to face.  We do not want to go by brute force
because It's gonna take way too long.

With my prehistoric face to face LSP:

Command: 3DF->CURVA_NIVEL
Selecciona las 3DFACE
Select objects: Specify opposite corner: 3073 found
Select objects:
Cronometro a 00:00
Tiempo para la funcion "c:3df->curva_nivel" : 20.985 seg.

With a intervale of 0.25 m in Z.-> 33000 contour level lines (separates )
Greetings.
Title: Re: Triangulation (re-visited)
Post by: dgorsman on May 21, 2011, 03:21:18 AM
Tracing contours between triangles may be an inefficient process, since triangles that span multiple elevation levels will have to be processed multiple times.  It may be more efficient to process all the contour segments for each triangle, then process those segments into buckets (for lack of a better term) for later processing.  Optimization would likely involve a pre-sort of the triangles based on the lowest (bottom to top) or highest (top to bottom).  Or, for even more optimization creation of the contour segments as each triangle is created.
Title: Re: Triangulation (re-visited)
Post by: ymg on May 21, 2011, 09:14:25 AM
Quote
Tracing contours between triangles may be an inefficient process, since triangles that span multiple elevation levels will have to be processed multiple times.  It may be more efficient to process all the contour segments for each triangle, then process those segments into buckets (for lack of a better term) for later processing.  Optimization would likely involve a pre-sort of the triangles based on the lowest (bottom to top) or highest (top to bottom).  Or, for even more optimization creation of the contour segments as each triangle is created.

Conrec by Bourke again http://paulbourke.net/papers/conrec/ (http://paulbourke.net/papers/conrec/) which dates back to 1987 advocates such a way. Once your done you still
need to unite all those line segment into polylines.  It is certainly feasible, but again in view of futher processing I strongly believe that we should
setup the triangulation in a way that permits us to traverse the TIN efficiently, be it for contour profile or anything else.

Whether we use a DCEL, a Winged Edge or a Quad edge (See Wikipedia for all these), the reason we triangulate is to organize the Data that is the point list in a way that we'll make it easy to find neighbor and to find in which face any point belongs to.

The way we are doing it we are keeping coordinates for points, coordinates again for 3dface and we are destructing the most important which is the
edges list.  The smart way is you keep your point with a pointer to the edge list.  The edge list is only pointer to the point list. The triangle is only pointer to the edge list.

ymg

Title: Re: Triangulation (re-visited)
Post by: ymg on May 21, 2011, 09:25:55 AM
Quote
What is more difficult is to trace that contour from face to face.  We do not want to go by brute force
because It's gonna take way too long.

With my prehistoric face to face LSP:

Your prehistoric lisp is actually doing the right thing, and probably follows what dgorsman is talking about and what Conrec
advocates.  But you still need to join all these line segments into polylines  and label them.

Now, If we want to do something else after contouring we are back to square one.  We need to find intersection to all the polylines which are already an interpolation.  However If we have  a good TIN you profile from the triangle.

Post your routine if you please. 

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 21, 2011, 04:36:39 PM
Hello
If you put the polilines in Z, then you can filter to select all under a certain height. It is easy then weld all.

In my LSP

( defun c:3df->SEG_curva_nivel ( / )   : 3dfaces a segment of polilines in deteminate Z

( defun c:seg_lwpoly->LWPL_curva_nivel ( / ) : Weld segments in the same Z

Greetings. :-)
Title: Re: Triangulation (re-visited)
Post by: ymg on May 21, 2011, 09:11:31 PM
Quote
If you put the polilines in Z, then you can filter to select all under a certain height. It is easy then weld all.

In my LSP

( defun c:3df->SEG_curva_nivel ( / )   : 3dfaces a segment of polilines in deteminate Z

Muy muchas gracias, Sofito.

I will look at it tomorrow, I am traveling at the moment.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 22, 2011, 03:57:23 PM
Hello:
A friend from another forum CAD, sent me this picture. And indeed, proving Evgeniy Elpanov triangulation, has discovered 2 errors.
A vertex is unbound and 2 triangles overlap.
I have reviewed a bit the code. With the changes I've made ​​and works well is the case of conflict. If another user wants to try, please let us know if the correction is adequate.
Thanks ... Greetings from Madrid.
Partial initial code:
Code: [Select]
(defun c:test (/ I L S tmp a b )
 (princ (strcat "\n select points"))
 (if (setq i 0
           s (ssget '((0 . "POINT")))
     ) ;_  setq
  (progn (repeat (sslength s)
          ( setq tmp (cdr (assoc 10 (entget (ssname s i)))))
          ( setq l (cons (list (/ (fix (* (car tmp ) 1000))  1000. )
                               (/ (fix (* (cadr tmp) 1000))  1000. )
                               (caddr tmp)
                         )
                    l)
                i (1+ i)
          ) ;_  setq
         ) ;_  repeat
         ( setq l ( remove_doubles l ) )
         ( setq  l  (vl-sort              
               l                         
               (function (lambda (a b) (>= (car a) (car b))))                            
                    ) ;_  vl-sort
         )
         ( eea-delone-triangulate  l )

  ) ;_  progn
 ) ;_  if
) ;_  defun
(defun eea-delone-triangulate (  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))
         i1 ( length l )
         i  1
         i1 (/ i1 100.)
         i2 0         
         x2 (caar l)
         y1 (cadar l)
         y2 y1
   ) ;_  setq   
   (while l
    (setq a  (fix (caar l))
          a1 (list (car l))         
      ;    l  (cdr l) ; <<<<<<<<<<<<<<<<<<<<<<< !!!!
    ) ;_  setq
.......
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 22, 2011, 04:03:56 PM
Unfortunately, your changes do not affect the algorithm. You try not to understand the algorithm, but only to correct the error.

ps. On Friday, I come from business trips and show you my code and change ...
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 27, 2011, 04:28:44 PM
program with the corrections:

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
 ;; edit 20.05.2011
 ;; 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 a2 (car l))
    (if (<= (cadr a2) y1)
     (setq y1 (cadr a2))
     (if (> (cadr a2) y2)
      (setq y2 (cadr a2))
     )
    )
    (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)
                           (cons (list (+ (car a2) (car a3)) (+ (cadr a2) (cadr a3)) a2 a3)
                                 (cons (list (+ (car a3) (car a1)) (+ (cadr a3) (cadr a1)) a3 a1) l1)
                           ) ;_  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)))
      (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))
       ) ;_  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))))
             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."))
 (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)))))))
  (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


look at the video, all my edits:
Title: Re: Triangulation (re-visited)
Post by: yarik on May 27, 2011, 04:44:11 PM
Thanks Elpanov, works great
Title: Re: Triangulation (re-visited)
Post by: ymg on May 28, 2011, 01:59:52 PM
Quote
program with the corrections:

You are still missing triangles on the convex hull Evgenyi.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 28, 2011, 03:04:57 PM
You are still missing triangles on the convex hull Evgenyi.

ymg

show the test point cloud
Title: Re: Triangulation (re-visited)
Post by: ymg on May 28, 2011, 04:42:13 PM
Code: [Select]
((1652.17 356.759 446.623) (1666.15 431.163 -353.053) (1688.64 379.861 -372.616) (1708.17 888.849 489.959)
 (1763.96 799.643 117.206) (1811.9 678.149 -387.295) (1818.56 140.657 -256.432) (1883.13 226.078 -79.1498)
 (1888.23 124.665 122.761) (1900.26 864.281 -41.0016) (1950.15 730.671 -164.785) (1979.73 671.496 -19.8523)
 (2031.64 260.656 -497.925) (2069.42 69.732 -278.069) (2071.19 123.139 -183.401) (2096.73 383.737 -280.053)
 (2173.55 135.927 283.044) (2241.51 1048.67 47.7767) (2298.4 460.399 211.447) (2304.6 871.301 -156.27)
 (2441.41 517.957 -411.649) (2455.6 695.636 -390.896) (2462.35 249.15 99.3225) (2585.43 387.857 201.498)
 (2591.77 477.032 100.238) (-17456.9 -419.739 0.0) (2121.97 20138.0 0.0) (21700.8 -419.739 0.0))

Same point cloud as on reply #16 in this thread, Evgenyi

You still miss two small triangle on the outside


ymg
Title: Re: Triangulation (re-visited)
Post by: VovKa on May 28, 2011, 07:06:52 PM
program with the corrections:
can you provide a test points list which will show the difference between old and new versions?
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 29, 2011, 01:18:57 AM
can you provide a test points list which will show the difference between old and new versions?

Code: [Select]
(foreach p '((2516.06 458.26 0.0)
             (2520.33 508.832 0.0)
             (2563.56 516.802 0.0)
             (2562.18 470.832 0.0)
             (2594.14 435.458 0.0)
             (2614.14 445.705 0.0)
             (2611.41 488.645 0.0)
            )
 (entmakex (list '(0 . "point") (cons 10 p)))
)
Title: Re: Triangulation (re-visited)
Post by: VovKa on May 29, 2011, 04:19:20 AM
now i see, thanx
Title: Re: Triangulation (re-visited)
Post by: ymg on May 29, 2011, 10:29:26 PM
In order  to generate contours, I have modified Piazza's Triangulate to work with pointers into ptlst instead of coordinates.

That is the trianglelst is now stored as '(((4 0 2) nil) ((4 2 3) nil) ((3 2 1) T) ((2 0 1) T)) where (nth 4 ptlst) are the coordinates of the first vertex of
the first  triangle.

From trianglelst I now generates an edges list edgelst of the form:
     '((4 0) (0 2) (2 4) (4 2) (2 3) (3 4) (3 2) (2 1) (1 3) (2 0) (0 1) (1 2))

Again we have pointers into ptlst.

With this I can follow triangle from one to the other by finding the reverse of an edge.
That is edge (0 2) is neighbor with edge (2 0).  If there is no reverse it means that we are on
the outside (convex hull) of the triangulation.

We also know that a contour crosses a triangle on two edges except in the case where one
or both points defining an edge are exactly on the level of the contour.

Anyway, here is the modified routine :
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                              ;
;;                                                                        ;
;; Changes by ymg :                                                      ;
;;                                                                        ;
;;Modified FINDSUPERTRIANGLE      19/05/2011                              ;
;;and PURGETRIANGLELST            19/05/2011                              ;
;;                                                                        ;
;;Removed recursion in NTH_SUBST  19/05/2011                              ;
;;                                                                        ;
;;Reverted to original GETCIRCIRCUMCIRCLE routine and changed it's       ;
;;name to GETCIRCUMCIRCLE         22/05/2011                              ;
;;                                                                        ;
;;Modified so that trianglelst and edgelst are now list of                ;
;;indices into ptlst.             25/05/2011                              ;
;;                                                                        ;
;;Removed EQUALMEMBER             25/05/2011                              ;
;;                                                                        ;
;;For Contour generation:                                                ;
;;                                                                        ;
;;Added GETEDGELST                28/05/2011                              ;
;;Added GETCROSSEDEGDGE   28/05/2011                              ;
;;Added MAKE_CONTOUR              28/05/2011                              ;
;;       ;

(defun C:TRIANGULATE (/ ss1   nv      
      i j   k
      circle    pt   flag
      ptlst     edgelst   trianglelst
      oldcmd    oldsnap   supertriangle
      triangle  start     think-cnt
      intv      zmax      zmin
     )

   (setq OLDCMD (getvar "CMDECHO"))
   (setvar "CMDECHO" 0)
   (command ".UNDO" "GROUP")
   (setq OLDSNAP (getvar "OSMODE"))
   (setvar "OSMODE" 0)
   
   (princ "\nSelect points...")
   (setq ss1 (ssget '((0 . "POINT"))))
   (setq start (car (_VL-TIMES))
THINK-CNT 0
   )                               
   (setq ptlst (getptlist ss1))
   (setq ptlst (xsort ptlst))
   (setq nv (length ptlst))
   
   (setq supertriangle (findsupertriangle ptlst))     
   (setq ptlst (append ptlst supertriangle))         
   (setq supertriangle (list nv (1+ nv) (+ 2 nv)))   
 
   (setq trianglelst (list(list supertriangle nil))) 
 
   (setq i 0) 
   
   (while (< i nv)
      (THINKING (strcat "Processing TIN - " (itoa (/ (* i 100) nv)) "%    ") think-cnt)
      (setq pt (nth i ptlst))
      (setq edgelst nil)
      (setq j 0)
      (while  (and trianglelst (setq triangle (car(nth j trianglelst))))
   (setq flag T)
   (if (not (cadr (nth j trianglelst)))
      (progn
(setq circle (GETCIRCUMCIRCLE triangle ptlst))         
(cond
   ((< (+ (caar circle) (cadr circle)) (car pt)) 
       (setq trianglelst (nth_subst j (list (car (nth j trianglelst)) T) trianglelst))
       
   )
   ((isinside pt circle)
       (setq edgelst (addtriangleedges triangle edgelst)
     trianglelst (nth_del j trianglelst)
     flag nil
       )
   )

      )
   )
   (if flag (setq j (1+ j)))
      )
     
      (setq edgelst (removedoublyedges edgelst))
      (setq trianglelst (addnewtriangles i edgelst trianglelst))
      (setq i (1+ i))
   )
   
   (setq trianglelst (purgetrianglelst trianglelst supertriangle))
   (repeat 3
      (setq ptlst (vl-remove (last ptlst) ptlst))
   )
   (foreach triangle (mapcar 'car trianglelst)
      (drawtriangle triangle ptlst)
   )
   
   ;; Now we set-up to trace the contour                                      ;
   
   (setq edgelst (GETEDGELST trianglelst)
            zmax (apply 'max (mapcar 'caddr ptlst))
    zmin (apply 'min (mapcar 'caddr ptlst))
    intv 1
    zmin (+ (fix zmin) intv)
    zmax (fix zmax)
   )
   
   (MAKE_CONTOUR zmin zmax intv edgelst ptlst)

   (setvar "OSMODE" OLDSNAP)
   (setq OLDSNAP nil)
   (command ".UNDO" "END")
   
   (princ (strcat "\r TIN Complete - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) start) 1000.) 2 4) " secs."))
   (princ)
   
   (setvar "CMDECHO" OLDCMD)
   (princ)
)






;; XSORT       ;
;;       ;
;;  Original Shell Sort function replaced with VLISP sort (much quicker :-)   ;
;;                                                                            ;

(defun XSORT ( PTLST /)
  (vl-sort PTLST (function (lambda (a b) (< (car a) (car b)))))
)

; NTH_DEL       ;
;       ;
; delete the n item in the list (by position, not by value!!)       ;
;       ;

(defun NTH_DEL (N LST / l)
   (repeat n
      (setq l (cons (car lst) l)
    lst (cdr lst)
      )
   )
   (append (reverse l) (cdr lst))
)

;(defun NTH_DEL (index lista)
; (setq i -1)
; (apply 'append (mapcar '(lambda (x)
;                                (if (= (setq i (+ 1 i)) index)
;                                    nil
;                                    (list x)
;                                )
;                         )
;                lista
;               )
; )
;)


; NTH_SUBST       ;
;       ;
; Replace the index element in the list with new element.       ;
;                       ;
;       ;

(defun NTH_SUBST (n new lst / l)
   (cond
      ((minusp n) lst)
      ((zerop n)(cons new (cdr lst)))
      (t  (repeat n
             (setq l   (cons (car lst) l)
           lst (cdr lst)
     )
          )
          (append (reverse l) (list new) (cdr lst))
      )
   )
)
;Need to check that one for speed
;(defun NTH_SUBST (n new lst / i)
; (setq i -1)
; (mapcar '(lambda (x)
;           (if (= (setq i (+ 1 i)) n)
;            new
;            x
;           )
;          )
;  lst
; )
;)
     

; GETPTLIST       ;
;       ;
; sset -> list (p1 p2 p3 ... pn)       ;
;       ;

(defun GETPTLIST (ss1 / i pt ptlst)
   (setq i 0)
   (if (not (zerop (sslength ss1)))
      (progn
(repeat (sslength ss1)
     (setq ptlst (cons (cdr (assoc 10 (entget (ssname ss1 i)))) ptlst)
               i (1+ i)
     )
)
      )
   )
   ptlst
)


; FINDSUPERTRIANGLE       ;
;       ;
; Search the supertriangle that contain all points in the data set      ;
;       ;
(defun FINDSUPERTRIANGLE (ptlst / xmax   xmin   ymax   ymin
  dmax   xmid   ymid   x1     x2     x3
  y1     y2     y3   
)
   (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))
dmax (max (- xmax xmin) (- ymax ymin))
xmid (* (+ xmax xmin) 0.5)
ymid (* (+ ymax ymin) 0.5)
   x1 (- xmid (* dmax 20.0)) ;modified ymg
   y1 (- ymid dmax)

   x2 xmid
   y2 (+ ymid (* dmax 20.0)) ;modified ymg

   x3 (+ xmid (* dmax 20.0)) ;modified ymg
   y3 (- ymid dmax)

   )
   (list (list x1 y1 0.0) (list x2 y2 0.0) (list x3 y3 0.0))
)



;;          GETCIRCUMCIRCLE                                                   ;
;;       ;
;;    Find circle passing through three points                                ;
;;       ;

(defun GETCIRCUMCIRCLE (triangle ptlst / p1 p2 p3 p1x p2x p3x
                 p1y p2y p3y d xc yc rad)
 (setq p1 (nth (car triangle) ptlst)
       p2 (nth (cadr triangle) ptlst)
       p3 (nth (caddr triangle) ptlst)
       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)
  (< (distance pt (car circle)) (cadr circle))
)

; 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 ptlst /)
  (entmake (list (cons 0 "3DFACE")
(cons 8 "TIN")
(cons 10 (nth (car triangle)   ptlst))
(cons 11 (nth (caddr triangle) ptlst))
                 (cons 12 (nth (cadr triangle)  ptlst))
(cons 13 (nth (cadr triangle)  ptlst))
   )
  )
)

; REMOVEDOUBLYEDGES       ;
;       ;
; Test the edge queue to remove duplicates (warning CW & CCW!)       ;
;       ;

(defun REMOVEDOUBLYEDGES (edgelst / j k)

   (setq j 0)
   (while (< j (length edgelst))
      (setq k (1+ j))
      (while (< k (length edgelst))
(if
    (or (and (= (car (nth j edgelst)) (car (nth k edgelst)))
     (= (cadr (nth j edgelst)) (cadr (nth k edgelst)))
)
(and (= (car (nth j edgelst)) (cadr (nth k edgelst)))
     (= (cadr (nth j edgelst)) (car (nth k edgelst)))
)
    )
      (setq edgelst (nth_subst j nil edgelst)
    edgelst (nth_subst k nil 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)
   
   (repeat (length edgelst)
      (if (nth j edgelst)
  (setq triangle    (append (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 / a b j triangle)
   (setq j 0)
   (repeat (length trianglelst)
      (setq triangle (car (nth j trianglelst)))
      (if
         (apply 'or (mapcar '(lambda (a) (apply 'or (mapcar '(lambda (b) (= b a )) supertriangle)))
triangle
)
)


   (setq trianglelst (nth_del j trianglelst))
   (setq j (1+ j))
      )
   )
   trianglelst
   
)

;;           GETEDGELST       ;
;;       ;
;;   Create list of EDGES of all triangles in trianglelst.       ;
;;       ;

(defun GETEDGELST (trianglelst / edgelst neighlst)
   (setq edgelst nil)
   (foreach tr trianglelst
      (setq edgelst (cons (list (caar tr)(cadar tr)) edgelst)
    edgelst (cons (list (cadar tr)(caddar tr)) edgelst)
    edgelst (cons (list (caddar tr)(caar tr)) edgelst)
      )
   )   
)


;;      GETCROSSEDEGDGE       ;
;;       ;
;;Traverses the edges list and creates a list of edges that are crossed by a  ;
;;contour level.  It then follows contlst from neigbour to neighbor until it  ;
;;reaches the exterior of the TIN or the starting point.       ;
;;       ;
;;Returns crossedlst which contains list of edges.       ;
;;       ;

(defun GETCROSSEDEDGE (lev edgelst ptlst / e pos nxt contlst lwplst
        z1 z2 crossedlst)
   
   (setq contlst nil)
   
   (foreach e edgelst
      (setq z1 (caddr (nth (car e) ptlst))
    z2 (caddr (nth (cadr e) ptlst))
      )
      (if (= z1 lev) (setq z1 (+ z1 1e-8)))
      (if (= z2 lev) (setq z2 (+ z2 1e-8)))
     
      (if (or (< z1 lev z2)(> z1 lev z2))   
          (setq contlst (cons e contlst))
      )
     
   )
   (setq crossedlst nil)     
   (while contlst
      (setq pos 0)
      ;;Find an edge on the convex hull and start from it if none start at last pos in contlst
      (while (and (member (reverse (nth pos contlst)) contlst) (< pos (1- (length contlst))))
(setq pos (1+ pos))
      )
     
      (setq  lwplst (list (nth pos contlst))
    contlst (nth_del pos contlst)
                pos (- pos (rem pos 2))
             lwplst (cons (nth pos contlst) lwplst)
        nxt (reverse (car lwplst))
    contlst (nth_del pos contlst)
      )
   
      (while (and (setq pos (vl-position nxt contlst)) (not (member nxt lwplst)))
(cond
    ((> (length contlst) 1)
            (setq contlst (nth_del pos contlst)
                  pos (- pos (rem pos 2))
                         lwplst (cons (nth pos contlst) lwplst)
      nxt (reverse (car lwplst))
      contlst (nth_del pos contlst)
            )
    )
    (t (setq lwplst  (cons (nth pos contlst) lwplst)
     contlst (nth_del pos contlst)
       )
    )
)   
      )
      (setq crossedlst (cons lwplst crossedlst))
    )
    crossedlst
 ) 

;;              MAKE_CONTOUR               ;
;;       ;
;; Creates LWPolylines from the edgelst         ;
;;       ;

(defun MAKE_CONTOUR (zmin zmax intv edgelst ptlst /
     lwplst lev edge p1 p2 d pt tmp)

   (setq lev zmin)
   
   (repeat (fix (/ (- zmax zmin) intv))
   
   (setq lwplst (GETCROSSEDEDGE lev edgelst ptlst))
   
   (foreach pline lwplst
      (setq tmp nil)
      (foreach edge pline
(setq   p1 (nth (car  edge) ptlst)
           p2 (nth (cadr edge) ptlst)
                            r (/ (- l (caddr p1)) (- (caddr p2) (caddr p1)));Modified june 2013
                         p1 (list (car p1)(cadr p1))                                 ; Calculation of contour was incorrect
           p2 (list (car p2)(cadr p2))                   
          d (* (distance p1 p2) r)
         pt (polar p1 (angle p1 p2) d)
pt (list (car pt)(cadr pt))
        tmp (cons (cons 10 pt) tmp)
)
      )
      (setq tmp (cons (cons 38 lev) tmp)
    tmp (cons (cons 43 0.0) tmp)
    tmp (cons (cons 70 0) tmp)
    tmp (cons (cons 90 (length pline)) tmp)
    tmp (cons (cons 100 "AcDbPolyline") tmp)
    tmp (cons (cons 100 "AcDbEntity") tmp)
    tmp (cons (cons 8 "Contour") tmp)
    tmp (cons (cons 0 "LWPOLYLINE") tmp)
      )            
             
   (entmake tmp)      
   )
   (setq lev (+ lev intv))
   )   
)   
   


     
;;       THINKING       ;
;;       ;
;;     Standard progress spinner       ;
;;       ;

(defun THINKING (prmpt think-cnt)
   (setq think-cnt (1+ think-cnt))
   (princ (strcat "\r" (nth (rem think-cnt 4) '("\|" "\/" "\-" "\\")) prmpt))
)


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


It now need to be optimized, probably sorting the edges would help.

Here is the point cloud :

Code: [Select]
(  (1652.17 356.759 16.4)
   (1666.15 431.163 16.1)
   (1688.64 379.861 16.2)
   (1708.17 888.849 13.8)
   (1763.96 799.643 14.1)
   (1818.56 140.657 18.3)
   (1883.13 226.078 18.1)
   (1888.23 124.665 18.8)
   (1900.26 864.281 12.8)
   (1950.15 730.671 12.5)
   (1979.73 671.496 12.1)
   (2031.64 260.656 15.8)
   (2069.42 69.732 22.2)
   (2071.19 123.139 19.8)
   (2096.73 383.737 15.4)
   (2173.55 135.927 20.7)
   (2241.51 1048.67 15.1)
   (2298.4 460.399 12.3)
   (2304.6 871.301 12.0)
   (2441.41 517.957 11.5)
   (2455.6 695.636 10.0)
   (2462.35 249.15 18.3)
   (2585.43 387.857 9.8)
   (2591.77 477.032 9.5)
)

I have attached an image of the results, not too sure if it will work.  :?

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 12:48:29 AM
Holas:
Very small differences, but actually looks better. :-)
The best test is made ​​with points that are in one of the sought contour. So, if the new curve (polyline) boundary passes over the points, is that the process can be taken for good.
Greetings  :-)
PS. In the dwg , old solution on color 21, new en cyan / blue
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 02:27:20 AM
Hello all:
The slopes plan also improves, the changes are more evenly distributed.
Thanks Eugeny. :-)
(http://img42.imageshack.us/img42/4808/slopeplane.gif) (http://imageshack.us/photo/my-images/42/slopeplane.gif/)

Greetings.
Title: Re: Triangulation (re-visited)
Post by: ymg on May 30, 2011, 11:19:15 AM
Quote
Holas:
Very small differences, but actually looks better.
The best test is made ​​with points that are in one of the sought contour. So, if the new curve (polyline) boundary passes over the points, is that the process can be taken for good.

Holas, Sofito

I did contour your drawing to get an idea of speed. On a 1m. interval took about 15 sec. on my laptop.

However you end up with joined plines at the proper level.  Also I am sure that this could be optimized
and accelerated quite a bit.  Also the triangle list could be destroyed once the edges list is created.

What I did to make the contour pass over the point is to disturb the elevation of the point(s)
by 1e -8 so that the edge is included in the contlst.

I notice on your revised drawing that some of the color 21 contour turn on themselves which is a
big No...No

Will look at it some more.

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 01:30:17 PM
Hello:

I did contour your drawing to get an idea of speed. On a 1m. interval took about 15 sec. on my laptop.<<<< 6 seg in my PC....is reasonable , for me.

Also I am sure that this could be optimized <<<< There are 3 or 4 methods ( p.e.geometric general calculation ) that are not highly purified. If you have any interest in something really fast ( productive? ), please write me to e-mail private.

I notice on your revised drawing that some of the color 21 contour turn on themselves which is a
big
<<<< is a problem that I could not solve. It occurs when there is a flat surface. The contour, effectively, turns on itself. She surrounds the flat area is closed, leaving some segment no possibility of being soldiers. Any idea to resolve it?

Will look at it some more.<<<< It is a fascinating subject
Greetings.
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 03:28:47 PM
Hello again:
Eugeny: another "extravaganza" of your algorithm ... please watch the yellow circle in extreme of wellow arrow......

(http://img851.imageshack.us/img851/9218/eugenybugtriangulating.gif) (http://imageshack.us/photo/my-images/851/eugenybugtriangulating.gif/)
 
Best regards..
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 30, 2011, 04:45:18 PM
Eugeny: another "extravaganza" of your algorithm ... please watch the yellow circle in extreme of wellow arrow......

I confirm that the program does not give the best result. But I have no need to continue in the near future to work on this program. The program was created for the contest on this site. Problem solver, I have used in the abridged version. All necessary for me, the decisions necessary for me, I already have. I give permission to fully use its program to your goals, including changes to its complement and so on. I ask only to maintain the link in the code for me.

PS. If you need to modify this program, or you want to order from me, a new program for your internal purposes only - please contact me directly offline. I am always ready to consider any offers.

Best Regards, ElpanovEvgeniy
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 05:11:14 PM
Hello again:
Evgeny: thanks for the information and the special program.
Greetings.
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 30, 2011, 05:19:51 PM
I hope my words do not destroy the great atmosphere of communication and mutual learning on the forum!  :?
Title: Re: Triangulation (re-visited)
Post by: pkohut on May 30, 2011, 05:27:44 PM
I hope my words do not destroy the great atmosphere of communication and mutual learning on the forum!  :?


You've contributed a great piece of code ElpanovEvgeniy, Спасибо.
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 30, 2011, 05:31:12 PM
You've contributed a great piece of code ElpanovEvgeniy, Спасибо.

You wrote this code in another language and your code is much faster...
Title: Re: Triangulation (re-visited)
Post by: pkohut on May 30, 2011, 05:57:19 PM
You've contributed a great piece of code ElpanovEvgeniy, Спасибо.

You wrote this code in another language and your code is much faster...


It uses a different algorithm than Paul Bourke's to crank the speed way up, being written in C++ just pushes it over the top more.  :-) With the internal data structures (http://www.theswamp.org/index.php?topic=15784.msg310523#msg310523), support for contouring and break lines would just need a little code.

Like you, it was done for personal interest at the time and achieved its primary goal
Title: Re: Triangulation (re-visited)
Post by: VovKa on May 30, 2011, 06:22:29 PM
another one
http://dwg.ru/dnl/?id=1766
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 08:02:13 PM
Hello again:
Evgeny: the one that says what is right, there is no threat to destroy anything.
In Spanish  is said : it that telling the truth, nor commit sin, nor lying.
I attached a dwg that I think will help you test the algorithm, for the day you get back the program.
Greetings from Madrid.
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 30, 2011, 08:18:23 PM
Adorning with colors the triangulations, the 3DFACE and contours with color by slope (3dfaces ) or level ( polylines) .
Endless possibilities.
Greetings.

(http://img38.imageshack.us/img38/4733/coloreandolastriangulac.gif) (http://imageshack.us/photo/my-images/38/coloreandolastriangulac.gif/)

 
Title: Re: Triangulation (re-visited)
Post by: ymg on May 30, 2011, 08:56:31 PM
Evgenyi,

Too bad you let it go.

The speed of that piece of code is incredible!

However as I told you earlier it is written very tightly which makes it a little difficult to follow.

You are and remain a great contributor.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on May 30, 2011, 09:46:45 PM
Quote
Also I am sure that this could be optimized <<<< There are 3 or 4 methods ( p.e.geometric general calculation ) that are not highly purified. If you have any interest in something really fast ( productive? ), please write me to e-mail private.

I only do it for the fun of it.

My days as a surveyor are done.  With the speed that are obtainable these days and a GPS we could probably keep a dynamic update of a site as it is being graded. For example red zone you need to cut, blue you need to fill, green you are on grade. 

All this while sitting in the bulldozer or the grader.  I know such system do exist but are very pricey.

It is now within our capacity to actually do it.

Fascinating subject as you said.

Hasta Luego,

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 31, 2011, 03:18:30 AM
Hello:
YMG: For example red zone you need to cut, blue you need to fill, green you are on grade.  <<<< fascinating and suggestive idea ...
Especially for large supericies / volumes with large cuts and fills thick ...
Will always be an approximation work ... the landfill is compressed ... the GPS has a tolerance ... the virgin land is washed by the rain ... but it can save many hours of expensive machinery ...
Very good idea! Bravo.
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on May 31, 2011, 03:27:38 AM
I think that the game of chess or the development of algorithms (programs) with a strong opponent (friend), it is a great honor! This session, the fastest way to perfection! Therefore, I always try to participate in the topics "-= {Challenge} =- ***".

My code, simple! Use gradual debugging - drawing lines after each piece of code, the example of five points in one hour, you can see and understand the entire algorithm. To do this, do not need to understand the code, draw the visualization of each segment code in the drawing.
I have heard many times that my code difficult to read. Maybe so, but it just and effective. For example in this problem ...

Each time, laying "a strong -  strong  code", I do contribute to the training of successors. I try to give food to those who seek to overtake me! I'm not interested to give the code for easy copying and use. I earn my bread by writing programs and algorithms. I am manage a process of comprehensive automation of CAD, CAM and CAE.
Title: Re: Triangulation (re-visited)
Post by: ymg on May 31, 2011, 07:50:24 PM
Quote
Very good idea! Bravo.

Well, the idea did not originate with me.

But I know it beats the hell out of driving stakes in the ground to mark cut&fill  :-D

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 31, 2011, 07:52:02 PM
Holas
YMG: Another approach on the navigation above / below a MDT:
Using the contour levels like LWPOLYLINES...  can take advantage using the functions VLAX-CURVE-GETxxxxx. They are native Autocad code and therefore very very fast ...
A progressive sophistication might be, calculate the Z of a point (obtained from a GPS) on the original model and model desirable.
Comparison of Z, the current (real) with the original and desirable, can give much information on the debris (filled) to move and where.
I have done some testing ...with one MDT based on LWPOLYLINE. I post a micro-video... I fail to generate animated GIF, so herewith a AVI  :lol:
Greetings  :-)
Title: Re: Triangulation (re-visited)
Post by: ymg on May 31, 2011, 07:56:49 PM
Quote
I am manage a process of comprehensive automation of CAD, CAM and CAE.

Evgenyi,

I got back into triangulation and contouring as part of a hobby of mine which is writing GCode
for router machine.

There exist very good program like Vectric Aspire but they are so pricey!!

If you want to do 3d routing, thr problem is essantially contouring.

Cheers!

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on May 31, 2011, 09:02:33 PM
Quote
I have done some testing ...with one MDT based on LWPOLYLINE. I post a micro-video... I fail to generate animated GIF, so herewith a AVI 
Greetings 

Unfortunately, Sofito, the video did not show anything.

But don't be discouraged I am also video challenged.  :-D

ymg
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on May 31, 2011, 09:21:02 PM
Hello:
In my pc, I can see. Sorry.
Moving the cursor over the LWPOLYLINES, clicking a point, it is estimated the 2 polylines closest point is interpolated and the original. It displays the Z. .. Now I can do with 2 sets of polylines in the layer "original" and "futuro".  The lsp displays the Z and differences and time calculation
I attached a static image of the video. and a dwg with the new method of double layer/mdt.
Greetings.
(http://img546.imageshack.us/img546/7660/surfincontours.gif) (http://imageshack.us/photo/my-images/546/surfincontours.gif/)
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 01, 2011, 05:32:10 AM
Quote
I am manage a process of comprehensive automation of CAD, CAM and CAE.

Evgenyi,

I got back into triangulation and contouring as part of a hobby of mine which is writing GCode
for router machine.

There exist very good program like Vectric Aspire but they are so pricey!!

If you want to do 3d routing, thr problem is essantially contouring.

Cheers!

ymg

A very strange sentence!
Individually written program, will always be significantly more expensive than to sell in large quantities. Free software has its price.

I sometimes are participating in such projects if the projects are necessary to me as well. Now I write 3D kernel for the direct management of vertices, edges and faces 3dsolid. All this is necessary in a project of precast ferroconcrete. The project is done under a separate technology, which is uncommon and can not be maintained well-known programs. The problem of triangulation, very far from me...
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on June 01, 2011, 06:10:09 AM
  Now I write 3D kernel for the direct management of vertices, edges and faces 3dsolid. All this is necessary in a project of precast ferroconcrete. 
 
You've worked with solids "ACIS"? They are solid autocad. Are well documented.
very very interesting.
Greetings
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 01, 2011, 07:22:19 AM
 
You've worked with solids "ACIS"? They are solid autocad. Are well documented.
very very interesting.
Greetings

Yes, I use a simple command entmakex, entmode, entupd...
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on June 01, 2011, 01:22:18 PM
Eugeny:
I thought rather  directly access to 3DSOLID with the  "SAT" files .
Giles opened the door, with a wonderful ACISDECODE subfunction.
From there, you can get the points, edges, vertices, faces, normal, holes, etc, of the ACIS solids.
I have read that it is developed also for NET.
It is always a good thing that new developments are compatible with something already established.
I do not know if your head warm.
A greeting.  :-)
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 01, 2011, 02:15:12 PM
Eugeny:
I thought rather  directly access to 3DSOLID with the  "SAT" files .
Giles opened the door, with a wonderful ACISDECODE subfunction.
From there, you can get the points, edges, vertices, faces, normal, holes, etc, of the ACIS solids.
I have read that it is developed also for NET.
It is always a good thing that new developments are compatible with something already established.
I do not know if your head warm.
A greeting.  :-)

All this is not the decoding and modeling - Change, for example, add or remove a port complex shape, or simply add / remove a point on one edge. This is the easiest level, and then - the cross section, sections, species in its representation of geometry. The problem of automatically creating species and the cross section with a complete adaptability - changed item, change views, changed views, changed body.

ps. arx not give full access to the geometry. Using arx, you can not add or delete a single point in the body ...
Title: Re: Triangulation (re-visited)
Post by: ymg on June 01, 2011, 04:17:09 PM
Quote
A very strange sentence!
Individually written program, will always be significantly more expensive than to sell in large quantities. Free software has its price.

Of course developing custom solutions is always more expensive.

However If you are doing it as a hobby, buying every nice program out there can break the bank.

Plus programming is another hobby.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 01, 2011, 05:42:42 PM
Of course developing custom solutions is always more expensive.

However If you are doing it as a hobby, buying every nice program out there can break the bank.

Plus programming is another hobby.

ymg

I think you will release the project documentation.
I manage the release of automation software.

If you need to develop small programs or large system, you can contact me with offers. Until now, I work only in the domestic market, but is willing to work with other countries...
Title: Re: Triangulation (re-visited)
Post by: ymg on June 02, 2011, 09:39:27 PM
Evgenyi,

To go back to your triangulation, I did follow part of the code.

When you setup your first triangle in list s it is too small.

The vertex of that triangle are inside the circumcircle of other points on the outside.
In the litterature some authors advocate infinity for the so called supertriangle.

I have attached a picture showing the interference of your s triangle with some
of the points in the cloud.

If you fix that, I believe we have a Delaunay triangulator. 

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on June 02, 2011, 11:58:54 PM
Here I have revised your code:

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
   ;; edit 20.05.2011
   ;; 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 a2 (car l))
    (if (<= (cadr a2) y1)
       (setq y1 (cadr a2))
       (if (> (cadr a2) y2)
  (setq y2 (cadr a2))
       )
    )
    (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)(* 20 a1));changed ymg
       mi (- (car a)(* 20 a1));changed ymg
       s  (list (list ma (cadr a) 0)
(list mi (+ (cadr a)(* 20 a1)) 0);changed ymg
(list mi (- (cadr a)(* 20 a1)) 0);changed ymg
  ) ;_ supertriangle
       l  (list (cons x2 (cons a (cons (+ a1 a1) s))))
       ;mi (- (car a)(* 20 a1));added ymg
       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
  )
  (cons (list (+ (car a2) (car a3))
      (+ (cadr a2) (cadr a3))
      a2
      a3
)
(cons (list (+ (car a3) (car a1))
    (+ (cadr a3) (cadr a1))
    a3
    a1
      )
      l1
)
  ) ;_  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))
   )
  (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))
     ) ;_  strcat
  ) ;_  setvar
  (setq i i1)
       ) ;_  progn
    ) ;_  if
) ;_  while
(foreach a l (setq l2 (cons (cdddr a) l2)))

;purge triangle list
(setq
    l2 (vl-remove-if-not
  (function
     (lambda (a)
(and (< mi (caadr a) ma) (< mi (caaddr a) ma))
     )
  )
  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."
  )
   )
   (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)))))
  )
       )
      (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

I don't think anything I did will affect the speed.

Attached a gif of the result.

ymg

Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on June 03, 2011, 04:39:20 AM
Hello:
ymg :
Another way of navigating in 3d over a MDT ... for a XY point which are the points that cut the 3dsolids generates by 3DFACEs.
In this case, draw lines between points found in the vertical of 3 3dsolid.

Time for the example: 0.374973 seg. for 6 points founds....To improve... :-(

(http://img600.imageshack.us/img600/4368/sandwichentretableroyte.gif) (http://imageshack.us/photo/my-images/600/sandwichentretableroyte.gif/)
Greetings. :-)
Title: Re: Triangulation (re-visited)
Post by: SOFITO_SOFT on June 03, 2011, 07:35:36 AM
Oops...
after some small problems with sub-functions releases  :lol: :lol:....the LSP for navigating in 3d over a MDT ( several 3dsolids )
( c:ll ) ...for initialize and...
( c:pto_en_Z_solid ) for run...
Greetings  :-)


Title: Re: Triangulation (re-visited)
Post by: ymg on June 14, 2011, 07:40:16 PM
Here is Evgeniy's triangulation commented for analysis.

The problem with first triangle is corrected, but note that in certain case it could be too small.

I have renamed most variables and commented.

Code: [Select]
(defun c:test (/ i pl s)
   (princ (strcat "\n select points"))
   (if (setq i 0
     s (ssget '((0 . "POINT")))
       )
      (progn (repeat (sslength s)
(setq pl (cons (cdr (assoc 10 (entget (ssname s i)))) pl)
      i (1+ i)
)
     )
     (triangulate pl)
      )
   )
)


;;********************************************************;
;;                                                        ;
;; Written by  ElpanovEvgeniy                             ;
;; 17.10.2008                                             ;
;; edit 20.05.2011                                        ;
;; Program triangulate an irregular set of 3d points.     ;
;; Modified by ymg June 2011                              ;
;;********************************************************;

(defun triangulate (pl  /  a   b   c   i   i1   i2
    bb  sl al  el  tl  l   ma  mi     
    ti  tr x1  x2  y1  y2  p    r  cp
   )
   
   (if pl
      (progn
(setq ti (car (_VL-TIMES))
       i  1
       i1 (/ (length pl) 100.)
       i2 0
       pl (vl-sort pl
     (function (lambda (a b) (< (car a) (car b))))
  )
       bb (list (apply 'mapcar (cons 'min pl))
                (apply 'mapcar (cons 'max pl))
                  )
               ;Replaced code to get minY and maxY with 3d Bounding Box Routine;
       ;A bit slower but clearer. minZ and maxZ kept for contouring    ;
       
       x1 (caar bb)   ;minX
       x2 (caadr bb)  ;maxX
       y1 (cadar bb)  ;minY
       y2 (cadadr bb) ;maxY
       
       
)
(setq cp (list (/ (+ x1 x2) 2.0) (/ (+ y1 y2) 2.0)); Midpoint of points cloud and center point of circumcircle through supertriangle.
        r (* (distance cp (list x1 y1)) 20) ;This could still be too small in certain case. No harm if we make it bigger.
       ma (+ (car cp) r);ma is maxX of supertriangle
       mi (- (car cp) r);mi is minX of supertriangle
       sl (list (list ma (cadr cp) 0);list of 3 points defining the supertriangle
(list mi (+ (cadr cp) r) 0)
(list mi (- (cadr cp) r) 0)
  )
 
       al  (list (cons x2 (cons cp (cons (* 20 r) sl))))
      ;al is a work list that contains active triangles each item is a list that contains:      ;
      ;     item 0: Xmax of points in triangle.                         ;
      ;     item 1: List 2d coordinates of center of circle circumscribing triangle.  ;
      ;     item 2: Radius of above circle.                     ;
      ;     item 3: List of 3 points defining the triangle                          ;

       ma (1- ma);Reducing ma to prepare for removal of triangles having a vertex
       mi (1+ mi);common with supertriangle. Done once triangulation is completed.
)               ;Increasing mi for same reason.


;Begin insertion of points

(repeat (length pl)
   
    (setq  p (car pl);Get one point from point list
  pl (cdr pl);Remove point from point list
  el nil     ;Initialize edge list
    )
    (while al ;Active triangle list
       (setq tr (car al);Get one triangle from active triangle list.
     al (cdr al);Remove the triangle from the active list.
       )
       (cond
  ((< (car tr) (car p)) (setq tl (cons (cdddr tr) tl)));This triangle inactive. We store it's 3 vertex in tl (Final triangle list).
  ((< (distance p (cadr tr)) (caddr tr));p is inside the triangle.
   (setq tr (cdddr tr) ;Trim tr to vertex of triangle only.
  a (car tr)   ;  First point.
  b (cadr tr)  ;  Second point.
  c (caddr tr) ;  Third point.
         el (cons (list (+ (car a) (car b))   ;We create 3 new edges ab, bc and ca,
        (+ (cadr a) (cadr b)) ;and store them in edge list.
a
b
  )
  (cons (list (+ (car b) (car c))
      (+ (cadr b) (cadr c))
      b
      c
)
(cons (list (+ (car c) (car a))
    (+ (cadr c) (cadr a))
    c
    a
      )
      el
)
  )
    )
   
   )
  )
  (t (setq l (cons tr l)))   ;tr did not meet any cond so it is still active.
       )      ;we store it a temporary list.
    );Go to next triangle of active list.
   
    (setq al l    ;Restore active triangle list from the temporary list.
   l nil  ;Re-initialize the temporary list to prepare for next insertion.
 
  el (vl-sort el ;Sort the edges list on X and Y
        (function (lambda (a b)
     (if (= (car a) (car b))
        (<= (cadr a) (cadr b))
        (< (car a) (car b))
     )
  )
        )
     )
    )

    ;Removes doubled edges, form new triangles, calculates circumcircles and add them to active list.
    (while el
       (if (and (= (caar el) (caadr el))
(= (cadar el) (cadadr el))
   )
  (setq el (cddr el))
  (setq al (cons (getcircumcircle p (cddar el)) al)
el (cdr el)
  )
       )
    )
    ;Spinner to show progress
    (if (and (< (setq i (1- i)) 1) (< i2 100))
       (progn
  (setvar
     "MODEMACRO"
     (strcat
"     "
(itoa (setq i2 (1+ i2)))
" %    "
(substr
   "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
   1
   i2
)
(substr "..." 1 (- 100 i2))
     )
  )
  (setq i i1)
       )
    )
) ;Go to insert next point.
;We are done with the triangulation

(foreach tr al (setq tl (cons (cdddr tr) tl)));What's left in active list is added to triangle list

;Purge triangle list of any triangle that has a common vertex with supertriangle.
(setq
    tl (vl-remove-if-not
  (function
     (lambda (a)
(and (< mi (caadr a) ma) (< mi (caaddr a) ma))
     )
  )
  tl
       )
)


;Create a layer and Draw the triangulation
(or (tblsearch "LAYER" "TIN")
    (entmake (list
'(0 . "LAYER")
                '(100 . "AcDbSymbolTableRecord")
        '(100 . "AcDbLayerTableRecord")
        '(2 . "TIN")
        '(70 . 0)
        '(62 . 8)
        '(6 . "Continuous")
                        '(290 . 1)
        '(370 . -3)
                     )
             )
)

(setvar "CLAYER" "TIN")

(foreach tr tl
    (entmake (list (cons 0 "3DFACE")
   (cons 10 (car tr))
   (cons 11 (car tr))
   (cons 12 (cadr tr))
   (cons 13 (caddr tr))
     )
    )
)
      )
   )
   (setvar "MODEMACRO" "")
   (princ (strcat "\n "
  (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4)
  " secs."
  )
   )
   (princ)
)


;;*********************************************************;
;;    ;
;; Written by  ElpanovEvgeniy    ;
;; 17.10.2008    ;
;; Calculation of the centre of a circle and circle radius ;
;; for program triangulate    ;
;;    ;
;; Modified ymg june 2011 (renamed variables)    ;
;;*********************************************************;

(defun getcircumcircle (a el /  b c c2 cp r ang)
     
   (setq b (car el)
c (cadr el)
c2 (list (car c) (cadr c)) ;c2 is point c but in 2d
   )
   (if (not
  (zerop
     (setq ang (- (angle b c) (angle b a)))
  )
       )
      (progn (setq cp (polar c2
  (+ -1.570796326794896 (angle c a) ang)
  (setq r (/ (distance a c2) (sin ang) 2.0))
      )
    r (abs r)
     )
     (list (+ (car cp) r) cp r a b c)
      )
   )
)



ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 15, 2011, 01:44:27 AM
Hi ymg! :)

Now, the program still seems complicated?
Title: Re: Triangulation (re-visited)
Post by: ymg on June 15, 2011, 12:11:56 PM
Quote
Now, the program still seems complicated?

No, since it is basically the same algoritmn as Piazza's program.

However I have to agree with Vovka's comment about your variable naming style  :-D

The Divide and conquer algorithm could be faster than this one if written properly.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 15, 2011, 12:30:00 PM
No, since it is basically the same algoritmn as Piazza's program.

Can be a little more detail?
I first heard about this algorithm ...
For this program, I developed an algorithm by yourself!
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on June 15, 2011, 02:00:57 PM
Evgeniy is always Wonderful !!!
This is I learn your Excellent Programe twice , The first is GA for TSP .
Saw the scriptures , rise the praise .

Thanks ymg , Thanks Evgeniy .


Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 15, 2011, 02:14:33 PM
The first is GA for TSP .

I'm sorry, I do not remember what this program. 
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 15, 2011, 02:16:26 PM
TSP = Travelling salesman problem
and
GA  = genetic algorithm
??
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on June 16, 2011, 07:07:58 AM
You're right . I'm sorry for excessive omitted .
Title: Re: Triangulation (re-visited)
Post by: ymg on June 16, 2011, 01:08:41 PM
Quote
Can be a little more detail?
I first heard about this algorithm ...
For this program, I developed an algorithm by yourself!]

For a good explanation of the various way to triangulate go to this link http://www.cs.cmu.edu/~quake/tripaper/triangle2.html (http://www.cs.cmu.edu/~quake/tripaper/triangle2.html)

Divide and conquer algorithm is conceptually the same as a Merge Sort.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on June 16, 2011, 05:53:05 PM
The Divide and conquer algorithm could be faster than this one if written properly.

do you propose to continue the contest in Lisp?
It is a challenge?
Title: Re: Triangulation (re-visited)
Post by: ymg on June 17, 2011, 07:11:15 AM
Quote
do you propose to continue the contest in Lisp?
It is a challenge?

Not really a challenge, like I told you I am very rusted in Lisp and computer programming in general.

Besides the divide and conquer force you to retriangulate everything when you want to add points.

Whereby  with Lawson's insertions you simply added them to what is there.

More useful would be to devise a scheme to add constrain to what we got.

ymg
Title: Re: Triangulation (re-visited)
Post by: TopoWAR on July 24, 2011, 11:37:49 AM
ymg , hello teacher, I found a small problem in your code of triangulation, the matter is repeated points, the question is which is the most accurate and quick delete repeated elements in a list? Thank you.
Title: Re: Triangulation (re-visited)
Post by: ymg on July 26, 2011, 08:06:12 PM
Hi,

Don't know if I understand your question correctly.  But all algorithm for Delaunay triangulation do not allow duplicate point.

Our implementation in Lisp does not check for the presence of duplicate as it is.  The dataset is supposed to have been filtered before and any duplicate removed or slightly perturbed so as not to cause havoc.

Hopes this covers it, and sorry about the late answer.

ymg
Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on July 27, 2011, 01:58:44 AM
We can calculate the convex hull for the remaining points, you can create the opposite direction. Need to reflect on the algorithms...
Title: Re: Triangulation (re-visited)
Post by: ymg on July 02, 2013, 01:31:39 PM
I know the thread is getting older but....

I've modified Evgeniy triangulation to operate on index instead of coordinates.
Speed penalty is minor as far as I can tell.

With the list of index created, I calculate contours.  However speed is not that
great.  I believe the bottleneck to be in all those vl-remove.

Would appreciate any pointers as  to how to accelerate.

There is a lot of comments explaining the way I go about it.

There is a small routine gen which create a bunch of points for testing,


Code - Auto/Visual Lisp: [Select]
  1. (defun c:tin (/ i s)
  2.    (princ (strcat "\n select points"))
  3.    (if (setq i 0
  4.              pl nil
  5.              s (ssget '((0 . "POINT")))
  6.        )
  7.       (progn
  8.         (repeat (sslength s)
  9.              (setq pl   (cons (cdr (assoc 10 (entget (ssname s i)))) pl)
  10.                     i   (1+ i)
  11.              )
  12.         )
  13.         (triangulate pl)
  14.      )
  15.    )
  16. )
  17.  
  18.  
  19. ;;*****************************************************************************;
  20. ;;                                                                                                                                             ;
  21. ;; Structure of Program by  ElpanovEvgeniy                                                                                ;
  22. ;; 17.10.2008                                                                                                                            ;
  23. ;; edit 20.05.2011                                                                                                                     ;
  24. ;; Program triangulate an irregular set of 3d points.                                                                    ;
  25. ;; Modified and Commented by ymg June 2011.                                                                          ;
  26. ;; Modified to operate on index by ymg in June 2013.                                                                 ;
  27. ;; Contour Generation added by ymg in June 2013.                                                                     ;
  28. ;;*****************************************************************************;
  29.  
  30. (defun triangulate (pl  /   a   b   c   i   i1   i2  xmin ymin zmin
  31.                     bb  sl  pl  el  l   zm  z1   z2  xmax ymax zmax                
  32.                     ti  tr  np  n   j   r   cp   sl  pl   al   intv
  33.                     p   tl  vc  pc  e   cl  xl   pol ent  nxt
  34.                    )
  35.    
  36.    (if pl
  37.       (progn
  38.          (setq ti (car (_VL-TIMES));Initialize timer for Triangulation                            ;
  39.                 i  1
  40.                i1 (/ (length pl) 100.)
  41.                i2 0
  42.                ; Variables and Constant to Control Progress Spinner                               ;
  43.                tl nil
  44.                
  45.                pl (vl-sort pl
  46.                      (function (lambda (a b) (< (car a) (car b))))
  47.                   )
  48.                ; Sort points list on x coordinates                                                ;
  49.                
  50.                bb (list (apply 'mapcar (cons 'min pl))
  51.                         (apply 'mapcar (cons 'max pl))
  52.                   )
  53.                ;Replaced code to get the min and max with 3d Bounding Box Routine                 ;
  54.                ;A bit slower but clearer. zmin and zmax kept for contouring                       ;
  55.                
  56.                xmin (caar bb)      
  57.                xmax (caadr bb)      
  58.                ymin (cadar bb)      
  59.                ymax (cadadr bb)      
  60.                zmin (caddar bb)    
  61.                zmax (caddr(cadr bb))
  62.                
  63.                np (length pl) ;Number of points to insert                                         ;
  64.                
  65.                cp (list (/ (+ xmin xmax) 2.0) (/ (+ ymin ymax) 2.0))
  66.                ; Midpoint of points cloud and center point of circumcircle through supertriangle. ;
  67.                 r (* (distance cp (list xmin ymin)) 20)
  68.                ; This could still be too small in certain case. No harm if we make it bigger.     ;
  69.                
  70.                sl (list (list (+ (car cp) r) (cadr cp) 0)          
  71.                         (list (- (car cp) r) (+ (cadr cp) r) 0)    
  72.                         (list (- (car cp) r) (- (cadr cp) r) 0)
  73.                   )
  74.                ; sl list of 3 points defining the Supertriangle,                                  ;
  75.                ; I have tried initializing to an infinite triangle but it slows down calculation  ;
  76.                pl (append pl sl)
  77.                ;Vertex of Supertriangle are appended to the Point list                            ;
  78.                sl (list np (+ np 1)(+ np 2))
  79.                ;sl now is a list of index into point list defining the supertriangle              ;
  80.                
  81.                al  (list(list xmax cp r sl))
  82.               ;Initialize the Active Triangle list                                                ;
  83.               ; al is a  list that contains active triangles defined by 4 items:                  ;
  84.               ;     item 0: Xmax of points in triangle.                                           ;
  85.               ;     item 1: List 2d coordinates of center of circle circumscribing triangle.      ;
  86.               ;     item 2: Radius of above circle.                                               ;
  87.               ;     item 3: List of 3 indexes to vertices defining the triangle                   ;
  88.                
  89.                n -1
  90.               ; n is a counting index into Point List                                             ;
  91.          )              
  92.  
  93.          
  94.          ;Begin insertion of points
  95.          
  96.          (repeat np
  97.            
  98.             (setq  n (1+ n)     ; Increment Index into Point List                                 ;
  99.                    p (nth n pl) ; Get one point from point list                                   ;
  100.                   el nil        ; el list of triangles edges                                      ;
  101.             )                   ;                                                                 ;
  102.             (repeat (length al) ; Loop to go through Active triangle list                         ;
  103.                (setq tr (car al); Get one triangle from active triangle list.                     ;
  104.                      al (cdr al); Remove the triangle from the active list.                       ;
  105.                )
  106.                (cond
  107.                   ((< (car tr) (car p)) (setq tl (cons (cadddr tr) tl)))
  108.                   ;This triangle inactive. We store it's 3 vertex in tl (Final triangle list).    ;
  109.                  
  110.                   ((< (distance p (cadr tr)) (caddr tr)) ; p is inside the triangle.                  ;
  111.                    (setq tr (cadddr tr)          ; Trim tr to vertex of triangle only.   ;
  112.                           a (car tr)                             ;  Index of First point.                    ;
  113.                           b (cadr tr)                           ;  Index of Second point.                ;
  114.                           c (caddr tr)                         ;  Index of Third point.                   ;
  115.                          el (cons (list a                            ; ((a b)(b c)(c a) (. .)(. .).....)      ;
  116.                                         b
  117.                                   )
  118.                                   (cons (list b
  119.                                               c
  120.                                         )
  121.                                         (cons (list c
  122.                                                     a
  123.                                               )
  124.                                               el
  125.                                         )
  126.                                   )
  127.                             )
  128.                            
  129.                    )
  130.                   )
  131.                   (t (setq l (cons tr l)))
  132.                   ;tr did not meet any cond so it remain active. We store it in the swap list     ;
  133.                );End cond
  134.              
  135.             );End repeat, go to next triangle of active list.
  136.            
  137.            
  138.             (setq al l    ;Restore active triangle list from the temporary list.                  ;
  139.                    l nil  ;Clear the swap list to prepare for next insertion.                     ;
  140.             )
  141.            
  142.             ;Removes doubled edges, calculates circumcircles and add them to al                   ;
  143.             (while el
  144.                (if (or (member (reverse (car el)) el)
  145.                        (member (car el) (cdr el))
  146.                    )
  147.                    (setq el (vl-remove (reverse (car el)) el)
  148.                          el (vl-remove (car el) el)
  149.                    )
  150.                    (setq al (cons (getcircumcircle n (car el)) al)
  151.                          el (cdr el)
  152.                   )
  153.                )
  154.             )
  155.            
  156.             ; Neat Spinner to show progress does not work too good with Window7                   ;
  157.             (if (and (< (setq i (1- i)) 1) (< i2 100))
  158.                (progn
  159.                   (setvar
  160.                      "MODEMACRO"
  161.                      (strcat
  162.                         "     "
  163.                         (itoa (setq i2 (1+ i2)))
  164.                         " %    "
  165.                         (substr
  166.                            "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||"
  167.                            1
  168.                            i2
  169.                         )
  170.                         (substr "..." 1 (- 100 i2))
  171.                      )
  172.                   )
  173.                   (setq i i1)
  174.                )
  175.             )
  176.            
  177.          ) ;End repeat np, Go to insert next point
  178.          
  179.        
  180.          ;We are done with points insertion. Any triangle left in al is added to tl               ;
  181.          
  182.          (foreach tr al (setq tl (cons (cadddr tr) tl)))
  183.          
  184.          ;Purge triangle list of any triangle that has a common vertex with supertriangle.        ;
  185.          (setq
  186.             tl (vl-remove-if-not
  187.                   (function
  188.                      (lambda (a)
  189.                         (and (< (car a) np)(< (cadr a) np)(< (caddr a) np))
  190.                      )
  191.                   )
  192.                   tl
  193.                )
  194.          )
  195.  
  196.  
  197.          ;Create a layer and Draw the triangulation                                               ;
  198.          (or (tblsearch "LAYER" "TIN")
  199.             (entmake (list
  200.                         '(0 . "LAYER")
  201.                         '(100 . "AcDbSymbolTableRecord")
  202.                         '(100 . "AcDbLayerTableRecord")
  203.                         '(2 . "TIN")
  204.                         '(70 . 0)
  205.                         '(62 . 8)
  206.                         '(6 . "Continuous")
  207.                         '(290 . 1)
  208.                         '(370 . -3)
  209.                      )
  210.              )
  211.          )
  212.          
  213.          (setvar "CLAYER" "TIN")
  214.          
  215.          (foreach tr tl
  216.             (entmake (list (cons 0 "3DFACE")
  217.                            (cons 10 (nth (car tr)   pl))
  218.                            (cons 11 (nth (car tr)   pl))
  219.                            (cons 12 (nth (cadr tr)  pl))
  220.                            (cons 13 (nth (caddr tr) pl))
  221.                      )
  222.             )
  223.          )
  224.       )
  225.    )
  226.    (setvar "MODEMACRO" "")
  227.    
  228.    (princ (strcat "\n     TIN Completed - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) " secs."))
  229.    (princ (strcat "\n         Generated   " (itoa (length tl)) " 3DFACES"))
  230.    (princ "\n")
  231.  
  232.  
  233.    ; Begin Calculation of Contour line                                                            ;
  234.    ;                                                                                              ;
  235.    ; The routine follows each contour from edge to edge resulting in Connected LWPOLYLINE.        ;
  236.    ;                                                                                              ;
  237.    ; 1. Sort triangle list tl on the maximum z value of vertices of each triangles.               ;
  238.    ; 2. Create el, a list containing all edges of all triangles. ((a b)(b c) (c a).......)        ;
  239.    ; 3. For each desired contour level l, we traverse el and  create list cl containing all       ;
  240.    ;    edges crossed by level l. At this point cl contains 2 edges per triangle.                 ;
  241.    ;    As we go through el, we can destroy any triplets of edges whose max z value is  below     ;
  242.    ;    the current contour level, thus limiting traversal for next level.                        ;
  243.    ; 4. We now process cl, first locating an edge with no neighbour if any and follow from edge   ;
  244.    ;    to edge until either we've closed the polyline or reached the convex hull.                ;
  245.    ;    As we complete each polyline list xl is formed.                                           ;
  246.    ; 5. We entmake each element of list xl thus completing.                                       ;
  247.    ;    Step 4 and 5 could be combined but it is easier to follow contour in index form           ;
  248.    ;                                                                                              ;
  249.    ;                                                                                              ;
  250.    ; An alternate way to do this would be compute all al segment between two edges joining with   ;
  251.    ; with a line for all contour level and at end Join everything together.                       ;
  252.    ;                                                                                              ;
  253.    
  254.    (setq tl (vl-sort tl
  255.                 (function
  256.                     (lambda (a b)
  257.                         (< (max (caddr (nth (car a) pl))
  258.                                 (caddr (nth (cadr a) pl))   ;Gotta be a more concise way         ;
  259.                                 (caddr (nth (caddr a) pl))  ;to write this probably with apply.  ;
  260.                            )
  261.                            (max (caddr (nth (car b) pl))
  262.                                 (caddr (nth (cadr b) pl))
  263.                                 (caddr (nth (caddr b) pl))
  264.                            )
  265.                         )
  266.                     )
  267.                 )
  268.             )
  269.    )
  270.  
  271.    ; Setup for Contouring                                                                         ;
  272.    (setq   ti (car (_VL-TIMES))   ; Re-initialize timer for Contouring                            ;
  273.          intv 1                   ; Interval between contour                                      ;
  274.          zmin (+ (fix zmin) intv) ; zmin was calculated during triangulation                      ;
  275.          zmax (fix zmax)          ; z2 also calculated at beginning                               ;
  276.             l zmin                ; Initialize Contour level                                      ;
  277.            el nil                 ; el will be list of all edges                                  ;
  278.            vc 0                   ; Vertices Count                                                ;
  279.            pc 0                   ; LWPOLYLINE Count                                              ;
  280.    )
  281.  
  282.    ; Extract all triangle edges from tl and form list el                                          ;
  283.    (foreach tr tl
  284.      (setq el (cons (list (car tr)(cadr tr)) el)
  285.            el (cons (list (cadr tr)(caddr tr)) el)
  286.            el (cons (list (caddr tr)(car tr)) el)
  287.      )
  288.    )
  289.                  
  290.    
  291.    (repeat (+(fix (/ (- zmax zmin) intv)) 1) ;Main Loop through each contour level                ;
  292.      
  293.      (setq cl nil   ; cl will be list of all edges crossed at current level l                     ;
  294.             j 0     ; Index into edge list el                                                     ;
  295.            zm 1e-8  ; zmax value for a given triplets of edges                                    ;
  296.      )
  297.      (repeat (length el)
  298.        (setq  e (nth j el)
  299.              z1 (caddr (nth (car e) pl))  ; Get elevation of edges from the point list.           ;
  300.              z2 (caddr (nth (cadr e) pl))
  301.              zm (max zm z1 z2)  
  302.               j (1+ j)
  303.        )
  304.        
  305.        (if (and (= (rem j 3) 0) (< zm (+ l intv))) ; Reduce size of el on zmax criteria.          ;
  306.              (setq  j (- j 3)
  307.                    el (vl-remove (nth j el) el)
  308.                    el (vl-remove (nth j el) el)
  309.                    el (vl-remove (nth j el) el)
  310.                    zm 1e-8
  311.              )
  312.        )
  313.              
  314.        (if (= z1 l) (setq z1 (- z1 1e-8))); If vertex is equal to l we disturb                    ;
  315.        (if (= z2 l) (setq z2 (- z2 1e-8))); the z value a little.                                 ;
  316.          
  317.        (if (or (< z1 l z2)(> z1 l z2))    
  318.           (setq cl (cons e cl))           ; Edge is added to Crossed List                         ;
  319.        )
  320.      );end foreach e
  321.      
  322.      ; cl now contains all edges where all contours at level l passes.                            ;
  323.  
  324.      (setq xl nil)      
  325.  
  326.      (while cl
  327.         (setq n 0)
  328.         ;;Find in cl an edge that has no reverse hence on the convex hull and start from it       ;
  329.         ;;if none is found we start with the last edge in cl                                      ;
  330.         (while (and (member (reverse (nth n cl)) cl) (< n (1- (length cl))))
  331.            (setq n (1+ n))
  332.         )
  333.      
  334.         (setq pol (list (nth n cl))
  335.                cl (vl-remove (nth n cl) cl)
  336.                 n (- n (rem n 2))
  337.               pol (cons (nth n cl) pol)
  338.               nxt (reverse (nth n cl))
  339.                cl (vl-remove (nth n cl) cl)
  340.                pc (1+ pc)
  341.                vc (+ 2 vc)
  342.         )
  343.        
  344.        (while (and (setq n (vl-position nxt cl)) (not (member nxt pol)))
  345.          
  346.            (setq  cl (vl-remove (nth n cl) cl)
  347.                    n (- n (rem n 2))
  348.                  pol (cons (nth n cl) pol)
  349.                  nxt (reverse (car pol))
  350.                   cl (vl-remove (nth n cl) cl)
  351.                   vc (1+ vc)
  352.            )
  353.          
  354.        );end while
  355.        
  356.        (setq xl (cons pol xl))
  357.        
  358.      );end while cl
  359.      
  360.      (foreach p xl
  361.            (setq ent nil)
  362.            (foreach e p
  363.               (setq p1 (nth (car e) pl)
  364.                     p2 (nth (cadr e) pl)
  365.                      r (/ (- l (caddr p1)) (- (caddr p2) (caddr p1)))
  366.                     p1 (list (car p1)(cadr p1))
  367.                     p2 (list (car p2)(cadr p2))
  368.                      d (* (distance p1 p2) r)
  369.                     pt (polar p1 (angle p1 p2) d)
  370.                    ent (cons (cons 10 pt) ent)
  371.               )
  372.            )
  373.            (setq ent (cons (cons 38 l) ent)
  374.                  ent (cons (cons 43 0.0) ent)
  375.                  ent (cons (cons 70 0) ent)
  376.                  ent (cons (cons 90 (length p)) ent)
  377.                  ent (cons (cons 100 "AcDbPolyline") ent)
  378.                  ent (cons (cons 100 "AcDbEntity") ent)
  379.                  ent (cons (cons 8 "Contour") ent)
  380.                  ent (cons (cons 0 "LWPOLYLINE") ent)
  381.  
  382.                  
  383.            )
  384.  
  385.            (entmake ent)
  386.      );end foreach p
  387.      
  388.      (setq l (+ l intv))
  389.      
  390.    );end repeat
  391.  
  392.    (princ (strcat "\n CONTOUR Completed - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) " secs."))
  393.    (princ (strcat "\n         Generated   " (itoa pc) " LWPOLYLINE."))
  394.    (princ (strcat "\n             Total   " (itoa vc) " Vertices."))     
  395.    (princ)
  396.  
  397. );end defun triangulate
  398.  
  399.  
  400. ;;*****************************************************************************;
  401. ;;                                                                             ;
  402. ;; Written by  ElpanovEvgeniy                                                  ;
  403. ;; 17.10.2008                                                                  ;
  404. ;; Calculation of the centre of a circle and circle radius                     ;
  405. ;; for program triangulate                                                     ;
  406. ;;                                                                             ;
  407. ;; Modified ymg june 2011 (renamed variables)                                  ;
  408. ;; Modified ymg June 2013 to operate on Index                                  ;
  409. ;;*****************************************************************************;
  410.  
  411. (defun getcircumcircle (a el /  b c c2 cp r ang vl)
  412.      
  413.    (setq p (nth a pl)
  414.          b (nth(car el) pl)
  415.          c (nth(cadr el) pl)
  416.         c2 (list (car c) (cadr c)) ;c2 is point c but in 2d                    :
  417.         vl (list a (car el) (cadr el))
  418.    )
  419.    (if (not
  420.           (zerop
  421.              (setq ang (- (angle b c) (angle b p)))
  422.           )
  423.        )
  424.       (progn (setq cp (polar c2
  425.                           (+ -1.570796326794896 (angle c p) ang)
  426.                           (setq r (/ (distance p c2) (sin ang) 2.0))
  427.                       )
  428.                     r (abs r)
  429.              )
  430.              (list (+ (car cp) r) cp r vl)
  431.       )
  432.    )
  433. )
  434.  
  435. ;Generate a bunch of points for testing on layer named points.                 ;
  436. (defun c:gen ( / n)
  437.        
  438.         (setq n (getint "\nNumber of points: "))
  439.         (while (> n 0)
  440.                 (entmake
  441.                     (list
  442.                        '(0 . "POINT")
  443.                        '(8 . "Points")
  444.                         (cons 10 (list (* 1000 (ran)) (* 750 (ran)) (* 10 (ran))))
  445.                     )
  446.                 )
  447.                 (setq n (1- n))
  448.         )
  449. )
  450. ; Random number generator                                                      ;
  451. (defun ran ()
  452.     (setq seed (if seed (rem (+ (* seed 15625.7) 0.21137152) 1) 0.3171943))
  453. )
  454.  
  455.  
  456.  
Title: Re: Triangulation (re-visited)
Post by: ymg on July 02, 2013, 05:37:44 PM
I am replying to myself here.

Upon close examination the biggest speed penalty is on the search for edges on the convex hull.
The list cl is scanned until the end at every edges.

I will modify the routine.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on July 03, 2013, 10:03:57 PM
I've modified the routine and am getting much better speed now.

Also modified the entmake part and insure that closed polylines are built accordingly.

Would still appreciate any speed tips or better ways to write some of the segment.

Next revision will implement smoothing of contour.  As you know simply splining
the polylines leads to contour crossing each other.

There is a way that was patented back in 1992 by IBM, I guess the patent has ran out now.
The method was developped by a fellow name Albert H. J.  Christensen and involves
inserting parabola at vertices of the polylines when a chosen angular threshold is reached.

If you do a search on "Contour Smoothing by an Eclectic Procedure" you will find the document.

Another thing we would need is a way to insert constraints (Breaklines) into the triangulation.

Meanwhile find below the revised code for contouring as well as the full lisp including
the triangulation as an attachment.

ymg

Code: [Select]
(grread t); Dummy just to force a screen update and avoid that Autocad hangs while calculating
 
   (princ (strcat "\n     TIN Completed - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) " secs."))
 
   (princ (strcat "\n         Generated   " (itoa (length tl)) " 3DFACES"))
   (princ "\n")


   ; Begin Calculation of Contour line                                                             
   ;                                                                                               
   ; The routine follows each contour from edge to edge resulting in Connected LWPOLYLINE.         
   ;                                                                                               
   ; 1. Sort triangle list tl on the maximum z value of vertices of each triangles.               
   ; 2. Create el, a list containing all edges of all triangles. ((a b)(b c) (c a).......)         
   ; 3. For each desired contour level l, we traverse el and  create list cl containing all       
   ;    edges crossed by level l. At this point cl contains 2 edges per triangle.                 
   ;    As we go through el, we can destroy any triplets of edges whose max z value is  below     
   ;    the current contour level, thus limiting traversal for next level.                         
   ; 4. We now process cl, following from edge to edge until either we've closed the polyline     
   ;    or reached the convex hull. If we reached the convex hull, we insure  that the beginning 
   ;    of the polyline is also on the convex hull.  Otherwise we reverse the polyline and         
   ;    continue processing.  Process continues until cl is empty.                           
   ;    As we complete each polyline list xl is formed.                                           
   ; 5. We entmake each element of list xl.
   ; 6. We go back to step 3 until all levels are completed and el is empty.
   ;                                       
   ;    Step 4 and 5 could be combined but it is easier to follow contour in index form           
   ;                                                                                               
   ;                                                                                               
   ; An alternate way to do this would be compute all al segment between two edges joining with   
   ; with a line for all contour level and at end Join everything together.                       
   ;                                                                                               


   ; Setup for Contouring                                                                         
   (setq   ti (car (_VL-TIMES))   ; Re-initialize timer for Contouring                             
     i  1
    i1 (/ (length pl) 100.)
    i2 0
         intv 1                   ; Interval between contour                                       
zmin (+ (fix zmin) intv) ; zmin was calculated during triangulation                       
zmax (fix zmax)          ; z2 also calculated at beginning                               
            l zmin                ; Initialize Contour level                                       
   el nil                 ; el will be list of all edges                                   
   vc 0                   ; Vertices Count                                                 
   pc 0                   ; LWPOLYLINE Count                                               
   )
 
   (setq tl (vl-sort tl
(function
    (lambda (a b)
        (< (max (caddr (nth (car a) pl))
(caddr (nth (cadr a) pl))   ;Gotta be a more concise way           
(caddr (nth (caddr a) pl))  ;to write this probably with apply.   
   )                                ;Help!!!                               
   (max (caddr (nth (car b) pl))
(caddr (nth (cadr b) pl))
(caddr (nth (caddr b) pl))
   )
)
    )
        )
            )
   )
 
 
   
   ; Extract all triangle edges from tl and form list el                                           
   (foreach tr tl
     (setq el (cons (list (car tr)(cadr tr)) el)
   el (cons (list (cadr tr)(caddr tr)) el)
   el (cons (list (caddr tr)(car tr)) el)
     )
   )

   
   
   (repeat (+(fix (/ (- zmax zmin) intv)) 1) ;Main Loop through each contour level                 
     
     (setq cl nil   ; cl will be list of all edges crossed at current level l                     
            j 0     ; Index into edge list el                                                     
           zm 1e-8  ; zmax value for a given triplets of edges                                     
     )
     (repeat (length el)
       (setq  e (nth j el)
     z1 (caddr (nth (car e) pl))  ; Get elevation of edges from the point list.           
     z2 (caddr (nth (cadr e) pl))
     zm (max zm z1 z2) 
      j (1+ j)
       )
       
       (if (and (= (rem j 3) 0) (< zm (+ l intv))) ; Reduce size of el on zmax criteria.           
             (setq  j (- j 3)
                   el (vl-remove (nth j el) el)
                   el (vl-remove (nth j el) el)
                   el (vl-remove (nth j el) el)
           zm 1e-8
             )
       )
             
       (if (= z1 l) (setq z1 (- z1 1e-8))); If vertex is equal to l we disturb                     
       (if (= z2 l) (setq z2 (- z2 1e-8))); the z value a little.                                 
       
       (if (or (< z1 l z2)(> z1 l z2))   
          (setq cl (cons e cl))           ; Edge is added to Crossed List                         
       )
     );end foreach e
     
     ; cl now contains all edges where all contours at level l passes.                             

     (setq xl nil)     

     (while cl
       ;We Initialize a Polyline                                                             
       (setq pol (list (cadr cl)(car cl)) ; We go reverse as we will cons the polyline             
             nxt (reverse (car pol))      ; nxt will be our next edge                             
      cl (cddr cl)                ; Remove first two edges from cl                         
       )
       (if (not (member nxt cl))          ;The previous edge was on convex hull                   
    (setq pol (reverse pol)       ;We reverse our Polyline                                 
  nxt (reverse(car pol))  ;and adjust our next edge                               
            )
       )
       
       (while (setq n (vl-position nxt cl))

           (setq  cl (vl-remove (nth n cl) cl)
                   n (- n (rem n 2))
                 pol (cons (nth n cl) pol)
  cl (vl-remove (nth n cl) cl)
   )
   (if (member nxt pol)
      (setq nxt nil)
      (setq nxt (reverse (car pol)))
   )

   (if (not (vl-position nxt cl))
      (setq pol (reverse pol)
            nxt (reverse (car pol))
      )
   ) 
       );end while
       
       (setq xl (cons pol xl))
       
     );end while cl
     
     (setq pc (+ pc (length xl)))
     (foreach p xl
           (setq ent nil
                  vc (+ vc (length p))
           )
           (if (equal (car p) (reverse (last p)))
         (setq isclosed 1
              p (cddr p)
)     
         (setq isclosed 0)
   )
           (foreach e p
(setq p1 (nth (car e) pl)
       p2 (nth (cadr e) pl)
        r (/ (- l (caddr p1)) (- (caddr p2) (caddr p1)))
       p1 (list (car p1)(cadr p1))
       p2 (list (car p2)(cadr p2))
                d (* (distance p1 p2) r)
       pt (polar p1 (angle p1 p2) d)
              ent (cons (cons 10 pt) ent)
       
                )
           )
           (setq ent (cons (cons 38 l) ent)
         ent (cons (cons 43 0.0) ent)
         ent (cons (cons 70 isclosed) ent)
         ent (cons (cons 90 (length p)) ent)
         ent (cons (cons 100 "AcDbPolyline") ent)
         ent (cons (cons 100 "AcDbEntity") ent)
         ent (cons (cons 8 "Contour") ent)
         ent (cons (cons 0 "LWPOLYLINE") ent)


           )

           (entmake ent)
           
     );end foreach p
     
     (setq l (+ l intv))
     
   );end repeat

   (setvar "MODEMACRO" "")
 
   (princ (strcat "\n CONTOUR Completed - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) " secs."))
   (princ (strcat "\n         Generated   " (itoa pc) " LWPOLYLINE."))
   (princ (strcat "\n             Total   " (itoa vc) " Vertices."))  
   (princ)
 
);end defun triangulate



Title: Re: Triangulation (re-visited)
Post by: ElpanovEvgeniy on July 04, 2013, 04:51:32 AM
Hello ymg!
See my version of the smoothing algorithm (http://www.theswamp.org/index.php?topic=42846.0;all)
Title: Re: Triangulation (re-visited)
Post by: ymg on July 04, 2013, 10:59:36 AM
Hello Evgeniy,

I've look into your algorythm.  Sure did appreciate the subject of the study.
Gives a whole new meaning to Model Space  :-)

Your process as far as I can tell is patch based (As opposed to line smoothing)
and applicable to matrix. 

If you look into Christensen's paper,

   http://www.asprs.org/a/publications/pers/2001journal/april/2001_apr_511-517.pdf (http://www.asprs.org/a/publications/pers/2001journal/april/2001_apr_511-517.pdf)

you will get some of the pluses and minuses of this method.

He then goes on proposing a so called eclectic method which actually is kind of an hybrid (line and patch)
solutions.

This is what I propose to implement.

Hope everything is well with you.


ymg

Title: Re: Triangulation (re-visited)
Post by: ymg on July 04, 2013, 12:54:10 PM
I inadvertently introduce a bug in the last posting of tin.lsp,    :ugly:
as I checked for isclosed polyline the list should be reduced by one element
instead of two.

Apologies !!

ymg

Code: [Select]
(foreach p xl
           (setq ent nil
                  vc (+ vc (length p))
           )
           (if (equal (car p) (reverse (last p)))
         (setq isclosed 1
              p (cdr p) ;Was wrong (cddr p)
)     
         (setq isclosed 0)
   )
Title: Re: Triangulation (re-visited)
Post by: MP on July 04, 2013, 12:57:46 PM
I inadvertently introduce a bug in the last posting of tin.lsp,    :ugly:

Unacceptable, 3 lashes with a wet hypotenuse.
Title: Re: Triangulation (re-visited)
Post by: ymg on July 20, 2013, 07:07:31 PM
Here is Contour with Christensen's Eclectic Method for smoothing implemented.

I still have not implemented the variable smoothing that he describe in his paper.

Would appreciate comments, bug(s) report and suggestion for speeding up.

Code: [Select]
; Contour by ymg   (July 2013)                                                                     
;                                                                                                 
; The routine follows each contour from edge to edge resulting in Connected LWPOLYLINE.           
; Basic method is per:                                                                             
;      http://www.originlab.com/www/helponline/Origin/en/UserGuide/Creating_Contour_Graphs.html   
;                                                                                                 
; Contours are then smoothed using Christensen's Eclectic Method.                                 
; For details see:                                                                                 
;      http://www.asprs.org/a/publications/pers/2001journal/april/2001_apr_511-517.pdf             
;      http://www.google.com.mx/patents/US5333248                                                 
;                                                                                                 
; General Flow:                                                                                   
;                                                                                                 
; 1. Sort triangle list tl on the maximum z value of vertices of each triangles.                   
;                                                                                                 
; 2. Create el, a list containing all edges of all triangles. ((a b)(b c) (c a).......)           
;                                                                                                 
; 3. For each desired contour level l, we traverse el and  create list cl containing all           
;    edges crossed by level l. At this point cl contains 2 edges per triangle.                     
;    As we go through el, we can destroy any triplets of edges whose max z value is  below         
;    the current contour level, thus limiting traversal for next level.                           
;                                                                                                 
; 4. We now process cl, following from edge to edge until either, we've closed the polyline       
;    or reached the convex hull. If we reached the convex hull, we insure  that the beginning     
;    of the polyline is also on the convex hull.  Otherwise we reverse the polyline and           
;    continue processing until cl is empty.                                         
;    As we complete each polyline list xl is formed and contains a series of list tracing all the 
;    polylines for a given level.                                                                 
;                                                                                                 
; 5. We entmake each element of list xl, smoothing vertices of polylines when conditions are right.
;    Variable mang controls if an edge can be smoothed or not preventing crossing of contours.     
;    Smoothing is controlled by two variables lang and lres which limits the number of segments   
;    tracing the parabolas of an edge.                                                             
;                                                                                                 
; 6. We go back to step 3 until all levels are completed and el is empty.                         
;                                                                                                 
;                                                                                                 
;                                                                                                 

 
(defun contour (tl pl / a a1 a2 a3 b bb c cl cn1 cn2 e el ent hfac i i1 i2 intv isclosed
        j l lang lres mang n nxt p1 p2 p3 p4 pc pl pol prv seg sm ti
        v1 v2 v3 vcr vcs xl z1 z2 zm zmax zmin
       )
 
   ; Setup for Contouring                                                                         
   (setq    ti (car (_VL-TIMES))   ; Re-initialize timer for Contouring                           
     i  1
    i1 (/ (length pl) 100.)
    i2 0
    bb (list
                             (apply 'mapcar (cons 'min pl))
                (apply 'mapcar (cons 'max pl))
                         )                           ; bb, bounding box of Point List.                               
  zmin (caddar bb)            ;  zmin, minimum z of Point List.                               
  zmax (caddr(cadr bb))    ; zmax, maximum z of Point List.                               
  intv 1                             ; Interval between contour.                                     
                zmin (+ (fix zmin) intv)  ; zmin, first level for contour.                               
  zmax (fix zmax)             ; zmax; last level for contour.                                 
        l zmin                      ; Initialize Current Contour level.                             
   intv 1                            ; Interval between contour.                                     
     el nil                           ; el, will be list of all edges in triangle list.               
   vcr 0                             ; Vertices Count for Raw Contours.                             
   vcs 0                             ; Vertices Count for Smoothed Contours.                         
    pc 0                              ;  LWPOLYLINE Count.                                             
  lres 5                              ; Linear resolution for smoothing.                             
  lang (dtr 5)                      ; Angular resolution for smoothing.                             
  mang (* lang 1.5)           ; Min angle controls if an edge is suitable for smoothing.     
  hfac 0.5                         ; Degree of smoothing, Max is 0.5, Min would be 1.             
   )
 
   (setq tl (vl-sort tl
(function
    (lambda (a b)
        (< (max (caddr (nth (car a) pl))
(caddr (nth (cadr a) pl))   ;Gotta be a more concise way           
(caddr (nth (caddr a) pl))  ;to write this probably with apply.   
   )                                ;Help!!!                               
   (max (caddr (nth (car b) pl))
(caddr (nth (cadr b) pl))
(caddr (nth (caddr b) pl))
   )
)
    )
        )
            )
   )
 
 
   
   ; Extract all triangle edges from tl and form edges list el                                     
   (foreach tr tl
     (setq el (cons (list (car tr)(cadr tr)) el)
   el (cons (list (cadr tr)(caddr tr)) el)
   el (cons (list (caddr tr)(car tr)) el)
     )
   )

   
   
   (repeat (+(fix (/ (- zmax zmin) intv)) 1) ;Main Loop through each contour level                 
     
     (setq cl nil   ; cl will be list of all edges crossed at current level l                     
            j 0     ; Index into edge list el                                                     
           zm 1e-8  ; zmax value for a given triplets of edges                                     
     )
     (repeat (length el)
       (setq  e (nth j el)
     z1 (caddr (nth (car e) pl))  ; Get elevation of edges from the point list.           
     z2 (caddr (nth (cadr e) pl))
     zm (max zm z1 z2) 
      j (1+ j)
       )
       
       (if (and (= (rem j 3) 0) (< zm (+ l intv))) ; Reduce size of el on zmax criteria.           
             (setq  j (- j 3)
                   el (vl-remove (nth j el) el)
                   el (vl-remove (nth j el) el)
                   el (vl-remove (nth j el) el)
           zm 1e-8
             )
       )
             
       (if (= z1 l) (setq z1 (- z1 1e-8))); If vertex is equal to l we disturb                     
       (if (= z2 l) (setq z2 (- z2 1e-8))); the z value a little.                                 
       
       (if (or (< z1 l z2)(> z1 l z2))   
          (setq cl (cons e cl))           ; Edge is added to Crossed List                         
       )
     );end foreach e
     
     ; cl now contains all edges where all contours at level l passes.                             

     (setq xl nil)     

     (while cl
       ;We Initialize a Polyline                                                             
       (setq pol (list (cadr cl)(car cl)) ; We go reverse as we will cons the polyline             
             nxt (reverse (car pol))      ; nxt will be our next edge                             
      cl (cddr cl)                ; Remove first two edges from cl                         
       )
       (if (not (member nxt cl))          ;The previous edge was on convex hull                   
    (setq pol (reverse pol)       ;We reverse our Polyline                                 
  nxt (reverse(car pol))  ;and adjust our next edge                               
            )
       )
       
       (while (setq n (vl-position nxt cl))

           (setq  cl (vl-remove (nth n cl) cl)
                   n (- n (rem n 2))
                 pol (cons (nth n cl) pol)
  cl (vl-remove (nth n cl) cl)
   )
   (if (member nxt pol)
      (setq nxt nil)
      (setq nxt (reverse (car pol)))
   )

   (if (not (vl-position nxt cl))
      (setq pol (reverse pol)
            nxt (reverse (car pol))
      )
   ) 
       );end while
       
       (setq xl (cons pol xl))
       
     );end while cl
     
     (setq pc (+ pc (length xl)))
     
     (foreach p xl
        (setq ent nil)
        (if (equal (car p) (reverse (last p)))
         (setq isclosed 1
              p (append p (list(cadr p)))
)
(setq isclosed 0
                            ent (list (cons 10 (clv l (car p) pl)))
)
        )
        (while (> (length p) 2)
           
   (setq  v1 (clv l (car p) pl)
          v2 (clv l (cadr p) pl)
   )
     
   (setq  v3 (clv l (caddr p) pl)
prv (car p)
nxt (caddr p)
           p (cdr p)
          p1 (nth (caar p) pl)
          p3 (nth (cadar p) pl)
  p2 (nth (car prv) pl)
  p4 (nth (car nxt) pl)
   )
           (if (or (equal p2 p1) (equal p2 p3))
         (setq p2 (nth (cadr prv) pl))
   )
   (if (or (equal p4 p1) (equal p4 p3))
         (setq p4 (nth (cadr nxt) pl))
   )

   ; Need to test here if edge p1 p3 is suitable for smoothing                             
   ; Unclear to me, what mang should if I set it to same value as lang, we still have     
           ; Contours crossing each other.  With (* lang 1.5) no crossing!!                       
   (cond
     ((and (or (< (defl p1 p2 p3) mang) (< (defl p1 p4 p3) mang))
           (> (/ (abs (- (caddr p1) (caddr p3))) intv) 1)
      )
                   (setq sm (list v2 )); Unsuitable we set v2 as the vertices on this edge.       
     )
     (t (setq cn1 (mapcar '(lambda (a b c) (/ (+ a b c) 3.0)) p1 p2 p3)
      cn2 (mapcar '(lambda (a b c) (/ (+ a b c) 3.0)) p1 p3 p4)
       a1 (cond
     ((inters cn1 p1 v1 v2))
     ((inters cn1 p3 v1 v2))
    ;((inters cn1 p2 v1 v2))
  )

       a2 v2
                       a3 (cond
     ((inters cn2 p1 v2 v3))
     ((inters cn2 p3 v2 v3))
    ;((inters cn2 p4 v2 v3))
  )
       sm (smooth a1 a2 a3); Probably would be better not to call a function here 
       ); end setq
     ); end t   
   ); end cond
 
   (foreach e sm
      (setq ent (cons (cons 10 e) ent))
   )
       
        ); end while p
       
        (if (= isclosed 0) (setq ent (cons(cons 10 v3) ent)))

        (setq seg (length ent)
      vcs (+ vcs seg)
      ent (cons (cons 38 l) ent)
      ent (cons (cons 43 0.0) ent)
      ent (cons (cons 70 isclosed) ent)
      ent (cons (cons 90 seg) ent)
      ent (cons (cons 100 "AcDbPolyline") ent)
      ent (cons (cons 100 "AcDbEntity") ent)
      ent (cons (cons 8 "Contour") ent)
      ent (cons (cons 0 "LWPOLYLINE") ent)


        )

        (entmake ent)
         
     );end foreach p
     
     (setq l (+ l intv))
     
   );end repeat

   (setvar "MODEMACRO" "")
   
   (princ (strcat "\n CONTOUR - Elapsed time: " (rtos (/ (- (car (_VL-TIMES)) ti) 1000.) 2 4) " secs, "
                   (itoa pc) " LWPOLYLINES, "
                   (itoa vcs) " Vertices."
          )
   )
   (princ)
 
);end defun contour



; Degree to Radian Conversion                                                                     
(defun dtr (a) (* pi (/ a 180.0)))


; Given:   lev, (level of contour being traced).                                                   
;           ed, (edge defined by a list of two index into the point list.,                         
;           pl, Point List                                                                         
;                                                                                                 
; Returns: Point where contour cross the edge.                                                     

(defun clv (lev ed pl /  r d p1 p2)
(setq p1 (nth (car ed) pl)
       p2 (nth (cadr ed) pl)
        r (/ (- lev (caddr p1))
             (- (caddr p2) (caddr p1))
  )
       p1 (list (car p1)(cadr p1))
       p2 (list (car p2)(cadr p2))
        d (* (distance p1 p2) r)
)
(polar p1 (angle p1 p2) d)
)



; Smooth an edge defined by 3 points                                                               
; Returns a list of points starting at a1 and ending at a3                                         
; which traces two parabolas joining a1 to a3                                                     
; Global variables lang and lres limits number of segments on parabolas.                           
; Global variable hfac controls how much smoothing we apply  (Not implemented yet)                 

(defun smooth (a1 a2 a3 / cp f ps tmp)
   (setq ps (subdiv a1 a2 a3)
          f nil
   )
(while ps
     (setq cp (car ps)
   ps (cdr ps)
     )
    (while (and (> (min (distance (car cp) (cadr cp))
                (distance (cadr cp) (caddr cp))
   )
   lres
                        )
(< (defl (car cp) (cadr cp) (caddr cp))
                           (- pi lang)
                        )
                   )
   
   (setq tmp (subdiv (caddr cp)(cadr cp)(car cp))
cp (car tmp)
ps (cons (cadr tmp) ps)
                   )
     
            ) ; end while lres lang
   
    (setq f (cons (car cp) f))
 
) ; end while ps

        (cons a1 f)
)

; Subdivide a pair defined by 3 points  (a1 b1 c1)                                                 
; Return a list of 6 points    (a1 d1 e1 e1 d2 c1)                                                 

(defun subdiv (a1 b1 c1 / m1 m2 m3)
     (setq m1 (mid a1 b1)
   m3 (mid b1 c1)
   m2 (mid m1 m3)
     )
     (list (list c1 m3 m2)(list m2 m1 a1))
)


; Midpoint of two points                                                                           
(defun mid (a b)
     (mapcar '(lambda (a b) (/ (+ a b) 2.0)) a b)
)

; Given 3 points, return the angle smaller than pi between a b and c b                             
(defun defl (a b c / d)
 (abs (- (setq d (- (angle b a)(angle b c))) (* (fix (/ d pi)) (+ pi pi))))
)

;Generate a bunch of points for testing on layer named points.                                     
(defun c:gen ( / n rangex rangey rangez)

(setq n (getint "\nNumber of points: ")
      rangex 5000          ; Extent in X for the points                                   
      rangey (/ rangex 1.6); Extent in Y * Golden Ratio                                   
      rangez 10            ; Extent in Z                                                   
)     
(while (> n 0)
        (entmake
                    (list
                       '(0 . "POINT")
                       '(8 . "Points")
                        (cons 10 (list (* rangex (ran)) (* rangey (ran)) (* rangez (ran))))
                    )
                )
(setq n (1- n))
)
)
; Random number generator                                                                         
(defun ran ()
    (setq seed (if seed (rem (+ (* seed 15625.7) 0.21137152) 1) 0.3171943))
)




Still to do:


To test, downlad tins.lsp attached to this post.
Appload it to Autocad.
Issue command gen and creates a bunch of point.
Issue command tin and select all points created on previous step.

Triangulation and Smoothed Contours will be created,
on layer "TIN" and layer "Contour".

ymg
Title: Re: Triangulation (re-visited)
Post by: snownut2 on July 21, 2013, 03:42:20 PM

Added Dialog Interface for Selecting Major & Minor Colors, Major & Minor Contour Intervals & Smoothing Factor
(all contained within LISP so there is only one file)

This also works in Bricscad 13......

See attached file for code.

I have also found a bug, the contour iteration is not completing all contours.  I have attached a dwg file that will not complete all contours.


Title: Re: Triangulation (re-visited)
Post by: ymg on July 22, 2013, 02:12:12 PM
Snowmut,

Thanks for the dcl for contour.

Quote
I have also found a bug, the contour iteration is not completing all contours

Bug is confirmed, need to do a test or disturb position of point when a point is exactly on Contour level.
In your example be bug at point # 3 which is exactly at elevation 100.

I'll add the necessary and repost.

Thanks for the testing.

ymg
Title: Re: Triangulation (re-visited)
Post by: snownut2 on July 22, 2013, 05:34:14 PM
Found another bug of sorts, the contours where not starting out at an elevation that was a multiple of the interval.  I have fixed that along with the contour layer placement (related) that added yesterday. 

See attached file (seems the entire code is to long to place in this box)

Revised code, added some error trapping in DCL and activated smoothing in DCL, also added some layer control to make sure appropriate layers are on.

Revised to;
       Accept either Blocks or Points
       Added TIN Name option to Dialog

Updated Dialog Interface to;
      Selection for Blocks or Point Objects to Contour (radio buttons)
      Added Error checking to ensure all required items have been done
     
I have also incorporated all the latest updates from ymg's latest file.
     

Title: Re: Triangulation (re-visited)
Post by: ymg on July 23, 2013, 03:57:42 PM
Find attached tinsd.lsp which fixes the bugs found so far.

Also implemented variable smoothing as per Christensen's

A dialogue box compliment of Snownut2 is included for setting
various parameters for contouring.  Some of this seems to be
from LeeMac.

I need to get going on adding constraints to the triangulation.

Again would appreciate suggestion for speeding up, bug report
as well as better way to write some section.

ymg

P.S. Take the latest file below
Title: Re: Triangulation (re-visited)
Post by: snownut2 on July 23, 2013, 04:06:10 PM
Indeed the LM:Popup function is included in the code, it has not been modified from Lee's code.

Bruce
Title: Re: Triangulation (re-visited)
Post by: ymg on July 23, 2013, 11:04:13 PM
Here is two images of the output from tinsd.

Only a 100 points and 20 contours.

On these images I have superimposed output with max smoothing,
min smoothing and no smoothing.

Not sure the variable smoothing is worth the trouble.

ymg

(http://contour.png)

(http://contour no tin.png)
Title: Re: Triangulation (re-visited)
Post by: snownut2 on July 24, 2013, 12:27:37 AM
I am not surprised at what you are calling minor differences, if the lines where shifted/relocated much more then the contour accuracy would suffer.  Typically as a rule in the civil world you cannot be more then 1/2 the contour interval off at any location in the contoured area.  That is on the outside you would like to see the contours closer than that, if the lines where adjusted anymore with the smoothing then you could be outside of that accuracy range.

Just as a thought, if you where to adjust the lres variable (line resolution) at the same time the smoothing factor is adjusted, making the resolution smaller as the smoothing factor is increased, the results are definitely more noticeable.

Regarding triangulation constraints, it appears that the very large angle >160deg triangles lay at the outer edges of the TIN.  These triangles cause the contour lines to take a sharp turn and head off into an un-natural direction.  How about not allowing triangles with any angle greater than 160deg. (or some other limit that makes sense)
Title: Re: Triangulation (re-visited)
Post by: ymg on July 24, 2013, 01:05:59 PM
Quote
Just as a thought, if you where to adjust the lres variable (line resolution) at the same time the smoothing factor is

That suggestion makes perfect sense.

The second one about triangle could prove to be more involved.
However if we implement constraints with an outer boundary
that problem would go away.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on July 24, 2013, 02:40:13 PM

Here is a case that tinsd does nor handle properly.

Has to do when a contour has only one edge.

Quote
Unacceptable, 3 lashes with a wet hypotenuse.

Will test and correct.  By now the wet hypotenuse
of MP is used up. So I will use the opposite side.   :pissed:

ymg

Title: Re: Triangulation (re-visited)
Post by: ymg on July 24, 2013, 04:22:53 PM
Bug is corrected.

Find below revised file.


Code - Auto/Visual Lisp: [Select]
  1. (foreach p xl
  2.         (setq ent nil)
  3.         (if (equal (car p) (reverse (last p)))
  4.                  (setq isclosed 1
  5.                               p (append p (list(cadr p)))
  6.                  )
  7.                  (setq isclosed 0
  8.                             ent (list (cons 10 (clv l (car p) pl)))
  9.                  )
  10.         )
  11.        
  12.        (if (< (length p) 3); Added to handle the case where polyline is of lenght 2.              
  13.            (setq  v3 (clv l (cadr p) pl))
  14.         )
  15.         (while (> (length p) 2)
  16.                
  17.            (setq  v1 (clv l (car p) pl)
  18.                   v2 (clv l (cadr p) pl)
  19.                   v3 (clv l (caddr p) pl)
  20.                  prv (car p)
  21.  
  22. etc. etc.
  23.  

(http://bug2.png)
Title: Re: Triangulation (re-visited)
Post by: ymg on August 30, 2013, 06:37:31 PM
I have updated the triangulation so that an adjacency list or neighbour list of triangle is generated.

Also modified so that we generate Voronoi Diagram plus the Delaunay triangulation.

The neigbour list now permits us to locate ourself efficiently in the tin via a new routine called Triloc.
This is actually Lawson's walk who invented the procedure.  Since then many different scheme have been devised
to accelerate it when you have huge triangulation.

Here's a paper showing some of them:  Walking Location Algorithms (http://graphics.zcu.cz/Download-document/165-Walking-Location-Algorithms)  by Roman Soukal.

Another paper Walking in a triangulation  (http://ftp://ftp-sop.inria.fr/geometrica/pion/publis/Walking_in_a_triangulation_socg_2001.pdf) by Olivier Devillers analyzes the speed of the  various schemes.

Here are the added function to the triangulation:

Code - Auto/Visual Lisp: [Select]
  1. ;;****************************************************************************;
  2. ;; (getneighbour pl tl)                                                                                                                   ;
  3. ;;                                                                                                                                                 ;
  4. ;; Parameters: pl, Points List                                                                                                        ;
  5. ;;             tl, Triangle List (index into pl)                                                                                       ;
  6. ;;    Returns: nl, List of neigbours for each edges of each triangles.                                               ;
  7. ;;                                                                                                                                                 ;
  8. ;; Should be done by the triangulation.                                                                                         ;
  9. ;;****************************************************************************;
  10.  
  11. (defun getneighbour (pl tl / nl pos tmp)
  12.  
  13.              (setq el nil)
  14.              (foreach tr tl
  15.                 (setq el (cons (list (car tr)(cadr tr)) el)
  16.                       el (cons (list (cadr tr)(caddr tr)) el)
  17.                       el (cons (list (caddr tr)(car tr)) el)
  18.                 )
  19.              )
  20.              (setq el (reverse el))
  21.  
  22.             (setq tmp nil)
  23.             (foreach e el
  24.                 (setq pos (vl-position (reverse e) el))
  25.                 (if pos
  26.                    (setq tmp (cons (/ pos 3) tmp)); Integer Division          ;
  27.                    (setq tmp (cons nil tmp))
  28.                 )
  29.             )
  30.            
  31.             (setq nl nil)
  32.             (while tmp
  33.                 (setq nl (cons (list (caddr tmp)(cadr tmp)(car tmp)) nl)
  34.                      tmp (cdddr tmp)
  35.                 )
  36.             )
  37.             nl
  38. )
  39.  
  40.  
  41. ;;****************************************************************************;
  42. ;; (triloc p)                                                                                                                                   ;
  43. ;;                                                                                                                                                ;
  44. ;; Locates triangle which encloses point p using Lawson's Walk.                                                    ;
  45. ;;                                                                                                                                                ;
  46. ;; Given p a point, Returns Index in tl of triangle containing the point.                                          ;
  47. ;; If outside the triangulation Return is nil.                                                                                   ;
  48. ;;                                                                                                                                                ;
  49. ;; Point List pl and Neigbour List nl are defined outside this routine.                                              ;
  50. ;; by ymg  August 2013                                                                                                               ;
  51. ;;****************************************************************************;
  52.  
  53. (defun triloc (p / v v1 v2 v3)
  54.  
  55.     ; Negative return, point is on right of v1->v2                                                                            ;
  56.     ; Positive return, point is on left of v1->v2                                                                                ;
  57.     ; 0 return, point is smack on the vector.                                                                                    ;
  58.     (defun onside (p v1 v2 /)
  59.         (setq v1 (list (car v1)(cadr v1))); We want 2d points                                                             ;
  60.         (apply '- (mapcar '* (mapcar '- v2 v1) (reverse (mapcar '-  p v1))))
  61.     )
  62.  
  63.     (setq found nil)
  64.     (if (not tn) (setq tn (/ (length tl) 2)))
  65.     (while (and (not found) tl)
  66.         (setq  v (nth tn tl)
  67.               v1 (nth (car v) pl)
  68.               v2 (nth (cadr v) pl)
  69.               v3 (nth (caddr v) pl)
  70.         )      
  71.         (cond
  72.            ((minusp (onside p v1 v2))(setq tn (car (nth tn nl))))
  73.            ((minusp (onside p v2 v3))(setq tn (cadr (nth tn nl))))
  74.            ((minusp (onside p v3 v1))(setq tn (caddr (nth tn nl))))
  75.            (t (setq found t))      ; We are inside triangle tn                ;
  76.         )
  77.         (if (not tn)(setq found t)); We are outside the triangulation         ;
  78.     )
  79.     tn ;
  80. )
  81.  
  82.  

The body of the triangulation has been modified to make use of the getneighbour function.
Note that contour has been removed from the attachment as I am currently doing a major
rewrite to make use of the neighbour list and hopefully gain some speed.

I have also put together a demoz routine to demonstrate what can be done with triloc and the neighbour list.
It is heavily inspired (leeched  :?)  by LeeMac's grtext.lsp  (http://www.lee-mac.com/grtext.html)

Here is the code to that one followed by an animated gif showing it in action.

Code - Auto/Visual Lisp: [Select]
  1. ;;****************************************************************************;
  2. ;; Demoz                                                                                                                                     ;
  3. ;; Demonstrate Lawson's Walk, to locate a triangle in a triangulation                                             ;
  4. ;; As you hover around the screen, the triangle under the screen is                                               ;
  5. ;; highlighted and the coordinates of the cursor position including Z, plus                                      ;
  6. ;; the triangle number are written to the screen.                                                                           ;
  7. ;;                                                                                                                                                 ;
  8. ;; Inspired by LeeMac's Grtext.lsp                                                                                                 ;
  9. ;;                                                                                                                                                 ;
  10. ;; by ymg        August 2013                                                                                                          ;
  11. ;;****************************************************************************;
  12.  
  13. (defun c:demoZ ( / *error* pnt str col ent etxt p prevtxt prevhatch
  14.                    prevtn scl t1 t2 t3 tr txt txtp xof yof)
  15.  
  16.     (defun *error* ( m ) (princ m) (redraw) (princ))
  17.  
  18.     (toglay '("Points" "Voronoi"))
  19.      
  20.     (while (= 5 (car (setq pnt (grread nil 13 0))))
  21.         (redraw)
  22.         (setq p (cadr pnt))
  23.         (setq tn (triloc p))
  24.         (if tn
  25.           (progn
  26.             (setq tr (nth tn tl)
  27.                   t1 (nth (car tr) pl)
  28.                   t2 (nth (cadr tr) pl)
  29.                   t3 (nth (caddr tr) pl)
  30.                    p (getz p t1 t2 t3)
  31.             )    
  32.             (if (not (= tn prevtn))
  33.                (progn
  34.                   ; Built with MakeEntmake.lsp by Cab at TheSwamp             ;
  35.                   (setq ent (list (cons 0 "HATCH")
  36.                                   (cons 100 "AcDbEntity")  
  37.                                   (cons 8 "Tin")
  38.                                   (cons 62 1)
  39.                                   (cons 440 33554534)
  40.                                   (cons 100 "AcDbHatch")
  41.                                   (cons 10 (list 0.0 0.0 0.0))
  42.                                   (cons 210 (list 0.0 0.0 1.0))      
  43.                                   (cons 2 "SOLID")
  44.                                   (cons 70 1)
  45.                                   (cons 71 0)
  46.                                   (cons 91 1)
  47.                                   (cons 92 1)
  48.                                   (cons 93 3)
  49.                                   (cons 72 1)
  50.                                   (cons 10 (list (car t1)(cadr t1) 0.0))
  51.                                   (cons 11 (list (car t2)(cadr t2) 0.0))
  52.                                   (cons 72 1)
  53.                                   (cons 10 (list (car t2)(cadr t2) 0.0))
  54.                                   (cons 11 (list (car t3)(cadr t3) 0.0))
  55.                                   (cons 72 1)
  56.                                   (cons 10 (list (car t3)(cadr t3) 0.0))
  57.                                   (cons 11 (list (car t1)(cadr t1) 0.0))
  58.                                   (cons 97 0)  
  59.                                   (cons 75 1)
  60.                                   (cons 76 1)
  61.                                   (cons 98 1)
  62.                                   (cons 10 (list 0.0 0.0 0.0))
  63.                                   (cons 450 0)
  64.                                   (cons 451 1)
  65.                                   (cons 460 0.0)
  66.                                   (cons 461 0.0)
  67.                                   (cons 452 0)
  68.                                   (cons 462 1.0)
  69.                                   (cons 453 2)
  70.                                   (cons 463 0.0)
  71.                                   (cons 63 5)
  72.                                   (cons 421 255)
  73.                                   (cons 463 1.0)
  74.                                   (cons 63 2)
  75.                                   (cons 421 16776960)
  76.                                   (cons 470 "LINEAR")    
  77.                             )
  78.                   )
  79.                   (if prevhatch (entdel prevhatch))
  80.                   (entmake ent)
  81.                   (setq prevhatch (entlast)
  82.                         prevtn tn
  83.                   )
  84.                   (redraw)
  85.                )
  86.             )
  87.           )
  88.           (if prevhatch (progn (entdel prevhatch)(setq prevhatch nil)))
  89.         )
  90.         ; Leeched from LeeMac (LM:GrText)                                     ;
  91.         ; http://www.lee-mac.com/lisp/html/GrTextDemo.html                    ;
  92.         (setq  col 2
  93.                xof 30
  94.                yof -30
  95.                scl (/ (getvar 'viewsize) (cadr (getvar 'screensize)))
  96.                  p (trans p 1 2)
  97.                str (mapcar 'rtos (trans p 1 0))
  98.               txtp (list (+ (car  p) (* xof scl)) (+ (cadr p) (* yof scl)))
  99.                txt (strcat  "T#: " (if tn (itoa tn) "Outside") "\\PX= " (car str) "\\PY= " (cadr str) "\\PZ= " (caddr str))
  100.         )
  101.         (setq etxt (list (cons 0 "MTEXT")
  102.                          (cons 100 "AcDbEntity")
  103.                          (cons 67 0)
  104.                          (cons 8 "Tin")
  105.                          (cons 62 col)
  106.                          (cons 100 "AcDbMText")
  107.                          (cons 10 txtp)
  108.                          (cons 40 30)
  109.                          (cons 41 360.0)
  110.                          (cons 46 0.0)
  111.                          (cons 72 5)
  112.                          (cons 1 txt)
  113.                          (cons 7 "Standard")
  114.                          (cons 11 '(1.0 0.0 0.0))
  115.                          (cons 42 237.0122)
  116.                          (cons 43 181.39154)
  117.                          (cons 50 0)
  118.                          (cons 73 1)
  119.                          (cons 44 1.0)
  120.                   )
  121.         )
  122.         (if prevtxt (entdel prevtxt))
  123.         (entmake etxt)
  124.         (setq prevtxt (entlast))
  125.     )
  126.     (if prevtxt (entdel prevtxt))
  127.     (if prevhatch (progn (entdel prevhatch)(setq prevhatch nil)))
  128.     (toglay '("Points" "Voronoi"))
  129.     (redraw)
  130.    
  131.    
  132.   (princ)
  133. )
  134.  
  135.  
  136.  
  137.  
  138. ;;****************************************************************************;
  139. ;; (getz p t1 t2 t3)                                                                                                                        ;
  140. ;; Given point p and triangle defined by points t1, t2, t3                                                                ;
  141. ;; Returns: (x y z) where z is on face of triangle.                                                                           ;
  142. ;;                                                                                                                                                  ;
  143. ;; By ymg  August 2013                                                                                                                ;
  144. ;;****************************************************************************;
  145.  
  146. (defun getz (p t1 t2 t3 /  v1 v2)
  147.      
  148.    ; Calculating a normal vector with gile's functions in line.               ;
  149.    ; Note that I do not calculate the unit vector for the normal n.           ;
  150.    (setq v1 (mapcar '- t2 t1)
  151.          v2 (mapcar '- t3 t1)
  152.           n (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2)))
  153.                   (- (* (caddr v1) (car v2)) (* (car v1) (caddr v2)))  
  154.                   (- (* (car v1) (cadr v2)) (* (cadr v1) (car v2)))    
  155.             )
  156.    )
  157.    (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- t1 p)))(caddr n)))
  158.  )
  159.  
  160. ;; (toglay lst)                                                                                               ;
  161. ;; Toggles (on/off) a list of Layer names.                                                      ;
  162. ;; From a Fishing Story by Cab at TheSwamp >-=((((°>                              ;
  163. ;;                                                                                                                ;
  164. ;; by ymg               August 2013                                                                 ;
  165.  
  166. (defun toglay (lst / entlist l)
  167.     (foreach l lst
  168.          (and (setq entlist (tblobjname "LAYER" l))
  169.               (setq entlist  (entget entlist))
  170.               (entmod (subst (cons 62 (- (cdr (assoc 62 entlist))))
  171.                              (assoc 62 entlist)
  172.                              entlist
  173.                       )
  174.               )
  175.          )
  176.     )
  177. )
  178.  
  179.  
  180.  
  181.  

To run the demo, load the attachment demoz.lsp, Issue command GEN and enter number of points desired
in the triangulation.

Isue command VOR, the Voronoi Diagram and the Delaunay Triangulation will generate.

Issue command DEMOZ, as you hover in the drawing, the triangle number plus coordinates X,Y and Z of point under cursor
is displayed and the currently occupied triangle is hatched in solid Red.

Sorry!, about the long post.

ymg

(http://demoz.gif)

 
Title: Re: Triangulation (re-visited)
Post by: Topographer on November 26, 2013, 05:08:41 PM
hi Ymg ,I have a question.I am typing Tin and the lisp run and makes tins
when i write contours the lisp not working.Which command should write to make contours ?
Title: Re: Triangulation (re-visited)
Post by: ymg on November 26, 2013, 10:54:17 PM
pedroantonio,

The version you have has a bug whereby the majcont variable is not initialized.

Make sure you change the setting for major contour in the dialog boxe
before running then it should work.

There will be a revised version very soon with many improvements.


ymg

You can also add in the code the following line majcnt 5.0:

Code - Auto/Visual Lisp: [Select]
  1. (setq dcl_id(load_dialog dclfiletemp)
  2.         intvL   '("  1" "  2" "  5" " 10")
  3.         majcntL '("  5" " 10" " 50" "100")
  4.         hfacL   '("  Max" "  -20%" "  -40%" "  -60%" "  -80%" " Min" "None")
  5.         hfacTL  '(0.5 0.6 0.7 0.8 0.9 1.0 5)   
  6.         majcolor 10
  7.         mincolor 14
  8.         intv 1
  9.         majcnt 5.0 ;add this line
  10.         hfac 0.5
  11.   )
  12.  
Title: Re: Triangulation (re-visited)
Post by: Topographer on November 27, 2013, 04:21:13 AM
Thank you ymg . Where i can find the new version ?
Title: Re: Triangulation (re-visited)
Post by: ymg on November 27, 2013, 05:52:19 AM
pedroantonio,

Quote
Thank you ymg . Where i can find the new version ?

The new version will be available here when I publish it.

Meanwhile just add the majcont line.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on November 29, 2013, 02:27:33 AM
Ymg i try to do this change in the code

Code: [Select]
intvL   '(" 0.02" " 0.1" " 0.2" "  1" "  2" "  5" " 10")
majcntL '(" 0.1" " 0.5" " 1" "  5" " 10" " 50" "100")
hfacL   '("  Max" "  -20%" "  -40%" "  -60%" "  -80%" " Min" "None")

but the code didn't work for this value !!

Code: [Select]
intvL   '(" 0.02" " 0.1" " 0.2"......................)
majcntL '(" 0.1" " 0.5" " 1" ..........................")
hfacL   '("  Max" "  -20%" "  -40%" "  -60%" "  -80%" " Min" "None")"

Sometimes is important if you  have the ability to give values ​​ smaller than 1m for (MNR)

can you tell me why ?
Title: Re: Triangulation (re-visited)
Post by: ymg on November 29, 2013, 04:05:39 AM
pedroantonio,

The next version does fractionnal contour.

You need to change more than just the list for this to happen.

If you have an immediate need I can send you an advanced copy.

However there are still quite a few things to change before issuing it at large.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on November 29, 2013, 05:44:36 AM
yes if you can ...

Thanks
Title: Re: Triangulation (re-visited)
Post by: ymg on February 09, 2014, 04:39:14 AM
Here is an update to the TIN program.

It is still very much a work in progress, so I would appreciate
bug reports and enhancement  suggestion.

The Contouring part has been accelerated a little, and now work with
fractionnal contours.

You can also call CONT to do contours on a Selection of existing 3DFACE.

Express Tools is required, as I make use of the Progress Bar.

Still missing Breaklines handling for the triangulation,
and Highlighting of Depression Contours for the contouring section.


ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on February 22, 2014, 11:02:07 PM
Here is an update to the TIN program.

It is still very much a work in progress, so I would appreciate
bug reports and enhancement  suggestion.

The Contouring part has been accelerated a little, and now work with
fractionnal contours.

You can also call CONT to do contours on a Selection of existing 3DFACE.

Express Tools is required, as I make use of the Progress Bar.

Still missing Breaklines handling for the triangulation,
and Highlighting of Depression Contours for the contouring section.
ymg
Great Program , Yang !
The last version is higher effect .
I has some suggest :
1. Just your description , "You can also call CONT to do contours on a Selection of existing 3DFACE." , I hope It can .
2. Smoothing effect may also need to improve some , just like the pic. I post .
Title: Re: Triangulation (re-visited)
Post by: ymg on February 22, 2014, 11:28:17 PM
chlh_jd,

Quote
"You can also call CONT to do contours on a Selection of existing 3DFACE." , I hope It can

Thanks for the feedback.

Agree that my wording is not too clear.  What was meant by that is in the case where you already have
a selection of 3dfaces created by another program, you can use command "CONT" and do the contouring.

Normal way would be a set of points that you triangulate and then contours.

I do not understand exactly what you mean from the attached image.  Could you clarify?

I've been playing around a bit with it and manage to accelerate it some. Shall post it soon.
We still need to handle the marking of "Depression Contour" currently not available.

I highly recommend Christensen's paper, it is well worth the reading.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on February 23, 2014, 09:17:39 AM
chlh_jd,

Quote
"You can also call CONT to do contours on a Selection of existing 3DFACE." , I hope It can

Thanks for the feedback.

Agree that my wording is not too clear.  What was meant by that is in the case where you already have
a selection of 3dfaces created by another program, you can use command "CONT" and do the contouring.
I'm sorry for Just finished reading your code in *gress5.lsp and foud add cont command .
Quote
Normal way would be a set of points that you triangulate and then contours.

I do not understand exactly what you mean from the attached image.  Could you clarify?
Sorry for my poor English ,  I mean that , the zone I mark , smoothed contours less than ideal .
I've seen Christensen's paper , perhaps no suit for the case I tesed .

Quote
I've been playing around a bit with it and manage to accelerate it some. Shall post it soon.
We still need to handle the marking of "Depression Contour" currently not available.

I highly recommend Christensen's paper, it is well worth the reading.

ymg
Not long ago, I also wrote a contour program is based on existing 3DFACES, and not for smoothing; and I'v lost in smooth annoying.
Looking forward to your masterpiece .  :-)

Title: Re: Triangulation (re-visited)
Post by: ymg on February 23, 2014, 10:38:13 AM
I believe that this zone would get crossing contours unless you do it this way.

In other word this is inherent to the method.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on March 05, 2014, 12:20:36 AM
I believe that this zone would get crossing contours unless you do it this way.

In other word this is inherent to the method.

ymg

Maybe this zone must add point .

I agree you just I can't found better .  :-)

In THis true case , it take some wrong .
Title: Re: Triangulation (re-visited)
Post by: ymg on March 05, 2014, 05:56:15 AM
chlh_jd,

Yes this is a bug introduced in the latest round of optimization.

Trying to locate the culprit and will post a revision.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on March 06, 2014, 12:28:37 AM
chlh_jd,

Here is a revision with the bugs in contour generation removed.

Bear in mind, that this is still a work in progress, so the routine does  not clean themselves nicely on completion.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 07, 2014, 10:38:39 AM
Hello ymg nice job but , I find a bug in your code and i dont know why.

I write TIN ---> make TINS
I write CONT ----> open the contour window and gives me this error

Code - Auto/Visual Lisp: [Select]
  1. Select objects: Specify opposite corner: 71 found
  2. Select objects:
  3. Error: bad argument type: fixnump: nil

here is my dwg file ... Any options ?

Thanks
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 07, 2014, 04:40:42 PM
with the  ContourTest.dwg  the lisp working fine but with my dwg is not working .Can you tell me what i am doing wrong ?
Thanks
Title: Re: Triangulation (re-visited)
Post by: CAB on March 07, 2014, 05:21:28 PM
I think he may be on the other side of the world so you may need to give some time for a response.  8)
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on March 08, 2014, 01:43:20 AM
chlh_jd,

Here is a revision with the bugs in contour generation removed.

Bear in mind, that this is still a work in progress, so the routine does  not clean themselves nicely on completion.

ymg
Thanks Ymg , it did fix the bug in the case "intv=0.25 majcnt=1.0"; but error in the case  "intv=0.1 majcnt=0.5".

And some suggestion in 'contour function:
1. defined lambda function will short codes
Code - Auto/Visual Lisp: [Select]
  1. (defun f1  ()
  2.     (setq  z1  (caddr p1)
  3.            z2  (caddr p2)
  4.            z1 (if (= z1 z) (setq z1 (- z1 1e-08)) z1)
  5.            z2 (if (= z2 z) (setq z2 (- z2 1e-08)) z2)
  6.            p1  (list (car p1) (cadr p1)))
  7.       (polar p1
  8.              (angle p1 p2)
  9.              (* (distance p1 p2) (/ (- z z1) (- z2 z1)))))
  10.  
to
Code - Auto/Visual Lisp: [Select]
  1. (if (equal (car p) (reverse (last p)))
  2.         (setq isclosed 1
  3.               p        (append p (list (cadr p))))
  4.         (setq isclosed 0
  5.               p1       (nth (caar p) pl)
  6.               p2       (nth (cadar p) pl)            
  7.               ent      (list (cons 10
  8.                                    (f1)));_
  9.               )
  10.         )
  11. ...
  12. (while (> (length p) 2)
  13.         (setq p1  (nth (caar p) pl)
  14.               p2  (nth (cadar p) pl)         
  15.               v1  (f1)
  16.               p1  (nth (caadr p) pl)
  17.               p2  (nth (cadadr p) pl)        
  18.               v2  (f1)
  19.               p1  (nth (caaddr p) pl)
  20.               p2  (nth (cadr (caddr p)) pl)          
  21.               v3  (f1)
  22.               prv (car p)
  23.               nxt (caddr p)
  24.               p   (cdr p)
  25.               p1  (nth (caar p) pl)
  26.               p3  (nth (cadar p) pl)
  27.               p2  (nth (car prv) pl)
  28.               p4  (nth (car nxt) pl)
  29.               )
  30.   ...
  31.  
2. in the 'cl while loop , (setq n ...) before (while n ...) will reduce one time vl-position .
Code - Auto/Visual Lisp: [Select]
  1.       (setq n (vl-position nxt cl));_reduce 'vl-position one time , edited by GSLS(SS)
  2.       (while n
  3.         (setq cl  (vl-remove nxt cl)
  4.               n   (- n (rem n 2))
  5.               m   (nth n cl)
  6.               pol (cons m pol)
  7.               cl  (vl-remove m cl));_reduce 'nth one time
  8.         (if (vl-position nxt pol)
  9.           (setq nxt nil)
  10.           (setq nxt (reverse (car pol))))
  11.         (and (not (setq n (vl-position nxt cl)))
  12.              (setq pol (reverse pol)
  13.                    nxt (reverse (car pol))
  14.                    n   (vl-position nxt cl)
  15.                    ))
  16.         )
  17.  
3. define lambda function to replace vl-remove for build cl list loop will improve speed for huge list ; It also can be added position determine for > (/ (length lst) 2) , and use reverse to speed up .
Code - Auto/Visual Lisp: [Select]
  1.   ;; remove-nth for contour
  2.   (defun f2  (i lst / fst)
  3.     (setq fst nil)
  4.     (repeat (rem i 4)
  5.       (setq fst (cons (car lst) fst)
  6.             lst (cdr lst)))
  7.     (repeat (/ i 4)
  8.       (setq fst (vl-list* (cadddr lst)
  9.                           (caddr lst)
  10.                           (cadr lst)
  11.                           (car lst)
  12.                           fst)
  13.             lst (cddddr lst)))
  14.     (append
  15.       (reverse fst)
  16.       (cdddr lst)))
  17. ;;--------------
  18. ...
  19. (while (< j (length el))
  20.       (setq e  (nth j el)
  21.             z1 (nth (car e) zl)
  22.             z2 (nth (cadr e) zl)
  23.             zm (max zm z1 z2)
  24.             j  (1+ j))
  25.  ;_ Reduce size of el on zmax criteria.                                  ;
  26.       (if (and (= (rem j 3) 0) (< zm (+ z intv)))
  27.         (setq j  (- j 3)
  28.               el (f2 j el);_remove-nth for contour
  29.               zm -1e19))
  30.       (if (= z1 z)
  31.         (setq z1 (- z1 1e-8)))          ; If on Interval we disturb         ;
  32.       (if (= z2 z)
  33.         (setq z2 (- z2 1e-8)))          ; the z value a little.             ;
  34.       (if (or (< z1 z z2) (> z1 z z2))
  35.         (setq cl (cons e cl))           ; Edge is added to Crossed List     ;
  36.         )      
  37.       )
  38. ...
  39.  
4. don't shift 'clayer between major with minor ,just use (cons 8 ...)
Code - Auto/Visual Lisp: [Select]
  1. ...
  2. (setq seg       (length ent)
  3.             vcs (+ vcs seg)
  4.             ent (vl-list* (cons 0 "LWPOLYLINE")
  5.                           (cons 100 "AcDbEntity")
  6.                           (cons 100 "AcDbPolyline")
  7.                           (cons 8
  8.                                 (if (zerop (rem z majcnt))
  9.                                   "Contour Major"
  10.                                   "Contour Minor"))                      
  11.                           (cons 90 seg)                  
  12.                           (cons 70 isclosed)
  13.                           (cons 38 z)
  14.                           ent)
  15.             )
  16. ...
  17.  
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 08, 2014, 05:37:28 AM
I try the prof command but i can't change the scale .

And what ocd command do?

Can you give me some instructions for them ?

Thanks
Title: Re: Triangulation (re-visited)
Post by: ymg on March 09, 2014, 08:09:55 AM
Topographer,

Sorry for late reply.

In your drawing there is an incorrect triangle, so the sorting routine for the triangle list
crash on it.

I will check the get_tin routine to see if there is a bug there, or add something to prevent this.

Meanwhile use tin on the point and you will get the contour.

The prof routine is for making a profile from a set of 3dface, or a triangulation or yet again from
a 3dpolyline.  However it is not finished, so not ready for production.

OCD is a routine i got from AUGI forum, will be use to trim triangulation.  This is also not ready

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on March 09, 2014, 08:20:03 AM
chlh_jd,

The previous version of the program had function like that.

But in an effort to squeeze every bit of speed I`ve put everything inline
thus gaining a tiny bit of speed by not calling a function.  On a 1000 points
triangulation this is about 50000 calls.

Second suggestion make sense to me.

Setting the current layer also saves me as many  calls to cons  as there are vertices in the contours.

Thanks for the helping hands.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 09, 2014, 12:29:15 PM
Thank you ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on March 10, 2014, 04:26:04 AM
Topographer,

The problem was in the GET_TIN routine.

Find a revision attached below.

By the way OCD stands for Outside Contour Delete.
So given a boundary polyline, it will trim everything that is outside
polyline.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 10, 2014, 05:04:05 AM
Thank you  ymg

I have two comments

1) OCD delete everething out of the polyline. This is a problem because  delete all my measure points, block,TINS , etc. Can  you fix it to delete only contours ?
It is good to have all the TINS .I forget  delete things in layer off
2) Can you add breaklines is important for correct TINS.

Nice job keep working we need you !!!!!

Thanks
Title: Re: Triangulation (re-visited)
Post by: ymg on March 10, 2014, 01:47:49 PM
Topographer,

As stated earlier OCD did not originate with me.

It is in the program because we will need a  functionnality of that type
sooner or later.  Right now I did not even look at it although I know we
will need to modify it some.

Breaklines are certainly a needed additions.  I have done some exploratory
work on it.  Still not ready for prime time though.

Understand that I do this as a hobby, so I sometimes get fed up of looking
at it and move on to something else.  I let it cook then come back to it.

There is also some more work to be done with the contouring.  Namely we
need to highlight the depression contour.  There too I've done a bit of work
but not ready.  Another recurring problem is roundoff error with the fractionnal
contour.  So far still haven't got 100% reliability.

ymg
Title: Re: Triangulation (re-visited)
Post by: Topographer on March 10, 2014, 01:56:44 PM
I understand , Thank you.
Title: Re: Triangulation (re-visited)
Post by: ymg on March 11, 2014, 10:59:09 AM
chlh_jd,

Quote
Thanks Ymg , it did fix the bug in the case "intv=0.25 majcnt=1.0"; but error in the case  "intv=0.1 majcnt=0.5".

I had to add a fuzz factor when we compare current contour z to see if the points z1 and z2 are equal to it.

I have also increased the disruption to 1e-03 from 1e-08.

Seems to have cured it!, so find attached the revision.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on March 12, 2014, 09:09:17 AM
Hello Ymg ,
I'm sorry , did not read your every version .
Thank you for your reply and fixing the bug .
I'll take much time to test the new version after some time ,
Very pleased to test it and help some . :-)
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on March 12, 2014, 09:31:21 AM
Hi Ymg,
In my head , Longer code still requires a larger computer memory , so shorting the codes is necessary for reduce computer memory usage , if we have enough computer memory  , we also can build other ways list for contour routine  to replace a large number of 'nth' . So we have to compare the speed after complied the file .
In the contour routine , make a zl will speed up some .
Code: [Select]
...
(setq ti   (car (_VL-TIMES)) ; Re-initialize timer for Contouring
zl   (mapcar (function caddr) pl) ;_only for zlist                               
zmin (apply (function min) zl) ;_improve speed by GSLS(SS) 2014.3.2             
zmax (apply (function max) zl)
zmin (* (fix (/ zmin intv)) intv) ;_zmin, first level for contour.   (* (+ (fix (/ zmin intv)) 1) intv)
zmax (* (fix (/ zmax intv)) intv) ;_ zmax; last level for contour.               
z    zmin ; Initialize Current Contour level.             
...                 
)
  (setq tl (vl-sort tl
    (function
      (lambda (a b)
(< (max (nth (car a) zl)
(nth (cadr a) zl)
(nth (caddr a) zl))
   (max (nth (car b) zl)
(nth (cadr b) zl)
(nth (caddr b) zl)))))))
...
Title: Re: Triangulation (re-visited)
Post by: ymg on March 13, 2014, 01:19:00 AM
chlh_jd,

Thanks, for the sort idea.

On a 2000 triangles list I do gain about 0.25 seconds (Non-Compiled)

Quote
In my head , Longer code still requires a larger computer memory

Very right, I am trading space for speed.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on March 14, 2014, 04:49:59 AM
Yet another bug in the contour generation routine removed.

Also implemented some of chlh_jd's suggestions.

Hopefully everything's OK now, so I can start work on
Depression Contours.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on March 15, 2014, 02:02:36 PM
Ymg,
The new version test OK in contour routine . Thank you very much .
I'm really looking forward to deal with discontinuous triangulation which just you said no support .
Title: Re: Triangulation (re-visited)
Post by: motee-z on March 15, 2014, 06:21:05 PM
excellent job for countor lines but it will be perfect if you take into consider break lines.
(break lines is a 3d polyline so 3d face not intersect with it because it represent a cliff in land)
thanks
Title: Re: Triangulation (re-visited)
Post by: ymg on March 16, 2014, 12:55:41 AM
motee-z,

I know what is a breakline.

I have done some preliminary work on it.

But we also have to take into considerations boundaries of the tin
plus holes into the triangulation.

Conceptually all we need to do is to erase all the triangles intercepted
by the breaklines, then re-triangulate the hole created.

My thinking would be a CCW polyline would represent a boundary where
the tin is inside. Clockwise poly would be a hole.

But there is a lot of catch in there.

ymg

Title: Re: Triangulation (re-visited)
Post by: teknomatika on March 31, 2014, 10:32:12 AM
ymg,
is a fantastic tool.
However, PROF function is not working. Always returns the same error.

Quote
Command: prof

Select a Linear Entity:

Error: bad argument type: fixnump: nil

I am using a lwpolyline as a linear entity. Not sure what am I doing wrong.
Attached the sample file
Title: Re: Triangulation (re-visited)
Post by: Pad on March 31, 2014, 11:20:46 AM
ymg,

brilliant work, amazing.

P
Title: Re: Triangulation (re-visited)
Post by: ymg on March 31, 2014, 11:25:31 AM
teknomatica,

Prof is not completed as explained in a previous post.

I am currently working on it.

Try it with the attached, you should be able to get the profile
from a tin at least.

There is still quite a bit of work to be done on the UI.

Also, I am not satisfied with the current selection method
the getneighbour routine is way too slow.

So far, I have come up with a new method.  Now I need
to complete the UI for every case. (Tin, Contour, 3DFACE and 3DPOLY)

Also working on the surface selection interface.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on March 31, 2014, 11:32:15 AM
Pad,

Thank you , for the encouragement.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on April 10, 2014, 05:56:29 AM
hi, ymg
The contour routine still get error  in the case of the dwg I post ,when intv = 0.1 majintv=0.5 .
Even the Tin routine don't get right result in this drawing .
Title: Re: Triangulation (re-visited)
Post by: ymg on April 10, 2014, 01:25:30 PM
chlh_jd,

In your drawing you had at least 20 illegal 3dfaces shaped as a line.

Don 't know how you generated them.  Once these are removed everything works.

You also claim that the Tin did not work.  I could not check as you did not provide
the points.

ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on April 10, 2014, 05:46:48 PM
I did regenerate your points.

And you are bound to have troubles when the points are nearly on top
of each other.

Same problem with your circle, reduce the density of the points.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on April 11, 2014, 07:50:47 AM
Ymg,
I'm sorry for not provide the points .  :oops:
Thank you for check this case .
In realistic terrain , there are some steep and cliffs , I wish the contour routine can support this case .
You have also mentioned, can not handle discontinuous 3DFACE , I wish it can .

Therefore, can we change an idea to deal with the  contour problem; For example, the first generation of each line segment at every elevations, and then combine them, and then fillet them .
Title: Re: Triangulation (re-visited)
Post by: ymg on April 11, 2014, 09:27:15 AM
chlh_jd,

I thik I found the culprit.  There were some duplicates (10 of them to be exact)
points that did not get filtered out.

modify this line at beginning of program:

Code: [Select]
(setq pl (vl-sort pl (function (lambda (a b) (< (cadr a) (cadr b))))) ;; Sort on Y Coordinate     
                pl (vl-sort pl (function (lambda (a b) (< (car  a) (car  b))))) ;; Sort on X Coordinate     
                pl (vl-remove nil                                               ;; Remove Duplicate Points 
                      (mapcar '(lambda (a b) (if (not (equal (cdr (reverse a)) (cdr (reverse b)) 1e-03)) a))
                            pl (append (cdr pl) (list (car pl)))                             
                      )
                   )
  )

Adding a fuzz factor on the comparison in the vl-remove line.

Seems to cure it.


ymg
Title: Re: Triangulation (re-visited)
Post by: ymg on April 11, 2014, 09:51:01 AM
chlh_jd,

contour ca run on disjointed set of 3dfaces.

In the picture below both sets were contoured
in one go.
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on April 13, 2014, 03:27:44 AM
Ymg,
Thank you for fix this case .
some suggest :
1. Another use case may can be : First select norm ,just like (210 1.0 0.0 0.0 ) , not only for (210 0. 0. 1.) ; and then trans all 3dfaces points , and then make contours .
   
Code: [Select]
  (contour  pl  tl intv  majcnt majcolor    mincolor hfac  norm)2. In the "get_ss" function , will be error , if the 3dface is manual change ( because of the osnap accuracy ), it can be solve by two way :
   2.1 tolerance index come in , such as (equalmember p pl 1e-6) replace  (vl-position p pl).
   2.2 fix or round coordinate before construct the pointlist , such as
     
Code: [Select]
(setq pl (cons (mapcar (function (lambda (a) (/ (fix (* a 1e6)) 1e6))) p) pl)3. In the contour function , some position may get error :
  3.1 build xl part
   
Code - Auto/Visual Lisp: [Select]
  1. (while (setq n (vl-position nxt cl))     
  2.            (setq  cl (vl-remove nxt cl)
  3.                    n (- n (rem n 2))
  4.                    m (nth n cl);_cl may be nil after first 'vl-remove .
  5.                  pol (cons m pol)
  6.                   cl (vl-remove m cl)
  7.            )
  8.            (if (vl-position nxt pol)
  9.               (setq nxt nil)
  10.               (setq nxt (reverse (car pol)))
  11.            )
  12.          
  13.            (if (not (vl-position nxt cl))
  14.               (setq pol (reverse pol)
  15.                     nxt (reverse (car pol))
  16.               )
  17.            )  
  18.        )
  19.    
   3.2 foreach for xl part
 
Code - Auto/Visual Lisp: [Select]
  1.     (foreach p xl ;_p may be has nil item , just like ((17 9) nil)
  2.         (setq ent nil)
  3.       ...)
  4.  
  3.2 in the cond t condtion
 
Code - Auto/Visual Lisp: [Select]
  1.     (foreach p xl
  2.             ...
  3.        (while (> (length p) 2)
  4.             ...
  5.           (cond (...)
  6.                    (...)
  7.                    (t (setq cn1 (list (/ (+ (car p1) (car p2) (car p3)) 3.)
  8.                            (/ (+ (cadr p1) (cadr p2) (cadr p3)) 3.))
  9.                  cn2 (list (/ (+ (car p1) (car p3) (car p4)) 3.)
  10.                            (/ (+ (cadr p1) (cadr p3) (cadr p4)) 3.))
  11.                  ;_ (cn1 p1 v1 v2) or (cn1 p3 v1 v2) may be collinear 4 points , how deal it ?
  12.                  a1  (cond
  13.                        ((inters cn1 p1 v1 v2))
  14.                        ((inters cn1 p3 v1 v2))
  15.                        )
  16.                  a3  (cond
  17.                        ((inters cn2 p1 v2 v3))
  18.                        ((inters cn2 p3 v2 v3))
  19.                        ))
  20.                           ...
  21.                        )
  22.                    ...
  23.                   )
  24.                 ...
  25.                )
  26.  
Title: Re: Triangulation (re-visited)
Post by: ymg on April 13, 2014, 01:40:57 PM
chlh_jd,

not sure I am following you on all this.

Originally the program would not check for duplicate points.

Maybe I could add something in "gettin" to detect illegal
3dfaces.  However all these checks do take time.

All in all there will always be a way where you can throw
the program off if only due to limited precision of floating
point operations.

If you feed it garbage you will get garbage.

ymg
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on April 14, 2014, 04:09:38 AM
You're right .
In the wrong direction to go right path indeed make road twists and turns .
Title: Re: Triangulation (re-visited)
Post by: chlh_jd on April 17, 2014, 10:17:15 AM
Can we use the autocad contour method ?
The output functions from acgex17.dll in ACAD2008 seems seems to have it .
Code: [Select]
?getContour@AcGeCurveBoundary@@QBEXAAHPAPAPAVAcGeEntity3d@@PAPAPAVAcGeCurve2d@@PAPAH3@Z
?getContour@AcGeImpCurveBoundary@@QBEXAAHPAPAPAVAcGeImpEntity3d@@PAPAPAVAcGeImpCurve3d@@PAPAH3@Z
?getContours@AcGeExternalBoundedSurface@@QBEXAAHAAPAVAcGeCurveBoundary@@@Z
?getContours@AcGeImpExternalBoundedSurface@@QBEXAAHAAPAVAcGeImpCurveBoundary@@@Z
Title: Re: Triangulation (re-visited)
Post by: ymg on April 17, 2014, 11:39:02 AM
chlh_jd,

Really don't know.  Worth a try!
Title: Re: Triangulation (re-visited)
Post by: hanhphuc on April 24, 2014, 03:41:25 AM
ymg,
is a fantastic tool.
However, PROF function is not working. Always returns the same error.

Quote
Command: prof

Select a Linear Entity:

Error: bad argument type: fixnump: nil

I am using a lwpolyline as a linear entity. Not sure what am I doing wrong.
Attached the sample file

hi teknomatika, maybe one of the endpoints (line entity's) is outside the TIN boundary?
Title: Re: Triangulation (re-visited)
Post by: hanhphuc on April 24, 2014, 05:03:06 AM

hi chlh_jd, this arx sample used by surveyor jobs can be applied in function: getz?
if we hide the first line (geomcal)** then we still can use ymg's getz, because ARX not loaded :-)
ymg

Code: [Select]
(if C:cal (princ "\nGeomcal loaded..")(arxload "geomcal"));** <-- this is optional

(if C:cal ;if loaded
 
(defun getz (p t1 t2 t3 / ptt ) ;defun opted by hanhphuc
(setq
ptt (list (car p)(cadr p)(1+(caddr p))))
(cal "ilp(p,ptt,t1,t2,t3)")
); getz
 
(defun getz .... code unchanged ) ;defun by ymg
       
) ;c:cal
Title: Re: Triangulation (re-visited)
Post by: ribarm on April 24, 2014, 06:08:41 AM
No need for geomcal.arx

Code - Auto/Visual Lisp: [Select]
  1. ;|(defun _ilp ( p1 p2 t1 t2 t3 / unit v^v Coplanar-p ptinsidetriangle-p n e p1e p2e p1ed p2ed p1p p2p v1 v2 )
  2.  
  3.   (defun unit ( v )
  4.     (mapcar '(lambda ( x ) (/ x (distance '(0.0 0.0 0.0) v))) v)
  5.   )
  6.  
  7.   (defun v^v ( v1 v2 )
  8.     (
  9.       (lambda ( a b )
  10.         (mapcar '(lambda ( a1 a2 b1 b2 ) (- (* a1 a2) (* b1 b2)))
  11.                  a (cdr b) b (cdr a)
  12.         )
  13.       )
  14.       (list (cadr v1) (caddr v1) (car v1) (cadr v1))
  15.       (list (cadr v2) (caddr v2) (car v2) (cadr v2))
  16.     )
  17.   )
  18.  
  19.   (defun Coplanar-p ( p1 p2 p3 p4 )
  20.     (
  21.       (lambda ( n1 n2 )
  22.         (equal (v^v n1 n2) '(0.0 0.0 0.0) 1e-8)
  23.       )
  24.       (v^v (mapcar '- p1 p2) (mapcar '- p1 p3))
  25.       (v^v (mapcar '- p1 p2) (mapcar '- p1 p4))
  26.     )
  27.   )
  28.  
  29.   (defun ptinsidetriangle-p ( pt p1 p2 p3 )
  30.     (if
  31.       (and
  32.         (Coplanar-p pt p1 p2 p3)
  33.         (not
  34.           (or
  35.             (inters pt p1 p2 p3)
  36.             (inters pt p2 p1 p3)
  37.             (inters pt p3 p1 p2)
  38.           )
  39.         )
  40.         (not
  41.           (or
  42.             (> (+ (distance pt p1) (distance pt p2)) (+ (distance p3 p1) (distance p3 p2)))
  43.             (> (+ (distance pt p2) (distance pt p3)) (+ (distance p1 p2) (distance p1 p3)))
  44.             (> (+ (distance pt p3) (distance pt p1)) (+ (distance p2 p3) (distance p2 p1)))
  45.           )
  46.         )
  47.       )
  48.       T
  49.       nil
  50.     )
  51.   )
  52.  
  53.   (setq n (unit (v^v (mapcar '- t2 t1) (mapcar '- t3 t1))))
  54.   (setq e (last (trans t1 0 n)))
  55.   (setq p1e (last (trans p1 0 n)))
  56.   (setq p2e (last (trans p2 0 n)))
  57.   (setq p1ed (abs (- p1e e)))
  58.   (setq p2ed (abs (- p2e e)))
  59.   (setq p1p (mapcar '+ p1 (setq v1 (mapcar '(lambda ( x ) (* p1ed x)) n))))
  60.   (setq p2p (mapcar '+ p2 (setq v2 (mapcar '(lambda ( x ) (* p2ed x)) n))))
  61.   (if (not (equal e (last (trans p1p 0 n)) 1e-8)) (setq p1p (mapcar '- p1 v1)))
  62.   (if (not (equal e (last (trans p2p 0 n)) 1e-8)) (setq p2p (mapcar '- p2 v2)))
  63.   (setq p (inters p1 p2 p1p p2p nil))
  64.   (if (ptinsidetriangle-p p t1 t2 t3)
  65.     p
  66.     nil
  67.   )
  68.  
  69. )|;
  70.  

Maybe this shorter variant - I think it's the same, but check - I would use this one :

Code - Auto/Visual Lisp: [Select]
  1. ;|
  2. (defun _ilp ( p1 p2 t1 t2 t3 / v^v Coplanar-p ptinsidetriangle-p n e p1n p2n p1p p2p )
  3.  
  4.   (defun v^v ( v1 v2 )
  5.     (
  6.       (lambda ( a b )
  7.         (mapcar '(lambda ( a1 a2 b1 b2 ) (- (* a1 a2) (* b1 b2)))
  8.                  a (cdr b) b (cdr a)
  9.         )
  10.       )
  11.       (list (cadr v1) (caddr v1) (car v1) (cadr v1))
  12.       (list (cadr v2) (caddr v2) (car v2) (cadr v2))
  13.     )
  14.   )
  15.  
  16.   (defun Coplanar-p ( p1 p2 p3 p4 )
  17.     (
  18.       (lambda ( n1 n2 )
  19.         (equal (v^v n1 n2) '(0.0 0.0 0.0) 1e-8)
  20.       )
  21.       (v^v (mapcar '- p1 p2) (mapcar '- p1 p3))
  22.       (v^v (mapcar '- p1 p2) (mapcar '- p1 p4))
  23.     )
  24.   )
  25.  
  26.   (defun ptinsidetriangle-p ( pt p1 p2 p3 )
  27.     (if
  28.       (and
  29.         (Coplanar-p pt p1 p2 p3)
  30.         (not
  31.           (or
  32.             (inters pt p1 p2 p3)
  33.             (inters pt p2 p1 p3)
  34.             (inters pt p3 p1 p2)
  35.           )
  36.         )
  37.         (not
  38.           (or
  39.             (> (+ (distance pt p1) (distance pt p2)) (+ (distance p3 p1) (distance p3 p2)))
  40.             (> (+ (distance pt p2) (distance pt p3)) (+ (distance p1 p2) (distance p1 p3)))
  41.             (> (+ (distance pt p3) (distance pt p1)) (+ (distance p2 p3) (distance p2 p1)))
  42.           )
  43.         )
  44.       )
  45.       T
  46.       nil
  47.     )
  48.   )
  49.  
  50.   (setq n (v^v (mapcar '- t2 t1) (mapcar '- t3 t1)))
  51.   (setq e (last (trans t1 0 n)))
  52.   (setq p1n (list (car (trans p1 0 n)) (cadr (trans p1 0 n)) e))
  53.   (setq p2n (list (car (trans p2 0 n)) (cadr (trans p2 0 n)) e))
  54.   (setq p1p (trans p1n n 0))
  55.   (setq p2p (trans p2n n 0))
  56.   (setq p (inters p1 p2 p1p p2p nil))
  57.   (if (ptinsidetriangle-p p t1 t2 t3)
  58.     p
  59.     nil
  60.   )
  61.  
  62. )|;
  63.  

Both codes aren't correct... My update can be found here :
www.cadtutor.net/forum/showthread.php?89154-Solids-intersection-and-something-else...&s=293246d38c1367703f834a7d79f7944d&p=610836#post610836

M.R.
Title: Re: Triangulation (re-visited)
Post by: ymg on April 24, 2014, 09:21:26 AM
I have gone through 5 different version of getz before settling on the last one.

Here they are:

Code - Auto/Visual Lisp: [Select]
  1. (defun getz0 (p t1 t2 t3 /  v1 v2)
  2.  
  3.    (setq v1 (mapcar '- t2 t1)
  4.          v2 (mapcar '- t3 t1)
  5.           n (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2)))
  6.                   (- (* (caddr v1) (car v2)) (* (car v1) (caddr v2)))  
  7.                   (- (* (car v1) (cadr v2)) (* (cadr v1) (car v2)))    
  8.             )
  9.    )
  10.    (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- t1 p)))(caddr n)))
  11.  )
  12.