Joint contrast enhancement and turbulence mitigation method

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Patent Application Ser. No. 62/020,235 filed Jul. 2, 2014, entitled “Joint Contrast Enhancement and Turbulence Mitigation Method”, the content of which is fully incorporated by reference herein.

FEDERALLY SPONSORED RESEARCH AND DEVELOPMENT

The Joint Contrast Enhancement and Turbulence Mitigation Method is assigned to the United States Government. Licensing inquiries may be directed to Office of Research and Technical Applications, Space and Naval Warfare Systems Center, Pacific, Code 72120, San Diego, Calif., 92152; telephone (619) 553-5118; email: ssc_pac_t2@navy.mil. Reference Navy Case No. 103222.

BACKGROUND

A common problem for imaging in the atmosphere is fog and atmospheric turbulence. Previous research has focused on either the fog or turbulence aspects, but not both. A system and method are needed that can incorporate both the atmospheric turbulence and fog aspects for applications such as real-time video processing.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a diagram illustrating the processing for an embodiment of the Joint Contrast Enhancement and Turbulence Mitigation Method.

FIGS. 2 and 3 show diagrams illustrating embodiments of a system configured to perform the methods in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method.

FIGS. 4-7 show input images and output images of objects in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method.

FIG. 8 shows a graph illustrating turbulence mitigation metrics for various values of K in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method.

FIG. 9 shows a graph illustrating the peak-to-signal noise ratio for various values of K in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method.

FIG. 10 shows a flowchart of an embodiment of a method in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method.

DETAILED DESCRIPTION OF SOME EMBODIMENTS

Reference in the specification to “one embodiment” or to “an embodiment” means that a particular element, feature, structure, or characteristic described in connection with the embodiments is included in at least one embodiment. The appearances of the phrases “in one embodiment”, “in some embodiments”, and “in other embodiments” in various places in the specification are not necessarily all referring to the same embodiment or the same set of embodiments.

Some embodiments may be described using the expression “coupled” and “connected” along with their derivatives. For example, some embodiments may be described using the term “coupled” to indicate that two or more elements are in direct physical or electrical contact. The term “coupled,” however, may also mean that two or more elements are not in direct contact with each other, but yet still co-operate or interact with each other. The embodiments are not limited in this context.

As used herein, the terms “comprises,” “comprising,” “includes,” “including,” “has,” “having” or any other variation thereof, are intended to cover a non-exclusive inclusion. For example, a process, method, article, or apparatus that comprises a list of elements is not necessarily limited to only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Further, unless expressly stated to the contrary, “or” refers to an inclusive or and not to an exclusive or.

Additionally, use of the “a” or “an” are employed to describe elements and components of the embodiments herein. This is done merely for convenience and to give a general sense of the invention. This detailed description should be read to include one or at least one and the singular also includes the plural unless it is obviously meant otherwise.

The embodiments disclosed herein disclose a joint turbulence mitigation and fog removal method that can recover an object image o(x) fast enough for near real-time performance. The embodiments of the method are based upon analysis in turbulence mitigation and include the fog model. The embodiments of the method perform well for most atmospheric conditions and are efficient for near-real time processing.

Imaging systems positioned near the ocean often suffer in performance, perceptually and objectively, because of atmospheric turbulence, fog, sun-glare, camera motion from wind buffeting and many other adverse weather conditions. For long distance imaging, the most prominent are camera motion, fog, atmospheric turbulence and blur (from optics and atmosphere). An example of a specific imaging environment is an optical system that is observing targets of interest with the optical axis being parallel with the horizon (over-the-horizon viewing). The environment itself will have fog or haze, wind, and heat that causes eddy currents, which are observed as turbulence in an imaging system.

For each wavelength (color), the atmospheric effects in the imaging system are represented with
i_k(x)=h_k(x)*(t_k({tilde over (x)}_k)o_k({tilde over (x)}_k)+(1−t_k({tilde over (x)}_k))a)+w_k(x)  (Eq. 1)

