Super–Resolution Restoration and Image Reconstruction for Passive Millimeter Wave Imaging

According to blackbody radiation theory [1], all substances at a finite absolute temperature will radiate electromagnetic energy. Passive millimeter-wave (PMMW) imaging system forms images by detecting the millimeter-wave radiation energy from the scene and utilizing the differences of the radiation intensity[2,3]. Although such imaging has been performed for decades (or more, if one includes microwave radiometric imaging), new sensor technology in the millimeter-wave regime has enabled the generation of PMMW imaging at video rates and has renewed interest in this area. Clouds and fog are effectively transparent to millimeter-wave and the cold sky is reflected by metallic objects on the ground making PMMW images similar to infrared (IR) and visible images. Due to being able to perform well under adverse weather conditions, PMMW imaging offers advantages over IR and visible imaging. It is widely used in airport security, scene monitoring, plane blind landing, medical diagnosis and environmental detection, et al [4-8].

However, the obtained images usually have the inherent problem of poor resolution, which is caused by limited aperture dimensions and the consequent diffraction limit. Images acquired from practical sensing operations usually suffer from poor resolution due to the finite size limitations of the antenna, or the lens, and the consequent imposition of diffraction limits. The fundamental operation underlying the sensing operation is the "low-pass" filtering effect due to the finite size of the antenna lens. The image at the output of the imaging system is a lowpass filtered version of the original scene. There is no useful signal beyond the cut-off frequency in the measured data, and the information lost by the imaging system are the fine details, i.e. high-frequency spectral components. In order to restore the details and improve the resolution of the image, some methods of image restoration will be needed. As is well-known that the problem of image restoration is a inverse problem and inverse problem is always singular or ill-posed. Traditional image restoration methods based on de-convolution approaches principally try to restore the information of the pass band and eliminate the effect of additive noise components. Therefore, these methods have merely limited resolution enhancement capabilities. Greater resolution improvements can only be achieved through a class of more sophisticated algorithms, called super-resolution algorithmsor image reconstruction algorithms, before the PMMW imaging can be employed.Some studies have www.intechopen.com Image Restoration -Recent Advances and Applications 260 indicated that the cost of an imager increases as (1/Resolution) raised to the power 2.5. Hence, a possible two-fold improvement in resolution by super resolution processing, roughly translates into a cost reduction of an imager by more than 5 times.
Super-resolution algorithms can be classified into iterative and non-iterative algorithms. Iterative algorithms are generally the preferred approach due to their numerous advantages and also since the iteration can be terminated once a solution of a reasonable quality is achieved. Non-iterative algorithms nvolves convolution operations in the spatial domain, direct inverse methods, regularizedpseudo inverse techniques [9].The existing iterative super-resolution methods include , MAP method [12], steepest descent, conjugate gradients [13] and projection on convex set (POCS) method [14].
This chapter considers the general problem of super-resolution restoration and image reconstruction. While our focus will be on application to PMMW imaging. This chapter is based on work presented in [15], portions of which appeared in [13,[15][16][17][18][19]. To solve the inherent problem of poor resolution which is caused by limited aperture dimensions and the consequent diffraction limit, this chapter presents system model, analyses the theoretical research results and design specifically for PMMW imaging. Firstly, we estimate the PSF of the PMMW imaging system by a variational Bayesian blind restoration algorithm. Secondly, we focus on mainly four algorithms, including Conjugate-Gradient algorithm (CG), Adaptive Projected Landweber super-resolution algorithm(APL), and Undedicated Steerable Pyramid Transform Projected Landweber algorithm (USPTPL), and Two-step Projection Iteration Thresholding (tw-PIT) supper resolution using compressed sensing architecture. Finally, we have verified the system model and super-resolution algorithms by experiment in different plane using differentsystem.

PMMW image formulation model
Passive millimeter-wave (PMMW) imaging is a method of forming images through the passive detection ofnaturally occurring millimeter-wave radiation from a scene.According toblackbody radiation theory [1], all substances with a finite absolute temperature will radiate electromagnetic energy. The radiated energy spectral intensitycan be described as a Brightness Temperature B (J.K-1) is Boltzmann's Constant; c is velocity of light and the T represents the absolute temperature. The brightness temperature can be assigned a gray level to generate a millimeter wave image. For the convenience of analyzing, we introduce a simplified focal plane model to calculate the imaging process, as illustrated in Fig.1 and Fig.2. The noise power can be received by sensor in the data Plane or in the image plane. Thus the millimeter wave image can be formed in the data plane or in the image plane. Fig.1 shows the imaging process of PMMW Focal Plane Array (FPA) imaging. The process is Space-variant by the non-uniformity of the antenna beam and inconsistencies of channels. Fig.2 shows the principle of imaging in different planes. System will obtain the higher www.intechopen.com Super-Resolution Restoration and Image Reconstruction for Passive Millimeter Wave Imaging 261 resolution and rate near to video when lens or reflectors are used in focal plane real-time imaging.
In Fig.2, the propagation process from object-plane to data-plane is two-dimensional (2D) spatial Fourier transformation (FT), while the propagation process from data-plane to focalplane is 2D spatial inverse Fourier transformation (IFT). FPA system forms images in the image-plane by receiving object-radiation energy, without 2D FT and IFT in post-data processing. Thus, the system is high real-time.
According to Fresnel and Fraunhofer diffraction theory [20], the incoherent image is given by nuv is the noise of image plane,   S  denotes the nonuniformity of the antenna beams and inconsistencies of channels.
Where 0 S is the gain consistency of the receiver-channels. This model only considers the non-uniform effects of array antenna beams.
The imaging system is linear space-invariant when the non-uniformity of antenna beams and the inconsistencies of channels can be ignored together. Fig.2 can be shown as a simplified convolution model Where the function (,) huv is the PSF of the imaging, which determines the radiant energy distribution in the image plane from a point source of radiant energy located in the object plane.
In matrix and vector form, we can write the problem as gH fn   Where Gxy , (,) Fxy , (,) Nxy respectively.
The so-called image restoration problem is solving for (,) fxy based on the above mathematic model, namely the inverse problem of equation (5). And the inverse problem is usually bizarre or morbid. According to equation (6), there is: From the equation (6) we know, when using the wrong or inaccurate PSF, the solution of ˆ( ,) Fuv would be wrong. And under the effect of system noise, using the inaccurate PSF can seriously influence the imaging quality. So on condition that it's unable to accurately determine the PSF of the imaging system, the way which according to the obtained blurred image, using blind processing method to estimate more accurate PSF, and then utilizing the classic image restoration algorithm to recover the image, can improve the passive millimeter wave imaging quality.
As theˆ(,) Huv is zero outside the cut-off frequency of the imaging system. Only from the frequency domain point of view, to restore high-frequency components outside the cutoff frequency appears to be impossible.
The analytical continuation theory is the theoretical basis for to achieve super-resolution. Analytic continuation theory includes two aspects, 1) The Fourier transform of any airspace bounded function is analytic functions. 2) For any analytic function, as long as we can accurately know the part information in a limited range , we can uniquely determine the entire function. Given the values of analytic function in a range, the overall reconstruction of the function is called analytic continuation [21].

