System and methods of regularized optimization for matrix factorization and image and video reconstruction

FIELD OF THE INVENTION

p-0002The present invention relates generally to image and video processing. More particularly, the present invention is directed to a system and methods of optimizing Amplitude-Modulation Frequency-Modulation (“AM-FM”) demodulation for processing stationary and non-stationary image and video content.

p-0003AM-FM demodulation is useful in a variety of contexts and applications including, for example, characterization and classification of image and video from imaging modalities such as electron microscopy, spectral and hyperspectral devices, ultrasound, magnetic resonance imaging (“MRI”), positron emission tomography (“PET”), histology, color and monochrome images, molecular imaging, radiographs (“X-rays”), computer tomography (“CT”), and others. The specific applications are in fingerprint identification, detection and diagnosis of retinal disease, malignant cancer tumors, cardiac image segmentation, atherosclerosis characterization, brain function, histopathology specimen classification, characterization of anatomical structure tracking such as carotid artery walls and plaques or cardiac motion and as the basis for computer-aided diagnosis to name a few.

BACKGROUND OF THE INVENTION

p-0004Image and video processing are forms of signal processing. Signal processing allows a set of characteristics or parameters related to the image or video to be obtained. Signal processing including analog signal processing, discrete time signal processing, and digital signal processing, which may involve a one-dimensional (“1D”), two-dimensional (“2D”) or three-dimensional (“3D”) input signal to which signal processing techniques are applied.

p-0005Signal processing techniques include transform-based processing such as discrete or integral transforms which were implemented prior to AM-FM processing. As an example, a 1D analysis of transform-based processing includes the use of short-time Fourier Transform (“STFT”) for non-stationary signals. When using STFT, the fast Fourier Transform (“FFT”) of different time intervals of the signals is used to determine the frequency and phase content. Thus, the STFT is a convenient 2D representation that provides frequency content information at different time intervals. A disadvantage is that the STFT cannot be effectively generalized to images and videos. For example, using STFT for images would produce a four-dimensional (“4D”) representation and using STFT for video would produce a six-dimensional (“6D”) representation.

p-0006The discrete Wavelet Transform (“DWT”) has also been used for transform-based image processing. Unlike Fourier Transforms, Wavelet Transforms are based on specific functions defined at different scales and durations. Thus, the DWT is a space-frequency representation of the input signal and it is related to harmonic analysis is as in Fourier Transform. While FFT uses equally spaced frequency division, DWT uses logarithmic divisions of the frequency. A disadvantage is that DWT does not measure frequency content directly.

p-0007The development of accurate methods for estimating amplitude-modulation frequency-modulation image decompositions is of great interest due to is potentially significant impact on image analysis applications including in the areas of signal, image and video processing. Applications in signal processing include speech signal analysis. Image processing applications include shape from shading, image pattern analysis, image interpolation, fingerprint classification, image retrieval in digital libraries, image segmentation, and damaged image texture repairs. Applications in video processing include cardiac image segmentation, motion estimation, and motion reconstruction, to name a few.

p-0008A number of techniques exist to reconstruct an image from its AM-FM representation in terms of amplitude, phase and frequency functions. Reconstruction of an image involves estimating or computing the amplitude, phase and frequency components of the signals emerging from each filter channel and using these components to create an AM-FM representation that best approximate the original image signal. Generally, the more components and channels used, the more information is recovered, the better the image signal is restored and the better the image is regenerated. If every component of every channel is used, this will yield to the best reconstruction of the original image. However, such approaches lead to very redundant representations that can lead to very inefficient applications in image analysis. Thus, the goal of an efficient reconstruction process is select few channels and components that best approximates the image signal.

p-0009The AM-FM Dominant Component Analysis (“DCA”) and Channelized Component Analysis (“CCA”) are methods used that consist of applying a filterbank to the Hilbert-transformed image, and then applying AM-FM demodulation of each bandpass filtered image. Using DCA, every pixel delivers estimated modulating functions corresponding to the AM-FM component that is locally dominant at that pixel. Using CCA, a filterbank partitions the image into components on a spatially global basis. Each resulting AM-FM component is restricted to lie in a single channel over the entire image domain. With CCA, the number of components in the computed image model is necessarily equal to the number of channels in the filterbank. AM-FM reconstructions based on the CCA use a reasonably small number of locally coherent components. In contrast, those based on the DCA only use one component—the estimates from the channel with the maximum amplitude estimate. A disadvantage of DCA and CCA are that they are known to produce noticeable visual artifacts.

