四参数函数及应用-转点及全图坐标系转换

wzg356

问题的提出，老老帖了，最近觉得这类资料少，算是补充
http://bbs.mjtd.com/forum.php?mo ... =%CB%C4%B2%CE%CA%FD
已知这四参数
纵坐标平移参数X=59.279033
横坐标平移参数Y=59.191856
旋转角度T=0.0001852315256（度）
伸缩尺度K=0.999985766528
可以写一个这样的程序来进行转换吗？因为转出来后几公分偏差，不知道写是不是正确的
(defun c:xzb()
  (setq sicanshu (ssget) )
  (command"move" sicanshu "" "0,0" "59.279033,59.191856" "")
  (command"_scale" sicanshu "" "0,0" "0.999985766528" "")
  (command"_rotate" sicanshu "" "0,0" "0.0001852315256" "")
)

;写个函数验证一下
;四参数获取，公共点对数量需对应
(defun pts24cs (oldps newps / ps1 ps2 p01 p02 ls1 ls2 ag1 ag2 k rt  dx dy)        
        (setq ps1 (mapcar '(lambda (a) (mapcar '+ (list 0.0 0.0)a)) oldps)
                ps2 (mapcar '(lambda (a) (mapcar '+ (list 0.0 0.0)a)) newps)
        );转二维               
        (while(and(setq p01(car ps1)  ps1(cdr ps1))
                        (setq p02(car ps2)  ps2(cdr ps2))
                )
                (setq ls1(cons (mapcar '(lambda(a)(distance p01 a)) ps1) ls1))
                (setq ag1(cons (mapcar '(lambda(a)(angle p01 a)) ps1) ag1))
                (setq ls2(cons (mapcar '(lambda(a)(distance p02 a)) ps2)ls2))
                (setq ag2(cons (mapcar '(lambda(a)(angle p02 a)) ps2)ag2))
        );所有可能的两点组合的距离及方位角
        (setq ls1 (apply 'append ls1) ag1(apply 'append ag1))
        (setq ls2(apply 'append ls2) ag2(apply 'append ag2))
        (setq k(/ (apply '+ (mapcar '/ ls2 ls1)) (length ls1)));平均k
        (setq rt(/ (apply '+(mapcar '- ag2 ag1)) (length ag1)));平均旋转弧度
        (setq oldps (mapcar '(lambda (a)(Rot2D a rt)) oldps));老点旋转后
        (setq oldps (mapcar '(lambda (a)(mapcar '* (list k k)a)) oldps));乘k        
        (setq dxys(mapcar '(lambda(a b)(mapcar '- a b)) newps oldps))
        (setq dx (/ (apply '+(mapcar 'car dxys)) (length (car dxys))));平均dx
        (setq dy (/ (apply '+(mapcar 'cadr dxys)) (length (cadr dxys))));平均dy
        (list dx dy rt k)
)
;;; 旋转向量到指定角度        ;;; 输入: 一个向量和指定的角度
;;; 输出: 被旋转后的向量
(defun Rot2D (v a / c s x y)
        (setq c (cos a) s (sin a))
        (setq x (car v) y (cadr v))
        (list (- (* x c) (* y s)) (+ (* x s) (* y c)))
)

;应用及验证
;曾用过的点对-公共已知点 ， 最少2个点对
(setq oldps (list(list 5087067.129 532083.741)(list 5087050.578 532228.286));老坐标系
      newps (list(list 5087679.709 609542.699)(list 5087664.971 609687.453));新坐标系
)
;现场实测第三点对-公共已知点
(setq oldp (list 5086983.205 532095.119) newp (list 5087595.926 609555.129))
;以上三个点对因保密原因，关键地理信息去掉了
(setq 4cs(pts24cs oldps newps));取得四参数
(setq newp(mapcar '+ (list (car 4cs)(cadr 4cs))(mapcar '* (list (cadddr 4cs)(cadddr 4cs)) (Rot2D oldp (caddr 4cs)))))
;转换第三个点
(list (rtos (car newp)2 5)(rtos (cadr newp)2 5))=>("5087595.92692" "609555.12984")与现场一致
;与其他软件计算结果也基本一致

;验证一下图面操作
(defun c:yz()         
        (setq oldps (list(list 5087067.129 532083.741)(list 5087050.578 532228.286))
      newps (list(list 5087679.709 609542.699)(list 5087664.971 609687.453))
    );坐标点对-公共点，老-新坐标系数据
    (setq 4cs(pts24cs oldps newps));;取得四参数(list dx dy rt k)
    (foreach cs 4cs (princ (strcat "\n"(rtos cs 2 10))))
        (setq en(entmakex(list '(0 . "line") (cons 10 (car oldps)) (cons 11 (cadr oldps)))))
        ;老坐标画的线，黑        
        (command"_rotate" en "" "0,0" (* (/ (caddr 4cs) pi)180));弧度转为度
        (command"_scale" en "" "0,0" (cadddr 4cs))
        (command"move" en "" "0,0" (list (car 4cs)(cadr 4cs)))
        (entmakex(list '(0 . "line") '(62 . 1)(cons 10 (car newps)) (cons 11 (cadr newps))))
        ;新坐标画的线，红
        ;两线重合        
)
;说明依据四参数，可以整图转换的
；需要注意的是，转换图形的时候，command “move”是最后一步

zhaoxt

很棒，感谢分享，学习一下