Variational bayesian estimate the PSF of the imaging system
In passive millimeter wave imaging signal processing, using accurate point spread function (PSF) is important for getting high quality restored image. In the process of imaging, the PSF is decided by antenna beam, atmospheric transmission, system noise etc. So the real PSF can not be substituted by a simple model. We can't use a simple model to replace the real PSF.
And the PSF will change due to the specific imaging environment so that it is difficult to acquire the exact PSF. The common methods to obtain the PSF include direct measurement and model parameter estimation. In order to improve the quality of image restoration, a variational Bayesian blind restoration algorithm for PMMW imaging is proposed in this section, which combines the posterior probability model of the PMMW imaging to obtain accurate point spread function by variational Bayesian estimation The variational Bayesian estimation algorithm is based on Bayesian framework [22]. The idea of this algorithm is according to the known priori information and assumptions establish a posteriori distribution model of the imaging system first. And then under a certain criterion, use the variational [23,24] method to obtain a optimal distribution which approximates with the posteriori distribution, so that reducing the complexity of the model and make the problem more analytic. Finally work out the estimation under the restriction of the cost function.
According to the Bayesian theorem, the posteriori distribution of passive millimeter wave imaging system is decided by the prior distribution of original image, the prior distribution of PSF, noise probability distribution and likelihood probability distribution. Because of the prior distribution of original image is established on its gradient value of the statistical distribution, according to the linear relationship of the above imaging model we can get is the gradient value of each point. The prior distribution of original image can be expressed by c-th order zero mean mixture gaussian distribution, the mathematical expressed as 1 ,, Where c v is the variance of the c-th Gaussian distribution, c  is the weight of the c-th Gaussian distribution.
Also, assuming that every point value of PSF is independent , Combine the equation (5), (6), (7), (8), the posteriori distribution of the passive millimeter wave imaging system is 2 2 1 ,, In order to minimize KL D , the minimum value can be deserved by partial derivate the equation (13) Obtain the parametric expression of () QH from the form of the known () PH . Factor () QH into KL D , and get the update equations.Use the gradient descent method to get the ( ) opt QH .Get the expected value of ( ) opt QH so that get the estimation of the PSF.
For PMMW imaging, the power pattern of the antenna is similar to gaussian function, so at the beginning of the iteration we could use the gaussian function as its initial value, and then through variational bayesian estimation to get more accurate PSF. That can reduce the iteration times, and estimate the PSF faster.
According to the above analysis and derivation, the process of variational bayesian blind restoration algorithm as shown in figure 3. First input blurred image and compute its gradient value. Second set the initial value of variational bayesian iteration algorithm. Then begin the variational bayesian estimation through the known priori information, when the K-L distance less than the threshold, stop the iteration and output the estimated PSF. At last use the Lucy-Richardson algorithm to restore the image, get the millimeter image of the scene.   . Let the experimental PSF convolute with the original image to get the blurred image which is showed in fig4-c. Utilize the variational bayesian estimate algorithm to get the PSF, the result is showed in fig4-d. Compare the estimated PSF with experimental one, we can find that the shapes are roughly same, only a few details are different.
In order to illustrate the correct of the variational bayesian estimation, we restore the image through the Lucy-Richardson algorithm using experimental PSF and estimated PSF separately. Each restoration process iterates 20 times, the recovery results are showed in fig5.

