空间后交.txt

很好的双像后方交会

    GetIO CDg1.FileName, a0, a1, a2, b0, b1, b2
    '改化左片像点坐标
    For i = 1 To n
        tempX = Xr(i): tempY = Yr(i)
        Xr(i) = a0 + a1 * tempX + a2 * tempY
        Yr(i) = b0 + b1 * tempX + b2 * tempY
        'txtShow.Text = txtShow.Text & vbCrLf & x1(i) & "  " & y1(i) & "  " & "x=" & Str(Format(x(i), "0.000000")) & " , y=" & Str(Format(y(i), "0.000000"))
    Next i
End Sub

'保存前方交会计算结果的过程
Private Sub mnuSave_Click()
    CDg1.Filter = "Text Files(*.TXT)|*.txt|All Files(*.*)|*.*"
    CDg1.DialogTitle = "保存计算结果"
    CDg1.FileName = ""
    CDg1.Action = 2
    If CDg1.FileName = "" Then Exit Sub
    Open CDg1.FileName For Output As #1
        Write #1, txtShow.Text
    Close #1
    'Shell "C:\WINdows\NOTEPAD.EXE " & CDg1.FileName, vbNormalFocus
End Sub

'空间前方交会计算过程space intersection
Private Sub mnuSpcIntersec_Click()
    '求方向余弦（旋转矩阵）和基线分量
    subGetR fai_L, omg_L, kap_L, R_L
    subGetR fai_R, omg_R, kap_R, R_R
    '求基线分量
    Bx = XsR - XsL: By = YsR - YsL: Bz = ZsR - ZsL
    
    '逐点进行空间前方交会
    Dim i#, X_L#, Y_L#, Z_L#, X_R#, Y_R#, Z_R#      '循环变量，像空辅助坐标
    Dim NL#, NR#                    '点投影系数
    ReDim X(1 To n), Y(1 To n), Z(1 To n)
    txtShow.Text = txtShow.Text & "空间前方交会计算结果：" & vbCrLf
    For i = 1 To n
        '求像空辅助坐标
        X_L = R_L(1, 1) * Xl(i) + R_L(1, 2) * Yl(i) - R_L(1, 3) * f  '左片
        Y_L = R_L(2, 1) * Xl(i) + R_L(2, 2) * Yl(i) - R_L(2, 3) * f
        Z_L = R_L(3, 1) * Xl(i) + R_L(3, 2) * Yl(i) - R_L(3, 3) * f
        X_R = R_R(1, 1) * Xr(i) + R_R(1, 2) * Yr(i) - R_R(1, 3) * f  '右片
        Y_R = R_R(2, 1) * Xr(i) + R_R(2, 2) * Yr(i) - R_R(2, 3) * f
        Z_R = R_R(3, 1) * Xr(i) + R_R(3, 2) * Yr(i) - R_R(3, 3) * f
        
        '求点投影系数
        NL = (Bx * Z_R - Bz * X_R) / (X_L * Z_R - X_R * Z_L)
        NR = (Bx * Z_L - Bz * X_L) / (X_L * Z_R - X_R * Z_L)
        
        '求地面摄影测量坐标
        X(i) = (NL * X_L + NR * X_R + XsL + XsR) / 2
        Y(i) = (NL * Y_L + NR * Y_R + YsL + YsR) / 2
        Z(i) = (NL * Z_L + NR * Z_R + ZsL + ZsR) / 2
        
        '求地面测量坐标：把地面摄影测量坐标换回地面测量坐标——把x,y互换回来
        Dim temp#               '交换辅助变量
        temp = X(i): X(i) = Y(i): Y(i) = temp
        
        '显示计算结果
        txtShow.Text = txtShow.Text & "第" & Str(i) & "个待测点的地面坐标：" & X(i) & " , " & Y(i) & " , " & Z(i) & vbCrLf
    Next i
End Sub

