(defun C:FC (/)
  ;;不同长度的弧度拟合误差分析
  (defun dif_arc_length(	
				  /		pnt_lst
				  result_lm	result_c
				  result_t	classic_param
				  taubin_param	classic_mse
				  taubin_mse
				 )
    (defun draw_circle (param color /)
      (entmakex	(list '(0 . "CIRCLE")
		      '(8 . "fit_line")
		      (cons 62 color)
		      (cons 10 (car param))
		      (cons 40 (cadr param))
		)
      )
    )
    (defun range (min_ max_ step / i range_lst)
      (setq i	      1
	    range_lst (list min_)
      )
      (repeat (fix (/ (- max_ min_) (float step)))
	(setq range_lst	(cons (+ (* step i) min_) range_lst)
	      i		(1+ i)
	)
      )
      (reverse range_lst)
    )
    (defun draw_tk (gt_cir head / pnt)
      (defun txt_entmake (layer rotate height point text / ent)
	(entmakex (list	'(0 . "TEXT")
			(cons 8 layer)
			(cons 50 rotate)
			(cons 40 height)
			(cons 10 point)
			(cons 1 text)
		  )
	)
	(setq ent (vlax-ename->vla-object (entlast)))
	(vla-put-Alignment ent acAlignmentMiddle)
	(vla-put-textalignmentpoint ent (vlax-3d-point point))
      )
      ;;(command "_insert" "C:\\Users\\sgidi\\Desktop\\papertk.dwg"
	       ;;(car gt_cir) "" "" "" )
      (setq pnt(mapcar '+ (car gt_cir) (list 0 -15.0))
	    pnt(list (car pnt) (cadr pnt) 0.0))
      (txt_entmake "0" 0 1 pnt head)
    )
    ;(draw_tk (list '(0 0) 10) "fasdf")
    
    (defun perform_compare(file	lst	    gt_cir	dr
				/	    i		head
				time	    f_txt	stpnt
				delta	    st_delta	pnt_lst
				width_lst   classic_time
				taubin_time LM_time	classic_param
				taubin_param		LM_param
				classic_str taubin_str	LM_str
			       )
      (setq f_txt(open file "w") i 0 st_delta(list 0 0))
      (if (zerop (- (caddr (car lst)) (caddr (last lst))))
	(setq time nil)
	(setq time T)
	)
      (foreach var lst
	(setq width_lst(list (* (cadr gt_cir) (- 1 (car var))) (* (cadr gt_cir) (+ 1 (car var))))
	      stpnt(mapcar '+ (last var) st_delta)
	      arc_len(/ (*  pi (cadr var)) 180.0))
	(setq pnt_lst (sample_circle_sector (car gt_cir) width_lst (list 0 arc_len) (caddr var)));随机生成点
	(setq head (strcat "弧长," (rtos (cadr var) 2 0) ","
			   "点宽度,"(rtos (car var) 2 5) ","
			   "待点个数," (rtos (caddr var) 2 0)","
			   "起点坐标," (rtos (cadr stpnt) 2 2)","))
	(setq st_time(getvar "MILLISECS")
	      classic_param (circle_fit_one_kpnt pnt_lst stpnt)
	      classic_time(getvar "MILLISECS")
	      taubin_param  (our_circle_fit_one_end pnt_lst stpnt 10)
	      taubin_time(getvar "MILLISECS")
	      LM_param	    (LM_method_one_end pnt_lst stpnt 50)
	      LM_time(getvar "MILLISECS")
	)
	;;计算残差
	(if time
	  (setq	classic_str (strcat (rtos (- classic_time st_time) 2 0) ",")
		taubin_str  (strcat (rtos (- taubin_time classic_time) 2 0)  ",")
		LM_str	    (strcat (rtos (- LM_time taubin_time) 2 0) ",")
	  )
	  (setq	classic_str (cal_mse pnt_lst classic_param gt_cir T)
		taubin_str  (cal_mse pnt_lst taubin_param gt_cir T)
		LM_str	    (cal_mse pnt_lst LM_param gt_cir T)
	  )
	)
	(write-line (strcat head classic_str taubin_str LM_str) f_txt)
	;;计算结果
;;;	(setq result_c (cons (cons arc_len classic_mse) result_c)
;;;	      result_t (cons (cons arc_len taubin_mse) result_t)
;;;	      result_L (cons (cons arc_len LM_mse) result_L)
;;;	      pnt_lst  nil
;;;	)
	;(rem 0 11)
	(if (and dr (zerop (rem i dr)))
	  (progn
	    (draw_pnt pnt_lst)
	    (draw_circle classic_param 5)
	    (draw_circle taubin_param 1)
	    (draw_circle LM_param 3)
	    (draw_circle gt_cir 4)
	    (draw_tk gt_cir head)
	    (setq delta(list (+ 4.0 (cadr gt_cir) (cadr gt_cir)) 0)
		  gt_cir(list (mapcar '+ (car gt_cir) delta) (cadr gt_cir))
		  st_delta (mapcar '+ st_delta delta))
	  )
	)
	(setq i(1+ i))
      )
      (close f_txt)
      )
   
    (defun dup_lst(lst loc others / new_lst new_el)
      (setq new_lst nil)
      (foreach var lst
	(cond ((equal loc 1) (setq new_el(list var (car others) (cadr others) (caddr others))))
	      ((equal loc 2) (setq new_el(list (car others) var (cadr others) (caddr others))))
	      ((equal loc 3) (setq new_el(list (car others) (cadr others) var (caddr others))))
	      ((equal loc 4) (setq new_el(list (car others) (cadr others) (caddr others) var))))
	(setq new_lst(cons new_el new_lst))
	)
      (reverse new_lst)
      )
      ;(dup_lst (list 1 2 3 4 5 6) 1 (list 0.05 22 (list 0 8)))
    ;;不同长度的弧长条件下，各种方法的误差
    (setq gt_cir(list '(0 0) 10)
	  st_pnt'(10 0))
    ;(s(open arc_txt "w")
    ;(close arc_txt)
    (setq arc_txt "弧长0.04.txt"
	  arc_lst(append (range 1 89 1) (range 90 179 5) (range 180 360 10))
	  arc_lst(dup_lst arc_lst 2 (list 0.04 800 st_pnt)))
    (perform_compare arc_txt arc_lst gt_cir nil)
    (setq gt_cir(list '(0 0) 10)
	  st_pnt'(10 0))
    (setq width_txt "\宽度45.txt"
	  width_lst(append (range 0.0001 0.005 0.0005) (range 0.005 0.020 0.001) (range 0.02 0.18 0.002))
	  width_lst(dup_lst width_lst 1 (list 45 800 st_pnt)))
    (perform_compare width_txt width_lst gt_cir nil)
    (setq gt_cir(list '(0 0) 10)
	  st_pnt'(10 0))
    (setq speed_txt "速度.txt"
	  speed_lst(range 100 20000 100)
	  speed_lst(dup_lst speed_lst 3 (list 0.03 45 st_pnt)))
    (perform_compare speed_txt speed_lst gt_cir nil)
    (setq stpnt_txt "起点误差.txt"
	  stpnt_lst(mapcar '(lambda(x) (mapcar '* (list 1 (1+ x)) st_pnt)) (range -0.06 0.06 0.001))  
	  stpnt_lst(dup_lst stpnt_lst 4 (list 0.05 45 800)))
    (perform_compare stpnt_txt stpnt_lst gt_cir nil)
    (draw_compare data st_pnt gt_cir nil)
    
    (setq stpnt_txt "D:\\Program\\circle_fit\\ffff.txt"
	  gt_cir(list '(0 0) 10)
	  st_pnt'(10 0))
    (get_circle_param  stpnt_txt (list 0.02 45  800 st_pnt) gt_cir)
  
    )
   ;;写入圆参数数据到txt中
    (defun get_circle_param (txt	 lst	     gt_cir
			     /		f_txt width_lst   stpnt
			     arc_len	 head	     classic_param
			     taubin_param	     LM_param
			    )
      (defun param_str (param /)
	(strcat	(rtos (car (car param)) 2 5)
		","
		(rtos (cadr (car param)) 2 5)
		","
		(rtos (cadr param) 2 5)
		","
	)
      )
      (setq f_txt(open txt "w"))
      (setq width_lst (list (* (cadr gt_cir) (- 1 (car lst)))
			    (* (cadr gt_cir) (+ 1 (car lst)))
		      )
	    stpnt     (last lst)
	    arc_len   (/ (* pi (cadr lst)) 180.0)
      )
      (setq pnt_lst (sample_circle_sector
		      (car gt_cir)
		      width_lst
		      (list 0 arc_len)
		      (caddr lst)
		    )
      )					;随机生成点
      (foreach var pnt_lst
	(setq
	  line (apply
		 'strcat
		 (mapcar '(lambda (x) (strcat (rtos x 2 3) ",")) var)
	       )
	)
	(write-line (strcat "point," line) f_txt)
      )
      (setq head (strcat "弧长,"
			 (rtos (cadr lst) 2 0)
			 ",\n"
			 "点宽度,"
			 (rtos (car lst) 2 5)
			 ",\n"
			 "待点个数,"
			 (rtos (caddr lst) 2 0)
			 ",\n"
			 "起点坐标,"
			 (rtos (cadr stpnt) 2 2)
			 ","
		 )
      )
      (write-line head f_txt)
      (setq classic_param (circle_fit_one_kpnt pnt_lst stpnt)
	    taubin_param  (our_circle_fit_one_end pnt_lst stpnt 10)
	    LM_param	  (LM_method_one_end pnt_lst stpnt 50)
      )
      (write-line
	(strcat "cparam," (param_str classic_param))
	f_txt
      )
      (write-line
	(strcat "oparam," (param_str taubin_param))
	f_txt
      )
      (write-line (strcat "lparam," (param_str LM_param)) f_txt)
      (close f_txt)
    )
 
  ;;计算残差
  (defun cal_mse(pntlst circle_param gt_circle str / cen rad n mse cen_diff rad_diff)
    (setq cen(car circle_param)
	  rad(cadr circle_param)
	  n(length pntlst)
	  mse(apply '+ (mapcar '(lambda(x) (/ (expt (-(distance x cen) rad) 2) (float n))) pntlst))
	  cen_diff(distance (car gt_circle) cen)
	  rad_diff(- (cadr gt_circle) rad))
    (setq result(list mse cen_diff rad_diff))
    (if str
      (setq result(apply 'strcat (mapcar '(lambda(x) (strcat (rtos x 2 10) ",")) result))))
    result
    )
  ;(setq file(open "E:/CAD_Program/地铁中轴线绘制程序/data/fff.txt" "w"))
  ;(close file)
  ;;写入列表文件
  (defun write_txt(lst / txt f_txt var)
    (setq txt(getfiled "选择要编辑的文件" "D:/" "txt" 0)
	  f_txt (open txt "w"))
    (foreach var lst
      (setq var(apply 'strcat (mapcar '(lambda(x) (strcat (rtos x 2 10) ",")) var)))
      (write-line var f_txt)
      )
    (close f_txt)
    )
  ;;************************************************************************
  ;;**                    start 扇形区域均匀采样方法                      **
  ;;************************************************************************
  ;(draw_pnt (sample_circle_sector (list 0 0) (list 1 1) (list 0 2)  100))
  (defun sample_circle_sector (cen r_sector theta_sector num / r_rand theta_rand)
    (setq r_rand     (random_lst num r_sector T);(abs(apply '- r_sector));
	  theta_rand (random_lst num theta_sector nil);(abs(apply '- theta_sector));
	  p_lst nil
    )
    (mapcar '(lambda (r th)
	       (list (+ (car cen) (* r (cos th)))
		     (+ (cadr cen) (* r (sin th)))
		     ;(+ (cadr cen) (* r (sin th)))
		     0.0
	       )
	     )
	    r_rand
	    theta_rand
    )
  )
  (defun rnd (/ modulus multiplier increment random)
    (if	(not seed)
      (setq seed (getvar "DATE"))
    )
    (setq modulus    65536
	  multiplier 25173
	  increment  13849
	  seed	     (rem (+ (* multiplier seed) increment) modulus)
	  random     (/ seed modulus)
    )
    random
  )
;;;  (defun random_seed( /)
;;;    (* (rem (getvar "cputicks") 1e6) 1e-6)) not false random
  (defun random_lst (n range sqrt_m / ed sum n_lst len num)
    (setq ed t len(abs(apply '- range)))
    (while ed
      (setq n_lst nil
	    sum	0
      )
      (repeat n
	(setq num   (rnd)
	      sum   (+ num sum)
	      n_lst (cons num n_lst)
	)
      )
      (setq ed nil)
      (if (< (abs (- (/ sum (float n)) 0.5)) 0.01)
	(setq ed nil)
      )
    )
    (if sqrt_m
      (setq n_lst(mapcar '(lambda(x) (sqrt (+ (* x (* len (apply '+ range))) (* (car range) (car range))))) n_lst))
      (setq n_lst(mapcar '(lambda(x) (+ (* x len) (car range))) n_lst)))
    n_lst
  )
  ;;************************************************************************
  ;;**                    end   扇形区域均匀采样方法                      **
  ;;************************************************************************

  ;;************************************************************************************
  ;;**             最小化目标为以下函数的圆弧拟合方法（一个端点约束)                  **
  ;;**                                                                                **
  ;;**            f(xc, yc) =  (sqrt((x - xc)^2 + (y - yc)^2) - r)^2                  **
  ;;**                                                                                **
  ;;************************************************************************************

  ;;************************************************************************
  ;;**                         start LM几何拟合                           **
  ;;************************************************************************
  (defun LM_method_one_end(pnt_lst stpnt	max_iter /	  lam
			       mark	i	 n	  avg
			       epsilon	stpnt_av pnt_av	  new_par
			       old_par	jg_m	 lam_mark dx_dy
			       dR	new_cen	 new_r
			      )
    (defun cal_sigma(pntlst circle_par / )
      (apply '+ (mapcar '(lambda(x) (* (- (cadr circle_par) (circle_sqrt x (car circle_par)))
				       (- (cadr circle_par) (circle_sqrt x (car circle_par))))) pntlst))
      )
    (defun circle_sqrt(xy ab )
      (sqrt (apply '+ (mapcar '* (mapcar '- xy ab)(mapcar '- xy ab)))))
    (defun xyztojb(ipnt stpnt param / J_xc J_yc  g_xc g_yc st_sqrt i_sqrt)
      (setq st_sqrt (circle_sqrt stpnt (car param))
	    i_sqrt  (circle_sqrt ipnt (car param))
	    J_xc    (+ (/ (- (car (car param)) (car stpnt)) st_sqrt)
		       (/ (- (car ipnt) (car (car param))) i_sqrt)
		    )
	    J_yc    (+ (/ (- (cadr (car param)) (cadr stpnt)) st_sqrt)
		       (/ (- (cadr ipnt) (cadr (car param))) i_sqrt)
		    )
	    g_xc    (* J_xc (- st_sqrt i_sqrt))
	    g_yc    (* J_yc (- st_sqrt i_sqrt))
      )
      (list (* J_xc J_xc)
	    (* J_xc J_yc)
	    (* J_yc J_yc)
	    g_xc g_yc)
      )
    (defun cal_dxdy(J11 J12 J22 G1 G2 lam / det dx dy)
      (setq det(- (* J11 J22 (1+ lam) (1+ lam))
		  (* J12 J12))
	    dy(/ (-(* G2  J11 (1+ lam))
		   (* G1 J12))
		 det)
	    dx(/ (-(* G1 J22 (1+ lam))
		   (* G2 J12))
		 det))
      (list dx dy)
      )
    (setq mark T i 0
	  lam 0.001
	  epsilon 3.e-8
	  n	 (length pnt_lst)
	  avg	 (mapcar '/
			 (apply 'mapcar (cons '+ pnt_lst))
			 (list (float n) (float n))
		 )
	  stpnt_av(mapcar '- stpnt avg)
	  pnt_av(mapcar '(lambda (a) (mapcar '- a avg)) pnt_lst)
	  ;new_par(circle_fit pnt_lst)
	  ;new_par(list (car new_par) (distance (car new_par) stpnt_av))
	  new_par(list '(0 0) (distance '(0 0) stpnt_av))
    )
    (while (and mark (< i max_iter))
      (setq old_par new_par
	    jg_m (apply
		   'mapcar
		   (cons '+
			 (mapcar
			   '(lambda (a)
			      (xyztojb a stpnt_av old_par
			      )
			    )
			   pnt_av
			 )
		   )
		 );;计算雅可比矩阵
	    jg_m(mapcar '(lambda(x) (/ x n)) jg_m)
	    lam_mark T)
      (while lam_mark
	(setq dx_dy(cal_dxdy (car jg_m) (cadr jg_m) (caddr jg_m) (cadddr jg_m) (last jg_m) lam)
	      dR(sqrt (apply '+ (mapcar '* dx_dy dx_dy))))
	(if (< (/ (apply '+ (mapcar 'abs (cons dR dx_dy))) (1+ (cadr old_par))) epsilon)
	  (setq lam_mark nil mark nil)
	  (progn
	    (setq new_cen(mapcar '- (car old_par) dx_dy)
		  new_r(distance new_cen stpnt_av)
		  new_par(list new_cen new_r)
		  new_sigma(cal_sigma pnt_av new_par)
	    )
	    (if (apply 'or (mapcar '> (mapcar 'abs new_cen) (list 1.e+6 1.e+6)))
	      (setq lam_mark nil mark nil)
	      (if (< new_sigma (cal_sigma pnt_av old_par))
		(setq lam(* lam 0.01) ;;减小步长靠近极值点
		      lam_mark nil mark T)
		(setq lam(* lam 10)   ;;增加步长跳出单调增加区域
		      lam_mark T mark T))
	      )
	    )
	  )
	(setq i(1+ i))
	)
      (setq i(1+ i))
      )
    (setq cen (mapcar '+ (car new_par) avg)
	  radian  (vl-sort
		    (mapcar '(lambda (x) (angle cen x)) (cons stpnt pnt_lst))
		    '(lambda (e1 e2) (< e1 e2))
		  )
    )
    (list cen (cadr new_par) (car radian) (last radian))
    )
  ;;************************************************************************
  ;;**                         end   LM几何拟合                           **
  ;;************************************************************************
  
  ;;************************************************************************************
  ;;**                  最小化目标为以下函数的圆弧拟合方法（一个端点约束)             **
  ;;**                                                                                **
  ;;**                            ((x - xc)^2 + (y - yc)^2 - r^2)^2                   **
  ;;**            f(xc, yc) =     ---------------------------------                   **
  ;;**                                ((x - xc)^2 + (y - yc)^2)                       **
  ;;**                                                                                **
  ;;**   这种方法计算速度也很快无需进行矩阵运算，特征方程只需要4-6次迭代就可以求得    **
  ;;************************************************************************************

  ;;************************************************************************
  ;;**               start       一个端点约束的圆弧拟合                   **
  ;;************************************************************************
  (defun our_circle_fit_one_end(pnt_lst st_pnt max_iter  /	p_pnt   n
			    xy_to_xyz	    avg	    xyz	    M_s
			    M_s	    M_xy    M_xx    M_yy    M_xz
			    M_yz    M_zz    M_z	    cov_xy  var_z
			    A3	    A2	    A1	    A0	    A22
			    A33	    mark    i	    iter_x  iter_y
			    delta_y x_new   y_new   DET	    x_c
			    y_c	    rad	    radian
			   )
    ;;max_iter为特征方程求解最大迭代次数，一般4到6次收敛，可设置10或更大
    (defun one_end_xy_to_xyz (xy avg stpnt / z zpnt xy-st)
      (setq xy	  (mapcar '- xy avg)
	    zpnt    (apply '+ (mapcar '* xy xy))
	    stpnt (mapcar '- stpnt avg)
	    z	  (- (apply '+ (mapcar '* xy xy))
		     (apply '+ (mapcar '* stpnt stpnt))
		  )
	    xy-st (mapcar '- xy stpnt)
      )
      (list (car xy-st) (cadr xy-st) z zpnt)
    )
    (defun square_root (a b c / sqrt_b2_4ac)
      (setq sqrt_b2_4ac (sqrt (- (* b b) (* 4.0 a c))))
      (list (/ (+ b sqrt_b2_4ac) -2.0 a)
	    (/ (- b sqrt_b2_4ac) -2.0 a)
      )
    )
    (defun scale(lst / min_value)
      (setq min_value(apply 'min (mapcar 'abs lst)))
      (mapcar '(lambda(x) (/ x min_value)) lst)
      )
    (setq n	 (length pnt_lst)
	  avg	 (mapcar '/
			 (apply 'mapcar (cons '+ pnt_lst))
			 (list (float n) (float n))
		 )
	  xyz	 (mapcar '(lambda (a) (one_end_xy_to_xyz a avg stpnt)) pnt_lst)
	  M_s	 (apply
		   'mapcar
		   (cons '+
			 (mapcar
			   '(lambda (a)
			      (list (* (car a) (cadr a))
				    ;;mxy
				    (* (car a) (car a))
				    ;;mxx
				    (* (cadr a) (cadr a))
				    ;;myy
				    (* (car a) (caddr a))
				    ;;mxz
				    (* (cadr a) (caddr a))
				    ;;myz
				    (* (caddr a) (caddr a))
				    ;;mzz
			      )
			    )
			   xyz
			 )
		   )
		 )
	  p_pnt  (mapcar '- stpnt avg)
	  p_pnt  (list (car p_pnt) (cadr p_pnt) (apply '+ (mapcar '* p_pnt p_pnt)))
	  ;p_pnt  (mapcar '(lambda(x) (* x n)) p_pnt)
	  ;M_s	 (mapcar '(lambda(x) (/ x n)) M_s)
	  ;M_s    (scale M_s)
	  M_xy	 (car M_s)
	  M_xx	 (cadr M_s)
	  M_yy	 (caddr M_s)
	  M_xz	 (- (cadddr M_s) (* 2.0 (car p_pnt)));(cadddr M_s)
	  M_yz	 (- (nth 4 M_s) (* 2.0 (cadr p_pnt)));(nth 4 M_s);
	  M_zz	 (nth 5 M_s)	  
	  ;N_z    (nth 6 M_s)
	  M_z	 (+ M_xx M_yy)    
	  cov_xy (- (* M_xx M_yy) (* M_xy M_xy))
	  A3	 (* 4.0 (caddr p_pnt))
	  A2	 (- (* -4.0 M_z (caddr p_pnt)) M_zz)
	  A1	 (+ (* M_zz M_z)
		    (* 4.0 cov_xy (caddr p_pnt))
		    (* -1.0 M_xz M_xz)
		    (* -1.0 M_yz M_yz)
		 )
	  A0	 (+ (* M_xz (- (* M_xz M_yy) (* M_yz M_xy)))
		    (* M_yz (- (* M_yz M_xx) (* M_xz M_xy)))
		    (* -1.0 cov_xy M_zz)
		 )
	  A22	 (+ A2 A2)
	  A33	 (+ A3 A3 A3)
	  mark	 T
	  i	 0
	  iter_x 0.0
	  iter_y A0
    )
    (while (and mark (< i max_iter))
      (setq delta_y (+ A1 (* iter_x (+ A22 (* iter_x A33)))))
      (if (< (abs delta_y) 0.00000000001)
	(setq mark nil)
	(progn
	  (setq x_new (- iter_x (/ iter_y delta_y)))
	  (if (< (abs (- x_new iter_x)) 0.00000001)
	    (setq mark nil)
	    (progn
	      (setq y_new
		     (+	A0
			(* x_new (+ A1 (* x_new (+ A2 (* x_new A3)))))
		     )
	      )
	      (if (>= (abs y_new) (abs iter_y))
		;;(< (abs(- y_new iter_y)) 0.00000001)
		(setq mark nil)
		(setq iter_x x_new
		      iter_y y_new
		)
	      )
	    )
	  )
	)
      )
      (setq i (1+ i))
    )
    ;;(setq iter_x(cadr (square_root A33 A22 A1))
    ;;	iter_y (+ A0 (* x_new (+ A2 (* A3 x_new)))))
    (setq DET (+ (* iter_x iter_x)
		 (* -1.0 iter_x M_z)
		 cov_xy
	      )
	  x_c (/ (- (* M_xz (- M_yy iter_x))
		    (* M_yz M_xy)
		 )
		 DET
		 2.0
	      )
	  y_c (/ (- (* M_yz (- M_xx iter_x))
		    (* M_xz M_xy)
		 )
		 DET
		 2.0
	      )
	  x_c (+ x_c (car avg))
	  y_c (+ y_c (cadr avg))
	  rad (distance (list x_c y_c) (list (car stpnt) (cadr stpnt)))
    )
    (setq
      radian (vl-sort
	       (mapcar '(lambda (x) (angle (list x_c y_c) x)) (cons stpnt pnt_lst))
	       '(lambda (e1 e2) (< e1 e2))
	     )
    )
    (list (list x_c y_c) rad (car radian) (last radian))
  )
  ;;************************************************************************ 
  ;;**               end         一个端点约束的圆弧拟合                   **
  ;;************************************************************************

  ;;************************************************************************************
  ;;**                 最小化目标为以下函数的圆弧拟合方法(一个端点约束)               **
  ;;**                                                                                **
  ;;**             f(xc, yc, r) = ((x - xc)^2 + (y - yc)^2 - r^2)^2                   **
  ;;**                                                                                **
  ;;**     这种方法为一种经典代数拟合方法，优点是计算速度快，但是拟合圆弧可能会偏小*****
  ;;************************************************************************************

  ;;************************************************************************
  ;;**                  start 一个端点约束的圆弧拟合                      **
  ;;************************************************************************
  (defun circle_fit_one_kpnt (  pnt_lst	kpnt /     to_xy n	  n_list
			      dif   sum	  a	b     m	    s	  cv
			      cu    cen	  rad	radian
			     )
    (defun to_xy (xy pq / xp yq xy2pq2)
      (setq xp	   (- (car pq) (car xy))
	    yq	   (- (cadr pq) (cadr xy))
	    xy2pq2 (- (apply '+ (mapcar '* xy xy))
		      (apply '+ (mapcar '* pq pq))
		   )
      )
      (list (* xp xp)
	    ;;Suu
	    (* xp yq)
	    ;; Suv
	    (* yq yq)
	    ;;Svv
	    (* -1.0 xp xy2pq2)
	    ;;Suuu
	    (* -1.0 yq xy2pq2)
	    ;;Suvv
      )
    )
    (setq n	 (length pnt_lst)
	  n_list (list (float n) (float n) (float n))
	  avg	 (mapcar '/
			 (apply 'mapcar (cons '+ pnt_lst))
			 n_list
		 )
	  sum	 (apply
		   'mapcar
		   (cons '+
			 (mapcar
			   '(lambda (a)
			      (to_xy (mapcar '- a avg) (mapcar '- kpnt avg))
			    )
			   pnt_lst
			 )
		   )
		 )
	  a	 (list (car sum)
		       (cadr sum)
		       (/ (cadddr sum) 2.)
		 )
	  b	 (list (cadr sum)
		       (caddr sum)
		       (/ (last sum) 2.)
		 )
	  m	 (/ (car a) (car b))
	  s	 (mapcar '- a (mapcar '(lambda (x) (* x m)) b))
	  cv	 (/ (caddr s) (cadr s))
	  cu	 (/ (- (caddr a) (* (cadr a) cv)) (car a))
	  cen	 (mapcar '+ (list cu cv) avg)
	  rad	 (distance cen (list (car kpnt) (cadr kpnt)))
    )
    (setq radian
	   (vl-sort
	     (mapcar '(lambda (x) (angle cen x)) (cons kpnt pnt_lst))
	     '(lambda (e1 e2) (< e1 e2))
	   )
    )
    (list cen rad (car radian) (last radian))
  )
  ;;************************************************************************
  ;;**                  end   一个端点约束的圆弧拟合                      **
  ;;************************************************************************
  
  ;;(circle_fit (list '(12.0 23.0) '(23.0 32.0) '(11.0 32.0) '(13.0 24.0)))
  (defun draw_res (obj pnt_lst h ro threold / i avg_res res)
    ;;计算残差并绘制到点上 h为字高， ro 为旋转
    (setq i	  0
	  obj	  (vlax-ename->vla-object obj)
	  avg_res 0
    )
    (foreach pnt pnt_lst
      (setq pnt	    (list (car pnt) (cadr pnt) 0.0)
	    res	    (distance (vlax-curve-getclosestpointto obj pnt t) pnt)
	    res	    (* 100 res)
	    avg_res (+ avg_res res)
      )
      ;;小于指定阈值变绿色，否则红色
      (if (< res threold)
	(entmakex (list	'(0 . "TEXT")
			'(8 . "Residuals")
			(cons 10 pnt)
			(cons 40 h)
			(cons 50 ro)
			(cons 62 3)
			(cons 1 (rtos res 2 3))
		  )
	)
	(entmakex (list	'(0 . "TEXT")
			'(8 . "Residuals")
			(cons 10 pnt)
			(cons 40 h)
			(cons 50 ro)
			(cons 62 1)
			(cons 1 (rtos res 2 3))
		  )
	)
      )
    )
    (/ avg_res (float (length pnt_lst)))
  )
  (defun draw_pnt (pnt_lst / i avg_res res)
    ;;计算残差并绘制到点上 h为字高， ro 为旋转
    (foreach pnt pnt_lst
      (entmakex (list	'(0 . "POINT")
			'(8 . "Points")
			(cons 10 pnt)
		  )
	))
  )

  (vl-load-com)
  (setq	CMDECHO	(getvar "CMDECHO")
	OSMODE	(getvar "OSMODE")
  )
  (setvar "CMDECHO" 0)
  (setvar "OSMODE" 0)
  (if (not (tblsearch "layer" "Residuals"))
    (command "_layer" "m" "Residuals" "c" 8 "" "")
  )
  (if (not (tblsearch "layer" "fit_line"))
    (command "_layer" "m" "fit_line" "c" 8 "" "")
  )
  (main_fun)
  (princ "拟合完成！")
  (setvar "CMDECHO" CMDECHO)
  (setvar "OSMODE" OSMODE)
)