The ideal image of the object without any atmospheric artifacts at frame k and pixel position xε² is o_k(x)ε. The observed image with the dominant artifacts of blur, turbulence, motion and fog is i_k(x)ε. The contrast degradation from fog is represented with the terms transmission t(x)ε and constant scalar airlight aε. The image noise that is independent of the atmosphere at frame k is w_k(x). The airlight is the light imaged when the object at x is infinitely far away from the imaging system. An example of a is the brightness observed of fog at the edge of the horizon.

The amount of light decay over distance from the atmospheric scattering medium is spatially varying and depends on the distance from the sensor to the object d(x)ε and scattering medium parameter βε, where the decay is characterized as
t_k(x)=e^−βdk(x).  (Eq. 2)

The global blur from the atmosphere and optical system is modeled with the blur kernel h (x). The observed motion from turbulence and camera motion is
{tilde over (x)}_k=Δ_g,k+Δ_t,k(x),  (Eq. 3)
with global motion represented with Δ_g,k and local motion from the object and turbulence is Δ_t,k(x). The blur, h(x), is modeled independent of local motion {tilde over (x)}. Although the blur caused by the exposure time is dependent on the turbulent local motion, the motion averaged together is modeled with the kernel h_k(x).

The performance of an optical system imaging through the atmosphere for short and long exposure settings has been sufficiently presented in the literature. A Fried seeing parameter r_o is used to describe the scale of the turbulence. If an optical system has a diameter smaller than r_o, then the turbulence is minimized. The modulation transfer function (MTF) of the optical system (including the atmosphere) degrades as the diameter increases above r_o. For the embodiments disclosed herein, the diameter of the optical system is assumed to be larger than r_o, which is a very common scenario.

The structure constant, C_n², from Kolmogorov theory is a parameter that “characterizes the strength of index-of-refraction fluctuations”. The structure constant is a function of altitude but is assumed constant in all horizontal viewing angles when an optical system is in an ocean environment.

An assumption is made that the imaging system is under isoplanatic conditions (spatially invariant blur) because the optical system is positioned close to the heat source (turbulent generator) which is most common for ground based optics in astronomy and imaging systems positioned near the ocean. Obtaining the image o_k(x) without knowledge of the blur h, transmission t, local motion {tilde over (x)}, or noise statistics is a difficult problem. Many methods have been proposed that account for some knowledge of the physics of the scene, however no method exists that accounts for both turbulence and fog. The embodiments of the method disclosed herein address how to incorporate information from the contrast enhancement process into the turbulence mitigation process.

FIG. 1 shows a diagram 10 illustrating the processing for an embodiment of the Joint Contrast Enhancement and Turbulence Mitigation Method. Given a foggy turbulent image at frame k i_k(x), single image noise removal (denoising) is applied at box 20 and passed into a single image defogging algorithm at box 30. The defogging algorithm estimates the transmission and t_max and enhances the contrast image, which then goes into the image alignment algorithm at box 40 to estimate the global motion vector {circumflex over (Δ)}_g,k. The enhanced image is then aligned with {circumflex over (Δ)}_g,k and locally aligned with an optical flow method at box 50 which estimates {circumflex over (Δ)}_t,k(x). The globally and locally aligned contrast enhanced image is then temporally averaged with K frames at box 60 and atmospheric blur is removed at box 70 to produce the estimated image of the object ô_k(x). A turbulence mitigation metric (TMM) is used at box 80 to measure performance by measuring the turbulence between the input and output image sequences. Feedback 82 regarding the TMM may be provided to the frame averaging step at box 60 to automatically change the number of frames to average K.

Contrast enhancement does not improve the image alignment performance when image noise is present. The result in terms of error is the same if enhancement is before or after global tracking Analysis does show, however, that reducing image noise can improve the image alignment performance. Additionally, removing noise before fog removal is a useful approach to improve transmission {circumflex over (t)} and airlight â estimation. Accordingly, some embodiments of the method disclosed herein remove noise, enhance the contrast and then estimate the global motion {circumflex over (Δ)}_g. A data dependent noise removal algorithm may be used because the amount of noise energy that will be enhanced will depend on the transmission. A median filter may also be used because it can preserve edges reasonably well and it is data dependent thus adapts to the image signal.

