曲线的转弯半径和曲率

在下面的这个帖子中讨论了椭圆的曲率和转弯半径
http://bbs.mjtd.com/thread-62980-1-1.html
现在我把这个主题深化一下,讨论一下曲线的两个函数:
vlax-curve-getSecondDeriv
vlax-curve-getFirstDeriv
这两个函数是什么意思呢?
我们考察AutoCAD里面的曲线类,主要是圆,椭圆,弧和样条曲线,多段线由这几种组合而成。
椭圆和样条曲线实际上都是由参数形成,因此,对于这类曲线,它们每点的坐标可以由参数方程表达:
譬如椭圆 x=a*cos(t); y=b*sin(t);
样条曲线也有方程,假设样条曲线的参数方程为: X= f(t);
Y=g(t);
因此可以对参数方程求导,得到每一点的切线矢量,曲线上每一点对应于一个参数 t0 ,
这个切线矢量的 的X值就是 f(t)在t0处的一阶导数,Y值就是g(t)在t0处的一阶导数,
即( f'(t0), g'(t0),0)
因而我们就理解了vlax-curve-getFirstDeriv 函数返回值的意义,
对于vlax-curve-getSecondDeriv的意义类似,只不过这次换成了二阶导数。
即( f''(t0), g''(t0),0)