Super-Resolution Restoration and Image Reconstruction for Passive Millimeter Wave Imaging 267
The image restoration effect is assessed by sum of squared differences (SSD) between recovery image and original scene, the formula is as follows: M, N are two images. The SSD error is smaller the recovery image is more approximate to the original one, namely the restoration effect is better. The SSD error of fig3-a is 179, the SSD error of fig3-b is 184. We can find that the SSD errors of the images which restored by two kind of PSF are proximate, so verifying the correct of the estimated PSF.  The PSF that deserved through variational bayesian estimation, showed in fig7-a, and then use the Lucy-Richardson algorithm through 20 iteration to get the recovery image, showed in fig7-b. And for PMMW image processing, we can also get the estimated PSF through the parameter model method that according to the diffractiFon cut-off characteristics of the image system, utilizing the spectrum of the blurred image to get the parameters of the gaussian-form PSF. The result showed in fig7-c. Then uses the same restoration algorithm to get the recovery image.

Fig. 7. The restoration of two methods
According to the results of the above experiment, the PSF estimated by variational bayesian method is more approximate to the antenna power pattern, the recovery image is more clearer. Also compute the SSD error between recovery image and original image. For fig7-b the SSD error is 878, and for fig7-d the SSD error is 1088. We can find that the SSD error is smaller and the recovery image is more approximate to the original scene when using the variational bayesian blind restoration algorithm rather than the parameter model method.
The third experiment uses the measured data that from PMMW imaging system to verify the algorithms availability. We use a single-channel scanning radiometer for imaging. Due to the limitation of scanning angle in pitch direction, the PMMW image isstitched together.Fig8-a is the optical image of the scene, fig8-b is the obtained millimeter wave image which is blurred, fig8-c is the estimated PSF, and fig8-d is the recovery image from variational bayesian blind restoration algorithm. We can find that the recovery image's contour and details are clearer, the imaging quality is improved effectively.

Super-resolution restoration and image reconstruction
In order to improve the resolution of the image, some methods of image restoration will be needed. The image restoration is an inverse problem in general, which is always ill posed. Traditional de-convolution approaches restore the information of the pass band and eliminate the effect of additive noise components. Therefore, these methods have only limited resolution enhancement capabilities. Greater resolution improvements can only be achieved through a class of more sophisticated algorithms, called super-resolution algorithms, including Lucy-Richardson algorithm, Conjugate-gradient (CG) algorithm,Adaptive Projected Landweber (APL) super-resolution algorithm, Undecimated Steerable Pyramid Transform Projected Landweber (USPTPL) algorithm and Two-step Projection Iteration Thresholding (tw-PIT) supper resolution using compressed sensing architecture algorithm. In this section,wepresentthe image restoration algorithms.

Fig. 8. Measured data restoration
Stochastic super-resolution image restoration, typically a Bayesian approach, provides a flexible and convenient way to model a priori knowledge concerning the solution. Bayesian estimation methods are used when the a posteriori probability density function (PDF) of the original image can be established. Using maximum-a-posteriori estimation, we can obtain an exact solution of the formula.

Conjugate-gradient algorithm
The Conjugate-gradient (CG) algorithm is an effective Krylov subspace method of solving an unconstrained large-scale optimization problem, which is equivalent to solve the quadratic minimization problem [25,26].
The target function can be denoted as And the gradient function is The k f is the k-th iteration estimate of the original scene, it generates a descent direction k d , n~ is conjugate to all previous search directions: 12 , ...... n dd dwith respect to matrix T HH ; that is TT nk dHH d, k.=0,1,…n-1. In other words, conjugate gradient algorithm deals with problem of 1D searching: where k  denotes optimal size at searching directions, which is decided by Because the standard conjugate gradient algorithm is a linear restoration method, it has only limited capability of super-resolution, that is, spectral extrapolation. Furthermore, during each iteration, we can not guarantee the nonnegative constraint of estimated images. An interesting feature of the CG method is that it can be modified in order to take into account the additional priori-information about the solution. The information can be expressed as a number of closed convex sets.
An example of a constraint which is rather natural in many problems of image restoration is the positivity of the solution. The constraint can also be projected as a closed convex set. Thus we can impose the constraint that the image is nonnegative on the iteration. Because the projection operation is non-linear operations, it introduces frequencies beyond the passband. Thus, the modified CG algorithm has the capacity of spectrum extrapolation. It can be shown as where c P denotes the projection operator on the constraint set C . Then the projection operator is given by The relationship between W, W -1 and Discrete Fourier Transform (DFT) is as follows: In image processing, n is selected in the 2-power, such as 28-2128, etc. The amount of calculation of image direct iteration restoring is very great at this time. However, H is Toeplitz matrix(BTTB), so the matrix H (BTTB) can be continued to circulant matrix structure (BCCB) by the above mentioned matrix-vector computing and Fast Fourier Transform (FFT). Thus, the CG algorithm has quick convergence.

