US9521326B2 - Misregistration correction using non-pinhole camera model and non-linear scan equations - Google Patents

Abstract

A method of misregistration correction in a line scanning imaging system includes: generating a model of scan motion over a focal plane of the imaging system, using a coupled system of non-linear scan equations with constant coefficients; estimating programmed motion positions across a plurality of detector junction overlap regions via a state transition matrix solution to the scan equations; at each detector junction overlap region, measuring actual motion positions via image correlation of overlapping detectors; generating differences between the actual motion positions and the estimated programmed motion positions; estimating updates to the constant coefficients based on the generated differences; generating corrections from the estimated updates to remove unwanted motion; and applying the updates to the constant coefficients.

Description

RELATED APPLICATION

This application claims the benefit of U.S. Provisional Patent Application No. 61/820,914, filed in the U.S. Patent and Trademark Office on May 8, 2013, the entire contents of which are incorporated herein by reference.

BACKGROUND

1. Technical Field

This disclosure relates to generation of images of remote objects and, more particularly, to generation of corrected images generated from a scanning sensor on a moving platform.

2. Discussion of Related Art

Using line scanning, an electro-optical sensor can acquire image data for a remote object while in motion. Support data associated with each acquisition by the sensor contains measurement histories for position, velocity, attitude angles and rates, and a priori error variances. Accurate support data is important to the image formation process, as errors can manifest themselves as geometric artifacts.

Misregistration correction (MRC) is a component in the image formation chain of a line scanning electro-optical sensor. MRC removes the geometric artifacts introduced by unprogrammed motion, which is unwanted motion not commanded as part of the image acquisition process. It can arise from vibrations, sensor attitude and position knowledge errors, etc.

SUMMARY

According to one aspect, a method of misregistration correction in a line scanning imaging system is provided. The method includes: generating a model of scan motion over a focal plane of the imaging system, using a coupled system of non-linear scan equations with constant coefficients; estimating programmed motion positions across a plurality of detector junction overlap regions via a state transition matrix solution to the scan equations; at each detector junction overlap region, measuring actual motion positions via image correlation of overlapping detectors; generating differences between the actual motion positions and the estimated programmed motion positions; estimating updates to the constant coefficients based on the generated differences; generating corrections from the estimated updates to remove unwanted motion; and applying the updates to the constant coefficients.

In some exemplary embodiments, estimating updates is performed using least-squares estimation. In some exemplary embodiments, the image correlation comprises normalized cross-correlation. In some exemplary embodiments, the image correlation comprises lag product cross-correlation. In some exemplary embodiments, the image correlation comprises least squares cross-correlation.

In some exemplary embodiments, the model of scan motion is generated over a predetermined time interval,

In some exemplary embodiments, the scan equations comprise a set of differential equations with constant coefficients. In some exemplary embodiments, the differential equations are first-order differential equations.

In some exemplary embodiments, the method further includes computing the scan equation coefficients using sensor platform parameters. In some exemplary embodiments, the sensor platform parameters comprise at least one of sensor position, velocity, attitude angles and rates. In some exemplary embodiments, the method further includes approximating errors in the scan equation coefficients from measurement errors in the platform parameters. In some exemplary embodiments, the method further includes generating a covariance matrix of focal plane coordinate errors. In some exemplary embodiments, the method further includes estimating platform parameter errors from the differences between the actual motion positions and the estimated programmed motion positions.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is further described in the detailed description which follows, in reference to the noted plurality of drawings by way of non-limiting examples of embodiments of the present disclosure, in which like reference numerals represent similar parts throughout the several views of the drawings.

FIG. 1 includes a schematic diagram of a focal plane scanning geometry/configuration.

FIG. 2 is a schematic functional flow diagram of a detector junction correlation process.

FIG. 3 includes a schematic functional block diagram of a detector junction measurement process.

FIG. 4 is a schematic diagram of traditional pinhole camera projection equations.

FIG. 5 includes a schematic logical flow diagram of generating correction terms to be added to the traditional pinhole camera model, according to some exemplary embodiments.

FIG. 6 includes a schematic diagram of an approach to working around a singularity issue, according to some exemplary embodiments.

FIG. 7 is a schematic diagram of an error perturbation model, according to some exemplary embodiments.

FIG. 8 includes a schematic block diagram of an MRC process, according to some exemplary embodiments.

FIG. 9 includes a schematic logical flow diagram of a line-by-line MRC computational process, according to some exemplary embodiments.

DETAILED DESCRIPTION

Some previous electro-optical line scanning sensors were modeled as pinhole cameras. Unknown amounts of system error are introduced using this approximation. As a result, new approaches have been developed for geometric corrections that do not adhere to the traditional model. An example is misregistration correction (MRC), which affects image quality. MRC removes geometric artifacts introduced by unprogrammed motion, which is unwanted motion not commended as part of the image acquisition process. It can arise from vibrations, sensor attitude and position knowledge errors, etc. An innovative new approach for performing MRC is described herein. The approach described herein does not assume a pinhole camera model. It simultaneously estimates and corrects unprogrammed motion while enabling estimation of sensor knowledge errors in position, velocity, attitude angles, and attitude angle rates.

According to some exemplary embodiments, projections are formulated via bivariate polynomials curve fit to ray traced variations from a pinhole model. Differentiating these projection equations yields non-linear scan equations. A first-order Taylor series approximation reduces the scan equations to a coupled system of first-order differential equations with constant coefficients. Programmed motion positions across each detector junction are estimated via a state transition matrix solution. The actual motion positions are measured via correlation of the overlapping detectors and differenced with the predicted locations. A least-squares fit estimates updates to the scan coefficients based on these error residuals, to remove the unwanted motion.

According to some exemplary embodiments, error propagation is formulated via perturbations to the state transition matrix and forcing vector, which does not require that motion errors be estimated separately by category (rectilinear, oscillatory, and random); a robust and computationally efficient algorithm results. The rigorous error propagation model employed inherently characterizes statistical consistency, which enables identification of junction measurement “outliers,” and reinforces robustness. The formulation of scan coefficient updates enables estimation of knowledge errors in sensor position, velocity, attitude angles, and attitude angle rates. Accuracy is enhanced by an iteration process based on these knowledge errors.

According to some exemplary embodiments, a covariance model relates scan coefficient errors to focal plane coordinate errors for non-linear scan equations. According to these embodiments, a line scanning electro-optical sensor that does not follow a pinhole model is considered. Given all the support data values necessary to define the scan (sensor position and velocity, attitude angles and rates) and the accuracy with which these parameters are known (mean errors, covariance), the error statistics (means, covariance) in the projected location of a ground point's image at some small time increment later in the scan are computed, via closed-form expressions. These statistics are performance-important for the image formation chain; specifically, the formulation of an effective geometric correction process.

According to some exemplary embodiments, the propagated coordinate errors are derived in terms of first-order error perturbations to the scan coefficients appearing in the state transition matrix and forcing vector. These scan coefficient errors arise from un-modeled vibrations and knowledge errors in the parameters employed to compute the coefficients. These scan coefficient errors are then determined via the detector junction measurement residuals. Next, closed-form approximations for the junction measurement residuals are derived from the measurement errors in the platform support parameters. Combining these steps completes the solution.

According to some exemplary embodiments, this new approach of the present disclosure produces an accurate expression for predicting the errors in the projected focal plane scan coordinates of a non-pinhole camera system. This is a major improvement over the known method of Monte Carlo simulation, which is computationally cumbersome. Also, this accurate solution enables on-the-fly computation of performance quality measures for band-to-band registered exploitation products. Such metrics are in demand by the user community. Also, the approach described herein enables on-the-fly statistical consistency and outlier checks for the junction measurements.

According to some exemplary embodiments, platform error estimation for non-linear scan equations is also provided. Unknown amounts of system error are introduced by modeling electro-optical line scanning sensors as pinhole cameras. Accurate support data, including measurement histories for position, velocity, attitude angles and rates, and a priori error variances, is important for removing geometric artifacts in the images captured by these sensors. Projected focal plane coordinate errors are estimated during image formation via a correlation process in the detector junction overlap regions across the array. The problem solved is that of estimating knowledge errors in the support data from these measured coordinate errors; thus, improving image quality.

According to some exemplary embodiments, one first accurately approximates the resulting focal plane coordinate errors arising from given random errors in the scan coefficients, as would be measured by the detector junction correlation process. These scan coefficient errors arise because the sensor support data employed in their computation contains measurement errors. One then derives expressions relating the scan coefficient errors to the measurement errors in the platform parameters (position, velocity, attitude angles and rates). These two intermediate expressions are then combined into a least-squares estimator relating measured detector junction coordinate residuals to platform parameter knowledge errors.

According to some exemplary embodiments, this new approach produces accurate estimates for platform parameter knowledge errors on a line-by-line basis. The improved accuracy of the platform parameters will positively impact image quality. The error estimates obtained from the initial execution of this approach yields error estimates that can be applied to improve the original platform parameter measurement histories. Subsequently, the updated histories can be iterated through another pass of junction measurements to further tighten overall knowledge accuracy. The approach of the exemplary embodiments facilitates on-the-fly statistical consistency and outlier checks for the junction measurement process.

Previous electro-optical line scanning systems are based on a post-correction pinhole camera model. Specifically, for these systems it was assumed that either corrections for optical distortion, atmospheric refraction, etc., had been successfully executed prior to formulation of the focal plane scan equations, or these correction terms were insignificant. The resulting linear scan equations were employed to implement a multi-stage misregistration correction (MRC) approach. This disclosure describes an original method for implementing a MRC approach for an electro-optical line scanning sensor, which cannot be modeled as a pinhole camera.

The MRC approach of the present disclosure treats total error without the need for separately estimating linear, oscillatory and random errors. Its formulation is mathematically rigorous and computationally robust. The rigorous error propagation model employed in the approach described herein characterizes statistical consistency, which enables identification of junction measurement “outliers,” and reinforces robustness. Its formulation of scan coefficient updates enables estimation of knowledge errors in sensor position, velocity, attitude angles, and attitude angle rates via an inversion process. Accuracy is further enhanced by an iteration process based on these knowledge errors.

According to the present disclosure, an approach is presented for mathematically modeling the non-linear variations from a pinhole model present in ray trace data characterizing the sensor optics via bivariate polynomial curve fits. These polynomial variation models are employed to derive the focal plane scan equations, which define the focal plane trajectory traced by the image of a specific point in object space during the scan. The resulting scan equations are a pair of coupled 1st-order non-linear differential equations.