To reduce the complexity of the algorithm such that the processing speed is near real-time, a two dimensional median filter may be used for each color channel for fast single image denoising. The denoising step should not be too aggressive in reducing spatial frequencies because important details can be lost and are unable to be recovered during the frame averaging and deblurring steps. If noise removal is too aggressive, it effectively imposes premature bandwidth limitation and high frequency regions will be lost before the bandwidth recovery process (deblurring). As an example, a very small window 3×3 may be used for the median filtering step for denoising.

Some embodiments of the disclosed method may employ the adaptive correlation filter, Minimum Output Sum of Squared Error (MOSSE) to estimate the global motion parameter Δ_g. This filter is fast and efficient and works well with turbulent video because the phase correlation update averages out the zero mean phase perturbations induced by turbulence. The local motion induced by turbulence must also be compensated in order to have a sharper image after frame averaging. Some embodiments involve enhancement before optical flow in order to improve the motion estimation. Enhancing before motion compensation reduces the number of possible motion vectors for each pixel. For video coding, intra-frame coding is used more when enhancing first therefore details are preserved whereas inter-frame coding is used more when there is low contrast which results in loss of details.

To remove fog effectively and efficiently for each frame and have an estimate of transmission, the method disclosed in U.S. patent application Ser. No. 14/501,400, titled “Single Image Contrast Enhancement Method Using the Adaptive Wiener Filter” of Gibson et al., may be used. The Wiener defogging method is 50 to 100 times faster than existing methods and can operate at real-time speeds for frames of sizes 720×480 and larger.

A modification to the Wiener defog method may be made in order to automate the defogging process. One step in removing fog with the Wiener defogging method is the proper choice of size for the smoothing window Ω in the spatial domain. The optimal Ω must be heuristically chosen for each image. This can be a cumbersome task for real-time video enhancements. The choice of window size is a common problem in single image defogging. A choice of two window sizes may be performed to properly refine the transmission estimate, but yet preserve depth discontinuities. The window should be sufficiently large to estimate a smooth transmission map that best reflects the scene depth.

Some embodiments of the disclosed method may automatically select |Ω| based on the statistics of the frames. Given a foggy color image iε^ψ, fog may be removed by estimating the transmission t and recovering the true non-foggy image oε^ψ. The number of color channels is denoted with ψ. While disclosure is made with regard to a 3-channel color image, the embodiments of the method can be used for grey images (ψ=1). As an assumption is made that the transmission is independent of wavelength, the image model (without motion, turbulence and blur) is the atmospheric dichromatic model
i(x)=t(x)o(x)+(1−t(x))a+w(x).  (Eq. 4)

The main idea in estimating the transmission with an adaptive local Wiener filter is that the observation model is composed of the veiling v=(1−t) and texture n_txt. The observation, the dark pixel measurement for each pixel, is

d ⁡ ( x ) = min c ∈ ( r , g , b ) ⁢ i c ⁡ ( x ) , ( Eq . ⁢ 5 )
where i_c(x) is the c^th color channel of the image i(x). For grey images, d(x)=i(x)ε. The dark pixel measure in (Eq. 5) takes advantage of the atmospheric dichromatic model in (Eq. 5) by assuming at least one color component is possibly dark, therefore exposing the atmospheric veiling: d(x)˜(1−t (x))a_c+w_c(x) for min_cεr,g,bo_c(c)˜0. Thus, the model at each pixel is
d(x)=v(x)+n_txt  (Eq. 6)
where v is the true veiling and n_txt is the textured noise. The observation is (Eq. 5) and the model is (Eq. 6). The goal in the Wiener defog method is to filter out the texture from the observation dby treating the texture n_txt as noise but preserve edges from depth discontinuities. For images with a large amount of texture (e.g. forest floor with leaves), the size of the window must be sufficiently large to filter out the texture variations. The choice of Ω arises from the local moment estimators

