Tuesday, April 17, 2012

Photogrammetry Details







One of the aspects not mentioned in the "Turning Points" description is the relative costs of photogrammetric rectification packages. I found them from $1000, most at $5000 and many to $25,000. As stated the public domain packages are time intensive to acquire sufficient expertise to utilize efficiently.


Shown here are the first two(2) figures.


This is Photo 147 cropped showing Highway 278 at Berlin

RESOLVING THE INTERNAL PHOTO RELATIONSHIPS
begin program output:
6
 11.836        76.037        152.4         669834.458    463903.566    256.253
-68.496        65.061        152.4         668211.814    463681.834    250.363
-19.952       -79.102        152.4         669225.487    460785.233    234.913
 68.597       -97.614        152.4         670994.677    460440.17     230.545
-70.431       -4.206         152.4         668181.048    462277.202    231.57
 56.594        4.138         152.4         670742.478    462467.424    243.797

 0             0             0             0             0             5000
*********** end of Input file ************************************************
omega =        0             phi =         0            kappa =        0
XL =           0             YL =          0             ZL =          5000
****************** BEGIN ITERATION ONE ***************************************

Comment: There will be 10 iterations, one pass for each control point--------------------------------------------------
Rotation matrix for Iteration:             1
 1             0             0
 0             1             0
 0             0             1
Rotation Matrix Check
 1             0             0
 0             1             0
 0             0             1
Pt Number      1
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21507.6024100164
vya =         -14827.5620870508
xa1 =          21519.4384100164
ya1 =          14903.5990870508
 11.8359999999993
 76.0370000000003
Pt Number      2
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21509.1830364198
vya =         -14812.9415719018
xa1 =          21440.6870364198
ya1 =          14878.0025719018
-68.4959999999992
 65.0609999999979
Pt Number      3
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21423.5411094538
vya =         -14816.2242202659
xa1 =          21403.5891094538
ya1 =          14737.1222202659
-19.9520000000011
-79.101999999999
Pt Number      4
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21371.9216283967
vya =         -14810.2136383235
xa1 =          21440.5186283967
ya1 =          14712.5996383235
 68.5969999999979
-97.6139999999996
Pt Number      5
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21425.6342252125
vya =         -14778.6801109338
xa1 =          21355.2032252125
ya1 =          14774.4741109338
-70.4310000000005
-4.20600000000013
Pt Number      6
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -21435.5823531119
vya =         -14814.4127257785
xa1 =          21492.1763531119
ya1 =          14818.5507257785
 56.594000000001
 4.13800000000083
---------------------
array VectorOut, sub Solve()
 82152821.8583712
-119136061.605371
-426201.315633722
-6966.59328362713
-2845.77156829018
 208.883605432399
 0
 0
 0
 0
---------------------
Sub gaussj
residuals for iteration      1
-0.01013194865169 -1.30050869153706E-02 -8.68642904038754E-03 -673688.904051919  217178.576640219 -2329.78383791447 
omega =        0.01013194865169            phi =         1.30050869153706E-02       kappa =        8.68642904038754E-03
XL =           673688.904051919            YL =         -217178.576640219            ZL =          7329.78383791447

Comment: This is the end of iteration 1,  shown are the residuals----------------------------------
Comment: There are 8 iterations similar to the above.
Comment: Begin the last iteration:
---------------------
Rotation matrix for Iteration:             10
 0.999878738904567           8.73254449169687E-03       -1.28938029034804E-02
-8.60001326212023E-03        0.999909985332769           1.02985923170277E-02
 0.012982575187713          -1.01864566224728E-02        0.999863835151054
Rotation Matrix Check
 0.999878738904567           8.73254449169687E-03       -1.28938029034804E-02
-8.60001326212023E-03        0.999909985332769           1.02985923170277E-02
 0.012982575187713          -1.01864566224728E-02        0.999863835151054
Pt Number      1
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =          0.043210130460256
vya =          6.03101648868778E-02
xa1 =          11.7927898707966
ya1 =          75.9766895748996
 11.8360000012569
 76.0369997397865
Pt Number      2
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -3.40662949737127E-03
vya =         -0.048007763032459
xa1 =         -68.4925933702779
ya1 =          65.1090075062389
-68.4959999997752
 65.0609997432065
Pt Number      3
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -8.35992118613723E-02
vya =         -8.55772187611816E-02
xa1 =         -19.8684007895645
ya1 =         -79.0164230437273
-19.9520000014258
-79.1020002624885
Pt Number      4
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =          6.81194962058081E-03
vya =          5.40036079379968E-02
xa1 =          68.5901880507898
ya1 =         -97.6680038737459
 68.5970000004104
-97.6140002658079
Pt Number      5
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =          0.062432048006606
vya =          1.43609441158232E-02
xa1 =         -70.4934320492948
ya1 =         -4.22036120131746
-70.4310000012882
-4.20600025720164
Pt Number      6
formulating b matrix
linerarized Collinearity condition for current parameters
vxa =         -2.52434608293936E-02
vya =          4.2791432117453E-03
xa1 =          56.6192434617323
ya1 =          4.13372059434486
 56.5940000009029
 4.1379997375566