p-0010Optimizing the quality of an AM-FM reconstruction image is important due to the potentially significant impact on various applications. Thus, there is demand for high quality reconstructions in both stationary and non-stationary processing for use in a variety of contexts and applications. The present invention satisfies this demand.

SUMMARY OF THE INVENTION

p-0011The present invention provides a system and methods of high quality reconstructions, both visually and quantitatively, when compared to standard reconstructions using various prior art techniques such as Dominant Component Analysis (“DCA”) and Channelized Component Analysis (“CCA”). The present invention is based on a Regularized Optimization (“RO”) to attain a small number of locally coherent components and simultaneously enforce a piecewise smooth constrain for one or more amplitude functions. In one embodiment, the small number of locally coherent components and piecewise smooth constrain for one or more amplitude functions is based on the estimates from the CCA. Image content from image signal components is obtained from processing a Hilbert-transformed image through a filter bank. More particularly, the present invention provides a Regularized Optimization (“RO”) method of reconstructing an image by applying an amplitude-modulation frequency-modulation (“AM-FM”) demodulation process to an image.

p-0012Although the present invention is discussed herein with respect to two-dimensional (“2D”) images, signals, and digital videos, it is contemplated the present invention can be extended to any dimensional (“ND”) images, signals and digital videos including three-dimensional (“3D”).

p-0013The AM-FM representation of images permits non-stationary image content to be modeled in terms of amplitude and phase functions using the following equation:

p-0014<maths id="MATH-US-00001" num="00001"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>b</mi><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>L</mi></munderover><mo>⁢</mo><mrow><mrow><msub><mi>α</mi><mi>n</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>⁢</mo><mrow><mi>cos</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><msub><mi>φ</mi><mi>n</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> Where b(ξ):<img id="CUSTOM-CHARACTER-00001" he="3.56mm" wi="3.56mm" file="US08908992-20141209-P00001.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />→<img id="CUSTOM-CHARACTER-00002" he="3.13mm" wi="2.46mm" file="US08908992-20141209-P00002.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" /> is the input image, ξ=(ξ₁,ξ₂)ε<img id="CUSTOM-CHARACTER-00003" he="3.56mm" wi="3.56mm" file="US08908992-20141209-P00001.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />, Lε<img id="CUSTOM-CHARACTER-00004" he="3.13mm" wi="3.56mm" file="US08908992-20141209-P00003.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />, α_n:<img id="CUSTOM-CHARACTER-00005" he="3.56mm" wi="3.56mm" file="US08908992-20141209-P00001.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />→[0,∞)⁻, and φ_n:<img id="CUSTOM-CHARACTER-00006" he="3.56mm" wi="3.56mm" file="US08908992-20141209-P00001.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />→<img id="CUSTOM-CHARACTER-00007" he="3.13mm" wi="2.46mm" file="US08908992-20141209-P00002.TIF" alt="custom character" img-content="character" img-format="tif" orientation="portrait" inline="no" />. The interpretation of Equation (1) suggest that the L AM-FM component images α_n(ξ)·cos(φ_n(ξ)), model the essential image modulation structure, the amplitude functions α_n(ξ) model stow-changing image intensity variations, and the FM components cos(φ_n(ξ)) capture cast-changing image intensity variations.

p-0015It is contemplated that Equation (1) can also be interpreted as a separation of texture—FM components cos(φ_n(ξ))—from piecewise smooth content—amplitude functions α_n(ξ)—in an image.

p-0016The AM-FM Dominant Component Analysis (“DCA”) and Channelized Component Analysis (“CCA”) consist of applying a collection of filterbanks—bandpass filters—to the original input image. The AM-FM demodulation of each bandpass filtered image provides estimations of instantaneous amplitude (“IA”) functions α_n(ξ), instantaneous phase (“IP”) functions φ_n(ξ), and instantaneous frequency (“IF”) functions ω_n(ξ)=∇φ_n(ξ).

p-0017The goal using CCA is to obtain a reasonably small number of locally coherent components such as modeling the input image as in Equation (1). The goal using DCA is to select the estimates from the channel with the maximum amplitude estimate using one component—the dominant component—to model the input image.

p-0018According to the present invention, the minimum of the following:

p-0019<maths id="MATH-US-00002" num="00002"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><mi>J</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>a</mi><mo>,</mo><mi>ζ</mi></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mrow><mfrac><mn>1</mn><mi>p</mi></mfrac><mo>⁢</mo><msubsup><mrow><mo></mo><mrow><mrow><mi>f</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>a</mi><mo>,</mo><mi>ζ</mi></mrow><mo>)</mo></mrow></mrow><mo>-</mo><mi>b</mi></mrow><mo></mo></mrow><mi>p</mi><mi>p</mi></msubsup></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>a</mi></msub><mo>⁢</mo><mrow><mi>T</mi><mo>⁡</mo><mrow><mo>(</mo><mi>a</mi><mo>)</mo></mrow></mrow></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>ζ</mi></msub><mo>⁢</mo><msub><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mn>1</mn></msub></mrow></mrow></mrow><mo>,</mo><mrow><mrow><mi>s</mi><mo>.</mo><mi>t</mi><mo>.</mo><mi>a</mi></mrow><mo>≥</mo><mn>0</mn></mrow><mo>,</mo><mrow><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mo>≤</mo><mn>1</mn></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> attains a small number of locally coherent components and simultaneously enforces a piecewise smooth constrain for α_n(ξ).

p-0020As used herein a_n, ζ_nand b are one-dimensional (“1D”) vectors representing the two-dimensional (“2D) instantaneous amplitude function α_n(ξ), the two-dimensional function cos(φ_n(ξ)) and the image b(ξ). In Equation (2), a equals

p-0021<maths id="MATH-US-00003" num="00003"><math overflow="scroll"><mrow><mo>[</mo><mtable><mtr><mtd><mrow><msubsup><mi>a</mi><mn>1</mn><mi>T</mi></msubsup><mo>,</mo></mrow></mtd><mtd><mrow><msubsup><mi>a</mi><mn>2</mn><mi>T</mi></msubsup><mo>,</mo></mrow></mtd><mtd><mi>…</mi></mtd><mtd><mrow><msup><mrow><msubsup><mi>a</mi><mi>L</mi><mi>T</mi></msubsup><mo>]</mo></mrow><mi>T</mi></msup><mo>,</mo><mi>ζ</mi></mrow></mtd></mtr></mtable></mrow></math></maths><br /> equals

p-0022<maths id="MATH-US-00004" num="00004"><math overflow="scroll"><mrow><mo>[</mo><mtable><mtr><mtd><mrow><msubsup><mi>ζ</mi><mn>1</mn><mi>T</mi></msubsup><mo>,</mo></mrow></mtd><mtd><mrow><msubsup><mi>ζ</mi><mn>2</mn><mi>T</mi></msubsup><mo>,</mo></mrow></mtd><mtd><mrow><msup><mrow><mrow><mi>…</mi><mo>⁢</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>,</mo><mstyle><mspace width="0.6em" height="0.6ex" /></mstyle><mo>⁢</mo><msubsup><mi>ζ</mi><mi>L</mi><mi>T</mi></msubsup></mrow><mo>]</mo></mrow><mi>T</mi></msup><mo>,</mo></mrow></mtd></mtr></mtable></mrow></math></maths><br /> and f(a,ζ) equals

p-0023<maths id="MATH-US-00005" num="00005"><math overflow="scroll"><mrow><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>L</mi></munderover><mo>⁢</mo><mrow><mrow><mi>diag</mi><mo>⁡</mo><mrow><mo>(</mo><msub><mi>a</mi><mi>n</mi></msub><mo>)</mo></mrow></mrow><mo>*</mo><mrow><msub><mi>ζ</mi><mi>n</mi></msub><mo>.</mo></mrow></mrow></mrow></math></maths><br /> The TV regularization generalization to vector-valued images with coupled channels is