μ ^ v ⁡ ( x ) = 1  Ω ⁡ ( x )  ⁢ ∑ j ∈ Ω ⁡ ( x ) ⁢ d ⁡ ( x ) , ( Eq . ⁢ 7 )

σ ^ v 2 ⁡ ( x ) = 1  Ω ⁡ ( x )  ⁢ ∑ j ∈ Ω ⁡ ( x ) ⁢ d 2 ⁡ ( x ) - μ ^ v 2 ⁡ ( x ) , ( Eq . ⁢ 8 )
where |Ω(x)| is the number of pixels within the local sample window positioned at a pixel location x.

The approach to estimating the optimum window size |Ω| is to choose a size sufficient enough to filter scene texture noise from the observation {circumflex over (μ)}_v in (Eq. 7) to a desired variance σ_e,Ω²,

σ e , Ω 2 = σ ^ n txt 2  Ω  ⟹  Ω  = σ n txt 2 σ e , Ω 2 ( Eq . ⁢ 9 )
As an example σ_e,Ω² may be equal to 5e⁻⁶. To determine |Ω|, the variance of the texture is estimated. Taking the maximum variance from the local variance measurement {circumflex over (σ)}_v²(x) in (Eq. 8) is not practical since the higher variance measurements are from the depth discontinuities and edges in the scene. The goal in smoothing d is to preserve these edges. One approach is to estimate the texture noise in an image. The approach is to measure the average variance in all color channels from the observed image i(x),

σ n txt 2 = 1 3 ⁢ 1 M ⁢ ∑ c = 1 3 ⁢ ∑ ∀ x ⁢ σ ^ i , c 2 ⁡ ( x ) , ( Eq . ⁢ 10 )
where

μ ^ i , c 2 ⁡ ( x ) = 1  Φ ⁡ ( x )  ⁢ ∑ y ∈ Φ ⁡ ( x ) ⁢ i c ⁡ ( y ) , ( Eq . ⁢ 11 )

σ ^ i , c 2 ⁡ ( x ) = 1  Φ ⁡ ( x )  ⁢ ∑ y ∈ Φ ⁡ ( x ) ⁢ i c 2 ⁡ ( y ) - μ ^ i , c 2 ⁡ ( x ) , ( Eq . ⁢ 12 )
and where M is the total number of pixels in an image. A sample window size |Φ| is selected initially in order to estimate the local mean and variance of the image. As an example, a sample window Φ=32×32 may be used for all image sequences. The sample window size |Ω| is estimates using (Eqs. 9 and 10) with σ_e,Ω²=5e⁻⁶ for images scaled to pixel values to [0,1].

The Wiener defogging method estimates the fog-free image ô and transmission {circumflex over (t)},
ô,{circumflex over (t)}←WienerDefog(i,Ω,t₀,d_p,1),  (Eq. 13)
where Ω is the |Ω|^1/2×|Ω|^1/2 size sample window, t₀ is the minimum possible transmission value used for numerical conditioning, d_p is the weight used to change the strength of the enhancement and the last parameter, 1, indicates the single-step approach is used. Typical values for t₀ and d_p are 1e⁻² and 0.9, respectively.

It is observed that the transmission estimation is invalid when there is a significant variation in fog causing the model described by (Eq. 6) to be invalid because the veiling and noise are not statistically independent. The texture variation n_txt is dependent on the amount of fog. Thus, the first step naively assumes the model is valid in order to remove fog and to provide the ability to re-estimate the local moments o_n² and o_v² that are indeed statistically uncorrelated in the 2-step approach. It also has been observed that images with significant depth discontinuities (e.g. buildings in near foreground with sky as backdrop), the 2-step approach is needed. The one-step approach may be sufficient for most of the embodiments disclosed herein because the noise is not significantly dependent on fog.

