Build an average line from a family of points

January 17, 2013

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

January 17, 2013

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

January 17, 2013

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

Agreed, look into the Least Squares Method.

January 17, 2013

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.

January 17, 2013

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)
)

January 17, 2013

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.

January 17, 2013

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

January 18, 2013

Lee Mac can be on arc or circle

January 18, 2013

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

January 18, 2013

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.

January 18, 2013

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.

January 23, 2013

@ 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

January 23, 2013

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.

January 23, 2013

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?

January 23, 2013

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.

October 29, 2013

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

October 29, 2013

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

November 7, 2013

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")

November 7, 2013

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)) )   )

November 7, 2013

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

Sign In

Build an average line from a family of points

Recommended Posts

MSasu

Top Posters In This Topic

Popular Days

Top Posters In This Topic

Popular Days

Popular Posts

Stefan BMR

Posted Images

Organic

Lee Mac

MSasu

Lee Mac

MSasu

Lee Mac

motee-z

Stefan BMR

MSasu

Lee Mac

Tyke

Lee Mac

Tyke

Tyke

ymg3

Lee Mac

ymg3

Guest

ymg3

Join the conversation

Stefan BMR

Tutorials

Forums

Activity