« previous next »
Pages: [1] |   Go Down

Author Topic: -={Challenge}=-Universal least squares method  (Read 2608 times)

0 Members and 1 Guest are viewing this topic.

chlh_jd

  • Guest
-={Challenge}=-Universal least squares method
« on: February 26, 2014, 01:14:54 AM »
Hi all , here's challenge to let gurus participation  :-)
The topic as the title .
Here's my code
Code - Auto/Visual Lisp: [Select]
  1. (defun solve-odeq (arg_mat res_vector / mat)
  2.   ;; Universal least squares method
  3.   ;; arg_mat    -- Parameters Matrix of Overdetermined Equations
  4.   ;; res_vector -- Right term of Overdetermined Equations
  5.   ;; by GSLS(SS) 2013.11.22
  6.   (setq mat ([*] ([trp] arg_mat) arg_mat)
  7.         mat ([inv] mat))
  8.   (mxv ([*] mat ([trp] arg_mat)) res_vector))
  9.  
  10. ;; Use function ---------------------------------------------
  11. ;; DotProduct
  12. (defun vxv(v1 v2)(apply(function +)(mapcar(function *)v1 v2)))
  13. ;; Matrix Vector Multiply
  14. (defun mxv(m v)(mapcar(function(lambda (r)(vxv r v)))m))
  15. ;; Transpose matrix
  17. ;; Matrix is multiplied by a matrix
  18. (defun [*](m q)(mapcar(function (lambda(r)(mxv ([trp] q)r)))m))
  19. ;; To get Unit diagonal matrix
  20. (defun [I] (d / i r n m)
  21.   (setq i d)
  22.   (while (<= 0 (setq i (1- i)))
  23.     (setq n d r nil)
  24.     (while (<= 0(setq n(1- n)))
  25.       (setq r (cons(if(= i n)1.0 0.0)r)))
  26.     (setq m(cons r m))))
  27. ;;; gile-inverse-matrix (gile) 2009/03/17
  28. ;; uses the gauss-jordan elimination method to calculate the inverse matrix of any dimension square matrix
  29. ;; argument : a square matrix
  30. ;; return : the inverse matrix (or nil if singular)
  31. ;; ([inv] '((1 2 3) (2 4 5) (3 5 6)));_([inv] '((2 1) (5 3)))
  32. (defun [inv] (mat / col piv row res)
  33.   (setq mat (mapcar (function (lambda (x1 x2)
  34.                  (append  x1  x2))) mat([I] (length mat))))
  35.   (while mat
  36.     (setq col (mapcar(function(lambda (x)(abs (car x))))mat))
  38.       (setq mat (append (cdr mat) (list (car mat)))))  
  39.     (if (equal (setq piv (caar mat))0.0 1e-14)
  40.       (setq mat nil  res nil)      
  41.       (setq piv (/ 1.0 (caar mat))
  42.             row (mapcar(function(lambda (x)(* x piv)))(car mat))
  43.             mat (mapcar (function(lambda (r / e)
  44.                      (setq e (car r))(cdr (mapcar(function(lambda (x n)
  45.                                (- x (* n e))))r row))))          
  46.                   (cdr mat))
  47.             res (cons (cdr row)
  48.                       (mapcar(function(lambda (r / e)
  49.                            (setq e (car r))(cdr(mapcar(function(lambda (x n)
  50.                                      (- x (* n e))))r row))))res)))))
  51.   (reverse res))
  52.  

EDIT:KDUB -> CODE FORMAT code=cadlisp-7
« Last Edit: February 26, 2014, 01:32:46 AM by chlh_jd »
Logged

chlh_jd

  • Guest
