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 :
;*****************************************************************************;
; 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 :
( (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