解一元二次方程lisp_晓东CAD家园-论坛-算法与数据结构-Lisp通用最小二乘法-通用最小二乘法在代码和原理上相对比较简单，麻烦的是每次应用都得把高次元参数转换成低次参数，这里写出我的写法和大家探...

;;;-------------------------------------------------------------

(defun LS-Pts->Ellipse (lst  /    n arg_mat   res_vec   a  b

c    d    e f    res  xc   yc   aa  bb

ang  co   si co2  si2)

;; by GSLS(SS)

;; Pointset return Ellipse by Least Square method

;; 2014-01-02

(if (> (setq n (length lst)) 4)

(progn

(mapcar (function (lambda (p / x y)

(setq x (car p)

y (cadr p))

(setq arg_mat (cons (list (* x y) (* y y) x y 1.)

arg_mat)

res_vec (cons (* -1. x x) res_vec))))

lst)

(setq res (solve-odeq arg_mat res_vec))

(mapcar (function set) (quote (a b c d e f)) (cons 1. res))

(if (= a c)

(setq ang _pi2)

(setq ang (* 0.5 (atan (/ b (- a c))))))

(if (and (/= f 0) (/= (setq delt (- (* 4. c) (* b b))) 0))

(progn

(setq xc  (/ (- (* b e) (* 2. c d)) delt)

yc  (/ (- (* b d) (* 2. e)) delt)

co  (cos ang)

si  (sin ang)

co2 (* co co)

si2 (* si si)

bb  (sqrt

(abs (/ (- (* (- (* f co2)

(* (+ (* xc xc)

(* b xc yc)

(* c yc yc))

co2))

(- (* 2 xc si2)

(* 2 yc co si)

(* (+ xc xc (* b yc)) si2)))

(* (- (* f si2)

(* (+ (* xc xc)

(* b xc yc)

(* c yc yc))

si2))

(- (* 2 xc co2)

(* -2 yc co si)

(* (+ xc xc (* b yc)) co2))))

(+ (* 2 xc co2)

(* 2 yc co si)

(- (* (+ xc xc (* b yc)) co2))))))

aa  (*

(sqrt

(abs (- (/ (+ (* 2 xc co2)

(* 2 yc co si)

(* -1 (+ xc xc (* b yc)) co2))

(- (* 2 xc si2)

(* 2 yc co si)

(* (+ xc xc (* b yc)) si2))))))

bb)

)

(list (list xc yc) (polar (list 0 0) ang aa) (/ bb aa)))))

))

;;;-------------------------------------------------------------

;; E.G.

;; For Ellipse

(defun c:test  (/ l r p10 a b p11 d40)

(prompt

"\nFit Ellipse though Least-square method , the points you selected must be at least 5 !!!")

(setq l (mapcar (function (lambda (x) (cdr (assoc 10 (entget x)))))

(vl-remove-if-not

(function (lambda (x) (eq (type x) 'ENAME)))

(mapcar 'cadr (ssnamex (ssget '((0 . "POINT")))))))

l (mapcar (function (lambda (p)

(trans p 0 1 t)))

l))

(if (and (> (length l) 4) (setq r (LS-Pts->Ellipse l)))

(progn

(setq p10 (car r)

p11 (cadr r)

d40 (caddr r))

(entmake

(append

'(

(000 . "ELLIPSE")

(100 . "AcDbEntity")

(100 . "AcDbEllipse"))

(list (cons 10 (trans p10 1 0 t))

(cons 11 (trans p11 1 0))

(cons 40 d40))

(list (cons 210 (trans '(0.0 0.0 1.0) 1 0 t)))))))

(princ)

)