Re: -={Challenge}=-Universal least squares method
« Reply #1 on: February 26, 2014, 01:31:00 AM »
E.G.
Thanks Lee Mac's program : Five point ellipse . http://www.lee-mac.com/5pointellipse.html  It's really a very concise function .
Here I use least-square method to solve it .
Code - Auto/Visual Lisp: [Select]
  1. (defun LS-Pts->Ellipse  (lst  /    n    arg_mat   res_vec   a    b
  2.                          c    d    e    f    res  xc   yc   aa   bb
  3.                          ang  co   si   co2  si2)
  4.   ;; by GSLS(SS)
  5.   ;; Pointset return Ellipse by Least Square method
  6.   ;; 2014-01-02
  7.   (if (> (setq n (length lst)) 4)
  8.     (progn
  9.       (mapcar (function (lambda (p / x y)
  10.                           (setq x (car p)
  11.                                 y (cadr p))
  12.                           (setq arg_mat (cons (list (* x y) (* y y) x y 1.)
  13.                                               arg_mat)
  14.                                 res_vec (cons (* -1. x x) res_vec))))
  15.               lst)
  16.       (setq res (solve-odeq arg_mat res_vec))
  17.       (mapcar (function set) (quote (a b c d e f)) (cons 1. res))
  18.       (if (= a c)
  19.         (setq ang _pi2)
  20.         (setq ang (* 0.5 (atan (/ b (- a c))))))
  21.       (if (and (/= f 0) (/= (setq delt (- (* 4. c) (* b b))) 0))
  22.         (progn
  23.           (setq xc  (/ (- (* b e) (* 2. c d)) delt)
  24.                 yc  (/ (- (* b d) (* 2. e)) delt)
  25.                 co  (cos ang)
  26.                 si  (sin ang)
  27.                 co2 (* co co)
  28.                 si2 (* si si)
  29.                 bb  (sqrt
  30.                       (abs (/ (- (* (- (* f co2)
  31.                                        (* (+ (* xc xc)
  32.                                              (* b xc yc)
  33.                                              (* c yc yc))
  34.                                           co2))
  35.                                     (- (* 2 xc si2)
  36.                                        (* 2 yc co si)
  37.                                        (* (+ xc xc (* b yc)) si2)))
  38.                                  (* (- (* f si2)
  39.                                        (* (+ (* xc xc)
  40.                                              (* b xc yc)
  41.                                              (* c yc yc))
  42.                                           si2))
  43.                                     (- (* 2 xc co2)
  44.                                        (* -2 yc co si)
  45.                                        (* (+ xc xc (* b yc)) co2))))
  46.                               (+ (* 2 xc co2)
  47.                                  (* 2 yc co si)
  48.                                  (- (* (+ xc xc (* b yc)) co2))))))
  49.                 aa  (*
  50.                       (sqrt
  51.                         (abs (- (/ (+ (* 2 xc co2)
  52.                                       (* 2 yc co si)
  53.                                       (* -1 (+ xc xc (* b yc)) co2))
  54.                                    (- (* 2 xc si2)
  55.                                       (* 2 yc co si)
  56.                                       (* (+ xc xc (* b yc)) si2))))))
  57.                       bb)
  58.                 )
  59.           (list (list xc yc) (polar (list 0 0) ang aa) (/ bb aa)))))
  60.     ))
  61.  
