单像空间后方交会编程实现

文中有详细注释，根据《摄影测量学（第三版）》（武汉大学出版）进行编写

import numpy as np
from numpy import *
from math import*
from numpy.linalg import inv
print('坐标数据并保存为CSV格式\n')   # 读取数据为矩阵形式
original_data1 = np.loadtxt(open('Coordinat_Date.csv'), delimiter=",", skiprows=0)
# 输入计算参数，以数组的形式
print("请依次输入比例尺，主距，x0,y0,迭代次数，均以逗号隔开")
A1 = input(" 比例尺，主距，x0，y0，迭代次数: ").split()
B = list(map(int, A1))
print(B)
m = B[0]
f1 = B[1]
x01, y01 = B[2], B[3]
num1 = B[4]
# 建立后方交会的类，用于面向对象编程
class Resection():# 初始化参数def __init__(self):self.original_data = original_data1
# 读取数据为矩阵形式self.f = f1self.x0, y0 = x01, y01self.num = num1self.xy = []self.XYZ = []self.fi, self.w, self.k = 0, 0, 0
# 《摄影测量学》P.78 一般相片倾角小于3°所以外方位元素近似取φ，ω，κ=0,书中原话
# 读取影像坐标，存到xy列表，相应地面点坐标存到XYZ列表self.xuanzhuan = []self.Xs0 = []self.Ys0 = []self.Zs0 = []self.rotate = []self.diedai = 0
# 建立转换方法函数     # 《摄影测量学》P.74def coordinate(self, m):for i in range(len(self.original_data)):self.xy.append([self.original_data[i][0]/1000, self.original_data[i][1]/1000])self.XYZ.append([self.original_data[i][2], self.original_data[i][3], self.original_data[i][4]])
# 定义系数矩阵A，常数项矩阵LA = np.mat(np.zeros((len(self.xy*2), 6)))L = np.mat(np.zeros((len(self.xy*2), 1)))
# 将xy和XYZ列表转化为矩阵self.xy = np.mat(self.xy)self.XYZ = np.mat(self.XYZ)XYZ_CHA = np.mat(np.zeros((len(self.xy), 3)))    # 便于推到偏导数建立的矩阵sum_X = 0sum_Y = 0
# Xs0 Ys0 取四个角上控制点坐标的平均值   Zs0=H=mf，# 《摄影测量学》P.78上方有详细说明for i in range(len(self.original_data)):sum_X = self.original_data[i][2]+sum_Xsum_Y = self.original_data[i][3]+sum_Yself.Xs0 = 0.25*sum_Xself.Ys0 = 0.25*sum_Yself.Zs0 = m * self.fwhile(self.diedai<self.num):
#旋转矩阵a1 = cos(self.fi)*cos(self.k)-sin(self.fi)*sin(self.w)*sin(self.k)a2 = (-1.0) * cos(self.fi) * sin(self.k) - sin(self.fi) * sin(self.w) * cos(self.k)a3 = (-1.0) * sin(self.fi) * cos(self.w)b1 = cos(self.w) * sin(self.k)b2 = cos(self.w) * cos(self.k)b3 = (-1.0) * sin(self.w)c1 = sin(self.fi) * cos(self.k) + cos(self.fi) * sin(self.w) * sin(self.k)c2 = (-1.0) * sin(self.fi) * sin(self.k) + cos(self.fi) * sin(self.w) * cos(self.k)c3 = cos(self.fi) * cos(self.w)self.xuanzhuan = np.mat([[a1, a2, a3], [b1, b2, b3], [c1, c2, c3]])for i in range(len(self.XYZ)):XYZ_CHA[i, 0] = self.XYZ[i, 0]-self.Xs0XYZ_CHA[i, 1] = self.XYZ[i, 1]-self.Ys0XYZ_CHA[i, 2] = self.XYZ[i, 2]-self.Zs0self.XYZ_ = self.xuanzhuan.T*XYZ_CHA.Tfor i in range(len(self.XYZ)):
# 系数阵：A[i*2, 0] = -self.f/(self.Zs0-self.XYZ[i, 2])A[i*2, 1] = 0A[i*2, 2] = -self.xy[i, 0]/(self.Zs0-self.XYZ[i, 2])A[i*2, 3] = -self.f*(1+pow(self.xy[i, 0], 2)/pow(self.f, 2))A[i*2, 4] = -(self.xy[i, 0]*self.xy[i, 1])/self.fA[i*2, 5] = self.xy[i, 1]A[i*2+1, 0] = 0A[i*2+1, 1] = -self.f/(self.Zs0 - self.XYZ[i, 2])A[i*2+1, 2] = -self.xy[i, 1]/(self.Zs0 - self.XYZ[i, 2])A[i*2+1, 3] = -(self.xy[i, 0]*self.xy[i, 1])/self.fA[i*2+1, 4] = -self.f*(1+pow(self.xy[i, 1], 2)/pow(self.f, 2))A[i*2+1, 5] = -self.xy[i, 0]# 常数项：L[i * 2, 0] = self.xy[i, 0]+self.f*(self.XYZ_[0, i]/self.XYZ_[2, i])L[i * 2 + 1, 0] = self.xy[i, 1]+self.f*(self.XYZ_[1, i]/self.XYZ_[2, i])# 结果：Result = ((A.T * A).I)*A.T*L  # 《摄影测量学》P.74-5-6self.Xs0 += Result[0]self.Ys0 += Result[1]self.Zs0 += Result[2]self.fi += Result[3]self.w += Result[4]self.k += Result[5]self.diedai = self.diedai+1# 像点坐标的改正值V = A * Result - Lsigma0 = sqrt((V.T*V)/(2*4-6)) print("单位权中误差：")print(sigma0)Q = inv(A.T * A)for i in range(1,6):D = Q[i, i]# 外方位元素的精度sigma = sqrt(D) * sigma0print("对角线上外方位元素精度：")print(sigma)a1 = cos(self.fi)*cos(self.k)-sin(self.fi)*sin(self.w)*sin(self.k)a2 = (-1.0) * cos(self.fi) * sin(self.k) - sin(self.fi) * sin(self.w) * cos(self.k)a3 = (-1.0) * sin(self.fi) * cos(self.w)b1 = cos(self.w) * sin(self.k)b2 = cos(self.w) * cos(self.k)b3 = (-1.0) * sin(self.w)c1 = sin(self.fi) * cos(self.k) + cos(self.fi) * sin(self.w) * sin(self.k)c2 = (-1.0) * sin(self.fi) * sin(self.k) + cos(self.fi) * sin(self.w) * cos(self.k)c3 = cos(self.fi) * cos(self.w)self.rotate = np.mat([[a1, a2, a3], [b1, b2, b3], [c1, c2, c3]])# 实例化对象
T = Resection()
t = T.coordinate(m)    # 以比例尺作为变量代入，这里的变量代入是为了调用函数而已，变量也可以是其他的参数值
print('原始数据如下(x,y,X,Y,Z):\n', original_data1)
print('计算结果\n', T.Xs0, '\n', T.Ys0, '\n', T.Zs0, '\n')
print('旋转矩阵\n', T.rotate)
print('迭代次数为：', T.diedai)

