The ability to view an object is limited by turbulence existing in the line of sight between the viewer and the object. Such turbulence may exist in the atmosphere, including a body of air, or in other intervening mediums, such as fluids. As is known, atmospheric turbulence is primarily caused by non-uniform heating and wind, and includes both naturally occurring and man-made turbulence, such as the atmospheric turbulence in the wake of a jet engine.
Passive incoherent imaging techniques such as phase diversity and wavelength diversity have long been used to compensate for the negative effects of atmospheric turbulence on conventional imaging systems. Both of these diversity schemes traditionally use iterative 2-D Fourier transforms in a sequential error-reduction methodology that is typically slow and relegated to post-processing applications.
Atmospheric coherence length, r0, is a measure of the effect of atmospheric turbulence on optical propagation. A larger value implies weaker turbulence. For optical imaging systems having apertures that are significantly larger than the atmospheric coherence length, r0 (Fried parameter), the turbulent atmosphere usually is the most significant contributor to the imaging system's loss of spatial resolution. Atmospheric turbulence introduces aberrations that can significantly degrade the performance of an optical imaging system. The degradation in imaging performance depends on a variety of factors, such as the operating wavelength(s), the relative size of the imaging system's entrance pupil diameter as compared to the atmospheric coherence length, r0, the angular instantaneous field of view (IFOV) as compared to the isoplanatic angle, θ0, sampling effects, signal-to-noise issues, and system effects such as imperfect optics, fixed pattern noise and platform jitter.
Given a circular aperture, the upper bound on achievable resolution is the so-called diffraction limit,
                                          Δ            ⁢                                                  ⁢                          x              diff                                =                      1.22            ⁢                                          λ                _                            D                        ⁢            z                          ,                            (        1        )            where λ is the center wavelength of the illuminating light, D is the diameter of the entrance pupil of the imaging system (for example, the diameter of a telescope, or camera lens), and z is the distance between the imaging system's entrance pupil and the object/target to be imaged. For example, an object located 2 km from a camera with a 4½ inch lens would have a diffraction-limited resolution of 1.059 cm at a center wavelength of 500 nm. This means that object features less than this size will not be resolved and object features larger than this size will be resolved if diffraction-limiting imaging conditions are attained.
Unfortunately, the atmosphere severely degrades the ability to achieve the diffraction-limited resolution indicated in equation (1). Instead, the conventionally attainable maximum resolution of an imaging system that is looking through atmospheric turbulence is given by,
                                          Δ            ⁢                                                  ⁢                          x              atm                                =                      1.22            ⁢                                          λ                _                                            r                0                                      ⁢            z                          ,                            (        2        )            where the r0 expression in the denominator of equation (2) is the atmospheric coherence length, or Fried parameter. For imaging over horizontal paths, recent experimental values for r0 have ranged from 1 to 4 cm leading to a “loss” of resolution of 11.52 (for an r0 of 1 cm) to 2.88 (for an r0 of 4 cm) with respect to the diffraction-limited resolution given in equation (1). Stated another way, by compensating for the effects of atmospheric turbulence, a maximum increase in resolution between 2.88 and 11.52 can be expected for the center wavelength and entrance-pupil diameter specified above. These potential effects are even more dramatic when a telescopic imaging system is used. For instance, if an eight-inch telescope is attached to the camera or video camera at optical wavelengths, the maximum increase in spatial resolution for the same values of r0 above jumps to 20.32 (for an r0 of 1 cm) or 5.08 (for an r0 of 4 cm).
As can be seen by dividing equation (2) by equation (1), the potential increase in spatial resolution of an imaging system with entrance-pupil diameter D is given by,
                    R        =                                            Δ              ⁢                                                          ⁢                              x                atm                                                    Δ              ⁢                                                          ⁢                              x                diff                                              =                                    D                              r                0                                      .                                              (        3        )            
With knowledge of the diameter of the imaging system's entrance pupil and the value of r0 for an illuminating wavelength, equation (3) can be used to determine the maximum achievable increase in resolution, neglecting system effects and assuming full compensation of atmospheric turbulence.
The value of r0 at a given wavelength can be scaled to another wavelength by,
                                          r            0                          λ              2                                =                                    r              0                              λ                1                                      ·                                          (                                                      λ                    2                                                        λ                    1                                                  )                                            6                5                                                    ,                            (        4        )            where r0λ1 is the value of r0 at wavelength λ1 and r0λ2 is the value of r0 at the wavelength of interest λ2. For example, if the value of r0 at the illuminating wavelength of 500 nm is 1 cm, then the value of r0 at 1.06 μm is about 2.46 cm. Once the new value of r0 is known, equation (3) can be used to determine the upper bound on the resolution increase of the optical system under consideration. For color images, the resolution increase is different at each wavelength and can be estimated by using the red (R), green (G), and blue (B) wavelengths fundamental to a RGB color scheme.
