近景摄影测量空间后方交会python

首先得有像点坐标和对应物方点坐标，顺序和数量没要求

与航空摄影不同，这里的坐标系需要转换，故读进来的XYZ有变化。

误差方程式：V=At+CX+DX-L

import numpy as np
from math import *# 读数据方法2
s = np.loadtxt('right.txt')
# 记得手动删除像点坐标文件里的重复点，
# 懒得写了
q = np.loadtxt('控制点坐标.txt')
# 点号,x,y
Pxy = []
for i in range(len(s)):if s[i][0] > 100:Pxy.append([s[i][0], s[i][1], s[i][2]])
print(len(Pxy))# X,Y,Z，与Pxy相对应
ground_list = []
# 用来存点号的
temp = [0]*len(q)
for i in range(len(q)):temp[i] = q[i][0]
temp_2 = []
for i in range(len(Pxy)):position = temp.index(Pxy[i][0])ground_list.append([q[position][2], q[position][3], -q[position][1]])
# print(ground_list)# x,y
image_list = []
for i in range(len(Pxy)):image_list.append([Pxy[i][1], Pxy[i][2]])# 初值设定
# 内外方位元素(9)+畸变系数(4)
phi = 0.0
omega = 0.0
kappa = 0.0
Xs = 4500
Ys = -150
Zs = -800
# 单位mm
f = 36.0
x0 = 0.0
y0 = 0.0
k1 = 0.0
k2 = 0.0
p1 = 0.0
p2 = 0.0# 迭代次数
num_ite = 0# 方程内的矩阵
rotate = np.mat(np.zeros((3, 3)))#旋转矩阵
A = np.mat(np.zeros((len(image_list)*2, 9+4)))
# LL:V=(AX+CX2+DXad-L)=(AX-L),L
L = np.mat(np.zeros((len(image_list)*2,1)))
# 存改正数的
delta = [0]*(9+4)
# 存中误差
m = [0]*(9+4)# 开始处理
while (True):# 旋转矩阵a1 = cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa)a2 = (-1.0) * cos(phi) * sin(kappa) - sin(phi) * sin(omega) * cos(kappa)a3 = (-1.0) * sin(phi) * cos(omega)b1 = cos(omega) * sin(kappa)b2 = cos(omega) * cos(kappa)b3 = (-1.0) * sin(omega)c1 = sin(phi) * cos(kappa) + cos(phi) * sin(omega) * sin(kappa)c2 = (-1.0) * sin(phi) * sin(kappa) + cos(phi) * sin(omega) * cos(kappa)c3 = cos(phi) * cos(omega)for i in range(0, len(image_list)):# 单位转换(像素->毫米)x = image_list[i][0]*5.1966/1000y = image_list[i][1]*5.1966/1000X = ground_list[i][0]Y = ground_list[i][1]Z = ground_list[i][2]rr = (x-x0)*(x-x0)+(y-y0)*(y-y0)Xp = a1*(X-Xs) + b1*(Y-Ys) + c1*(Z-Zs);Yp = a2*(X-Xs) + b2*(Y-Ys) + c2*(Z-Zs);Zp = a3*(X-Xs) + b3*(Y-Ys) + c3*(Z-Zs);A[i * 2,0] = 1.0*(a1*f + a3*(x-x0))/ZpA[i * 2,1] = 1.0*(b1*f + b3*(x-x0))/ZpA[i * 2,2] = 1.0*(c1*f + c3*(x-x0))/ZpA[i * 2,3] = (y-y0)*sin(omega) - ((x-x0)*((x-x0)*cos(kappa) - (y-y0)*sin(kappa))/f + f*cos(kappa))*cos(omega)A[i * 2,4] = -f*sin(kappa) - (x-x0)*((x-x0)*sin(kappa) + (y-y0)*cos(kappa))/fA[i * 2,5] = y-y0A[i * 2,6] = (x-x0)/fA[i * 2,7] = 1A[i * 2,8] = 0A[i * 2,9] = -(x-x0)*rrA[i * 2,10] = -(x-x0)*rr*rrA[i * 2,11] = -rr-2*(x-x0)*(x-x0)A[i * 2,12] = -2*(x-x0)*(y-y0)A[i * 2 + 1,0] = 1.0*(a2*f + a3*(y-y0))/ZpA[i * 2 + 1,1] = 1.0*(b2*f + b3*(y-y0))/ZpA[i * 2 + 1,2] = 1.0*(c2*f + c3*(y-y0))/ZpA[i * 2 + 1,3] = -(x-x0)*sin(omega) - ((y-y0)*((x-x0)*cos(kappa)-(y-y0)*sin(kappa))/f-f*sin(kappa))*cos(omega)A[i * 2 + 1,4] = -f*cos(kappa) - (y-y0)*((x-x0)*sin(kappa) + (y-y0)*cos(kappa))/fA[i * 2 + 1,5] = -(x-x0)A[i * 2 + 1,6] = (y-y0)/fA[i * 2 + 1,7] = 0A[i * 2 + 1,8] = 1A[i * 2 + 1,9] = -(y-y0)*rrA[i * 2 + 1,10] = -(y-y0)*rr*rrA[i * 2 + 1,11] = -2*(x-x0)*(y-y0)A[i * 2 + 1,12] = -rr-2*(y-y0)*(y-y0)detax = (x-x0)*(k1*rr+k2*rr*rr)+p1*(rr+2*(x-x0)*(x-x0))+2*p2*(x-x0)*(y-y0)detay = (y-y0)*(k1*rr+k2*rr*rr)+p2*(rr+2*(y-y0)*(y-y0))+2*p1*(x-x0)*(y-y0)L[i*2, 0] = x + f*Xp/Zp-x0+detaxL[i*2+1, 0] = y + f*Yp/Zp-y0+detayAT = A.T #转置矩阵ATA = AT*A#矩阵相乘ATAInv = ATA.I#求逆矩阵ATAAT = ATAInv*ATdelta = ATAAT*L# 添加改正数Xs += delta[0]Ys += delta[1]Zs += delta[2]phi += delta[3]omega += delta[4]kappa += delta[5]f += delta[6]x0 += delta[7]y0 += delta[8]k1 += delta[9]k2 += delta[10]p1 += delta[11]p2 += delta[12]num_ite+=1# 限差(看具体情况设置限差)if (fabs(delta[3]) < 1e-6) and (fabs(delta[4]) < 1e-6) and (fabs(delta[5]) < 1e-6):breakif num_ite>60:print("Error:overtime")break'''
后续处理，各类残差、中误差
'''print("迭代次数:%f\n"%num_ite)
rotate[0,0] = cos(phi)*cos(kappa) - sin(phi)*sin(omega)*sin(kappa)
rotate[0,1] = (-1.0)*cos(phi)*sin(kappa) - sin(phi)*sin(omega)*cos(kappa)
rotate[0,2] = (-1.0)*sin(phi)*cos(omega)
rotate[1,0] = cos(omega)*sin(kappa)
rotate[1,1] = cos(omega)*cos(kappa)
rotate[1,2] = (-1.0)*sin(omega)
rotate[2,0] = sin(phi)*cos(kappa) + cos(phi)*sin(omega)*sin(kappa)
rotate[2,1] = (-1.0)*sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa)
rotate[2,2] = cos(phi)*cos(omega)AX = A*delta
# 残差
v = AX - Lvv = 0
for it in v:vv += it[0,0]*it[0,0]
m0 = sqrt(vv/(2*len(Pxy)-13))print('---------------------')
# 存残差的
print_v = np.zeros((len(Pxy), 3))
for i in range(0,len(Pxy)):print_v[i][0] = Pxy[i][0]print_v[i][1] = v[2*i][0]print_v[i][2] = v[2*i+1][0]# print("%d   %f   %f"%(Pxy[i][0], v[2*i][0], v[2*i+1][0]))
# 可以排个序来挑出较大残差，这里没有
for i in range(0, len(Pxy)):print("点号:%d  vx:%f  vy:%f\n"%(print_v[i][0], print_v[i][1], print_v[i][2]))
print('---------------------')print("单位权中误差m0=%f\n"%m0)
# print(ATAInv)
for n in range(9+4):m[n] = m0 * sqrt(abs(ATAInv[n,n]))print("\n结果：\n")
print("Xs=%f  m=%f"%(Xs, m[0])+"  单位：mm\n")
print("Ys=%f  m=%f"%(Ys, m[1])+"  单位：mm\n")
print("Zs=%f   m=%f"% (Zs, m[2])+"  单位：mm\n")
print("phi=%f  m=%f"% (phi, m[3])+"  单位：弧度\n")
print("omega=%f    m=%f"% (omega, m[4])+"  单位：弧度\n")
print("kappa=%f    m=%f"%(kappa, m[5])+"  单位：弧度\n")
print("f=%f    m=%f"%(f, m[6])+"  单位：mm\n")
print("x0=%f    m=%f"%(x0, m[7])+"  单位：mm\n")
print("y0=%f    m=%f"%(y0, m[8])+"  单位：mm\n")
print("k1=%e    m=%e\n"%(k1, m[9]))
print("k2=%e    m=%e\n"%(k2, m[10]))
print("p1=%e    m=%e\n"%(p1, m[11]))
print("p2=%e    m=%e\n"%(p2, m[12]))
print("\n旋转矩阵R为：\n")
print(rotate)# 写文件
with open('后方交会结果.txt','w') as fo:fo.writelines("\n像点残差：\n")for i in range(0, len(Pxy)):fo.writelines("点号:%d  vx:%f  vy:%f\n"%(print_v[i][0], print_v[i][1], print_v[i][2]))fo.writelines("   \n")fo.writelines("单位权中误差m0=%f\n"%m0)fo.writelines("          \n")fo.writelines("Xs=%f  m=%f"%(Xs, m[0])+"  单位：mm\n")fo.writelines("Ys=%f  m=%f"%(Ys, m[1])+"  单位：mm\n")fo.writelines("Zs=%f   m=%f"% (Zs, m[2])+"  单位：mm\n")fo.writelines("phi=%f  m=%f"% (phi, m[3])+"  单位：弧度\n")fo.writelines("omega=%f    m=%f"% (omega, m[4])+"  单位：弧度\n")fo.writelines("kappa=%f    m=%f"%(kappa, m[5])+"  单位：弧度\n")fo.writelines("   \n")fo.writelines("f=%f    m=%f"%(f, m[6])+"  单位：mm\n")fo.writelines("x0=%f    m=%f"%(x0, m[7])+"  单位：mm\n")fo.writelines("y0=%f    m=%f"%(y0, m[8])+"  单位：mm\n")fo.writelines("     \n")fo.writelines("k1=%e    m=%e\n"%(k1, m[9]))fo.writelines("k2=%e    m=%e\n"%(k2, m[10]))fo.writelines("p1=%e    m=%e\n"%(p1, m[11]))fo.writelines("p2=%e    m=%e\n"%(p2, m[12]))

像点坐标：点号，x，y

控制点坐标：点号，xyz

结果：

近景摄影测量空间后方交会python相关推荐

近景摄影测量 matlab,近景摄影测量课件.ppt
工业与工程摄影测量学(1) 测绘学院徐进军航空航天研究所4-308 学完本课后应该: 掌握近景摄影测量的特点(精度和适用场合) 掌握近景摄影测量实际作业方法和过程掌握解析近景摄影测量的数据处理方 ...
单片空间后方交会程序设计（代码共享）
单片空间后方交会程序设计 1 目的用程序设计语言(VC或者VB)编写一个完整的单片空间后方交会程序,通过对提供的试验数据进行计算,输出像片的外方位元素并评定精度.本实验的目的在于让学生深入理解单片空 ...
空间后方交会前方交会 MFC实现 CSU 摄影测量学
空间后方交会前方交会 MFC 实现 CSU 摄影测量学空间后方交会前方交会 MFC 实现 CSU 摄影测量学摄影测量学基础一. 实验目的二.实验内容与要求三.设计与实现: 3.1设计思路 3 ...
单向空间后方交会C++代码实现
摄影测量作业空间后方交会的代码各位小伙伴要用的话只需修改main函数里边的像点地面点坐标,摄影比例尺分母,焦距等这些已知元素就行. #include<iostream> #inclu ...
空间后方交会c++程序和matlab（可直接运行）
本文不详细说明空间后方交会的原理,只着重说明空间后方交会的程序,并附带一个样例. 样例来源:<摄影测量学>(第二版)武汉大学出版社,张剑清,潘励,王树根. 空间后方交会的误差方程式: 可以 ...
python爬虫新闻网页的浏览量转载量,Python爬取新闻网标题、日期、点击量
最近接触Python爬虫,以爬取学校新闻网新闻标题.日期.点击量为例,记录一下工作进度目前,感觉Python爬虫的过程无非两步: Step1.获取网页url(利用Python库函数import ur ...
使用几何光学实现空间相对定位(python+opencv)
我从2019年3月份开始学习python,在有一定的基础后,我看到学校有一个物理实验竞赛:北京联合大学第十二届物理实验竞赛,其中有一个题目是空间定位,即利用物理原理,自行搭建实验装置,实现物体的空间定 ...
利用卡口数据绘制路段基本图（出入量法）——Python交通数据分析
目录一.卡口数据简介二.路段基本图一.卡口数据简介数据名称含义 DEVICEID 设备ID TRAVELID 车辆ID,可分辨车辆 hpzl 车辆类型,1表示大型车,2表示小汽车 SJ 时间 ...
python广告刷量_用python实现刷点击率的示例代码
背景同事的老爸参加微信的一个活动,需要刷点击率,因此,写了一个程序助之. 准备微信活动也是有真实地址的. 通过mitmproxy(man in the middle proxy)的方式,可以获取微 ...
python qq空间上传_QQ空间的Python接口
QQ空间说说接口这是一个可以用来访问QQ空间说说详细信息的Python模块,能够为用户解析出有用的信息. 用法首先要通过传入cookies创建一个Qzone对象,其次调用它的emotion_lis ...

近景摄影测量空间后方交会python

近景摄影测量空间后方交会python相关推荐

最新文章

热门文章