WHHIWIWINWIWWWWHWIHHNHWHI 109 058 __THS ”5'3 F, manna? Michigan 5mm University This is to certify that the dissertation entitled A MODEL FOR THE CORRECTION OF THE GEOMETRIC DISTORTION OF MULTISPECTRAL SCANNER DATA presented by Henry Richard Gee has been accepted towards fulfillment of the requirements for PhD degree in Mathematics . / ~ ~ \ ' .. 1' Major professor l 2 l Tien Yien Li Date Ju y 6, 985 MS U i: an Affirmative Action/Equal Opportunity Institution 0-12771 MSU LlBRARlES RETURNING MATERIAL§: Place in book drop to remove this checkout from your record. filfl§§_wili be charged if book is returned after the date stamped below. A I‘lODEL FOR THE CORRECTION OF THE GEOMETRIC DISTORTION OF NULTISPECTRAL SCANNER DATA 39 Henry Richard Gee A DISSERTATION Submitted to flichigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics I 985 ABSTRACT A MODEL FOR THE CORRECTION OF THE GEOMETRIC DISTORTION OF HULTISPECTRAL SCANNER DATA 39 Henry Richard Bee The purpose of this work is to solve a problem suggested by Dr. David Coupland, in a talk given at Michigan State University. A multispectral scanner, mounted in an aircraft, collects data from a certain geographical location, and the objective is to print a geometrically accurate picture of the region using the data. The solution presented herein, involves modelling the location and orientation of the aircraft using cubic spline functions. The coefficients of these spline functions are determined by obtaining a set of “match points“, plugging these into the modelling equations, and solving the resulting nonlinear system by a generalized Newton's method, to obtain a "best fit" for the data. An efficient algorithm for solving this problem is discussed in this work, and a full FORTRAN program listing is given in the appendix. ACKNOWLEDGEMENTS I want to express my gratitude to all of the members of my guidance committee for their help during my years at Michigan State University. and for all of the interest they have shown. I especially want to thank Professors Tien Yien Li and Richard Hill for their guidance and patience in the preparation of this dissertation. TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES Chapter I: INTRODUCTION l.l= Multispectral Scanners 1.2: Description of the Geometric Problem Chapter 2: CONSTRUCTION OF A MODEL Chapter 3: THE NUMERICAL MODEL 3.1: The Use of Cubic Splines 3.2: Iterative Solutions Chapter 4: THE ALGORITHM 4.l: Input and Initialization 4.2: Computation of F and J 4.3: Solving the Linear Least Squares Problem 4.4: Choosing Step Size and Updating 4.5: Variations Chapter 5: NUMERICAL RESULTS APPENDIX LIST OF REFERENCES TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES Chapter 1: INTRODUCTION 1.1: Multispectral Scanners 1.2: Description of the Geometric Problem Chapter 2: CONSTRUCTION OF A MODEL Chapter 3: THE NUMERICAL MODEL 3.1: The Use of Cubic Splines 3.2: Iterative Solutions Chapter 4: THE ALGORITHM 4.1: Input and Initialization 4.2: Computation of F and J 4.3: Solving the Linear Least Squares Problem 4.4: Choosing Step Size and Updating 4.5: Variations Chapter 5: NUMERICAL RESULTS APPENDIX LIST OF REFERENCES iii iv iv \OQU‘IN 13 17 I7 19 21 24 27 29 35 43 LIST OF TABLES TABLE 1: The Use of Weights TABLE 2: The Use of Weights (continued) TABLE 3: Variations in the Number of Intervals LIST OF FIGURES FIGURE 1: Variations in the Viewing Area FIGURE 2: The Orientation of the Aircraft FIGURE 3: The Plane of the Scanner iv 31 32 33 Chapter 1: INTRODUCTION 1.1: Multispectral Scanners In this section multispectral scanners (M88) and their operation will be described. Although much of this section will apply to MSS's in general, specific information will be given about the AAD31268 Digital 1153 System produced by Daedalus Enterprises, Inc, since this is the model that collected the data that is used in the section on numerical results. A multispectral scanner is a data—collecting device which may be mounted in an aircraft. As the aircraft flies over a region, the 1183, by means of a rotating mirror, scans the ground below from right to left, and measures the intensities of several bands of radiation, generally in the range of .3 to 14.0 pm. The mirror is mounted on a rotating, horizontal shaft at an angle of 45° to the shaft. Light or infrared radiation emanating from the ground strikes the mirror and is reflected horizontally into some focusing optics. The infrared bands are split off using dichroic filters, and the visible light is separated by a prism. The intensity of several bands is measured and stored on magnetic tape. Later the tape can be read, and an image of the region can be produced. The mirror rotates at a constant angular speed and the sampling is done at a constant rate. The samples are recorded according to their scan line and pixel number. In the case of the Daedalus scanner each line consists of 716 pixels, and the mirror rotates through an angle of 1.5 radians in going from pixel 0 through pixel 715 in each row. One of the advantages of this type of imaging, over ordinary 1 a .i! r- 0 I d s L I ,,i i. a \ ‘ ‘ t it a " in? .. . I“ ‘ V ~ 4 L a 3 ’ . . . l 1 v a ‘ O I . Hit: , :u_ -. s ‘- »r . 3 ~.: :..... , ii ’iiIV ': unit: '3‘, .1 ri‘ iLI. Hit. L: c~ inlnh 3r. 1' Pin . "5 I I' ' :iD a . .f . .‘ _ e u ,“ l. i i; ‘ii it , a, it libti'lt. tit. m s e : v. . ' . . i t ~‘iir )« ' '='i' "ii"' " ' ' I.. .. 4 ' ‘ E. y . .~." . g.‘ .1 t 1“ .‘j i it it): - '- '- ..n} 5' ' - ‘.; , o ‘; _1i':§ii. ot'? . . i 1.‘ I " | 5‘ ' 'L': if. ‘ "I i ?r ”all: 11’» .1 ‘u‘ :iiz’e’ul :. :2' ‘ Mi: ‘ is air-i: an? ”'31”- n pm i. ..., ... '2'er . ii‘i‘il’. 'iim . nun : m ... y D . ‘ I “' , ' ‘.. I § , i. n . w... mi‘i ,in up: I ".ti‘ [1 312:..- l '3! ”x- it}: . it” v ‘ ‘:({ .g-‘ '1. ?| ,"' a b f'l‘ g: :“ g: |p '1! . . I. b v v I ' .‘ 3 g T ‘ ' .' z ' I " i’ L' . , “ 1' ‘ ' ' ' ‘i w "t l' 'vL .‘-:.L- ‘.‘1 _. ' ‘ .J ”I?“ I‘. 3‘2... " ' , .. ' " I. ' I ' “‘ ‘ ’7 ’l ‘_| "f ' ‘ a ‘ I . V. ‘ a l J‘ I K 1. i? .3 "3i ‘ " “I ‘ i ‘. s-I'; p' ‘9 t ' :3: ‘ u i t 3 . fl) 1 1‘ l , ‘ O 0 I t l I. ‘. . _ i . . .Ti‘ « photography is that the data, being in digital form, can be processed in a variety of ways, so that one might, for example, ’ highlight certain tunes of vegetation, or certain geological formations. A second advantage is that the geometry of the image may be easily altered, by changing the function used to plot the pixels, as compared to a photograph which, once taken, is not readily corrected. This second point is the subject of this work. 1.2: The Geometric Problem This paper deals with the geometric accuracy of the image that can be produced Distortion of this image arises in a variety of ways. First, the size and shape of the area in view at any instant is not constant, but rather it is a function of the direction and distance from the scanner to the object being viewed. The field of view for the Daedalus model is within a cone of angle 2.5 mrad. If the scanner is directed vertically downward then the field of view is nearly circular, but as the scan beam moves away from vertical, the field becomes elliptical. For example (See FIGURE 1), since pixel number 0 in any row is of an area at an angle of .75 radians to the right of vertical FIGURE 1: Variations in the Viewing Area w ‘ .“2:L- iii "vii"- :'3=Jt.".L‘U“"t’3 ' " .i‘w in wimmvi , .« .3 . _ . r 1": I i;. ‘- .. _ . ‘ " . ‘ s , ‘ . " i l. “ f." l I: .. ‘ 1-‘ 311”” 5;.) ’ ~i. 1".- "l L' I! I ‘J t. IIL‘ . . a ‘ . n. ‘ , .3 c .u_ H. ‘ ‘o .‘ o '\ . g ,p ‘ t_. “ I , I . . crawl u a ‘J t“. Pi , ’. "I I. . ft'blhi O , ‘ I.. V. a § .‘ A . ', v‘ ‘ f: * ' 1 ' I. . .’ iI‘ ._1. ‘ {‘- ',:: I ,‘ :1_ i’ .I III ‘:~ T]; Lditiyftt‘“. t , '. | . . a :a I . . 1‘. ; Lt? 02“, a ‘ It ‘1 . 1")..." :HI . 3‘ I'le‘ I I.) f"( ,glft“ t "3“ " . : v {at A 1 . ' Iiii" '- l ‘i 4“? "I " .J "L ‘1. ’ “’ . ‘-..‘o‘ , r I ' i I ‘ ‘. . .. : a: _ , ~ . - - w » . (I ”' l‘r.2.} :, “ I I . 1 ‘ . a T' 'i : a iv. “'1. i. , O . " - ' ' ‘ v i '.‘ ‘ . ~ .. . it ~ ‘. .i 1‘: - . wet. I ‘ . . : .- -.ii.i~-- 11: it . ’ l 1‘\T l.‘ , ~ A v| - a I 4 V t. ‘I r I b w .. i ‘ . _ . . . F . 1‘ . ‘ a \ . ’ ' .‘ ‘:"- . 'o . :— I‘. '3 'a‘ N.’. ‘ ~ .1 'i V i: ‘l . .4 a. i .‘. . 9 ‘ . ‘ ‘ D ' . u , ‘ ‘ . V! t. 5. ‘ t . v . i 1‘: ' - . 'a n i“._' | ‘ i'. i t . . i a f . ‘ l ‘1 ' : ' ‘ . ' ‘l ' I’ -..:'m_:: 2-1: ’ , ,~ ~. « .‘ , I S I ‘ .‘,’, l , 'l i 1 ‘ .tg‘ . a . . from the scanner, the view is of an ellipse having a focal width given by 2h sin.00125 I cos.75 feet, whereas pixel number 358 is a view of an area directly beneath the aircraft, so that it represents a circle with diameter 2h sin .00125 feet, where h is the altitude in feet of the scanner. Second, the forward motion of the aircraft perpendicular to the direction of scan, the yaw of the aircraft, and the constant scan rate cause the scan lines to be somewhat S-shaped. The forward motion and the yaw cause a shear and the constant scan rate causes the pixels to be spread out near the edges and be compressed near the center of each line. Finally turbulence in the atmosphere may cause some erratic motion of the aircraft union will distort the image. Most liSS‘s have a roll compensation device to partially correct for this. The Daedalus accomplishes this by means of a gyroscope, which keeps track of the vertical. If the aircraft has rolled to the right for example, then the scanner will recognize this and start the sampling for that row sooner and stop sooner correspondingly. The methods currently in use for making the necessary correction to this distortion usually involve assumptions such as flat earth, and level, straight flight path. One views the set of pixels as a rectangular array, with 716 columns and an arbitrary number of rows. A set of “match points“ can be obtained: through visual inspection of an uncorrected image, one may be able to associate certain pixels to certain points on the map. It is then possible to derive a bilinear transformation, which will map a rectangular array of data onto an 113V”) IIII'Vvi 1131-51 13 t)!‘:‘::1 "1- 111?") qr.- "_ -: ‘-.."3:'- «‘1 "1413;, _ .-.§‘ ”ii-'51 115 If:- wuv i- 1).." ".tslim'.” fun: :ri~.i"'l"\. ,l-'~‘.l= Liz-1.. \ ;' ‘."u.':i- " Ln! iiirw SUIT} s ....=ir.~.:-iqsi :« mm o. JIE—‘a‘ifl'tlb em nimmsu pl‘w 1:..- can: 3111 to last .11 about“: an! at ii 3131114 .1331 65:100. me 11'. wail-imp I;- "twine-3:: slit of 'Ir-IUHDTI'SQ'VSQ 116131”; 3:11 to nonlnm bwewwl am ,bnoosc' 9161 ”6313 instance suit has Jimmie .mi is won sdi need to notiaswib brie nouom b1awinl sdl pagoda-8 teaiwsmoe ad 01 estiil nose am sauce 01 eisxiq sill «Anus.- siai i‘lea 111818110.) anti one incur. a seven way 1:11 doe-'9 in wins) 3111 wean baeasaqmoo mi pm. uf“LJIitJ 3111 teen iuo L-os'iqe an ‘iilll notiom 3136119 smog scues yam eisuiqamms an! m SitiSlUiI‘lUi ulicm? lIo-i L- Welt {3811 icon 3011111? +1111 I'Ti_-I:..'IIJ iliw 1mm 1115mm aim to aulebaefl 911T an" 11‘.’ !.J'i'l‘il.\’- [1116111511 01 sown) nosisensqnm.- 3.11 to 'I'Jt'.“l cit-33.1 113111.- .EiiJONU’w e lo ens-«m 9:! anti eastzilqniuvie ‘1113 113111 ,siqmcsa ml iityi'i au.‘ 2:1 1.1.: 1- i can “mum's 3:11 11 its-NIH. tannin; urn 15111 '11." prismatic 3111 Melt. tins anti ’i‘alflglrilv'l H w isnne 1r ”IIJIIIDilL-ulltilii') "tattoo..- out: tutu ._.1 11011311211023 mot-cameo hit unriem ~..-i can an pz'nu'nuu «.‘UI'TII'L‘TI ‘1111 halo ‘111163 1311 EC- 'iJii. *i.i.i:3'l‘::=.'r: triuvul Lil'iun- l-"t" 'le- ~-:"‘1.' . , _1v't’3€ r. ('1. Elam} 1:) 1-27: fivl .‘waw sill! ailbu 1111.1.1 Hiring :1 awn m 'istiuum gnarl-”Ii to 1'6 gm; r-imulo - il l :J N ‘ub‘l h. 1- : l‘U TM '1 ‘ . i; It:- I‘lL' "J'SIIILIII IN -:1 i 1'. ‘.--I '.' :" :. ".1. '1': '1" - “3H,"; "31.. ' :1 '2' " i m law: “1511:: statues-1:1 [-1 521123 .1: yh-‘n 'riu spen- 051 la 1 “our tfz‘PlItilt-i E- 'Jwiiil’: n :--_'. an; H.311? :- 1! Jim: -" ‘:- cum» "uh-2‘ .I:. ”1:1. LL'JI} h. L'i. “ ‘i 1|. ll‘."'l. ' ,.-, . '. “51'. i , (7“. .. . ‘ .. I ._ H! I arbitrary quadrilateral, that will best fit this set of data. Commonly, the entire scan is broken into several segments, and this solution is applied to each segment. One then attempts to properly match up the corners of the various quadrilaterals. In the next chapter a theoretical model for the correction of the distortion is constructed, and in Chapter 3, a numerical model is discussed in which we also use a set of match points and try to obtain a best fit. Chapter 4 describes an algorithm for the implementation of this model, and the last chapter gives the results of some numerical experiments using this algorithm. ,pmummoil 5150 :21 13c «25:11 1:: teed list-x 15:11 IE‘i'flBli‘tbsKLttn 1.. mil die a: uZL-lIUIOII 5:111 We ,cfelamuu lewsvae 91m mil-tr. s: 1167*. sums .3111 '31" cu itriinm ylieqt-g q u? ._-'..ns.i;e was? 'ilU insult-3.: ‘1'163 ‘11 I": ‘ Lin sinisteiu vueup cum‘inv ".1111 tr- clan-:7: am to {IOTITPSYTUL' mil ‘iol lamlii i' ZITITS‘USIII a vii-gosh 1w: mil .1! I‘lfljt'I' l. nuns-.11: u .5 zuiqed‘i :ll 1-1: ,beimi'nial . 4 ~3lllIdlII‘ 1116111005 on late ainroq doibm 1c 19.; c ecu 0815 u slots at: L'*i:.‘del<°.lll in notieinsms'qnu 9111 "IL"? mnirtogie n: asrii‘i‘asht .«'-ieii) .’ 1.sd a Minimum exam: 1-: dime-3': slit csvu. :'.:'qr‘fl‘:~ 1-:2’ ' ‘ imc zebziri 2.1.1.1 iiifliiiuL-ffi ,1“) l." insanisqr; Chapter 2: Construction of a Model Data is recorded on tape according to scan line number, which is denoted by u, and pixel number, denoted by v. The objective is to produce a model that given a pair (u,v) will determine the map coordinates of the point that was scanned, which is herein denoted by" (x,y). To be accurate this model must also take into account the topography of the region, so we assume that the function h(x,y), which gives the elevation of the surface being scanned at the point (x,y), is known. Then the problem is to find functions f, and f2 such that x = fl(u,v,h(x,y)) and y = f2(u,v,h(x,y)). (2-1) Consider a fixed (x,y,z)-coordinate system, with the x,y-plane being horizontal at sea level. For definiteness, take the x-axis to be in the intended direction of flight, and the z-axis to be directed vertically upward. Let x(t), V(t), 2(t) give the position of the scanner relative to this (x,y,z)-coordinate system at time 1. Also let p and a be functions of twhich give the orientation of the scanner at time t as follows (See FIGURE 2). In general an aircraft is subject to 3 rotations; roll, pitch and yaw. Due to the roll compensation discussed above, ignore the first ‘2’ T2 \ ‘P FIGURE 2: The Orientation of the Aircraft .4. vi a a. 6 of these and use 9 and e to model the other two, so that 9 will represent a rotation of the x,2-plane and a will represent a rotation of the x,y-plane. We can transform the standard coordinate vectors 1,3,! in the x,y,z-coordinate system into a set of coordinates 121,3. relative to the scanner by cose -sine cosp 0 sins» cosecosp -sine cosesinp ”13.1311th sine case 0 0 l 0 =1 ml sinecosp cosa sine sing O 0 l -sin90cosp -sinp O c039 Therefore, at any time t the scanner is pointing in a direction that lies in the plane determined by I and t Let a be the angle between the direction of scan and the vector 4. Then a: is in the interval [’%,%1 (for the Daedalus this is [-.75,.75]), and given pixel number v we have eta-oi... + v Act, mere so: = 20:.“ ,p is the constant scan angle between two successive pixels (so: = Lill-"6 =.00209497Zl for the Daedalus). THUS, the direction OT 3080 13 -cos« case sing + sintx sine V: -cos -—‘l [9 c><= c;l Lu rc: <3 c> <31 'o'ob'o —. <3 <3, c: i-l columns 10 {—1 8in ‘2d 80d 8" 8 a" I Biz ° 55 F6 -! 0 + H a. z: 0 + n a. .35 ‘6 '6 " 0 0 0 +fl2 To + 3' 2 :yzg : for i=l,2,...N-l "4' "24 "34 . nnrl4 N—i, v'c> c><31 — ...+ "V"1 n" , ai+l .1 lg? c: c: c>1 [cu i+2 +Hn ON 013) (310 and ' ‘ (Bis 81+] 3 8Hz 3 "' aN-l 3 8no an: 8N2 3N3 )T - k3 Claim: 1 1'11 32d? ((J-Il3-j3w3 ml. 0 1 -2jd (6j-3)d2 o o 1 -3d 0 o o o 3mg]; True by inspection for jsl. Assume it is true for j=k. F —kd kzdz ((k-1)3-k3)d3 1 -d d2 -d3 Then "1mg 0 1 -2kd (alt-3M2 o 1 —2d 3d? 0 o 1 -3a 0 o 1 -3a 0 o o o o o o o 1 (-d-kd) (d2+2kd2+k2d2) (-d3-3kd3-3k2d3) ._. o 1 (-2d-2kd) (3d2+6kd2) o o 1 -3d {2 o o o -(k+l)d (k+l)2d2 (k3—(k+1)3)d3 —. . o 1 -2(k+l)d (6(k+l)-3)d 2 o 1 -3d 0 o o 0 Therefore, using the matrix l1 as shown above we can express the aiJ-‘s in terms of a”, 323»"'»°N-l 3 and 5N0» am, 3N2: 8N3, which leaves N+3 unknown a's. The bij's are related in the same wag as are the ciJ's, dij's, and en's, leaving us 5N+l5 unknowns. The accuracy of approximation by cubic splines has been investigated by many researchers. De Boor [5] p68, and Hall [6] give the following result: 'a K." I 3”. E,- 1 3'.‘Ij.g ‘l_-' .1‘, ‘ i. ‘1 ‘ l D on. 4 .— o. 9 I . I 11" 11 ,. .7- ‘6‘ .. 17’ 1 '1" ‘_ 1". :1" 1 w! : Slt'!‘.- 1.91. "all it 1.‘ o . l 12 m If tcc4ia,til, and ii is a partition of [a,b] with mesh(lT) a a, then there is a cubic spline function scczla,bl with knots in n, such that llf-sll s c WWII 6“- The norm here is the uniform norm on Cla,b]. Hall gives the value of c as 5’384° Suppose we are given uo, v0, and ho. Let 00,90) be obtained from (2.4) using exact functions X ,Y, Z, 9, e, and let (xo*,go*) be obtained from (2.4) using cubic spline approximations X”, Y”, 2", 9", a” with the accuracy provided in the theorem. Then lxo-xo'l = ltxtto)+(ho-ztto))vx(to))- (x*(to)+(ho-z*(toiivx’kto»I s [IX-it'll + Rho-200))“x(to)-vx*(to))l + ll(Z-2")Vx"ll We assume that Vx‘ is bounded, mich is equivalent to the vector V in FIGURE 3 being bounded away from the horizontal, so we have llx-x‘ll + lltz-z*)vx*|| = otd“). The term [(ho-ZU0))(Vx(to)-Vx*(to))l is also 0(d4). To see this note that vx = -cos tx cos a sin 9 «t sin or sin e -cos a cosy so that the assumption that Vx is bounded implies that H9" < 11/2 and ”at" = «m < 11/2, This implies that both avx ._. -cosoz cose + sintx sinp sine 3” cosot c0329 l3 and avx = sintx cose + cosot sinp sine 3° cosot costp are bounded, and therefore llvx-vxfll = 0(d4). Thus we have ‘xo-xo'l = EKG“) and similarly for lyo-yo"l. The term I(ho-Z(t0))(Vx(to)-Vx*(to))l will most likely be the dominant one since Z(to)-ho is the altitude of the aircraft above the level ho. 3.2: Iterative Solutions In this section we discuss two well known techniques for determining the coefficients of the cubic splines. Suppose that we are given a set of match points: (x yy J)<->(u J,vJ-), jsi,2, ---m, and corresponding values for the surface elevation, h J=h(x J,y j)° Then, plugging each of these pairs into equations (2.4) gives 2 equations in a,b,c,d,e, and in this way we obtain a set of 2m equations in Sit+i5 unknowns. Write this system as F(a,b,c,d,e) s 0. (3.5) To simplify notation write He) in place of F(a,b,c,d,e), as long as no confusion will result. In practice the match points are obtained empirically by printing an uncorrected image of the region and identifying some recognizeable features. Because of this the points are subject to some inaccuracy, and therefore the model we produce using them is also. It is possible I . I 0 s I i l i. {-7 :33 (Wt ‘i ‘t' .' 73'!“ ii i l“ v‘ ‘ 0‘ ' O \! (V's ." il' . .«Y - .1. 331-1% x’19 ll 14 to improve the situation by having 2m ) SN+IS. Then we may not have an exact solution to F(a)=0, but we can solve for that a mich minimizes IIF(a)II2, and thereby lessen the effect of the innaccurate points. In this research two iterative schemes for solving this minimization problem will be investigated, the first is to approximate F by a linear function,say F(a)~F(ao) + J(ao)ea where J(ao) is the Jacobian matrix of F at a0, and so = a — a0, and solve for that sea which minimizes I|F(ao)+J(ao)AaII2 . This latter minimization is accomplished by solving the equation J(ao) aao = -F(ao) (3.6) as a linear least squares problem, giving aao = -J( ao)*F(ao) (3.7) where Jlsor = (Jtao)TJ(ao))"J(ao)T is the so-called pseudo-inverse of J(ao). Here, and throughout we will assume that J(ao) has full column rank. The object of this is that one then sets a, = a0 + see and hopes that ||F(°1)H2 ‘ ||F(°o)llz : (3-3) \mich is not necessarily the case, unless F is linear. However, if this is not the case, it is possible to decrease the length of aao in such a way that (3.8) holds, since: Liam); -J(ao)*F(ao) is a descent direction for H F(a)ll2 at so, if J(80)T F(ao) ilo. 2mm; Conside'r IIF(80)I|22 = F(ao)TF(ao), which has a directional derivative in the direction u of DuF(aO)TF(aO) = (J(ao)u)TF(ao) + F(a°)T(J(ao)u) = 2uTJ(aO)TF(80) bl svsri ion them av MM .3448 < mi Eli'H‘t/i‘li pd i".itibUHC’ mil svmqmt oi assn‘nmm {hiriw s isrit 'iol avioe I‘lfiJ a»: 'lJl'5 ,illztsii navitiiciu were no .iiiloll eit‘i'luzl’anrli ‘ilii in halls} :iti! H‘iniei l_li!"“1“|1'U 1111f: .hiiu)?“ ”Oilfiflminim suit gl'iluh'iu "lli' ""ill|‘lii‘iv' ‘}'vl!i'j"i‘i)l (fa-"V "TJ'WIS‘j’fi‘. 5.111 iii wouml 13 yd i sleilll'u’nqqs oi .21 lit-wit ml: ,bstepdeowm oil iii-e m’ii-imq i;f’vf,l1i}ii'.uli.ii 6A(DB)L i (palietel‘l 6 s=sAhne, .,.l 03 )r i in -1-iit.m lik’ldUJE-L slit .31 (unit. WIDTN’ .wll -§;I$/~.(|,5,‘L‘(I)U;IH :u' tantrum ii'izriw .6». lost ”ml fiVlflt'f bnr. l} 1'..«=?i‘1;‘v'-: “if“ priming all i infiziiiiniiiil 1 ' limit? i'i'illél :ll Hills! (55.?) (”(5)31- = ”8A (013 ‘i. ,_»~.-‘ 2;; ,"izsiiiulq ~ l v .v‘ u. >' (\ ‘ (Ulifbi’fl'ti); = ,‘hr‘. 1' \ i .4 '- v ' l - f , i . ‘08” Mr]; )- tab)‘ 1 ‘._n HmLJ 1”» 4w "" “ill/ind? i":!':8 .‘- 1-1? (”ti ‘1, :0 :J._ :5.--.'1i.- J!" ““T be? title “.1“ -t 1. ii, 'l ‘ all i ’* int: i :‘ifil‘lti i111}, ”1.5-. i ”L ~. it. “,3?! “ill: '1" . i-v Exil'fn >-‘.1 "it? .\‘) A ‘b‘vi' “ {i I} "if “*3 "'.‘l“:' T - ‘5, 1- ‘ --_5“ 1; --~ ~11 "it 't , . ' tat. :i‘w " ' i . ' l” .413 1 » v 1; ‘5‘ 1‘. .._, z'nli‘t,t). .. s 1‘ I i l .J!’ "' I ' " I ' ',‘i i- ’ll t“)'.:‘ . H : t .‘ Letting u = -J(ao)*F(ao), we see that (-J( ao)*F(ao))TJ(a°)TF(aO)=-((J( saints,»-lJtaoilrtaoiihtaoilrtao) = -<.l(a,)T Ftaolil ( J(ao)TJ(ao))"J(a°)TF(ao) < 0, unless ..l(ao)T F(ao) = 0, since (J(ao)TJ(ao))“ is positive definite. So the directional derivative in the direction given by -J(ao)*F(ao) is negative. Given this result, we conclude that there is some co>0 so that for a, = a0 + cocao with sea given by (3.7), we have IIF(a,)II2 < IIF(ao)II2. Continuing in this manner to find a2,a3,~--, we have an iterative scheme, and we can guarantee that IIF( ak)I|2 decreases to some non-negative limit. This scheme is sometimes referred to as a generalized Newton's method in the literature, since in the case that J(ao) is square nonsingular, we have J(ao)*=J(ao)", and so the above becomes Newton's method. Ben—Israel (See [2]) gives the following: Ihmmms Let f: E"-> Em, uocE", and r > 0, such that f¢C1(B(uo,r)). Let N, N be positive constants such that for all u, v in 8(uo,r) with u-v cRg(JT(v)), we have IIJ(v)(u-v)-f(u)+f(v)ll slillu—v“ II(J*(u)-J*(v))f(u)ll s NHu-vll and H IIJ*(u)II + N .t k < l for all ucB(uo,r) and IIJ*(uo)f(uo)|| ( (l-k)r. Then the sequence generated by "n+1 = un - J*(un)f(un) n = 0,l,' ' ' converges to a solution of JT(u)f(u) = 0, which lies in B(uo,r). 16 This illustrates that, theoretically, local convergence can occur without the decrease in step size discussed above. One should note here that quadratic convergence is not given in this theorem, as it is with Newton's method, and that numerically the scheme may behave poorly so that the result of the theorem does not hold. An alternate scheme mich is also investigated here is the negative gradient, or steepest descent method. Since -J(a°)TF(ao) is also a descent direction for ”H a)”2 at so, there is some c°>0 so that if we let .I 3 a0 - coJ(a°)TF(ao) then we have IIF(a')II2 ( IIF(ao)II2. Therefore, this method will also generate a sequence, “U such that ”H °k)”2 will decrease. In the following two chapters, the implementation of these schemes is discussed and some numerical results are compare d. ' I I1 5: IJ I " 1 (it!) ’. lllU" I 'ith'ti' n ,;l! .t -! 215' k", : '1ti2.l . 3.7“ 12*in 51.3 r, stic- iii-I'M.- - in .D - 1- : 1‘. sl" . . '. i"~ q}: I I l" C u f a "‘ IOIQ'U D "t' .l‘l’ . . ."' 1" ‘. )i"! “:| ,7 0 .‘ ' II '3: -.f in. «is . w " I e"' III 1 .‘H" Chapter 4- The Algorithm This description of the algorithm will be divided into four parts: 4.l: Input and Initialization 4.2: Computation of F and J 4.3: Solving the Linear Least Squares Problem 4.4: Choosing Step Size and Updating The final section of this chapter deals with possible variations to this algorithm. 4.1: Input and Initialization The first task is to set the parameters of the problem, N, «m, Act , P=pixels per scan line, L=scan lines per second, and an initial approximation of the speed s, altitude H, and the yaw 90 of the aircraft in flight. In order to initialize the cubic spline functions, assume level, straight line flight in the direction of the x-axis, at the given speed, altitude and yaw X(t)=st gives em = isd, an=s, a,2=ai3=0 Y(t)=0 gives bifo j:0,l,2,3 Z(t)=H gives ci0=H, and c“: i2=ci3'° 9(t)=0 gives difo j=o,l,2,3 and e(t)=eo gives eio=°o: and en=ei2= i3=0, i=l,2,- ' ' ,N. Next, input the set of match points uk,vk,xk,yk,hk, k=l,2,...,m, where hk=h(xk,yk) is the elevation at (xk,yk). At this point we need to determine the interval containing each match point, and so If t I (UK+VK*2%/flp)/L C [ti'l’til l7 define Ak = t-tl- and set xx = X(t) = aio+an Ak+8iz Akz+ai3 A? = °i0*°il AK Yk = 0 2k 3 H 9k z 0 9k a so Now compute the scan angle for each match point GKBM4‘VKA“, and for the purpose of later computing J(a) it is also desirable at this time to compute "k 'fl . fl and fl .. ah. an. ad.. ae.. l.) ‘J U l.) ‘3 These are in fact, all the same. For example, (see equation (3.3)): Xk = (1 AK Akz Aka) .i (4.1) 3 (1 AK Akz Ak3) Aifl 2 3 . and, Yk :- (i 4k AK AK ) AIb, etc.. Therefore, the row vector q(k)= (i Ak akz Ak3)Ai gives the partial derivatives of xk, as well as those of Yk,2k, 9k» and BK. This calculation can be made fairly efficiently by noting the form of A,- from equation (3.4), and observing that: (l s s2 3’) "3"4 = (33-(3-1)3)d3+(6(J+l)-3)d ZA-Sd s2 - (l s s2 s3) 1134 +6uzts -j d) j-l,2,-- ' 1 x \i V- 't ‘ . ll '- 1,. .06 .. ‘6 'lilUt') 11.22r1m .117" ’1'? 'lw‘l v'lL." it)? :livi.'tqiilt".t my 1-:‘T 3.. i} Hit-2%“ «)rif " 3. (b)! " "1' l L)» ~ '» Juliunrlf Win) .Ttli‘ihfilfl mi ,- - -H ‘t If I l A“ z ‘ .; ('i.) I' , f\ _ Li ’ \ 1 .i ‘. ' -- i , ,, , . ' ,d V t » ‘ '. 1‘ ' it"? 1 '1 a ’ 2 I g \ .1117 ii. i’ ”mi r1( A (0., .I “ “3.57.1 ' "l' -."’-’)"‘“'g' l t i . . , \ \ ~ ,_ xv , y r i.« i '1)‘ ¢ . ' u' '. . l iii-1 "“ . ‘i‘ tztijlff‘ .i: ll: ' 63-.) 1. tit- / t . ' - .ri. ‘ .r ‘ ,. , _‘ l. r i i \r i, . 1 = . .t <_ -,_ . . . were the subscript k has been suppressed. Thus, to compute q = (l a a7- A3) Ai; set qJ- = 0 , j=l,2,"‘,i-l q, = 43 qi+1= (A'd) 3‘A3 qi+i=qi+ 60'2“ - .id) . rinse-2 flu = ' qN+1= A ' (N-i)d qu+2= qu+12 qN*3= (IN-1+ 6d2(A ' (ti-I)!” 4.2: Computation of F and J The computation of the Jacobian at each step is a fairly large part of the total computation, but can be accomplished quite efficiently using the partial derivatives computed above. As an example consider the calculation of the partial derivatives of the ri-m equation, which is either, Fn = xk-Xk+(hk-Zk) M with k = ("m/2 0,. Fn = yk-Yk+(hk-Zk) vgk with k = "/2 depending on whether n is odd or even. Then, aFn = 3F“ axk 3813 an 38.0 ab” aVk fin- aFri = 3':n 32,, BC” 32k 3C“ .Jt\ .- i4 4- , ... I \e .s l . l i . .1 MW . luv. l a 9 f. i I a Q 0 ‘1. 1' s I l s a 4 I .V .II t, I- h o «1‘ .e' e le 39k alle is = aFn s, RU ask EU and as discussed above the second factor in each of these is the same constant, and therefore does not need to be re-computed at each step. The partial derivatives of Fn with respect to Xk, and Yk are -l and 0, or 0 and -i respectively, depending on whether it is odd or even. So we only have three factors to compute for each n, those being the partial derivatives of Fn with respect to 2k, 9k» and ek. These are computed now along with Fn- For K=i,2, I'll, given Xk, Xk, yk, Yk’ hk’ 2.0 9k, 8k, «k, set: ka= -cos txkcos 8k sin 'Pk + sin otk sin 6k -cos “k cospk (4.2) VU": -cos otk sin 9k sin 9k -sin “k cos 9k -cos “k cos ’k and set F2k-l‘ "k — Xk - (h.< - 2k ) ka (4.3) F2k=gk'yk-(hk'zk)vgk and aFZk-l = M 32k 36 - 2k - v9k 21 aF2k_‘ _._ -costxk cossk + Sifltxk sing),< sinek (hk'zk) 39k 6031ka c0329),( 52k = -sine kcosotk - sinotk sinpk cosek (hk‘zk) 39k cos “k cos 29k aF2k_, _._ sintxk cosek + cosotk sinpk sine,< (hk'zk) 33 k cos txk cospk 3':th = sintit.< sine.< - costxk sinpk cosh.< (hk'zk) 39 k cos etk cospk Using these, along with the vectors q(k) computed above we can now assemble the Jacobian matrix. In doing this we may arrange the rows and columns in such a way that that the Jacobian will have a block upper-triangular form, Mitch will be advantageous in solving the linear least squares problem. To do this, first arrange the columns of J according to the following order on the unknowns: a13 biS °13 ‘113 ‘313. 823 b23""31l-ls 3N0 “no 6no “no eno am‘ ' ‘ ‘3113.- Next, order the match points according to increasing uk, with the corresponding arrangement of the equations obtained from (2.4). 4.3: Solving the Linear Least Squares Problem It is not desirable to solve for one which minimizes II F(ao) + J(ao) ea ”2 by calculating no = —[J( ao)TJ(ao)]"J(ao)TF(ao) since a lot of calculation is required (J(ao) is 2mx5N+i5), and if _ v. . l v t i i .. .. 1 . . I \ . ’l . i . . l l J . i . . r: .. . ' . . I ‘ r o i ‘P 1' r i . . _ U ' n I e ‘0 f (I r I I ' . O as . . . 2. . 1. . it . -l. l. ah .l ‘9 22 J(ao)TF(ao) is nearly 0 this may behave poorly (See [ill,p202). Consider the general problem; Solve AX = 8 in the least square sense, where A is mxn, X is nxi, and B is mxl. Factor A as A = 0R where 0 is an mxn matrix with 070:1“, the nxn identity matrix, and R is nxn upper triangular. If A has full column rank then R is nonsingular, and it is well known that the least square solution of AX = B is the unique solution of RX = 073, which is easily obtained by back substitution. This method of solving the least square problem is advantageous for two reasons. First, the matrix J(a) as discussed above has a form which is block upper triangular and so the amount of work required in factoring J(a) is not nearly as much as that of factoring a full matrix. Second, if J(a)TF(a) is nearly 0, then the problem is ill conditioned. Since we seek a minimum for H He) ”2, at which point, necessarily J(a)TF(a) = 0 , this problem cannot be completely avoided. In this case the solution method discussed here is advantageous because 0 is orthogonal, so RTR = RTOTOR = AIA. Therefore Cond(R) =r/;(RTR) 9((RTR)") =fpthThl pttATAr‘), whereas Cond(ATA) = MATA) pttATAr‘) = Cond(R)2 Suppose that at stage 3 we have computed A J where A=0jA-, and A J has only zeroes below the principal diagonal in the first j—l columns, and .= .1 let l3J OJ 8. Let AJ+1=EmJ‘ ' ' EJ+1JAJ’ where .y’ .t (I x... will .“‘-t . 43 .2; . l l" 1. t’i . T . —\é’i .. (.2 .0. i 10. ii .' ' .vl cl .5]. I “in, a p. o o . o1 13-1 = . 0. CJ'J' cm 0 ; o o EIJ ' - ‘ I1-3-1 . . ° °ll '°JJ ° 0 o ‘ Im-i-l o o o , L ._ lJ _Jg' JJ _gg’ IJ IJ , K K and where an and a” are the ijtm- and ”“1 elements respectively, of Ei-lj "' Ej+lj AJ Then each EU is orthonormal, and A3” has all zeroes below the principal diagonal in the first jcolumns. Set ...T . .T... .T "3+1 ' OJ EJ+1J EJ+2J Ema and we have A a: 03+) A)”. Repeating this for each j-l,2,---n, we obtain finally, A ' "n+1 An+i where A has the form R [.3 and R is nxn, upper triangular. We may partition amt as on” = l 0 = U l with 0 being mxn, and we have the decomposition A = 0R. In implementing this method, the matrix 0 is not explicitly computed, but rather the rows of A and B are operated on so as to sucessively reduce the columns of A to the desired form. These row operations being equivalent to multiplication by one of the Eij's. The following simple FORTRAN routine accomplishes this: 24 DO 30 K=l,NC DO 20 I=K+I,NR IF(A(I,K).LE.TOL)GO TD 20 SIM K,K)/SQRT(A(K,K)“2+ A( I,K)“2) T: A( I,K)/SQRT(A(K,K)I*2+ A( I,K)Ii2) 00 l =K, A(l<,.l)=s:u+v:l A(T,J)=SsV-UIT V=B( I) 8(K)IS*U+ V!“ B( I )=5IV- U‘IT 20 CONTINUE 30 CONT INUE At this point the array A(I,J) contains the matrix R and several zero rows, and array 8(1) contains OTB in entries Isl,2,°"n. We now solve the equation RX s 078 by back substitution for X. At this time we may also compute the condition number of R as Cond(R) = IMX{A(I,I)=I=1,2, "'n} HINMI,I)=I=L2, "n} 4.4: Choosing Step Size and Updating Having solved the linear least squares problem, we obtain ea from which we must calculate AXk for k=l,2,'",m, with which to update Xk. There are at least two ways to accomplish this. First, Xk is the value of X(t) = gaiju'ti” for some i, at the point t = (uktrvksZotu/fiPflL determined by the k-M:l match point. Using the notation of Section 4.l above: xk = °i0*'il Alttiliz Altzft'is 41:3 = (I Ak Akz 133) a,- So we may compute aai by so, = Aiaa (see equation (3.3)) i i)fi)TCuL3\{l,t)n" i .\ -.'.:i 1’ vii ._ ll! 1..” xiii .,i)a=tl [M 3 l )l’: _ft. )1 _- 'l‘)‘l l'liil ‘tlt.'.'. ll; I { tgl! .li: ' I '_ 2" If‘Ltt: 1t" 3 i do 'ltl .fl‘u ;.i.' 1”?) i"! U ill ' at; .l r. '1 . 3.. .. ililiw' .l', .l 1* i =‘iw (I '4. v\. at 25 and then set AX = (1 AK Akz Ak3 )Aai. Alternatively, we may make use of q(k) = (1 AK Akz 15,30 A,, which has been computed and stored as discussed above, and set an = q(k)aa for each k. To count the number of multiplications required by the second method, recall that the vector q(k) has i-i leading zero elements, so that if we take N/2 as an average value for i then we must perform N+3-(N/2-l)=N/2+4 multiplications to compute q(kMa. Thus, this method requires a total of approximately m(N/2+4) multiplications. The first method described above, requires that one multiply Aiaa for each i=l,2,“-,N, which would take about 4(N/2+4)N multiplications, and then multiply one of these by (l Ak Akz Aks) for each k=l,2,~--,m. So this method uses approximately 4(N/2+4)N+3m multiplications. Comparing these two we can see that the second method may be preferable if m(N/2+4)< 4(N/2+4)N+3m that is, if m(N/2+l)< 4(N/2+4)N. In the numerical experiments which are discussed in the last chapter, we are working with about 60 match points, so the condition becomes 60(N/2+l)< 4(N/2+4)N or 2N2-l4N-60 ) 0 which is satisfied for N > to. In that last chapter we will compare results for various values of N, some larger and some smaller than to. For now, let us proceed to describe the algorithm using the second method for updating. For k=i,2,---,m , set 26 AXk qkaa qkab szk :1 qkAc Apk :2 qkad ask := qkae .. ll AYk For n=0,l,"-,l5, set c=2’", and For k=l,2,"-m set - Xk(n) := Xk + c AXk Yk(n) := Yk + c Mk Zk(n) = 2k + C Azk 9km) -.= 9k + c A9k ek(n) 1: 9k + c ask and compute FZk-l and F2k according to (4.2) and (4.3). FNORl‘l(n) = FNORI‘l(n) + F2k-1**2 + sz“? Choose the n for Vnich FNORfl(n) is minimum, and if this value of FNORli is greater than that obtained in the previous iteration, then STOP. Otherwise, proceed to update: For k=l,2,-"m, set Xk== Xk(n) Ykfl Yk(n) 2km 2km) 9k“ 9km) aka ek(n) a + The D + 2'"Ab c + That: a + 2'"sd II a + 2'"ee and also 0508’ At this point, to continue the iteration, return to the routine for 27 computing the Jacobian matrix. Section 4.5: Variations The algorithm for solving this problem by the steepest descent method is identical to this except that the subroutine that solves the linear least square problem is replaced by a subroutine to compute the product x = ATB. Next, note here that the system (3.5), obtained from equations (2.4), is linear in a, b, and c, so that if we fix d = 110 and e = e0, there is a unique solution to min HF(a,b,c,do,eo)II So that it may be desireable to solve this smaller problem as an initial step before proceeding with the full algorithm. Another variation, suggested by numerical experience, is to pose some additional endpoint constraints on the splines. We introduce the so-called natural boundary conditions: x‘(0) = X“(T) = o, (4.4) and likewise for Y', Z", 9', and a“. In addition it is advantageous to add conditions that will control 9 and 9. Neither of these should be large, nor should they very rapidly, and the conditions that have worked best are 9'“) = 0 and 8'“) = 0 t: l'T/N, i=0,l,"‘N. (4.5) If these derivatives are close to zero at the nodes, then the splines cannot behave too badly in between. The additional equations obtained from (4.4) and (4.5) are combined with the equations that constitute the system (3.5), giving us a system ls -l I t. -41: 5 1.. Itil‘ ‘. I .1. e la) a {1‘31 4.1 |’- [:4 {if i =1! ' i'lrlj’ll 7, 12%;- 3 3.3‘11‘ . L 9...: I HI .V'l a 28 of 2m+10+2(N+1) equations in the 5N+15 unknowns to be solved in the least square sense. This is a reasonable approach since we expect these conditions to be nearly satisfied, but not necessarily exactly satisfied. Numerical experience shows that some additional constraints of this type are necessary, since without these the solution may exhibit wild behavior at the ends, such as an aircraft plunging from an altitude of 50,000 feet, to 10,000 feet below sea level (15,000 feet below ground level) in the final seconds of the scan, and spinning several times in the process. A final variation to be discussed here, is the possibility of using weights. Given a diagonal matrix W = diag(w1 w2--~wm), one may solve WAX = W8 in place of AX = 8. According to probability theory (See [10], p77), if one wishes to solve the equation AX = B where the elements b,- of B are measured quantities subject to errors that are independent, normally distributed, and have mean 0, then taking the wi's to be inversely proportional to the variances of these errors, the solution of WAX = B in the least square sense gives the most probable values for the elements of X. In the problem that we are considering here we have no knowledge of these variances, so some numerical studies will be done using various weights for the equations in each of the three groups, those obtained from (2.4), from (4.4), and from (4.5). Chapter 5: Numerical Results For the purpose of numerical experimentation, data was obtained from GeoSpectra Corporation of Ann Arbor, Michigan. The data consisted of 62 match points, giving row and pixel numbers for data points, and corresponding longitude and latitude for map points. This author transformed the longitude and latitude coordinates into a rectangular system in units of 1000 feet, and also determined the corresponding surface elevation from a US. Geological Survey map of the area, with contour intervals of 40 feet. Finally the aircraft yaw, so, was approximated by visual inspection of the plotted match points. The scan was taken from an aircraft flying at an altitude of approximately 40,000 feet above sea level. At this altitude, the field of view for the scanner when it is pointing straight down (See Figure 1, p3), is a circle of diameter approximately 100 feet. The computations described below were performed on the CDC Cyber 750 at liichigan State University. Programming was done by this author. using the FORTRAN programming language, and a full listing of the program is given in the appendix following this chapter. In this discussion no will be the average over a set of match points, of the Euclidean distance between a given point (x,y), and a point (x',y") obtained by use of the model under discussion. The values given for HD will be in units of feet. Upon entering the data and running the program in various ways it was observed that certain points made quite a large contribution to the value of liD. Two points seemed to be particularly bad, so computations were performed with a data set reduced to 60 points. The result was a 29 .' . a \ ”Ti ..- Ii 1),]? mitt-m . i. , -: Mi ‘. ‘ : ‘ - M ii . 1'5 iii?‘ ~ , 21w: - Ali“. 11:1. 31".? ‘ilJ"? :l~ -114.) It). 1 113-“ 95f . til": i; 01.: is, r. if: 1.11'ig‘fl'i’ t-i ‘ 1:51 titli' lii"-.*1ttilllrli1.111111 lit - it 1‘11: triul'ilil'i-iq '31: will .tlLi iii? it". 1:110): t pm. :11 new» 1. 3} 3111113111 ! "~13 limit i ,'.?1)""} i "l L, {if} ‘ vii 1‘ 1“ ~ l1.<_t:it"lf lint 11.11 t. till». ,-.ii'-.i.~t.tiim.: .14,“ thug-.15.: 1‘1‘1'1'. first) .311. Mirna. tum 11,131 "i" l (Jil‘liiflt. 'Ili‘f ffl' .\ ‘ 1‘) ' 1 1" 'ii‘lVg. ”Flt-’3 ' ' .“. ‘ O t ‘0 I " i.‘ V :1 .‘ 1 . I i1 i 01'“; 1. ' ‘ it '9 1;." 31.: : 3.11.?" “it. I t‘.“'qt.‘( ‘1' I: .1'. sir! ‘é' I.o f )l .l..~.. ! 18,! 3.1." xii .li Hg. .trfl' '11 ‘:_'1l 1. ‘. l 3 23% t .H‘,‘ ”'11 .111 '-ii {1' ;‘.:1:t171~.-'1 . - :-‘ :- 3-33411‘1 L15}""EI :1." n. 0’. ‘I‘1E‘i1l v l 1'. it) ' I 1;": { I‘e" 01 .I a‘vll J -')‘1 (”it it“, til. 41211111114111)“ :1 ii i 'i {'1 i \ "1 I it} )1) 6 limitliltl‘! t Vié-“H‘JHI ‘11).) 1 1 i;_-. . 1 m; 1 '5‘ 11.1}. I; 1! .‘r.it l; ,: “i 1’ '. 131:. “1:1“ it ; - .. : - . 11.. 113’. . &:2; 5181111 Lil‘qulb [1,“ t l (H :1 .. , ll)! '3 3‘1 30 decrease in the value of ND by about 16.5%. For this reason it was decided that the two points would be removed from the data set for all of the following computations. One should be cautious in removing a point from the data set since the removal of any point would likely cause some reduction in the value of ND. For example, if this were carried to extremes, and the 60 points were removed and only the two remained, then one could easily construct a model for which no = 0. This would not necessarily be a good model. Example 1: The use of weights In this example we will experiment with various weights as discussed in Chapter 4. We have 3 groups of equations, those obtained from (2.4), the boundary conditions arising from (4.4), and the conditions obtained from (4.5). We restrict our investigation here to assigning the same weight to each equation in any group. Let us assign the weight of i to the equations in the first group, W1 will be the weight assigned to the second group and W2 will be the weight given the equations arising from (4.5). In TABLE 1 we see the results using various values of W, and W2. Column 4 gives the condition number of the matrix R obtained from the OR decomposition of the Jacobian matrix as described in Section 4.3. The last column gives the number of iterations used, where the stopping criterion was that IIF(a)Il2 does not decrease after the step size has been halved 15 times. 31 TABLE 1: The Use of Weights W1 W2 flO Cond( R) ‘ of iterations 1 .1 70.2 809.58 6 1 i .0 70.3 797.09 1 5 .1 10.0 71.4 767.98 9 1.0 .1 66.7 431.60 26 1.0 I .0 65.9 377.18 10 1.0 l 0.0 73.0 466.79 2 10.0 .1 67.2 558.37 26 10.0 I .0 66.0 455.94 1 0 I 0.0 1 0.0 73.0 500.80 3 This data seems to suggest that the best model is produced with W, in the range [1.0,10.0), and W2 in [.l,i.0]. Also it was observed that after a few iterations, very little progress was made in most cases, so that one might wish to choose a different stopping criterion. The difficulty in making comparisons here is due to the fact that 110 is not an absolute test of the quality of a model. The true test is how well the model will plot all of the points in a data set obtained during a scan, not only the relatively few given match points. In order to analyse this situation more quantitatively, consider an imaginary flight path given by: X(t) = 575t + 12 - 100 sin(t/7 - 3) V(t) = 1200 sin(t/4) ~300 cos(t/12) 2(1) = 35000 + 1000 cos(t/5+l) 9(1) = .04 sin(t/3) +03 cos(t/7) e(t) = .14 + .05 cos(t/3 - 5) for tc[0,34]. Two sets of 60 match points each were created using a random number generator to select row and pixel numbers, and ground elevation. The first of these sets wil be used to construct a model by the algorithm described in Chapter 4, and then the second set will be used to test 32 that model. In the following, 110, and 1102 will refer to the value of 110 (as discussed above) obtained from the first and second data set respectively. Repeating the experiment with weights as above, we obtain TABLE 2. The number of iterations before stopping was 3 in all cases except where W: .1, W2=10.0, in which case 4 iterations were required. The condition number at the stopping time, and the ratio 1102/1101 are also given for comparison. TABLE 2: The Use of Weights (continued) W1 W2 1101 CODd(R) "D2 1102/1101 .001 .001 7.15 286.70 53.93 7.538 .01 .01 7.87 157.27 42.55 5.407 .1 .1 7.69 43.761 13.90 1.809 1 1.0 9.52 43.733 16.14 1.695 1 10.0 229.49 42.871 280.24 1 .221 1.0 .1 17.91 27.409 35.99 2.010 1.0 1.0 9.03 27.566 17.20 1.905 1 .0 10.0 42.23 25.395 78.24 1.853 10.0 .1 24.95 142.44 48.87 1.959 10.0 1.0 1 1.84 143.08 20.21 1.707 10.0 10.0 44.30 144.73 73.58 1.661 100.0 100.0 438.12 1461.6 506.95 1.157 The first thing that one will notice in this is that the condition number seems heavily dependent on W1 and nearly independent of W2, and that it does not seem to reflect the accuracy of a model. Compare for example lines 4 and 5. Also, it seems clear that one should avoid large values for W2 in particular. The optimal values for both W, and W2 here seem to be between .1 and 1.0, with the best results coming when they are the same. 33 Example 2: Variation in the number of intervals As discussed in Section 3.1, the accuracy of a model is 0(d4), where d is the length of the intervals, so increasing the number of intervals should improve the result. 0n the other hand, with more intervals come more unknowns, and the possibility that the Jacobian matrix will have less than full column rank. Let us examine this situation now. The calculations for TABLE 3 were performed using weights w1=w2=.s, and varying the number N, of intervals. Computations were made for both data sets, and the additional data set was used for comparison. Also given in the table are the number of iterations, and the condition number at the stopping time. TABLE 3: Variations in the Number of Intervals Actual Scan Hypothetical Scan N 11D “ I TER CONN R) 1101 HDZ ‘ I TER (JOHN R) 2 1 02.05 6 l 48.51 675.24 854.53 7 268.68 4 76.1 9 7 352.08 226.25 286.46 7 507.81 6 66.73 16 384.85 29.23 49.44 3 671.88 8 66.09 1 1 367.27 12.41 22.38 3 825.59 1 0 67.01 1 5 650.90 6.91 1 4.1 4 3 880.59 12 68.22 20 602.30 5.02 1 1.1 0 3 995.23 14 67.26 1 1 1976.90 4.63 l 1.29 3 l 1 19.4 In both cases, for values of 11 greater than 14, there were intervals with no match points. The best choice for N appears to be 8 for the actual scan data, whereas for the hypothetical case it is 12 or 14, depending on whether one considers HDZ or HD1. all. 7 34 Example 3: Comparison with the Steepest Descent liethod In this example we compare the results above with results obtained using the steepest descent method for iterating. We started with the initial approximation discussed in Section 4.1, and fixing 9 and 9, updated X, Y, and 2 using the generalized Newton's method. As discussed above this is a reasonable first step since the equations are linear in the coefficients of X, Y, and Z. We then proceeded to iterate using the steepest descent method. Partitioning the scan into 10 intervals, using the actual scan data, we obtained a value for MD of 70.25 after 25 iterations, compared with "0:67.01 after 15 iterations using the generalized Newton's method. In the case of the hypothetical data, using 10 intervals, after 16 iterations, 1101:39088 and H02=S3l.l9, compared with 6.91 and 14.14 respectively, after 3 iterations using the generalized Newton's method. In conclusion, the scheme presented herein, appears to work quite well, with rapid convergence. The values obtained for "D are quite reasonable considering that the field of view at all times has diameter at least 100 feet, as discussed above. The discrepancy between the results using the "real" data and the other may be only partially due to the input error. Another factor may be that all of the given match points were in the right half of the region being scanned. APPENDIX 35 APPENDIX: Program listing 90 66 666 100 PROGRAN NSSPT DOUBLE PRECISION ABC(108) CONNON 11PS ,211P( 20) C011110N /BONE/U(62) ,V(62), X(62) ,Y(62) ,Z(62) C0111‘10N IBFOR/ REAL OHEG(6Z) ,PHI(62) ,KAPP(6Z) DATA 0866/62I.0/,P11I162s.0/,KAPP/62a.1 4/ DATA ABC/108*0.0/ OPEN<22, FILE-' A ) READ S,T CALL ABIN(N, 1100, T) REWIND 2 I'M-0.0001 THNAX=0.7S P8716 "00:3 NC-NO0I(N+ 3) DO 7 1:1 READ(2,1)UH(I),V(1),X(1)V,(l),Z(I) CONTINUE NTCaO CALL EXIN(P,N,T,TH11AX,011EG,PHI,KAPPJlOD) CALL JACOB(ONEG,PNI,XAPP,NOD,N) NR=2i11PS+6 CALL CUEAHMNRIS NC ,CHG) CALL ABUP(N, 1101) CALL EXUP(P, N, T, THliAX, 011EG, PHI ,XAPP ,1100) 1100=5 NC=5IN+15 =2I(11PS+1100+N+3) CALL JACOB(011EO,PHI,KAPP,HOD,N) PRINT 100 CALL CUEAHR