Equation (4) can be obtained from
                                          r            0                    =                      0.185            ⁢                                          (                                                      λ                    2                                                                              ∫                      0                      z                                        ⁢                                                                                            c                          n                          2                                                ⁡                                                  (                          ξ                          )                                                                    ⁢                                                                                          ⁢                                              ⅆ                        ξ                                                                                            )                                            3                5                                                    ,                            (        5        )            where the term in the denominator is the integral of the atmospheric structure constant along the optical path [0, z]. Very often in imaging applications, the Fried parameter is parameterized. In this case, the integral in expression (5) can be estimated from the parameter r0.
In order to attain an increase in resolution approaching the limits given by equation (3), the effects of atmospheric turbulence must be de-convolved from the aberrated image. FIG. 1 shows the simulated effects of atmospheric turbulence for a target that is simulated to be 1.3 km away from the optical imaging system. The center wavelength was 550 nm and the distance between the simulated object and the 8-inch telescopic imaging system was 1.3 km. The Fried coherence length was 2 cm. The imaging system itself was assumed to be noise-free to isolate and demonstrate the atmospheric turbulence effects on the resulting image.
Three basic ways to estimate and remove the effects of atmospheric turbulence are adaptive optics systems, post-processing atmospheric turbulence compensation systems and hybrids of the two.
Adaptive optics systems are hardware-based systems that can correct atmospheric-turbulence effects in real-time, for example, at rates faster than 30 Hz. Adaptive optics systems are often cumbersome, require extensive hardware, alignment, and expertise, and are often expensive. They are also predominantly designed for fixed sites and are not typically man-portable, rugged or covert.
Post-processing atmospheric turbulence compensation systems are largely implemented in software (software dominant) but traditionally are very slow (not real-time), and can involve additional hardware. Hybrid methods are a cross between the two and generally relax some of the hardware processing requirements and then correct for the loss in performance using post-processing methods.
Many of the existing software-dominant atmospheric turbulence compensation (deconvolution) approaches such as wavelength diversity, phase diversity, multi-frame blind deconvolution and other post-processing approaches have been predominantly focused on phase-only corrections (near-field turbulence approximation) and have been iterative and slow. There are some exceptions. In “Phase-diversity correction of turbulence-induced space-variant blur”, R. Paxman and others (“Paxman”) described use of phase diversity for reconstructing the unaberrated object brightness through distributed turbulence. (Optics Letters, Vol. 19, No. 16, pp. 1231-1233, 1994). In “Automated Video Enhancement from a Stream of Atomospherically-distorted Image, the Lucky-Region Fusion Approach”, M. Aubailly, M. Vorontsov and others (“Aubailly”) described methods for fusing segments of “lucky” image frames to de-blur images obtained over horizontal optical paths. (Proc. SPIE 7463, 74630C, 2009). In “Special-purpose hardware for real-time compensation of atmospheric effects in long-range imaging” and “Reconfigurable device for enhancement of long-range imagery”, F. Ortiz, C. J. Carrano and others (“Carrano”) generalized the traditional speckle-imaging technique that is and has been used to great effect in astronomy and other near-field turbulence application areas. (Atmospheric Optical Modeling, Measurement, and Simulation II, Proc. SPIE Vol. 6303, 2006, and Airborne Intelligence, Surveillance, Reconnaissance (ISR) Systems and Applications IV, Proc. SPIE Vol. 6546, 2007, respectively).
The approach of Carrano improves the speed of the traditional speckle imaging-based post-processing methods. Carrano's approach parameterized the Korff transfer function in terms of the atmospheric coherence length, r0 to estimate the object magnitude and used the bispectrum technique to estimate the object phase without requiring the presence of a reference source such as a star or laser guide-star. The approach was applied to a horizontal path turbulence scenario with observable improvement in clarity for optical path lengths of 0.5 km to 10 km. Processing was done on a 256 by 256 pixel by 100 image frame data cube using a 1.7 GHz Pentium IV processor. It originally took about 10 seconds to process this data cube, which represented an isoplanatic patch from a 1280 by 1024 image. To recover the full image, 20 similar data cubes must be processed resulting in 200 seconds of processing time for their non-optimized algorithm. Carrano described continuing to improve processing times to approximately 1 second using faster processing equipment.
Nonetheless, there are significant limitations to the Carrano approach. For instance, the Carrano approach could require hundreds of frames of images or more to provide the parameters needed for it to work. The frames require additional processing power and time. The method further requires that the frames capture the same image—necessitating that the subject being viewed is not moving. Limitations such as these make the approach problematic. For example, consider atmospheric turbulence compensating applications that require real-time image processing of uncooperative targets, such as covert surveillance of people for military or law enforcement purposes.
The following is an overview of the optical-systems model used in general-purpose incoherent optical imaging systems.
For many incoherent imaging applications, a linear, shift-invariant imaging model is appropriate, with the optical systems model given by,i({right arrow over (x)})=o({right arrow over (x)})*|hi({right arrow over (x)})|2,   (6)where o({right arrow over (x)}) is the 2-D pristine object brightness function, |hi({right arrow over (x)})|2 is the imaging system's point spread function (PSF), i({right arrow over (x)}) is the atmospherically degraded (e.g. “blurry”) image due to atmosphere and optical imaging system effects, and {right arrow over (x)} is a 2-D position vector in the image plane. The asterisk represents 2-D spatial convolution.
By taking the 2-D Fourier transform of both sides of equation (6), the frequency-space equivalent of equation (6) is given byI({right arrow over (f)})=O({right arrow over (f)})H({right arrow over (f)}),   (7)where I({right arrow over (f)}) is the image spectrum, O({right arrow over (f)}) is the object spectrum, H({right arrow over (f)}) is the optical transfer function (OTF) and {right arrow over (f)} is a 2-D spatial frequency variable. Equations (6) and (7) apply at each spatial coordinate {right arrow over (x)} and at each spatial frequency {right arrow over (f)}.
The PSF can be related to the optical transfer function by,
                                          H            ⁡                          (                              f                →                            )                                =                                    𝔍              ⁢                              ⌊                                                                                                                        h                        i                                            ⁡                                              (                                                  x                          →                                                )                                                                                                  2                                ⌋                                                                    𝔍                ⁡                                  [                                                                                                                                    h                          i                                                ⁡                                                  (                                                      x                            →                                                    )                                                                                                            2                                    ]                                                                              f                  →                                =                0                                                    ,                            (        8        )            