According to the exemplary embodiments, an accurate approximation is then introduced that reduces the non-linear differential equations to a set of linear differential equations, and which preserves the significant non-pinhole characteristics of the camera model. The solution to these approximate scan equations yield an error propagation model, which describes the flow of knowledge errors in the sensor position, velocity, attitude angles and rates into coordinate errors in the predicted position of an image point. These coordinate errors are measured via a correlation process in the detector junction overlap regions.

According to some exemplary embodiments, the junction measurements serve as inputs into a least-squares process, which is based on the error flow model. The least-squares process yields error estimates to the nominal scan coefficients, which are consistent with the measured coordinate residuals. It is these scan coefficient errors that are employed to affect MRC. An inversion process derives estimates of knowledge errors in the sensor position, velocity, attitude angles and rates from the scan coefficient error updates. The MRC process feeds back the knowledge errors into the a priori platform support data and then iterates the MRC solution.

A typical traditional MRC approach is now described as a reference point for the detailed description of the approach of the exemplary embodiments, which begins in section 1 below.

In terms of nomenclature, the term “programmed motion” used herein refers to commanded motion that is planned as part of the sensor's image acquisition process. Specifically, it is the motion of the image of a selected stationary ground point across the focal plane during a commanded scan. The terms “unprogrammed motion,” “image motion error” or “focal plane coordinate errors” refer to the motion of the image of a selected stationary ground point across the focal plane, which differs from the commanded motion, and can occur due to attitude errors, vibrations, or a number of other causes.

Many prior electro-optical line scanning sensor systems have been based on a post-correction pinhole camera model. Specifically, these systems assume that either corrections for optical distortion, atmospheric refraction, etc., have been successfully executed prior to formulation of the focal plane scan equations, or these correction terms were insignificant. The resulting linear scan equations were employed to implement a multi-stage MRC approach.

The scan equations were derived by direct differentiation of the pinhole camera projection equations with respect to time. All pinhole camera model parameters were allowed to vary with time except for the coordinates of the ground point of interest, camera focal length, and the camera principal point coordinates. The resulting scan equations were a coupled system of linear, first-order differential equations with six time varying coefficients. Ideally, one would supply the time histories of the sensor position, velocity, attitude angles, and attitude angle rates, compute the scan six time varying coefficients, and solve the system of differential equations.

In general these differential equations were found to be a difficult, if not impossible, problem to solve in closed form without simplifying assumptions. The simplifying assumption that was adopted is to consider a time increment sufficiently small to ensure that one can accurately consider the computed scan coefficients to be constant over this interval. One then solves the resulting equations and utilizes this solution until the time constraint has been reached. One then repeats the process until the entire time period has been covered increment-by-increment.

FIG. 1 includes a schematic diagram of a focal plane scanning geometry/configuration. The scan configuration 10 includes a forward array 12 of sensor modules 14 and a trailing array 16 of sensor modules 14. Each sensor module 14 includes a plurality of sensors 28 disposed linearly. The sensor modules 14 of the forward array 12 are staggered with respect to the sensor modules 14 of the trailing array 16, such that an overlap 20 of sensor modules is created. The forward array 12 and trailing array 16 are offset by an offset 18. Through appropriate processing of received sensor data, the forward array 12 and trailing array 16 combine to form a synthetic array 26, as shown in phantom in FIG. 1. During scanning image data collection, an object 24 is first detected by the forward array 12, and, due to the motion of the sensor platform, the object 24 is then detected by trailing array 16. The object 24 follows a scan trajectory 22 along the arrays 12 and 16. Previous approaches involved electro-optical line scanning sensors with offset and overlapping detector placements, as illustrated in FIG. 1.

The solution to the scan equations was employed to predict the object crossing points of the leading (forward) array 12 and the trailing array 16 of the scan trajectory 22 of an image point. The correlated detector responses of the forward and trailing arrays within each overlap region 20 determined the actual shift in positions of the crossing points. Variations between the predicted and the measured shifts would be used to estimate the impact of any unplanned motion present, and correct for it during formation of the synthetic image array.

One exemplary detector junction correlation process is illustrated in FIG. 2, which is a schematic functional flow diagram of a detector junction correlation process. The correlation process relies on the fact that the leading and trailing arrays “see” the same image but at different times. The correlator can be a normalized cross correlator, a lag product cross correlator, a least squares cross correlator, or other known cross correlator. The image correlator is used to determine where a selected portion of the image viewed by the leading array 12 best aligns with the image seen by the trailing array 16. This provides a measured set of focal plane coordinates that can be compared with coordinates that were predicted via the scan equation solution using support data parameters.

In general, the highest measurable spatial frequency of a motion disturbance is a function of line rate and array separation. FIG. 3 includes a schematic functional block diagram of the detector junction measurement process 50. According to the process 50, the difference between the measured and predicted focal plane coordinates feed the MRC process. Referring to FIG. 3, scan coefficients are computed 54, using camera attitude and position data. Using the computed scan coefficients, programmed trajectory is also computed 56. Using the computed programmed trajectory, the trailing array crossing point is estimated 58. In parallel, also using camera attitude and position data as well as the junction overlap data, junction overlaps are correlated 60. Using the correlation of the junction overlaps, the trailing array crossing point is measured 62. The differences or residuals between the trailing array crossing estimate 58 and the measured trailing array crossing point are computed 64. As noted above, these differences or residuals are fed to the MRC process 52.

Thus, in previous approaches, the focal plane coordinate error residuals obtained from the detector junction measurements were exploited to estimate the amount of unprogrammed motion due to unwanted linear, oscillatory, and random motion. The quantified amount of unprogrammed motion in these categories was then “corrected” when forming the synthetic image array: The unwanted displacements were subtracted out when mapping the synthetic array interpolation locations.

In the following sections numbered 1 through 7, preferred embodiments which overcome drawbacks and disadvantages of the previous approaches are described in detail.

1 Non-Pinhole Camera Variation Models

FIG. 4 is a schematic diagram of traditional pinhole camera model and projection equations. The pinhole camera model yields linear scan equations. However, it is unknown how much error is introduced into the system by using this traditional simplified model. It would be beneficial to construct a system of non-linear scan equations to evaluate the errors introduced by the traditional pinhole camera model. A non-linear model of the projection equations is numerically determined from ray trace data in accordance with the following.

FIG. 5 includes a schematic logical flow diagram of generating correction terms to be added to the traditional pinhole camera model, according to some exemplary embodiments. Referring to FIG. 5, ray traced points 302 through the optics are compared with locations 304 predicted from a pinhole model. The position residuals are then least-squares curve fit 306 to bivariate polynomial expressions to obtain correction terms 308 to be added to the traditional pinhole model. The order of the polynomial fit required to accurately model the variation from a pinhole camera may vary from design to design but the principle is the same. Once an accurate polynomial representation has been found; the numerical coefficients of the polynomial fit characterize the system's variation from a pinhole camera. In designs that allow for focus variation, it may be necessary to generate a set of polynomials for several different focus ranges.

The MRC approach described herein will work for polynomials of any degree. However, in order to avoid unnecessary abstraction and to make the mathematics easier to follow, the following quadratic bivariate polynomials are featured:
Δx=a ₀₀ +a ₀₁ y+a ₀₂ y ² +a ₁₀ x+a ₁₁ xy+a ₂₀ x ²  Eq. 1-1a
Δx=b ₀₀ +b ₀₁ y+b ₀₂ y ² +b ₁₀ x+b ₁₁ xy+b ₂₀ x ²  Eq. 1-1b
This quadratic model provides six coefficients to represent the variation from a pinhole camera in each of two coordinates; for a total of twelve real-valued coefficients.
Hence, the non-pinhole camera projection equations are:

$\begin{array}{cc}x-x_p+f\frac{L}{N}+a_{00}+a_{01}y+a_{02}{y}^2+a_{10}x+a_{11}\mathrm{xy}+a_{20}{x}^2=0& \mathrm{Eq}.1\text{-}2a\\ y-y_p+f\frac{M}{N}+b_{00}+b_{01}y+b_{02}{y}^2+{\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}x+b_{11}\mathrm{xy}+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}{x}^2=0& \mathrm{Eq}.1\text{-}2b\end{array}$
2 Non-Linear Scan Equations

The scan equations are derived by direct differentiation of the non-pinhole camera projection equations with respect to time. All the parameters of the non-pinhole camera model are allowed to vary with time except for the coordinates of the ground point of interest, camera focal length, the polynomial coefficients of the variation models, and the camera principal point coordinates.

Differentiating equation 1-2 with respect to time yields:

$\begin{array}{cc}\frac{dx}{dt}+\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]+a_{01}\frac{dy}{dt}+2a_{02}y\frac{dy}{dt}+a_{10}\frac{dx}{dt}+a_{11}\frac{dx}{dt}y+a_{11}x\frac{dy}{dt}+a_{20}2x\frac{dx}{dt}=0& \mathrm{Eq}.2\text{-}1a\\ \frac{dy}{dt}+\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]+b_{01}\frac{dy}{dt}+2b_{02}y\frac{dy}{dt}+{\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}\frac{dx}{dt}+b_{11}\frac{dx}{dt}y+b_{11}x\frac{dy}{dt}+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}2x\frac{dx}{dt}=0& \mathrm{Eq}.2\text{-}1b\end{array}$
These expressions may be immediately simplified to:

$\begin{array}{cc}\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]+\frac{dy}{dt}\left(a_{01}+2a_{02}y+a_{11}x\right)+\frac{dx}{dt}\left(1+a_{10}+a_{11}y+a_{20}2x\right)=0& \mathrm{Eq}.2\text{-}2a\\ \left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]+\frac{dy}{dt}\left(1+b_{01}+2b_{02}y+b_{11}x\right)+\frac{dx}{dt}\left({\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}+b_{11}y+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}2x\right)=0& \mathrm{Eq}.2\text{-}2b\end{array}$
Rearranging terms yields:

$\begin{array}{cc}\left(a_{01}+2a_{02}y+a_{11}x\right)\frac{dy}{dt}+\left(1+a_{10}+a_{11}y+a_{20}2x\right)\frac{dx}{dt}=-\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]& \mathrm{Eq}.2\text{-}3a\\ \left(1+b_{01}+2b_{02}y+b_{11}x\right)\frac{dy}{dt}+\left({\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}+b_{11}y+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}2x\right)\frac{dx}{dt}=-\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]& \mathrm{Eq}.2\text{-}3b\end{array}$
Equations 2-3a and 2-3b may be rewritten in matrix form as:

