Topic: Fit Circle to a Selection Set of Points

[ymg]

Probably been done before but I could not find it in the usual sources.

Would appreciate if somebody with access to C3D would check it, and let me know.

For Derivation of formulas see: Least-Squares Circle Fit

ymg

Code - Auto/Visual Lisp: [Select]

1. ;; fitcircle    by ymg                                                        ;
2. ;; Fit a Circle to a Selection Set of Points                                  ;
3. ;;                                                                            ;
4. ;; For Derivation of Formulas see:                                            ;
5. ;; http://www.dtcenter.org/met/users/docs/write_ups/circle_fit.pdf            ;
6. ;;           by R. Bullock                                                    ;
7.  
8.  
9. (defun c:fc (/ *AcadDoc* *error* a avg b cen cu cv dif diff errl m n p pl rad res s ss sum x)
10.    (vl-load-com)
11.    (defun *error* (msg)
12.         (mapcar 'eval errl)
13.         (if (and msg (not (wcmatch (strcase msg) "*BREAK*,*CANCEL*,*EXIT*")))
14.            (princ (strcat "\nError: " msg))
15.         )
16.         (and *AcadDoc* (vla-endundomark *AcadDoc*))
17.         (princ)
18.    )
19.      
20.    (setq errl '("CMDECHO" "DIMZIN")
21.          errl (mapcar (function (lambda (a) (list 'setvar a (getvar a)))) errl)
22.    )
23.      
24.      
25.    (or *AcadDoc*
26.       (setq *AcadDoc* (vla-get-activedocument (vlax-get-acad-object)))
27.    )
28.    
29.    (princ (strcat "\n Select Points"))
30.    (if (setq  n 0
31.              pl nil
32.              ss (ssget '((0 . "POINT")))             
33.        )
34.       (progn
35.          (vla-startundomark *AcadDoc*)     
36.          (setvar 'CMDECHO 0)
37.          (setvar 'DIMZIN 0)
38.          (repeat (sslength ss)
39.                  (setq pl (cons (cdr (assoc 10 (entget (ssname ss n)))) pl)
40.                         n (1+ n)
41.                  )
42.          )
43.          (setq  avg  (mapcar '/  (apply 'mapcar (cons '+ pl)) (list n n n))
44.                 dif  (mapcar '(lambda (a) (mapcar '- a avg)) pl)
45.                 sum  (apply
46.                         'mapcar
47.                            (cons '+
48.                               (mapcar
49.                                  '(lambda (a)
50.                                       (list (* (car a) (car a))
51.                                             (* (car a) (cadr a))
52.                                             (* (cadr a) (cadr a))
53.                                             (* (car a) (car a) (car a))
54.                                             (* (car a) (cadr a) (cadr a))
55.                                             (* (car a) (car a) (cadr a))
56.                                             (* (cadr a) (cadr a) (cadr a))
57.                                       )
58.                                   )
59.                                   dif
60.                               )
61.                            )
62.                      )
63.                 ;; sum contains:  (sUU sUV sVV sUUU sUVV sUUV sVVV)        ;
64.                    a (list (car sum) (cadr sum) (/ (+ (cadddr sum) (car (cddddr sum))) 2))
65.                    b (list (cadr sum) (caddr sum) (/ (+ (cadr (cddddr sum))(caddr (cddddr sum))) 2))
66.                    m (/ (car a) (car b))
67.                    s (mapcar '- a (mapcar '(lambda (x) (* x m)) b))
68.                   cv (/ (caddr s) (cadr s))
69.                   cu (/ (- (caddr a) (* (cadr a) cv)) (car a))
70.                  cen (mapcar '+ (list cu cv) avg)
71.                  rad (sqrt (+ (* cu cu) (* cv cv) (/ (+ (car sum) (caddr sum)) n)))
72.                    n (1- n)
73.          )
74.          (entmakex (list '(0 . "POINT") (cons 10 cen)))
75.          (entmakex (list '(0 . "CIRCLE") (cons 10 cen) (cons 40 rad)))
76.          (while (> n -1)
77.              (setq   p (cdr (assoc 10 (entget (ssname ss n))))
78.                    res (- rad (distance p cen))
79.                      n (1- n)
80.              )      
81.              (entmakex (list '(0 . "TEXT") '(8 . "Residuals") (cons 10 p) (cons 40 3.0) '(50 . 1.5708) (cons 1 (rtos res 2 3))))
82.          )
83.       )
84.    )
85.    (*error* nil)
86. )
87.  
88. (princ "\nFit a Circle to Selection Set of Points, Type FC to run")
89.  

[edata]

Thanks in advance.

When the number of points is 1 of the code will go Eorror.

Code: [Select]

(if (setq  n 0
     pl nil
     ss (ssget '((0 . "POINT")))      
       )

ep.

Code: [Select]

(if (and (setq  n 0
     pl nil
     ss (ssget '((0 . "POINT"))))
  (if ss (if (> (sslength ss) 1) t) t))
  (progn
    ;;progn code

    ))


[ymg]

When the number of points is 1 of the code will go Eorror.

You are right, but no need for a double if. To do a fit you need at least 3 points.

You will also get funny results if all your points are colinear.

Code - Auto/Visual Lisp: [Select]

1. (if (setq  n 0
2.              pl nil
3.              ss (ssget '((0 . "POINT")))
4.              go (> (sslength ss) 2)
5.        )
6.       (progn
7.           ...
8.           ...
9.  

ymg

[Lee Mac]

Code - Auto/Visual Lisp: [Select]

1. (if (setq  n 0
2.              pl nil
3.              ss (ssget '((0 . "POINT")))
4.              go (> (sslength ss) 2)
5.        )
6.       (progn
7.           ...
8.  

(sslength nil) will error 

[ymg]

(sslength nil) will error 

So... will pulling the plug! 

ymg

[Lee Mac]

(sslength nil) will error 

So... will pulling the plug! 

If that is indeed your attitude, why bother to include any error trapping at all...?

[Krushert]

If that is indeed your attitude, why bother to include any error trapping at all...?

Is that like going commando?      

[ymg]

There U go, Lee.

Did it on the batteries. 

But on a serious note, would dearly like to see a comparison of results with Civil 3d.

ymg

Code: [Select]

(princ (strcat "\n Select Points"))
   (if (and (setq ss (ssget '((0 . "POINT"))))
            (> (sslength ss) 2)
       )     
      (progn
         (vla-startundomark *AcadDoc*)    
         (setvar 'CMDECHO 0)
         (setvar 'DIMZIN 0)
         
         (setq  n 0  pl nil)       
(repeat (sslength ss)
         (setq pl (cons (cdr (assoc 10 (entget (ssname ss n)))) pl)


[ymg]

Here is a slightly better version of the Fit circle routine.

This one does a Non Linear regression and thus really minimize the distance to the circle.

The linearized form sumitted before is used to obtain an approach position.

Then the non linear regression is applied until correction is smaller than 1e-06,
and make use of Evgenyi's implementation of Gaussian Elimination.

ymg

Code - Auto/Visual Lisp: [Select]

1. ;; fitcircle    by ymg                                                        ;
2. ;; Fit a Circle to a Selection Set of Points                                  ;
3. ;;                                                                            ;
4. ;; For Derivation of Formulas see:                                            ;
5. ;; http://www.dtcenter.org/met/users/docs/write_ups/circle_fit.pdf            ;
6. ;;           by R. Bullock                                                    ;
7.  
8.  
9. (defun c:fc (/ *AcadDoc* *error* a avg b cen cu cv dif diff errl m n p pl rad res s ss sum x)
10.    (vl-load-com)
11.    (defun *error* (msg)
12.         (mapcar 'eval errl)
13.         (if (and msg (not (wcmatch (strcase msg) "*BREAK*,*CANCEL*,*EXIT*")))
14.            (princ (strcat "\nError: " msg))
15.         )
16.         (and *AcadDoc* (vla-endundomark *AcadDoc*))
17.         (princ)
18.    )
19.      
20.    (setq errl '("CMDECHO" "DIMZIN")
21.          errl (mapcar (function (lambda (a) (list 'setvar a (getvar a)))) errl)
22.    )
23.      
24.      
25.    (or *AcadDoc*
26.       (setq *AcadDoc* (vla-get-activedocument (vlax-get-acad-object)))
27.    )
28.    
29.    (princ (strcat "\n Select Points"))
30.    (if (setq  n 0
31.              pl nil
32.              ss (ssget '((0 . "POINT")))
33.              go (> (sslength ss) 2)
34.        )
35.       (progn
36.          (vla-startundomark *AcadDoc*)     
37.          (setvar 'CMDECHO 0)
38.          (setvar 'DIMZIN 0)
39.          (repeat (sslength ss)
40.                  (setq pl (cons (cdr (assoc 10 (entget (ssname ss n)))) pl)
41.                         n (1+ n)
42.                  )
43.          )
44.          (setq  avg  (mapcar '/  (apply 'mapcar (cons '+ pl)) (list n n n))
45.                 dif  (mapcar '(lambda (a) (mapcar '- a avg)) pl)
46.                 sum  (apply
47.                         'mapcar
48.                            (cons '+
49.                               (mapcar
50.                                  '(lambda (a)
51.                                       (list (* (car a) (car a))
52.                                             (* (car a) (cadr a))
53.                                             (* (cadr a) (cadr a))
54.                                             (* (car a) (car a) (car a))
55.                                             (* (car a) (cadr a) (cadr a))
56.                                             (* (car a) (car a) (cadr a))
57.                                             (* (cadr a) (cadr a) (cadr a))
58.                                       )
59.                                   )
60.                                   dif
61.                               )
62.                            )
63.                      )
64.                 ;; sum contains:  (sUU sUV sVV sUUU sUVV sUUV sVVV)        ;
65.                    a (list (car sum) (cadr sum) (/ (+ (cadddr sum) (car (cddddr sum))) 2))
66.                    b (list (cadr sum) (caddr sum) (/ (+ (cadr (cddddr sum))(caddr (cddddr sum))) 2))
67.                    m (/ (car a) (car b))
68.                    s (mapcar '- a (mapcar '(lambda (x) (* x m)) b))
69.                   y0 (/ (caddr s) (cadr s))
70.                   x0 (/ (- (caddr a) (* (cadr a) y0)) (car a))
71.                   r0 (sqrt (+ (* x0 x0) (* y0 y0) (/ (+ (car sum) (caddr sum)) n)))
72.                   x0 (+ x0 (car avg))
73.                   y0 (+ y0 (cadr avg))
74.                    n (1- n)
75.          )
76.  
77.          (setq corr 9999999999
78.                tol 1e-06
79.              astep 0
80.            areport "Circular Non-Linear Regression Results \n"
81.          )
82.          (princ (strcat aReport  "\n"
83.                                  "  Initial: " (itoa aStep) "\n"
84.                                  "       Xo: " (rtos x0 2 8)  "\n"
85.                                  "       Yo: " (rtos y0 2 8)  "\n"
86.                                  "       R : " (rtos r0 2 8) "\n\n")
87.          )      
88.          (while (and (< astep 21) (> corr tol))
89.              (setq astep (1+ astep)
90.                    a00 0.0 a01 0.0 a02 0.0 a03 0.0
91.                    a10 0.0 a11 0.0 a12 0.0 a13 0.0
92.                    a20 0.0 a21 0.0 a22 0.0 a23 0.0
93.              )
94.              (foreach p pl
95.                   (setq x (car p) y (cadr p))
96.                   (setq bott (sqrt(+ (* x x) (* -2.0 x x0) (* y y)(* -2.0 y0 y) (* x0 x0)(* y0 y0)))
97.                           k1 (sqrt (+ (*(- x x0)(- x x0)) (* (- y y0) (- y y0))))
98.                            a (/ (- x0 x) bott)
99.                            b (/ (- y0 y) bott)
100.                            c -1.0
101.                            k (- 0.0 (- k1 r0))
102.                          a00 (+ a00 (* a a))  ; A*Sum(A**2)+B*2Sum(AB) +C*Sum(AC) +Sum(AK) =0
103.                          a01 (+ a01 (* a b))  ; A*Sum(AB) +B*2Sum(B**2)+C*Sum(BC) +Sum(BK) =0
104.                          a02 (+ a02 (* a c))  ; A*Sum(AC) +B*2Sum(BC) +C*Sum(C**2)+Sum(-K) =0
105.                          a03 (+ a03 (* a k))
106.                          a11 (+ a11 (* b b))
107.                          a12 (+ a12 (* b c))
108.                          a13 (+ a13 (* b k))
109.                          a22 (+ a22 (* c c))
110.                          a23 (+ a23 (* -1. k))
111.                   )
112.              )
113.              (setq a10 a01 a20 a02 a21 a12)
114.              (setq  MatA (list  (list a00 a01 a02 a03)
115.                                 (list a10 a11 a12 a13)
116.                                 (list a20 a21 a22 a23)
117.                          )       
118.                     MatR (Gauss MatA)
119.                  
120.                        c (/ (last (last MatR)) (caddr(last MatR)))
121.                        b (/ (- (last (cadr MatR)) (* (caddr (cadr MatR)) c)) (cadr(cadr MatR)))
122.                        a (/ (- (last (car  MatR)) (* (cadr (car MatR)) b) (* (caddr (car MatR)) c)) (car(car MatR)))
123.                     corr (min (max (abs a)(abs b)(abs c)) corr)
124.                       x0 (+ x0 a)
125.                       y0 (+ y0 b)
126.                       r0 (+ r0 c)
127.            
128.              )
129.              (setq aReport (strcat  "\n"
130.                               "Iteration: " (itoa aStep) "\n"
131.                               "       Xo: " (rtos x0 2 8)  "\n"
132.                               "       Yo: " (rtos y0 2 8)  "\n"
133.                               "       R : " (rtos r0 2 8) "\n\n"
134.                               "  Matrix : " (rtos a00 2 3) "  " (rtos a01 2 3) "  " (rtos a02 2 3) "  " (rtos a03 2 3) "\n"
135.                               "           " (rtos a10 2 3) "  " (rtos a11 2 3) "  " (rtos a12 2 3) "  " (rtos a13 2 3) "\n"
136.                               "           " (rtos a20 2 3) "  " (rtos a21 2 3) "  " (rtos a22 2 3) "  " (rtos a23 2 3) "\n\n"
137.                               "  Delta X: " (rtos a 2 8) "\n"
138.                               "  Delta Y: " (rtos b 2 8) "\n"
139.                               "  Delta R: " (rtos c 2 8) "\n\n"
140.                            )
141.              )
142.              (princ aReport)
143.              (setq aReport "")
144.          )
145.          (entmakex (list '(0 . "POINT") (cons 10 (list x0 y0))))
146.          (entmakex (list '(0 . "CIRCLE") (cons 10 (list x0 y0)) (cons 40 r0)))
147.          (foreach p pl
148.              (setq   ;p (cdr (assoc 10 (entget (ssname ss n))))
149.                    res (- r0 (distance p (list x0 y0)))
150.                      n (1- n)
151.              )      
152.              (entmakex (list '(0 . "TEXT") '(8 . "Residuals") (cons 10 p) (cons 40 3.0) '(50 . 1.5708) (cons 1 (rtos res 2 3))))
153.          )
154.       )
155.    )
156. )
157.  
158. ;*****************************************************************************;
159. ;; gauss                                                                      ;
160. ;; Implementation Gaussian Elimination     by ElpanovEvgeniy                  ;
161. ;; For text:  1x+2y+3z=2                                                      ;
162. ;;           10x+1y+8z=17                                                     ;
163. ;;               7z+2y=5                                                      ;
164. ;; (gauss '((1.0 2.0 3.0 2.0) (10.0 1.0 8.0 17.0) (0.0 2.0 7.0 5.0)))         ;
165. ;; =>                                                                         ;
166. ;; ((1.0 2.0 3.0 2.0) (0.0 -19.0 -22.0 -3.0) (0.0 0.0 4.68421 4.68421))       ;
167. ;;****************************************************************************;
168.  
169. (defun gauss (lst)
170.  
171.   (if (car lst)
172.     (if (zerop (caar lst))
173.       (if (vl-every (function zerop) (mapcar (function car) lst))
174.         (if (cdr lst)
175.           (cons
176.             (car lst)
177.             (mapcar
178.               (function (lambda (x) (cons 0. x)))
179.               (gauss (mapcar (function cdr) (cdr lst)))
180.             )
181.           )
182.           lst
183.         )
184.         (gauss
185.           (cons
186.             (mapcar
187.               (function +)
188.               (car lst)
189.               (car (vl-remove-if
190.                      (function (lambda (x) (zerop (car x))))
191.                      (cdr lst)
192.                    )
193.               )
194.             )
195.             (cdr lst)
196.           )
197.         )
198.       )
199.       (cons
200.         (car lst)
201.         (mapcar
202.           (function (lambda (x) (cons 0. x)))
203.           (gauss
204.             (mapcar
205.               (function
206.                 (lambda (x / i)
207.                   (setq i (/ (car x) (caar lst)))
208.                   (mapcar
209.                     (function -)
210.                     (cdr x)
211.                     (mapcar (function (lambda (x1) (* x1 i)))
212.                             (cdar lst)
213.                     )
214.                   )
215.                 )
216.               )
217.               (cdr lst)
218.             )
219.           )
220.         )
221.       )
222.     )
223.   )
224. )
225.  
226.  

[jarekw90]

Your Lisp is great, but I've one proposition to update.
Could you modify your Lisp to maked circle in the average elevation of selected points? It would be very usefull.