【越飞越高讲堂14】椭圆论

hhh454

真的是大师级的，受教了，万分感谢

扶风

我就喜欢挺这种贴，抱着学无止境的态度追逐

ynhh

向大师反应一个问题：
椭圆长短轴长为260X90这个尺寸，在双数等分周长时，CAD2010以下版本（如04，08）多数情况都会出错，但略改尺寸，如改为260X106等等或在单数等分时，则不出错。以上问题，在CAD2010版本中，260X90的椭圆在双数等分时又正常。
请教大师，看看是不是CAD版本的原因？

chlh_jd

感谢高飞大师的总结！
Lee Mac 大师写了个5点拟合椭圆的程序:http://www.lee-mac.com/5pointellipse.html
我写了个最小二乘法拟合椭圆源码，欢迎大家批评指正：
;;最小二乘法拟合椭圆

2. 
   
3. ;;;-------------------------------------------------------------
   
4. (defun LS-Pts->Ellipse (lst  /    n arg_mat   res_vec   a  b
   
5.     c    d    e f    res  xc   yc   aa  bb
   
6.     ang  co   si co2  si2)
   
7.   ;; by GSLS(SS)
   
8.   ;; Pointset return Ellipse by Least Square method
   
9.   ;; 2014-01-02
   
10.   (if (> (setq n (length lst)) 4)
    
11.     (progn
    
12.       (mapcar (function (lambda (p / x y)
    
13.      (setq x (car p)
    
14.     y (cadr p))
    
15.      (setq arg_mat (cons (list (* x y) (* y y) x y 1.)
    
16.            arg_mat)
    
17.     res_vec (cons (* -1. x x) res_vec))))
    
18.        lst)
    
19.       (setq res (solve-odeq arg_mat res_vec))
    
20.       (mapcar (function set) (quote (a b c d e f)) (cons 1. res))
    
21.       (if (= a c)
    
22. (setq ang _pi2)
    
23. (setq ang (* 0.5 (atan (/ b (- a c))))))
    
24.       (if (and (/= f 0) (/= (setq delt (- (* 4. c) (* b b))) 0))
    
25. (progn
    
26.    (setq xc  (/ (- (* b e) (* 2. c d)) delt)
    
27.   yc  (/ (- (* b d) (* 2. e)) delt)
    
28.   co  (cos ang)
    
29.   si  (sin ang)
    
30.   co2 (* co co)
    
31.   si2 (* si si)
    
32.   bb  (sqrt
    
33.         (abs (/ (- (* (- (* f co2)
    
34.            (* (+ (* xc xc)
    
35.           (* b xc yc)
    
36.           (* c yc yc))
    
37.        co2))
    
38.         (- (* 2 xc si2)
    
39.            (* 2 yc co si)
    
40.            (* (+ xc xc (* b yc)) si2)))
    
41.      (* (- (* f si2)
    
42.            (* (+ (* xc xc)
    
43.           (* b xc yc)
    
44.           (* c yc yc))
    
45.        si2))
    
46.         (- (* 2 xc co2)
    
47.            (* -2 yc co si)
    
48.            (* (+ xc xc (* b yc)) co2))))
    
49.          (+ (* 2 xc co2)
    
50.      (* 2 yc co si)
    
51.      (- (* (+ xc xc (* b yc)) co2))))))
    
52.   aa  (*
    
53.         (sqrt
    
54.    (abs (- (/ (+ (* 2 xc co2)
    
55.           (* 2 yc co si)
    
56.           (* -1 (+ xc xc (* b yc)) co2))
    
57.        (- (* 2 xc si2)
    
58.           (* 2 yc co si)
    
59.           (* (+ xc xc (* b yc)) si2))))))
    
60.         bb)
    
61.   )
    
62.    (list (list xc yc) (polar (list 0 0) ang aa) (/ bb aa)))))
    
63.     ))
    