where the symbol, ℑ[●], denotes taking the 2-D Fourier transform of the expression inside the brackets. The optical transfer function is seen to be the 2-D Fourier transform of the PSF and then normalized so that the maximum value of the OTF is 1.
The relationship between the optical system's impulse response hi({right arrow over (x)}) and the generalized pupil function (GPF) is given by,hi({right arrow over (x)})=ℑ−1[W({right arrow over (x)})],   (9)where ℑ−1[●] represents taking the 2-D inverse Fourier transform of the expression inside the brackets, andW({right arrow over (x)})=A({right arrow over (x)})ejφ({right arrow over (x)}),   (10)is the generalized pupil function (GPF). The function A({right arrow over (x)}) is an aperture function that has a value of 1 inside the clear aperture of the telescope and 0 outside of it. The function φ({right arrow over (x)}) is the atmospherically induced phase aberration realized at spatial coordinate {right arrow over (x)} in the entrance pupil of the imaging system. In the near-field turbulence approximation model, such as if one is looking up through the atmosphere into space from a ground-based telescope, the amplitude variations are assumed negligible within an isoplanatic patch, and are set to 1. It can be seen that there are two components that make up equation (10)—one due to the imaging system's aperture, and the other due to the effects of atmospheric turbulence.
Traditional “diversity-based” post-processing atmospheric turbulence compensation methods involve inserting a known diversity into the imaging system by some artifice and then simultaneously capturing the original image and the diversity image. The OTFs of the aberrated image and the diversity image are related to each other and a suitable error metric is used to select an error minimized OTF. Ideally, the error minizimed OTF is the OTF that produces the lowest local and global error. Depending on the approach taken, in relating the OTF to the diversity OTF, often the resulting error metric will be a function of the aberrated image spectrum, the diversity image spectrum (both measured) and analytical expressions of the OTF and diversity OTF that are both functions of entrance pupil phase estimates or phase difference estimates. For instance, one phase-diversity post-processing atmospheric turbulence compensation method uses a predetermined additive phase term in the expression of the phase diversity generalized pupil function,Wpd({right arrow over (x)})=A({right arrow over (x)})ej(φ({right arrow over (x)})+φpd({right arrow over (x)})),   (11)where the subscript pd denotes the phase diversity method was used. The expression Wpd({right arrow over (x)}) is the phase diversity generalized pupil function and, as shown, has a predetermined phase diversity φpd({right arrow over (x)}) added to the unknown atmospheric turbulence phase φ({right arrow over (x)}) at every entrance pupil spatial coordinate {right arrow over (x)}. Often a predetermined quadratic phase factor can be introduced in the phase diversity image by slightly defocusing the diversity image.
Another diversity method is one by which an image is captured simultaneously at two different narrow-band wavelengths centered at λ1 and λ2. The wavelength diversity generalized pupil function is then given by,
                                          W            wd                    ⁡                      (                          x              →                        )                          =                              A            ⁡                          (                              x                →                            )                                ⁢                                    ⅇ                              j                ⁡                                  (                                                                                    λ                        1                                                                    λ                        2                                                              ⁢                                          ϕ                      ⁡                                              (                                                  x                          →                                                )                                                                              )                                                      .                                              (        12        )            In the traditional diversity-based atmospheric turbulence compensation methods, a diversity OTF is generated by using equations (8) through (10) in reverse order and substituting the appropriate diversity generalized pupil function from equation (11) or (12) depending on which diversity method one is using for equation (10).