程序中的Coordinat_Date.csv格式如下：

-86.15,-68.99,36589.41,25273.32,2195.17
-53.4,82.21,37.63108,31324.51,728.69
-14.78,-76.63,39100.97,24934.98,2396.5
10.46,64.43,40426.54,30319.81,757.31

可以根据自己的需求使用不同的检验代码：表中数据每一行都以x,y,X,Y,Z的顺序排列

坐标数据并保存为CSV格式请依次输入比例尺，主距，x0,y0,迭代次数，均以逗号隔开比例尺，主距，x0，y0，迭代次数: 1200 1 1 1 20
[1200, 1, 1, 1, 20]
单位权中误差：
0.03808031015889729
对角线上外方位元素精度：
13.002322081468744
对角线上外方位元素精度：
205.3896722827482
对角线上外方位元素精度：
0.0346622511350548
对角线上外方位元素精度：
0.04088792697435863
对角线上外方位元素精度：
0.5470752369098285
原始数据如下(x,y,X,Y,Z):[[-8.615000e+01 -6.899000e+01  3.658941e+04  2.527332e+04  2.195170e+03][-5.340000e+01  8.221000e+01  3.763108e+01  3.132451e+04  7.286900e+02][-1.478000e+01 -7.663000e+01  3.910097e+04  2.493498e+04  2.396500e+03][ 1.046000e+01  6.443000e+01  4.042654e+04  3.031981e+04  7.573100e+02]]
计算结果[[5814720.85507374]] [[-6723902.61979151]] [[791.63275503]] 旋转矩阵[[-0.94997338  0.30254713  0.07756161][-0.06446153  0.05306299 -0.99650842][-0.30560642 -0.95165621 -0.03090578]]
迭代次数为： 20Process finished with exit code 0