<br /> ;;;-------------------------------------------------------------<br /> (defun LS-Pts->Ellipse (lst  /    n arg_mat   res_vec   a  b<br />     c    d    e f    res  xc   yc   aa  bb<br />     ang  co   si co2  si2)<br />   ;; by GSLS(SS)<br />   ;; Pointset return Ellipse by Least Square method <br />   ;; 2014-01-02<br />   (if (> (setq n (length lst)) 4)<br />     (progn<br />       (mapcar (function (lambda (p / x y)<br />      (setq x (car p)<br />     y (cadr p))<br />      (setq arg_mat (cons (list (* x y) (* y y) x y 1.)<br />            arg_mat)<br />     res_vec (cons (* -1. x x) res_vec))))<br />        lst)<br />       (setq res (solve-odeq arg_mat res_vec))<br />       (mapcar (function set) (quote (a b c d e f)) (cons 1. res))<br />       (if (= a c)<br /> (setq ang _pi2)<br /> (setq ang (* 0.5 (atan (/ b (- a c))))))<br />       (if (and (/= f 0) (/= (setq delt (- (* 4. c) (* b b))) 0))<br /> (progn<br />    (setq xc  (/ (- (* b e) (* 2. c d)) delt)<br />   yc  (/ (- (* b d) (* 2. e)) delt)<br />   co  (cos ang)<br />   si  (sin ang)<br />   co2 (* co co)<br />   si2 (* si si)<br />   bb  (sqrt<br />         (abs (/ (- (* (- (* f co2)<br />            (* (+ (* xc xc)<br />           (* b xc yc)<br />           (* c yc yc))<br />        co2))<br />         (- (* 2 xc si2)<br />            (* 2 yc co si)<br />            (* (+ xc xc (* b yc)) si2)))<br />      (* (- (* f si2)<br />            (* (+ (* xc xc)<br />           (* b xc yc)<br />           (* c yc yc))<br />        si2))<br />         (- (* 2 xc co2)<br />            (* -2 yc co si)<br />            (* (+ xc xc (* b yc)) co2))))<br />          (+ (* 2 xc co2)<br />      (* 2 yc co si)<br />      (- (* (+ xc xc (* b yc)) co2))))))<br />   aa  (*<br />         (sqrt<br />    (abs (- (/ (+ (* 2 xc co2)<br />           (* 2 yc co si)<br />           (* -1 (+ xc xc (* b yc)) co2))<br />        (- (* 2 xc si2)<br />           (* 2 yc co si)<br />           (* (+ xc xc (* b yc)) si2))))))<br />         bb)<br />   )<br />    (list (list xc yc) (polar (list 0 0) ang aa) (/ bb aa)))))<br />     ))<br />

通用最小二乘法函数

2. 
   
3. (defun solve-odeq  (arg_mat res_vector / mat)
   
4.   ;; Universal least squares method
   
5.   ;; arg_mat    -- Parameters Matrix of Overdetermined Equations
   
6.   ;; res_vector -- Right term of Overdetermined Equations
   
7.   ;; by GSLS(SS) 2013.11.22
   
8.   (setq mat ([*] ([trp] arg_mat) arg_mat)
   
9. mat ([inv] mat))
   
10.   (mxv ([*] mat ([trp] arg_mat)) res_vector))
    
11. ;; Use function ---------------------------------------------
    
12. ;; DotProduct
    
13. (defun vxv (v1 v2) (apply (function +) (mapcar (function *) v1 v2)))
    
14. ;; Matrix Vector Multiply
    
15. (defun mxv (m v) (mapcar (function (lambda (r) (vxv r v))) m))
    
16. ;; Transpose matrix
    
17. (defun [trp] (a) (apply (function mapcar) (cons (function list) a)))
    
18. ;; Matrix is multiplied by a matrix
    
19. (defun [*] (m q) (mapcar (function (lambda (r) (mxv ([trp] q) r))) m))
    
20. ;; To get Unit diagonal matrix
    
21. (defun [I]  (d / i r n m)
    
22.   (setq i d)
    
23.   (while (<= 0 (setq i (1- i)))
    
24.     (setq n d
    
25.    r nil)
    
26.     (while (<= 0 (setq n (1- n)))
    
27.       (setq r (cons (if (= i n)
    
28.         1.0
    
29.         0.0)
    
30.       r)))
    
31.     (setq m (cons r m))))
    
32. ;;; gile-inverse-matrix (gile) 2009/03/17
    