A common error metric is then applied at every point in the image spectrum as a way to determine when the OTF estimate is accurate enough. One such error metric is described by Gonsalves and others (“Gonsalves”) in “Phase Retrieval from Modulus Data” (J. Opt. Soc. Am., Vol. 66, No. 9, September 1976) and “Wavefront Sensing by Phase Retrieval” (SPIE. Vol. 207, Applications of Digital Image Processing III, 1979),
                                          E            ⁡                          (                              f                →                            )                                =                                                                                                                                I                      ⁡                                              (                                                  f                          →                                                )                                                              ⁢                                                                                            H                          ^                                                d                                            ⁡                                              (                                                  f                          →                                                )                                                                              -                                                                                    I                        d                                            ⁡                                              (                                                  f                          →                                                )                                                              ⁢                                                                  H                        ^                                            ⁡                                              (                                                  f                          →                                                )                                                                                                                        2                                                                                                                                    H                      ^                                        ⁡                                          (                                              f                        →                                            )                                                                                        2                            +                                                                                                                                    H                        ^                                            d                                        ⁡                                          (                                              f                        →                                            )                                                                                        2                                                    ,                            (        13        )            The carat symbol ^ above the OTF and diversity OTF indicate that these quantities are estimated. This and other known error equations need modification according to the teachings of this disclosure to work generally using the wavelength diversity methodology when large wavelength separations are employed between the primary image and the simultaneously captured diversity image.
Object brightness (e.g., Watts/cm2) has posed a problem for prior wavelength diversity compensation methods. For example, the error metric in equation (13) assumes that the object brightness and the diversity object brightness are the same—hence they divide out in the derivation process of this known equation. This is true for the phase diversity approach but not always true for the wavelength diversity approach. In wavelength diversity, the object brightness changes as a function of wavelength. As an example, looking at the red, green and blue components of an RGB image, a given object's brightness can be significantly different at each of the red, green and blue color components. This means that equation (13) does not remain valid except for cases where the different wavelength image (e.g. diversity image) is close to the original image. Equation (13) cannot be used for the wavelength diversity method in the general case where wavelength separations between the image and diversity image may be large. For example, this prior art technique is not expected to be reliable for separations of greater than approximately 10% or 40 nanometers.
Others have noted this shortcoming. For example, in “Parallel multiframe blind deconvolution using wavelength diversity”, H. R. Ingleby and others (“Ingleby”) described an alternate error metric that supposedly mitigated the object brightness issue, proposing the following error metric (Proceedings of the SPIE, Volume 5562, pp. 58-64, October 2004):
                                                        E              m                        ⁡                          (                              a                kn                            )                                =                                    ∑              v                        ⁢                          [                                                                    ∑                                          n                      =                      1                                        N                                    ⁢                                                                                                          I                        mn                                                                                    2                                                  -                                                                                                                                                    ∑                                                      n                            =                            1                                                    N                                                ⁢                                                                              H                            mn                            *                                                    ⁢                                                      I                            mn                                                                                                                                      2                                                                              ∑                                              n                        =                        1                                            N                                        ⁢                                                                                                                    H                          mn                                                                                            2                                                                                  ]                                      ⁢                                  ⁢                  and          ,                                    (        16        )                                                      O            m                    =                                    ∑                              n                =                1                            N                        ⁢                                          H                mn                *                            ⁢                                                I                  mn                                /                                                      ∑                                          n                      =                      1                                        N                                    ⁢                                                                                                          H                        mn                                                                                    2                                                                                      ,                            (        17        )            where the “m” index sums over wavelengths, the “n” index sums over the number of frames (e.g. snapshots taken at different time intervals), the “v” index sums over spatial frequency locations in the entrance pupil of the imaging system, and the “k” index runs over the Zernike modes. The “a” coefficient is the estimated weight on the “kth” Zernike mode of the nth frame. The object estimate for the mth wavelength component is then given by equation (17), which follows from application of linear systems theory. Equation (16) has the unfortunate characteristic, however, that it goes to zero for any arbitrary choice of Hm when the number of frames is identically one (e.g. for a single snapshot). Therefore, equation (16) is unsuitable for WD systems that only have one frame of data (e.g. a single snapshot).