'根据输入的旋转角计算旋转矩阵：
'输入旋转角fai、omg、kap
'              |a1 a2 a3 | |R11 R12 R13|
'输出旋转矩阵R=|b1 b2 b3 |=|R21 R22 R23|
'              |c1 c2 c3 | |R31 R32 R33|
Public Sub subGetR(fai#, omg#, kap#, R)
    Dim cosFai#, sinFai#, cosOmg#, sinOmg#, cosKap#, sinKap#
    cosFai = Cos(fai): sinFai = Sin(fai)
    cosOmg = Cos(omg): sinOmg = Sin(omg)
    cosKap = Cos(kap): sinKap = Sin(kap)
    R(1, 1) = cosFai * cosKap - sinFai * sinOmg * sinKap
    R(1, 2) = -cosFai * sinKap - sinFai * sinOmg * cosKap
    R(1, 3) = -sinFai * cosOmg
    R(2, 1) = cosOmg * sinKap
    R(2, 2) = cosOmg * cosKap
    R(2, 3) = -sinOmg
    R(3, 1) = sinFai * cosKap + cosFai * sinOmg * sinKap
    R(3, 2) = -sinFai * sinKap + cosFai * sinOmg * cosKap
    R(3, 3) = cosFai * cosOmg
End Sub

'空间后方交会的通用过程
Public Sub subSpaceResection(fai#, omg#, kap#, Xs#, Ys#, Zs#, X#(), Y#(), esp#)
    Dim dX#(1 To 6)             '空间后方交会中的未知数向量，对应6个外方位元素的改正数
    Dim A#(1 To 8, 1 To 6)       '误差方程的系数矩阵
    Dim L#(1 To 8)              '误差方程的常数向量
    Dim R#(1 To 3, 1 To 3)      '旋转矩阵
    Dim i%, bLoop As Boolean
    
    
    Do
        subGetR fai, omg, kap, R    '计算旋转矩阵
        '组成误差方程的系数矩阵
        For i = 1 To 4
'            A(2 * i - 1, 1) = -f / H * RU: A(2 * i - 1, 2) = 0: A(2 * i - 1, 3) = -X(i) / H * RU
'            A(2 * i - 1, 4) = -f * (1 + X(i) * X(i) / (f * f)): A(2 * i - 1, 5) = -X(i) * Y(i) / f: A(2 * i - 1, 6) = Y(i)
'            A(2 * i, 1) = 0: A(2 * i, 2) = -f / H * RU: A(2 * i, 3) = -Y(i) / H * RU
'            A(2 * i, 4) = -X(i) * Y(i) / f: A(2 * i, 5) = -f * (1 + Y(i) * Y(i) / (f * f)): A(2 * i, 6) = -X(i)
            A(2 * i - 1, 1) = (R(1, 1) * f + R(1, 3) * X(i)) / H * RU: A(2 * i - 1, 2) = (R(2, 1) * f + R(2, 3) * X(i)) / H * RU: A(2 * i - 1, 3) = (R(3, 1) * f + R(3, 3) * X(i)) / H * RU
            A(2 * i - 1, 4) = Y(i) * Sin(omg) - (X(i) * (X(i) * Cos(kap) - Y(i) * Sin(kap)) / f + f * Cos(kap)) * Cos(omg)
            A(2 * i - 1, 5) = -f * Sin(kap) - X(i) * (X(i) * Sin(kap) + Y(i) * Cos(kap)) / f
            A(2 * i - 1, 6) = Y(i)
            A(2 * i, 1) = (R(1, 2) * f + R(1, 3) * Y(i)) / H * RU: A(2 * i, 2) = (R(2, 2) * f + R(2, 3) * Y(i)) / H * RU: A(2 * i, 3) = (R(3, 2) * f + R(3, 3) * Y(i)) / H * RU
            A(2 * i, 4) = -X(i) * Sin(omg) - (X(i) * (X(i) * Cos(kap) - Y(i) * Sin(kap)) / f - f * Sin(kap)) * Cos(omg)
            A(2 * i, 5) = -f * Cos(kap) - Y(i) * (X(i) * Sin(kap) + Y(i) * Cos(kap)) / f
            A(2 * i, 6) = -X(i)
        Next i
        
        '计算误差方程的常数项
        Dim Tx#, Ty#, Tz#, t#
'        Debug.Print "to show the L matrix:"
        For i = 1 To 4
            Tx = Xtp(i) - Xs: Ty = Ytp(i) - Ys: Tz = Ztp(i) - Zs
            t = R(1, 3) * Tx + R(2, 3) * Ty + R(3, 3) * Tz
            L(2 * i - 1) = X(i) + f * (R(1, 1) * Tx + R(2, 1) * Ty + R(3, 1) * Tz) / t
            L(2 * i) = Y(i) + f * (R(1, 2) * Tx + R(2, 2) * Ty + R(3, 2) * Tz) / t
'            Debug.Print L(2 * i - 1)
'            Debug.Print L(2 * i)
        Next i
        
        '解算误差方程：调用有关的函数
        Dim P#(1 To 8, 1 To 8)
        For i = 1 To 8
            P(i, i) = 1
        Next i
        InAdjust A, P, L, dX    '调用间接平差通用过程解算
        
        '收敛判断：0.0004对应于0.1秒
        bLoop = False
        For i = 1 To 3
            If Abs(dX(i)) > esp Then bLoop = True
        Next i
        For i = 4 To 6
            If Abs(dX(i)) > esp * 1000 Then bLoop = True
        Next i
        
        '计算外方位元素的结果
        fai = fai + dX(1): omg = omg + dX(2): kap = kap + dX(3)
        Xs = Xs + dX(4): Ys = Ys + dX(5): Zs = Zs + dX(6)
    Loop While bLoop
End Sub

'空间前方交会的通用过程
Public Sub subSpaceIntersection()

End Sub

'获取定向参数
Public Sub GetIO(strFileName As String, a0#, a1#, a2#, b0#, b1#, b2#)
    Dim strTemp As String, iTemp1 As Integer, iTemp2 As Integer
    
    Open CDg1.FileName For Input As #1
        While Not EOF(1)
            Line Input #1, strTemp
            If InStr(strTemp, "a0=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                a0 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "a0=" & Str(sa0)
            End If
            If InStr(strTemp, "a1=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                a1 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "a1=" & Str(sa1)
            End If
            If InStr(strTemp, "a2=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                a2 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "a2=" & Str(sa2)
            End If
            If InStr(strTemp, "b0=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                b0 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "b0=" & Str(sb0)
            End If
            If InStr(strTemp, "b1=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                b1 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "b1=" & Str(sb1)
            End If
            If InStr(strTemp, "b2=") > 0 Then
                iTemp1 = InStr(strTemp, "="): iTemp2 = InStr(strTemp, ",")
                b2 = Val(Mid(strTemp, iTemp1 + 1, iTemp2 - iTemp1))
                strTemp = Right(strTemp, Len(strTemp) - iTemp2)
                'txtShow.Text = txtShow.Text & vbCrLf & "b2=" & Str(sb2)
            End If
        Wend
    Close #1
End Sub