33. ;; uses the gauss-jordan elimination method to calculate the inverse matrix of any dimension square matrix
    
34. ;; argument : a square matrix
    
35. ;; return : the inverse matrix (or nil if singular)
    
36. ;; ([inv] '((1 2 3) (2 4 5) (3 5 6)));_([inv] '((2 1) (5 3)))
    
37. (defun [inv]  (mat / col piv row res)
    
38.   (setq mat (mapcar (function (lambda (x1 x2)
    
39.     (append x1 x2)))
    
40.       mat
    
41.       ([I] (length mat))))
    
42.   (while mat
    
43.     (setq col (mapcar (function (lambda (x) (abs (car x)))) mat))
    
44.     (repeat (vl-position (apply (function max) col) col)
    
45.       (setq mat (append (cdr mat) (list (car mat)))))
    
46.     (if (equal (setq piv (caar mat)) 0.0 1e-14)
    
47.       (setq mat nil
    
48.      res nil)
    
49.       (setq piv (/ 1.0 (caar mat))
    
50.      row (mapcar (function (lambda (x) (* x piv))) (car mat))
    
51.      mat (mapcar (function
    
52.      (lambda (r / e)
    
53.        (setq e (car r))
    
54.        (cdr (mapcar (function (lambda (x n)
    
55.            (- x (* n e))))
    
56.       r
    
57.       row))))
    
58.    (cdr mat))
    
59.      res (cons (cdr row)
    
60.         (mapcar
    
61.    (function
    
62.      (lambda (r / e)
    
63.        (setq e (car r))
    
64.        (cdr (mapcar (function (lambda (x n)
    
65.            (- x (* n e))))
    
66.       r
    
67.       row))))
    
68.    res)))))
    
69.   (reverse res))
    

<br /> (defun solve-odeq  (arg_mat res_vector / mat)<br />   ;; Universal least squares method<br />   ;; arg_mat    -- Parameters Matrix of Overdetermined Equations<br />   ;; res_vector -- Right term of Overdetermined Equations<br />   ;; by GSLS(SS) 2013.11.22<br />   (setq mat ([*] ([trp] arg_mat) arg_mat)<br /> mat ([inv] mat))<br />   (mxv ([*] mat ([trp] arg_mat)) res_vector))<br /> ;; Use function ---------------------------------------------<br /> ;; DotProduct<br /> (defun vxv (v1 v2) (apply (function +) (mapcar (function *) v1 v2)))<br /> ;; Matrix Vector Multiply<br /> (defun mxv (m v) (mapcar (function (lambda (r) (vxv r v))) m))<br /> ;; Transpose matrix<br /> (defun [trp] (a) (apply (function mapcar) (cons (function list) a)))<br /> ;; Matrix is multiplied by a matrix<br /> (defun [*] (m q) (mapcar (function (lambda (r) (mxv ([trp] q) r))) m))<br /> ;; To get Unit diagonal matrix<br /> (defun [I]  (d / i r n m)<br />   (setq i d)<br />   (while (<= 0 (setq i (1- i)))<br />     (setq n d<br />    r nil)<br />     (while (<= 0 (setq n (1- n)))<br />       (setq r (cons (if (= i n)<br />         1.0<br />         0.0)<br />       r)))<br />     (setq m (cons r m))))<br /> ;;; gile-inverse-matrix (gile) 2009/03/17<br /> ;; uses the gauss-jordan elimination method to calculate the inverse matrix of any dimension square matrix<br /> ;; argument : a square matrix<br /> ;; return : the inverse matrix (or nil if singular)<br /> ;; ([inv] '((1 2 3) (2 4 5) (3 5 6)));_([inv] '((2 1) (5 3)))<br /> (defun [inv]  (mat / col piv row res)<br />   (setq mat (mapcar (function (lambda (x1 x2)<br />     (append x1 x2)))<br />       mat<br />       ([I] (length mat))))<br />   (while mat<br />     (setq col (mapcar (function (lambda (x) (abs (car x)))) mat))<br />     (repeat (vl-position (apply (function max) col) col)<br />       (setq mat (append (cdr mat) (list (car mat)))))<br />     (if (equal (setq piv (caar mat)) 0.0 1e-14)<br />       (setq mat nil<br />      res nil)<br />       (setq piv (/ 1.0 (caar mat))<br />      row (mapcar (function (lambda (x) (* x piv))) (car mat))<br />      mat (mapcar (function<br />      (lambda (r / e)<br />        (setq e (car r))<br />        (cdr (mapcar (function (lambda (x n)<br />            (- x (* n e))))<br />       r<br />       row))))<br />    (cdr mat))<br />      res (cons (cdr row)<br />         (mapcar<br />    (function<br />      (lambda (r / e)<br />        (setq e (car r))<br />        (cdr (mapcar (function (lambda (x n)<br />            (- x (* n e))))<br />       r<br />       row))))<br />    res)))))<br />   (reverse res))<br />