The Wiener defogging method estimates air light where it is assumed that if the sky is present in an image (t=0), the brightest pixels are the actual air light color a. A practical implementation is to pick the 0.1% darkest pixels in the transmission estimate {circumflex over (t)} to estimate the airlight a=(a₁, a₂, a₃)^T,

a c = 1000 M ⁢ ∑ y ⁢ i c ⁡ ( y ) : t ^ ⁡ ( y ) ≤ T , ( Eq . ⁢ 14 )
where T represents a threshold chosen such that there are M/1000 pixels that have a value lower than T.

The variance of temporally averaging (denoted with angle brackets) a defogged image sequence that has already been aligned with motion compensation is

var ⁡ [ 〈 defog ⁡ ( i ⁡ ( x ) ) 〉 ] ≈ 1 K ⁢ ( σ o 2 ⁡ ( x ) + σ ω 2 μ t 2 ⁡ ( x ) ) . ( Eq . ⁢ 15 )
The noise variance is reduced with K but amplified with <t(x)>_K²=μ_t²(x). In order to reduce the noise power, the averaging parameter K must depend on the time averaged transmission μ_t(x). One approach is to choose a large K value, but this is not ideal for most scenes where a smaller K is a better choice to prevent excessive blur.

The fewer the number of frames used for averaging, the less blur will occur if the frames are well aligned. A global averaging method helps to reduce flickering and achieve low complexity however adaptive to the signal noise and contrast. An approach to solve the problem is to provide a method for determining the lower bound for K given a desired o_e,K². This is achieved by using the maximum transmission value t_max estimated from the image sequence. Therefore, the minimum number of frames to average for an image sequence with M pixels per frame can be determined with a contrast dependent measure and averaged o_o²,

K min = 1 o e , K 2 ⁢ ( 1 M ⁢ ∑ x ⁢ σ ^ o 2 ⁡ ( x ) + σ ω 2 ^ t max 2 ) , ( Eq . ⁢ 16 )
with t_max:t_max>{circumflex over (t)}(x)∀x. As an example, o_e,K²=2×10⁻⁵, but other values may be used as determined by a user.

After the image i(x) has been defogged, motion compensated, and frame averaged, the next step is to remove the blur caused by the atmosphere and image alignment errors (deblur). In some embodiments, a Wiener Deconvolution method may be used for this purpose. In order to remove the blur, a kernel must be known. One approach is to apply a blind deconvolution method. For isoplanatic conditions, the blur kernel is not exactly known but the structure is. There are two different forms of the optical transfer function (OTF), or blur kernel, due to turbulence: one for long exposure and the second for short exposure. The exposure time is in relation to the atmospheric turbulence. Short exposure measurements effectively freeze the turbulence during the image capture integration time. Only one instance of turbulence is captured in a short exposure frame. Long exposure is an integration of multiple turbulent instances in one frame capture. In embodiments of the method herein, the long exposure model may be used since the frames are averaged together.

Given an OTF that represents the blur from turbulence, H (u, λ), where u is the spatial frequency and λ is the optical wavelength (color channel), the long exposure OTF is

H LE ⁡ ( u , λ ) = e - 1 2 ⁢ 6.88 ⁢ ( λ ⁢ ⁢ l r ⁢ ⁢ o ⁢  u  ) 5 / 3 ( Eq . ⁢ 17 )
The optic focal length is l. The parameter α is for near (α=1) and far field (α=0.5). From experiments, the far field imaging is addressed, therefore α=0.5. The Fried parameter is r₀. Modifying the ratio l/r₀ for both short and long exposure is sufficient for situations when the optical system and seeing parameters are not known, and in fact is a common situation because r₀ is typically unknown.