Simulation and experiments
To verify the super-resolution capability of the CG algorithm, representative results of restoring blurred images is shown in Figs. 10a-c. The simulation image is a 256 256  synthetic image composed of a series of concentric disks with the background (dark) at intensity value zero and the disks (bright) at intensity value one.
The simulation image and its spectrum are shown in Fig.10a and Fig.11a, respectively. To simulate the blurring caused by a diffraction limited imaging sensor, we blur the ideal image by the PSF of a circular aperture antenna which is equivalent to a low-pass filter. Zero-mean Gaussian noise is added to the blurred image to get the noisy blurred image. The blurred image and its spectrum are shown in Fig.4b Fig.9 and Fig.10, we can see that the super resolution capability (spectralextrapolation) of the CG algorithm is improvement remarkable. Furthermore, the CG algorithm can effectively reduce rings effects.
To evaluate the capacity of the recovery algorithm, the simplest common assessment criteria is to calculate the 2 L norms of the deviation between the original image and restore image, that is the Mean Square Error (MSE), which is expressed as where f denotes the ideal image, k f  represents the k-th restoration result, 2  is 2 L norms. Fig. 11 plots the MSE of the CG algorithms versus iteration numbers. Obviously, the CG algorithm has a fast convergence of the MSE. Generally, the MSE will be below 0.065 when the iterations number is 6. Clearly, the CG algorithm has the good convergence performance.
In the second experiment, a small number of 32x32 gun images are collected by 91.5GHz mechanically scanned mono-channel radiometer with horn antenna in the lab and around surroundings, as shown in Fig.12. The captured PMMW image is shown in Fig.12b. The restored image is shown in Fig.12c. From Fig.12, we can see that the super-resolution performance and spatial resolution are enhanced by the algorithm. Experiment results demonstrate that the CG super-resolution algorithm has fast convergent rate and the good spectral-extrapolation capacity. The CG algorithm improves the spatial resolution and reduce the ringing effects which are caused by regularizing the image restoration. Thus, the CG algorithm can be used in image restoration and PMMWI to enhance the super-resolution performance and eliminate most of the effects of blurring.

Adaptive projected Landweber super-resolution algorithm
It is well-known that the problem of image restoration is the computation of f , given the image data g and the PSF h . Some iterative methods have been introduced to solve equation (5) [28]. Thebasic feature is that the number of iterations plays the role of a regularization parameter because semi-convergence holds true in the case of noisy images. Iterative methods of image restoration have an advantage over one-step methods in that the partial solution may be examined at each step of the iteration and any constraints on the solution can be enforced at that time.

Landweber algorithm
The Landweber method (successive-approximation method) [29] is the simplest iterative regularizing algorithm to solve linear ill-posed problems. Because standard Landweber algorithm is a linear restoration method, it has only limited capability of super-resolution (spectral extrapolation). Furthermore, during each iteration, we cannot guarantee the nonnegative constraint of estimated image. Other disadvantages of the Standard Landweber method are slow convergence and the difficulty of choosing proper update parameter.
Landweber algorithm is the simplest iterative method for approximating the least-square solutions of equation (5). It is characterized by the equation Where T H is the transpose of the blurring operator H ; superscript k denotes the th k iteration;  is a relaxation parameter controlling the convergence, in order to guarantee the convergence of k f , the value of  are given by Where 1  is the largest singular value of matrix H [30]. The initial guess 0 f is usually set to 0.

Projected Landweber super-resolution algorithm
Because standard Landweber algorithm is a linear restoration method, it has only limited capability of super-resolution, that is, spectral extrapolation. Furthermore, during each iteration, we cannot guarantee the nonnegative constraint of estimated image.
An interesting feature of the Landweber method is that it can be modified in order to take into account additional priori information about the solution. Fundamental of reliable estimation of high frequencies is the utilization of a priori known information during the processing iteration. In fact, since image restoration is inherently an ill-posed problem, the quality of restoration and the extent for achievable super-resolution depend on the accuracy and the amount of a priori information. As shown there, many physical constraints on the unknown object can be expressed by requiring that it belong to some given closed and convex sets. While efforts at the modeling of constraint sets and the use of these in projection-based set theoretic image recovery constitutes a popular direction for current research, it seems that the idea of combining the strong point of Landweber schemes and that of the projection-based methods has not been paid much attention to. The modified Landweber method, also called the Projected Landweber Algorithm[31] is as follows Where C P denotes the projection operator onto the constraint set C .
An example of a constraint which is rather natural in many problems of image restoration is the positivity of the solution. We can impose the constraint that the image is nonnegative on the iteration. Because this constraint is non-linear operation, it introduces frequencies beyond the passband. The super-resolution performance of the Projected Landweber algorithm can be verified from subsequent experiments. Then the projection operator is given by Other convex and closed set C include finite support constraint, band limited constraint and spatial limited constraint [27].

Adaptive projected Landweber Super-resolution alogithm
One disadvantage of the Landweber method is slow convergence and the difficulty to choose proper update parameter. If  is too large, the iterative process may diverge. If  is too small, the iterative process would be slow. Therefore, it is necessary to determine a suitable  for the iterative algorithm. Lie Liang and Yuanchang Xu proposed modification of this method, Adaptive Landweber method [32], in which constant  is calculated adaptively in each iteration. The Adaptive Landweber algorithm has a better result and faster convergence than standard Landweber algorithm. Instead of using a constant update parameter, the Adaptive Landweber method computes the update parameter at each iteration and chooses the maximum of the computed parameter and the preset constantparameter to use in the next iteration.
Then we proposed a hybrid algorithm that attempt to combine the strong points of both Projected Landweber scheme (simplicity of execution, Super-resolution, etc.) and the Adaptive adjustments relax parameter  (faster convergence rate, lower mean square error, etc.). For a brief description, each cycle of this "Adaptive Projected Landweber algorithm" consists of executing the three steps: Step 1. Implement standard Landweber algorithm with initial relax parameter 0  and Step 2. Implement projection onto the convex set C Step 3. Implement relax parameter 1 k   updating algorithm Where k f  denotes the first order derivative of k f . The algorithm hence adaptively updates the relax parameter  to speed up convergence and obtain improved object estimation after each cycle of implementation. A flow-chart depicting the various step is shown in Fig. 13.