$\begin{array}{cc}\left(\begin{array}{cc}\left(1+a_{10}+a_{11}y+a_{20}2x\right)& \left(a_{01}+2a_{02}y+a_{11}x\right)\\ \left({\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}+b_{11}y+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}2x\right)& \left(1+b_{01}+2b_{02}y+b_{11}x\right)\end{array}\right)\left(\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right)=\left(\begin{array}{c}-\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\\ -\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\end{array}\right)& \mathrm{Eq}.2\text{-}4\end{array}$
Thus, assuming the indicated matrix is invertible, one may write that:

$\begin{array}{cc}\left(\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right)={\left(\begin{array}{cc}\left(1+a_{10}+a_{11}y+a_{20}2x\right)& \left(a_{01}+2a_{02}y+a_{11}x\right)\\ \left({\mathrm{<figure-callout\; id="10"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{10}+b_{11}y+{\mathrm{<figure-callout\; id="20"\; label="b"\; filenames="US09521326-20161213-D00001.png"\; state="\{\{state\}\}">b</figure-callout>}}_{20}2x\right)& \left(1+b_{01}+2b_{02}y+b_{11}x\right)\end{array}\right)}^{-1}\left(\begin{array}{c}-\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\\ -\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\end{array}\right)& \mathrm{Eq}.2\text{-}5\end{array}$

Determining whether or not the matrix is invertible would require specification of particular values for the polynomial coefficients. However, a quick glance at the matrix reveals that its determinate is a quadratic bivariate polynomial. Setting this expression equal to zero will result in a quadratic plane curve. The near linear, short-scan trajectories of interest here are likely only to cross this forbidden curve in at most two places. FIG. 6 includes a schematic diagram of an approach to working around this singularity issue, according to some exemplary embodiments. Referring to FIG. 6, one can “work around” these singularity problem points by stopping the projected scan trajectory at a stop point just short of the singular point and extending the solution past it by a simple linear segment whose slope matches the scan trajectory's first derivative at the stop point, i.e., a tangent line. The solution is then picked up after a small time increment, i.e., across an exclusion band, on the other side.

For a 2×2 matrix, the inverse is given by:

$\begin{array}{cc}{\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)}^{-1}=\frac{1}{\mathrm{ad}-\mathrm{bc}}\left(\begin{array}{cc}d& -b\\ -c& a\end{array}\right)& \mathrm{Eq}.2\text{-}6\end{array}$
Consequently,

$\begin{array}{cc}\left(\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right)=\left(\frac{1}{E\left(x,y\right)}\right)\left(\begin{array}{cc}D\left(x,y\right)& -B\left(x,y\right)\\ -C\left(x,y\right)& A\left(x,y\right)\end{array}\right)\left(\begin{array}{c}-\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\\ -\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]\end{array}\right)& \mathrm{Eq}.2\text{-}7a\\ A\left(x,y\right)=\left(1+a_{10}+a_{11}y+a_{20}2x\right)& \mathrm{Eq}.2\text{-}7b\\ B\left(x,y\right)=\left(a_{01}+2a_{02}y+a_{11}x\right)& \mathrm{Eq}.2\text{-}7c\\ C\left(x,y\right)=\left(b_{10}+b_{11}y+b_{20}2x\right)& \mathrm{Eq}.2\text{-}7d\\ D\left(x,y\right)=\left(1+b_{01}+2b_{02}y+b_{11}x\right)& \mathrm{Eq}.2\text{-}7e\\ E\left(x,y\right)=A\left(x,y\right)D\left(x,y\right)-B\left(x,y\right)C\left(x,y\right)& \mathrm{Eq}.2\text{-}7f\end{array}$
Now, let

$\begin{array}{cc}P\left(t\right)=-\left(f\frac{L}{N}\right)\left[\frac{1}{L}\frac{dL}{dt}-\frac{1}{N}\frac{dN}{dt}\right]& \mathrm{Eq}.2\text{-}8a\\ Q\left(t\right)=-\left(f\frac{M}{N}\right)\left[\frac{1}{M}\frac{dM}{dt}-\frac{1}{N}\frac{dN}{dt}\right]& \mathrm{Eq}.2\text{-}8b\end{array}$
Equation 2-7a then becomes

$\begin{array}{cc}\left(\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right)=\left(\begin{array}{c}\frac{D\left(x,y\right)}{E\left(x,y\right)}P\left(t\right)-\frac{B\left(x,y\right)}{E\left(x,y\right)}Q\left(t\right)\\ -\frac{C\left(x,y\right)}{E\left(x,y\right)}P\left(t\right)+\frac{A\left(x,y\right)}{E\left(x,y\right)}Q\left(t\right)\end{array}\right)& \mathrm{Eq}.2\text{-}9\end{array}$

Thus the scan equation 2-9 is a couple system of first order, non-linear differential equations. Practically speaking, such equations can rarely be solved in closed form. Consequently, an accurate approximation is required that can (1) capture the significant non-pinhole character of the projection equations and (2) that leads to a soluble system. Such an approximation is provided in the next section.

3 An Approximate Solution