测试函数

2. 
   
3. ;; E.G.
   
4. ;; For Ellipse
   
5. (defun c:test  (/ l r p10 a b p11 d40)
   
6.   (prompt
   
7.     "\nFit Ellipse though Least-square method , the points you selected must be at least 5 !!!")
   
8.   (setq l (mapcar (function (lambda (x) (cdr (assoc 10 (entget x)))))
   
9.     (vl-remove-if-not
   
10.       (function (lambda (x) (eq (type x) 'ENAME)))
    
11.       (mapcar 'cadr (ssnamex (ssget '((0 . "POINT")))))))
    
12. l (mapcar (function (lambda (p)
    
13.          (trans p 0 1 t)))
    
14.     l))
    
15.   (if (and (> (length l) 4) (setq r (LS-Pts->Ellipse l)))
    
16.     (progn
    
17.       (setq p10 (car r)
    
18.      p11 (cadr r)
    
19.      d40 (caddr r))
    
20.       (entmake
    
21. (append
    
22.    '(
    
23.      (000 . "ELLIPSE")
    
24.      (100 . "AcDbEntity")
    
25.      (100 . "AcDbEllipse"))
    
26.    (list (cons 10 (trans p10 1 0 t))
    
27.   (cons 11 (trans p11 1 0))
    
28.   (cons 40 d40))
    
29.    (list (cons 210 (trans '(0.0 0.0 1.0) 1 0 t)))))))
    
30.   (princ)
    
31.   )
    

<br /> ;; E.G.<br /> ;; For Ellipse<br /> (defun c:test  (/ l r p10 a b p11 d40)<br />   (prompt<br />     "\nFit Ellipse though Least-square method , the points you selected must be at least 5 !!!")<br />   (setq l (mapcar (function (lambda (x) (cdr (assoc 10 (entget x)))))<br />     (vl-remove-if-not<br />       (function (lambda (x) (eq (type x) 'ENAME)))<br />       (mapcar 'cadr (ssnamex (ssget '((0 . "POINT")))))))<br /> l (mapcar (function (lambda (p)<br />          (trans p 0 1 t)))<br />     l))<br />   (if (and (> (length l) 4) (setq r (LS-Pts->Ellipse l)))<br />     (progn<br />       (setq p10 (car r)<br />      p11 (cadr r)<br />      d40 (caddr r))<br />       (entmake<br /> (append<br />    '(<br />      (000 . "ELLIPSE")<br />      (100 . "AcDbEntity")<br />      (100 . "AcDbEllipse"))<br />    (list (cons 10 (trans p10 1 0 t))<br />   (cons 11 (trans p11 1 0))<br />   (cons 40 d40))<br />    (list (cons 210 (trans '(0.0 0.0 1.0) 1 0 t)))))))<br />   (princ)<br />   )<br />

看天的小树

楼主的胸怀让人敬佩

zzcollan

大师啊，真正的大师

zst1978

highflybir  是中国公开研究和发表很多椭圆相关高精尖知识的大师

tanjurun

highflybir  的胸怀，所公开的内容没有留一手，并很热心帮助明经上的朋友，热心免费提供源代码
向 highflybir  致敬

tanjurun

真的是大师级的，受教了，万分感谢

tomonkey239

tanjurun 发表于 2021-1-19 08:53
highflybir  的胸怀，所公开的内容没有留一手，并很热心帮助明经上的朋友，热心免费提供源代码
向 highfly ...

的确，值得钦佩。。。。