Simulation and experiments
In order to verify the super-resolution capability of this algorithm mentioned above, two images have been adopted in simulation experiments, one is a 256x256 synthetic image composed of a series of concentric disks, and the other is an 32x32 gun image captured by 91.5 GHz mechanically scanned mono channel radiometer.
In the first experiment, the ideal image and its spectrum are shown in Fig.14a and Fig.16a, respectively. For simulating the blurring caused by a diffraction limited imaging sensor, we blur the ideal image by the PSF of a low-pass filter that simulates a sensor with a circular aperture of diameter 8 pixels. The blurred image and its spectrum are shown in Fig.14b and Fig.15b respectively. the simulation results of these algorithms are shown Fig.14c-14e, and their spectrum are shown in Fig.15c-15e, respectively. It is clear that from Fig.14 and Fig.15 the results obtained by the Adaptive Porjected Landweber algorithm are better than the results obtained by the Porjected Landweber algorithm and the standard Landweber algorithm. The super-resolution capabilities (spectrum extrapolation) of these projectionbased methods are obvious. However, the convergence of the standard Landweber method is slow compared to the Adaptive Porjected Landweber algorithm and the Porjected Landweber algorithm, mainly because the standard Landweber algorithm is a linear method, and it has hardly super-resolution capability. In the second experiment, the visual image is shown in Fig.17a. The blurred image is shown in Fig.17b acquired by the PMMW radiometer. The restored images are shown in Fig.17c-17e respectively. The Adaptive Projected Landweber algorithm gives better result than standard Landweber algorithm. The results obtained by the Adaptive Projected Landweber algorithm are very similar to (slightly better than) the results obtained by the Projected Landweber algorithm. Furthermore, Gibbs rings of standard Landweber algorithm aggravate as the iteration increases. The APL algorithm, which iteratively applies a cycle of Projected Landweber algorithm followed by a relax parameter adaptive adjustment, combines the strong points of the two approaches and hence possesses a number of implementation benefits. The adaptive update parameter aims to emphasize speed and stability. Experiment results demonstrate that the results of Adaptive Projected Landweber algorithm are better than those of standard Landweber algorithm and Projected Landweber algorithm. The Adaptive Projected Landweber Super-resolution algorithm has lower MSE and produces sharper images. These constraints and adaptive character speed up the convergence of the Landweber estimation process. The superior restoration of the object features observed in the image domain and the significant extrapolation of spatial frequencies observed in the spectral domain lead to the conclusion that the Adaptive Prejected Landweber super-resolution algorithm can be used for restoration and super-resolution processing for PMMW imaging.

Undecimated steerable pyramid transform projected Landweber algorithm
Super resolution algorithms have two tasks: restoring the spectrum components within the passband (by reversing the effects of convolution with the point spread function of imaging system) and re-create the lost frequencies due to the imposition of sensor diffraction limits by spectral extrapolation. Recently, a more sophisticated spectrum decomposition technique is adopted. Using multi-scale technique, an image is decomposed into a hierarchical manner where each level corresponds to a reduced-resolution approximation of the image. It is equivalent to a filter bank that decomposes an image into different frequency components. By such decomposition, one can restore the passband firstly and then extrapolate high frequency components stage by stage. The most commonly used multi-scale methods are based on the Pyramid transform (such as Laplacian Pyramid, steerable Pyramid), the wavelets transform, the contourlets transform and so on [33][34][35][36]. The multi-scale transforms mentioned above belong to linear transform, but they are not shift invariant due to the down sampling. The lack of shift-invariance is a problem for many applications such as image restoration and image denoising because it causes pseudo-Gibbs phenomena around singularities [33]. Undecimated multi-scale methods can avoid such problem thus it has been introduced in several studies [34][35][36].

Undecimated pyramid transform
The Undecimated Pyramid transform (UPT) uses filter bank ( 10 H , 11 H ) in the analysis part and ( 10 G , 11 G ) in the synthesis part. The ideal frequency response of the building block of the UPT is shown in Fig. 18(a). The perfect reconstruction condition is given as 10 10 11 11 1 HG HG   Filters of the UPT do not need to be orthogonal or bi-orthogonal and this lack of the need for orthogonality or bi-orthogonality is beneficial for design freedom. To perform the multi-scale decomposition, we construct non-subsampled pyramids by iterated non-subsampled filter banks. For the next level, we upsample all filters by 2 in both dimensions. Therefore, they also satisfy the perfect reconstruction condition. The cascading of the analysis part is shown in Fig. 18 In above mentionProjection Landweber algorithm, essentially the property of semiconvergence indicates that the Landweber algorithm is a regularization method. Suppose a perturbed linear equation is given as below Where A is an mn  system matrix, which describes the system geometry, x is an 1 n  vector of the image pixels and b  is an 1 m  vector of the measured data with perturbation  . Let k x  be the solution of equation (35)  Known from equation (36), the Landweber stops more quickly if it is used to restore an image with a bigger  . In the restoration of a badly contaminated image, the Landweber may have not plenty iterative steps to restore or re-create the frequency components that we need. Unfortunately, signal-to-noise ratio of PMMW image is quite low due to limited integral time and bandwidth. Because of this, the PL algorithms can not provide a satisfied resolution improvement in most practical applications. The practical performance of the PL algorithm is shown in Fig. 23(c).
The denoising technique can avoid the fast termination of the Landweber algorithm. Because PSF of a PMMW imaging system is approximately band-limited, a low passband filter can perform the function. By such pre-procession, the high-frequency component of  is attenuated thus the Landweber algorithm avoid the adverse effect of the amplification of high-frequency noise.

