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

0 Members and 1 Guest are viewing this topic.

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #90 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!

chlh_jd

  • Guest
Re: Triangulation (re-visited)
« Reply #91 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 .



ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #92 on: June 15, 2011, 02:14:33 PM »
The first is GA for TSP .

I'm sorry, I do not remember what this program. 

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #93 on: June 15, 2011, 02:16:26 PM »
TSP = Travelling salesman problem
and
GA  = genetic algorithm
??

chlh_jd

  • Guest
Re: Triangulation (re-visited)
« Reply #94 on: June 16, 2011, 07:07:58 AM »
You're right . I'm sorry for excessive omitted .

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #95 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

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

ymg
« Last Edit: June 16, 2011, 04:14:44 PM by CAB »

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #96 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?

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #97 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

TopoWAR

  • Newt
  • Posts: 135
Re: Triangulation (re-visited)
« Reply #98 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.
Thanks for help

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #99 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

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #100 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...

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #101 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.  
« Last Edit: July 02, 2013, 01:48:00 PM by ymg »

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #102 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

ymg

  • Swamp Rat
  • Posts: 725
Re: Triangulation (re-visited)
« Reply #103 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



« Last Edit: July 03, 2013, 10:19:10 PM by ymg »

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1569
  • Moscow (Russia)
Re: Triangulation (re-visited)
« Reply #104 on: July 04, 2013, 04:51:32 AM »