the preferred embodiment is best understood by referring to the figures wherein like numerals and symbols and used for like and corresponding references to the various drawings . the method of the preferred embodiment permits determining the attitude ( yaw , pitch and roll ) of an imaging sensor relative to its attitude at a reference time . inputs are images at a first time and a second time and relative sensor position . the relative sensor position may be provided by an inertial navigation system ( ins ) or other source . the method meets the needs of passive ranging methods and systems where accurate estimates of attitude are crucial to the missile or platform performance and are difficult to obtain , particularly when the sensor and ins are not co - located . for many of these systems , however , it is not difficult to obtain estimates of relative translational motion to the required accuracy . the preferred embodiment is best understood when considering an example in which certain assumptions hold . the first assumption is that the sensor is in or on an aircraft , such as a missile , flying over terrain . the example also assumes that image features are automatically evaluated or selected for robust tracking and that image trackers continually relocate features within the field - of - view ( fov ). furthermore , the example assumes that new features are introduced to replace those lost through motion . finally , the example uses an assumption that features derive from stationary objects in the scene . with these assumptions , the explanation of the preferred embodiment proceeds as follows . feature locations are reported by the trackers as points in the image , although the features themselves are distributed objects . interpolation techniques can provide sub - pixel resolution , but the example of the preferred embodiment assumes tracker outputs are reported as pixel locations . the image of the tracking sensor typically is in a world coordinate system , as depicted in fig1 . understanding the preferred embodiment becomes simple , however , in a canonical coordinate system that is obtainable by transforming world coordinates . fig2 shows such a canonical coordinate system . in fig1 the center - of - focus ( cof ) is located at point p . a line from p to the fixation point f passes through the center of the imaging array and is normal to it . this determines the direction of the unit normal to the array , denoted n . the center of the imaging array is located at p + ηn , where η is the focal length . positions within the imaging plane are defined by the 2 - dimensional coordinate system ( u 2 , u 3 , where from the perspective of the observer , u 2 is associated with horizontal and u 3 , with vertical . the canonical coordinate system of fig2 may be obtained by a translation and rotation of the world coordinate system of fig1 . this transformation leaves invariant distances such as focal length . note that fig1 and 2 employ similar nomenclature , with primes distinguishing quantities in the canonical system . in fig2 the cof is at the origin and the fixation point is on the x 1 axis at x 1 &# 34 ;=| f - p |. axes u 2 &# 34 ;- u 3 &# 34 ; are offset from , but aligned with , x 2 &# 34 ;- x 3 &# 34 ;. the transformation from world to canonical coordinates consists of certain simple steps . the first step is to translate the origin to p . the system is then rotated by a yaw angle φ 3 about the x 3 - axis so that the fixation point is directly below the x 1 &# 39 ;- axis in the resulting intermediate coordinate system x 1 &# 39 ;- x 2 &# 39 ;- x 3 . next , the system is rotated by pitch angle φ 2 about the x 2 &# 39 ;- axis so that the fixation point falls on the x 1 &# 34 ;- axis in the resulting intermediate coordinate system x 1 &# 34 ;- x 2 &# 39 ;- x 3 . the system is then rotated by a roll angle φ 1 about the x 1 &# 34 ;- axis so that the positive u 2 &# 34 ;- axis aligns with the positive x 2 &# 34 ;- axis of the resulting final coordinate system x 1 &# 34 ;- x 2 &# 34 ;- x 3 &# 34 ;. where 3 - vectors x and x &# 34 ; are locations in world and canonical coordinates , respectively . if x represents the location of a feature , it is by assumption , constant . however , its mapping into canonical coordinates , x &# 34 ;, is time - varying due to sensor motion . the 3 - by - 3 rotation matrix s is given by : ## equ1 ## where c i = cos φ i and s i = sin φ i . an important property of s is orthonormality , that is where s t denotes the matrix transpose of s . this property makes it easy to transform between reference frames . the euler angles φ 2 and φ 3 may be related to the unit normal n in world coordinates using the special case of 1 : ## equ2 ## to solve this , transpose and multiply on the right by s . from this expression , identify : where ## equ3 ## these results are useful when θ 1 and n are used to specify attitude . the tracking sensor system assigns to each feature within the fov a pixel location in the imaging array . this mapping is independent of choice of coordinate system . the physical center of this pixel is taken as the location ( u 2 &# 34 ;, u 3 &# 34 ;) of the feature in the image plane . it is appropriate to assume the imaging array is rectangular with n r rows and n c columns . pixels are identified by 2 - tuples of the form ( j , i ) where j is the column index and i is the row index . the remaining details of this mapping are contained in ## equ4 ## where δ 2 is the column spacing and δ 3 is the row spacing . fig3 and 4 are illustrations of very simple imaging arrays . referring now to fig5 by similar triangles : ## equ5 ## note that ( η , u 2 &# 34 ;, u 3 &# 34 ;) and ( x 1 &# 34 ;, x 2 &# 34 ;, x 3 &# 34 ;) may have different units of measure provided each system is self - consistent . as an expedient , we take η as the unit of measure in the image plane , such that ## equ6 ## if u &# 34 ; in equation ( 14 ) is the mapping of x &# 34 ; in object space , by equations ( 12 ) and ( 13 ) the effect of motion on images collected at two times , 0 and 1 ( in appropriate units ) is an important parameter of the preferred embodiment . the derivation is simplified by picking the coordinate system to coincide with canonical coordinates at time 0 , yielding using equations ( 1 ) and ( 15 ), at time 1 the following holds : the first row of equation ( 17 ) defines x 1 &# 34 ;( 1 ) and may be used to eliminate it in the second and third rows , with the result where s i i ( 1 ) denotes the i th row of s ( 1 ). substituting from equations ( 15 ) and ( 16 ) for x , by eliminating x 1 between equations ( 22 ) and ( 23 ) and simplifying , ## equ7 ## where u 2 and u 3 are shorthand for u 2 &# 34 ;( 1 ) and u 3 &# 34 ;( 1 ), respectively . hence , the determinant in equation ( 25 ) vanishes when motion parameters p ( 1 ), φ 1 ( 1 ), φ 2 ( 1 ) and φ 3 ( 1 ) are consistent with the observed feature locations u ( 0 ) and u ( 1 ). however , the solution of equation ( 25 ) is underdetermined , since there are three unknowns . this is overcome by adjoining equations for other features . adopting a least squares approach , it is practical to define the nonnegative cost function for k features ( k ≧ 3 ) as ## equ8 ## where δ k is the determinant for feature k . the method of the preferred embodiment , thus , is to determine the euler angles which minimize the cost function j . the problem of minimizing equation ( 26 ), however , is nonlinear and must be solved iteratively . newton &# 39 ; s method or other gradient - based approaches can be applied as analytical expressions for the derivatives can be derived . however , efficient procedures are available that avoid the need for differentiation , for example , powell &# 39 ; s method as described in m . j . d . powell , &# 34 ; an efficient method of finding the minimum of a function of several variables without calculating derivatives ,&# 34 ; the computer journal , 7 , 155 - 162 , 1964 , and g . r . walsh , methods of optimization , john wiley and sons , 1975 , pp . 129 - 138 , provides an applicable optimization technique . an illustrative embodiment uses a two - stage search that the following discussion describes . testing of the preferred embodiment has been performed with synthetic data . the minimum of the cost function j was sought by searching the three - dimensional solution space about a nominal set of euler angles , such as might be provided by an ins . experiments were repeated for several levels of measurement noise . the resolution and fov of the sensor will affect performance . the sensor model of the tests assumed an array of 480 rows and 960 columns of light sensitive elements or pixels . each pixel is square and subtends an angle of 0 . 1 milliradian . this corresponds to a pixel dimension of 10 microns and a focal length of 10 centimeters , for example . the horizontal fov is 96 milliradians or 5 . 5 degrees . the narrow fov makes it difficult to maintain a feature within the fov for long . the pose of the sensor at times 0 and 1 also affects performance . the locations of the cof at times 0 and 1 were chosen as follows : ## equ9 ## attitude was specified as a roll angle φ 1 and the normal n to the imaging array . at time 0 , roll was assumed to be 2 . 5 degrees and the normal , ## equ10 ## at time 1 , roll was assumed to be 5 . 0 degrees . the normal was not specified directly , but was computed by requiring that the fixation point be identical to that at time 0 . in this calculation , fixation points were defined at zero attitude . computing n ( 1 ) in this manner increases the probability that a feature will be visible at both locations . a contributor to measurement error is quantization of ( u 2 &# 34 ;, u 3 &# 34 ;) to the discrete pixel location ( j , i ). while sub - pixel interpolation may reduce this error , the tests of the preferred embodiment do not model sub - pixel interpolation . to model other tracker errors , gaussian white noise was added to the feature locations before quantization . the noise had zero mean . the variance was experimentally controlled . feature locations in the imaging array at time 0 were picked stochasticly using a gaussian random number generator . the nature of the distribution caused feature candidates to cluster near the center of the array . a ray from the cof was struck through the center of the pixel and projected until it reached an attitude h , another random variable which modeled the topography of the terrain . this point was then back - projected to the cof at time 2 and the intersection with the imaging plane noted . the location of the intersection was corrupted with noise and quantized to a pixel location . if the location was within the fov , the feature was accepted . if not , it was discarded and a new candidate was processed . candidates were processed until eight features had been defined . since the coordinate system at time 0 is not canonical , the derivations of the preceding section must be modified . fortunately , it is only necessary to replace s ( 1 ) with s ( 1 ) s t ( 0 ) and p ( 1 ) with s ( 0 )[ p ( 1 )- p ( 0 )]. table 1______________________________________results ( deg ) with resolution = 1 . 0 deg . rms noise ( pixels ) euler angle 0 . 2680 0 . 6867 1 . 7340______________________________________roll true 5 . 00 5 . 00 5 . 00 estimated 4 . 10 4 . 20 4 . 40 scatter 0 . 46 0 . 60 0 . 60pitch true 70 . 10 70 . 10 70 . 10 estimated 70 . 00 70 . 00 70 . 00 scatter 0 . 04 0 . 04 0 . 05yaw true 93 . 80 93 . 80 93 . 80 estimated 92 . 50 92 . 60 92 . 60 scatter 0 . 47 0 . 42 0 . 46______________________________________ a two - stage search was employed to seek the minimum value of the cost function j . resolution of the first stage was specified interactively and was either 1 . 0 degree or 0 . 1 degree . this is the resolution referred to in the titles of the tables 1 through 4 , herein . the second stage resolution was 10 times the first stage resolution , i . e ., 0 . 1 degree and 0 . 01 degree , respectively . the range searched was ± 5 resolution elements about a nominal . to allow for possible entrapment at local minima , the euler angles associated with the 15 lowest values of j were stored and used as nominal values in the second stage . results of the second stage searches were sorted to extract the five results with the smallest value of j . the medians of the euler angles in this group were taken as the estimates . as a measure of confidence , the rms deviation of the individual estimates from the median was computed . table 1 lists these measurements as the &# 34 ; scatter .&# 34 ; table 1 summarizes the results obtained with a first stage resolution of 1 . 0 degree . note that good estimates of pitch are obtained throughout , while roll and yaw have slightly lower correlations with their true values . the scatter is a good qualitative indicator of accuracy . table 2 gives values of the cost function j for the true euler angles and those yielding the minimum value of j . for the lowest level of measurement noise , the true euler angles yielded a smaller value for j that the minimizing values . this indicates that the search did not uncover the neighborhood of the true minimum . to see if this could be corrected by increasing the resolution of the search , resolution was increased by 10 and the experiment was then repeated . the results are shown in table 3 . note that errors in roll and yaw are greatly diminished . additionally , the low scatter suggests reliable estimates . values of the cost function of this search are given in table 4 . as expected , values for the minimizing euler angle are smaller , especially the one for the lowest measurement noise . table 2______________________________________j ( 10 . sup . 6 ) with resolution = 1 . 0 deg . rms noise ( pixels ) euler angles 0 . 2680 0 . 6867 1 . 7340______________________________________true 1 . 4363 4 . 2261 7 . 3928minimizing 2 . 3221 3 . 6228 4 . 5824______________________________________ table 3______________________________________results ( deg ) with resolution = 1 . 0 deg . rms noise ( pixels ) euler angle 0 . 2680 0 . 6867 1 . 7340______________________________________roll true 5 . 00 5 . 00 5 . 00 estimated 4 . 96 4 . 86 5 . 25 scatter 0 . 04 0 . 06 0 . 05pitch true 70 . 10 70 . 10 70 . 10 estimated 70 . 10 70 . 10 70 . 10 scatter 0 . 00 0 . 00 0 . 00yaw true 93 . 80 93 . 80 93 . 80 estimated 93 . 75 93 . 75 93 . 75 scatter 0 . 01 0 . 03 0 . 01______________________________________ table 4______________________________________j ( 10 . sup . 6 ) with resolution = 1 . 0 deg . rms noise ( pixels ) euler angles 0 . 2680 0 . 6867 1 . 7340______________________________________true 1 . 4363 4 . 2261 7 . 3928minimizing 1 . 3736 3 . 4680 4 . 3660______________________________________ fig6 shows a block diagram of a tracking sensor 30 in accordance with the present invention . the tracking sensor 30 includes a memory 32 and a processor 34 . in accordance with the present invention , the tracking sensor 30 is operable to perform the operations shown in the flow diagram in fig7 . as shown in fig7 those operations includes storing a reference image at block 40 , sensing a second image at block 42 , mapping the reference image feature locations and the sensed image feature locations at block 44 , determining the relative position at block 46 , correlating the sensed image feature locations with the reference image feature locations at block 48 , and then calculating the tracking sensor 30 attitude using the correlations at block 50 . in summary , the preferred embodiment illustrates the use of image manipulation to improve the sensor pose determinations by more accurately producing sensor attitude measurements . by providing an exact , discrete time formulation of the attitude determination problem , the preferred embodiment enhances methods of passive ranging . although the present invention and its advantages have been described in detail , it should be understood that various changes , substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims .