Author Topic: The minimum circle contain set of points, Thanks for Menzi's function  (Read 4123 times)

0 Members and 1 Guest are viewing this topic.

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Yesterday, talk with some friend about The minimum circle contain set of points. After some time of consideration, I write down the following code.
Code: [Select]
;;; ========================================================================
;;; Thanks to Menzi's great code for find the convex hull of points        ;
;;; ========================================================================
;;; Some of the following codes are writen by QJCHEN                       ;
;;; Civil engineering Department, South China University of Technology     ;
;;; http://autolisper.googlepages.com                                      ;
;;; Purpose: To find the minimum radius circle that contain a set of points;
;;;          in z=0 plane                                                  ;
;;; Algorithm: First, use Menzi's great code to find the convex hull of the;
;;;           points. Then, in the new point list, repeat each three points;
;;;           get the minimum radius circle that contain each three points ;
;;;           then the maximum radius cirlce of this set is the minimum    ;
;;;           radius circle that contain the points.                       ;
;;; The second step algorithm is my own thought, and havent been proved    ;
;;;           strictly, but it seems right, I will find a more strict way  ;
;;; The efficiency of the second step is slowly but after all it can run   ;
;;; Version: 0.1                                                           ;
;;; Platform: R2000 and after                                              ;
;;; Original post :www.Theswamp.org                                        ;
;;; 2006.06.14                                                             ;
;;; Thanks to Menzi's great code about the convex hull which I get from    ;
;;; Google group                                                           ;
;;; His website is http://www.menziengineering.ch/                         ;
;;; ========================================================================

