谢谢Stefan!
您的代码产生了良好的结果。 @李
您的代码适用于三维点还是仅适用于二维点?我在一批3D点上试用了它,结果很奇怪,而Stefan的结果似乎很管用。
经过进一步检查,Stefan的代码生成了一条穿过三维点的连接线,但它只是一个二维实体。
有没有可能仅仅通过点簇创建一条与XLINE或光线平行的线?
嗨,小家伙,
由于我的代码使用简单线性回归,它将只适用于2变量系统,即2D点;当然,更高维度的线性回归是可能的,但要复杂一些。
谢谢李,我想可能是这样的。我有来自3D点云的点子集,大约200个点,需要获得点的最佳3D拟合线。有什么想法可以实现吗? 我的LISP确实很生疏,所以我在VBA中实现了它,并在点簇中画了一条线。我现在打算做进一步的线性回归,但使用点的X和Z坐标来获得通过点的线的梯度。三维线性回归所需的数学知识超出了我的能力。 在原始问题的背景下,我认为Stefan的解决方案优于
对直线应用正交回归,并最小化点到直线的距离。
Lee的解决方案(杰出的实现)是另一个beast,它将沿Y坐标到直线的距离缩小,当X被视为独立变量时,它是完全正确的。
该问题的解决方案首次发表在R.J.ADCOCK出版的《分析员》第5卷第2期(1878年3月)第53-54页。
这是本文的链接:最小二乘法中的一个问题
它也称为戴明回归或正交回归。
ymg公司
这是一张来自维基百科的图片,展示了斯特凡的解决方案:
http://upload.wikimedia.org/wikipedia/commons/thumb/8/81/Total_least_squares.svg/220px-Total_least_squares.svg.png
以下是李正在最小化的内容:
非常有趣且信息丰富的比较ymg-感谢您提供的信息! 泰克,
这是3D中的正交回归。
可能已经太晚了,但我正忙于其他事情。
代码还可以,但我仍然需要添加一个残差表来分析结果。
公式推导见:Jean Jacquelin的三维线性回归
ymg公司
;; fl3d, Orthogonal Regression in 3D by ymg ;
;; ;
;; For Derivation of Formulas See Following Paper: ;
;; 3-D LINEAR REGRESSIONby 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 (setqn 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)
)
)
(setqavg(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")
我使用这个lisp代码
5 Prodromosm,
我知道这一个,它来自Cadalyst,做普通最小二乘法。
这并不理想,因为您希望将一条线拟合到您测量的点。
斯特凡的那部更好。
新的用于三维拟合,也适用于二维线,
然而,可能对舍入更加敏感。我没有检查。
ymg公司
页:
1
[2]