Undecimated Pyramid Projected Landweber (UPPL) super-resolution algorithm
The basic strategy of image super-resolution based on the UPT is to use a band-limit frequency selection rule to construct a multilevel Pyramid-like restoration model from the Pyramid representations of the original data. The UPT of an image can be described as collection of low-or band-pass copies of an original image. For a band-limited image, the decomposition of UPT can implement the function of image denoising.
The framework of the UPPL algorithm is shown in Fig. 19. First the UPPL algorithm decomposes a PMMW images by the UPT. Sub images   12 1 ,,, j yy y   acquired by the decomposition correspond to the frequency bands shown in Fig. 18(b). Then the UPPL restores these frequency components from the lowest frequency band to the highest. Because the UPT has the capability to attenuate noise in higher frequency bands, the UPPL is able to improve the restoration of lower frequency bands by providing the PL algorithm plenty iterate steps. In each stage the initial guess of the current PL iterations is the restored result of the last stage. So the restoration of current frequency band is always based on a sufficient restoration of lower frequency components. The UPPL algorithm is summarized as below: Step 1. Compute the UPT of the input image for J levels   12 1 ,,, J yy y   .
Step 3. Implement j K PL iterations with initial update parameter j Step 4. Repeat step 3 to the highest scale of 1 J  .
In each stage, the UPPL restores a frequency band by an independent PL algorithm, which has its own update parameter. Known from equation (36), a big update parameter accelerates the amplification of the data error || || kk xx   thus it decreases the number of iterative steps. Because of this, the UPPL enjoys the flexibility of controlling the convergence speed of different frequency bands. The principle for choosing a suitable j  is that the number of iterative steps should be sufficient but not excessive. Commonly the concrete value of j  is acquired by a number of experiments.

Simulation and experiments
In comparison with PL algorithm, two images have been adopted in simulation experiments, one is a 256x256 synthetic image composed of a series of disks, and the other is an 32x32 gun image captured by 91.5 GHz mechanically scanned mono channel radiometer. In the first experiment, the ideal image and its spectrum are shown in Fig. 20(a) and Fig.  20(b), respectively. For simulating the blurring caused by a diffraction limited imaging sensor, we blur the ideal image by the PSF of a low-pass filter that simulates a sensor with a circular aperture of diameter 16 pixels. Zero-mean Gaussian noise was added to the blurred image to get the observed noisy blurred image at 30 dB Blurred Signal-to-Noise Ratio (BSNR). The blurred image and its spectrum are shown in Fig. 20(c) and Fig. 20(d), respectively.
In the experiment, scale number of UPT is 3. The update parameters are 1 1.5   , 2 1.2   , 3 1.0   and 4 0.4   respectively. The simulation results of the PL algorithm and UPPL algorithm are shown in Fig. 21(a) and 21(c) after 200 iterations, and their spectrum are shown in Fig. 21(b) and 21(d), respectively. It is clear that the results obtained by the UPPL algorithm are better than the results of the PL algorithm. The super-resolution capability (spectrum extrapolation) of the UPPL algorithm is more obvious. In the second experiment, the visual image is shown in Fig. 23(a). The blurred image is shown in Fig. 23(b) acquired by a 91.5 GHz mechanically scanned radiometer. We use the same experimental parameter as the first experiment. The restored images are shown in Fig.  23(c) and 23(d) respectively. It is clear that a significant resolution improvement is achieved by the UPPL algorithm.
A reasonable frequency decomposition scheme, the UPT is also presented. Experiment results demonstrate that the result of UPPL is better than that of the PL algorithm. The UPPL algorithm has lower MSE and produces sharper images. The effectiveness of the UPPL algorithm for practical PMMW images is also validated. Undecimated steerable pyramid decomposition of the image generates a number of different sub-band frequency images. As these details and the differences between the useful signals and noise in the approximate sub-band images, we can use different relaxation parameters and iterations in order to suppress the noise. The UPPL algorithm combines the advantages of the PL algorithms and multi-level pyramid recovery methods. The acquired image are first decomposed by non-sampling pyramid transform into some sub-images. Then we operate the super-resolution process to each grade low-pass image using the PL algorithm. Since each sub-image contains different frequency content, we select different relaxation parameters and iterations in the super-resolution processing. The algorithm improves the ability of spectral extrapolation, and makes the restored image more sharpen.

