Build an average line from a family of points

[Guest]

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

 

I am use block attribiut points and Stefan lisp not working to me ..........

[Guest]

Salut Mircea

 

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

 

can you convert it to support attribute blocks. Here is the block i use

Point.dwg

[ymg3]

prodromosm,

 

Just change the filter in the ssget line to the following and it will select either

both your point block and any autocad points.

 

(setq ss (ssget '((-4 . "<OR")
                         (0 . "POINT")
                         (-4 . "<AND")(0 . "INSERT")(2 . "Point")(-4 . "AND>")
                     (-4 . "OR>")))
)       

 

If you change this way, only your Insert "Point" will be selected.

 

(setq ss (ssget '((-4 . "<AND")(0 . "INSERT")(2 . "Point")(-4 . "AND>")))) 

 

ymg

Edited November 7, 2013 by ymg3

[Guest]

If i need this lisp to support more blocks with different names can i write !!!!!!!!!!!

(setq ss (ssget '((-4 . "<OR")
                         (0 . "POINT")
                         (-4 . "<AND")(0 . "INSERT")(2 . "Point")(2 . "block2")(2 . "block3")(2 . "block 4")(-4 . "AND>")
                     (-4 . "OR>")))
)

[MSasu]

May just use:

(2 . "Point,block2,block3,block 4")

[Guest]

thanks MSasu

[MSasu]

You're welcome! For more options you may want to check the help of WCMATCH function.

[Stefan BMR]

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

ymg, I've made some tests in 2D and your routine gives a zero length line in average point.

This is my routine for 3d points, based on fact that seeking line has direction of an eigenvector of symmetrical matrix:

((Ixx Ixy Ixz)
(Iyx Iyy Iyz)
(Izx Izy Izz))

My eigenvector function uses a numerical method, hence it is an approximation.

;vectors&matrix
(defun vxv (u v) (apply '+ (mapcar '* u v)))                                ;dot_product
(defun vxs (v s) (mapcar '(lambda (a) (* a s)) v))                        ;scale vector
(defun mod (v)   (sqrt (vxv v v)))                                                  ;vector's module
(defun mxv (m v) (mapcar '(lambda (u) (vxv u v)) m))                        ;matrix x vector
(defun m+m (m n) (mapcar '(lambda (a b) (mapcar '+ a b)) m n))                        ;add matrix
(defun m-m (m n) (mapcar '(lambda (a b) (mapcar '- a b)) m n))                        ;matrix subs.
(defun mxs (m s) (mapcar '(lambda (a) (vxs a s)) m))                         ;scale matrix
(defun trm (m)   (apply 'mapcar (cons 'list m)))                        ;transpose matrix
(defun mxm (m n) (mapcar '(lambda (a) (mapcar '(lambda (b) (vxv a b)) (trm n))) m)) ;matrix product
(defun unv (v)   (vxs v (/ 1. (mod v))))                            ;unit vector

;eigenvector coresponding to greatest eigenvalue
(defun eigenvec (m / u u0 ev)
 (repeat (length m) (setq u0 (cons 1.0 u0)))
 (setq u (mxv m u0)
       ev (car u0)
       )
 (while
   (not (equal ev (setq ev (/ (car u) (car u0))) 1e-12))
   (setq u0 u
         u  (unv (mxv m u0))
   )
 )
 u
)

(defun C:TEST ( / ss i l o dx dy dz Ixx Ixy Ixz Iyy Iyz Izz n)
 (if
   (setq ss (ssget '((0 . "POINT"))))
   (progn
     (repeat (setq i (sslength ss))
       (setq l (cons (cdr (assoc 10 (entget (ssname ss (setq i (1- i)))))) l))
     )
     (setq o   (vxs (apply 'mapcar (cons '+ l)) (/ 1.0 (length l)))
           l   (mapcar '(lambda (p) (mapcar '- p o)) l)
           dx  (mapcar 'car l) dy  (mapcar 'cadr l) dz  (mapcar 'caddr l)
           Ixx (vxv dx dx)     Iyy (vxv dy dy)      Izz (vxv dz dz)
           Ixy (vxv dx dy)     Iyz (vxv dy dz)      Ixz (vxv dx dz)
           n   (eigenvec (list (list Ixx Ixy Ixz)
                               (list Ixy Iyy Iyz)
                               (list Ixz Iyz Izz)
                         )
               )
     )
     (entmake (list '(0 . "XLINE") '(100 . "AcDbEntity") '(100 . "AcDbXline") (cons 10 o) (cons 11 n)))
   )
 )
 (princ)
)

[ymg3]

Hi Stefan,

 

My function will indeed not work if all the points are on the same plane.

 

I like your code, very clean!!

 

Thanks,

 

ymg

Edited November 17, 2013 by ymg3

[satishrajdev]

Can anyone provide code for Best fit Plane 3D?

[Stefan BMR]

Can anyone provide code for Best fit Plane 3D?

I can only guess that this is the solution.

The plan is defined by 2 perpendicular Xlines.

;Plan approximation of 3D points
;Stefan M. 28.06.2016

;vectors&matrix
(defun vxv (u v) (apply '+ (mapcar '* u v)))                                        ;dot_product
(defun vxs (v s) (mapcar '(lambda (a) (* a s)) v))                                  ;scale vector
(defun mod (v)   (sqrt (vxv v v)))                                                  ;vector's module
(defun mxv (m v) (mapcar '(lambda (u) (vxv u v)) m))                                ;matrix x vector
(defun m+m (m n) (mapcar '(lambda (a b) (mapcar '+ a b)) m n))                      ;add matrix
(defun m-m (m n) (mapcar '(lambda (a b) (mapcar '- a b)) m n))                      ;matrix subs.
(defun mxs (m s) (mapcar '(lambda (a) (vxs a s)) m))                                ;scale matrix
(defun trm (m)   (apply 'mapcar (cons 'list m)))                                    ;transpose matrix
(defun mxm (m n) (mapcar '(lambda (a) (mapcar '(lambda (b) (vxv a b)) (trm n))) m)) ;matrix product
(defun unv (v)   (vxs v (/ 1. (mod v))))                                            ;unit vector

;eigenvector coresponding to greatest eigenvalue
(defun eigenvec (m / u u0 ev)
 (repeat (length m) (setq u0 (cons 1.0 u0)))
 (setq u (mxv m u0)
       ev (car u0)
       )
 (while
   (not (equal ev (setq ev (/ (car u) (car u0))) 1e-12))
   (setq u0 u
         u  (unv (mxv m u0))
   )
 )
 u
)

(defun C:TEST ( / a b i l o p ss x)
 (if
   (setq ss (ssget '((0 . "POINT"))))
   (progn
     (repeat (setq i (sslength ss))
       (setq l (cons (cdr (assoc 10 (entget (ssname ss (setq i (1- i)))))) l))
     )
     (setq o   (vxs (apply 'mapcar (cons '+ l)) (/ 1.0 (length l)))
           l   (mapcar '(lambda (p) (mapcar '- p o)) l)
           a   (average_line l)
           b   (trans (average_line (mapcar '(lambda (x) (2dpoints (trans x 0 a t))) l)) a 0 t)
     )
     (entmake (list '(0 . "XLINE") '(100 . "AcDbEntity") '(100 . "AcDbXline") (cons 10 o) (cons 11 a)))
     (entmake (list '(0 . "XLINE") '(100 . "AcDbEntity") '(100 . "AcDbXline") (cons 10 o) (cons 11 b)))
   )
 )
 (princ)
)

(defun 2dpoints (p) (list (car p) (cadr p) 0.0))

;Linear approximation of 3D points
(defun average_line (l / dx dy dz Ixx Ixy Ixz Iyy Iyz Izz n)
 (setq
   dx  (mapcar 'car l) dy  (mapcar 'cadr l) dz  (mapcar 'caddr l)
   Ixx (vxv dx dx)     Iyy (vxv dy dy)      Izz (vxv dz dz)
   Ixy (vxv dx dy)     Iyz (vxv dy dz)      Ixz (vxv dx dz)
   n   (eigenvec
         (list
           (list Ixx Ixy Ixz)
           (list Ixy Iyy Iyz)
           (list Ixz Iyz Izz)
         )
       )
 )
)

[satishrajdev]

Seems like great approximation but not 100% accurate. Z for for 4 points can't be same...

[Stefan BMR]

Seems like great approximation but not 100% accurate. Z for for 4 points can't be same...

show me the points.

[satishrajdev]

Best Fit Plane 3D.dwg

 

Please find attachment.

[Stefan BMR]

You didn't sent me 4 points and you didn't explain the approximation error you are talking about.

 

I've ran some tests and here are the conclusions:

- for 3 points in 3D, the plan contain each point, as it should be.

- for a symmetric cloud of points on a ellipsoid surface, the plan passes through the main axis of the ellipsoid.

 

It seems the result is just fine.

[satishrajdev]

I think I didn't explain it properly... Can you make modifications in your code where I can get line/polyline instead of xline?

 

I tried it from my end... But failed with it.... So, If once I get the length of line/polyline I can easily convert them into Plane (3D face).