Patent Publication Number: US-11024001-B2

Title: System and methods for attaining optimal precision stereoscopic direction and ranging through air and across refractive boundaries using minimum variance sub-pixel registration

Description:
BACKGROUND OF THE INVENTION 
     The present invention relates to the precise estimation of direction and range to remote objects submersed in a liquid, such as water, or a plasma, as well as remote, atmospheric objects residing in air, by means of a stereoscopic camera system. 
     Miniaturization of stereo-camera systems, for example in the field of smart-phones and robotics, has created a trend of progressively smaller camera separation distances requiring the ability to resolve ever minute parallax angles. Moreover, the ability to accurately estimate far-field object distances from digital images is sure to have a wide reaching impact on both real-time and offline applications for law enforcement and criminal forensics, security and surveillance, surveying, aerial drones and machine vision in general, but it requires accurate estimates of slim parallax angles. And notwithstanding the availability of high-resolution cameras, in order to accurately measure ever smaller parallax angles, modern stereo-camera systems must have the capability to perform sub-pixel registration, as is known in the art. 
     Sub-pixel registration, which estimates the pixel offset between images is also required when estimating far-field object distances, where the parallax angle narrows substantially and covers only a small number of pixels with the sub-pixel component comprising a non-negligible percentage of the offset between stereo images. 
     The problem of estimating the parallax angle from two stereo digital images is familiar to one trained in the art. The prior art discloses the use of peak detection on the spatial correlation between stereo images. Furthermore, the prior art discloses that more accurate peaks can be obtainable by the use of phase correlation, for which the 2-D Discrete Fourier Transform (DFT) is applied to both images, followed by the Schur (i.e. Hadamard) Product of the two transforms (one conjugated), and the result is magnitude normalized (to isolate the phase) and inverse transformed. Additionally, the prior art discloses a variant of phase correlation based on the 1-D Discrete Fourier Transform of horizontal image rows with averaging across the vertical dimension, which is substantially more computationally efficient than the 2-D approach and is the technique employed by the present invention. To address the problem of estimating the sub-pixel registration, the prior art discloses various techniques to compute an interpolated correlation peak based on an existing integer pixel offset already produced by the spatial or phase correlation approaches. Techniques disclosed in the prior art include Lagrange polynomial interpolation, linear and non-linear interpolation or regression, and polyphase interpolation, typically applied in the spatial domain to the correlation peak value, its nearest neighbors, and the associated spatial indices. Techniques such as these are heuristic in nature and make assumptions about the behavior of the phase correlation function between pixels. For example, the Lagrange polynomial interpolation technique assumes the phase correlation function is inherently polynomial, while the linear interpolation technique assumes a linear variation of the correlation function. 
     One shortcoming of such interpolation techniques is that they produce estimates based on a single measurement of the non-integer pixel offset. Thus, such techniques do not produce a consistent estimate based on an arithmetic mean of multiple measurements, wherein the estimate&#39;s error variance could decay as the number of measurements is increased, thus producing unbiased, high-precision estimates. Accordingly, such prior art interpolation techniques produce estimates that are often unpredictably biased relative to the true non-integer pixel offset and lack precision. 
     Another shortcoming of the prior art is lack of an optimal strategy for phase noise suppression. Phase noise limits the precision achievable with phase correlation approaches. Some prior art approaches to suppressing phase noise are based on truncating upper frequencies of the cross-power spectrum, from which phase correlation derives. Other approaches provide spectrum weighting of the normalized cross-power spectrum. Such approaches are heuristic and tend to distort the cross-power spectrum phase response in a non-predictable manner, which contributes to the bias in their associated sub-pixel registration estimates. 
     A final shortcoming of prior art approaches is the lack of quantifying and reporting the estimation error variance. For example, in the case of the polyphase interpolation approach, the reliability of the sub-pixel registration estimate is given as the percentage of power in the main peak of the phase correlation. While such a reliability value conveys the level of confidence in a given sub-pixel registration estimate, it does not indicate the precision of that estimate in units of pixels, as might be inferred from an estimation error variance. Characterization of error variance is a cornerstone of estimation theory and is essential to effective localization and decision-making. 
     BRIEF SUMMARY OF THE INVENTION 
     The goal of the present invention is to provide a system and methods capable of optimal precision direction and range determination through air and across a refractive boundary separating air from a liquid or plasma, using stereoscopic digital images captured by a stereo-camera system having dual left and right optical sensors, separated at a known distance, and operating in the visible or infrared spectrum. While achieving this stated goal, the present invention offers remedies to the shortcomings heretofore identified with respect to the prior art in the areas of precision, phase noise suppression, and error variance reporting. 
     The present invention achieves improved precision in sub-pixel registration by employing a globally minimum variance linear phase estimation method for determining precise estimates of the parallax angle. The minimum variance method disclosed by the present invention goes beyond conventional linear phase estimation approaches by selecting the optimal frequencies to use in the estimation process, in a manner that minimizes the affects of phase noise. 
     In addition, the present invention tracks measurement and estimation variances as they propagate through the system in order to provide a comprehensive error variance analysis of all estimated quantities, thereby informing the user of the precision for each estimate reported, wherein precision information is reported on the same display illustrating the subject under focus, such as a viewfinder or smart-phone display. 
     A first aspect of the invention uses a left and right stereoscopic camera system to make an optimal precision, sub-pixel resolution estimate of the parallax angle at the subject subtended by the dual optical sensor pair in order to triangulate an estimate of the line-of-sight (LOS) range to the subject under focus. To achieve optimal precision with robustness to noise, the invention first makes a consistent estimate of the spectral transfer function relating the transforms of the left and right images. It is assumed at the outset that the spectral transfer function is a frequency response having a phase response which is substantially linear with respect to frequency, with a slope equaling the sub-pixel registration, or fractional sample offset between left and right subject images. Next, the spectral transfer function estimate and its associated error variance are used to derive the signal-to-noise-ratio (SNR) for each frequency, and the SNR is used to predict the corresponding phase noise power levels, or phase error variances, which are then sorted in ascending order. Then, the present invention estimates the phase response of the spectral transfer function estimate, where the frequencies included in the estimate are selected efficiently and optimally to minimize the predicted error variance of the estimated slope of a linear phase dependence.  FIG. 1  illustrates the predicted minimum variance  101  achieved by the present invention.  FIG. 2  illustrates the high degree of phase linearity over accepted optimal frequencies  200 , in contrast to rejected frequencies  201 . The estimated slope is taken as the sub-pixel registration estimate. Additionally, the invention provides a means of calculating and reporting the precision of the sub-pixel estimate and the line-of-sight range estimate on which it depends. A second aspect of the invention uses the line-of-sight range estimate to calculate the horizontal range and altitude of an atmospheric subject. A third aspect of the invention employs the line-of-sight range estimate to calculate the horizontal range, slant-range, depth, and depression angle of a subject submersed in a liquid or plasma by compensating for refraction according to Snell&#39;s Law with respect to the liquid&#39;s or plasma&#39;s temperature and pressure; in both the second and third aspects, azimuth is calculated as an offset within the field-of-view relative to a known heading. Moreover, the second and third aspects provide means for deriving a precision for each estimated quantity based on the error variances of the line-of-sight range estimate and relevant measurements.  FIG. 3  provides an overview of the methods used by the stereo-camera system disclosed in the present invention. 
     In further regard to the first aspect of the invention, and in contrast to the prior art, the present invention discloses a consistent estimate of the true non-integer pixel offset and reports the error variance of the non-integer pixel offset estimate, indicating the precision of that estimate. A consistent estimate is defined as one whose expected value equals the random parameter being estimated, while having an error variance that tends to zero as the number of measurements approaches infinity. The present invention is capable of achieving such an estimate by exploiting two facets of dual camera systems. First, stereo-camera systems are aligned with the lateral baseline they are mounted on, which means the phase difference affects only the horizontal dimension, relative to the dual camera system&#39;s local reference frame. This fact suggests that the M vertically stacked rows (i.e. horizontal vectors) of pixels can provide multiple measurements of the cross power spectrum between left and right images that can be averaged to produce a frequency response having a reduced-noise phase response. Second, in a stereo-camera system, the left image is right-shifted relative to the right image, creating a linear phase response whose slope magnitude corresponds to the non-integer pixel offset separating the images. So, the reduced-noise phase response produced by averaging the cross power spectrum should be substantially linear, allowing for the regression of linear phase with respect to radian frequency. The slope of the phase response is determined by dividing the real-valued covariance, between the phase angle and radian frequency, by the mean-squared radian frequency in the sample. 
    
    
     
       BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING 
         FIG. 1 : Predicted Error Variance of Sub-Pixel Registration Estimate 
         FIG. 2 : Phase Response of Spectral Transfer Function 
         FIG. 3 : Overview 
         FIG. 4A : Integer Pixel Processing 
         FIG. 4B : Sub-pixel Processing 
         FIG. 4C : Alignment of Left/Right Overlapping Segments 
         FIG. 4D : Signal Model in the Complex Plane 
         FIG. 4E : Spectral Transfer Function Error Analysis in the Complex Plane 
         FIG. 4F : Integer Pixel Processing Equations 
         FIG. 4G : Sub-Pixel Processing Equations 
         FIG. 4H : SNR Estimate Equations 
         FIG. 4I : Phase Error Variance Estimate Equations 
         FIG. 4J : Sub-Pixel Error Minimization Equations 
         FIG. 4K : Sub-Pixel Estimate Equations 
         FIG. 4L : Sub-Pixel Error Analysis Equations 
         FIG. 4M : Cross-Correlator Computation Structure 
         FIG. 4N : Recurrence Relation Computation Structure 
         FIG. 5A : Optical Sensor 
         FIG. 5B : Parallax Angle Estimation Equations 
         FIG. 6A : Line-of-Sight Range Estimate Equations 
         FIG. 6B : Line-of-Sight Range Estimate Geometry 
         FIG. 7A : Estimate Equations (Atmospheric Subject) 
         FIG. 7B : Estimation and Measurement Error Variance Equations (Atmospheric Subject) 
         FIG. 7C : Altitude Estimate Error Variance Equations (Atmospheric Subject) 
         FIG. 7D : Range Estimate Error Variance Equations (Atmospheric Subject) 
         FIG. 7E : Problem Definition (Atmospheric Subject) 
         FIG. 8A : Problem Definition (Submersed Subject) 
         FIG. 8B : Application of Snell&#39;s Law (Submersed Subject) 
         FIG. 8C : Estimate Equations (Submersed Subject) 
         FIG. 8D : Estimation and Measurement Error Variance Equations (Submersed Subject) 
         FIG. 8E : Depth Estimate Error Variance Equations (Submersed Subject) 
         FIG. 8F : Range Estimate Error Variance Equations (Submersed Subject) 
         FIG. 8G : Slant-Range Estimate Error Variance Equations (Submersed Subject) 
         FIG. 8H : Depression Angle Estimate Error Variance Equations (Submersed Subject) 
         FIG. 8I : Estimation Error Variance Computation Structure (Submersed Subject) 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The system of the present invention comprises six methods as indicated in the overview shown in  FIG. 3 . 
     The purpose of the first method shown in Overview  FIG. 3 , Integer Pixel Registration  300 , is to determine the integer horizontal pixel offset of a subject captured in a left stereo image relative to a right stereo image. Integer Pixel Registration  300  informs method Sub-Pixel Registration  301  how to align subject pixels shared in common between left and right, and discard pixels that do not appear in both subjects, thereby reducing the effective noise level presented to Sub-Pixel Registration  301 . Method Integer Pixel Registration  300  comprises ten steps as illustrated in  FIG. 4A . 
     The first step of method Integer Pixel Registration  300  is Stereo Image Capture  400 , which simultaneously captures left and right images using a stereo-camera having corresponding left and right optical sensors separated at a known distance, and operating in the visible or infrared spectrum. 
     The second step of method Integer Pixel Registration  300  is Subject Selection Buffer  401 , which allows a user or automated agent to select from the right image, the smallest rectangular buffer that contains the subject, so that a substantial portion of the selection buffer reflects the line-of-sight range to the subject. The dimensions of the rectangle are defined to have M pixel rows and N pixel columns, comprising an M×N subject image. Similarly, the right subject selection buffer location and dimensions are then used to select the corresponding left subject. 
     The third step of method Integer Pixel Registration  300  is RGB to Luminance Conversion  402 , which converts Red-Green-Blue color triplet values to scalar luminance values for both left and right subject selection buffers. 
     The fourth step of method Integer Pixel Registration  300  is Subject Storage Buffer  403 , which comprises storing the M×N subject image in a storage means for later use by subsequent steps in method Integer Pixel Registration  300  as well as subsequent steps in method Sub-Pixel Registration  301 . 
     The fifth step of method Integer Pixel Registration  300  is Multiply M Rows by N-pt Hanning Window  404 . In this step, each row of both left and right subject images are multiplied by an N sample Hanning window for the purpose of reducing spectral smearing in the following step N-pt DFT of M Rows  405 . 
     The sixth step of method Integer Pixel Registration  300  is N-pt DFT of M Rows  405 , which comprises computing the one-dimensional N sample Discrete Fourier Transform (DFT) for each row of the left and right M×N subject images, wherein rows are treated as observations indexed by the lower-case letter i, as shown in  FIG. 4F  Equation  440 . The DFT is computed at N equally spaced frequencies which are indexed by the lower-case letter k, as shown in  FIG. 4F  Equation  439 . 
     The seventh step of method Integer Pixel Registration  300  is Estimate Cross-Power Spectrum (X)  406 , which comprises estimating the cross-power spectrum between the left and right subject images in order to determine the relative phase-shift between left and right subject images as a function of frequency. Due to noise, the true relative phase-shift between left and right subject images cannot be determined exactly and must be estimated. 
     To illustrate,  FIG. 4D  shows the complex plane with real-axis is  426 , imaginary-axis  427 , and a number of signal vectors showing the relationship between DFT samples for left and right subject images for a typical frequency k. Vector R_i(k)  428  is the observed right subject DFT sample for the i-th of M rows. Vector  429  is the theoretical left subject DFT sample for the i-th row, which is related to R_i(k)  428  by rotation over an unknown, negative phase angle  430  (measured counter-clockwise), reflecting the fact that the left subject image is right-shifted relative to the right subject image. However, the theoretical left subject DFT sample  429  is not observable due to the addition of complex noise vector z_i(k)  431 . In the presence of complex noise vector z_i(k)  431 , L_i(k)  432  is observed instead.  FIG. 4F  Equation  441  shows that L_i(k)  432  equals the sum of complex noise vector z_i(k)  431  and R_i(k)  428  phase-shifted by true phase angle  430 , which is linearly proportional to the pixel offset P, between left and right subject images, and wherein it is assumed that P is a natural number defined on the interval [0, N) for the purposes of method Integer Pixel Registration  300 , as shown in  FIG. 4F  Equation  444 . The result is a phase error  433  between the observed phase-shift and the unknown true phase-shift  430 . 
     It is assumed that R_i(k)  428  and z_i(k)  431  are statistically independent as shown in  FIG. 4F  Equation  442 . Furthermore, it is assumed that z_i(k)  431  from different rows are independent and identically distributed according to a complex normal distribution as shown in  FIG. 4F  Equation  443 , wherein the complex mean is zero and the variance, sigma squared subscript z, is split evenly between real and imaginary components. Under the stated statistical assumptions, the present invention seeks to estimate the cross-power spectrum between L_i(k) and R_i(k). 
     The cross-power spectrum X(k) is defined as the expected value of the product L_i(k) times the complex conjugate of R_i(k), as shown in  FIG. 4F  Equation  445 , where complex conjugation is denoted by asterisk. Multiplying  FIG. 4F  Equation  441  by the complex conjugate of R_i(k) amounts to the squared magnitude of R_i(k) times the phase-shift plus complex noise. By approximating the expectation as the sample mean taken over M rows, as indicated by sigma hat subscript LR in  FIG. 4F  Equations  446  and  447 , the present invention can produce a consistent estimate of the cross-power spectrum with an error variance that decays inversely as the number of row observations M increases. It should be noted that the resulting cross-power spectrum estimate represents a reduction from M dimensions to one dimension. 
     To facilitate sample mean estimation of the expectation of the product of a first variable, such as L_i(k), and the complex-conjugate of a second variable, such as R_i(k), the present invention discloses a first computation structure shown in  FIG. 4M . Each variable is assumed to have an observation sequence of length L, which are ordered such that opposing observations of first  49901  and second  49902  variables are in direct correspondence with respect to a given temporal, frequency or spatial index, such as the row index i. This first computation structure can also be used to calculate mean square value, wherein the first and second variable are duplicate copies for all observations, and variance, with the additional assumption of zero-mean variables. Also, the first computation structure can also operate on real-valued variables, wherein complex-conjugation of the second variable does not apply. Furthermore, multiply  49903  and accumulate  49904  operations can execute in parallel for limited segments of the input data. The final element of the first computation structure performs sample mean normalization by the observation length (L) of the observation sequences. As such, the present invention employs the first computation structure to implement  4 F Equation  447 , wherein the observation length (L) is equal to M, the number of horizontal rows in the left and right subject images. 
     The eighth and ninth steps of method Integer Pixel Registration  300  are Normalize Cross-Power Spectrum  407 , followed by N-pt IDFT to get Phase Corr, comprising dividing the cross-power spectrum estimate by its magnitude across all frequencies k to produce an all-pass response that approximates the phase-shift factor, and then taking the one-dimensional N sample Inverse Discrete Fourier Transform (IDFT) to produce a spatially shifted phase correlation sequence, as shown in  FIG. 4F  Equation  448 . 
     The tenth step of method Integer Pixel Registration  300  is 3-pt Sliding Window Peak Detector (Integer Pixel)  409 , comprising the search for the pixel offset that yields the peak sum of any three consecutive samples having adjacent circular (modulo N) indices, as shown in  FIG. 4F  Equation  449 . The resulting pixel offset reflects the group delay of the normalized cross-power spectrum and is taken as the integer pixel registration relating the left and right subject images, hereafter denoted P subscript int, as indicated in  FIG. 4F  Equation  449 . 
     The purpose of the second method shown in Overview  FIG. 3 , Sub-Pixel Registration  301 , is to estimate the sub-pixel registration utilizing the integer pixel registration estimated in method  300 , by step 3-pt Sliding Window Peak Detector (Integer Pixel)  409 , along with the left and right subject images stored by step Subject Storage Buffer  403 . Method Sub-Pixel Registration  301  comprises eleven steps. 
     The first step of method Sub-Pixel Registration  301 , as shown in  FIG. 4B , is Align L/R Overlapping Segments (Mx{N−P})  410 . The purpose of this first step is to discard subject image samples that are not shared, and do not correlate, between left and right images due to the integer pixel offset (P) separating them, and which contribute to noise that lowers the precision of the sub-pixel registration estimate. 
     The first step  410  begins by receiving the M×N left and right subject images from step  403  and the integer pixel registration from step  409 .  FIG. 4C  shows the procedure of aligning left and right subject images. All M rows  423  are retained for both left and right subject images. The column extent for the right subject image consists of the number of columns ({N−P})  421  that overlap with the left subject image plus the number (P) of non-overlapping columns  422 . Similarly, the column extent for the left subject image consists of the number ({N−P}) of columns  425  that overlap with the right subject image plus the number (P) of non-overlapping columns  424 . 
     Note that the overlapping portion of the left subject image is right-shifted P pixels relative to the left edge of the left subject image, while the overlapping portion of the right subject image is aligned with the left edge of the right subject image. 
     The first step  410  of method Sub-Pixel Registration  301  is completed by discarding P non-overlapping left and right pixel columns,  424  and  422  respectively, and retaining {N−P} left and right overlapping pixel columns,  425  and  421  respectively. 
     The second step of method Sub-Pixel Registration  301  is Multiply M Rows by {N−P} pt Hanning Window  411 . In this step, each row of both left and right subject images are multiplied by a {N−P} sample Hanning window for the purpose of reducing spectral smearing in the forthcoming DFT. 
     The third step of method Sub-Pixel Registration  301  is Zero-Pad P Samples for M Rows (M×N)  412 . In this step, P pixel columns to the right of each subject image are padded with zeros to restore the full count of N columns in left and right subject images. 
     The fourth step of method Sub-Pixel Registration  301  is N-pt DFT of M Rows  413 , which comprises computing the one-dimensional N sample Discrete Fourier Transform (DFT) for each row of the left and right M×N subject images, wherein rows are treated as observations indexed by the lower-case letter i, as shown in  FIG. 4G  Equation  440 . The DFT is computed at N equally spaced frequencies which are indexed by the lower-case letter k. However, because left and right subject images are real-valued, only non-negative frequency DFT samples are used in subsequent processing since negative-frequency DFT samples are redundantly conjugate symmetric, with the further exclusion of zero and half the sampling frequency, whose DFT samples are always real-valued and devoid of any phase information. The number of frequencies whose DFT samples are retained for subsequent processing is N/2−1, as shown in  FIG. 4G  Equation  450 . 
     The fifth step of method Sub-Pixel Registration  301  is Estimate L/R Spectral Transfer Function (H)  414 , which comprises estimating the spectral transfer function relating the left and right subject images in order to determine the relative phase-shift between left and right subject images as a function of frequency. Due to noise, the true sub-pixel relative phase-shift between left and right subject images cannot be determined exactly and must be estimated. According to the fifth step, the present invention makes an estimate of the spectral transfer function. 
     To illustrate the basic elements defining the spectral transfer function,  FIG. 4D  shows the complex plane with real-axis is  426 , imaginary-axis  427 , and a number of signal vectors that show the relationship between DFT samples for left and right subject images for a typical frequency k. Vector R_i(k)  428  is the observed right subject DFT sample for the i-th of M rows. Vector  429  is the theoretical left subject DFT sample for the i-th row, which is related to R_i(k)  428  by rotation over an unknown, negative phase angle  430  (measured counter-clockwise), reflecting the fact that the left subject image is right-shifted relative to the right subject image. However, the theoretical left subject DFT sample  429  is not observable due to the addition of complex noise vector z_i(k)  431 . In the presence of complex noise vector z_i(k)  431 , L_i(k)  432  is observed instead.  FIG. 4G  Equation  441  shows that L_i(k) 432 equals the sum of complex noise vector z_i(k)  431  and R_i(k)  428  phase-shifted by true phase angle  430 , which is linearly proportional to the pixel offset P, between left and right subject images, and wherein it is assumed that P is a real number defined on the interval (−1, 1) for the purposes of method Sub-Pixel Registration  301 , as shown in  FIG. 4G  Equation  451 . The result is a phase error  433  between the observed phase-shift and the unknown true phase-shift  430 . 
     It is assumed that R_i(k)  428  and z_i(k)  431  are statistically independent as shown in  FIG. 4G  Equation  442 . Furthermore, it is assumed that z_i(k)  431  from different rows are independent and identically distributed according to a complex normal distribution as shown in  FIG. 4G  Equation  443 , wherein the complex mean is zero and the variance, sigma squared subscript z, is split evenly between real and imaginary components. Under the stated statistical assumptions, the present invention seeks to estimate the spectral transfer function between left and right subject images as the slope of the estimate of L_i(k) with respect to R_i(k). 
     The true spectral transfer function takes the form of a unity-magnitude complex number whose phase angle is proportional to P, the sub-pixel registration offset, as shown in  FIG. 4G  Equation  452 . To illustrate the problem of estimating the spectral transfer function,  FIG. 4E  shows the complex plane with real-axis is  426 , imaginary-axis  427 , and the relationship between the true spectral transfer function and its estimate for a typical frequency k in the i-th of M rows. The spectral transfer function is the true, unknown unity magnitude phase-shift vector  434  that multiplies R_i(k) in order to produce L_i(k) in the absence of noise. The true associated phase-shift angle  430  is unknown, but drawn as a negative angle (measured counter-clockwise) to reflect that the left subject image is right-shifted relative to the right subject image. Arc  438  is the unit circle that circumscribes the true phase-shift vector  434 . Due to the addition of complex noise vector  435 , the observed phase-shift vector  436  has a phase-shift angle error  437 . By modeling the spectral transfer function as the slope of the estimate of L_i(k) with respect to R_i(k), the present invention seeks to estimate the phase-shift angle  430  as well as the error variance of that estimate given the statistics of phase-shift angle error  437 , the later of which is derived from the signal-to-noise ratio of the spectral transfer function estimate. 
     To complete the fifth step of method Sub-Pixel Registration  301 , Estimate L/R Spectral Transfer Function (H)  414 , the slope of the L_i(k) estimate with respect to R_i(k) is obtained by first noting that multiplying  FIG. 4G  Equation  441  by the complex conjugate of R_i(k), or R_i*(k), and taking the expectation of both sides yields  FIG. 4G  Equation  453 . Then, solving for spectral transfer function H(k) produces the cross-power spectrum E[L_i(k).R_i*(k)] divided by the variance E[R_i(k).R_i*(k)], where the dividend and divisor expectations have consistent estimates in their sample means as shown in  FIG. 4G  Equations  454 ,  447 , and  455 . The first computation structure disclosed in the present invention implements  FIG. 4G  Equations  447  and  455 , wherein the observation length (L) is equal to M, the number of horizontal rows in the left and right subject images. 
     The sixth step of method Sub-Pixel Registration is Estimate SNRs  415 . As mentioned above, the signal-to-noise-ratio (SNR) is used to derive the error variance of the phase-shift angle estimate.  FIG. 4H  Equation  456  defines for the i-th of M rows, the estimate L_i(k) hat as the product R_i(k) times the spectral transfer function estimate H(k) hat.  FIG. 4H  Equation  457  defines for the i-th of M rows, the complex error, z_i(k) hat, between the observed and estimated left subject DFT samples, L_i(k) and L_i(k) hat respectively.  FIG. 4H  Equation  458  defines the error variance, sigma squared subscript z hat, of z_i(k), and is implemented by the first computation structure disclosed in the present invention, wherein the observation length (L) is equal to M.  FIG. 4H  Equation  459  estimates the H(k) hat error variance by summing over each of M random inputs to H(k) hat, the product of the error variance of that input times a scale factor equaling the magnitude squared of the partial derivative of H(k) hat with respect to that input. In light of  FIG. 4G  Equation  441 , the M random inputs are simply z_i(k) with error variance sigma squared subscript z hat, and the partial derivative with respect to z_i(k) is equivalent to the partial derivative with respect to L_i(k) using the chain rule, yielding  FIG. 4H  Equation  460 .  FIG. 4H  Equation  460  demonstrates that the error variance of H(k) hat decays inversely as the number of row observations M increases, making H(k) hat a consistent estimate of the spectral transfer function. Finally, the SNR is the magnitude squared of the estimate H(k) hat divided by the H(k) hat error variance, as shown in  FIG. 4H  Equation  461 . 
     The seventh step of method Sub-Pixel Registration  301  is Estimate Phase Err Variances  416 . Recall from  FIG. 4E  that H(k) hat will suffer phase angle error  437  due to the presence of complex noise vector  435 . The purpose of the seventh step is to estimate the variance of phase angle error  437  using SNR calculated in the sixth step. 
     Given the estimate H(k) hat and the error variance estimated in step six, the true spectral transfer function H(k) follows a complex normal distribution, in rectangular coordinates (x, y), with complex mean H(k) hat and bi-variate variances equal to half the error variance of H(k) hat, as shown in  FIG. 4I  Equations  462 ,  463 , and  464 . The marginal density of the phase error, theta, of H(k) with respect to H(k) hat can be found by making a change of variables to polar coordinates (radius, theta) and integrating over radius from  0  to infinity, yielding  FIG. 4I  Equation  465 , wherein capital phi denotes the Gauss Error function. 
     Under simplifying assumptions of high SNR, and small phase error theta,  FIG. 4I  Equation  465  approaches a normal distribution as shown in  FIG. 4I  Equation  472 . Furthermore, under these assumptions, and according to definitions  FIG. 4I  Equation  462  and  463 , the phase error variance, sigma squared subscript theta hat, for a frequency k can be estimated as one over twice SNR for that frequency, as shown in  FIG. 4I  Equation  473 . 
     The eighth and ninth steps of method Sub-Pixel Registration  301  select the optimal frequencies to use when estimating the sub-pixel registration as the slope of the spectral transfer function phase response. As shown above, there is an inverse relationship between SNR and phase error variance. So, frequencies with the lowest spectral transfer function SNR should be omitted from the sub-pixel registration estimation, since their respective phase estimates have unreliably high error variances. 
     However, if too many frequencies are omitted, then the estimate will lack sufficient data to effectively cancel estimation noise. So, there is clearly an optimal set of frequencies that minimizes the error variance of the sub-pixel registration estimate. 
     To select the optimal set of frequencies, the present invention first sorts the frequencies by ascending phase error variance, and then from all Q length sets of sorted frequencies spanning the lowest error variance up to the Q-th lowest error variance, selects the Q that minimizes the predicted sub-pixel registration error variance. Note that the Q length sets will in general not have contiguous frequencies. 
     The eighth step of method Sub-Pixel Registration  301  is Sort Freqs by Ascending Phase Err Variance  417 . The purpose of the eighth step is to sort frequency indices (k) such that corresponding phase error variances, defined by  FIG. 4I  Equation  473 , increase as a reference index (q) is incremented, and a mapping function (s) from q to k is provided. The sorting and index mapping functions are defined in  FIG. 4J  Equation  477 , wherein the frequency (k) with the highest phase error variance is mapped to q=N/2−1, the total number of frequencies. 
     The ninth step of method Sub-Pixel Registration  301  is Minimize Predicted Sub-Pixel Err Variance  418 . The purpose of the ninth step is to select the optimal set of frequencies k that minimize the predicted sub-pixel registration error variance. 
     Recall that the spectral transfer function H(k) takes the form of a unity-magnitude complex number whose phase angle is proportional to P, the sub-pixel registration offset, as shown in  FIG. 4J  Equation  452 . Then, the observed phase angle theta(k) is defined by  FIG. 4J  Equation  474 , wherein it is assumed that the epsilon(k) are independent and identically distributed normal random variables that are statistically independent from omega(k), the radian frequency, as shown in  FIG. 4J  Equations  475  and  476 . 
     The sub-pixel registration estimate (P hat) is given by  FIG. 4J  Equation  479 ,  480 , and  481  and determines how the associated error variance is calculated.  FIG. 4J  Equation  482  estimates the predicted sub-pixel registration error variance, or sigma squared subscript P tilde, by summing over each of Q inputs to P hat, the product of the error variance of that input times a scale factor equaling the magnitude squared of the partial derivative of P hat with respect to that input, which is given by  FIG. 4J  Equation  483 . Note that the predicted sub-pixel registration error variance is a function of Q and can therefore be minimized with respect to Q. 
     Because predicted sub-pixel registration error variance of consecutive values of Q share common elements in their numerator and denominator, it is possible to construct a recurrence relation for efficient computation of the error variance, as shown in  FIG. 4J  Equations  484 ,  485 , and  486 . The recurrence relation is the ratio of partial sum A(Q) divided by the square of partial sum B(Q). The present invention discloses a second computation structure, shown in  FIG. 4N , that implements this recurrence relation. The second computation structure comprises a first computation unit, wherein partial sum A(Q−1)  49908 , provided by delay element  49907 , and its increment, provided by  49906 , are added together by summing element  49909  to produce partial sum A(Q). Similarly, the second computation structure comprises a second computation unit, wherein partial sum B(Q−1)  49912 , provided by delay element  49911 , and its increment, provided by  49910 , are added together by summing element  49913  to produce partial sum B(Q). 
     The ninth step of method Sub-Pixel Registration  301  is completed by evaluating the recurrence relation for Q between 1 and N/2−1 and determining the optimal value Q subscript opt that minimizes the predicted sub-pixel registration error variance, as shown in  FIG. 4J  Equation  487 . 
       FIG. 1 . shows the sorted phase error variance estimates 100 plotted on the right axis as a function of sample size Q.  FIG. 1  also shows the predicted sub-pixel registration error variance plotted on the left axis as a function of sample size Q. The predicted minimum error variance  101  is achieved at optimal sample size Q opt  102 . For the example shown, the predicted minimum error variance  101  is roughly 6E-5, the predicted root-mean-square error is roughly 0.008 pixels, and Q opt is 52. Note that where a predicted error variance at sample size Q is an outlier with respect to the other Q−1 variances in the sample, the Q-th and higher error variances are rejected, which has the effect of moderating the spread between small and large phase error variances, thus lowering the kurtosis of the sample, and preventing heavy-tailed phase errors. 
     With regard to complexity, letting n represent the N/2−1 frequencies under consideration, sorting requires O(n log n) iterations, while the recurrence relation requires O(n) iterations, making the present invention an efficient way to select the frequencies that achieve the minimum predicted sub-pixel registration error variance. By comparison, to exhaustively search for the lowest predicted sub-pixel registration error variance over every possible sample of Q frequencies, for all possible sample sizes Q, has complexity exceeding n!/Q!/(n−Q)! iterations, wherein the worst case is Q=n/2. 
     The tenth step of method Sub-Pixel Registration  301  is Estimate Sub-Pixel Registration  419 . The purpose of the tenth step is to produce the sub-pixel registration estimate. First, the phase angles of the spectral transfer function estimate are calculated at each of the Q opt frequencies s subscript q, as shown in  FIG. 4K  Equation  478 . Note that while the arc-tangent function is expensive computationally, it is only executed on Q opt samples. Next, the phase angles, theta of s subscript q, are used to estimate the sub-pixel registration estimate (P hat), as shown in  FIG. 4K  Equations  479 ,  488 , and  489 . The first computation structure disclosed in the present invention implements  FIG. 4K  Equation  488 , with first input variable being the real-valued phase angle and the second input variable being the corresponding real-valued radian frequency with an observation length (L) equal to Q opt. Similarly, the first computation structure disclosed in the present invention implements  FIG. 4K  Equation  489 , the mean square radian frequency, with first and second input variables being the real-valued radian frequency with an observation length (L) equal to Q opt. 
     The eleventh step of method Sub-Pixel Registration  301  is Estimate Actual Sub-Pixel Err Variance  420 . The purpose of the eleventh step is to estimate the error variance of the sub-pixel registration estimate. Using the sub-pixel registration estimate, P hat, the estimate of the phases, theta hat, is estimated as shown in  FIG. 4L  Equation  490 . Next, the phase estimation error, epsilon hat of s subscript q, is found by subtracting the phase estimates from the observed phases, as shown in  FIG. 4L  Equation  491 . Then, the error variance of the phase estimates is simply the sample variance of the phase errors, as shown in  FIG. 4L  Equation  492 ; this sample variance is computed by the first computation structure disclosed in the present invention, wherein the first and second input variables are the real-valued phase errors with an observation length (L) equal to Q opt. 
     Finally,  FIG. 4L  Equation  493  estimates the actual sub-pixel registration error variance, or sigma hat squared subscript P, by summing over each of Q opt inputs to P hat, the product of the phase error variance times the magnitude squared of a scale factor equaling the partial derivative of P hat with respect to the phase error. Note that in light of  FIG. 4J  Equation  474 , the partial derivative of P hat with respect to the phase error (epsilon) is equivalent to the partial derivative with respect to the phase (theta), using the chain rule. Substituting  FIG. 4K  Equations  479 ,  488 , and  489  into  FIG. 4L  Equation  493  yields  FIG. 4L  Equation  494 . 
     For subsequent methods, the sub-pixel registration estimate P hat will be referred to as P subscript sub, as shown in  FIG. 4L  Equation  495 . 
       FIG. 2  illustrates the result of employing the present invention to estimate the sub-pixel registration by plotting the spectral transfer function phase as a function of radian frequency. The Q opt frequencies  200  accepted for sub-pixel registration estimation are shown as triangles, while the rejected frequencies  201  are shown as plus symbols. It is clear that the present invention is capable of isolating the phase samples  200  with low error variance from phase samples  201  with high error variance such that the minimum sub-pixel registration error variance is achieved. Note that accepted frequencies show a strong linear phase relationship with a slope less than unity, indicating the sub-pixel registration. For the example shown, the actual error variance is roughly 4E-5, while the actual root-mean-square error is roughly 0.006 pixels, showing a comparable order of magnitude to the predicted quantities shown in  FIG. 1 . 
     Sub-pixel registration precision is driven by the error variance, which improves with increases in subject image dimensions M and N as well as SNR. Recall that the predicted phase error variances given by  FIG. 4I  Equation  473  are proportional to the estimated spectral transfer function error variances given by  FIG. 4H  Equation  460 , which decay inversely as the number of rows M increases. Moreover, it is assumed that predicted phase error variances  473 , are positively correlated with the actual phase error variances given by  FIG. 4L  Equation  492 . Hence, the sub-pixel registration error variance, given by  FIG. 4L  Equation  494 , decays inversely as the number of rows M increases. Similarly, predicted phase error variance  473 , and hence sub-pixel registration error variance, decays inversely as SNR increases, as shown in  FIG. 4I  Equation  473 . Furthermore,  FIG. 4L  Equation  494  shows that the sub-pixel registration error variance decays inversely as Q opt increases. Q opt can be increased by increasing the number of columns N, which increases the frequency resolution and, hence, the number of frequencies with a low predicted phase error variance. 
     The purpose of the third method shown in Overview  FIG. 3 , Parallax Angle Estimation  302 , is to estimate the parallax angle between left and right subject images given the sub-pixel registration estimate already estimated, in addition to standard optical sensor parameters. Additionally, the third method estimates the error variance of the parallax angle estimate.  FIG. 5A  illustrates the optical sensor along with standard optical sensor parameters comprising diagonal dimension  500 , horizontal pixel resolution  502 , and vertical pixel resolution  503 . 
     The first step of method Parallax Angle Estimation  302  is to derive the horizontal dimension D subscript x  501  of the optical sensor from the standard optical sensor parameters diagonal dimension D  500 , horizontal pixel resolution rho subscript x  502 , and vertical pixel resolution rho subscript y  503 , as shown in  FIG. 5B  Equation  504 . 
     The second step of method Parallax Angle Estimation  302  is to derive the horizontal field-of-view FOV subscript x from the horizontal dimension D subscript x  501  of the optical sensor and the standard optical parameter focal length, as shown in  FIG. 5B  Equation  505 . 
     The third step of method Parallax Angle Estimation  302  is to derive angular resolution per pixel rho subscript psi from the horizontal field-of-view FOV subscript x and the horizontal pixel resolution rho subscript x  502 , as shown in  FIG. 5B  Equation  506 . 
     The fourth step of method Parallax Angle Estimation  302  is to derive the parallax angle between left and right subject images from the angular resolution per pixel rho subscript psi  506  and the sum of integer pixel P subscript int and sub-pixel P subscript sub registrations, as shown in  FIG. 5B  Equation  507 . Recall that integer pixel registration was defined in  FIG. 4F  Equation  449 , and sub-pixel registration was defined in  FIG. 4L  Equation  495 . 
     The fourth step of method Parallax Angle Estimation  302  is completed by estimating the error variance of parallax angle psi as the product of sub-pixel registration error variance  FIG. 4L  Equation  494  times a scale factor equaling the magnitude squared of the partial derivative of parallax angle psi with respect to sub-pixel registration P subscript sub, as shown in  FIG. 5B  Equations  508  and  509 . 
     The fourth method shown in Overview  FIG. 3 , Line-of-Sight Range Estimation  303 , comprises estimating a line-of-sight range estimate based on the parallax angle estimate  FIG. 5B  Equation  507  and the dual optical sensor spacing (delta), as shown in  FIG. 6A  Equation  600 . Additionally, the fourth method comprises estimating the error variance of the line-of-sight range estimate as shown in  FIG. 6A  Equations  601  and  602 . 
       FIG. 6B  illustrates the line-of-sight range estimate geometry. The left optical sensor is located at position  603  and the right optical sensor is located at position  604 . The left and right optical sensors are separated by a distance (delta)  605 . The subject, whose line-of-sight range is to be estimated relative to the left optical sensor  603 , is located at position  606 . Angle  607  is the parallax angle (psi) subtended by the dual optical sensors relative to the subject position  606 . Direction angle  608  (beta) measures the subject  606  direction relative to an axis collinear with left and right optical sensor positions  603  and  604 , respectively, within a plane formed by left and right optical sensors  603  and  604  and subject  606 . Thus, the line-of-sight range (r subscript LOS)  609  relative to the left optical sensor  603 , is given by  FIG. 6A  Equation  600 . 
     The error variance of the line-of-sight range estimate is obtained by summing the product of error variances for parallax angle (psi) and subject direction angle (beta) times the corresponding squared partial derivatives of line-of-sight range with respect to parallax angle (psi) and subject direction angle (beta), as shown in  FIG. 6A  Equations  601  and  602 . 
     The fifth method shown in Overview  FIG. 3  is Atmospheric Subject: Altitude &amp; Range  304 . The purpose of the fifth method is to estimate the altitude and horizontal range of a subject residing in air based on known measurements of the dual optical sensor&#39;s height above ground (h), depression/elevation angle (phi) and the estimate of the line-of-sight range (r subscript LOS) obtained in the fourth method  303 , to compute error variances for all estimated quantities, and to report all information on a display for review by a user. The display can be a camera viewfinder, a smart-phone display, or a generic computer display used for illustrating a subject under focus. 
       FIG. 7E  illustrates the problem definition. A stereo-camera  712  with dual optical sensors is positioned at height (h)  713  above zero-altitude plane  714 . The stereo-camera  712  captures left and right images of subject  715 , which is at altitude (a)  716 . Relative to the stereo-camera, the subject is positioned at a horizontal range (r)  717 , line-of-sight range (r_LOS)  718 , and depression/elevation angle (phi)  719 . 
     The set of possible depression/elevation angles (phi) is defined by  FIG. 7A  Equation  700 . The altitude estimate is given by  FIG. 7A  Equation  701 . And the horizontal range is given by  FIG. 7A  Equation  702 . 
     Additionally, the fifth method estimates error variances for altitude and horizontal range estimates in terms of the error variances of the underlying measurements and prior estimates on which they are based. The error variance for line-of-sight range is given by  FIG. 7B  Equation  602 , the error variance for height,  FIG. 7B  Equation  703 , is determined by user behavior, and the error variance for depression/inclination angle,  FIG. 7B  Equation  704 , is determined by the error variance of a gyroscope output, sigma squared subscript G, that provides the depression/inclination angle measurement. 
     The altitude estimate is derived from three inputs: the line-of-sight range estimate (r subscript LOS) and the measurements height (h) and depression/elevation angle (phi). So, the error variance of the altitude estimate is obtained by the summing over all inputs, the error variance for each input times a scale factor equaling the squared partial derivative of altitude with respect to said input, as given by  FIG. 7C  Equations  705 ,  706 ,  707 , and  708 . 
     The horizontal range is derived from two inputs: line-of-sight range estimate and depression/elevation angle. So, the error variance of the horizontal range estimate is obtained by summing over all inputs, the error variance for each input times a scale factor equaling the squared partial derivative of horizontal range with respect to said input, as given in  FIG. 7D  Equations  709 ,  710 , and  711 . 
     The sixth method shown in Overview  FIG. 3  is Submersed Subject: Depth, Range &amp; D/E Angle  305 . The purpose of the fifth method is to estimate subject depth, horizontal range, slant-range, and depression angle for a subject submersed in a liquid such as water, or a plasma, to compute error variances for all estimated quantities, and to report all information on a display for review by a user. The display can be a camera viewfinder, a smart-phone display, or a generic computer display used for illustrating a subject under focus. 
       FIG. 8A  illustrates the problem definition for the case where the subject is submersed in liquid water. A stereo-camera  800  with dual optical sensors is positioned at height (h)  801  above water surface  802 . The stereo-camera  800  captures left and right images of subject  803 , which is submersed at a depth (d)  804 . Relative to the stereo-camera, the subject is positioned at a horizontal range (r)  805 , slant-range (r_s)  806 , and depression angle (phi)  807 . 
       FIG. 8B  illustrates the application of Snell&#39;s Law. The stereo-camera  800  is oriented along a transmission range (r_t)  808  passing through air and making transmission angle (gamma_t)  809  with respect to the vertical axis at the point where a virtual subject image ray leaves water surface  802 , wherein gamma_t  809  is between 0 and gamma_c  816 , the critical angle, as shown in  FIG. 8C  Equation  813 . Beyond the critical angle  816 , total internal reflection occurs at water surface  802 , and the subject is not visible to the stereo-camera system. The critical angle and the ratio of incidence (n_i) index of refraction (of water) to transmission (n_t) index of refraction (of air) are defined in  FIG. 8C  Equation  816 . Stereo-camera  800  is oriented facing down at the (pi/2) complementary angle of transmission angle (gamma_t)  809 , the apparent depression angle. The virtual subject image ray appears to travel toward the surface on virtual path  810 , but actually travels on refracted path  811 , wherein both paths have length equaling incidence range (r_i). Note that virtual path  810  and transmission range (r_t)  808  are collinear, as indicated by the 180 degree arc. The real path  811  makes angle (gamma_i)  812  with respect to the vertical axis at the point where refracted subject image ray leaves the water. 
     The first step of method Submersed Subject: Depth, Range &amp; D/E Angle  305  is to use line-of-sight estimate (r_LOS) and measurements height (h), transmission angle (gamma_t), indices of refraction for incidence (n_i) in water and for transmission (n_t) in air, to estimate the quantities transmission range (r_t)  808  with  FIG. 8C  Equation  814 , incidence range (r_i)  810  with  FIG. 8C  Equation  815 , incidence angle (gamma_i)  812  with  FIG. 8C  Equation  817 , depth (d)  804  with  FIG. 8C  Equation  818 , horizontal range (r)  805  with  FIG. 8C  Equation  819 , slant-range (r_s)  806  with  FIG. 8C  Equation  820 , and depression angle (phi)  807  with  FIG. 8C  Equation  821 . 
     The second step of method Submersed Subject: Depth, Range &amp; D/E Angle  305  is to use error variances for line-of-sight range (r_LOS) estimate  FIG. 8D  Equation  602 , and error variances for measurements height (h)  FIG. 8D  Equation  703 , transmission angle (gamma_t)  FIG. 8D  Equation  822 , and ratio of incidence to transmission index of refraction (n)  FIG. 8D  Equations  823 - 824 , to estimate error variances for estimates depth (d), horizontal range (r), slant-range (r_s) and depression angle (phi). Note that the error variance for height (h)  FIG. 8D  Equation  703  is determined by user behavior and the error variance for transmission angle (gamma_t)  FIG. 8D  Equation  822  is derived from a gyroscope error variance, sigma squared subscript G, associated with a gyroscope that provides a measurement of the apparent depression angle, which is the (pi/2) complementary angle of transmission angle (gamma_t). Also, note that the error variance for ratio of incidence to transmission indices of refraction (n)  823  is itself an estimate based on a temperature dependence and barometric pressure dependence along with associated error variances. So, the error variance of estimate (n)  823  is obtained by summing the error variance of a temperature (T) measurement times the squared partial derivative of n  823  with respect to temperature (T) plus the error variance of a barometric pressure (B) measurement times the squared partial derivative of n  823  with respect to barometric pressure (B), as shown in  FIG. 8D  Equation  824 . 
     The error variance of the depth (d) estimate is obtained by summing over all inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), the error variance for each input times the square of a scale factor equaling the partial derivative of depth (d) with respect to said input, as shown in  FIG. 8E  Equation  825 ; wherein it is assumed that each input is statistically independent from all other inputs. Using the chain rule, the partial derivatives of depth (d) with respect to inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), are given in  FIG. 8E  Equations  826 ,  827 ,  828 , and  829 , respectively. The present invention discloses a third computation structure, shown in  FIG. 8I , for efficiently obtaining the error variance of the depth (d) estimate by computing first, second, third and fourth scale factors  846  using computation unit  845  for inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively; wherein the squares of said first, second, third and fourth scale factors  846  are multiplied by error variances for line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively, and summed by depth error variance computation unit  847 . 
     The error variance of the horizontal range (r) estimate is obtained by summing over all inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), the error variance for each input times the square of a scale factor equaling the partial derivative of horizontal range (r) with respect to said input, as shown in  FIG. 8F  Equation  830 ; wherein it is assumed that each input is statistically independent from all other inputs. Using the chain rule, the partial derivatives of horizontal range (r) with respect to inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), are given in  FIG. 8F  Equations  831 ,  832 ,  833 , and  834 , respectively. The third computation structure, shown in  FIG. 8I , further comprises a computation unit  848  that computes fifth, sixth, seventh and eighth scale factors  849  for inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively; wherein said fifth, sixth, seventh and eighth scale factors  849  are derived from said first, second, third and fourth scale factors  846 , respectively; wherein the squares of said fifth, sixth, seventh and eighth scale factors  849  are multiplied by error variances for line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively, and summed by horizontal range error variance computation unit  850 . 
     The error variance of the slant-range (r_s) estimate is obtained by summing over all inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), the error variance for each input times the square of a scale factor equaling the partial derivative of slant-range (r_s) with respect to said input, as shown in  FIG. 8G  Equation  835 ; wherein it is assumed that each input is statistically independent from all other inputs. Using the chain rule, the partial derivatives of slant-range (r_s) with respect to inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), are given in  FIG. 8G  Equations  836 ,  837 ,  838 , and  839 , respectively. The third computation structure, shown in  FIG. 8I , further comprises a computation unit  851  that computes ninth, tenth, eleventh and twelfth scale factors  852  for inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively; wherein said ninth, tenth, eleventh and twelfth scale factors  852  are derived from said first, second, third and fourth scale factors  846 , respectively; wherein said ninth, tenth, eleventh and twelfth scale factors  852  are also derived from said fifth, sixth, seventh and eighth scale factors  849 , respectively; wherein the squares of said ninth, tenth, eleventh and twelfth scale factors  852  are multiplied by error variances for line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively, and summed by slant-range error variance computation unit  853 . 
     The error variance of the depression angle (phi) estimate is obtained by summing over all inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), the error variance for each input times the square of a scale factor equaling the partial derivative of depression angle (phi) with respect to said input, as shown in  FIG. 8H  Equation  840 ; wherein it is assumed that each input is statistically independent from all other inputs. Using the chain rule, the partial derivatives of depression angle (phi) with respect to inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), are given in  FIG. 8H  Equations  841 ,  842 ,  843 , and  844 , respectively. The third computation structure, shown in  FIG. 8I , further comprises a computation unit  854  that computes thirteenth, fourteenth, fifteenth and sixteenth scale factors  855  for inputs line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively; wherein said thirteenth, fourteenth, fifteenth and sixteenth scale factors  855  are derived from said first, second, third and fourth scale factors  846 , respectively; wherein said thirteenth, fourteenth, fifteenth and sixteenth scale factors  855  are also derived from said fifth, sixth, seventh and eighth scale factors  849 , respectively; wherein the squares of said thirteenth, fourteenth, fifteenth and sixteenth scale factors  855  are multiplied by error variances for line-of-sight range (r_LOS), height (h), transmission angle (gamma_t), and index of refraction ratio (n), respectively, and summed by depression angle error variance computation unit  856 . 
     Thus, the present invention discloses substantial improvements in stereo-camera direction and ranging systems and methods, in comparison to the prior art, by providing optimal sub-pixel registration precision using a minimum variance linear phase estimation method, along with comprehensive error analysis and reporting. And notwithstanding that the above proceedings describe, to one of ordinary skill, the present best mode of practicing the present invention, additional modes and intended uses may become apparent to those of ordinary skill, based on the subject matter presented in the above proceedings. As such, the present invention should not be limited by the preferred embodiment set forth above, but should include all embodiments falling within the limitations, in spirit and scope, of the following claims.