Two-step projection iteration thresholding supper-resolution using compressed sensing architecture
For passive millimeter wave image, it has certain structure, which can be sparse decomposed in a particular base. Therefore, we can use the sparse prior information of PMMW images for the image reconstruction process. Because the noise is not sparse, the supper-resolution based on sparse prior information can suppress noise. The low-pass effect of the PMMW imaging system make the space resolution very poor. The tw-PIT algorithm uses the sparse prior information and the non-negative finite value information of PMMW images, which can separate the noise from the signals in the iteration process. The algorithm is very effective when the image noise is existed which can also has good super-resolution performance.

Projection Iteration Threshold (PIT) supper-resolution
In the compressed sensing(CS) theory, we use sparsity to describe a signal's feature or structure. Passive millimeter wave(PMMW) images have a certain structure so that they can have sparse representation on some basis. Therefore, during the processing of superresolution, we can use the prior that the PMMW images can be sparse represented to reconstruct them. And in image restoration, the introduction of l1-norm optimization can proform a more effective way to recover a original image.
In the CS theory, we take advantage of the signal's sparsity by bringing p-norm restrict condition in solving the objective function. So far, using iterative way to solve nonlinear reconstruct problem has achieved remarkable results in the field of CS. So in the field of image restoration, we can also bring in p-norm restriction to obtain images' sparse prior information. The process of image degradation as shown in section 2, for the model of PMMW imaging system which is supposed space-invariable, the operator K is simplified a convolution process. And the solution of this problem requires minimizing the difference between the optimal solution and the real one: Unlike classical regularization way, we add regularization to the sparse prior information of image, where the constraint is not quadric, but the lp-norm( 12 p   ) of signal f. It is here that the introduction of p-norm make the solution of objective function have sparsity.
The objective function to be optimized is: Where   is the orthogonal basis which f could be sparse represented with, and ( ) ww    is the weight of the coefficients in transform-domain. When K is a identitymatrix, and   is a wavelet basii function, the objective function mentioned above turns to denoising via wavelet transform based on Besov priori information .
Then, the corresponding variational equation is: The nonlinear equation shown above which involved symbolic function is a tricky one in practice. Define constant C satisfies: In equation (43), when  w0 , it means the objective function doesn't include any sparse priori information. Therefore, the algorithm is degraded to a process of landweber iteration algorithm.
When p=1, the variational equation of the objective function is expressed as: For the process of symbolic function, when f  >0, the solution is: otherwise, when the equation do not satisfy the two conditions mentioned above, which is In conclusion, the iterating thresholding algorithm based on sparse prior is described as follow Where ,1 w S  is defined as follow: www.intechopen.com The sparse prior iterating thresholding algorithm will achieve a better effect in image restoration, which is introducted l 1 -norm restricted condition and the sparse priori information. And in the PMMW imaging system, the non-negative limitation of the image can be used as a prior in image super-resolution process. So we add that in our algorithm, and we have the sparse prior super-resolution algorithm restricted by l 1 -norm -projected iterating thresholding(PIT) algorithm. The iterative formula is shown as follows: Where the definition of operator u S is soft threshold operation, and  is the coefficient of the image's sparse representation on a certain orthogonal basis function.
Passive millimeter wave images can be sparsely represented. Assuming the orthogonal basis is 12 [,, . . . , ] N     , the sparsity of image is expressed as follows: Where, f is the scene signal,  is the sparse coefficient in orthogonal basis.
A flow chart depicting the various steps of Projected Iterating Thresh-holding (PIT) superresolution algorithm is shown as Fig. 24. After the above steps, one iteration is finished. Depending on the specific situation, the algorithm can perform multiple iterations until reaching design requirement. The PIT algorithm taking advantage of the sparse prior of the image, can have a better recovery performance. Because of the isotropy, the noise has no sparsity in transform-domain. Since we keep projecting to the l1-ball in every iteration, the algorithm can effectively eliminate the impact of the noise's amplification. We will verify it in the later simulation.
In one iteration, there are two convolutions of images, two wavelet transform operations, two plus-minus operation and thresholding operation. The convolutions of images can perform in frequency-domain, which could reduce calculations by using FFT, and wavelet transform could achieve quickly by using Mallat algorithm.
For a brief description, each cycle of this "tw-PIT supper-resolution algorithm" consists of executing the four steps: Step 1. The iterative principal process Where H is the PSF of system, g is the Captured PMMW image,  is the sparse coefficient in orthogonal basis  .
Step 2. Soft threshold procedure:  is the given parameters.
Step 3. Updating algorithm two-step projection iteration The designation "two-step" stems from the fact that depends on both k f and 1 k f  , rather than only on k f .
Step 4 Tw-PIT supper-resolution algorithm can finish one iteration though the above four steps. The algorithm converges very fast. The convergence is much faster with 1  is much smaller.
But it does not guarantee optimal performance of the reconstruction algorithm.