Exemplary desired initial conditions are selected: a point in the focal plane (x₀,y₀) at time t₀. There exists a time increment and associated neighborhood of N(x₀,y₀) for which one may approximate to the first order (bivariate Taylor's series):

$\begin{array}{cc}P\left(t\right)\approx P\left(t_0\right)t\in \left[t_0,t_0+\mathrm{\Delta \tau }\right]& \mathrm{Eq}.3\text{-}1a\\ Q\left(t\right)\approx Q\left(t_0\right)t\in \left[t_0,t_0+\mathrm{\Delta \tau }\right]& \mathrm{Eq}.3\text{-}1b\\ N\left(x_0,y_0\right)=\left\{\left(x,y\right)|\sqrt{{\left(x-x_0\right)}^2+{\left(x-x_0\right)}^2}\le R_0\right\}& \mathrm{Eq}.3\text{-}1c\\ \frac{D\left(x,y\right)}{E\left(x,y\right)}\approx \frac{D\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}+\left[\frac{\partial }{\partial x}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(x-x_0\right)+\hspace{1em}\left[\frac{\partial }{\partial y}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(y-y_0\right)& \mathrm{Eq}.3\text{-}1d\\ \frac{B\left(x,y\right)}{E\left(x,y\right)}\approx \frac{B\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}+\left[\frac{\partial }{\partial x}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(x-x_0\right)+\hspace{1em}\left[\frac{\partial }{\partial y}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(y-y_0\right)& \mathrm{Eq}.3\text{-}1e\\ \frac{C\left(x,y\right)}{E\left(x,y\right)}\approx \frac{C\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}+\left[\frac{\partial }{\partial x}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(x-x_0\right)+\hspace{1em}\left[\frac{\partial }{\partial y}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(y-y_0\right)& \mathrm{Eq}.3\text{-}1f\\ \frac{A\left(x,y\right)}{E\left(x,y\right)}\approx \frac{A\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}+\left[\frac{\partial }{\partial x}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(x-x_0\right)+\hspace{1em}\left[\frac{\partial }{\partial y}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]\left(y-y_0\right)& \mathrm{Eq}.3\text{-}1g\end{array}$
Substituting Equations 3-1a through 3-1g into Equation 2-9 and setting

$\begin{array}{cc}{\Gamma }_1=\left[\frac{\partial }{\partial x}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)-\left[\frac{\partial }{\partial x}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)& \mathrm{Eq}.3\text{-}2a\\ {\Gamma }_2=\left[\frac{\partial }{\partial y}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)-\left[\frac{\partial }{\partial y}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)& \mathrm{Eq}.3\text{-}2b\\ {\Gamma }_3=\frac{D\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}P\left(t_0\right)-\frac{B\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}Q\left(t_0\right)-\left[\frac{\partial }{\partial x}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)x_0+\left[\frac{\partial }{\partial x}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)x_0-\left[\frac{\partial }{\partial y}{\left(\frac{D\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)y_0+\left[\frac{\partial }{\partial y}{\left(\frac{B\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)y_0& \mathrm{Eq}.3\text{-}2c\\ {\Gamma }_4=-\left[\frac{\partial }{\partial x}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)+\left[\frac{\partial }{\partial x}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)& \mathrm{Eq}.3\text{-}2d\\ {\Gamma }_5=-\left[\frac{\partial }{\partial y}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)+\left[\frac{\partial }{\partial y}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)& \mathrm{Eq}.3\text{-}2e\\ {\Gamma }_6=-\frac{C\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}P\left(t_0\right)+\frac{A\left(x_0,y_0\right)}{E\left(x_0,y_0\right)}Q\left(t_0\right)+\left[\frac{\partial }{\partial x}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)x_0-\left[\frac{\partial }{\partial x}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)x_0+\left[\frac{\partial }{\partial y}{\left(\frac{C\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]P\left(t_0\right)y_0-\left[\frac{\partial }{\partial y}{\left(\frac{A\left(x,y\right)}{E\left(x,y\right)}\right)}_{\left(x_0,y_0\right)}\right]Q\left(t_0\right)y_0& \mathrm{Eq}.3\text{-}2f\end{array}$
Allows the scan equation 2-9 to be rewritten in the form:

$\begin{array}{cc}\left(\begin{array}{c}\frac{dx}{dt}\\ \frac{dy}{dt}\end{array}\right)=\left(\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_{1}x+{\Gamma }_{2}y+{\Gamma }_3\\ {\Gamma }_{4}x+{\Gamma }_{5}y+{\Gamma }_6\end{array}\right)=\left(\begin{array}{cc}{\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1& {\Gamma }_2\\ {\Gamma }_4& {\Gamma }_5\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{\Gamma }_3\\ {\Gamma }_6\end{array}\right)& \mathrm{Eq}.3\text{-}2g\end{array}$
Here Γ₁, . . . Γ₆ are constant coefficients. Hence, the approximate scan equations constitute a coupled system of first order differential equations with constant coefficients. One may now rewrite the scan equations in matrix form as:

$\begin{array}{cc}\frac{d}{dt}\stackrel{\rightharpoonup }{s}\left(t\right)={\Lambda }_{2\times 2}\stackrel{\rightharpoonup }{s}\left(t\right)+\stackrel{\rightharpoonup }{F}& \mathrm{Eq}.3\text{-}3a\\ \stackrel{\rightharpoonup }{s}\left(t\right)=\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right)& \mathrm{Eq}.3\text{-}3b\\ {\Lambda }_{2\times 2}=\left(\begin{array}{cc}{\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1& {\Gamma }_2\\ {\Gamma }_4& {\Gamma }_5\end{array}\right)& \mathrm{Eq}.3\text{-}3c\\ \stackrel{\rightharpoonup }{F}=\left(\begin{array}{c}{\Gamma }_3\\ {\Gamma }_6\end{array}\right)& \mathrm{Eq}.3\text{-}3d\end{array}$

It is noted that each of the scan coefficients Γ₁, . . . Γ₆ depends on all the polynomial coefficients of the variation models, and thus these scan coefficients reflect the significant contributions of the non-pinhole projection equations.

There are numerous mathematical techniques for solving such a system including: matrix methods, the Laplace transform, method of elimination, etc. Employing matrix methods one can show that the solution to equation 3-3 is:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\Lambda \left(t-t_0\right)}\stackrel{\rightharpoonup }{s}\left(t_0\right)+{\int }_{t_0}^t{e}^{\Lambda \left(t-\zeta \right)}\stackrel{\rightharpoonup }{F}d\zeta & \mathrm{Eq}.3\text{-}4\end{array}$
It is noted that the function e^{Λ(t−t 0 )} yields a matrix upon evaluation, which in linear systems theory is referred to as the state transition matrix. One can represent this function of a matrix in series form as:

$\begin{array}{cc}{e}^{\Lambda t}=\sum _{p=0}^{\infty }\frac{{\Lambda }^p{t}^p}{p!}& \mathrm{Eq}.3\text{-}5\end{array}$
As the exponential is an entire function in the complex plane, then the series in equation 3-5 converges for all complex 2×2 matrices.

Evaluation of the exponential function of the matrix Λ may be accomplished in closed form by employing the spectral resolution theorem:

Spectral Resolution Theorem:

If A is a member of the set of n×n complex matrices, ƒ is a function defined on the spectrum of A, ƒ_k ^(j) is the value of the j^th derivative of ƒ at the eigenvalue λ_k (k=1, 2 . . . , s, j=0, . . . m_k−1), and m_k is the index of λ_k, then there exists matrices Z_kj, independent of ƒ, such that:

$\begin{array}{cc}F\left(A\right)=\sum _{k=1}^s\sum _{j=1}^{m_k}{f}_k^{\left(j-1\right)}{Z}_{\mathrm{kj}}& \mathrm{Eq}.3\text{-}6\end{array}$
Moreover, the matrices Z_kj are linearly independent n×n complex matrices and commute with A.
Recall that the index m_k of an eigenvalue λ_k of the matrix Λ is the power of the term (λ−λ_k)^{m k} appearing in the minimal annihilating polynomial of Λ.
Successful application of this theorem requires that one first:

• • Determine the eigenvalues of the matrix Λ.
  • Compute the matrices Z_kj.
    The eigenvalues are determined via the characteristic equation:

$\begin{array}{cc}\det \left(\begin{array}{cc}{\Gamma }_1-\lambda & {\Gamma }_2\\ {\Gamma }_4& {\Gamma }_5-\lambda \end{array}\right)=0& \mathrm{Eq}.3\text{-}7\end{array}$
The resulting characteristic equation is quadratic:
λ²−(Γ₁+Γ₅)λ+(Γ₁Γ₅−Γ₂Γ₄)=0
One finds the eigenvalues to be:

$\begin{array}{cc}{\mathrm{<figure-callout\; id="1"\; label="\lambda "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\lambda </figure-callout>}}_1=\frac{\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1+{\Gamma }_5\right)+\sqrt{{\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1+{\Gamma }_5\right)}^2-4\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1{\Gamma }_5-{\Gamma }_2{\Gamma }_4\right)}}{2}& \mathrm{Eq}.3\text{-}8a\\ {\lambda }_2=\frac{\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1+{\Gamma }_5\right)-\sqrt{{\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1+{\Gamma }_5\right)}^2-4\left({\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1{\Gamma }_5-{\Gamma }_2{\Gamma }_4\right)}}{2}& \mathrm{Eq}.3\text{-}8b\end{array}$

This results in three distinct cases depending on the sign of the discriminant D:
D=(Γ₁+Γ₅)²−4(Γ₁Γ₅−Γ₂Γ₄)  Eq. 3-9

• • Case 1: D>0 Two distinct real eigenvalues
  • Case 2: D<0 Two distinct complex conjugate eigenvalues
  • Case 3: D=0 One repeated real eigenvalue

The function e^Λt is known for each of these cases:

$\begin{array}{cc}\mathrm{Case}1\text{:}D>0\mathrm{Two}\mathrm{distinct}\mathrm{real}\mathrm{eigenvalues}{e}^{\Lambda \left(t-t_0\right)}=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{cc}\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}& {\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\\ {\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)& \begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\end{array}\right)& \mathrm{Eq}.3\text{-}10a\\ \alpha =\frac{1}{2}\left({\Gamma }_1+{\Gamma }_5\right)& \mathrm{Eq}.3\text{-}10b\\ \beta =\frac{1}{2}\sqrt{{\left({\Gamma }_1+{\Gamma }_5\right)}^2-4\left({\Gamma }_1{\Gamma }_5-{\Gamma }_2{\Gamma }_4\right)}& \mathrm{Eq}.3\text{-}10c\\ \varepsilon =\frac{1}{2}\left({\Gamma }_1-{\Gamma }_5\right)& \mathrm{Eq}.3\text{-}10d\\ \mathrm{Case}2\text{:}D<0\mathrm{Two}\mathrm{distinct}\mathrm{complex}\mathrm{conjugate}\mathrm{eigenvalues}{e}^{\Lambda \left(t-t_0\right)}=\frac{e^{\alpha \left(t-t_0\right)}}{\tilde{\beta }}\left(\begin{array}{cc}\begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)+\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}& {\Gamma }_2\sinh \left(\tilde{\beta }\left(t-t_0\right)\right)\\ {\Gamma }_4\sinh \left(\tilde{\beta }\left(t-t_0\right)\right)& \begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)-\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}\end{array}\right)& \mathrm{Eq}.3\text{-}11a\\ \tilde{\beta }=\frac{1}{2}\sqrt{-{\left({\Gamma }_1+{\Gamma }_5\right)}^2+4\left({\Gamma }_1{\Gamma }_5-{\Gamma }_2{\Gamma }_4\right)}& \mathrm{Eq}.3\text{-}11b\\ \mathrm{Case}3\text{:}D=0\mathrm{One}\mathrm{repeated}\mathrm{real}\mathrm{eigenvalue}{e}^{\Lambda \left(t-t_0\right)}=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{cc}1+\varepsilon \left(t-t_0\right)& {\Gamma }_2\left(t-t_0\right)\\ {\Gamma }_4\left(t-t_0\right)& 1+\varepsilon \left(t-t_0\right)\end{array}\right)& \mathrm{Eq}.3\text{-}12\end{array}$

Hence, the approximate solution is by case:

$\begin{array}{cc}\mathrm{Case}1\text{:}D>0\mathrm{Two}\mathrm{distinct}\mathrm{real}\mathrm{eigenvalues}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{cc}\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}& {\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\\ {\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)& \begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\end{array}\right)\stackrel{\rightharpoonup }{s}\left(t_0\right)+{\int }_{t_0}^t\frac{e^{\alpha \left(t-\zeta \right)}}{\beta }\left(\begin{array}{cc}\begin{array}{c}\beta \cosh \left(\beta \left(t-ϛ\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-ϛ\right)\right)\end{array}& {\Gamma }_2\sinh \left(\beta \left(t-ϛ\right)\right)\\ {\Gamma }_4\sinh \left(\beta \left(t-ϛ\right)\right)& \begin{array}{c}\beta \cosh \left(\beta \left(t-ϛ\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-ϛ\right)\right)\end{array}\end{array}\right)\stackrel{\rightharpoonup }{F}d\zeta & \mathrm{Eq}.3\text{-}13\end{array}$
Hence, expanding terms

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{c}\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+{\int }_{t_0}^t\left(\frac{e^{\alpha \left(t-\zeta \right)}}{\beta }\right)\left(\begin{array}{c}\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-ϛ\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-ϛ\right)\right)\end{array}\right]{\Gamma }_3+\left[{\Gamma }_2\sinh \left(\beta \left(t-ϛ\right)\right)\right]{\Gamma }_6\\ \left[{\Gamma }_4\sinh \left(\beta \left(t-ϛ\right)\right)\right]{\Gamma }_3+\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-ϛ\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-ϛ\right)\right)\end{array}\right]{\Gamma }_6\end{array}\right)d\zeta & \mathrm{Eq}.3\text{-}14\end{array}$
Observe that one may rewrite the integral in equation 3-14 to obtain:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{c}\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+{\int }_0^{t-{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}\left(\frac{e^{\alpha z}}{\beta }\right)\left(\begin{array}{c}\left[{\mathrm{\beta \Gamma }}_3\cosh \left(\beta z\right)\right]+\left[\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)\sinh \left(\beta z\right)\right]\\ \left[{\mathrm{<figure-callout\; id="6"\; label="\beta \Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\beta \Gamma </figure-callout>}}_6\cosh \left(\beta z\right)\right]+\left[\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)\sinh \left(\beta z\right)\right]\end{array}\right)dz& \mathrm{Eq}.3\text{-}15\end{array}$
Direct integration establishes that:

$\begin{array}{cc}\int e^{\alpha z}\cosh \left(\beta z\right)dz=\frac{1}{2}\int e^{\left(\alpha +\beta \right)z}dz+\frac{1}{2}\int e^{\left(\alpha -\beta \right)z}dz=\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}+\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}& \mathrm{Eq}.3\text{-}16a\\ \int e^{\alpha z}\sinh \left(\beta z\right)dz=\frac{1}{2}\int e^{\left(\alpha +\beta \right)z}dz-\frac{1}{2}\int e^{\left(\alpha -\beta \right)z}dz=\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}-\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}& \mathrm{Eq}.3\text{-}16b\end{array}$
Substituting equations 3-16a and 3-16b into 3-15 yields:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{c}\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+\left(\begin{array}{c}{{\Gamma }_3\left[\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}+\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}\right]}_0^{t-t_0}+{\frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\beta }\left[\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}-\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}\right]}_0^{t-t_0}\\ {{\Gamma }_6\left[\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}+\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}\right]}_0^{t-t_0}+{\frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\beta }\left[\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}-\frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}\right]}_0^{t-t_0}\end{array}\right)& \mathrm{Eq}.3\text{-}17\end{array}$
Consequently, the solution to Case 1 (D>0 two distinct real eigenvalues) is:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\beta }\left(\begin{array}{c}\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)+\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sinh \left(\beta \left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sinh \left(\beta \left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\beta \cosh \left(\beta \left(t-t_0\right)\right)-\\ \varepsilon \sinh \left(\beta \left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+\left(\begin{array}{c}\begin{array}{c}{\Gamma }_3\left[\frac{e^{\left(\alpha +\beta \right)\left(t-t_0\right)}}{2\left(\alpha +\beta \right)}+\frac{e^{\left(\alpha -\beta \right)\left(t-t_0\right)}}{2\left(\alpha -\beta \right)}\right]+\frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\beta }\left[\begin{array}{c}\frac{e^{\left(\alpha +\beta \right)z}}{2\left(\alpha +\beta \right)}-\\ \frac{e^{\left(\alpha -\beta \right)z}}{2\left(\alpha -\beta \right)}\end{array}\right]-\frac{{\Gamma }_3}{2}\left[\frac{1}{\left(\alpha +\beta \right)}+\frac{1}{\left(\alpha -\beta \right)}\right]-\\ \frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\beta }\left[\frac{1}{\left(\alpha +\beta \right)}+\frac{1}{\left(\alpha -\beta \right)}\right]\end{array}\\ \begin{array}{c}{\Gamma }_6\left[\frac{e^{\left(\alpha +\beta \right)\left(t-t_0\right)}}{2\left(\alpha +\beta \right)}+\frac{e^{\left(\alpha -\beta \right)\left(t-t_0\right)}}{2\left(\alpha -\beta \right)}\right]+\frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\beta }\left[\begin{array}{c}\frac{e^{\left(\alpha +\beta \right)\left(t-t_0\right)}}{2\left(\alpha +\beta \right)}-\\ \frac{e^{\left(\alpha -\beta \right)\left(t-t_0\right)}}{2\left(\alpha -\beta \right)}\end{array}\right]-\frac{{\Gamma }_3}{2}\left[\frac{1}{\left(\alpha +\beta \right)}+\frac{1}{\left(\alpha -\beta \right)}\right]-\\ \frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\beta }\left[\frac{1}{\left(\alpha +\beta \right)}+\frac{1}{\left(\alpha -\beta \right)}\right]\end{array}\end{array}\right)& \mathrm{Eq}.3\text{-}18\end{array}$

