Python programming
GGE 4313 – AIRBORNE MAPPING SYSTEMS LAB #2 – Coding Collinearity Equations, Due Date: Feb 14th at 11:59 PM Lab Outline During this project, students will get familiarized with the fundamentals of Python programming required in mathematical models in photogrammetry. The students also learn how to code collinearity equations to estimate ground and image coordinates and determine Exterior Orientation (EO) parameters. Learning Objectives • Learning Python programming fundamentals including defining new variables, writing expressions, using built-in functions, creating and using lists, and Python control instructions • Learning some more advanced programming concepts namely symbolic mathematics in python and related functions • Understanding and coding non-linear equations in photogrammetry • Underdoing and codding collinearity equations in photogrammetry Data The data consists of 11 GCPs that are visible in 2 stereoscopic images. You can download the table from D2L under Lab#2 materials (an excel file named “Data”). Deliverables 3 Jupyter notebook files including your codes for the algorithms described in part1, part2, and part3. (You need to have a file for each part separately) A pdf file including screenshots of your results for each part (screenshots of the results that you get by running your code. You need to have ONLY ONE pdf file for all three parts) Part 1 (20%): • In part1 we assume that we know all required parameters and we simply want to use collinearity equations for obtaining image coordinates (image 1) from object coordinates (use equation 1 in part 2). • Considering the following IO and EO parameters and collinearity equations write a code that calculates the image coordinates (x, y) of GCPs 7, 8 and 9 by using their ground coordinates (X, Y, Z). (angels are in degrees) • f= 0.008604233 (m), x0= 0.006416659 (m), y0= -0.00427777248(m) • ?? = 2548878.98, ??= 7464989.11 , ?? = 61.8, omega = 0.09, phi = 0.04, kappa = 21.76 Coding Collinearity Equations Winter 2021 Part 2 - Single Photo Resection (50%): In part2 we assume that we only have IO parameters and we also have the image and ground coordinates of some control points (use GCPs 1 to 6). Now we need to estimate EO parameters by solving a system of non-linear equations using least square optimization. Since collinearity equations are non-linear, we use Newton-Raphson Method to make them linear. To do that we need to follow these steps: Step 1: For all GCPs (we have n GCPs) estimate the approximate image coordinates using the following initial values for EOPs which is called computed image coordinates (we already know the true values of the computed image coordinates but this estimation is part of the parameter optimization process). You use the following collinearity equation for this estimation. The computed image coordinates are denoted by ?1̅̅ ̅ , ?1̅̅ ̅, ?2̅̅ ̅, ?2̅̅ ̅, etc. Equation (1) Hint : Use the code in Appendix for simplicity. Step 2: Compute the partial derivatives of the collinearity equation (Equation 1) with respect to the EOPs. By repeating the computation for all data points, we form matrix A. Example: From step1 we have >>> x = x0 + (-f) * (M[0,0]*(X-X0) + M[0,1]*(Y-Y0) + M[0,2]*(Z-Z0)) / (M[2,0]*(X- X0) + M[2,1]*(Y-Y0) + M[2,2]*(Z-Z0)) we calculate the first element of A ( ??1 ??0 ) by first calculating ?? ??0 by using sym.diff(): >>> diff_x_X0= sym.diff(x,X0) You Should expect a result with the following structure (values might be different): (sin(kappa)*sin(omega) - sin(phi)*cos(kappa)*cos(omega))*(-0.008604233*(X - Xl)*cos(kappa)*cos(phi) + 0.008604233*(Y - Yl)*sin(kappa)*cos(phi) - 0.008604233*(Z - Coding Collinearity Equations Winter 2021 Zl)*sin(phi))/((X - Xl)*(sin(kappa)*sin(omega) - sin(phi)*cos(kappa)*cos(omega)) + (Y - Yl)*(sin(kappa)*sin(phi)*cos(omega) + sin(omega)*cos(kappa)) + (Z - Zl)*cos(omega)*cos(phi))**2 + 0.008604233*cos(kappa)*cos(phi)/((X - Xl)*(sin(kappa)*sin(omega) - sin(phi)*cos(kappa)*cos(omega)) + (Y - Yl)*(sin(kappa)*sin(phi)*cos(omega) + sin(omega)*cos(kappa)) + (Z - Zl)*cos(omega)*cos(phi)) Then we calculate ??1 ??0 by substituding the EO variables with initial values: Val_x_X0= diff_x.subs ([(omega,omega_init), (phi,phi_init), (kappa, kappa_init), (Xl,Xl_init),(Yl, Xl_init),(Zl, Xl_init)]) The followings are the initial values for the EO parameters (?????). Image 1: ?? = 2548880, ??= 7464992, ?? = 61, omega = 0, phi = 0, kappa = 21.0 Image 2: ?? = 2548870, ??= 7465010, ?? = 62, omega = 0, phi = 0, kappa = 21.0 Step 3: Compute corrections values of the EO parameters ∆? by using the following equation: where ∆Y is the difference between computed and true image coordinates computed as: ∆? = [ ?1 − ?1̅̅ ̅ ?1 − ?1̅̅ ̅ …… ?? − ??̅̅ ̅ ?? − ??̅̅ ̅] Step 4: Update EO parameters ξ: ???? = ????? + ∆? Step 5: Repeat steps 2 to 4 until convergence (∆? should be less than 0.001) Part 3 – Space Intersection By Collinearity (30%): In space intersection, two photos (in a stereo pair) with known exterior orientation parameters are used to determine the spatial position of an object point in the overlap area since the two rays to the same object point from the two photos must intersect at the point. Using the code you wrote in part2, calculate EO parameters of both images. Then write a code that estimates the ground coordinates of the GCPs 7 to 11 using their image coordinates. You need to use the collinearity equations for obtaining object coordinates : Coding Collinearity Equations Winter 2021 You need to use the initial ground coordinates that is provided in the data table. Appendix Collinearity equations: x = x0 + (-f) * (M[0,0]*(X-X0) + M[0,1]*(Y-Y0) + M[0,2]*(Z-Z0)) / (M[2,0]*(X-X0) + M[2,1]*(Y-Y0) + M[2,2]*(Z-Z0)) = y0 + (-f) * (M[1,0]*(X-X0) + M[1,1]*(Y-Y0) + M[1,2]*(Z-Z0)) / (M[2,0]*(X-X0) + M[2,1]*(Y-Y0) + M[2,2]*(Z-Z0)) Rotation Matrix : M = Array([[sym.cos(phi)*sym.cos(kappa), -sym.cos(phi)*sym.sin(kappa), sym.sin(phi)], [sym.cos(omega)*sym.sin(kappa)+sym.sin(omega)*sym.sin(phi)*sym.cos(kappa),sym.cos(omeg a)*sym.cos(kappa)-sym.sin(omega)*sym.sin(phi)*sym.sin(kappa), - sym.sin(omega)*sym.cos(phi)]]) [sym.sin(omega)*sym.sin(kappa)-sym.cos(omega)*sym.sin(phi)*sym.cos(kappa), sym.sin(omega)*sym.cos(kappa)+sym.cos(omega)*sym.sin(phi)*sym.sin(kappa), sym.cos(omega)*sym.cos(phi)])