(defun c:test (/ os ssp ssp1 ssp2 mmmax)
  (setq os (getvar "osmode"))
  (setvar "osmode" 0)
  (setvar "pdmode" 2)
  (setq ssp (ssget '((0 . "POINT"))))
  (setq ssp1 (ssgetpoint ssp))
  (setq ssp2 (MeGetConvexHull ssp1))
  (setq mmmax (getmin ssp2))
  (command "circle" (nth 1 mmmax) (nth 0 mmmax))
  (setvar "osmode" os)
)


;;;repeat the convex hull by each 3 points to find the circle
(defun getmin (lst / i n a j b k c mm lstmm mmmax)
  (setq i 0
mmmax (list nil nil)
  )
  (setq n (length lst))
  (repeat (- n 2)
    (setq a (nth i ssp2))
    (setq j (1+ i))
    (repeat (- (- n i) 2)
      (setq b (nth j ssp2))
      (setq k (1+ j))
      (repeat (- (- n j) 1)
(setq c (nth k ssp2))
(setq mm (3pcircle a b c))
(if (> (nth 0 mm) (nth 0 mmmax))
  (setq mmmax mm)
)        
(setq k (1+ k))
      )
      (setq j (1+ j))
    )
    (setq i (1+ i))
  )
  mmmax
)

;;; get the minimum circle of point a b c
(defun 3pcircle (a b c / ab bc ca index lst)
  (setq ab (distance a b))
  (setq bc (distance b c))
  (setq ca (distance a c))

  (if (and
(> (cosabc ab bc ca) 0)
(> (cosabc bc ab ca) 0)
(> (cosabc ca ab bc) 0)
      )
    (setq index 1)
    (setq index 0)
  )

  (if (= index 0)
    (setq lst (longestabc a b c ab bc ca))
  )
  (if (= index 1)
    (setq lst (longestabc1 a b c))
  )
  lst
)

;;; get the angle
(defun cosabc (a b c / d)
  (setq d (- (+ (* b b) (* c c)) (* a a)))
  d
)

;;; get the minimum circle of angle big than 90 degreed.
(defun longestabc (a b c ab bc ca / d r p)
  (setq d (max
    ab
    bc
    ca
  )
r (* 0.5 d)
  )
  (cond
    ((= d ab)
      (setq p (midpoint a b))
    )
    ((= d bc)
      (setq p (midpoint b c))
    )
    ((= d ca)
      (setq p (midpoint c a))
    )
  )
  (list r p)
)

;;; get the minimum circle of angle low than 90 degreed.
(defun longestabc1 (a b c)
  (command "circle" "3p" a b c)
  (setq d (entlast))
  (setq e (entget d))
  (setq p (cdr (assoc 10 e)))
  (setq r (cdr (assoc 40 e)))
  (entdel d)
  (list r p)
)
;;; midpoint function
(defun midpoint (p1 p2)
  (mapcar
    '(lambda (x)
       (/ x 2.)
     )
    (mapcar
      '+
      p1
      p2
    )
  )
)

(defun ssgetpoint (ss / i listpp a b c)
  (setq i 0
listpp nil
  )
  (if ss
    (progn
      (repeat (sslength ss)
(setq a (ssname ss i))
(setq b (entget a))
(setq c (cdr (assoc 10 b)))
(setq listpp (append
       listpp
       (list c)
     )
)
(setq i (1+ i))
      )
    )
  )
  listpp
)
;;;
;;; == Function MeGetConvexHull
;;; Calculates the convex hull from a 'cloud' of points.
;;; Arguments [Type]:
;;;   Lst = Point list [LIST]
;;; Return [Type]:
;;;   > Hull points [LIST]
;;; Notes:
;;;   - Algorithm by Graham
;;;
(defun MeGetConvexHull (Lst / FstCnt LstLen MinCnt NxtCnt TmpLst)
  (setq LstLen (length Lst)
MinCnt 0
FstCnt 1
  )
  (while (< FstCnt LstLen)
    (if (< (cadr (nth FstCnt Lst)) (cadr (nth MinCnt Lst)))
      (setq MinCnt FstCnt)
    )
    (setq FstCnt (1+ FstCnt))
  )
  (setq FstCnt 0)
  (while (< FstCnt LstLen)
    (if (and
  (equal (cadr (nth FstCnt Lst)) (cadr (nth MinCnt Lst)) 1E-8)
  (> (car (nth FstCnt Lst)) (car (nth MinCnt Lst)))
)
      (setq MinCnt FstCnt)
    )
    (setq FstCnt (1+ FstCnt))
  )
  (setq TmpLst (MeSwapList Lst 0 MinCnt)
TmpLst (vl-sort TmpLst (function (lambda (e1 e2)
   (< (MeCalcTheta (nth 0 TmpLst) e1)
      (MeCalcTheta (nth 0 TmpLst) e2)
   )
)
       )
       )
TmpLst (cons (last TmpLst) TmpLst)
NxtCnt 3
FstCnt 4
  )
  (while (<= FstCnt LstLen)
    (while (>= (MeGetCcw (nth NxtCnt TmpLst) (nth (1- NxtCnt) TmpLst)
(nth FstCnt TmpLst)
       ) 0
   )
      (setq NxtCnt (1- NxtCnt))
    )
    (setq NxtCnt (1+ NxtCnt)
  TmpLst (MeSwapList TmpLst FstCnt NxtCnt)
  FstCnt (1+ FstCnt)
    )
  )
  (mapcar
    '(lambda (l)
       (nth l TmpLst)
     )
    (MeRepeatedList 0 (1- NxtCnt) 1)
  )
)
;;;
;;; == Function MeSwapList
;;; Swaps 2 atoms in a list.
;;; Arguments [Type]:
;;;   Lst = List to swap [LIST]
;;;   Fst = Source position [INT]
;;;   Nxt = Target position [INT]
;;; Return [Type]:
;;;   > Swapped list [LIST]
;;; Notes:
;;;   None
;;;
(defun MeSwapList (Lst Fst Nxt / FstVal NxtVal TmpLst)
  (setq FstVal (nth Fst Lst)
NxtVal (nth Nxt Lst)
TmpLst (MeEditList 0 Fst NxtVal Lst)
  )
  (MeEditList 0 Nxt FstVal TmpLst)
)
;;;
;;; == Function MeCalcTheta
;;; Calculates a ordinal number between 2 points.
;;; Arguments [Type]:
;;;   Pt1 = First point [LIST]
;;;   Pt2 = Second point [LIST]
;;; Return [Type]:
;;;   > A value between 0 and 360 [REAL]
;;; Notes:
;;;   None
;;;
(defun MeCalcTheta (Pt1 Pt2 / X__Abs Y__Abs X__Dif Y__Dif TheVal)
  (setq X__Dif (- (car Pt2) (car Pt1))
Y__Dif (- (cadr Pt2) (cadr Pt1))
X__Abs (abs X__Dif)
Y__Abs (abs Y__Dif)
TheVal (if (equal (+ X__Abs Y__Abs) 0 1E-5)
0
(/ Y__Dif (+ X__Abs Y__Abs))
       )
  )
  (if (< X__Dif 0)
    (setq TheVal (- 2.0 TheVal))
    (if (< Y__Dif 0)
      (setq TheVal (+ 4.0 TheVal))
    )
  )
  (* 90.0 TheVal)
)
;;;
;;; == Function MeGetCcw
;;; Determines the direction of 3 points.
;;; Arguments [Type]:
;;;   Pt1 = First point [LIST]
;;;   Pt1 = Second point [LIST]
;;;   Pt3 = Third point [LIST]
;;; Return [Type]:
;;;   >  1 = ccw [INT]
;;;   > -1 = cw [INT]
;;;   >  0 = Colinear [INT]
;;; Notes:
;;;   None
;;;
(defun MeGetCcw (Pt0 Pt1 Pt2 / X1_Dif X1_Sqr X2_Dif X2_Sqr Y1_Dif Y1_Sqr
     Y2_Dif Y2_Sqr
)
  (setq X1_Dif (- (car Pt1) (car Pt0))
Y1_Dif (- (cadr Pt1) (cadr Pt0))
X2_Dif (- (car Pt2) (car Pt0))
Y2_Dif (- (cadr Pt2) (cadr Pt0))
X1_Sqr (* X1_Dif X1_Dif)
Y1_Sqr (* Y1_Dif Y1_Dif)
X2_Sqr (* X2_Dif X2_Dif)
Y2_Sqr (* Y2_Dif Y2_Dif)
  )
  (cond
    ((> (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif))
      1
    )
    ((< (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif))
      -1
    )
    ((or
       (< (* X1_Dif X2_Dif) 0)
       (< (* Y1_Dif Y2_Dif) 0)
     )
      -1
    )
    ((< (+ X1_Sqr Y1_Sqr) (+ X2_Sqr Y2_Sqr))
      1
    )
    (0)
  )
)
;;;
;;; == Function MeRepeatedList
;;; Fills a list with numbers from Srt to End, using increment of Inc.
;;; Arguments [Type]:
;;;   Srt = Start number [INT] or [REAL]
;;;   End = End number [INT] or [REAL]
;;;   Inc = Increment to use [INT] or [REAL]
;;; Return [Type]:
;;;   > List containting all numbers between and including Srt and End
;;; Notes:
;;;   If End is not a repeated of End - Srt, End will not be included
;;;
(defun MeRepeatedList (Srt End Inc / TmpVal RetVal)
  (setq TmpVal (- Srt Inc))
  (while (< TmpVal End)
    (setq RetVal (append
   RetVal
   (list (setq TmpVal (+ TmpVal Inc)))
)
    )
  )
  RetVal
)
;;;
;;; == Function MeEditList
;;; Delete, replace or append a list.
;;; Arguments [Type]:
;;;   Mde = Mode  1 append  [INT]
;;;         Mode  0 replace [INT]
;;;         Mode -1 delete  [INT]
;;;   Pos = Position in list [INT]
;;;   Val = New value [INT/REAL/STR/LIST]
;;;   Lst = List [LIST]
;;; Return [Type]:
;;;   > Modified list [LIST]
;;; Notes:
;;;   - Delete and replace, require the position (Pos) in list.
;;;   - Append and replace, require the new value (Val).
;;;   - :vlax-null items are not allowed in the list argument
;;;
(defun MeEditList (Mde Pos Val Lst / Countr TmpLst TmpVal)
  (if (= Mde 1)
    (reverse (cons Val (reverse Lst)))
    (progn
      (setq Countr -1
    TmpVal (if (= Mde -1)
     :vlax-null
     Val
   )
    TmpLst (mapcar
     '(lambda (l)
(if (= Pos (setq Countr (1+ Countr)))
  TmpVal
  l
)
      )
     Lst
   )
      )
      (vl-remove :vlax-null TmpLst)
    )
  )
)

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