那么如何求样条曲线或者椭圆的每一点的曲率及其半径呢?
根据参数方程的求曲率公式,若曲线由下面参数方程给出:
\[\left\{\begin{array}{l}x=\varphi(t) \\ y=\psi(t) \end{array}\right.\]
那么曲率如何计算?提示:
\[K=\frac{\left|\varphi^{\prime}(t) \psi^{\prime \prime}(t)-\varphi^{\prime \prime}(t) \psi^{\prime}(t)\right|}{\left[\varphi^{\prime 2}(t)+\psi^{\prime 2}(t)\right]^{3 / 2}}\]
下面是求值代码:

;;;=============================================================
;;;功能: 获取椭圆上一点处的曲率和离心半径(适用于单次利用此法)
;;;参数: 椭圆实体和椭圆上的一点
;;;返回: 此处离心圆圆心、离心半径及其曲率(离心率)
;;;说明: 如果要在CAD中几何作图,可以参考此贴:
;;;      http://bbs.mjtd.com/thread-62980-1-1.html
;;;=============================================================
(defun ELL:GetCurvature (en pt / obj a b px x y v1 v2 rad cen)
  (setq obj (vlax-ename->vla-object en))
  (setq a   (vla-get-MajorRadius obj))                          ;椭圆的半长轴
  (setq b   (vla-get-MinorRadius obj))				;椭圆的半短轴
  (setq pt  (vlax-curve-getclosestpointto en pt))		;保证此点在椭圆上
  (setq px  (vlax-curve-getParamAtPoint en pt))			;此点的椭圆参数
  (setq v1  (vlax-curve-getFirstDeriv en px))			;此点的一阶矢量
  (setq v2  (list (- (cadr v1)) (car v1) (caddr v1)))		;此点的切线矢量
  (setq x   (* a (sin px)))
  (setq y   (* b (cos px)))
  (setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b)))		;得到转弯半径
  (setq cen (polar pt (angle '(0 0 0) v2) rad))			;圆心
  (list cen rad)						;圆心及半径
)
;;;=============================================================
;;; 利用参数方程的求椭圆的曲率及其离心半径(适用于多次利用此法)
;;; 功能: 获取椭圆上一点处的转弯半径和离心圆圆心
;;; 参数: 椭圆实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun ELL:GetCurvature1 (en a b pt / px v1 v2 x y rad cen)
  (setq pt  (vlax-curve-getclosestpointto en pt))		;保证此点在椭圆上
  (setq px  (vlax-curve-getParamAtPoint en pt))			;此点的椭圆参数
  (setq v1  (vlax-curve-getFirstDeriv en px))			;此点的一阶矢量
  (setq v2  (list (- (cadr v1)) (car v1) (caddr v1)))		;此点的切线矢量
  (setq x   (* a (sin px)))
  (setq y   (* b (cos px)))
  (setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b)))		;得到转弯半径
  (setq cen (polar pt (angle '(0 0 0) v2) rad))			;圆心
  (list cen rad)						;圆心及半径
)
;;;=============================================================
;;; 一般平面曲线参数方程的曲率离心公式
;;; 功能: 获取曲线上一点处的离心半径和离心圆圆心
;;; 参数: 曲线实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun CUR:GetCurvature (en pt / ob px v1 v2 v3 x1 y1 x2 y2 cen rad d1 d2)
  (setq ob (vlax-ename->vla-object en))
  (setq pt (vlax-curve-getclosestpointto en pt))                ;保证此点在曲线上
  (setq px (vlax-curve-getParamAtPoint en pt))                  ;此点的曲线参数
  (setq v1 (vlax-curve-getFirstDeriv en px))                    ;此点的一阶矢量
  (setq v2 (vlax-curve-getSecondDeriv en px))			;此点的二阶矢量
  (setq v3 (list (- (cadr v1)) (car v1) (caddr v1)))            ;此点的切线矢量
  (setq x1 (car  v1))                                           ;一阶导数的 X值
  (setq y1 (cadr v1))						;一阶导数的 Y值
  (setq x2 (car  v2))						;二阶导数的 X值
  (setq y2 (cadr v2))						;二阶导数的 Y值
  (setq d1 (expt (+ (* y1 y1) (* x1 x1)) 1.5))
  (setq d2 (- (* x1 y2) (* x2 y1)))				;转弯内外的判定
  (if (/= d2 0)                                                 ;如果不为直线段
    (progn
      (setq rad (/ d1 d2))
      (if (vlax-method-applicable-p ob 'GetBulge)               ;如果为多段线(含圆弧)
	(if (< (vla-GetBulge ob (fix px)) 0)			;如果此段凸度小于0
	  (setq rad (- rad))
        )
      )
      (list (polar pt (angle '(0 0 0) v3) rad) (abs rad))	;圆心及半径
    )
  )
)

一些测试代码:

;;;=============================================================
;;;测试程序1: 获取曲线一点处的曲率和离心半径
;;;=============================================================
(defun c:tt1 (/ ent pnt ret)
  (setq ent (car (entsel "\n请选取曲线:")))
  (while (setq pnt (getpoint "\n点取一点"))
    (setq pnt (trans pnt 1 0))
    (setq pnt (vlax-curve-getclosestpointto ent pnt))
    (setq ret (CUR:GetCurvature ent pnt))
    (princ "\n离心半径是:")
    (princ ret)
    (and ret (apply 'Ent:Make_Circle ret))
  )
  (princ)
)
;;;=============================================================
;;;测试程序2: 获取椭圆一点处的曲率和离心半径,并比较两种方法
;;;=============================================================
(defun c:tt2 ( / ent obj a b pnt par 1st 2st r1 r2 r3)
  (if (and (setq ent (car (entsel "\n请选取椭圆:")))
	   (setq obj (vlax-ename->vla-object ent))
	   (vlax-property-available-p obj 'MajorRadius)		;这个地方应该加出错处理
  	   (setq a (vla-get-MajorRadius obj))
	   (setq b (vla-get-MinorRadius obj))
      )
    (while (setq pnt (getpoint "\n点取一点:"))
      (setq pnt (trans pnt 1 0))
      (setq pnt (vlax-curve-getclosestpointto ent pnt))
      (setq par (vlax-curve-getparamatpoint ent pnt))
      (setq 1st (vlax-curve-getFirstDeriv ent par))
      (setq 2st (vlax-curve-getSecondDeriv ent par))
      (setq r1  (distance '(0 0) 2st))                          ;这个secondDeriv并不意味着半径
      (setq r2  (cadr (ELL:GetCurvature1 ent a b pnt)))
      (setq r3  (cadr (CUR:getcurvature ent pnt)))
      (princ "\nRadius 1 is:")
      (princ r1)
      (princ "\nRadius 2 is:")
      (princ r2)
      (princ "\nRadius 3 is:")
      (princ r3)
      (UTI:Bench
        10000
        (list
	  (list 'ELL:GetCurvature1 ent a b pnt)
	  (list 'CUR:getcurvature ent pnt)
        )
      )
    )
  )
)
;;;=============================================================
;;;测试程序3: 动态演示曲线的离心半径
;;;=============================================================
(defun c:tt3 (/ CIR ENT LIN PNT PT0 RET 1ST 2ST CEN PAR PT2 VEC)
  (setq ent (car (entsel "\n请选取曲线:")))
  (setq cir (Ent:Make_Circle '(0 0 0) 1))
  (setq lin (ent:make_Line '(0 0 0) '(0 0 0)))
  (setq vec (ent:make_Line '(0 0 0) '(0 0 0)))
  (setq lin (vlax-ename->vla-object lin))
  (setq vec (vlax-ename->vla-object vec))
  (setq cir (vlax-ename->vla-object cir))
  (vla-put-color lin 1)
  (vla-put-color cir 3)
  (vla-put-color vec 6)
  (princ "\n按空格或者回车退出!")
  (while (= (car (setq pt0 (grread 5 0))) 5)
    (setq pnt (trans (cadr pt0) 1 0))
    (setq pnt (vlax-curve-getclosestpointto ent pnt))
    (setq par (vlax-curve-getparamatpoint ent pnt))
    (setq 1st (vlax-curve-getFirstDeriv ent par))
    (setq 2st (vlax-curve-getSecondDeriv ent par))
    (if (setq ret (CUR:GetCurvature ent pnt))
      (progn
	(setq pt2 (mapcar '+ pnt 2st))
	(setq cen (car ret))
        (vlax-put Cir 'Center cen )
        (vlax-put cir 'Radius (cadr ret))
        (vlax-put lin 'StartPoint cen)
        (vlax-put lin 'EndPoint pnt)
	(vlax-put vec 'StartPoint pt2)
        (vlax-put vec 'EndPoint pnt)
      )
    )
  )
  (vla-erase cir)
  (vla-erase lin)
  (vla-erase vec)
  (princ)
)
;;;=============================================================
;;;测试程序4: 由曲线的离心半径描绘其轨迹
;;;=============================================================
(defun c:tt4 (/ ent lst par px1 px2 pxn pnt Inf cen)
  (if (setq ent (car (entsel "\n请选取曲线:")))
    (progn
      (setq px1 (vlax-curve-getStartParam ent))
      (setq px2 (vlax-curve-getendparam ent))
      (setq pxn (* (- px2 px1) 0.02))
      (setq par px1)
      (while (<= par px2)
        (setq pnt (vlax-curve-getpointatparam ent par))
        (setq Inf (CUR:GetCurvature ent pnt))
        (if (setq cen (car Inf))
	  (setq lst (cons cen lst))
	)
        (setq par (+ par pxn))
      )
      (setq lst (reverse lst))
      (Ent:Make_LWPoly lst 1)
    )
  )
)

相关联代码:

;;;************************************************************;
;;;实体创建部分
;;;************************************************************;
;;;-------------------------------------------------------------
;;;创建一个点
;;;输入: 一个三维或者二维的点
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Point (p)
  (entmakex (list '(0 . "POINT") (cons 10 p)))
)
;;;-------------------------------------------------------------
;;;创建一个带颜色的点(此函数为测试或者其他用途)
;;;输入: 一个三维或者二维的点表和一个颜色号
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:MakePoint-1 (p c)
  (entmakex (list '(0 . "POINT") (cons 10 p) (cons 62 c)))
)
;;;-------------------------------------------------------------
;;;创建一条直线段
;;;输入: 两个三维或者二维的点
;;;输出: 线段实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Line (p q)
  (entmakeX (list '(0 . "LINE") (cons 10 p) (cons 11 q)))
)
;;;-------------------------------------------------------------
;;;创建一个由三条直线组成的三角形
;;;输入: 三个三维或者二维的点
;;;输出: 由三条直线组成的三角形
;;;-------------------------------------------------------------
(defun Ent:Make_Triangle (p1 p2 p3)
  (mapcar 'Ent:Make_Line (list p1 p2 p3) (list p2 p3 p1))
)
;;;-------------------------------------------------------------
;;;创建一个三维多段线
;;;输入: 三维的点集
;;;输出: 三维多段线实体
;;;-------------------------------------------------------------
(defun Ent:Make_Poly (pts Closed / e)
  (if Closed
    (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 9))))
    (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 8))))
  )
  (foreach p pts
    (entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
  )
  (entmake '((0 . "SEQEND")))
  (entlast)
)
;;;-------------------------------------------------------------
;;;创建轻多段线
;;;输入: 二维的点集
;;;输出: 轻多段线实体名
;;;-------------------------------------------------------------
(defun Ent:Make_LWPoly (pts closed /)
  (entmakeX
    (append
      '((0 . "LWPOLYLINE")
        (100 . "AcDbEntity")
        (100 . "AcDbPolyline")
       )
      (list (cons 90 (length pts)))                      	;顶点个数
      (mapcar (function (lambda (x) (cons 10 x))) pts)  	;多段线顶点
      (list (cons 70  closed))                  	        ;闭合的
    )
  )
)
;;;-------------------------------------------------------------
;;; 创建圆实体
;;; 输入: 圆心C和半径R
;;; 输出: 圆的实体名
;;;-------------------------------------------------------------
(defun Ent:Make_Circle (C R)
  (entmakex (list '(0 . "CIRCLE") (cons 10 C) (cons 40 R)))
)

程序1效果演示:

程序2效果演示:


发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注