Simulation and experiments
In order to verify the super-resolution capability of thetw-PITalgorithm mentioned above, two images have been adopted in simulation and experiments, one is a lena image, and the other is an 210×96 University of ElectronicScience and Technology of China (UESTC) library's image captured by 94 GHz mechanically scanned mono channel radiometer.
In the first experiment, the optical lena image and the blurred image are shown in Fig.25a and Fig.25b  It is clear from Fig.25 that the results obtained by the twPIT algorithm are better than the results obtained by the PL algorithm.
When the ideal image is available for comparison, various distance metrics can be readily postulated to compare the images and their spectra. Straightforward measure is the Mean Square Error (MSE), which given by Where f denotes the ideal image, k f  represents the th k restoration result, 2  is 2 L norms. Fig. 26 plots the MSE of the three algorithms.
It can be observed that the twPIT algorithm has a faster decrease of the MSE than the PL algorithm and PIT algorithm from Fig.27. Also, the MSE of the twPIT algorithm is lower than that of the other two algorithms. Clearly, the twPIT algorithm has the best convergence performance and lowest MSE. This is because the twPIT algorithm suppress the noise in an iterative process. The experiments show that the twPIT algorithm convergence faster, which is suitable for the required real-time applications, such as airport security and Concealed Weapons Detection.
In the second experiment, the visual image is shown in Fig.28a. The PMMW image is shown in Fig.28b, which is captured by a 3mm Passive millimeter wave focal plane linear array scanning imaging system. The restored images are shown in Fig.28c-28e, respectively. Fig.27 www.intechopen.com

290
is the PSF which is estimated by variational bayesian. In captured image, the image are assembled together by three times mechanically scanned mono-channel radiometer, as a single mechanical scan pitch direction field limit. The MSE is a global evaluation criterion for super-resolution algorithm. We can use the local image variance method to estimate the image noise level. This approach is based on the assumption, which the image is a large number of small pieces of uniform composition, and image noise to additive noise based. The captured passive millimeter wave images meet this requirement, which there are many flat background and noise is mainly additive noise, basically meet the above assumptions. Local variance calculated as follows [43]: 1. The image is divided into many small block, such as 4×4, 5×5. Then, we calculated the local mean and local variance of each block. The local variance is defined as: Where, (, ) ij indicates the location of the image. , PQ is the local variance calculation window size.
2. The equally spaced intervals are created between the maximum local variance and minimum local variance. The each local variance is enclosed into the appropriate interval. The local variance average of which contains the largest number of extents is used as the image noise variance. For the simulation image and the actual millimeter wave image, we set the block number is 100. In our experiments, the three kinds of super-resolution algorithm including PL, PIT and twPIT are executed for low resolution PMMW images.The symmetric extension technology is adopted in the image boundaries convolution for eliminating shock ringing. The restored result of PL, PIT and twPIT are shown in Fig.28c-28e after 20 iterations, respectively.
The calculated local noise variance is 55.3552, 24.2451 and 34.6997 by PL, PIT and twPIT algorithm, while the 2 PQ   and the window size of the local variance is 5×5.
From above the results, we can see that PL algorithm can effectively perform super-resolution processing, but the noise is magnified in the recovery process. PIT and twPIT algorithm can effectively suppress high frequency noise and make the flat region of original image maintain its flatness in the iterative process and maintain high-resolution of the images.
In the second experiment, the scene is relatively simple millimeter-wave imaging image processing. The Point Spread Function of the imaging system is estimated by the variational Bayesian method. For the PIT algorithm, the parameter  is chosen as 1. For tw-PIT algorithm, the parameter  is settled to 2, the parameter p is set to 0.5. All the algorithms are iterative 100 times, the experimental results as shown below: the covered the car visual image is shown in Fig. 29(a),the unobstructed visual image is shown in Fig. 29(b). The obtainedPMMW image is shown in Fig. 29(c), which is acquired by a W band mechanically scanned radiometer. The spectrum of the PMMW image is shown in Fig. 29(d). The result and its spectrum after100 iterations PL algorithm are shown in Fig. 29(e) and Fig. 29(f). The result and its spectrum after100 iterations PIT are shown in Fig. 29(g) and Fig.29 h). The result and its spectrum after100 iterations tw-PIT are shown in Fig. 29(i) and Fig. 29  In the third experiment, we process the obtained PMMW image is the UESTC Qingshuihe campus's Classroom Building, as shown in Fig.30a  The twPIT algorithm is introduced for the PMMW supper resolution process. The update equation depends on the two previous estimates (thus, the term two-step), rather than only on the previous one. This class contains and extends the iterative thresholding methods recently introduced. This algorithm uses the sparse prior information and the non-negative finite value information of PMMW images, which can separate the noise from the signals in the iteration process. Our algorithm is very effective when the image noise is existed, which can also has good super-resolution performance. The experimental results have shown that PL, PIT and twPIT algorithm can effectively be super-resolution processing, but the noise is magnified in the PL algorithm recovery process. The PIT and twPIT algorithm can effectively suppress high frequency noise smooth the region and maintain its flatness in an iterative process, while maintaining a high-resolution images.
And the twPIT algorithm can in fact be tuned to converge much faster than PL and PIT algorithm, specially in severely ill-conditioned problems. The twPIT algorithm has minimum MSE, which restore original signal more accurate. The reason is that noise is restrained in the iteration process. Sparse prior super-resolution algorithm has good MSE performance is mainly due to which can separate the noise in the iterative process and minimize the impact of noise for image restoration. Thus, the twPIT algorithm can be used in the occasion which requires PMMW imaging quality and real-time such as airports and Concealed Weapons Detection (CWD).
In the second experiment, due to the limited field of view, the single-channel scanning process has introduced a certain artifacts. So eliminate and abate the artifacts is the problem that need to be researched and solved in our follow-up work.

Acknowledgment
This work is supported by the National Natural Science Foundation of China.