Bob Wahr

  • Guest
I can't see myself ever having a need for this functionality which of course means that I will need it sometime in the next week, but it looks pretty nifty.

Jürg Menzi

  • Swamp Rat
  • Posts: 598
  • Oberegg, Switzerland
Yesterday, talk with some friend about The minimum circle contain set of points. After some time of consideration, I write down the following code.

祝賀, kool implementation of ConvexHull!
A computer's human touch is its unscrupulousness!
MENZI ENGINEERING GmbH
Current A2k16... A2k21 - Start R2.18

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
thank you, Menzi
 :-)
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

Joe Burke

  • Guest
Maybe some food for thought...

The code as posted errors when one or no points are selected.
And it does not work as expected in a UCS. Which IMO, is always
an important concern.

The following functions address those issues, with less code.
I think the net result is the same, though not well tested.

Code: [Select]
;Argument: a selection set
;Returns: a list of VLA objects
(defun SSVLAList (ss / obj lst i)
  (setq i 0)
  (if ss
    (repeat (sslength ss)
      (setq obj (vlax-ename->vla-object (ssname ss i))
            lst (cons obj lst)
            i (1+ i)
      )
    )
  )
  (reverse lst)
) ;end

;; Arguments - a point and a point list.
;; Returns - the point farthest from point.
(defun FarthestPoint (pt ptlst / x dist res)
  (setq x 0)
  (foreach p ptlst
    (setq dist (distance p pt))
    (if (> dist x)
      (setq x dist res p)
    )
  )
  res
) ;end