有不好的地方，还请各位看官指正批评！

单像空间后方交会编程实现相关推荐

双象空间前方交会代码_单像空间后方交会和双像解析空间后方-前方交会的算法程序实现...
单像空间后方交会和双像解析空间后方 - 前方交会的算法程序实现遥感科学与技术摘要:如果已知每张像片的 6 个外方位元素,就能确定被摄物体与航摄像片的关系.因此, 利用单像空间后方交会的方法,可以 ...
单像空间后方交会（C语言）
单像空间后方交会(C语言) 1 原理介绍 1.1 定义 1.2 基本思想 1.3 详细计算 1.4 精度评定 2 问题求解 2.1 问题重述 2.2 问题解读与说明 2.3 c语言求解实现代码 2.4 ...
解析摄影测量之单像空间后方交会（MATLAB）
目录一.题目二.理论基础三.MATLAB代码一.题目二.理论基础三.MATLAB代码 clc;clear; %输入初值 %像点坐标,单位统一化,以米为单位 imgPt_X=[-86.15, ...
单像后方交会、pnp问题迭代计算的数学原理
先提出几个问题: 1.为什么后方交会要迭代法? 2.这个求"改正数"的迭代法怎么保证收敛? 3.这个迭代法的精度分析? 4.单像后方交会与PNP问题有什么联系? 参考<数值分 ...
读芯术python教程答案_攻略Python的免费书单：走进编程，从这五本书开始
全文共1245字,预计学习时长5分钟图源:unsplash Python一向是数据科学家最青睐的编程语言,它的语法相对简单.易于学习.除了机器学习数据库之外,还有非常活跃的开发人员社区,维护着各种库 ...
用php实现一个简易的web表单生成器,网络编程PHP Web表单生成器案例分析
本文实例讲述了PHP Web表单生成器.分享给大家供大家参考,具体如下: 1.实例: 2. 需求分析在项目的实际开发中,经常需要设计各种各样表单.直接编写HTML表单虽然简单,但修改.维护相对麻烦. ...
核心编程第五版配套代码_攻略Python的免费书单：走进编程，从这五本书开始...
全文共1245字,预计学习时长5分钟图源:unsplash Python一向是数据科学家最青睐的编程语言,它的语法相对简单.易于学习.除了机器学习数据库之外,还有非常活跃的开发人员社区,维护着各种库 ...
北理工python慕课10次测验的单选题和编程题答案_20春-程序设计及应用(Python)-何俊-2_章节测验,期末考试,慕课答案查询公众号...
[判断题]液控单向阀与普通单向阀最大的区别是在控制油的作用下能够反向流动. [填空题]She learned the rules of the game step _____ step. [单选题]I ...
python---简单的图形编程
Python自带一套简单的图形开发工具 Turtle 小乌龟用来画一些简单的二维图形通过写代码来画画 import turtle 导入turtle工具 showturtle() 显示箭头指示刚开 ...
Matlab解算空间后方交会外方位元素
目录 1 问题描述 2 实现部分参考文献本问题为武汉大学摄影测量学教材课后习题,现在用MATLAB实现,供大家学习参考. 1 问题描述已知摄像机主距f=153.24mm,四对点的像点坐标 ...

单像空间后方交会编程实现

单像空间后方交会编程实现相关推荐

最新文章

热门文章