关于点到平面距离问题讨论

树櫴希德 · 发表于 2022-10-28 17:51

关于点到平面距离问题讨论

点到平面距离公式是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 (-(* (- 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 (* -1 (+(* a x1) (* b y1) (* c z1)) )
)
;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 height scale / pt pt1 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 70 0)
(cons 74 2)
)
)
;;;;;;;;;;;;;;;;;;;;;;;
;;;插入属性
;;;结束标志
(entmake '((0 . "SEQEND")))
(princ)
)
;;;;;;;;===========================================

(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  (-(* (- 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    (* -1 (+(* a  x1) (* b  y1) (* c  z1)) ) 
       
       ) 
;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 height  scale  / pt  pt1 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 70  0) 
              (cons 74 2) 
           ) 
   ) 
;;;;;;;;;;;;;;;;;;;;;;; 
;;;插入属性 
   
   
   ;;;结束标志 
   (entmake '((0 . "SEQEND"))) 
   (princ) 
) 
 
 
;;;;;;;;===========================================

树櫴希德 · 发表于 2023-4-25 22:35

树櫴希德发表于 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 16:40

本帖最后由树櫴希德于 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)))
(- (* (car v) (caddr u)) (* (car u) (caddr v)))
(- (* (car u) (cadr v)) (* (car v) (cadr u)))
)
) ;@[偏爱云~小吴]偏
(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 函数来自[x_s_s_1]生无可恋大神) "))
(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 pt p0 p1 p2 )))
(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)
)

;;反正弦函数(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))) 
    (- (* (car  v) (caddr u)) (* (car  u) (caddr v))) 
    (- (* (car  u) (cadr  v)) (* (car  v) (cadr  u))) 
  ) 
)  ;@[偏爱云~小吴]偏 
 
(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 函数来自[x_s_s_1]生无可恋大神) ")) 
  (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 pt  p0  p1  p2 ))) 
(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 22:43

本帖最后由 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|/......
(sqrt (+ (expt a 2) (expt b 2)(expt c 2) ))
) )

树櫴希德 · 发表于 2022-10-28 22:21

不知道是精度有问题还是公式错了本来该10.446的

树櫴希德 · 发表于 2022-10-30 09:46

tigcat 发表于 2022-10-29 22:43
楼主图片展示的似乎是一个平面上各点的线性插值
感觉就像是在平面上线性插值求标高.

分子加上abs是对的但现在差的很小如果改正是完全精确的那么待求点与另外三点可以共面在CAD里面可以拉伸extrude验证事实上不能差少少是不是vlax-curve-getpointatparam函数获取的三维点精度不够还是计算中公式精度不够？跟UNITS有没有关系？想不明白

tigcat · 发表于 2022-10-30 10:21

curve系列函数有缺陷，当坐标特别大时会出错。

gzxl · 发表于 2022-10-30 14:30

本帖最后由 gzxl 于 2022-10-30 14:32 编辑

印象中 (expt a 2) 要改成 (expt a 2.)，不然会影响精度。
会不会就不清楚了，毕竟现在不接触 lisp 了。

树櫴希德 · 发表于 2022-10-30 20:18

tigcat 发表于 2022-10-30 10:21
curve系列函数有缺陷，当坐标特别大时会出错。

估计是不知道换乘GETPOINT怎么样

树櫴希德 · 发表于 2022-10-30 20:46

gzxl 发表于 2022-10-30 14:30
印象中 (expt a 2) 要改成 (expt a 2.)，不然会影响精度。
会不会就不清楚了，毕竟现在不接触 lisp 了。

试了没有效果

树櫴希德 · 发表于 2023-4-25 21:18

本帖最后由树櫴希德于 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)))

账号		自动登录	找回密码
密码			注册

关于点到平面距离问题讨论

本帖子中包含更多资源

本帖子中包含更多资源

本帖子中包含更多资源