Case 2: D<0 Two distinct complex conjugate eigenvalues

Given the analogy between equations 3-10a and 3-11a, one may write immediately via analogy with equation 3-15:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\tilde{\beta }}\left(\begin{array}{c}\left[\begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)+\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)-\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+{\int }_0^{t-t_0}\left(\frac{e^{\alpha z}}{\tilde{\beta }}\right)\left(\begin{array}{c}\left[\tilde{\beta }{\Gamma }_3\cos \left(\tilde{\beta }z\right)\right]+\left[\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)\sin \left(\tilde{\beta }z\right)\right]\\ \left[\tilde{\beta }{\Gamma }_6\cos \left(\tilde{\beta }z\right)\right]+\left[\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)\sin \left(\tilde{\beta }z\right)\right]\end{array}\right)dz& \mathrm{Eq}.3\text{-}19\end{array}$
One can show that

$\begin{array}{cc}\int e^{\alpha z}\cos \left(\tilde{\beta }z\right)dz=\frac{e^{\alpha z}\left(\alpha \cos \left(\tilde{\beta }z\right)+\tilde{\beta }\sin \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}& \mathrm{Eq}.3\text{-}20a\\ \int e^{\alpha z}\sin \left(\tilde{\beta }z\right)dz=\frac{e^{\alpha z}\left(\alpha \sin \left(\tilde{\beta }z\right)-\tilde{\beta }\cos \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}& \mathrm{Eq}.3\text{-}20b\end{array}$
Substituting equations 3-20a and 3-20b into 3-19 yields:

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\tilde{\beta }}\left(\begin{array}{c}\left[\begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)+\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}\right]x_0+\left[{\Gamma }_2\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]x_0+\left[\begin{array}{c}\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)-\\ \varepsilon \sin \left(\tilde{\beta }\left(t-t_0\right)\right)\end{array}\right]y_0\end{array}\right)+\left(\begin{array}{c}{{\Gamma }_3\left[\frac{e^{\alpha z}\left(\alpha \cos \left(\tilde{\beta }z\right)+\tilde{\beta }\sin \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]}_0^{t-t_0}+{\frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\tilde{\beta }}\left[\frac{e^{\alpha z}\left(\alpha \sin \left(\tilde{\beta }z\right)-\tilde{\beta }\cos \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]}_0^{t-t_0}\\ {{\Gamma }_6\left[\frac{e^{\alpha z}\left(\alpha \cos \left(\tilde{\beta }z\right)+\tilde{\beta }\sin \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]}_0^{t-t_0}+{\frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\tilde{\beta }}\left[\frac{e^{\alpha z}\left(\alpha \sin \left(\tilde{\beta }z\right)-\tilde{\beta }\cos \left(\tilde{\beta }z\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]}_0^{t-t_0}\end{array}\right)& \mathrm{Eq}.3\text{-}21\end{array}$
Hence, one finds the solution for Case 2 (D<0 Two distinct complex conjugate eigenvalues) to be:

$\begin{array}{cc}s\left(t\right)=\frac{e^{\alpha \left(t-t_0\right)}}{\tilde{\beta }}(\begin{array}{c}\left[\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)+\mathrm{\varepsilon sin}\left(\tilde{\beta }\left(t-t_0\right)\right)\right]x_0+\left[{\Gamma }_2\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]y_0\\ \left[{\Gamma }_4\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right]x_0+\left[\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)-\mathrm{\varepsilon sin}\left(\tilde{\beta }\left(t-t_0\right)\right)\right]y_0\end{array})+\hspace{1em}\left(\begin{array}{c}\begin{array}{c}{\Gamma }_3\left[\frac{e^{\alpha z}\left(\mathrm{\alpha cos}\left(\tilde{\beta }\left(t-t_0\right)\right)+\tilde{\beta }\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]+\\ \frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\tilde{\beta }}\left[\frac{e^{\alpha z}\left(\mathrm{\alpha sin}\left(\tilde{\beta }\left(t-t_0\right)\right)-\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)\right)}{{\alpha }^2+{\beta }^2}\right]-\\ \frac{{\Gamma }_3}{{\alpha }_2+{\tilde{\beta }}^2}-\frac{\left({\mathrm{\varepsilon \Gamma }}_3+{\Gamma }_2{\Gamma }_6\right)}{\tilde{\beta }\left({\alpha }^2+{\tilde{\beta }}^2\right)}\end{array}\\ \begin{array}{c}{\Gamma }_6\left[\frac{e^{\alpha z}\left(\mathrm{\alpha cos}\left(\tilde{\beta }\left(t-t_0\right)\right)+\tilde{\beta }\sin \left(\tilde{\beta }\left(t-t_0\right)\right)\right)}{{\alpha }^2+{\tilde{\beta }}^2}\right]+\\ \frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\tilde{\beta }}\left[\frac{e^{\alpha z}\left(\mathrm{\alpha sin}\left(\tilde{\beta }\left(t-t_0\right)\right)-\tilde{\beta }\cos \left(\tilde{\beta }\left(t-t_0\right)\right)\right)}{{\alpha }^2+{\beta }^2}\right]-\\ \frac{{\Gamma }_6}{{\alpha }_2+{\tilde{\beta }}^2}-\frac{\left({\Gamma }_4{\Gamma }_3-{\mathrm{\varepsilon \Gamma }}_6\right)}{\tilde{\beta }\left({\alpha }^2+{\tilde{\beta }}^2\right)}\end{array}\end{array}\right)& \mathrm{Eq}.3\text{-}22\end{array}$

Case 3: D=0 One repeated real eigenvalue

Substituting equation 3-12 into 3-4 yields

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{cc}1+\varepsilon \left(t-t_0\right)& {\Gamma }_2\left(t-t_0\right)\\ {\Gamma }_4\left(t-t_0\right)& 1-\varepsilon \left(t-t_0\right)\end{array}\right)\stackrel{\rightharpoonup }{s}\left(t_0\right)+\underset{t_0}{\overset{t}{\int }}{e}^{\alpha \left(t-\zeta \right)}\left(\begin{array}{cc}1+\varepsilon \left(t-\zeta \right)& {\Gamma }_2\left(t-\zeta \right)\\ {\Gamma }_4\left(t-\zeta \right)& 1-\varepsilon \left(t-\zeta \right)\end{array}\right)\stackrel{\rightharpoonup }{F}d\zeta & \mathrm{Eq}.3\text{-}23\end{array}$
Hence, one may write that;

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{c}\left[1+\varepsilon \left(t-t_0\right)\right]x_0+{\Gamma }_2\left(t-t_0\right){\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">y</figure-callout>}}_0\\ {\Gamma }_4\left(t-t_0\right)x_0+\left[1-\varepsilon \left(t-t_0\right)\right]y_0\end{array}\right)+\underset{0}{\overset{t-{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}{\int }}{e}^{\alpha \left(t-\zeta \right)}\left(\begin{array}{c}\left[1+\varepsilon \left(t-\zeta \right)\right]{\Gamma }_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\left(t-\zeta \right)\\ {\Gamma }_3{\Gamma }_4\left(t-\zeta \right)+\left[1-\varepsilon \left(t-\zeta \right)\right]{\Gamma }_6\end{array}\right)d\zeta & \mathrm{Eq}.3\text{-}24\end{array}$
A change of variables yields;

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{c}\left[1+\varepsilon \left(t-t_0\right)\right]x_0+{\Gamma }_2\left(t-t_0\right){\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">y</figure-callout>}}_0\\ {\Gamma }_4\left(t-t_0\right)x_0+\left[1-\varepsilon \left(t-t_0\right)\right]y_0\end{array}\right)+\underset{0}{\overset{t-{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}{\int }}{e}^{\alpha z}\left(\begin{array}{c}\left[1+\varepsilon z\right]{\Gamma }_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_{2}z\\ {\Gamma }_3{\Gamma }_{4}z+\left[1-\varepsilon z\right]{\Gamma }_6\end{array}\right)dz& \mathrm{Eq}.3\text{-}25\end{array}$
Consequently,

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{c}\left[1+\varepsilon \left(t-t_0\right)\right]x_0+{\Gamma }_2\left(t-t_0\right){\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">y</figure-callout>}}_0\\ {\Gamma }_4\left(t-t_0\right)x_0+\left[1-\varepsilon \left(t-t_0\right)\right]y_0\end{array}\right)+\underset{0}{\overset{t-{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">t</figure-callout>}}_0}{\int }}{e}^{\alpha z}\left(\begin{array}{c}\left({\Gamma }_3+\left[{\mathrm{\varepsilon \Gamma }}_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\right]z\right)\\ \left({\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6+\left[{\Gamma }_3{\Gamma }_4-{\mathrm{\varepsilon \Gamma }}_6\right]z\right)\end{array}\right)dz& \mathrm{Eq}.3\text{-}26\end{array}$
Integrating by parts yields that:

$\begin{array}{cc}\int e^{\alpha z}\left[{\mathrm{\varepsilon \Gamma }}_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\right]z=\left[{\mathrm{\varepsilon \Gamma }}_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\right]\left(\frac{\alpha z-1}{{\alpha }^2}\right)e^{\alpha z}& \mathrm{Eq}.3\text{-}27a\\ \int e^{\alpha z}\left[{\Gamma }_3{\Gamma }_4-{\mathrm{\varepsilon \Gamma }}_6\right]z=\left[{\Gamma }_3{\Gamma }_4-{\mathrm{\varepsilon \Gamma }}_6\right]\left(\frac{\alpha z-1}{{\alpha }^2}\right)e^{\alpha z}& \mathrm{Eq}.3\text{-}27b\end{array}$
Substituting equations 3-27a and 3-27b into 3-26 yields the solution to Case 3 (D=0 One repeated real eigenvalue):

$\begin{array}{cc}\stackrel{\rightharpoonup }{s}\left(t\right)=e^{\alpha \left(t-t_0\right)}\left(\begin{array}{c}\left[1+\varepsilon \left(t-t_0\right)\right]x_0+{\Gamma }_2\left(t-t_0\right){\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="US09521326-20161213-D00004.png"\; state="\{\{state\}\}">y</figure-callout>}}_0\\ {\Gamma }_4\left(t-t_0\right)x_0+\left[1-\varepsilon \left(t-t_0\right)\right]y_0\end{array}\right)+\left(\begin{array}{c}\begin{array}{c}\left(\frac{{\Gamma }_3}{\alpha }+\left[{\mathrm{\varepsilon \Gamma }}_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\right]\left(\frac{\alpha \left(t-t_0\right)-1}{{\alpha }^2}\right)\right)e^{\alpha \left(t-t_0\right)}+\\ \left(-\frac{{\Gamma }_3}{\alpha }+\left(\frac{\left[{\mathrm{\varepsilon \Gamma }}_3+{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6{\Gamma }_2\right]}{{\alpha }^2}\right)\right)\end{array}\\ \begin{array}{c}\left(\frac{{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\alpha }+\left[{\Gamma }_3{\Gamma }_4-{\mathrm{\varepsilon \Gamma }}_6\right]\left(\frac{\alpha \left(t-t_0\right)-1}{{\alpha }^2}\right)\right)e^{\alpha \left(t-t_0\right)}+\\ \left(-\frac{{\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\alpha }+\left(\frac{\left[{\Gamma }_3{\Gamma }_4-{\mathrm{\varepsilon \Gamma }}_6\right]}{{\alpha }^2}\right)\right)\end{array}\end{array}\right)& \mathrm{Eq}.3\text{-}28\end{array}$

It should be noted that the constant scan coefficients employed in the approximate solution presented here depend on the coordinates X, Y, Z of the ground target point. For each solution, a set of ground coordinates is chosen to be used in computing Γ₁, . . . Γ₆. In the case where one wishes to specify a point in the focal plane instead, one may employ the non-pinhole camera projection equations to establish where that point intersects the earth's surface (usually accomplished via iteration). The resulting ground coordinates may then be employed to compute the scan coefficients used in the approximate solution.

4 Error Propagation Model and MRC

The solution to the image scan equations is propagated in time via the state transition matrix and the forcing vector. As the scan coefficients completely determine the state transition matrix and forcing vector, any errors in image motion propagation arise from errors in the scan coefficients (within the accuracy of the constant scan coefficients assumption). Thus, the error propagation model can be constructed as indicated in FIG. 7, which is a schematic diagram of an error perturbation model, according to some exemplary embodiments. FIG. 7 illustrates a state at time t₀ and two possible positions at some time t after t₀. One of the positions is defined by the nominal state transition matrix and forcing vector provided in FIG. 7, and the other position is defined by a perturbed state transition matrix and forcing vector, also provided in FIG. 7. The difference in position between the two points is the position error, which is defined by the fundamental error propagation equation, also provided in FIG. 7 and reproduced below as Eq. 4-1a.

Thus the fundamental error propagation equation is:

$\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{s}\left(t\right)=\left[e^{\left(\Lambda +\mathrm{\Delta \Lambda }\right)\left(t-t_0\right)}-e^{\Lambda \left(t-t_0\right)}\right]\stackrel{\rightharpoonup }{s}\left(t_0\right)+\hspace{1em}\left[\underset{t_0}{\overset{t}{\int }}{e}^{\left(\Lambda +\mathrm{\Delta \Lambda }\right)\left(t-\zeta \right)}\left(\stackrel{\rightharpoonup }{F}+\Delta \stackrel{\rightharpoonup }{F}\right)-e^{\Lambda \left(t-\zeta \right)}\stackrel{\rightharpoonup }{F}\right]d\zeta & \mathrm{Eq}.4\text{-}1a\\ \mathrm{\Delta \Lambda }=\left(\begin{array}{cc}{\mathrm{<figure-callout\; id="1"\; label="\Delta \Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Delta \Gamma </figure-callout>}}_1& {\mathrm{\Delta \Gamma }}_2\\ {\mathrm{\Delta \Gamma }}_4& {\mathrm{\Delta \Gamma }}_5\end{array}\right)& \mathrm{Eq}.4\text{-}1b\\ \Delta \stackrel{\rightharpoonup }{F}=\left(\begin{array}{c}{\mathrm{\Delta \Gamma }}_3\\ {\mathrm{\Delta \Gamma }}_6\end{array}\right)& \mathrm{Eq}.4\text{-}1c\end{array}$
Here Δ

(t) represents the total error in the focal plane coordinates arising from propagating forward, from the initial conditions

(t₀), in time by (t−t₀) in the presence of scan coefficient errors ΔΓ₁, . . . , ΔΓ₆. It is noted that this approach does not require the individual estimation of linear, oscillatory and random errors, as was the case with other previous applications.

According to exemplary embodiments, a first-order error propagation model is constructed by expanding the indicated state transition matrices in equation 4-1a in time, and retaining only the first order terms:

$\begin{array}{cc}{e}^{\left(\Lambda +\mathrm{\Delta \Lambda }\right)t}=\sum _{p=0}^{\infty }\frac{{\left(\Lambda +\mathrm{\Delta \Lambda }\right)}^p{t}^p}{p!}\approx I+\left(\Lambda +\mathrm{\Delta \Lambda }\right)\left(t-t_0\right)& \mathrm{Eq}.4\text{-}2a\\ e^{\Lambda \left(t-t_0\right)}=\sum _{p=0}^{\infty }\frac{{\left(\Lambda \right)}^p{\left(t-t_0\right)}^p}{p!}\approx I+\Lambda \left(t-t_0\right)& \mathrm{Eq}.4\text{-}2b\end{array}$
Substituting equations 4-2a and 4-2b into 4-1 and retaining only first order terms in the scan coefficient errors ΔΓ₁, . . . , ΔΓ₆ yields:

$\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{s}\left(t\right)=\left[\mathrm{\Delta \Lambda }\left(t-t_0\right)\right]\stackrel{\rightharpoonup }{s}\left(t_0\right)+\hspace{1em}\hspace{1em}\left[\Delta \stackrel{\rightharpoonup }{F}\left(t-t_0\right)\right]+\hspace{1em}\frac{1}{2}\left[\left(\mathrm{\Delta \Lambda }\right)\stackrel{\rightharpoonup }{F}+\Lambda \left(\Delta \stackrel{\rightharpoonup }{F}\right)\right]{\left(t-t_0\right)}^2& \mathrm{Eq}.4\text{-}3\end{array}$

It can be observed from equation 4-3 that the focal plane coordinate errors vary with propagation time. The longer the propagation time, the larger the resulting focal plane errors will be. This trend has significant practical impacts for line scanning sensors with slow line rates, or where the separation between the leading and trailing detector arrays 12, 16 is large.

Thus, according to the disclosure, an equation has been derived that predicts focal plane coordinate errors in terms of propagation time and scan coefficient errors. This equation is most useful if one has knowledge of the scan coefficient errors. The detector junction measurements can be utilized to derive the desired scan coefficients errors. Specifically, equation 4-3 will be rewritten to establish a measurement sensitivity matrix, which relates the errors in the scan coefficients to the detector junction measured coordinate differences (between the predicted and correlation-determined positions).

Equation 4-3 is expanded in terms of the individual scan coefficient errors ΔΓ₁, . . . , ΔΓ₆ and terms are rearranged to obtain:

$\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{s}\left(t\right)=H_{2\times 6}\Delta {\stackrel{\rightharpoonup }{\Gamma }}_{6\times 1}=\left(\begin{array}{cccccc}{H}_{11}& H_{12}& H_{13}& H_{14}& H_{15}& H_{16}\\ H_{21}& H_{22}& H_{23}& H_{24}& H_{25}& H_{26}\end{array}\right)\left(\begin{array}{c}{\mathrm{\Delta \Gamma }}_1\\ {\mathrm{\Delta \Gamma }}_2\\ {\mathrm{\Delta \Gamma }}_3\\ {\mathrm{\Delta \Gamma }}_4\\ {\mathrm{\Delta \Gamma }}_5\\ {\mathrm{\Delta \Gamma }}_6\end{array}\right)& \mathrm{Eq}.4\text{-}4a\\ H_{11}=\left[x_0\left(t-t_0\right)+{\Gamma }_1\frac{{\left(t-t_0\right)}^2}{2}\right]& \mathrm{Eq}.4\text{-}4b\\ H_{12}=y_0\left(t-t_0\right)+\frac{{\Gamma }_6}{2}{\left(t-t_0\right)}^2& \mathrm{Eq}.4\text{-}4c\\ H_{13}=\left(t-t_0\right)+\frac{{\Gamma }_1}{2}{\left(t-t_0\right)}^2& \mathrm{Eq}.4\text{-}4d\\ H_{14}=0& \mathrm{Eq}.4\text{-}4e\\ H_{15}=0& \mathrm{Eq}.4\text{-}4f\\ H_{16}={{\Gamma }_2\left(t-t_0\right)}^2/2& \mathrm{Eq}.4\text{-}4g\\ H_{21}=0& \mathrm{Eq}.4\text{-}4h\\ H_{22}=0& \mathrm{Eq}.4\text{-}4i\\ H_{23}=\frac{{\Gamma }_4}{2}{\left(t-t_0\right)}^2& \mathrm{Eq}.4\text{-}4j\\ H_{24}=x_0\left(t-t_0\right)+{\Gamma }_3\frac{{\left(t-t_0\right)}^2}{2}& \mathrm{Eq}.4\text{-}4k\\ H_{25}=y_0\left(t-t_0\right)+{\Gamma }_6\frac{{\left(t-t_0\right)}^2}{2}& \mathrm{Eq}.4\text{-}4l\\ H_{26}=\left(t-t_0\right)+{\Gamma }_5\frac{{\left(t-t_0\right)}^2}{2}& \mathrm{Eq}.4\text{-}4m\end{array}$

At this point one assumes the scan coefficient errors to be constant across an image line, and then employs equations 4-4a through 4-4m to execute a least-squares estimate of these scan errors given sufficient detector junction measurements Δ

_j(t) j=1 . . . m to over determine the system. It is noted that only the initial condition

_j(t₀) j=1 . . . , m and the scan coefficients ΔΓ_1j, . . . , ΔΓ_6j j=1 . . . , m vary from detector junction to junction. Thus the least squares solution becomes:

$\begin{array}{cc}\Delta {\stackrel{\rightharpoonup }{\Gamma }}_{6\times 1}={\left({\mathrm{<figure-callout\; id="6"\; label="H"\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">H</figure-callout>}}_{6\times 2m}^T{H}_{2m\times 6}\right)}^1{{\mathrm{<figure-callout\; id="6"\; label="H"\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">H</figure-callout>}}_{6\times 2m}^T\left(\begin{array}{c}\Delta {\stackrel{\rightharpoonup }{s}}_1\left(t\right)\\ ⋮\\ \Delta {\stackrel{\rightharpoonup }{s}}_m\left(t\right)\end{array}\right)}_{2m\times 1}& \mathrm{Eq}.4\text{-}5a\\ \Delta {\stackrel{\rightharpoonup }{s}}_j\left(t\right)=\left(\begin{array}{c}\Delta x_j\\ \Delta y_j\end{array}\right)& \mathrm{Eq}.4\text{-}5b\end{array}$

It is noted that equation 4-5a is only valid when the indicated matrix inversion exists. In the event that the matrix is not invertible for a particular line, this line would have its total motion error correction by interpolating the closest lines that yielded successful solutions.

FIG. 8 includes a schematic block diagram of the MRC process described in detail in the foregoing text. Given the estimate of the scan coefficient errors on a line-by-line basis, the MRC correction process of the disclosure can be implemented as depicted in FIG. 8. According to the process, total motion error is removed prior to synthetic array generation. Referring to FIG. 8, a least squares estimate of Δ

_6×1 at 140, is computed using camera attitude and position data and measured junction residuals. Using the least squares estimate, position errors are computed at 142. Scan coefficients are computed at 154, also using the camera attitude and position data. The programmed trajectory is computed at 156 using the computed scan coefficients. The position errors computed at 142 are subtracted at 164 from the programmed trajectory computed in 156 to generate the synthetic array.

5 Covariance Error Model

According to exemplary embodiments, an error model has been developed that relates a given set of scan coefficient errors to the resulting focal plane coordinate errors. Considering the scan coefficient errors to be random variables with a known joint probability distribution, expressions for the mean and covariance of the resulting focal plane coordinate errors are obtained. These expressions are referred to herein as the covariance model.

Define the expectation operation

with respect to the joint probability distribution p(ΔΓ₁, . . . , ΔΓ₆) via:

ƒ(ΔΓ₁, . . . ,ΔΓ₆)

=∫∫∫ƒ(ΔΓ₁, . . . ,ΔΓ₆)p(ΔΓ₁, . . . ,ΔΓ₆)d(ΔΓ₁) . . . d(ΔΓ₆)  Eq. 5-1

This expectation operator will now be employed to compute the desired mean and covariance of the focal plane coordinate errors. Hence, applying the expectation operator to equation 4-4a yields the mean error

Δ

(t)

to be:

Δ

(t)

=H _2×6

Δ

_6×1

  Eq. 5-2
Likewise one can compute the covariance

of the focal plane coordinate errors via:

≡

(Δ

(t)−

Δ

(t)

)(Δ

(t)−

Δ

(t)

)^T

  Eq. 5-3a

=H _2×6

Δ

_6×1Δ

^T _1×6

H ^T _6×2 −H _2×6

Δ

_6×1

Δ

^T _1×6

H ^T _6×2  Eq. 5-3b

=H _2×6 P _ΔΓ H ^T _6×2  Eq. 5-3c

Here P_ΔΓ is the 6×6 covariance matrix of the scan coefficient errors. A representative scan coefficient covariance is computed by employing a sufficient number of the line-to-line individual scan coefficient error vectors from the least squares process via a “moving window” approach.

While the covariance model does not appear directly in the MRC process equations, it is important to being able to execute “on the-fly” statistical consistence checks and detection junction measurement outlier rejection.

6 Platform Error Inversion Process

The estimated errors in the scan coefficients are related on a line-by-line basis to knowledge errors in the fifteen parameters of the platform support data (sensor position (3) and velocity (3), attitude angles (3) and rates (3), ground point location (3)). It is assumed that the underlying platform parameter knowledge errors are constant across an image line. One then employs all the valid junction measurements for a particular line and then inverts the first-order error model via an over determined least squares process as follows. Specifically, at each detector junction one has:

$\begin{array}{cc}\Delta \stackrel{\rightharpoonup }{s}\left(t\right)=H_{2\times 6}\Delta {\stackrel{\rightharpoonup }{\Gamma }}_{6\times 1}==H_{2\times 6}{\mathrm{<figure-callout\; id="6"\; label="D"\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">D</figure-callout>}}_{6\times 15}\Delta {\stackrel{\rightharpoonup }{\Psi }}_{15\times 1}& \mathrm{Eq}.6\text{-}1a\\ {\mathrm{<figure-callout\; id="6"\; label="D"\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">D</figure-callout>}}_{6\times 15}=\left(\begin{array}{ccccccccccccccc}\frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial \omega }& \frac{\partial {\Gamma }_1}{\partial \varphi }& \frac{\partial {\Gamma }_1}{\partial \kappa }& \frac{\partial {\Gamma }_1}{\partial \stackrel{.}{\omega }}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial \stackrel{.}{\varphi }}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial \stackrel{.}{\kappa }}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial X}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial Y}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial Z}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial X_c}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial Y_c}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial Z_c}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial {\stackrel{.}{X}}_c}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial {\stackrel{.}{Y}}_c}& \frac{\partial {\mathrm{<figure-callout\; id="1"\; label="\Gamma "\; filenames="US09521326-20161213-D00002.png,US09521326-20161213-D00006.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_1}{\partial {\stackrel{.}{Z}}_c}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\Gamma }_6}{\partial \omega }& \frac{\partial {\Gamma }_6}{\partial \varphi }& \frac{\partial {\Gamma }_6}{\partial \kappa }& \frac{\partial {\Gamma }_6}{\partial \stackrel{.}{\omega }}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial \stackrel{.}{\varphi }}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial \stackrel{.}{\kappa }}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial X}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial Y}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial Z}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial X_c}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial Y_c}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial Z_c}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial {\stackrel{.}{X}}_c}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial {\stackrel{.}{Y}}_c}& \frac{\partial {\mathrm{<figure-callout\; id="6"\; label="\Gamma "\; filenames="US09521326-20161213-D00008.png"\; state="\{\{state\}\}">\Gamma </figure-callout>}}_6}{\partial {\stackrel{.}{Z}}_c}\end{array}\right)& \mathrm{Eq}.6\text{-}1b\\ \Delta {\stackrel{\rightharpoonup }{\Psi }}_{15\times 1}={\left(\begin{array}{c}\mathrm{\Delta \omega }\\ ⋮\\ \Delta {\stackrel{.}{Z}}_c\end{array}\right)}_{15\times 1}& \mathrm{Eq}.6\text{-}1c\end{array}$
One can employ all the valid junction measurement for a particular line to write:

$\begin{array}{cc}\left(\begin{array}{c}\Delta {\stackrel{\rightharpoonup }{s}}_1\\ \Delta {\stackrel{\rightharpoonup }{s}}_2\\ ⋮\\ \Delta {\stackrel{\rightharpoonup }{s}}_M\end{array}\right)={\left(\begin{array}{c}{\hspace{0.17em}}_{1}H{\hspace{0.17em}}_{1}D\\ {\hspace{0.17em}}_{2}H{\hspace{0.17em}}_{2}D\\ ⋮\\ {\hspace{0.17em}}_{M}H{\hspace{0.17em}}_{M}D\end{array}\right)}_{2M\times 15}{\left(\begin{array}{c}\mathrm{\Delta \omega }\\ ⋮\\ \Delta {\stackrel{.}{Z}}_c\end{array}\right)}_{15\times 1}& \mathrm{Eq}.6\text{-}2a\\ \Delta {\stackrel{\rightharpoonup }{S}}_{2M\times 1}={\Theta }_{2M\times 15}{\mathrm{\Delta \Psi }}_{15\times 1}& \mathrm{Eq}.6\text{-}2b\end{array}$
The resulting least squares solution over m valid junctions is:

$\begin{array}{cc}\Delta {\stackrel{\rightharpoonup }{\Psi }}_{15\times 1}={\left\{{\left({\Theta }^T\Theta \right)}^{-1}{\Theta }^T\right\}}_{15\times 2M}\Delta {\stackrel{\rightharpoonup }{S}}_{2m\times 1}& \mathrm{Eq}.6\text{-}3a\\ {\Theta }_{2M\times 15}={\left(\begin{array}{c}{\hspace{0.17em}}_{1}H{\hspace{0.17em}}_{1}D\\ {\hspace{0.17em}}_{2}H{\hspace{0.17em}}_{2}D\\ ⋮\\ {\hspace{0.17em}}_{M}H{\hspace{0.17em}}_{M}D\end{array}\right)}_{2M\times 15}& \mathrm{Eq}.6\text{-}3b\end{array}$

Hence, one can estimate the platform parameter knowledge errors on a line-by-line basis and employ these estimates along with the a priori knowledge estimates to improve the knowledge of the platform parameters throughout an image. Indeed, one can iteratively feedback these improved platform position errors to iterate the MRC correction process, thus improving overall algorithm robustness and accuracy.

7 MRC Summary Description

According to the exemplary embodiments, the MRC approach described herein in detail includes certain assumptions upon which the approach is constructed. These assumptions are:

Assumption 1: The non-pinhole camera model variations can be accurately modeled via nth-order bivariate polynomials obtained from curve fitting the non-pinhole disparities obtained from a sufficient number of ray traces. This disclosure describes quadratic bivariate polynomials for mathematical brevity and proof of concept, as an exemplary illustration only.

Assumption 2: The nonlinear scan equations can be approximated about a spatial neighborhood of any chosen initial projection point and corresponding time increment via a first order Taylor series expansion. This results in a set of linear first order differential equations. The extent of the approximation is sufficient to accurately project from the leading to the trailing array.

Assumption 3: The resulting system of coupled first-order differential equations may be integrated over the associated time increment by holding the scan coefficients constant.

Assumption 4: The error propagation characteristics of the linearized scan equations can be determined via a first-order perturbation of the state transition matrix and the associated forcing vector via additive scan coefficient errors.

Assumption 5: While the scan coefficients obtained from the linearized scan equations will vary at each detector junction, it is assumed that a single set of scan coefficient errors can be obtained across all the detector junctions via a least squares process that is consistent with the measured junction coordinate residuals.

Assumption 6: The knowledge errors in the platform support parameters can be assumed constant across an individual scan line.

FIG. 9 includes a schematic logical flow diagram of the line-by-line MRC computational process, according to some exemplary embodiments. Referring to FIG. 9, it is noted that the entire process is iterative, line-by-line. Specifically, initial platform support parameters are used in step 202 to project an initial point on the leading array to the trailing array at each detector junction across a line. Next, in step 204, correlation at each detector junction is performed to establish actual projected point location. Next, at step 206, at each junction, a difference between projected and measured focal plane coordinates is computed. Next, at step 208, the coordinate differences computed at 206 are used to estimate, for example, by least-squares estimation, the scan coefficient errors for a line. Next, in step 212, the coordinate differences computed at 206 are used to estimate, for example, by least-squares estimation, errors in platform support parameters. If another iteration is required, in step 210, the estimated errors in platform support parameters are used to correct the errors, and flow returns to step 202. In contrast, after step 212, if iteration is complete, in step 214, the estimated scan errors are used to correct the projection from each position along the synthetic array to the actual detector locations. It is noted that the iteration branch depicted in FIG. 9 is optional based on the accuracy of the estimated errors in the support parameters. It should also be noted that in some applications the rate at which the platform position, velocity, attitude angles and rates are measured by on-platform equipment may not support a line-by-line implementation. In these cases one would compute new corrections every m lines, where m is chosen to be consistent with the interpolation accuracy supported by the time sampling granularity.

Various embodiments of the above-described systems and methods may be implemented in digital electronic circuitry, in computer hardware, firmware, and/or software. The implementation can be as a computer program product (i.e., a computer program tangibly embodied in an information carrier). The implementation can, for example, be in a machine-readable storage device and/or in a propagated signal, for execution by, or to control the operation of, data processing apparatus. The implementation can, for example, be a programmable processor, a computer, and/or multiple computers.

A computer program can be written in any form of programming language, including compiled and/or interpreted languages, and the computer program can be deployed in any form, including as a stand-alone program or as a subroutine, element, and/or other unit suitable for use in a computing environment. A computer program can be deployed to be executed on one computer or on multiple computers at one site.

Method steps can be performed by one or more programmable processors and/or controllers executing a computer program to perform functions of the invention by operating on input data and generating output. Method steps can also be performed by, and an apparatus can be implemented as, special purpose logic circuitry. The circuitry can, for example, be a FPGA (field programmable gate array) and/or an ASIC (application-specific integrated circuit). Modules, subroutines, and software agents can refer to portions of the computer program, the processor, the special circuitry, software, and/or hardware, e.g., a controller such as a microcontroller, that implements that functionality.

Processors suitable for the execution of a computer program include, by way of example, both general and special purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor receives instructions and data from a read-only memory or a random access memory or both. The essential elements of a computer are a processor for executing instructions and one or more memory devices for storing instructions and data. Generally, a computer can be operatively coupled to receive data from and/or transfer data to one or more mass storage devices for storing data, e.g., magnetic, magneto-optical disks, or optical disks.

Data transmission and instructions can also occur over a communications network. Information carriers suitable for embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices. The information carriers can, for example, be EPROM, EEPROM, flash memory devices, magnetic disks, internal hard disks, removable disks, magneto-optical disks, CD-ROM, and/or DVD-ROM disks. The processor and the memory can be supplemented by and/or incorporated in special purpose logic circuitry.

To provide for interaction with a user, the above described techniques can be implemented on a computer having a display device. The display device can, for example, be a cathode ray tube (CRT) and/or a liquid crystal display (LCD) monitor. The interaction with a user can, for example, be a display of information to the user and a keyboard and a pointing device, e.g., a mouse or a trackball, by which the user can provide input to the computer, e.g., interact with a user interface element. Other kinds of devices can be used to provide for interaction with a user. Other devices can, for example, be feedback provided to the user in any form of sensory feedback, e.g., visual feedback, auditory feedback, or tactile feedback. Input from the user can, for example, be received in any form, including acoustic, speech, and/or tactile input.

The above described techniques can be implemented in a distributed computing system that includes a back-end component. The back-end component can, for example, be a data server, a middleware component, and/or an application server. The above described techniques can be implemented in a distributing computing system that includes a front-end component. The front-end component can, for example, be a client computer having a graphical user interface, a Web browser through which a user can interact with an example implementation, and/or other graphical user interfaces for a transmitting device. The components of the system can be interconnected by any form or medium of digital data communication, e.g., a communication network. Examples of communication networks include a local area network (LAN), a wide area network (WAN), the Internet, wired networks, and/or wireless networks.

The system can include clients and servers. A client and a server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other.

Packet-based networks can include, for example, the Internet, a carrier internet protocol (IP) network, e.g., local area network (LAN), wide area network (WAN), campus area network (CAN), metropolitan area network (MAN), home area network (HAN), a private IP network, an IP private branch exchange (IPBX), a wireless network, e.g., radio access network (RAN), 802.11 network, 802.16 network, general packet radio service (GPRS) network, HiperLAN, and/or other packet-based networks. Circuit-based networks can include, for example, the public switched telephone network (PSTN), a private branch exchange (PBX), a wireless network, e.g., RAN, Bluetooth, code-division multiple access (CDMA) network, time division multiple access (TDMA) network, global system for mobile communications (GSM) network, and/or other circuit-based networks.

The computing system can also include one or more computing devices. A computing device can include, for example, a computer, a computer with a browser device, a telephone, an IP phone, a mobile device, e.g., cellular phone, personal digital assistant (PDA) device, laptop computer, electronic mail device, and/or other communication devices. The browser device includes, for example, a computer, e.g., desktop computer, laptop computer, with a World Wide Web browser, e.g., Microsoft® Internet Explorer® available from Microsoft Corporation, Mozilla® Firefox available from Mozilla Corporation. The mobile computing device includes, for example, a Blackberry®, iPAD®, iPhone® or other smartphone device.

Whereas many alterations and modifications of the disclosure will no doubt become apparent to a person of ordinary skill in the art after having read the foregoing description, it is to be understood that the particular embodiments shown and described by way of illustration are in no way intended to be considered limiting. Further, the subject matter has been described with reference to particular embodiments, but variations within the spirit and scope of the disclosure will occur to those skilled in the art. It is noted that the foregoing examples have been provided merely for the purpose of explanation and are in no way to be construed as limiting of the present disclosure.

While the present disclosure has been described with reference to exemplary embodiments, it is understood that the words that have been used herein, are words of description and illustration, rather than words of limitation. Changes may be made, within the purview of the appended claims, as presently stated and as amended, without departing from the scope and spirit of the present disclosure in its aspects.

Although the present disclosure has been described herein with reference to particular means, materials and embodiments, the present disclosure is not intended to be limited to the particulars disclosed herein; rather, the present disclosure extends to all functionally equivalent structures, methods and uses, such as are within the scope of the appended claims.

Whereas many alterations and modifications of the disclosure will no doubt become apparent to a person of ordinary skill in the art after having read the foregoing description, it is to be understood that the particular embodiments shown and described by way of illustration are in no way intended to be considered limiting. Further, the subject matter has been described with reference to particular embodiments, but variations within the spirit and scope of the disclosure will occur to those skilled in the art. It is noted that the foregoing examples have been provided merely for the purpose of explanation and are in no way to be construed as limiting of the present disclosure.

Claims (13)

The invention claimed is:

1. A method of misregistration correction in a line scanning imaging system, comprising:

providing the line scanning imaging system with a forward array of optical detectors and a trailing array of optical detectors for receiving optical signals from a scene being scanned and converting the optical signals into electrical signals, the forward array of optical detectors receiving the optical signals at a different time than the trailing array of optical detectors receiving the optical signals, the forward array of optical detectors and the trailing array of optical detectors defining a plurality of detector junction overlap regions;

receiving the electrical signals at a processor;

processing the received electrical signals with the processor, the processing comprising the steps of:

generating a model of scan motion over a focal plane of the imaging system, using a coupled system of non-linear scan equations with constant coefficients;

estimating programmed motion positions across the plurality of detector junction overlap regions via a state transition matrix solution to the non-linear scan equations;

at each detector junction overlap region, measuring actual motion positions via image correlation of overlapping detectors;

generating differences between the actual motion positions and the estimated programmed motion positions;

estimating updates to the constant coefficients based on the generated differences;

generating corrections from the estimated updates to remove unwanted motion; and

applying the updates to the constant coefficients; wherein

the steps of generating a model of scan motion, estimating programmed motion positions, measuring actual motion positions, generating differences, estimating updates, generating corrections and applying the updates are performed by means of a non-transitory computer-usable medium having computer-readable code that, when executed by the processor, cause the processor to perform the steps.

2. The method of claim 1, wherein estimating updates is performed using least-squares estimation.

3. The method of claim 1, wherein the image correlation comprises normalized cross-correlation.

4. The method of claim 1, wherein the image correlation comprises lag product cross-correlation.

5. The method of claim 1, wherein the image correlation comprises least squares cross-correlation.

6. The method of claim 1, wherein the model of scan motion is generated over a predetermined time interval.

7. The method of claim 1, wherein the non-linear scan equations comprise a set of differential equations with constant coefficients.

8. The method of claim 7, wherein the differential equations are first-order differential equations.

9. The method of claim 1, further comprising computing the non-linear scan equation coefficients using sensor platform parameters.

10. The method of claim 9, wherein the sensor platform parameters comprise at least one of sensor position, velocity, attitude angles and rates.

11. The method of claim 10, further comprising approximating errors in the non-linear scan equation coefficients from measurement errors in the platform parameters.

12. The method of claim 11, further comprising generating a covariance matrix of focal plane coordinate errors.

13. The method of claim 10, further comprising estimating platform parameter errors from the differences between the actual motion positions and the estimated programmed motion positions.

Images