(defun AveragePts (ptlist / lst)
 (if (= (length (car ptlist)) 2) ;2D point list
   (setq lst
     (list
       (apply '+ (mapcar '(lambda (x) (car x)) ptlist))
       (apply '+ (mapcar '(lambda (x) (cadr x)) ptlist))
     )
   )
   (setq lst ;3D point list
     (list
       (apply '+ (mapcar '(lambda (x) (car x)) ptlist))
       (apply '+ (mapcar '(lambda (x) (cadr x)) ptlist))
       (apply '+ (mapcar '(lambda (x) (caddr x)) ptlist))
     )
   )
 )
 (mapcar '(lambda (x) (/ x (length ptlist) 1.0)) lst)
)

;; Primary function ;;
(defun c:test2 ( / *error* osm ssp ext apt p1 p2 p3 ptlst)

  (defun *error* (msg)
    (cond
      ((not msg))
      ((wcmatch (strcase msg) "*QUIT*,*CANCEL*"))
      (T (princ (strcat "\nError: " msg)))
    )
    (setvar "osmode" osm)
    (princ)
  ) ;end error

  (setq osm (getvar "osmode"))

  (if
    (and
      (setq ssp (SSVLAList (ssget '((0 . "POINT")))))
      (< 1 (length ssp))
      (foreach x ssp
        (setq ptlst (cons (trans (vlax-get x 'Coordinates) 0 1) ptlst))
      )
      (setq apt (AveragePts ptlst))
      (setvar "osmode" 0)
    )
    ;then
    (cond
      ((= 2 (length ptlst))
        (command "circle" "2P" (car ptlst) (cadr ptlst))
      )
      (T
        (setq p1 (FarthestPoint apt ptlst)
              ptlst (vl-remove p1 ptlst)
              p2 (FarthestPoint apt ptlst)
              ptlst (vl-remove p2 ptlst) 
              p3 (FarthestPoint apt ptlst)
        )
        (command "circle" "3P" p1 p2 p3)
      )
    )
    ;else
    (princ "\nInsufficient data to draw circle. ")
  )
  (*error* nil)
) ;end

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Thank you, dear Joe Burke
the first code cant deal with the error of select less 2 points and cant be used in ucs. I will improve it later.

but it seems the solution of problem is not the average point and the longest point among them.
just for the following dwg picture.
the test2 give a smaller circle but cant contain all the points
test give a larger circle but contain all the points.

I need to think it over
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

Joe Burke

  • Guest
I see your point. I'll play with it some more over the weekend.

Regards

highflyingbird

  • Bull Frog
  • Posts: 415
  • Later equals never.

Here is my algorithm.
Code: [Select]
;;;******************************************************
;;; A routine for finding a smallest enclosing circle for
;;; a set of points (SEC problem)                       
;;; command: test                                       
;;; usage:  select some points ,then you will get the SEC
;;;******************************************************
(defun C:test (/ ss pts t0 x cen rad ptmax)
  ;; select the points
  (setq ss (ssget '((0 . "POINT"))))
  (setq pts (ssgetpoint ss))

  ;; get started
  (setq t0 (getvar "TDUSRTIMER"))
  (setq x (mincir pts))
  (princ "\nIt takes")
  (princ (* (- (getvar "TDUSRTIMER") t0) 86400))
  (princ "seconds")
  (if (null x)
    (alert "Invalid input!")
    (progn
      (setq cen   (car x)
    rad   (cadr x)
    ptmax (caddr x)
      )
      ;;draw the smallest circle and the radiu.
      (make-circle cen rad)
      (make-line cen ptmax)
    )
  )
  (princ)
)
;;;*****************************************************
;;;Function : find the SEC with a set of points         
;;;Arguments: pts ,the set of points                   
;;;Return: the center and radius of the SEC             
;;;*****************************************************
(defun mincir (pts / CEN CEN_R P1 P2 P3 PTMAX R rad X i)
  (cond
    ((= (length pts) 0)
     nil
    )
    ((= (length pts) 1)
     (alert "One point,radius is zero.")
     (list (car pts) 0 (car pts))
    )
    ((= (length pts) 2)
     (alert "Two points,radius is the half distance of these points.")
     (setq cen (mid (car pts) (cadr pts))
   rad (/ (distance (car pts) (cadr pts)) 2)
     )
     (list cen rad (car pts))
    )
    (t
     ;;where it works
     (setq p1 (car pts)
   p2 (cadr pts)
   p3 (caddr pts)
     )
     (setq cen_r (3pc p1 p2 p3))
     (setq ptmax (maxd-cir pts (car cen_r)))
     (setq i 0)
     (while (null (in1 ptmax (car cen_r) (cadr cen_r)))
       (setq cen_r (4pc p1 p2 p3 ptmax))
       (setq p1 (car (caddr cen_r))
     p2 (cadr (caddr cen_r))
     p3 (caddr (caddr cen_r))
       )
       (setq ptmax (maxd-cir pts (car cen_r)))
       (setq i (1+ i))
     )
     (list (car cen_r) (cadr cen_r) ptmax)
    )
  )
)
(defun make-circle (cen rad)
  (entmake
    (list
      '(0 . "circle")
      (cons 10 cen)
      (cons 40 rad)
      (cons 62 1)
    )
  )
)
(defun make-line (p q)
  (entmake
    (list
      '(0 . "LINE")
      (cons 10 p)
      (cons 11 q)
    )
  )
)
;;; Gather the coordinates of these points
(defun ssgetpoint (ss / i l a b c)
  (setq i 0)
  (if ss
    (repeat (sslength ss)
      (setq a (ssname ss i))
      (setq i (1+ i))
      (setq b (entget a))
      (setq c (cdr (assoc 10 b)))
      (setq l (cons c l))
    )
  )
  (reverse l)
)
;;; the midpoint of two points
(defun mid (p1 p2)
  (list
    (* (+ (car p1) (car p2)) 0.5)
    (* (+ (cadr p1) (cadr p2)) 0.5)
    (* (+ (caddr p1) (caddr p2)) 0.5)
  )
)
;;; whether a point is inside a circle or not
(defun in1 (pt cen r)
  (< (- (distance pt cen) r) 1e-8)
)
;;; whether some points is is inside a circle or not
(defun in2 (ptl cen r / pts pt)
  (setq pts ptl)
  (while (and (setq pt (car pts))
      (in1 pt cen r)
)
    (setq pts (cdr pts))
  )
  (null pt)
)
;;; Function: get the SEC with 3 points
;;; Arguments: pa pb pc ,the 3 points
;;; Return: the SEC for 3 points.
(defun 3pc (pa pb pc / D MIDPT)
  (cond
    ( (in1 pc (setq midpt (mid pa pb)) (setq d (/ (distance pa pb) 2)))
      (list midpt d (list pa pb pc))
    )
    ( (in1 pa (setq midpt (mid pb pc)) (setq d (/ (distance pb pc) 2)))
      (list midpt d (list pb pc pa))
    )
    ( (in1 pb (setq midpt (mid pc pa)) (setq d (/ (distance pc pa) 2)))
      (list midpt d (list pc pa pb))
    )
    (t
      (3pcircle pa pb pc)
    )
  )
)
;;; Function: Get the circle of a triangle
;;; Arguments: three points of this triangle
;;; Return: the center and radius
(defun 3PCirCle(P0 P1 P2 / X0 Y0 X1 Y1 X2 Y2 DX1 DY1 DX2 DY2 D 2D C1 C2 CE)
  (setq X0  (car P0)
Y0  (cadr P0)
X1  (car P1)
Y1  (cadr P1)
X2  (car P2)
Y2  (cadr P2)
DX1 (- X1 X0)
DY1 (- Y1 Y0)
DX2 (- X2 X0)
DY2 (- Y2 Y0)
  )
  (setq D (- (* DX1 DY2) (* DX2 DY1)))
  (if (/= D 0.0)
    (progn
      (setq 2D (+ D D)
    C1 (+ (* DX1 (+ X0 X1)) (* DY1 (+ Y0 Y1)))
    C2 (+ (* DX2 (+ X0 X2)) (* DY2 (+ Y0 Y2)))
    CE (List (/ (- (* C1 DY2) (* C2 DY1)) 2D)
     (/ (- (* C2 DX1) (* C1 DX2)) 2D)
       )
      )
      (list CE (distance CE P0) (list p0 p1 p2))
    )
  )
)
;;; Function: get the SEC with 4 points
;;; Arguments: pa pb pc ptmax,the 4 points
;;; Return: the SEC for 4 points.
(defun 4pc (p1 p2 p3 ptmax / pts mind minr r 4ps)
  (setq pts (list (3pc p1 p2 ptmax)
  (3pc p1 p3 ptmax)
  (3pc p2 p3 ptmax)
    )
  )
  (setq 4ps (list p1 p2 p3 ptmax))
  (setq minr 1e308)
  (foreach n pts
    (setq r (cadr n))
    (if (and (< r minr)
     (in2 4ps (car n) r)
)
      (setq mind n)
    )
  )
  mind
)

;;; Function: Get the farthest point from a center.
;;; Arguments: ptl,the points
;;;            cen,the center
;;; Return: the farthest point
(defun maxd-cir (ptl cen / pmax dmax d)
  (setq dmax 0.0)
  (foreach pt ptl
    (if (> (setq d (distance pt cen)) dmax)
      (setq dmax d
    pmax pt
      )
    )
  )
  pmax
)
I am a bilingualist,Chinese and Chinglish.