Build an average line from a family of points

[MSasu]

I have a family of points that resemble a straight line – those points were coming from a third-part application and, due to its low precision, they doesn’t lay on a line, even if they should – please refer to picture below. I’m trying right now to design an algorithm to rebuild a straight line from those points.

 

 

I will appreciate any suggestions/solutions on this matter.

Thank you very much for your attention

[Organic]

Linear regression could probably do it although it might not be easy to program.

[Lee Mac]

Linear regression could probably do it although it might not be easy to program.

 

Agreed, look into the Least Squares Method.

[MSasu]

Thank you for your suggestions, Gentlemen. I will investigate the math behind those theorems to see if I'm able to implement them in AutoLISP.

[Lee Mac]

Since you only have one dependent and one independent variable, you can use Simple Linear Regression...

 

(defun _lineartrendline ( l / a ax ay b f x y )
   (setq f '(( l ) (/ (apply '+ l) (length l)))
         x  (mapcar 'car  l)
         y  (mapcar 'cadr l)
         ax (f x)
         ay (f y)
         b  (/ (- (f (mapcar '* x y)) (* ax ay))
               (- (f (mapcar '* x x)) (* ax ax))
            )
         a  (- ay (* b ax))
   )
   (entmake
       (list
          '(000 . "RAY")
          '(100 . "AcDbEntity")
          '(100 . "AcDbRay")
           (list 10 0.0 a 0.0)
           (list 11 1.0 b 0.0)
       )
   )
)

Test code:

 

(defun c:test ( / i l s )
   (if (setq s (ssget '((0 . "POINT") (-4 . "*,*,=") (10 0.0 0.0 0.0))))
       (progn
           (repeat (setq i (sslength s))
               (setq l (cons (cdr (assoc 10 (entget (ssname s (setq i (1- i)))))) l))
           )
           (_lineartrendline l)
       )
   )
   (princ)
)

[MSasu]

Lee, I have to say, that is really EXTRAORDINARY! Thank you so much for your support!

 

 

 

 

PS. That is an interesting way to define a lambda function; this is the first time I see it.

[Lee Mac]

Lee, I have to say, that is really EXTRAORDINARY! Thank you so much for your support!

 

Many thanks Mircea!

 

PS. That is an interesting way to define a lambda function; this is the first time I see it.

 

I'm glad you like it, I sometimes use this construct for extra concision; though, note that it is not strictly a lambda function, but rather the literal equivalent of a defun-q expression, e.g.:

 

_$ (defun-q f ( x ) (* x x))
F
_$ f
((X) (* X X))
_$ (setq g '(( x ) (* x x)))
((X) (* X X))
_$ (equal f g)
T

[motee-z]

Lee Mac can be on arc or circle

[Stefan BMR]

Salut Mircea

 

Linear regression assure that the sum of differences between a series of values (e.g. y) and the approximation line is minimum.

It is working good if data meet some criteria.

For example, set your points around a vertical line, then use simple linear regression (like Lee's lisp. BTW Lee, your code is amazing).

 

Here is my lisp, using a formula from another (surprisingly) domain.

 

(defun C:TEST ( / vxv ss i n l o d dx dy a)
 
 (defun vxv (a b) (apply '+ (mapcar '* a b)))
 
 (if
   (setq ss (ssget '((0 . "POINT"))))
   (progn
     (repeat (setq i (sslength ss) n i)
       (setq l (cons (cdr (assoc 10 (entget (ssname ss (setq i (1- i)))))) l))
       )
     (setq o  (mapcar '/ (apply 'mapcar (cons '+ l)) (list n n))
           d  (mapcar '(lambda (a) (mapcar '- o a)) l)
           dx (mapcar 'car d)
           dy (mapcar 'cadr d)
           a  (* 0.5 (atan (* 2 (vxv dx dy)) (- (vxv dx dx) (vxv dy dy))))
           )
     (entmake
       (list
         '(0 . "XLINE")
         '(100 . "AcDbEntity")
         '(100 . "AcDbXline")
         '(62 . 1)
         (cons 10 o)
         (list 11 (cos a) (sin a) 0.0)
         )
       )
     )
   )
 (princ)
 )

Edited January 18, 2013 by Stefan BMR

[MSasu]

Multumesc mult Stefan pentru solutie! Thank you a lot Stefan for your solution.

 

Although the results of the two routines are slightly different on a scattered set of points, both give excellent result for the data that I have to deal with.

Thank you again Gentlemen for your valuable support.

[Lee Mac]

Linear regression assure that the sum of differences between a series of values (e.g. y) and the approximation line is minimum.

It is working good if data meet some criteria.

For example, set your points around a vertical line, then use simple linear regression (like Lee's lisp. BTW Lee, your code is amazing).

 

Here is my lisp, using a formula from another (surprisingly) domain.

 

Thank you Stefan!

Your code produces good results.

[Tyke]

@ Lee

 

Does your code work for 3D points or just for 2D points? I tried it out on a batch of 3D points and it produced some queer results, whilst that from Stefan seemed to work.

 

After further checking, the code from Stefan produces an XLINE through the 3D points but it is only a 2D entity.

 

Is it possible to create a line as apposed to an XLINE or a RAY just through the cluster of points?

Edited January 23, 2013 by Tyke
Update

[Lee Mac]

Does your code work for 3D points or just for 2D points? I tried it out on a batch of 3D points and it produced some queer results, whilst that from Stefan seemed to work.

 

Hi Tyke,

 

Since my code uses Simple Linear Regression, it will only apply to 2-variable systems, i.e. 2D points; Linear Regression of higher dimensions is certainly possible, however, somewhat more complex.

[Tyke]

Since my code uses Simple Linear Regression, it will only apply to 2-variable systems, i.e. 2D points; Linear Regression of higher dimensions is certainly possible, however, somewhat more complex.

 

Thanks Lee, I thought that might be the case. I have subsets of points from a 3D point cloud, approx 200 points, and need to get the best 3D fit line for the points. Any ideas how I could achieve that?

[Tyke]

My LISP is very rusty indeed so I have implimented it in VBA and have drawn a line just in the cluster of points. I intend now to do a further linear regression but using the X and Z coordinates of the points to get the gradient of the line through the points. The maths required for a 3D linear regression is beyond my capabilities.

[ymg3]

In the context of the original question I believe Stefan's solution is superior as it

does apply an orthogonal regression to the line and minimize the distance of the points to the line.

 

Lee's solution (Outstanding implementation) is a different beast and minize the distance along the Y ordinate to the line and is totally correct when the X is considered an independant variable.

 

Solution to this problem was first published in The Analyst, Vol. 5, No. 2 (Mar., 1878 ), pp. 53-54Published by R. J. ADCOCK.

Here's a link to this paper : A PROBLEM IN LEAST SQUARE

 

It is also known as Deming's regression or Orthogonal's Regression.

 

ymg

 

Here an image from Wikipedia showing Stefan's solution:

 

And here is what Lee is minimizing:

 

Edited November 6, 2013 by ymg3

[Lee Mac]

Very interesting & informative comparison ymg - thank you for the information!

[ymg3]

Tyke,

 

Here is your Orthogonal Regression in 3D.

 

It might be too late but was busy on something else.

 

Code is OK, but I still need to add a table of residuals for analysis of results.

 

For derivation of formulas see: 3-D LINEAR REGRESSION by Jean Jacquelin

 

ymg

 

;; fl3d, Orthogonal Regression in 3D   by ymg                                  ;
;;                                                                             ;
;; For Derivation of Formulas See Following Paper:                             ;
;;        3-D LINEAR REGRESSION  by Jean Jacquelin                             ;
;; http://www.scribd.com/doc/31477970/Regressions-et-trajectoires-3D           ;
;;                                                                             ;

(defun c:fl3d (/ *acaddoc* a avg b c0 c1 c2 cos2a cosa errl go k00 k01 k10 k11 k12 k22
               n p phi pl q r rho s2m sin2a sina ss sum sxx sxy sxz syy syz szz u v w xm ym zm)
  (vl-load-com)
  (defun *error* (msg)
     	(mapcar 'eval errl)
(if (and msg (not (wcmatch (strcase msg) "*BREAK*,*CANCEL*,*EXIT*")))
          (princ (strcat "\nError: " msg))
       )
(and *AcadDoc* (vla-endundomark *AcadDoc*))
       (princ)
  )
 
  (defun acos (z)  (atan (sqrt (- 1.0 (* z z))) z))
 
  (setq errl '("CMDECHO" "DIMZIN")
        errl (mapcar (function (lambda (a) (list 'setvar a (getvar a)))) errl)
  )
    
    
  (or *AcadDoc*
     (setq *AcadDoc* (vla-get-activedocument (vlax-get-acad-object)))
  )
  
  (princ (strcat "\n Select Points"))
  (if (setq  n 0
     pl nil
     ss (ssget '((0 . "POINT")))
            go (> (sslength ss) 2) 
      )
     (progn
        (vla-startundomark *AcadDoc*)	   
        (setvar 'CMDECHO 0)
        (setvar 'DIMZIN 0)
 (repeat (sslength ss)
         (setq pl (cons (cdr (assoc 10 (entget (ssname ss n)))) pl)
                n (1+ n)
         )
 )
        (setq  avg  (mapcar '/  (apply 'mapcar (cons '+ pl)) (list n n n))
               ;dif  (mapcar '(lambda (a) (mapcar '- a avg)) pl)
               sum  (apply
                       'mapcar
                          (cons '+
                             (mapcar
                                '(lambda (a)
                                     (list (* (car a) (car a))       ;  XX 
                                           (* (car a) (cadr a))      ;  XY 
                                           (* (cadr a) (cadr a))     ;  YY 
                                           (* (car a) (caddr a))     ;  XZ 
                                           (* (caddr a) (caddr a))   ;  ZZ 
                                           (* (cadr a) (caddr a))    ;  YZ 
                                     )
                                 )
                                 pl
                             )
                          )
                    )
               sum  (mapcar '/ sum (list n n n n n n))
               ;; avg contains:  (Xm Ym Zm)                                  ;
               ;; sum contains:  (sXX sXY sYY sXZ sZZ sYZ) Divided by n      ;
                  Xm (car avg) Ym (cadr avg) Zm (caddr avg)
                 Sxx (- (car sum)(* Xm Xm))
                 Sxy (- (cadr sum)(* Xm Ym))
                 Syy (- (caddr sum)(* Ym Ym))
                 Sxz (- (cadddr sum)(* Xm Zm))
                 Szz (- (car (cddddr sum))(* Zm Zm))
                 Syz (- (cadr (cddddr sum))(* Ym Zm))
                   a (/ (atan (/ (* Sxy 2.)(- Sxx Syy))) 2.)         
                Cosa (cos a)        Sina (sin a)
               Cos2a (* Cosa Cosa) Sin2a (* Sina Sina)            
                 K11 (+ (* (+ Syy Szz) Cos2a) (* (+ Sxx Szz) Sin2a) (* -2. Sxy Cosa Sina))
                 K22 (+ (* (+ Syy Szz) Sin2a) (* (+ Sxx Szz) Cos2a) (*  2. Sxy Cosa Sina))        
                 K12 (+ (* (* Sxy -1.)(- Cos2a Sin2a))(* (- Sxx Syy) Cosa Sina))
                 K10 (+ (* Sxz Cosa)(* Syz Sina))
                 K01 (+ (* Sxz Sina -1.)(* Syz Cosa))
                 K00 (+ Sxx Syy)           
                  c2 (* (+ K00 K11 K22) -1.)
                  c1 (+ (* K00 K11)(* K00 K22)(* K11 K22)(* K01 K01 -1.)(* K10 K10 -1.))
                  c0 (+ (* K01 K01 K11)(* K10 K10 K22)(* -1. K00 K11 K22))
                   p (- c1 (/ (* c2 c2) 3.))
                   q (+ (/ (* 2.0 c2 c2 c2) 27.)(/ (* -1. c1 c2) 3.) c0)
                   R (+ (/ (* q q) 4.)(/ (* p p p) 27.))         
        )
        (if (minusp R)
           (progn
               (setq rho (sqrt (/ (* -1. p p p) 27.))
                     phi (acos (/ q (* -2. rho)))
                     s2m (min (+ (/ c2 -3.)(* 2. (expt rho (/ 1. 3.))(cos (/ phi 3.))))
                              (+ (/ c2 -3.)(* 2. (expt rho (/ 1. 3.))(cos (/ (+ phi pi pi) 3.))))
                              (+ (/ c2 -3.)(* 2. (expt rho (/ 1. 3.))(cos (/ (+ phi pi pi pi pi) 3.))))
                         )
               )      
           )
           (setq s2m (+ (/ c2 -3.)(expt (+ (/ q -2.)(sqrt R)) (/ 1. 3.))(expt (- (/ q -2.)(sqrt R)) (/ 1. 3.))))
        )

        (setq a (+ (* (/ (* -1. K10)(- K11 s2m)) Cosa) (* (/ K01 (- K22 s2m)) Sina))
              b (+ (* (/ (* -1. K10)(- K11 s2m)) Sina) (* (/ (* -1. K01) (- K22 s2m)) Cosa))
              u (* (/ 1. (+ 1. (* a a) (* b b))) (+ (* (+ 1. (* b b)) Xm) (* -1. a b Ym) (* a Zm)))
              v (* (/ 1. (+ 1. (* a a) (* b b))) (+ (* -1 a b Xm) (* (+ 1 (* a a)) Ym) (* b Zm)))
              w (* (/ 1. (+ 1. (* a a) (* b b))) (+ (* a Xm) (* b Ym) (* (+ (* a a) (* b b)) Zm)))
        )      
        (entmakex (list '(0 . "LINE") (cons 10 avg) (cons 11 (list u v w))))
        
        
     ) 
  )
  (*error* nil)
)

(princ "\nFit a 3D Line to Selection Set of Points, Type FL3D to run")

[Guest]

I use this lisp code

 

;Tip1677:  LIN.LSP   Best-fit Line   (c)2001, Jorn Anke
; PURPOSE: You can select several points and the routine will calculate (using
;          Least Square Method) the line that fits best between the points, and 
;          then draw  the line.

(defun C:LIN ( / xmin  xmax )

; Select points

;;;  (print)
;;;  (princ
;;;    "Vil du gi inndata som punkter eller objekter ? (For objekter benyttes"
;;;  ) ;_ end of princ
;;;  (print)
;;;  (setq SVAR
;;;         (getstring
;;;           "innsettingspunktene til blokker, sirkler etc.) (P/O) <O> ?"
;;;         ) ;_ end of getstring
;;;  ) ;_ end of setq
;;;
;;;  (print)
;;;  (princ
;;;    "Do you vant to select points by snapping - or by selecting objects ?"
;;;  ) ;_ end of princ

 (print)
 (setq SVAR (getstring
              "Do you vant to select points or objects? (P/O) <O> ?"
            ) ;_ end of getstring
 ) ;_ end of setq

 (setq SVAR (strcase SVAR))
 (setq X 0)
 (setq I 0)
 (setq SX 0)
 (setq SY 0)
 (setq SXY 0)
 (setq SX2 0)
 (if (= SVAR "P")
   (progn ; Svar = Points
     (while X ; while X <> nil, e.g. X has got a value
       (setq PUNKT (getpoint "Select point, <Return> to end ...\n"))
       (setq X (car PUNKT)   Y (cadr PUNKT))
       (if X
         (progn
           (setq SX (+ SX X)    SY (+ SY Y)  
                 SXY (+ SXY (* X Y)) 
                 SX2 (+ SX2 (* X X)) )
           (if (= I 0)    (setq XMIN X   XMAX X)  ) ;_ end of if
           (if (> XMIN X) (setq XMIN X) ) ;_ end of if
           (if (< XMAX X) (setq XMAX X) ) ;_ end of if
           (setq I (+ I 1))
         ) ;_ end of progn
       ) ;_ end of if
     ) ;_ end of while
   ) ;_ end of progn
   (progn ; Svar = Objects
     (setq LISTE (ssget)    LISTELENGDE (sslength LISTE))
     (while (<= I (- LISTELENGDE 1))
       (setq ELEMENT (entget (ssname LISTE I)))
       (setq X (cadr (assoc 10 ELEMENT)))
       (setq Y (caddr (assoc 10 ELEMENT)))
       (setq SX (+ SX X))
       (setq SY (+ SY Y))
       (setq SXY (+ SXY (* X Y)))
       (setq SX2 (+ SX2 (* X X)))
       (if (= I 0)
         (progn
           (setq XMIN X)
           (setq XMAX X)
         ) ;_ end of progn
       ) ;_ end of if
       (if (> XMIN X)
         (setq XMIN X)
       ) ;_ end of if
       (if (< XMAX X)
         (setq XMAX X)
       ) ;_ end of if
       (setq I (+ I 1))
     ) ;_ end of while
   ) ;_ end of progn
 ) ;_ end of if
 (setq A (/ (- SXY (/ (* SX SY) I)) (- SX2 (/ (* SX SX) I))))
 (setq B (- (/ SY I) (* A (/ SX I))))
 (setq X1 XMIN)
 (setq Y1 (+ (* A X1) B))
 (setq X2 XMAX)
 (setq Y2 (+ (* A X2) B))
 (command "line" (list X1 Y1) (list X2 Y2) "")
) ;_ end of defun
;***** END RUTINE: LIN




; Linear Least square fit of 2D points 
(defun Lsq_pL (  pL  / xmin  xmax  kA  kB )
  (setq X 0   i 0   kX 0  kY 0  kXY 0  kX2 0
        XMiN  (car p)   XMAX  XMiN  )   
  (foreach p pL 
    (setq x   (car p)    kX (+ kX X) 
          y   (cadr p)   kY (+ kY Y)  
          kXY (+ kXY (* X Y))    i (+ i 1)
          kX2 (+ kX2 (* X X))  )
    (if (> XMiN X) (setq XMiN X) )    ; X min-max vaLs 
    (if (< XMAX X) (setq XMAX X) )  ) ; pL
  (setq kA (/ (- kXY (/ (* kX kY) i)) 
              (- kX2 (/ (* kX kX) i))))
  (setq kB (- (/ kY i) (* kA (/ kX i))))
  (List (List XMiN (+ (* kA XMiN) kB)) 
        (List XMAX (+ (* kA XMAX) kB)) )   ) 

[ymg3]

Prodromosm,

 

I know this one, it's from Cadalyst and does Ordinary Least Square.

 

This is not ideal in you want to fit a line to point that you survey.

The one by Stefan is better for that.

 

The new one is for fitting in 3D and will also work with a 2D line,

however could be more sensitive to roundoff. I have not checked.

 

ymg