;| ConvexHull3D.lsp - fast 3D Least Convex Hull implemented in AutoLISP CH3D calculates the least convex hull of a set of 3D point entities, and provides multiple output options: a 3DSOLID, facets, edges, struts, and/or joints. Duplicate points are ignored. Points internal to the hull, to an edge, or to a face of the hull are ignored because they do not contribute to the definition of the least convex hull. Polygonal faces with more than three sides are treated as a single facet, and are not subdivided into arbitrary triangular sections. Several methods have been used to optimize performance by removing from consideration points interior to the hull before final computation, resulting in execution speed several hundred times faster than without optimization. Even with optimization, computation time still increases 1) exponentially with large numbers of points and 2) with deviation of the cloud from a rectangular shape. Sample times with elliptical clouds on a 3.0 GHz CPU and AutoCAD 2011: points seconds ------ ------- 25 0.4 50 1.1 100 4.1 200 12.1 400 60.5 800 184.0 1600 2745.0 NOTE: this is a highly computationally-intensive routine, and runs fastest when AutoCAD is RUNNING ALONE, and when the AutoCAD window HAS FOCUS. If other applications are active and have focus, times can be several times as long. Note that internally CH3D uses a list of lists of 3D coordinates instead of a selection set, and so could readily be adapted to use points from a text file instead of actual AutoCAD point entities. Thanks to: Somebuddy for (remove_doubles) gile for the dot product functions Reini Urban for (remove) Martti Halminen for ways to reduce computational overhead and Nicole at Clarkson for suggesting the update. v 1.00 Sep 06 2008 - published on autodesk.autocad.customization v 2.00 Jan 11 2011 - completely rewritten - polygonal facets supported - output options: solid (3DSOLID) facets (regions) edges (lines) struts (extruded cylinders) joints (spheres) - provides a separate layer for each option - reports number of requested output results - reports volume of hull if solid option chosen v 2.01 Jan 15 2011 - modest performance improvement v 2.02 Jan 16 2011 - further performance improvement v 2.03 Jan 17 2011 - 4x performance improvement w/ dot products v 3.00 Jan 18 2011 - 45x(!) performance improvement by pruning interior points v 3.01 Jan 22 2011 - stability improvements - code modularized - optimized for different numbers of points v 3.02 Jan 24 2011 - performance improvement with large clouds. The larger the cloud, the greater the increase: 400 points: 2.5x 800 points: 4x Check for updates at www.realerthanreal.com/autolisp To generate clouds of points, try RandomPoints.lsp, also found there. (c) 2011 by Bill Gilliss bill [atttt] realerthanreal [dotttt] com Comments and suggestions always welcome. ===================================================================== |; (defun c:CH3D ( / ;; ===== system variables *cmdecho olderror *osmode *3dosmode *delobj *clayer *ucsname *gridmode *ucsicon *navvcubedisplay ;; ===== subroutine functions myerror prologue cleanup getcent getExtPoints date2sec startTimer endTimer setSubtract setCopy setPurge layernames remove_doubles dotProduct crossProduct normal3Points normalize ilp delDupEdges ;; ===== local variables, sets, lists, and entities facetLayer facetColor edgeLayer edgeColor strutLayer strutColor jointLayer jointColor solidLayer solidColor hullPoints hullFaces hullFacets hullEdges hullSolids dupEdges shortList facets edges struts joints diam solids solidVolume inputFlag pt ptlist n m a b c d p1 p2 p3 p4 fac cent intPoint jointPoints extPoints extLength skipFace allObjects strutCircle jointSphere centroidEntity dummyPt timer1 timer2 elapsed finalHullFaces finalHull pts listlength iterations ) ;; ================== PROGRAM FUNCTIONS ======================= ;; Do not change layer names. Feel free to change layer colors. (defun setLayerNames () (setq facetLayer "CH3D-facets" edgeLayer "CH3D-edges" strutLayer "CH3D-struts" jointLayer "CH3D-joints" solidLayer "CH3D-solids" facetColor "_green" edgeColor "_yellow" strutColor "_blue" jointColor "_red" solidColor 9 ) ) (defun myerror (msg) (princ msg) (cleanup) ) (defun prologue () (setq *cmdecho (getvar "cmdecho")) (setvar "cmdecho" 0) (setq olderror *error*) (setq *error* myerror) (setq *osmode (getvar "osmode")) (setvar "osmode" 16384) ;;off (setq *3dosmode (getvar "3dosmode")) (setvar "3Dosmode" 1) ;;off (setq *delobj (getvar "delobj")) (setvar "delobj" 0) (setq *clayer (getvar "clayer")) (setq *ucsname (getvar "ucsname")) (setq *gridmode (getvar "gridmode")) (setq *ucsicon (getvar "ucsicon")) (setq *navvcubedisplay (getvar "navvcubedisplay")) (vl-load-com) (arxload "geomcal") (command "._undo" "_begin") (command "._layer" "_unlock" "CH3D*" "") (setq ptlist (list) hullFaces (list) hullFacets (ssadd) hullEdges (ssadd) hullPoints (list) hullSolids (ssadd) jointPoints (list) extPoints (list) shortList (list) finalHullFaces (ssadd) ) ) (defun cleanup () (setq *error* olderror) (setvar "osmode" *osmode) (setvar "3dosmode" *3dosmode) (setvar "delobj" *delobj) (setvar "clayer" *clayer) (setvar "gridmode" *gridmode) (setvar "ucsicon" *ucsicon) (setvar "navvcubedisplay" *navvcubedisplay) (if (not edges) (command "._erase" hullEdges "")) (command "._erase" finalHullfaces "") (command "._undo" "_end") (setvar "cmdecho" *cmdecho) ) (defun getUserInput () (prompt "Select points for 3D convex hull: ") (setq inputFlag nil) (while (not inputFlag) (setq pts (ssget '((0 . "POINT")))) (if (< (sslength pts) 4) (alert "Need at least 4 unique points.") (setq inputFlag T) ) ) (princ "\nCH3D can provide any combination of a 3Dsolid, facets,") (princ "\nedges, cylindrical struts, and/or spherical joints.\n") (setq inputFlag nil) (while (not inputFlag) (princ "Do you want...") (initget 0 "Yes No") (setq solids (getkword "\nA solid? [Yes or No] : ")) (if (eq solids "Yes") (setq solids T inputFlag T) (setq solids nil)) (initget 0 "Yes No") (setq facets (getkword "Facets? [Yes or No] : ")) (if (eq facets "Yes") (setq facets T inputFlag T) (setq facets nil)) (initget 0 "Yes No") (setq edges (getkword "Edges? [Yes or No] : ")) (if (eq edges "Yes") (setq edges T inputFlag T) (setq edges nil)) (initget 0 "Yes No") (setq struts (getkword "Struts? [Yes or No] : ")) (if (eq struts "Yes") (setq struts T inputFlag T) (setq struts nil)) (initget 0 "Yes No") (setq joints (getkword "Joints? [Yes or No] : ")) (if (eq joints "Yes") (setq joints T inputFlag T) (setq joints nil)) (if (not inputFlag) (progn (vlr-beep-reaction) (princ "\nGotta choose something. ") ) ) ) (if struts (progn ;;leave strutdiam global for reuse (if (not (numberp strutdiam)) (setq strutdiam 1.0)) (initget 6) ;;no zero, no negative (setq diam (getreal (strcat "Enter strut diameter: <"(rtos strutdiam)">: ")) ) (if diam (setq strutdiam diam)) ) ) (if joints (progn ;;leave jointdiam global for reuse (if (not (numberp jointdiam)) (setq jointdiam 1.0)) (initget 6) (setq diam (getreal (strcat "Enter joint diameter: <"(rtos jointdiam)">: ")) ) (if diam (setq jointdiam diam)) ) ) (princ "\nCalculating...\n") );getUserInput ;; convert set of points to list for lower computational overhead (defun pointsToList () (setq n 0) (while (< n (sslength pts)) (setq pt (cdr (assoc 10 (entget (ssname pts n))))) (setq ptlist (cons pt ptlist)) (setq n (1+ n)) ) (setq ptlist (remove_doubles ptlist nil)) (if (> 4 (length ptlist)) (progn (alert "Need at least 4 unique points.") (quit) ) ) ) ;; Find approximate centroid of hull, to be used as a known interior ;; point when evaluating other points. Also, surface faces will be ;; lofted to this point to create 3Dsolids. Note that this is not the ;; actual mathematical centroid of the final convex hull. (defun getcent ( / n pt x y z cx cy cz) (setq cX 0 cY 0 cZ 0) (foreach pt ptlist (setq x (car pt) y (cadr pt) z (caddr pt) cX (+ cX x) cY (+ cY y) cZ (+ cZ z) ) ) (setq cX (/ cX (length ptlist)) cY (/ cY (length ptlist)) cZ (/ cZ (length ptlist)) ) (setq cent (list cX cY cZ)) ;; need ann interior dummy point at end of list to allow last actual ;; point to be evaluated (setq pt (nth 0 ptlist)) ;;could be any point in list (setq dummy (c:cal "(pt + cent)/2")) (setq ptlist (reverse (cons dummy (reverse ptlist))) listLength (length ptlist) ) (if solids (progn (entmake (list (cons 0 "POINT") (cons 10 cent))) (setq centroidEntity (entlast)) ;;need entity for LOFT ) ) );;getcent ;| Since points known to be on the hull are fewer and farther from the centroid than the interior points, they should (and do) disqualify potential surface faces much faster than using the entire cloud. This function finds such extreme points, those with minimum and maximum x, y, and z values relative to several different UCS's. The optimum number of UCS's to use proved to increase with the number of points, approximately proportional to the square of the log (base 10) of the number of points for irregularly shaped clouds with uniform distribution. |; (defun getExtPoints ( / orig pt pt2 x y z maxx minx maxy miny maxz minz maxxPt minxPt maxyPt minyPt maxzPt minzPt) (setvar "gridmode" 0) ;;to prevent screen flicker during UCS changes (setvar "ucsicon" 0) (setvar "navvcubedisplay" 0) (if (= (getvar "ucsname") "") (progn (command "._UCS" "_named" "_delete" "CH3D") (command "._UCS" "_named" "_save" "CH3D") (setq *ucsname "CH3D") ) ) (setq iterations ;;at least three UCS's (max 3 (fix (expt (/ (log (sslength pts)) (log 10)) 2))) ) (setq rotang (/ 90.0 iterations)) (repeat iterations (repeat iterations (repeat iterations (setq orig (trans (car ptlist) 0 1)) ;; any point in the cloud (mapcar '(lambda (x) (set x orig)) ;; initialize local maxima (list 'maxX 'minX 'maxY 'minY 'maxZ 'minZ) ) (foreach pt ptlist ;; WCS coords (setq pt2 (trans pt 0 1)) ;; UCS coords (setq x (car pt2) y (cadr pt2) z (caddr pt2)) (if (> x (car maxx)) (setq maxx pt2 maxxPt pt) (if (< x (car minx)) (setq minx pt2 minxPt pt) ) ) (if (> y (cadr maxy)) (setq maxy pt2 maxyPt pt) (if (< y (cadr miny)) (setq miny pt2 minyPt pt) ) ) (if (> z (caddr maxz)) (setq maxz pt2 maxzPt pt) (if (< z (caddr minz)) (setq minz pt2 minzPt pt) ) ) );foreach (foreach p (list maxxPt minxPt maxyPt minyPt maxzPt minzPt) ;;WCS coords (if (not (member p extPoints)) (setq extpoints (cons p extpoints)) ) ) (command "._ucs" "x" rotang) ) (command "._ucs" "y" rotang) ) (command "._ucs" "z" rotang) );repeat (setq extPoints (remove_doubles extPoints nil)) (setq extpoints (remove nil extpoints)) (setq extpoints (reverse (cons dummy (reverse extpoints)))) (setq extlength (length extpoints)) (if (not (eq *ucsname "CH3D")) (command "._UCS" "_named" "_restore" *ucsname) (command "._UCS" "_named" "_restore" "CH3D") ) (setvar "gridmode" *gridmode) (setvar "ucsicon" *ucsicon) (setvar "navvcubedisplay" *navvcubedisplay) );extpoints (defun makeFacetLayer () (if (tblsearch "layer" facetLayer) (setvar "clayer" facetLayer) (command "._layer" "_make" facetLayer "_color" facetColor facetLayer "") ) ) ;; ======== Build hull of 3-point faces ;; Output depends upon when function is called. ;; For each potential face (combination of three points), we have to ;; check that every other point in the list is on the same side of the ;; face as the centroid. Takes a while... (defun makeHull ( pointlist / listlength a b c d p1 p2 p3 p4) (command "._UCS" "_w") ;; so inner loop does not have to use [trans] (setq listLength (length pointlist)) (setq list1 pointlist) (while (cdddr list1) (setq p1 (car list1)) (setq list2 (cdr list1)) (while (cddr list2) (setq p2 (car list2)) (setq list3 (cdr list2)) (while (cdr list3) (setq p3 (car list3)) (if (inters p1 p2 p1 p3) ;;if 3 define a plane (progn (setq skipFace nil) (setq n 0) (while (< n listlength) (setq p4 (nth n pointlist)) (if (and (/= p4 p1) (/= p4 p2) (/= p4 p3) (not (inters p1 p2 p3 p4 nil));;and 4 not co-planar ) (progn (setq planeNormal (normal3Points p1 p2 p3)) (setq intPoint (ilp cent p4 p1 planeNormal)) (if (equal (distance p4 cent) (+ (distance p4 intPoint) (distance intPoint cent) ) 1e-6 ) (setq skipFace T) ) (if (and (= n (1- listLength)) (not skipFace)) (cond ( (not finalHull) (setq hullFaces (cons (list p1 p2 p3) hullFaces)) ) ( T (progn (entmake (list (cons 0 "3DFACE") (cons 10 p1) (cons 11 p2) (cons 12 p3) (cons 13 p3) ) ) (ssadd (entlast) finalHullFaces) (command "._redraw") ) ) ) );if ) );if (if skipface (setq n listlength) (setq n (1+ n))) );p4 ) ) (setq list3 (cdr list3)) );p3 (setq list2 (cdr list2)) );p2 (setq list1 (cdr list1)) );p1 );makeHull ;; Create short list of possible outer points by discarding inner ;; ones. For each point in the cloud, determine if it is outside at ;; least one face of the temporary hull. If so, add it to the short list ;; of outer points. Only these will be used to compute the final hull. ;; This only takes a fraction of a second, since the faces of the temp- ;; orary hull are known. (defun makeShortList () (foreach p4 ptlist (if (not (member p4 extpoints)) ;;i.e., not already an outer point (progn (setq n 0) (while (< n (length hullFaces)) (setq fac (nth n hullFaces) p1 (car fac) p2 (cadr fac) p3 (caddr fac) planeNormal (normal3Points p1 p2 p3) intPoint (ilp cent p4 p1 planeNormal) ) (if (equal (distance p4 cent) (+ (distance p4 intPoint) (distance intPoint cent) ) 1e-6 ) (progn (setq shortList (cons p4 shortList)) (setq n (length hullFaces)) ;;skip to next point ) (setq n (1+ n)) ;; proceed to next face ) );while n ) ) ) (setq shortlist (append shortlist extpoints)) ) ;; Convert 3Dfaces to facets [regions] and union them to ;; merge coplanar faces and eliminate superfluous edges (defun createFacets () (setq allObjects (ssget "X")) (setvar "delobj" 1) ;;to get rid of 3dfaces (command "._region" finalHullFaces "" ) (setq hullFacets (setSubtract (list (ssget "X") allObjects))) (command "._union" hullFacets "") (setq hullFacets (setPurge hullFacets)) ) (defun createSolids () (setvar "delobj" 0) ;;leave faces and centroidEntity (if (tblsearch "layer" solidLayer) (setvar "clayer" solidLayer) (command "._layer" "_make" solidLayer "_color" solidColor solidLayer "") ) (setq n 0) (while (< n (sslength finalHullFaces)) (command "._loft" (ssname finalHullFaces n) centroidEntity "" "_bulge" 0 "" ) (ssadd (entlast) hullSolids) (setq n (1+ n)) ) (command "._union" hullSolids "") (entdel centroidEntity) (setq solidVolume (rtos (/ (vlax-get-property (vlax-ename->vla-object (entlast)) 'volume)))) ) (defun createEdges () (if facets (command "._copy" hullFacets "" '(0 0 0) '(0 0 0)) ) (if (not (tblsearch "layer" edgeLayer)) (command "._layer" "_make" edgeLayer "_color" edgeColor edgeLayer "") ) (setq allObjects (ssget "X")) (setq n 0) (while (< n (sslength hullFacets)) (command "._explode" (ssname hullFacets n)) (setq n (1+ n)) ) (setq hullEdges (setSubtract (list (ssget "X") allObjects))) (command "._change" hullEdges "" "_properties" "_layer" edgelayer "") (delDupEdges) ) ;; ==== get rid of duplicate edges [OVERKILL might not be available] (defun delDupEdges ( / n m en1 ed1 start1 end1 en2 ed2 start2 end2) (setq dupEdges (ssadd)) (setq n 0) (while (< n (sslength hullEdges)) (setq en1 (ssname hullEdges n) ed1 (entget en1) start1 (cdr (assoc 10 ed1)) end1 (cdr (assoc 11 ed1)) ) (setq m (1+ n)) (while (< m (sslength hullEdges)) (setq en2 (ssname hullEdges m) ed2 (entget en2) start2 (cdr (assoc 10 ed2)) end2 (cdr (assoc 11 ed2)) ) (if (or (and (equal start1 start2) (equal end1 end2)) (and (equal start1 end2) (equal end1 start2)) ) (progn (ssadd en2 dupEdges) (setq m (sslength hullEdges)) ) ) (setq m (1+ m)) );while m (setq n (1+ n)) );while n (setq hullEdges (setSubtract (list hullEdges dupEdges))) (setq hullEdges (setPurge hullEdges)) (command "._erase" dupEdges "") ) ;; ======== create struts (defun createStruts () (if (not (tblsearch "LAYER" strutlayer)) (command "._layer" "_make" strutLayer "_color" strutColor strutLayer "") (setvar "clayer" strutlayer) ) (command "._circle" '(0 0 0) (/ strutdiam 2.0)) (setq strutCircle (entlast)) (setvar "delobj" 0) (setq n 0) (while (< n (sslength hullEdges)) (command "._sweep" strutCircle "" (ssname hullEdges n)) (setq n (1+ n)) ) (entdel strutCircle) ) ;; ========= create points and joints (defun createJoints () (setq n 0) (while (< n (sslength hullEdges)) (setq en1 (ssname hullEdges n) ed1 (entget en1) start1 (cdr (assoc 10 ed1)) end1 (cdr (assoc 11 ed1)) ) (setq jointPoints (cons start1 (cons end1 jointPoints))) (setq n (1+ n)) ) (setq jointPoints (remove_doubles jointPoints nil)) (if (not (tblsearch "LAYER" jointlayer)) (command "._layer" "_make" jointLayer "_color" jointColor jointLayer "") (setvar "clayer" jointlayer) ) (command "._sphere" '(0 0 0) (/ jointdiam 2.0)) (setq jointSphere (entlast)) (foreach pt jointPoints (command "._copy" jointSphere "" '(0 0 0) (trans pt 0 1)) ) (entdel jointSphere) ) (defun showResults () (setq elapsed (rtos (setq elapsed (- timer2 timer1)) 2 2)) (princ (strcat "\nTime in loop: " elapsed " seconds")) (if facets (princ (strcat "\nFacets created: " (itoa (sslength hullFacets))))) (if edges (princ (strcat "\nEdges created: " (itoa (sslength hullEdges ))))) (if struts (princ (strcat "\nStruts created: " (itoa (sslength hullEdges ))))) (if joints (princ (strcat "\nJoints created: " (itoa (length jointPoints))))) (if solids (princ (strcat "\nVolume of hull: " solidVolume))) (princ) ) ;;=============== UTILITY FUNCTIONS =========================== (defun date2sec ( / s seconds) (setq s (getvar "DATE")) (setq seconds (* 86400.0 (- s (fix s)))) ) (defun startTimer ( ) (setq timer1 (date2sec)) ) (defun endTimer ( ) (setq timer2 (date2sec)) ) (defun SetSubtract (setlist / n ct setname result) ;;by Bill Gilliss (setq result (SetCopy (nth 0 setlist))) (if (= (type result) 'VLA-OBJECT) (setq result (Set->AL result))) (setq n 1) (while (< n (length setlist)) (setq ct 0 setname (nth n setlist)) (if (= (type setname) 'VLA-OBJECT)(setq setname (Set->AL setname))) (repeat (sslength setname) (ssdel (ssname setname ct) result) (setq ct (1+ ct)) ) (setq n (1+ n)) ) result ) (defun SetCopy (setname / n result) ;;by Bill Gilliss (if (= (type setname) 'VLA-OBJECT) (setq setname (Set->AL setname))) (setq result (ssadd)) (setq n 0) (while (< n (sslength setname)) (ssadd (ssname setname n) result) (setq n (1+ n)) ) result ) (defun SetPurge (setname / n result) ;;by Bill Gilliss (setq result (ssadd)) (setq n 0) (while (< n (sslength setname)) (setq en (ssname setname n)) (if (entupd en) ;;can't update a non-entity (ssadd en result) ) (setq n (1+ n)) ) result ) ;; DotProduct (gile) ;; Returns the dot product (scalar product) of two vectors ;; ;; Arguments: two vectors (defun DotProduct (v1 v2) (apply '+ (mapcar '* v1 v2))) ;; CrossProduct (gile) ;; Returns the cross product of two vectors ;; ;; Arguments: two vectors (defun CrossProduct (v1 v2) (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2))) (- (* (caddr v1) (car v2)) (* (car v1) (caddr v2))) (- (* (car v1) (cadr v2)) (* (cadr v1) (car v2))) ) ) ;; Normalize (gile) ;; Returns the single unit vector of a vector ;; ;; Argument: a vecteur (defun Normalize (v) ((lambda (l) (if (/= 0 l) (mapcar (function (lambda (x) (/ x l))) v) ) ) (distance '(0 0 0) v) ) ) ;; Normal3points (gile) ;; Returns the normal vector of the plane defined by 3 points ;; ;; Arguments: three points (defun Normal3points (p0 p1 p2) (Normalize (CrossProduct (mapcar '- p1 p0) (mapcar '- p2 p0))) ) ;; IntersectLinePlane (gile) ;; Returns the intersection point between an unbounded line and a plane ;; ;; Arguments ;; p1 p2: two points defining the line ;; org: a point on the plane ;; nor: the plane normal vector (defun ilp (p1 p2 org nor / scl) (if (and (/= 0 (setq scl (DotProduct nor (mapcar '- p2 p1)))) (setq scl (/ (DotProduct nor (mapcar '- p1 org)) scl)) ) (mapcar (function (lambda (x1 x2) (+ (* scl (- x1 x2)) x1))) p1 p2 ) ) ) ;;; REMOVE - Non-destructive way to remove an item from a list/tree. ;;; Reini Urban ;;; recursive, double elements allowed (defun remove (item from) (cond ((atom from) from) ((equal (car from) item) (remove item (cdr from))) (T (cons (car from) (remove item (cdr from))))) ) (defun remove_doubles (lst fuzz) ;;by SomeBuddy (cond ((atom lst) lst) (T (cons (car lst) (remove_doubles (if fuzz (vl-remove-if '(lambda (x)(equal (car lst) x fuzz)) lst) (vl-remove (car lst) lst) ) fuzz ) ) ) ) ) ;; ================= MAIN PROGRAM STARTS HERE ======================= (prologue) (setLayerNames) (getUserInput) (pointsToList) (startTimer) ;;time all computational loops (getcent) (makeFacetLayer) (getExtPoints) (setq finalHull nil) (makeHull extPoints) ;;create temporary hull (makeShortList) (setq finalHull T) (makeHull shortList) ;;create final hull (endTimer) ;;stop timer before creating output entities (if solids (createSolids)) (if (or facets edges struts joints) (createFacets)) (if (or edges struts joints) (createEdges)) (if struts (createStruts)) (if joints (createJoints)) (cleanup) (showResults) );main defun (princ "\nEnter CH3D to run program.") (princ)