关于点到平面距离问题讨论
关于点到平面距离问题讨论点到平面距离公式是d=|Ax0+By0+Cz0+D|/√(A2+B2+C2)。
已知三点p1（x1,y1,z1），p2(x2,y2,z2)，p3(x3,y3,z3)，要求确定的平面方程
点到平面距离是指空间内一点到平面内一点的最小长度。特殊的，当点在平面内时，该点到平面的距离为0。公式中的平面方程为Ax+By+Cz+D=0
;A = (y2 - y1)*(z3 - z1) - (z2 -z1)*(y3 - y1);
;B = (x3 - x1)*(z2 - z1) - (x2 - x1)*(z3 - z1);
;C = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
;
;又D = -(A * x1 + B * y1 + C * z1)，所以可以求得D的值。
(defun vxs (e / i v lst)
(vl-load-com)
(setq i 0)
(while
(setq v (vlax-curve-getpointatparam e (setq i (1+ i))))
(setq lst (cons v lst))
)
(reverse lst)
)
(defun c:ddm (/ x1 x2 x3 y1 y2 y3 z1 z2 z3 a b c d dd pts pt1 )
;;;;;;;;;;;;;;;
(setq pts (vxs(car(entsel "\n请选择三角网中一个三角形3dpoly线："))))
(setq x1(car(car pts))
y1 (cadr(car pts))
z1 (nth 2 (car pts))
x2(car(cadr pts))
y2 (cadr(cadr pts))
z2 (nth 2 (cadr pts))
x3 (car(nth 2 pts))
y3(cadr(nth 2 pts))
z3 (nth 2 (nth 2 pts))
)
(setq a(-(* (- y2y1)(- z3z1)) (*(- z2 z1)(- y3y1)) )
b(-(* (- x3x1)(- z2z1)) (*(- x2x1)(- z3z1)) )
c(-(* (- x2x1)(- y3y1)) (*(- x3x1)(- y2y1)) )
d (* -1 (+(* ax1) (* by1) (* cz1)) )
)
;A = (y2 - y1)*(z3 - z1) - (z2 -z1)*(y3 - y1);
;B = (x3 - x1)*(z2 - z1) - (x2 - x1)*(z3 - z1);
;C = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
;
;又D = -(A * x1 + B * y1 + C * z1)，所以可以求得D的值。
;点到平面距离公式是d=|Ax0+By0+Cz0+D|/√(A2+B2+C2)
(setq pt1(getpoint "\n请选择待求点："))
(setq dd (/ (+(* a(car pt1))(* b(cadr pt1))(* c(nth 2 pt1)) d)
(sqrt (+ (expt a 2) (expt b 2)(expt c 2) ))
) )
;22、 点
;(entmake (list '(0 . "POINT") (cons 10 (list (car pt1) (cadr pt1)(- (nth 2 pt1)dd)))))
(mkgcd (list (car pt1) (cadr pt1)(- (nth 2 pt1)dd))(- (nth 2 pt1)dd)1)
(princ)
)
(defun mkgcd (inspt heightscale/ ptpt1 blkdef obj)
(setvar "CMDECHO" 0)
(command "layer" "m" "检查高程点" "c" "1" "" "L" "CONTINUOUS" """")
(if height
(setq height (rtos height 2 3))
(setq height "")
)
(regapp "SOUTH")
;;;检查字体 "HZ" 是否存在
(if (not (tblobjname "style" "HZ"))
(command "style" "HZ" "rs.shx,hztxt.shx" 0 1 0 "" "" "")
)
;;;检查是否存在高程点图块定义
(if (not (tblobjname "block" "GC2000"))
(progn
;13、entmake生成普通块
(defun emkblk ( pt name /)
(entmake (list '(0 . "block") (cons 2 name) '(70 . 0) (cons 10 pt)))
(entmake (list '(0 . "LWPOLYLINE") '(100 . "AcDbEntity") '(100 . "AcDbPolyline") (cons 90 4) (cons 10 (list (+ (car pt) 0.75)(+ (cadr pt) 1) ))(cons 10 pt) (cons 10 (list (- (car pt) 0.75)(+ (cadr pt) 1) ))
(cons 10 (list (+ (car pt) 4.25)(+ (cadr pt) 1) ))
))
(entmake '((0 . "ENDBLK")))
;(entmake (list '(0 . "INSERT") (cons 2 name) (cons 10 pt)))
)
(emkblk '(0 0) "GC2000")
)
)
;;;插入块
(entmake (list
'(0 . "INSERT")
'(100 . "AcDbEntity")
'(100 . "AcDbBlockReference")
'(66 . 1);;;属性跟随标志,1跟随，0不跟随
(cons 2 "GC2000")
(cons 10 inspt)
(cons 41 scale)
(cons 42 scale)
(cons 43 scale)
'(-3 ("SOUTH" (1000 . "202101")))
)
)
;;;插入属性
(entmake (list
'(0 . "ATTRIB")
'(100 . "AcDbEntity")
'(100 . "AcDbText")
(cons 10 (setq pt (polar inspt (* 0.5 PI) (* 2.25 scale))))
(cons 40 (* 2.0 scale))
(cons 50 0)
(cons 62 3)
(cons 41 0.8)
(cons 51 0)
(cons 1 height)
(cons 7 "HZ")
(cons 72 0)
(cons 11 pt)
'(100 . "AcDbAttribute")
(cons 2 "height")
(cons 700)
(cons 74 2)
)
)
;;;;;;;;;;;;;;;;;;;;;;;
;;;插入属性
;;;结束标志
(entmake '((0 . "SEQEND")))
(princ)
)
;;;;;;;;===========================================
树櫴希德 发表于 2023-4-25 21:18
;;反正弦函数(asin 值)
(defun asin(x)
(cond
((or (> x 1) (< x -1)) nil)
((= x 1) (* pi 0.5))
((= x -1) (* pi -0.5))
(t (atan (/ x (sqrt (- 1 (* x x))))))
)
)
;;反余弦函数(acos 值)
(defun acos(x)
(if (>= x 0)
(asin (sqrt (- 1 (* x x))))
(+ pi (- (asin (sqrt (- 1 (* x x))))))
)
) 本帖最后由 树櫴希德 于 2023-4-25 23:24 编辑
;;反正弦函数(asin 值)
(defun asin(x)
(cond
((or (> x 1) (< x -1)) nil)
((= x 1) (* pi 0.5))
((= x -1) (* pi -0.5))
(t (atan (/ x (sqrt (- 1 (* x x))))))
)
)
;;反余弦函数(acos 值)
(defun acos(x)
(if (>= x 0)
(asin (sqrt (- 1 (* x x))))
(+ pi (- (asin (sqrt (- 1 (* x x))))))
)
)
;;;=============================================
;;; 通用函数点到平面投影
;;;参数： pt------原始点
;;; pt0/1/2------面上三点
;;;返回值：投影点
(defun xty-G-pertoface(pt pt0 pt1 pt2 / nv e)
(setq nv (mat:vxv (mapcar (function -) pt1 pt0)
(mapcar (function -) pt2 pt0)))
(setq
e (/ (-(apply (function +) (mapcar (function *) nv pt0))
(apply (function +) (mapcar (function *) nv pt)))
(apply (function +) (mapcar (function *) nv nv))))
(mapcar (function +)
pt
(mapcar (function (lambda (x) (* x e))) nv)))
(defun mat:vxv ( u v )
(list
(- (* (cadr u) (caddr v)) (* (cadr v) (caddr u)))
(- (* (carv) (caddr u)) (* (caru) (caddr v)))
(- (* (caru) (cadrv)) (* (carv) (cadru)))
)
);@[偏爱云~小吴]偏
(defun xty-g-pertoface1 (pt pt0 pt1 pt2 / nv e)
(setq nv (mat:vxv (mapcar '- pt1 pt0) (mapcar '- pt2 pt0)))
(setq e (/ (- (apply '+ (mapcar '* nv pt0)) (apply '+ (mapcar '* nv pt))) (apply '+ (mapcar '* nv nv))))
(mapcar '+ pt (mapcar (function (lambda (x) (* x e))) nv))
)
; 53829.5,221425.0
;(apply'xty-G-pertoface (list(getpoint) (getpoint) (getpoint) (getpoint) ))
(defun PtPerFace (pt p0 p1 p2 / nv e)(setq nv (mat:vxv(mapcar '- p1 p0)(mapcar '- p2 p0))e(/(- (apply '+ (mapcar '* nv p0))(apply '+ (mapcar '* nv pt)))(apply '+ (mapcar '* nv nv))))(mapcar '+ pt(mapcar'(lambda(x)(* x e)) nv)))
;;
(defun c:ddpm1 ( /p0 p1 p2 pt p11 p12 juli gcc xcc jd xcd )
(alert (strcat "\n ddpm1""\n绿色点为原位重力方向相较于平面的点 " "\n 白色为待求点垂直于平面的点""\n 函数来自生无可恋大神) "))
(setq p0 (getpoint "\n请指定面上第1点："))
(setq p1 (getpoint "\n请指定面上第2点："))
(setq p2 (getpoint"\n请指定面上第3点："))
(setq pt (getpoint"\n请指定待求点："))
(setq p11 (apply'xty-G-pertoface (list ptp0p1p2 )))
(entmake (list '(0 . "POINT")(cons 62 7) (cons 10 p11)))
(setq juli (sqrt (+ (expt (-(car pt) (car p11)) 2.0) (expt (-(cadr pt) (cadr p11)) 2.0))) )
(setq gcc (- (last p11) (last pt) ) )
(setq xcc (sqrt (+(expt juli 2.0) (expt gcc 2.0)) ))
(setq jd (atan (/(abs gcc) juli) ))
(setq xcd (/ xcc (sin jd) ))
(if(>= gcc 0)
(setq xcd xcd)
(setq xcd (* -1 xcd))
)
(setq p12 (list (car pt ) (cadr pt ) (+ (last pt ) xcd) ))
(entmake (list '(0 . "POINT") (cons 62 3)(cons 10 p12)))
(if( < xcd 0)
(alert (strcat "\n待求点比平面高" (rtos (* -1 xcd)2 3) "米"))
(alert (strcat "\n待求点比平面低" (rtos (* 1 xcd)2 3) "米"))
)
(princ)
) 本帖最后由 tigcat 于 2022-10-29 23:01 编辑
楼主图片展示的似乎是一个平面上各点的线性插值
感觉就像是在平面上线性插值求标高.
如果划分成三角网,三角网三个顶点有三个坐标,从而可以确定一个平面方程Ax+By+Cz+D=0,在这个三角网内确定Xi,yi后都有唯一的Zi对应.
如果是怀疑公式在程序中出错,可以考虑换trans方式对比计算结果.
没有dwg,lisp程序中有一句比较可疑:
[*](setq dd (/ (+(* a(car pt1))(* b(cadr pt1))(* c(nth 2 pt1)) d);分子应该加上abs,保证是正值|Ax+By+Cz+D|/......
[*]http://bbs.mjtd.com/source/plugin/imc_colorcode/images/jssc_none.gif
[*]http://bbs.mjtd.com/source/plugin/imc_colorcode/images/jssc_none.gif (sqrt (+ (expt a 2) (expt b 2)(expt c 2) ))
[*]http://bbs.mjtd.com/source/plugin/imc_colorcode/images/jssc_none.gif ) )
不知道是精度有问题还是公式错了 本来该10.446的 tigcat 发表于 2022-10-29 22:43
楼主图片展示的似乎是一个平面上各点的线性插值
感觉就像是在平面上线性插值求标高.
分子加上abs是对的 但现在差的很小 如果改正是完全精确的 那么待求点与另外三点可以共面 在CAD里面可以拉伸extrude验证事实上不能 差少少是不是vlax-curve-getpointatparam函数获取的三维点精度不够 还是计算中公式精度不够？跟UNITS有没有关系？想不明白 curve系列函数有缺陷，当坐标特别大时会出错。 本帖最后由 gzxl 于 2022-10-30 14:32 编辑
印象中 (expt a 2) 要改成 (expt a 2.)，不然会影响精度。
会不会就不清楚了，毕竟现在不接触 lisp 了。
tigcat 发表于 2022-10-30 10:21
curve系列函数有缺陷，当坐标特别大时会出错。
估计是 不知道换乘GETPOINT怎么样 gzxl 发表于 2022-10-30 14:30
印象中 (expt a 2) 要改成 (expt a 2.)，不然会影响精度。
会不会就不清楚了，毕竟现在不接触 lisp 了。
试了 没有效果 本帖最后由 树櫴希德 于 2023-4-25 23:29 编辑
树櫴希德 发表于 2023-4-25 16:40
(setq p11 (apply'xty-G-pertoface (list(getpoint) (getpoint) (getpoint) (getpoint) )))
(entmake (list '(0 . "POINT") (cons 10 p11)))
页:
[1]
2