---------------------
array VectorOut, sub Solve()
-3.9193457756315E-06
-2.27558512645529E-07
 6.38136271862301E-08
-1.07887282052499E-09
-1.07936164293905E-09
 3.75699982703177E-11
 0
 0
 0
 0
---------------------
Sub gaussj
residuals for iteration      10

Comment: The following are the residuals for -10-.... notice values are small.
-6.58650037074729E-14 -1.05427665029161E-13 -1.77141193660531E-13  2.74643877656817E-10 -7.27011384398665E-08 -3.99021923776564E-11

Comment: Then, finally the 3 rotations and the coordinates of the Observation Point. Also in the actual output file, an ASCII file, these values will be on one line.
  Final Values: omega   phi  kappa XL  YL  ZL -------------------------
 1.01874913978596E-02        1.29829399116193E-02        8.60084414823962E-03        669650.413992352            462343.808977171            3305.40699175135
degrees
 0.583700260923186           0.743867662607736           0.492792069943922

Comment: End of Output File.
After the internal characteristics are determined for both photos, omega phi kappa XL YL ZL, we can then utilized those photo relationships to establish an unknown ground point.

At this milestone of the process, the redundant input can "shroud" a poor control point (or the photo coordinates), just as a single measurement can be hidden for a survey that has 100 observations. Just one can have high residuals and since the remainder have small residuals they "overwhelm" to poor point.
Shown is the stereo pair, the right photo has all rotation angles = 0. Notice the difference in ground coverage. Of course all parameters are exaggerated for clarity of the viewer.
Figure showing an unknown ground point and the ray intersection with both photo positives.


These are the OP values w/o adjustment for earth curvature: 

Final Values: omega   phi  kappa XL  YL  ZL -------------------------
 1.01874913978596E-02        1.29829399116193E-02        8.60084414823962E-03       

      XL 669650.41                  YL 462343.81           Z 3305.41
degrees
 0.583700260923186           0.743867662607736           0.492792069943922


These values are final OP with curvature adjustments:

 Final Values: omega   phi  kappa XL  YL  ZL -------------------------
 8.87707101416532E-03        1.35039640397676E-02        9.21384170508731E-03       

     XL  669652.98                  YL  462349.56            Z  3306.65
degrees
 0.50861870354959            0.773720146175117           0.527914242803126


The phi & kappa rotations show some variations.

data concerning atmospheric refraction:
See Analytic Aerotriangulation, 1962 Coast And Geodetic Survey, Technical Bulletin No 21, Dated July 1963
He states an average terrain height can assumed over the area involved.  The K coefficient is dependent upon the barometric pressure and the wavelength of light and other parameters. These can be simplified for the project area; thus for our applications we can normally assume both flight paths are completed in the same day. Paul Wolf also has a slightly different approach.
Harris, 1962.
                                                                              
Also note here that Harris used f = c, so the focal length "f" is the parameter "c" in the equation.


So here we see the external conditions being evaluated using the collinearity equations for two lines.













Wolf states:
"These are the two equations for a straight line in an image, corresponding to the two collinearity equations for a point in an image. The photogrammetric condition equations using straight lines will fail when the straight line in object space is in the same plane as the two (or more) image perspective centers (the epipolar condition). When this occurs, the projection planes for the corresponding straight lines from the images are coplanar. Since the planes do not intersect in a single line, the object line cannot be  determined, Consequently, in practice, straight lines that are parallel or are close to parallel to the epipolar plane should be avoided. "

Sample input for the solution of the unknown value. This routine is yet to be finished. The image adjustments are made from a spreadsheet:

9652.79465192302            2349.70930221402            3306.963    8.8332392387544E-03         1.34473361198426E-02        9.21546910637473E-03  152.4
9673.07190937109            4018.56307591905            3307.66     8.04394743779159E-03        8.42477365500256E-03        0.013275825020702     152.4
3
 56.6068          4.017     53.709       -79.604     10742.478     2467.424      0
  64.3287      50.4046    61.6548      -33.0248   10838           4634         235
  64.1585      49.7376    61.4927      -33.6192   10838           4634         235


For the above: You must have  an estimated value for the unknown point so that GAUSS will know where to start, i.e. 10838  4634  235.
Some Output:

Estimated point coordinates and iteration deltas, Iteration 1
XL :  10742.575       0.097 YL :   2467.669       0.245 ZL :    261.430     261.430
Image Errors Left   Photo   4.085     0.405          Right      3.801    -5.763

...........................................
Estimated Point coordinates and iteration deltas, Iteration 10
XL :  10742.574       0.000 YL :   2467.652      -0.000 ZL :    242.167      -0.000
Image Errors Left  Photo   0.074    -0.001         Right    -0.074     0.002


How would you check yourself?
Include a point, say #1 of the 3 above that has a known elevation and position. Or say one that has a known elevation. Very prudent issue.
Add to bibliography shown in "Turning Points":

9)Sprott, Julien, Numerical Recipes Routines and Examples in BASIC, Cambridge University Press, 1998

10) Press, William Et Al., Numerical Recipes in Fortran 77, Cambridge University Press, 2003

All contents Copyright 2012, J M Hillman



No comments:

Post a Comment