In traditional diversity-based atmospheric turbulence compensation methods, the process for estimating the OTF (and also the diversity OTF by analogy) includes:                1. Using a suitable basis set like the Zernike polynomials initially to generate an entrance pupil plane-phase estimate. This is typically done by initially guessing the phase—for example, assuming that all phase values are zero;        2. Forming the generalized pupil function with this entrance pupil phase “guess” using equation (10). Equations (11) or (12) are typically selected for estimating the diversity OTF;        3. Zero-packing the GPF for sampling reasons in preparation of generating an OTF estimate;        4. Forming the impulse response in accordance with equation (9);        5. Determining the PSF estimate by point-wise taking the magnitude squared of the result of step 4 above;        6. Using equation (8) to form the OTF estimate;        7. After forming both the OTF and diversity OTF estimate, applying an error metric such as the Gonsalvez error metric given in equation (13) and optionally one or more constraints known to one of ordinary skill in the art to determine the instantaneous error at each spatial location in the OTF. The constraints typically involve enforcing some physical aspect of the imaging problem such as positivity of the object brightness, applying conservation of energy principles, point spread function and image, and an experiment specific constraint such as including a pre-defined entrance pupil aperture size.        8. Summing the errors to determine the total summed error due to the initial entrance pupil plane phase estimation error;        9. Changing the weights on the Zernike polynomials in a methodical manner to come up with a new entrance pupil plane phase estimate;        10. Repeating steps 2 through 8 to generate a new error estimate;        11. Comparing the new error estimate to the old estimate and keeping the phase estimates associated with the lowest error; and        12. Continuing to execute steps 10 and 11 until the error is minimized and the best OTF estimate is obtained.        
Once the error-minimized OTF estimate is obtained, a Wiener filter can be generated that removes the effect of atmospheric turbulence. The OTF itself is generated by phase estimates that are due to atmospheric turbulence effects and aperture effects (for example, diffraction effects). If the effects of the atmospheric turbulence are mitigated by filtering them out using the Wiener filter, then the only remaining effect is that due to diffraction, and so the diffraction-limited result is obtained. To attempt to remove the effects of the aperture, super-resolution methods can be employed.
After the error-minimized OTF estimate is determined, the Wiener filter is given by,
                                                        H                              -                1                                      ⁡                          (                              f                →                            )                                =                                                    H                *                            ⁡                              (                                  f                  →                                )                                                    (                                                                                                              H                      ⁡                                              (                                                  f                          →                                                )                                                                                                  2                                +                α                            )                                      ,                            (        14        )            where the asterisk on the right side of equation (14) represents complex conjugation. Care must be taken for the case where the denominator of equation (14) approaches zero. A parameter a based on system noise is sometimes included in the denominator to prevent equation (14) from blowing up as H({right arrow over (f)}) approaches zero. As can be seen from equation (7), multiplying the image spectrum by the Wiener filter leads to an unaberrated object spectrum,O({right arrow over (f)})=I({right arrow over (f)})H−1({right arrow over (f)}),   (15)and the atmospheric turbulence free object brightness estimate is simply obtained by taking the 2-D inverse Fourier transform of equation (15).
Needed are embodiments of modified wavelength diversity compensation that overcome or mitigate one or more of the deficiencies of the prior art.