To illustrate that the turbulence has been mitigated in time, a Turbulence Mitigation Metric (TMM) is presented. The goal in turbulence mitigation is to not only recover the true structure (phase) of the scene but also to enforce temporal consistency in each color channel. Instead of using a subjective measure by displaying multiple frames of a video sequence, an objective measure is developed. The goal in designing the TMM is to make it an objective measure of turbulence mitigation, simple to construct, and easily accessible such that it can be used in any other turbulence mitigation method. One approach to developing a TMM is to utilize the motion estimated from optical flow and global estimations. This approach, however, is complex given that not every method uses optical flow. It is desired to only require an input image sequence and output image sequence.

The TMM is a measurement between the input image sequence i and the output image sequence o,

TMM ⁡ ( i , o ) =  i  TV ⁢ ⁢ 1 -  o  TV ⁢ ⁢ 1  i  TV ⁢ ⁢ 1 , ( Eq . ⁢ 18 )
where the TV/L1 norm, ∥f∥_TV1:^M1×M2×3×K→, is

 f  TV ⁢ ⁢ 1 = ∑ x , λ ⁢ ∑ k = 0 K - 1 ⁢  f ⁡ ( x , λ , k ) - f ⁡ ( x , λ , k - 1 )  ( Eq . ⁢ 19 )

where each image sequence is of spatial dimension M₁×M₂, has 3 color (λ) channels, and is of length K frames. This is simply a total variation (TV) measurement in the time domain similar to the anisotropic TV norm where the spatial weights are set to zero and the time weight factor is set to one. Additionally, (Eq. 19) accumulates the variation for all color channels.

If the turbulence is reduced, then the value of the TMM approaches unity. The TMM may be negative if the input sequence exhibits less turbulence than the output sequence. When contrast enhancement is involved, it is important to measure ∥f∥_TV1 immediately after contrast enhancement. This is because the color variances are smaller when fog is present and the color variations in time are smaller compared to the same image sequence with the contrast enhanced, which would produce erroneous negative TMM values.

FIG. 2 shows a diagram 100 illustrating an embodiment of a system 110 configured to perform the methods in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. System 110 includes a processor 112 operatively connected, such as via a wired or wireless connection, to a memory module 114 and display 116. Processor 112 may be a general purpose processing device or a processing device specifically programmed to perform the embodiments of the method disclosed herein. Memory module 114 is configured to store various types of data therein and may include computer-implementable instructions stored therein for performing the embodiments of the method, such as method 900, disclosed herein. Display 116 may comprise a general purpose display device as recognized by one having ordinary skill in the art. System 110 is operatively connected, such as via a wired or wireless connection, to an imaging device 120, which may comprise any type of image capture device known in the art, such as a video camera, smart phone, digital camera, or the like.

In operation, imaging device 120 may provide image data, either delayed or in real-time, to processor 112. Processor 112 may then execute computer-implementable instructions stored therein or stored within memory module 114, to remove any contrast reducing scattering medium present in the image data and cause a contrast enhanced image to be displayed to a user on display 140.

FIG. 3 shows a diagram 200 of illustrating an embodiment of a system 210 configured to perform the methods in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. System 210 includes a processor 212 operatively connected, such as via a wired or wireless connection, to a memory module 214 and imaging device 216. The components shown in FIG. 3 may be the same as or similar to the components shown in FIG. 2, with the difference being that system 210 includes imaging device 216, while system 210 is operatively connected to imaging device 220 and may include a separate display.

In operation, imaging device 216 may provide image data, either delayed or in real-time, to processor 212. Processor 212 may then execute computer-implementable instructions stored therein or stored within memory module 214, to remove any contrast reducing scattering medium present in the image data and cause a contrast enhanced image to be displayed to a user on imaging device 216.

Systems 100 and 200 may be configured with the appropriate software modules to perform the embodiments of the method, such as method 900, as discussed herein. As an example, processors 112 and 212 may have such modules stored therein, or, memory modules 114 and 214 may have such modules stored therein, which are accessed by processors 112 and 212, respectively. The term “module” generally refers to a software module. A module may be implemented as a collection of routines and data structures that performs particular tasks or implements a particular abstract data type. Modules generally are composed of two parts. First, a software module may list the constants, data types, variables, and routines that may be accessed by other modules or routines. Second, a module may be configured as an implementation, which may be private (i.e., accessible only to the module), and which contains the source code that actually implements the routines or subroutines upon which the module is based. Thus, the use of the term “module” herein, indicates reference to such software modules or implementations thereof.