p-0024<maths id="MATH-US-00006" num="00006"><math overflow="scroll"><mrow><mrow><mi>T</mi><mo>⁡</mo><mrow><mo>(</mo><mi>a</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mfrac><mn>1</mn><mi>q</mi></mfrac><mo>⁢</mo><mrow><msubsup><mrow><mo></mo><msqrt><mrow><mrow><munder><mo>∑</mo><mi>n</mi></munder><mo>⁢</mo><msup><mrow><mo>(</mo><mrow><msub><mi>D</mi><mi>x</mi></msub><mo>⁢</mo><msub><mi>a</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow><mo>+</mo><msup><mrow><mo>(</mo><mrow><msub><mi>D</mi><mi>y</mi></msub><mo>⁢</mo><msub><mi>a</mi><mi>n</mi></msub></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></msqrt><mo></mo></mrow><mi>q</mi><mi>q</mi></msubsup><mo>.</mo></mrow></mrow></mrow></math></maths><br /> Horizontal discrete derivative operators are represented by D_xand vertical discrete derivative operators are represented by D_y.

p-0025The present invention is useful in a variety of contexts and applications, and allows the identification of disease at different stages, such as retinal disease (diabetic retinopathy, age-related macular degeneration, glaucoma, etc.), pulmonary diseases (pneumoconiosis, lung nodules tumors, etc.), breast cancer, cellular abnormalities, or any pathological structure in a medical or biomedical image or video.

p-0026The present invention and its attributes and advantages will be further understood and appreciated with reference to the detailed description below of presently contemplated embodiments, taken in conjunction with the accompanying Figures.

BRIEF DESCRIPTION OF THE DRAWINGS

p-0027<figref idrefs="DRAWINGS">FIG. 1</figref> illustrates a block diagram of one embodiment of AM-FM reconstruction according to the present invention;

p-0028<figref idrefs="DRAWINGS">FIG. 2</figref> illustrates a block diagram of an exemplary computer system for implementing the methods according to the present invention;

p-0029<figref idrefs="DRAWINGS">FIG. 3</figref> illustrates an input image and reconstructions using the CCA method, the DCA method, and the method according to the present invention;

p-0030<figref idrefs="DRAWINGS">FIG. 4</figref> illustrates an input image and reconstructions using the CCA method, the DCA method, and the method according to the present invention; and

p-0031<figref idrefs="DRAWINGS">FIG. 5</figref> illustrates an input image and reconstructions using the CCA method, the DCA method, and the method according to the present invention.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

p-0032<figref idrefs="DRAWINGS">FIG. 1</figref> illustrates a block diagram of AM-FM reconstruction <b>100</b> according to the present invention. Specifically, <figref idrefs="DRAWINGS">FIG. 1</figref> illustrates the solution to the minimization problem:

p-0033<maths id="MATH-US-00007" num="00007"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><mi>J</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>a</mi><mo>,</mo><mi>ζ</mi></mrow><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mrow><mfrac><mn>1</mn><mi>p</mi></mfrac><mo>⁢</mo><msubsup><mrow><mo></mo><mrow><mrow><mi>f</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>a</mi><mo>,</mo><mi>ζ</mi></mrow><mo>)</mo></mrow></mrow><mo>-</mo><mi>b</mi></mrow><mo></mo></mrow><mi>p</mi><mi>p</mi></msubsup></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>a</mi></msub><mo>⁢</mo><mrow><mi>T</mi><mo>⁡</mo><mrow><mo>(</mo><mi>a</mi><mo>)</mo></mrow></mrow></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>ζ</mi></msub><mo>⁢</mo><msub><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mn>1</mn></msub></mrow></mrow></mrow><mo>,</mo><mrow><mrow><mi>s</mi><mo>.</mo><mi>t</mi><mo>.</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>a</mi></mrow><mo>≥</mo><mn>0</mn></mrow><mo>,</mo><mrow><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mo>≤</mo><mn>1</mn></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>

p-0034An input image signal is provided at step <b>110</b>, although it is contemplated that an input video may also be provided. Enforcing two constrains in the AM-FM reconstruction, a small number of locally coherent components and a piecewise smooth constrain for the amplitude functions b(ξ) in the following equation:

p-0035<maths id="MATH-US-00008" num="00008"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>b</mi><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>=</mo><mrow><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>L</mi></munderover><mo>⁢</mo><mrow><mrow><msub><mi>a</mi><mi>n</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>⁢</mo><mrow><mi>cos</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><msub><mi>φ</mi><mi>n</mi></msub><mo>⁡</mo><mrow><mo>(</mo><mi>ξ</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> is desired.

p-0036An extended analytic signal of the input image is computed at step <b>120</b> by applying a Hilbert transform to form a 2D extension analytic signal b(ξ) of the 1D analytic signal. The extended analytic signal is processed through the filterbank and data is interpreted at step <b>130</b>.

p-0037The functional J(a,ζ) is convex in a or ζ; but is not necessarily convex in both variables together. In one embodiment of the invention, the numerical value of the local minimum may be solved by an optimization procedure shown in step <b>140</b> for which updates are alternated for each independent variable.

p-0038In one embodiment, the optimization procedure shown in step <b>140</b> is summarized by the following Equation (3) and Equation (4), for k=1, 2, . . . :

p-0039<maths id="MATH-US-00009" num="00009"><math overflow="scroll"><mtable><mtr><mtd><mrow><msup><mi>ζ</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></msup><mo>=</mo><mrow><mrow><munder><mi>min</mi><mrow><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mo>≤</mo><mn>1</mn></mrow></munder><mo>⁢</mo><mrow><mfrac><mn>1</mn><mrow><mi>p</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>ζ</mi></mrow></mfrac><mo>⁢</mo><msubsup><mrow><mo></mo><mrow><mrow><mi>f</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><msup><mi>a</mi><mrow><mo>(</mo><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow><mo>)</mo></mrow></msup><mo>,</mo><mi>ζ</mi></mrow><mo>)</mo></mrow></mrow><mo>-</mo><mi>b</mi></mrow><mo></mo></mrow><mrow><mi>P</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>ζ</mi></mrow><mrow><mi>P</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>ζ</mi></mrow></msubsup></mrow></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>ζ</mi></msub><mo>⁢</mo><msub><mrow><mo></mo><mi>ζ</mi><mo></mo></mrow><mn>1</mn></msub></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>a</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></msup><mo>=</mo><mrow><mrow><munder><mi>min</mi><mrow><mi>a</mi><mo>≥</mo><mn>0</mn></mrow></munder><mo>⁢</mo><mrow><mfrac><mn>1</mn><msub><mi>p</mi><mi>α</mi></msub></mfrac><mo>⁢</mo><msubsup><mrow><mo></mo><mrow><mrow><mi>f</mi><mo>⁡</mo><mrow><mo>(</mo><mrow><mi>a</mi><mo>,</mo><msup><mi>ζ</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></msup></mrow><mo>)</mo></mrow></mrow><mo>-</mo><mi>b</mi></mrow><mo></mo></mrow><mrow><mi>p</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>α</mi></mrow><mrow><mi>p</mi><mo>⁢</mo><mstyle><mspace width="0.3em" height="0.3ex" /></mstyle><mo>⁢</mo><mi>α</mi></mrow></msubsup></mrow></mrow><mo>+</mo><mrow><msub><mi>λ</mi><mi>a</mi></msub><mo>⁢</mo><mrow><mi>T</mi><mo>⁡</mo><mrow><mo>(</mo><mi>a</mi><mo>)</mo></mrow></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths><br /> to solve Equation (2). Specifically, ζ^(k)of Equation (3) can be solved based on the non-negative quadratic programming optimization algorithm and the FOCUSS algorithm. a^(k)of Equation (4) can be solved using the vector valued IRN-NQP (vv-IRN-NQP) algorithm which is based on the non-negative quadratic programming optimization algorithm and on the iteratively reweighted norm (“IRN”) algorithm for vector-valued images.

p-0040Solving for ζ^(k)of Equation (3) and a^(k)of Equation (4) to ultimately solve for Equation (2) is strongly dependent on the accuracy of the initial instantaneous amplitude estimates, for example, a⁽⁰⁾is the instantaneous amplitude estimate from CCA. In one embodiment of the invention, the instantaneous amplitude estimate is obtained by the Quasi-Local Method (“QLM”) and is preferred over the Quasi-Eigen Approximately (“QEA”) method since the instantaneous estimates using QLM are less sensitive to perturbations—i.e., noise—than the instantaneous amplitude estimates of QEA. After optimization at step <b>140</b>, the image is reconstructed at step <b>150</b>.

p-0041<figref idrefs="DRAWINGS">FIG. 2</figref> illustrates an exemplary computer system <b>200</b>, or network architecture, that may be used to implement the methods according to the present invention. One or more computer systems <b>200</b> may carry out the methods presented herein as computer code. One or more processors, such as processor <b>204</b>, which may be a special purpose or a general-purpose digital signal processor, is connected to a communications infrastructure <b>206</b> such as a bus or network. Computer system <b>200</b> may further include a display interface <b>202</b>, also connected to communications infrastructure <b>206</b>, which forwards information such as graphics, text, and data, from the communication infrastructure <b>206</b> or from a frame buffer (not shown) to display unit <b>230</b>. Computer system <b>200</b> also includes a main memory <b>205</b>, for example random access memory (“RAM”), read-only memory (“ROM”), mass storage device, or any combination thereof. Computer system <b>200</b> may also include a secondary memory <b>210</b> such as a hard disk drive <b>212</b>, a removable storage drive <b>214</b>, an interface <b>220</b>, or any combination thereof. Computer system <b>200</b> may also include a communications interface <b>224</b>, for example, a modem, a network interface (such as an Ethernet card), a communications port, a PCMCIA slot and card, wired or wireless systems, etc.

p-0042It is contemplated that the main memory <b>205</b>, secondary memory <b>210</b>, communications interface <b>224</b>, or a combination thereof function as a computer usable storage medium, otherwise referred to as a computer readable storage medium, to store and/or access computer software and/or instructions.

p-0043Removable storage drive <b>214</b> reads from and/or writes to a removable storage unit <b>215</b>. Removable storage drive <b>214</b> and removable storage unit <b>215</b> may indicate, respectively, a floppy disk drive, magnetic tape drive, optical disk drive, and a floppy disk, magnetic tape, optical disk, to name a few.

p-0044In alternative embodiments, secondary memory <b>210</b> may include other similar means for allowing computer programs or other instructions to be loaded into the computer system <b>200</b>, for example, an interface <b>220</b> and a removable storage unit <b>222</b>. Removable storage units <b>222</b> and interfaces <b>220</b> allow software and instructions to be transferred from the removable storage unit <b>222</b> to the computer system <b>200</b> such as a program cartridge and cartridge interface (such as that found in video game devices), a removable memory chip (such as an EPROM, or PROM) and associated socket, etc.

p-0045Communications interface <b>224</b> allows software and instructions to be transferred between the computer system <b>200</b> and external devices. Software and instructions transferred by the communications interface <b>224</b> are typically in the form of signals <b>225</b> which may be electronic, electromagnetic, optical or other signals capable of being received by the communications interface <b>224</b>. Signals <b>225</b> are provided to communications interface <b>224</b> via a communications path <b>226</b>. Communications path <b>226</b> carries signals <b>225</b> and may be implemented using wire or cable, fiber optics, a phone line, a cellular phone link, a Radio Frequency (“RF”) link or other communications channels.

p-0046Computer programs are stored in main memory <b>205</b> and/or secondary memory <b>210</b>. Computer programs may also be received via communications interface <b>224</b>. Computer programs, when executed, enable the computer system <b>200</b>, particularly the processor <b>204</b>, to implement the methods according to the present invention. The methods according to the present invention may be implemented using software stored in a computer program product and loaded into the computer system <b>200</b> using removable storage drive <b>214</b>, hard drive <b>212</b> or communications interface <b>224</b>. The software and/or computer system <b>200</b> described herein may perform any one of, or any combination of, the steps of any of the methods presented herein. It is also contemplated that the methods according to the present invention may be performed automatically, or may be invoked by some form of manual intervention.

p-0047The invention is also directed to computer products, otherwise referred to as computer program products, to provide software to the computer system <b>200</b>. Computer products store software on any computer useable medium. Such software, when executed, implements the methods according to the present invention. Embodiments of the invention employ any computer useable medium, known now or in the future. Examples of computer useable mediums include, but are not limited to, primary storage devices (e.g., any type of random access memory), secondary storage devices (e.g., hard drives, floppy disks, CD ROMS, ZIP disks, tapes, magnetic storage devices, optical storage devices, Micro-Electro-Mechanical Systems (“MEMS”), nanotechnological storage device, etc.), and communication mediums (e.g., wired and wireless communications networks, local area networks, wide area networks, intranets, etc.). It is to be appreciated that the embodiments described herein can be implemented using software, hardware, firmware, or combinations thereof.

p-0048The computer system <b>200</b>, or network architecture, of <figref idrefs="DRAWINGS">FIG. 2</figref> is provided only for purposes of illustration, such that the present invention is not limited to this specific embodiment. It is appreciated that a person skilled in the relevant art knows how to program and implement the invention using any computer system or network architecture.

p-0049The performance of the present invention in terms of image reconstruction quality was compared with that of several alternative approaches, including Channelized Component Analysis (“CCA”), Dominant Component Analysis (“DCA”), Least-Squares Reconstructions (“LESHA” and “LESHAL”) and Multi-Scale Least-Squares Reconstructions (“MULTILES”). For the CCA and DCA methods, the instantaneous amplitude (“IA”) was computed using QLM, the instantaneous phase (“IP”) was computed using the QEA, the LESHA, LESHAL and MULTILES use QEA to estimate IA and IP. As discussed more fully below, the image reconstruction quality using Regularized Optimization (“RO”) was superior.

p-0050A filterbank covering the whole frequency spectrum consisting of one low-pass and one high-pass filter is used. Each separable channel filter has support over four quadrants. To maintain support over only two quadrants needed for the QEA method, Fast Fourier Transform (“FFT”) pre-filtering is used to remove support in two quadrants, for example, the two left quadrants or two right quadrants. Thus, each bandpass filter has frequency support in only two quadrants of the frequency spectrum so that, in effect, each channel filter operates over a single quadrant. The filters are designed using a min-max, equiripple approach. In a preferred embodiment, passband ripple is set at 0.017 dB and the stopband attenuation is set to 66.02 dB. Because the filterbank covers the entire spectrum, it can be expected that the instantaneous frequency will fall within the spectral support of one of the channel filters. It is assumed that local image coherency will force the instantaneous frequency estimate to fall within the passband of the dominant bandpass filter.

p-0051<figref idrefs="DRAWINGS">FIG. 3</figref>, <figref idrefs="DRAWINGS">FIG. 4</figref> and <figref idrefs="DRAWINGS">FIG. 5</figref> illustrate original input images and reconstructions using the CCA method, the DCA method, and the method according to the present invention. Specifically, <figref idrefs="DRAWINGS">FIG. 3</figref> illustrates the “Radial Chirp” image, <figref idrefs="DRAWINGS">FIG. 4</figref> illustrates the “Barbara” image and <figref idrefs="DRAWINGS">FIG. 4</figref> illustrates the “Lena” image. All images contain 512×512 pixels and the simulation is carried out using Matlab-only code on a 1.83 GHz Intel Dual core CPU with a 4G RAM.

p-0052Table 1 below is a comparison of the signal to noise ratio of the reconstructed images using DCA, CCA, LESHA, LESHAL, MULTILES and the RO method (Equation (2)+DCA) according to the present invention. Values with ( )* indicates that the reconstructed image has been normalized. For all cases, the RO approach has superior performance, especially for the gray scale photograph images Barbara and Lena (see <figref idrefs="DRAWINGS">FIG. 4</figref> and <figref idrefs="DRAWINGS">FIG. 5</figref>). The signal-to-noise ratio is greater than 22 dB, compared with the modest signal-to-noise ratio of less than 15.2 than of all other methods for all three images tested. In no event is the signal-to-noise ratio less than 15.

p-0053<tables id="TABLE-US-00001" num="00001"><table frame="none" colsep="0" rowsep="0" pgwide="1"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="1" colwidth="28pt" align="left" /><colspec colname="2" colwidth="238pt" align="center" /><thead><row><entry namest="1" nameend="2" rowsep="1">TABLE 1</entry></row></thead><tbody valign="top"><row><entry namest="1" nameend="2" align="center" rowsep="1" /></row><row><entry /><entry>SNR (dB)</entry></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="7"><colspec colname="1" colwidth="28pt" align="left" /><colspec colname="2" colwidth="49pt" align="center" /><colspec colname="3" colwidth="49pt" align="center" /><colspec colname="4" colwidth="28pt" align="center" /><colspec colname="5" colwidth="35pt" align="center" /><colspec colname="6" colwidth="42pt" align="center" /><colspec colname="7" colwidth="35pt" align="center" /><tbody valign="top"><row><entry>Image</entry><entry>DCA</entry><entry>CCA</entry><entry>LESHA</entry><entry>LESHAL</entry><entry>MULTILES</entry><entry>(2) + DCA</entry></row><row><entry namest="1" nameend="7" align="center" rowsep="1" /></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="9"><colspec colname="1" colwidth="28pt" align="left" /><colspec colname="2" colwidth="21pt" align="right" /><colspec colname="3" colwidth="28pt" align="left" /><colspec colname="4" colwidth="21pt" align="right" /><colspec colname="5" colwidth="28pt" align="left" /><colspec colname="6" colwidth="28pt" align="char" char="." /><colspec colname="7" colwidth="35pt" align="center" /><colspec colname="8" colwidth="42pt" align="center" /><colspec colname="9" colwidth="35pt" align="center" /><tbody valign="top"><row><entry>Radial</entry><entry>6.51</entry><entry>(14.43)*</entry><entry>3.21</entry><entry>(−2.63)*</entry><entry>0.51</entry><entry>13.56</entry><entry>13.56</entry><entry>15.41</entry></row><row><entry>Chirp</entry><entry /><entry /><entry /><entry /><entry /><entry /><entry /><entry /></row><row><entry>Barbara</entry><entry>0.92</entry><entry>(7.69)*</entry><entry>1.15</entry><entry>(9.76)*</entry><entry>10.48</entry><entry>12.69</entry><entry>12.69</entry><entry>22.39</entry></row><row><entry>Lena</entry><entry>0.83</entry><entry>(5.38)*</entry><entry>0.46</entry><entry>(6.58)*</entry><entry>14.97</entry><entry>15.16</entry><entry>15.16</entry><entry>24.40</entry></row><row><entry namest="1" nameend="9" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0054Table 2 is a comparison of the same reconstructed images using a structured similarity index (“SSIM”). SSIM measures the visual structural similarity between the reconstructed images with the original reference image. Again, the RO approach (Equation (2)+DCA) offers the highest SSIM rating, above 0.95 for all three pictures.

p-0055<tables id="TABLE-US-00002" num="00002"><table frame="none" colsep="0" rowsep="0"><tgroup align="left" colsep="0" rowsep="0" cols="2"><colspec colname="1" colwidth="35pt" align="left" /><colspec colname="2" colwidth="182pt" align="center" /><thead><row><entry namest="1" nameend="2" rowsep="1">TABLE 2</entry></row></thead><tbody valign="top"><row><entry namest="1" nameend="2" align="center" rowsep="1" /></row><row><entry /><entry>SSIM index [15]</entry></row></tbody></tgroup><tgroup align="left" colsep="0" rowsep="0" cols="7"><colspec colname="1" colwidth="35pt" align="left" /><colspec colname="2" colwidth="21pt" align="center" /><colspec colname="3" colwidth="21pt" align="center" /><colspec colname="4" colwidth="28pt" align="center" /><colspec colname="5" colwidth="35pt" align="center" /><colspec colname="6" colwidth="42pt" align="center" /><colspec colname="7" colwidth="35pt" align="center" /><tbody valign="top"><row><entry>Image</entry><entry>DCA</entry><entry>CCA</entry><entry>LESHA</entry><entry>LESHAL</entry><entry>MULTILES</entry><entry>(2) + DCA</entry></row><row><entry namest="1" nameend="7" align="center" rowsep="1" /></row><row><entry>Radial</entry><entry>0.929 </entry><entry>0.766</entry><entry>0.144</entry><entry>0.815</entry><entry>0.815</entry><entry>0.978</entry></row><row><entry>Chirp</entry><entry /><entry /><entry /><entry /><entry /><entry /></row><row><entry>Barbara </entry><entry>0.730</entry><entry>0.544</entry><entry>0.619</entry><entry>0.656</entry><entry>0.656</entry><entry>0.987</entry></row><row><entry>Lena</entry><entry>0.623</entry><entry>0.378</entry><entry>0.731</entry><entry>0.731</entry><entry>0.731</entry><entry>0.962</entry></row><row><entry namest="1" nameend="7" align="center" rowsep="1" /></row></tbody></tgroup></table></tables>

p-0056The visual quality of the reconstructed images is consistent with the quantitative measurements. <figref idrefs="DRAWINGS">FIG. 3</figref> is a side-by-side comparison of the original “Radial Chirp” image (<figref idrefs="DRAWINGS">FIG. 3A</figref>) with images reconstructed using the CCA method (<figref idrefs="DRAWINGS">FIG. 3B</figref>), the DCA method (<figref idrefs="DRAWINGS">FIG. 3C</figref>), and the RO method (<figref idrefs="DRAWINGS">FIG. 3D</figref>), respectively. Both the CCA and DCA methods create strong circular artifacts around the four edges. The contrast of the CCA reconstructed image is noticeably muted. The DCA method introduces discontinuous patches into the original continuous image. The RO method produced an image without visible artifacts, while maintaining a level of contrast and continuity comparable to the original image.

p-0057The same improved quality holds true for photographic images. <figref idrefs="DRAWINGS">FIG. 4</figref> is a side-by-side comparison of the original “Barbara” image (<figref idrefs="DRAWINGS">FIG. 4A</figref>) with images reconstructed using the CCA method (<figref idrefs="DRAWINGS">FIG. 4B</figref>), the DCA method (<figref idrefs="DRAWINGS">FIG. 4C</figref>), and the RO method (<figref idrefs="DRAWINGS">FIG. 4D</figref>), respectively. Both the CCA and the DCA methods create visibly blurry reconstructions with artifacts of lines crossing the image at different angles. Both images lost the texture detail of the checkered tablecloth and the creases on the pants. The RO reconstruction produces a clean image with the appropriate texture without visible artifacts. A similar result is achieved in the “Lena” image shown in <figref idrefs="DRAWINGS">FIG. 5</figref>.

p-0058While the disclosure is susceptible to various modifications and alternative forms, specific exemplary embodiments thereof have been shown by way of example in the drawings and have herein been described in detail. It should be understood, however, that there is no intent to limit the disclosure to the particular embodiments disclosed, but on the contrary, the intention is to cover all modifications, equivalents, and alternatives falling within the scope of the disclosure as defined by the appended claims.