Test
Code - Auto/Visual Lisp: [Select]
  1. (defun c:test  (/ l r p10 a b p11 d40)
  2.   (prompt
  3.     "\nFit Ellipse though Least-square method , the points you selected must be at least 5 !!!")
  4.   (setq l (mapcar (function (lambda (x) (cdr (assoc 10 (entget x)))))
  5.                   (vl-remove-if-not
  6.                     (function (lambda (x) (eq (type x) 'ENAME)))
  7.                     (mapcar 'cadr (ssnamex (ssget '((0 . "POINT")))))))
  8.         l (mapcar (function (lambda (p)
  9.                               (trans p 0 1 t)))
  10.                   l))
  11.   (if (and (> (length l) 4) (setq r (LS-Pts->Ellipse l)))
  12.     (progn
  13.       (setq p10 (car r)
  14.             p11 (cadr r)
  15.             d40 (caddr r))
  16.       (entmake
  17.         (append
  18.           '(
  19.             (000 . "ELLIPSE")
  20.             (100 . "AcDbEntity")
  21.             (100 . "AcDbEllipse"))
  22.           (list (cons 10 (trans p10 1 0 t))
  23.                 (cons 11 (trans p11 1 0))
  24.                 (cons 40 d40))
  25.           (list (cons 210 (trans '(0.0 0.0 1.0) 1 0 t)))))))
  26.   (princ)
  27.   )
  28.  
« Last Edit: March 08, 2014, 10:50:35 AM by chlh_jd »
Logged

chlh_jd

  • Guest
Re: -={Challenge}=-Universal least squares method
« Reply #2 on: February 26, 2014, 01:34:47 AM »
Thanks KDUB for Edit  :-)
Logged

chlh_jd

  • Guest
Re: -={Challenge}=-Universal least squares method
« Reply #3 on: February 26, 2014, 01:38:36 AM »
Fit Circle though Least-square method .
Code - Auto/Visual Lisp: [Select]
  1. (defun LS-Pts->Circle  (lst   nchsep      /     arg_mat     res_vec
  2.                         ps    pe    xs    ys    xe    ye    xc    yc
  3.                         a     b     c     d     res)
  4.   ;; Function:
  5.   ;; Pointset return Circle by Least Square method
  6.   ;; args :
  7.   ;;     lst --  Pointset , > 3 points
  8.   ;;     nchsep -- Is Not change the start and end point ? T or Nil
  9.   ;; by GSLS(SS)
  10.   ;; 2013-11-16
  11.   ;; 2014-01-01 Edited , Add constant the startpoint and endpoint mode .
  12.   (if (> (length lst) 2)
  13.     (if nchsep
  14.       (progn
  15.         (setq ps (car lst)
  16.               xs (car ps)
  17.               ys (cadr ps)
  18.               pe (last lst)
  19.               xe (car pe)
  20.               ye (cadr pe))
  21.         (mapcar (function (lambda (p / x y)
  22.                             (setq x (car p)
  23.                                   y (cadr p))
  24.                             (setq arg_mat (cons (list (- (+ x x) xs xe)
  25.                                                       (- (+ y y) ys ye))
  26.                                                 arg_mat)
  27.                                   res_vec (cons
  28.                                             (- (+ (* x x) (* y y))
  29.                                                (/2 (+ (* xs xs)
  30.                                                       (* ys ys)
  31.                                                       (* xe xe)
  32.                                                       (* ye ye))))
  33.                                             res_vec))))
  34.                 (butlast (cdr lst)))
  35.         (setq res (solve-odeq arg_mat res_vec))
  36.         (if (and (cadr res)
  37.                  (setq xc (car res)
  38.                   yc (cadr res))
  39.                  (setq d (/2    (+ (distance (list xc yc) ps)
  40.                            (distance (list xc yc) pe))))
  41.                  (setq c
  42.                    (car
  43.                      (vl-sort
  44.                        (c_int_c ps d pe d)
  45.                        (function (lambda (a b)
  46.                                    (< (distance a (list xc yc))
  47.                                       (distance b (list xc yc)))))))))
  48.           (list c d)))
  49.       (progn
  50.         (mapcar (function (lambda (p / x y)
  51.                             (setq x (car p)
  52.                                   y (cadr p))
  53.                             (setq arg_mat (cons (list x y 1.) arg_mat)
  54.                                   res_vec (cons (/2 (+ (* x x) (* y y))) res_vec))))
  55.                 lst)
  56.         (setq res (solve-odeq arg_mat res_vec))
  57.         (if (and (setq c (caddr res))
  58.                  (setq a (car res)
  59.                        b (cadr res))
  60.                  (> (setq d (+ (* a a) (* b b) c c)) 0))
  61.           (list (list a b) (sqrt d)))))))
  62.  
Test
Code: [Select]
(defun c:Test  (/ l r p10 d40)
  (prompt "\nFit Circle though Least-square method , the points you selected must be at least 3 !!!")
  (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) 2)(setq r (LS-Pts->Circle l nil)))    
    (progn
      (setq p10 (car r)    
    d40 (cadr r))
      (entmakex
(list (cons 0 "CIRCLE")
      (cons 10 (trans p10 1 0 t))
      (cons 40 d40)
      (cons 210 (trans (quote(0.0 0.0 1.0)) 1 0 t))))
)) 
  (princ)
  )
Logged

chlh_jd

  • Guest
Re: -={Challenge}=-Universal least squares method
« Reply #4 on: February 26, 2014, 01:45:38 AM »
Fit Line though least-square method .
Code - Auto/Visual Lisp: [Select]
  1. (defun LS-Pts->Line  (lst   /    n  arg_mat res_vec a b c d  res)
  2.   ;; Pointset return Line by Least Square method
  3.   ;; lst -- 2D Pointset
  4.   ;; kx+b=Y
  5.   ;; by GSLS(SS) 2013-11-16
  6.   (if (> (setq n (length lst)) 1)
  7.     (progn
  8.       (mapcar (function (lambda (p / x y)
  9.                           (setq x (car p)
  10.                                 y (cadr p))
  11.                           (setq arg_mat (cons (list x 1.) arg_mat)
  12.                                 res_vec (cons y res_vec)))) lst)
  13.       (solve-odeq arg_mat res_vec)            
  14.       ))
  15.   )
  16.  
Test
Code: [Select]
(defun c:TEST  (/ l r k b xa ya xb yb pa pb p10 p11)
  (prompt "\nFit Line though Least-square method，the points you selected must be at least 2 !!!")
  (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) 1)(setq r (LS-Pts->Line l)))    
    (progn
      (setq k (car r)
    b (cadr r))
      (setq l (vl-sort l (function (lambda (a b) (< (car a) (car b))))))
      (setq xa (caar l)
    ya (+ (* k xa) b)
    xb (car (last l))
    yb (+ (* k xb) b)
    pa (list xa ya)
    pb (list xb yb)
    p10 (pedal (car l) pa pb)
    p11 (pedal (last l) pa pb))     
      (setq p10 (trans p10 1 0 t)
    p11 (trans p11 1 0 t))
      (entmakex
(list (cons 0 "LINE")
      (cons 10 p10)
      (cons 11 p11)
      (cons 210 (trans (quote(0.0 0.0 1.0)) 1 0 t))))
)) 
  (princ)
  )
« Last Edit: March 08, 2014, 10:52:34 AM by chlh_jd »
Logged

dgpuertas

  • Newt
  • Posts: 81
Re: -={Challenge}=-Universal least squares method
« Reply #5 on: February 26, 2014, 11:57:25 AM »
Logged

chlh_jd

  • Guest
Re: -={Challenge}=-Universal least squares method
« Reply #6 on: February 27, 2014, 03:46:31 AM »
Thanks dgpuertas a lot ! :-)
Logged
Pages: [1] |   Go Up
« previous next »