FIGS. 4-7 show input images and output images of objects in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. FIG. 4 shows a diagram 300 with an input image 310 of a boat within a contrast reducing scattering medium and an output image 320 representing the processed image in accordance with the method disclosed herein. FIG. 5 shows a diagram 400 with an input image 410 of a star image within a contrast reducing scattering medium and an output image 420 representing the processed image in accordance with the method disclosed herein. FIG. 6 shows a diagram 500 with an input image 510 of a tower within a contrast reducing scattering medium and an output image 520 representing the processed image in accordance with the method disclosed herein. FIG. 7 shows a diagram 600 with an input image 610 of a chimney within a contrast reducing scattering medium and an output image 620 representing the processed image in accordance with the method disclosed herein.

FIG. 8 shows a graph 700 illustrating TMMs for various values of K in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. Line 710 represents the TMM for K=10, line 720 represents the TMM for K=44, and line 730 represents the TMM for K=90. FIG. 9 shows a graph 800 illustrating the peak-to-signal noise ratio (PNSR) for various values of K in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. Line 810 represents the PNSR for K=10, line 820 represents the PNSR for K=44, and line 830 represents the PNSR for K=90.

FIG. 10 shows a flowchart of an embodiment of a method 900 in accordance with the Joint Contrast Enhancement and Turbulence Mitigation Method. Method 900 disclosed herein may be implemented as a series of modules, either functioning alone or in concert with physical electronic and computer hardware devices. Such modules may be utilized separately and/or together locally and/or remotely to form a program product thereof, that may be implemented through recordable media. Method 900 may be implemented as a program product comprised of a plurality of such modules, which can be interactively displayed for a user on a display, such as display 116 or imaging device 216.

Some or all of the steps of method 900 may be stored on a non-transitory computer readable storage medium, wherein the steps are represented by computer-readable programming code. The steps of method 900 may also be computer-implemented using a programmable device, such as a computer-based system. Method 900 may comprise instructions that, when loaded into a computer-based system, cause the system to execute the steps of method 900. Method 900 may be computer-implemented using various programming languages, such as “Java”, “C”, or “C++”. Further, method 900 may be implemented within a system such as system 100 shown in FIG. 2 or system 200 shown in FIG. 3.

For illustrative purposes, method 900 will be discussed with reference to the steps being performed by processors 112 or 212, using instructions stored therein or memory modules 114 or 214, respectively, with the image input provided by imaging devices 120 or 216, respectively and the image output being displayed on display 116 or imaging device 216, respectively.

Additionally, while FIG. 9 shows one embodiment of method 900 to include steps 910-970, other embodiments of method 900 may contain fewer or more steps. Further, while in some embodiments the steps of method 900 may be performed as shown in FIG. 9, in other embodiments the steps may be performed in a different order, or certain steps may occur simultaneously with one or more other steps.

Method 900 begins at step 910, which involves generating a denoised image by performing a denoising operation on an input image of an object ô_k(x). The input image sequence, such as image 310, 410, 510, and 610, includes turbulence and a contrast reducing scattering medium therein, examples of which include fog, haze, water, or dust. In some embodiments, the denoising operation is performed by a median filter, such as a two-dimensional median filter.

Step 920 involves generating a contrast enhanced image from the denoised image using a contrast enhancement algorithm. In some embodiments, the contrast enhancement algorithm is a locally adaptive wiener filter. It should be noted that other contrast enhancement algorithms may be used including, but not limited to, histogram equalization, contrast limited histogram equalization, and the median dark channel prior method.

Step 930 involves estimating a global motion vector {circumflex over (Δ)}_g,k of the contrast enhanced image using an image alignment algorithm. In some embodiments, the image alignment algorithm is a minimum output sum of squared error adaptive correlation filter, such as is described by Bolme et al. Other image alignment algorithms may be used, including, but not limited to, phase correlation tracking Kanade-Lucas-Tomasi feature tracking, and feature tracking using SURF.

