摄影测量学 ：基于python的双像空间前方交会实验

 输入数据：

1. 模拟像片一对：左片号31，右片号32；
2. 航摄机主距：f=150mm；左右像片x0=0；y0=0。

左片外方位元素：

φ=-0°25'00″                              ω=-1°00'00″                                          k=-0°10'00″

Xs=103007.006117        Ys=139998.994849           Zs=4801.9989994  （m）

右片外方为元素：

φ=1°39'59″                              ω=-0°10'00″                                          k=0°40'00″

Xs=106002.023762        Ys=140005.002780            Zs=4797.009648  （m）

 

已知像点坐标如下表：

25号片			26号片
点号	X（mm）	Y(mm)	X(mm)	Y(mm)
1	0.864	2.652	-100.251	1.443
2	97.084	2.933	-4.425	0.320
3	1.082	-76.227	-102.695	-79.618
4	0.640	82.003	-99.904	81.754
5	95.207	-75.521	-5.334	-78.443
6	96.797	83.77	-3.509	79.587

• 实验步骤：

（1）读入同名点坐标；

（2）读入内外方位元素；

（3）在利用相应的外定向角元素组成旋转矩阵，将像空间坐标转化为像空辅助坐标（X，Y，Z），公式如下：

（4）再计算点投影系数利用公式：

（5）计算模型坐标，公式：

完整代码：

# #!/usr/bin/env python
# # -*- coding: UTF-8 -*-
# '''
# @Project ：Calculate_basin
# @File    ：Ex2.py
# @Author  ：Ryo
# @Date    ：2021/11/15 8:27
# '''
import numpy as np
import pandas as pd
import math

#Define a Rotation Array
def Rotation(phi,w,k):
    R=np.array([[1,-k,-phi],[k,1,-w],[phi,w,1]])
    return R
#import main parameters
f=150
# XS1,YS1,ZS1
leftlin=np.array([103007.006117,139998.994849,4801.9989994])
rightlin=np.array([106002.023762,140005.002780,4797.009648])

leftele=np.array([-25,-60,-10])/60*math.pi/180
rightele=np.array([60+39+59/60,-10,40])/60*math.pi/180
# 左片外方位角元素
# -0°25′00″,-1°00′00″,-0°10′00″
# 右片外方位角元素
# 1°39′59,-0°10′00″,0°40′00″
dataarr=np.array([[0.864,2.652,-100.251,1.443],
[97.084,2.933,-4.425,0.32],
[1.082,-76.227,-102.695,-79.618],
[0.64,82.003,-99.901,81.754],
[95.207,-75.521,-5.334,-78.443],
[96.797,83.77,-3.509,79.587]])


# B(Bx,By,Bz)
B=rightlin-leftlin
R1=Rotation(leftele[0],leftele[1],leftele[2])
R2 = Rotation(rightele[0], rightele[1], rightele[2])

lst=[('Name','X',"Y","Z")]
for i in range(dataarr.__len__()):
    # xy1=np.array(df.loc[i,['25X'	,'25Y']])
    # xy2=np.array(df.loc[i,['26X'	,'26Y']])
    xy1=dataarr[i][0:2]
    xy2=dataarr[i][2:4]
    print('Left (x,y)=',xy1,'right (x,y)=',xy2)

    # a1(X1,Y1,Z1) , a2(X2,Y2,Z2)
    xyz1=np.array(list(xy1)+[-f])
    xyz2=np.array(list(xy2)+[-f])

    a1=np.dot(R1,xyz1)
    a2=np.dot(R2,xyz2)
    print('left (X,Y,Z)=',a1,'right (X,Y,Z)=',a2)
    N1=(B[0]*a2[2]-B[2]*a2[0])/(a1[0]*a2[2]-a1[2]*a2[0])
    N2=(B[0]*a1[2]-B[2]*a1[0])/(a1[0]*a2[2]-a1[2]*a2[0])

    X=leftlin[0]+N1*a1[0]
    Z=leftlin[2]+N1*a1[2]
    Y=0.5*((leftlin[1]+N1*a1[1])+(rightlin[1]+N2*a2[1]))

    lst.append((i,X,Y,Z))
    print('name:',i,'  X=',X,'  Y=',Y,' Z=',Z,'\n')

'''Create GUI'''

from tkinter import *


class Table:

    def __init__(self, root):

        # code for creating table
        for i in range(total_rows):
            for j in range(total_columns):
                self.e = Entry(root, width=20, fg='black',
                               font=('Times New Roman', 20))

                self.e.grid(row=i, column=j)
                self.e.insert(END, lst[i][j])


#Create GUI
from tkinter import ttk
from tkinter import *

# find total number of rows and
# columns in list
total_rows = len(lst)
total_columns = len(lst[0])

# create root window
root = Tk()
t = Table(root)
root.mainloop()

root = Tk()  # 初始框的声明

# Table2
# [X,Y,Z]=inputlist[2:]
columns = ('x_left','y_left','x_right','y_right')
treeview2 = ttk.Treeview(root, height=18, show="headings", columns=columns)  # 表格

for i in columns:
    treeview2.column(i,width=75,anchor='center')
    treeview2.heading(i,text=i)


treeview2.pack(side=LEFT, fill=BOTH)

c1 = [1,2,3,4]

ipcode = [xy1[0], xy1[1],xy2[0],xy2[1]]
for i in range(4):  # 写入数据
    treeview2.insert('', i, values=(xy1[0], xy1[1],xy2[0],xy2[1]))



#text of Rotaion

text = Text(root, width=40, height=20)
# insert的第一个参数为索引;第二个为添加的内容
text.insert(1.,'Left_Rotation is\n')
text.insert(2.0, R1,'\n')
text.insert(7.,'Right_Rotation is=\n')
text.insert(8.0, R2)
# text.insert(3,'fgjsdfkjgierjgidfkgjol')
text.pack()

root.mainloop()  # 进入消息循环