Step 940 involves generating a contrast enhanced aligned image by globally aligning the contrast enhanced image to {circumflex over (Δ)}_g,k and locally aligning the contrast enhanced image with an optical flow method. In some embodiments, step 940 is performed by the image alignment algorithm. In some embodiments, the optical flow method is a dense optical flow block-based method so that every pixel has a motion vector associated with it. The motion vectors are used to warp the images to be aligned with each other.

Step 950 involves temporally averaging the contrast enhanced aligned image with K frames. In some embodiments, the step of temporally averaging the contrast enhanced aligned image comprises determining a lower bound for, K_min, according to the equation K_min=

1 o e , K 2 ⁢ ( 1 M ⁢ ∑ x ⁢ σ ^ o 2 ⁡ ( x ) + σ ^ ω 2 t max 2 ) ,
with t_max:t_max>{circumflex over (t)}(x) for all x, where M is the total number of pixels in the contrast enhanced aligned image, where t_max is the maximum transmission value, o_e,K² is a desired variance (i.e. 2×10^—5), {circumflex over (σ)}_o²(x) is the variance of the defogged image in time at pixel x, {circumflex over (σ)}_ω² is the noise variance in time, {circumflex over (t)}(x) is the estimated transmission at pixel location x. In some embodiments, the temporal averaging is performed by averaging each pixel through time.

In some embodiments, the number of frames K to average is automatically determined using feedback regarding the TMM determined at step 970. This determination may be performed by a user setting a desired TMM, e.g. 0.95 or other desired value, and the algorithm will either decrease or increase the number of frames K to reach the desired TMM value of 0.95 within a reasonable tolerance, which may also be user selected.

Step 960 involves generating an output image, such as images 320, 420, 520, and 620, of the object ô_k(x) by performing a deblur operation on the temporally averaged contrast enhanced aligned image. In some embodiments, the deblur operation involves a the Richardson-Lucy Deconvolution method.

Step 970 involves measuring turbulence using a turbulence mitigation metric (TMM). In some embodiments, step 970 comprises taking a measurement between an input image sequence i and an output image sequence o, where each of i and o are of spatial dimension M₁×M₂, have three color channels, and are of length K frames. In some embodiments, the TMM is measured using a sequence of input images that each have been denoised and defogged, such as shown in FIG. 1. In some embodiments, method 900 may further comprise the step of outputting the output image to a display, such as display 116 or imaging device 216. In some embodiments, the TMM may be used as a way to change the number of frames K, via feedback provided to step 950 regarding the determined TMM. If a desired TMM value is set, the algorithm could change K to meet that specific TMM.

Various storage media, such as magnetic computer disks, optical disks, and electronic memories, as well as non-transitory computer-readable storage media and computer program products, can be prepared that can contain information that can direct a device, such as a micro-controller, to implement the above-described systems and/or methods. Once an appropriate device has access to the information and programs contained on the storage media, the storage media can provide the information and programs to the device, enabling the device to perform the above-described systems and/or methods.

For example, if a computer disk containing appropriate materials, such as a source file, an object file, or an executable file, were provided to a computer, the computer could receive the information, appropriately configure itself and perform the functions of the various systems and methods outlined in the diagrams and flowcharts above to implement the various functions. That is, the computer could receive various portions of information from the disk relating to different elements of the above-described systems and/or methods, implement the individual systems and/or methods, and coordinate the functions of the individual systems and/or methods.

Many modifications and variations of the Joint Contrast Enhancement and Turbulence Mitigation Method are possible in light of the above description. Within the scope of the appended claims, the embodiments of the systems described herein may be practiced otherwise than as specifically described. The scope of the claims is not limited to the implementations and the embodiments disclosed herein, but extends to other implementations and embodiments as may be contemplated by those having ordinary skill in the art.