Image Restoration for Long-Wavelength Imaging Systems

Basically, the quality of an image can be evaluated on its spatial and spatial-frequency resolutions, image interpolation and superresolution are perhaps the way to respectively produce high spatial and spatial-frequency resolutions of images especially for a single downsampled image. For convenience, the term “hyper-resolution” used here represents the approach to enhancing both the spatial and the spatial-frequency resolutions of an image.


Introduction 1.1 Overview
Basically, the quality of an image can be evaluated on its spatial and spatial-frequency resolutions, image interpolation and superresolution are perhaps the way to respectively produce high spatial and spatial-frequency resolutions of images especially for a single downsampled image.For convenience, the term "hyper-resolution" used here represents the approach to enhancing both the spatial and the spatial-frequency resolutions of an image.
As known, the process of decimation or down-sampling is an effective way often used to reduce image sizes, thus, reducing the amount of information transmitted through the communication channels and the local storage requirements, while trying to preserve as much as possible the image quality.Conversely, the reverse procedure of this, referred to as interpolation or up-sampling, is useful in restoring the original high resolution image from its decimated version or for resizing or zooming a digital image.Decimation and interpolation are used for several purposes in many practical applications, such as progressive image transmission systems, multimedia applications, and so forth.A number of conventional interpolation techniques [Hou & Andrews, 1978;Jain, 1989;Keys, 1981] have been proposed to increase the spatial resolution of an image.Undoubtedly, these techniques degrade the quality of the magnified image.
Furthermore, images may be corrupted by degradation such as blurring distortion, noise, and blocking artifacts.These sources of degradation may arise during image capture or processing and have a direct bearing on visual quality.Various methods of restoration have been described in the literature; this diversity reveals the importance of the problem and its great difficulty.The purpose of image deconvolution or restoration is to recover degraded images by removing noise, highlighting image contrast, and preserving edge features of image.
Image superresolution was developed in 1950s to improve image quality and pilot research of this field is derived from the early work (Toraldo di Francia, 1952Francia, , 1955) ) where the term "superresolution" was used in the paper.Following that, clear definition, description and some of the obvious contribution to this field can be found in the work (Gerchberg, 1974;Hunt & Sementilli, 1992) in which their work, superresolution, was meant to seek to recover image processing techniques can improve the image quality of tomographic images between iterations of image reconstruction.
The technique using image restoration gradually becomes popular for an mm-wave or an NIR DOT imaging system; the difference of both imaging systems is that the former is postprocessing and the latter is inter-processing.

Remark
In this section, we have described a number of terms such as spatial resolution, spatialfrequency resolution, interpolation, restoration, superresolution, hyper-resolution, interprocessing, and post-processing.In addition, advantages and drawbacks of long wavelength imaging systems were addressed and general description of restoration algorithms was made.It is worth emphasizing that long wavelength imaging systems have the same problem to be dealt with so image restoration can be used to improve such an imaging system.Following this introduction, this chapter is organized as follows.Section 2 describes mathematical model of image formation; image restoration algorithms and further consideration on image restoration are explained in Sec. 3 and Sec. 4, respectively.Subsequently, Sec. 5 demonstrates related applications of image restoration.Finally, conclusion is drawn in Sec. 6.

Mathematical model of Image formation
In this section, imaging systems, image formation model, and forward problem and inverse problem are described in the following.

Common imaging system
Usually, the imaging process of a common imaging system is formed as follows.Suppose we have a scene of interest that is going to pass through a common imaging system where it has been corrupted by a linear blurring function and some additive noise.The blurring function h accounts for the imperfectness of the imaging system including optical lens or the human factors in shooting the images.Some typical examples are a diffraction-limited or defocused lens and camera motion or shaking during the exposure.The noise arises from the inherent characteristics of the recording media, e.g., electronic noise and quantization noise when the images are digitized (or discretized).
In practice, the available blurred image not only follows exactly the above description but also is constrained with the film size, in most cases the images have to be truncated at the boundaries.Instead, what is available now becomes a windowed blurred image where a rectangular window is usually accounting for the film aperture shape and size.One inherent problem with this is that many ringing artifacts are introduced into the restored image when the linear or nonlinear filter is applied directly to the truncated blurred image.

Medical imaging system
Here, we use NIR DOT imaging system as an example.Basically, an NIR DOT imaging system is composed of a measuring instrument associated with image reconstruction scheme for the purpose of reconstructing the NIR optical-property tomographic images of phantoms/tissue of interest.The reconstructed images reveal the NIR optical properties of tissue computed by using measured radiances emitted from the circumference of the object.A schematic diagram of the NIR DOT measuring system in the frequency domain is shown in Fig. 2.1.Fig. 2.1.Schematic diagram of NIR DOT measuring system in the frequency domain.

Image formation model
The image formation is modelled as where f is the original scene, h is the point-spread function (p.s.f.) of the imaging system,  is the convolution operator, n is the noise, and g is the corrupted image.Subsequently, the corrupted image is windowed due to the film size/support area and sampled for digitization.
Aliasing is arising, which causes different signals to become indistinguishable when sampled.It also refers to the distortion or artifact that results when the signal reconstructed from samples is different from the original continuous signal.

Forward problem & inverse problem
In a common imaging system, the image is formed as the above description in which finding an estimated original signal/image (f) is an inverse problem for a given corrupted signal/image (g) while the reverse process is a forward problem.In tomographic imaging, the reconstruction of optical-property images is done iteratively using a Newton method, requiring inversion of a highly ill-posed and ill-conditioned matrix.The goal of DOT is to estimate the distribution of the optical properties in tissue from non-invasive boundary measurements.For the purpose of determining the optical properties (the absorption coefficient and the diffusion/scattering coefficient) from measurement data, which is an inverse problem in DOT, a forward model is needed to describe the physical relation between the boundary measurements of tissue and the optical properties that characterize the tissue.

Forward problem in DOT
In general, such a forward model of NIR DOT that gives the description of this physical relation is the diffusion equation, where   ,  r is the photon density at position r and  is the light modulation frequency.

 
, S  r is the isotropic source term and c is the speed of light in tissue.a  and  denote the optical absorption and diffusion coefficients, respectively.In addition, the finite element method (FEM) and a Robin (type-III) [Brendel & Nielsen, 2009;Holboke et al., 2000] boundary condition are applied on Eq. (2.2) to solve this forward problem, i.e., calculating the photon density for a given set of optical property within the tissue.

Inverse problem in DOT
Owing to the non-linearity with respect to the optical properties, an analytic solution to the inverse problem in DOT is absent.Instead, the numerical way of obtaining the inverse solution is to iteratively minimize the difference between the measured diffusion photon density data, M Φ , around the tissue and the calculated model data, C Φ , from solving the forward problem with the current estimated optical properties.This data-model misfit difference is typically defined as follows, where M N is the number of measurements.
By means of the first order Taylor series to expand

Remark
In this section, we have explained a common imaging system which includes the operation of convolution, support area, sampling, and noise as well as a medical imaging system of which the optical-property images are formed with the reconstruction algorithm from 1D signals.

Image restoration algorithms
This section will discuss non-iterative, iterative and statistical methods; in addition, regularization is also used frequently in image restoration algorithms.More descriptions are explained in the following.
As known, the image degradation is basically modelled as where f is the original scene, h is the point-spread function (p.s.f.) of the imaging system,  is the convolution operator, n is the noise, and g is the corrupted image.
Generally, the non-linear iterative restoration algorithms (Archer & Titterington, 1995;Hunt, 1994;Meinel, 1986;Singh et al., 1986;Stewart & Durrani, 1986) to enhance image quality by restoring the high frequency spectrum of the corrupted images can be simply modelled as the following form: where the subscript n is the n-th iteration, Eqs. (

Non-iterative methods
Non-iterative restoration algorithms are described in this sub-section such as the inverse and Wiener filters usually recovering the spatial-frequencies below the diffraction limit.Filters in the Fourier domain are respectively given by the following expressions: However, Eq. (3.4) is not able to be directly implemented; usually, one uses a so called pseudo-inverse filter with a small constant as below.
Pseudo-inverse Filter = 1 Wiener filter is described as Eq.(3.6) in the following.
where H is the modulation transfer function (MTF) of p.s.f.; the superscript asterisk (*) denotes the complex conjugate; [ n / f ], the ratio of noise-to-signal. n and  f represent the power spectral densities for noise and the true images, respectively.Apparently, applying the Wiener filter to the restoration problem has to know the power spectral densities for the noise and the original image (or more precisely, their ratio).Unfortunately, this a priori knowledge is not available in most cases.Nevertheless, the noise-to-signal ratio (NSR), [ n / f ], is usually approximated by a small constant .In such a case, the Wiener filter becomes Wiener filtering achieves a compromise between the improvement obtained by boosting the amplitude of spatial-frequency coefficients up to the diffraction limit and the degradation that occurs because of the noise amplification of the inverse filtering.Noise propagation tends to be reduced by the convolution with p.s.f.; this has a smoothing effect in the result.This fact reveals that Wiener filtering is more immune to noise than inverse filtering.

Recursive wiener filter
This technique is briefly described here; further, a more detailed description of the implementation of this algorithm can be found in the literature [Kundur & Hatzinakos, 1998].
Briefly, such a recursive Wiener-like filtering operation in the Fourier domain can be expressed as Eqs.(3.8) and (3.9). .ˆˆ/ The real constant α represents the energy of the additive noise and is determined by prior knowledge of the noise contamination level, if available.The algorithm is run for a specified number of iterations or until the estimates begin to converge.The method is popular for its low computational complexity.The major drawback of the method is its lack of reliability.The uniqueness and convergence properties are, as yet, uncertain.

Lucy-Richardson method
The Richardson-Lucy algorithm, also known as Lucy-Richardson deconvolution, is an iterative procedure for recovering a latent image that has been blurred by a known point spread function.
The Richardson-Lucy (RL) algorithm has been widely used for the data from astronomical imaging.The RL algorithm (Richardson, 1972;Lucy, 1974) generates a restored image through an iterative method, which is derived using a Bayesian statistical approach to guess the original image (f ), to convolute it (f n-1 ) with the p.s.f.(h) and to compare the result with the real image (g).Usually the guessed image for the first iteration is the blurred image.It uses such an iterative approach:

Poisson MAP algorithm
The Poisson MAP superresolution algorithm begins with Bayes' law associated with Poisson models for the statistics of image and object to estimate the object by finding the maximum probability on the object (f) given the image (g).Mathematically, the Poisson MAP (Hunt & Sementilli, 1992) is given by where  represents a convolution; *, a correlation; n f , the restored signal/image; g is the blurred signal/image; h, p.s.f.; 0 f , the initial guess signal/image; subscript n, the iteration number.Here,

www.intechopen.com
Image Restoration for Long-Wavelength Imaging Systems 237 C can be regard as the correction term during the iterative restoration process.In terms of the operation of the Poisson MAP, it is an iterative algorithm where successive estimate of the restored image is obtained by multiplication of the current estimate by a quantity close to one.The quantity close to one is a function of the detected image divided by a convolution of the current estimate with p.s.f.. Indeed, one can replace the exponential in Eq. (3.12) by the first order approximation e x ~ 1+x because of low contrast in a blurred signal/image to achieve Eq. (3.13).
Equation (3.14) shows that the Poisson MAP superresolution is consistent with Eq. (3.2).Experience reveals that when implemented for simple point objects, the Poisson MAP algorithm is able to expand the bandwidth much more than done for more complex objects and the Poisson MAP superresolution algorithm requires hundreds of iterations for a final solution.

Improved P-MAP
Following that, the Poisson MAP can be improved by itself by operating upon the edge map with a re-blurring technique; that is, the g and f n-1 of the Poisson MAP are replaced by the corresponding gradients of the g ⊗ h and f n-1 along with the integrated p.s.f.(h ⊗ h).Mathematically, it is shown that Thus, the final hyper-resolved image f can be obtained by integrating (f n )'.The whole process of this improved Poisson MAP includes re-blurring, differentiation, restoration, integration, and then correction for a DC offset.More details concerning this algorithm can be found in the author's previous work [Pan, 2003].

Regularization
Regularization presents a very general methodology for image restoration.The main technique of a regularization procedure is to transform this ill-posed problem into a wellposed one.Roughly speaking, restoration problem with regularization comes down to the minimization problem [Chen et al., 2000;Landi, 2007].
In our real life, one cannot get the whole blurred and noisy images but only can get part of blurred and noisy images because of the limited support size.According to the part of blurred and noisy image, ones want to reconstruct an approximate true image by deconvolving the part of blurred and noisy image.Thus, noise (n) in general meaning should include both additive noise (n add ) and the effect of the limited support size (n limited ) at least.Normally, Q 1 is usually used with a true h which is, however, not known and optimal, whereas Q 2 is expected to be used with an ĥ , which is supposed to be optimal in practice.
Here, Q 2 is proposed for the purpose of reducing the error energy coming from noise and ringing artifacts while only Q 1 is considered.Thus, a new objective function combines Q 1 with Q 2 , and its regularization term is ; it is approaching to null when iteration is increasing.Finally, we define an objective function as Eq.(3.17) where is the regularization parameter and then minimize Eq. (3.17) with respect to f n-1 ; i. e. , where Note that h  in Eq. (3.25), normally, is equal to h but it is chosen as a user-guess p.s.f.when h is unknown.Here, h hp is chosen ash  , where a delta function and a Gaussian function adopted for h lp1 and h lp2 in numerical simulation, respectively.Equations (3.23)- (3.25) show that the restored signal/image can be obtained from the increment iteratively updated using the projection of the high frequency spectra of the increment.As discussed, h hp is defined as the difference of a delta function and a Gaussian function; in addition, an edge operator like a Laplacian operator defined as Eq.(3.26) is adopted for h hp in the following experimental verification.

Remark
In this section, we have established a framework of image restoration/superresolution including (pseudo) inverse filter, Wiener filter, recursive Wiener filter, Lucy-Richardson method, Poisson MAP algorithm, and improved P-MAP algorithm.Of restoring image quality and reducing ringing artifacts, the error-energy-reduction-based regularization algorithm has been proposed here for long-wavelength imaging systems as well.

Further consideration on image restoration
In this section, the topics of improvement of spatial resolution, rapid convergence, and inverse pitfall for image restoration are described.

Improvement of spatial resolution
Usually, hyper-resolution of a noisy image is considered as an interpolation followed with restoration/superresolution; generally, the procedure for processing noisy images is shown in Fig. 4.1(a), that is, noise removal, interpolation, and then superresolution, whereas the proposed scheme is dealing with interpolation and noise removal simultaneously, as shown in Fig. 4.1(b).

Rapid convergence
As known, restoration/superresolution or the reconstruction of optical-property images with an iteration procedure is usually computed off-line and computationally expensive.Most of studies, however, focused mainly on improving the spatial and spatial-frequency resolutions.If a real-time resolution processing is required, dedicated reconstruction hardwares or specialized computers are mandatory.Moreover, fast reconstruction algorithms should also be considered to reduce the computation load.It is worth emphasizing that our proposed method can reduce computation time with the regularization term which is designed on the viewpoint of the update characteristics in the iteration procedure but not utilizing any spatial/spectral a priori knowledge or constraints; some results can be found in the author's work show how to speed up the computation to find an inverse solution for reconstructing optical-property images by using regularization with an iteration domain technique; similarly, this proposed method is capable of being applied to image restoration/superresolution for other imaging systems.

Algorithm of rapid convergence
Image reconstruction tasks contain forward modeling and inverse problem.The forward computation consists in obtaining the intensity out of a subject under investigation for a given source, and the initial-guess (or iterated result) on scattering and absorption coefficients.The inverse computation is to compute the scattering and absorption coefficients for a known light source and measured intensities in an iterative manner.
Since we utilize cw light illumination or DC data, the physical process of NIR light illuminating through a highly-scattering medium can be approximated by the steady-state diffusion equation ,  is the vector composed of D k and  l , and  is the vector with differences between calculated intensities (Φ cal. ) and measured intensities (Φ meas.).Also, D k for k = 1, 2, …, K and l for l = 1, 2, …, L are the reconstruction parameters for the optical-property profile.The opticalproperty image reconstruction is actually a process of successively updating the distribution of optical coefficients so as to minimize the difference between measured intensities and computed ones from the forward process.More details can be found in [Paulsen and Jiang, 1995] where the Levenberg-Marquardt procedure was adopted to update the diffusion and absorption coefficients iteratively.
It is known that to solve Eq. (4.4) is an ill-posed problem.Tikhonov regularization is a method stabilizing the inverse problem through incorporating a priori assumptions to constraint the desired solution.It is able to convert an ill-posed problem into a well-posed one, and further to improve an ill-conditioned problem.The regularization term (penalty term) introduced in the process regularizes the problem and makes the update stable.It also strengthens the robustness of algorithm to noisy data with the adequate design of the regularization term.Generally, Tikhonov regularization is to optimize this ill-conditioned problem as where () is a constraint on the estimate , and E is a quantity confining the constraint to be an energy bound.Applying Lagrange optimization technique, we seek a solution to the constrained objective function where λ is referred to as the regularization parameter.A solution to Eq. (4.7) is given by 2( ) 0 and equivalently where Eq. (4.9) is a constrained estimate of , but becomes an unconstrained one when λ equals to zero.It is noted that the minus sign in Eq. (4.6), the objective function, corresponds to the regularization term proposed here as the term is constrained to an energy bound.

Constraints on the spatial domain
A constraint on the spatial domain can generally be expressed as If L is the identity matrix (I), a solution to Eq. (4.9) is given by On the other hand, if L is the discrete Laplacian matrix, substituting Eq. (4.10) into Eq.(4.9), the corresponding solution is Equation (4.11) is usually a primary inverse solution to optical-property image reconstruction, which is also Levenberg's contribution to the inverse problem; and Eq.(4.12) is a constrained inverse solution implemented to improve the quality of the reconstructed NIR DOT images, which is identical to Marquardt's work.

Constraints on the iteration domain
In NIR DOT, it is also crucial to accelerate the computation.But, up to now, speeding up the computation in the iteration domain has not been explored yet.Here we consider this issue through the use of a Lorentzian distributed function taking a natural logarithm computation as a constraint, i.e.where p is the calculated nodes in the subject under investigation and is a user defined positive parameter.As can be seen, The Lorentzian distributed function, as depicted in Fig. 4.4, is employed here owing to its following two characteristics: a. Lorentzian distributed function has a sharp peak with a long tail, describing the histogram distribution of Δχ, many of Δχ (~0) at its peak and a small rest of Δχ distributing along its long tail, and b. its histogram distribution can be further tuned with the parameter ( ) as iteration increasing.Related to the consideration in convergence, the updated quantity, Δχ, decreases, ranging from the peak to the tail, as the iteration increases whereas it has a smooth distribution in the beginning stage of iteration.
In addition, as the shape of the histogram would be affected, it is smooth with a big value of and sharp with a small value of .Thus, Lorentzian distributed function can characterize the nature of Δχ in the iterative process as the distribution from a smooth to a sharp distribution to be used as a constraint for the purpose of speeding up computation.) at various .As can be seen, it has a smooth distribution for a big and a sharp distribution as is small.

Inverse pitfall
The ill-posed nature of inverse problems means that any restoration or reconstruction algorithm will have limitations on what images it can accurately reconstruct and that the images degrade with noise in the data.When developing a restoration or reconstruction algorithm it is usual to test it initially on simulated data.Moreover, the restoration or reconstruction algorithm typically incorporates a forward solver.A natural first test is to use the same forward model to generate simulated data with no simulated noise and to then find that the simulated data can be recovered fairly well.If one is fortunate enough to have a good data collection system and phantom, and someone skilled enough to make some accurate measurements with the system, one could then progress to attempting to reconstruct images from experimental data.However, more often the next stage is to test further with simulated data and it at this stage that one must take care not to cheat and commit a so-called inverse pitfall or inverse crime.Simply to say, inverse pitfall or inverse crime arises from the reason of 'limited for infinite', e.g., limited support area for infinite scenery, finite elements for continuous zone, or given noise for unknown noise.The best practice is to use a forward model independent of an inverse model.For example, in the case of a finite element forward model one would use a much finer mesh while a coarse mesh is used in the inverse model.

Remark
In this section, we have proposed some extra points about image restoration.Interpolator with noise removal, design of regularization term for reducing computational burden, and inverse pitfall/crime have been illustrated and discussed.

Related application
In this section, application to a mm-wave imaging system or near infrared diffuse optical tomography using image restoration is demonstrated for post-processing or interprocessing.To verify the proposed method in the previous section (Sec.3.4), a computergenerated signal/image and an image of real scene were tested.

Post-processing: Application to a millimeter-wave imaging system [Pan, 2010]
A 1-D noiseless signal and a 2-D noisy image were used, originally blurred with a p.s.f. of Gaussian function plus additive white Gaussian noise.White Gaussian noise is defined with a zero mean and variance, σ 2 , specified by a blurred signal-to-noise ratio (BSNR).Recall that where M, N are the dimension of the processed image and i, j are the indexes of a pixel and X means the average value of X.In many practical situations, the blur is often unknown and little information is available about the true image; therefore, several h   of the Gaussian blur around the true σ h were tested in the following examples; f o and α are chosen to g and g , respectively.In this work, the stopping criterion is or 1% (for 2-D image).The mean square error (MSE) of the restored signal/image relative to the original signal/image is provided here for the evaluation of image quality, thus supporting the visual assessment.
The proposed algorithm was applied to a 1-D signal as well as both simulated and real atmospherically degraded images, one of a simulated blur and one of a real blur.The purpose of the simulation was to enable a comparative evaluation of the results given the original signal/image and to explain the algorithm characteristics.In the real-blur example shown here, a 256 × 400 pixel millimeter-wave image was tested and the image was captured at 94 GHz by the Defence Evaluation and Research Agency, Malvern, UK.
For a comparison purpose, non-iterative Gaussian filtering was used in the case of 1-D signal and the common Richardson-Lucy (RL) deconvolution method was implemented using a built-in MATLAB function deconvlucy in the cases of both 1-D signal and 2-D images.This RL method employs an iterative procedure to estimates the original signal/image, and therefore requires an initial guess of it as well.It is clear that this restored signal is considerably better than the blurred signal shown in Fig. 5.1(c) whereas the restored signal using the RL method reveals lots of ringing artifacts.Figure 5.1(f) shows that the result using the proposed algorithm with h hp equal to -h  ( h   =1.5) presents higher contrast and less ringing artifact than other two methods.Following the above discussion, Fig. 5.2 shows the iterations used by the RL method and the proposed algorithm satisfying with the stopping criterion.In the case of 1-D signal, our algorithm usually converges within fewer iterations than the RL method, the former using 34 iterations and the latter using 187 iterations.3 shows that the nature of our proposed method possesses the ability to reconstructing frequency spectrum beyond the diffraction limit, where a 1-D noiseless signal was used.  equal to 1.2, 1.5, and 1.8, respectively; and the MTFs of the original and the restored signals are depicted in Fig. 5.3(g)-(i).The restored signals in Fig. 5.3(e) and (f) display the performance of high resolution and the two peaks are separated in Fig. 5.3(d) even with a small h   .Compared with that of the original signal, high-frequency information of the restored signals was definitely generated beyond the diffraction limit as shown between the two dashed lines in Fig. 5.3(g)-(i), explaining that the proposed method possesses the highresolution ability.  =2.5 was used; the MSEs of these three results are 181.17,49.45, and 52.61, respectively.

Results for synthetically blurred signal and image
These three restored images demonstrate high quality but Fig. 5.4(c) still shows ringing artifact especially in the boundary of the image.In Fig. 5.4(d), simultaneously, the image quality can also be improved by reducing most of the ringing artifact and preserving more edge information.Also, it can be seen that our method with a Laplacian filter still works well, shown in Fig. 5.4(e).Corresponding to Fig. 5.4(c)-(e), Fig. 5.5 shows the iterations used by the RL method and the proposed algorithm where fewer iterations was used in the RL method than our algorithm, the former using 46 iterations and the latter two using about 200 iterations.It should be noted that the proposed algorithm is considerably more computationally expensive than the RL method.However, in our experiments we did not find any significant improvement but even more ringing artifacts when the RL method was employed for a further iteration number.For further inspection into our proposed algorithm, we investigated the effect of this algorithm using the high pass filter, -h  , with varied h   .

Results for a real degraded image
It is always expected that a novel algorithm can be implemented on a real image; Fig. 5.7(a) presents a real degraded image captured by an mm-wave imaging system.where the subscript n is the n-th iteration, "max" means the maximum value, and the superscript T denotes a transposition operation.One way to improve the convergence rate is using n    as the Type-1 soft prior and using    5.9 (a-f) shows the 2D reconstructions of phantoms with two and three inclusions, where slight discrepancy can be observed.Figure 5.9 (g-l) depicts their corresponding 1D circular transection profiles to reveal noticeable differences.Basically, there is a better separation resolution but a lower intensity owing to a highly suppressed signal by a hard prior rather than a soft prior.Additionally, Fig. 5.9 (m-o) exhibits good convergences obtained by using both soft and hard priors.

Image restoration applied to NIR DOT
The phantoms employed for justifying our proposed technique (Sec.3.4) incorporate two or three inclusions with various sizes, locations and separations, illustrated in Fig. 5.10, where R denotes radius in the unit of mm.Of the phantom, the background absorption ( a ) and reduced scattering ( ' s ) values are about 0.0025 mm -1 and 0.25 mm -1 , respectively, while the maximum absorption and reduced scattering for the inclusion are 0.025 mm -1 and 2.5 mm -1 , thereby assuming the contrast ratio of the inclusion to background 10:1, because high contrast results in much more overlapping effects than low contrast although a contrast of 2~10 were used throughout other published works.As depicted in Fig. 5.10, Case 1, 2, respectively, have two inclusions separated with a similar distance but different sizes.As the separation resolution of inclusions is examined, several (two or three) embedded inclusions are necessary, and different inclusion sizes are considered as well.For the convenience in discussion latter, we denote M0-4 as the reconstructions with the schemes using non-filtering, -g2 (σ 2 =1.5), g1-g2(σ 1 =0.75, σ 2 =1.5), wavelet (a dilated factor a=0.5), and Laplacian high-pass filter (HPF) in their 2D form, respectively.Currently, absorption-coefficient images are presented for our continuous wave image reconstruction algorithm.
In FEM-based image reconstruction, the homogeneous background ( a = 0.0025mm -1 , ' s = 0.25mm -1 ) was adopted as an initial guess.Thirty-iteration assignment was employed for each case as the normalized increasing rate, i.e. mean value of

Case 1
This case was designed as a phantom with three smaller inclusions.Several improved images were obtained by using appropriate filtering, as shown in Fig. 5.11(b-e) of 1D circular profiles passing through the centers of inclusions.Likewise, M2 resulted in worse resolved image than others with HP filtering.Negative artifacts occurred in each reconstructed image, as depicted in Fig. 5.11(g-j).It is well noted that M4 overestimated the inclusion amplitudes, which yields a higher inclusion-to-background contrast.

Case 2
In this highly challenging case, a phantom with two closest-separation inclusions was designed.As shown in Fig. 5.12(a-e), all reconstructed images underestimated inclusions, and offered relatively bad resolution for two separate inclusions.It is rather competitive for these employed filters.Based upon a quantitative comparison, as depicted in Fig. 5.12(i) and (j), M3 and M4 schemes demonstrate better resolution discrimination to separate bigger and closer inclusions in comparison of Case 1.
From the results of Case 1 and 2 for a phantom with inclusions of both small size and close separation, it can be concluded that the wavelet-like HP filtering (M3) demonstrates the best spatial-frequency resolution capability to the inclusions.
It evidently shows that the enhancement of reconstruction through the incorporation of our proposed HPF approach can effectively improve computed images.As illustrated above, the wavelet-like HP filtering schemes (M3, M4) further yields better results than the LPFcombined HP filtering schemes (M1, M2).In the aspects of sensitivity and stability of evaluation, M3 yielded results closest to the true absorption property than other schemes.However, M4 visually characterizes the inclusion-to-background contrast best.

Performance investigation
In terms of the optical properties within the inclusion and background, it is worth noted that the image reconstruction is not only pursuing qualitative correctness but also obtaining favorably quantitative information about the optical properties of either the inclusions or background.Parameters of interest such as size, contrast and location variations associated image quantification measures are most frequently investigated and discussed.Readers can refer to the research work [Pan et al., 2008].

Remark
In this section, we have demonstrated the performance of our proposed image restoration algorithms exactly applied in the imaging process for 'inter-processing' and to corrupted images for 'post-processing.'

Conclusions 6.1 Concluding remark
In this chapter, we have explained the background and the mathematical model of image formation and image restoration for long-wavelength imaging systems; as well, image restoration algorithms, further consideration on image restoration, and their related application have been described and demonstrated.In the meanwhile, a promising method to restore images has been proposed.As discussed in this chapter, the proposed algorithm was applied to both simulated and real atmospherically degraded images.Restoration results show significantly improved images.Especially, the restored millimeter-wave image highlights the superior performance of the proposed method in reality.The main novelty here is that error energy resulting from noise and ringing artifact is highly suppressed with the algorithm proposed in this chapter.Also, we have used such a resolution-enhancing technique with HP filtering incorporated with the FEM-based inverse computation to obtain highly resolved tomographic images of optical-property.
In addition, we have developed and realized the schemes for expediting NIR DOT image reconstruction through the inverse solution regularized with the constraint of a Lorentzian distributed function.Substantial improvements in reconstruction have been achieved without incurring additional hardware cost.With the introduction of constraints having a form of the Lorentzian distributed function, rapid convergence can be achieved owing to the fact that decreasing Δχ results in the increase of as the iteration process proceeds, and vice versa.It behaves like a criterion in the sense of a rapid convergence that the optimal iteration number is founded as seeking an inverse solution regularized with the Lorentzian distributed function.

Future work
It is anticipated that of regularizing mean square error (residual term) with error energy reduction and rapid convergence (a priori terms) an algorithm is explored to restore images effectively and efficiently.In addition, it is no doubt that image restoration for interdiscipline application is the focus in the future research.

Fig. 3 .
Fig. 3.1.A schematic diagram of forming a real image and proposing our algorithm.To develop the novel algorithm with regularization, we plot a schematic diagram, Fig.3.1, to show the mechanism of the concept proposed here and thus define the following functions, Equation(3.16).

. 1 ) 3 )
where()  S r and ()  r denote the source and the intensity, respectively, as well as ()a  r , cand() r D are the absorption coefficient and the diffusion coefficient, respectively.For solving Eq. (4.1), the boundary condition,       D n Flux  , and finite element method are employed.Thus, the following discrete equations can be obtained [Paulsen and Jiang, 1995] AC   , (4.2) where A and C are matrices dependent on the optical properties and the source-detection locations, respectively.The forward solution,  , can be explicitly evaluated by Eq. (4.2).Partially differentiating Eq. (4.2) with DWith an approximation to applying the Newton-Raphson method and ignoring higher order terms, we obtain

2 (
L can be the identity matrix (I) or the discrete Laplacian matrix[Pogue et al., 1999;Davis et al., 2007].

14 )
the requirement of Eq.(4.5).Performing the differentiation indicated in Eq. (4.9), we can obtain the solution in an iterative For further inspection in Eqs.(4.13) and (4.14), as known, a and D are generally searched in a range of [10 -3 :10 -1 ] mm -1 and mm, respectively; and thus Δχ is much smaller than a unit.It can be proven that even the use of the natural logarithm in the constraint Ψ(Δχ) still makes it a positive and finite value.The other reason to use ln is because the regularization term in Eq. (4.14) still remains in a form of the Lorentzian distributed function derived from the constraint associated with the Lorentzian distributed function in Eq. (4.13).

Figure 5 .
Figure 5.1(a) and (c) present an original signal containing 256 pixels and a blurred version of this signal, obtained by convolving it with a Gaussian function with h  equal to 1.5, Fig. 5.1(b), which approximates an atmospheric blur.Figure 5.1(d)-(f) show a comparison between the results obtained from the implementation of Gaussian filtering, the RL deconvolution method and our proposed algorithm, the MSEs of which are 188.29,210.23, and 184.50, respectively.The resulting Wiener-filtered restored signal (with = 0.001) is shown in Fig. 5.1(d).It is clear that this restored signal is considerably better than the blurred signal shown in Fig.5.1(c) whereas the restored signal using the RL method reveals lots of ringing artifacts.Figure5.1(f)shows that the result using the proposed algorithm with h hp equal to -h  ( h

Figure 5 .
1(d)-(f) show a comparison between the results obtained from the implementation of Gaussian filtering, the RL deconvolution method and our proposed algorithm, the MSEs of which are 188.29,210.23, and 184.50, respectively.The resulting Wiener-filtered restored signal (with = 0.001) is shown in Fig. 5.1(d).

Figure 5 .
Figure 5.3  shows that the nature of our proposed method possesses the ability to reconstructing frequency spectrum beyond the diffraction limit, where a 1-D noiseless signal was used.Figure 5.3(a)-(c) shows the original signal, p.s.f. and its modulation transfer function (MTF); the degraded (σ h = 1.5) and the restored signals are shown in Fig. 5.3(d)-(f) with h

Figure 5 .
Figure 5.3  shows that the nature of our proposed method possesses the ability to reconstructing frequency spectrum beyond the diffraction limit, where a 1-D noiseless signal was used.Figure 5.3(a)-(c) shows the original signal, p.s.f. and its modulation transfer function (MTF); the degraded (σ h = 1.5) and the restored signals are shown in Fig. 5.3(d)-(f) with h

Fig. 5 . 3 .
Fig. 5.3.Demonstration of the high resolution of the proposed algorithm.(a) Original signal, (b) Gaussian form (solid line) and MTF (dotted line) of the blurring function (σ h = 1.5),(c) blurred signal, (d)-(f) restored signals withh  and h   =1.2, 1.5, and 1.8, respectively, and (g)-(i) MTFs of the blurred (solid line) and the restored (dotted line) signals.Note that the region between two dashed lines is the high frequency beyond the diffraction limit.

Figure 5 .
Figure 5.4 represents an image (256 × 256) of clown which is a built-in image in MatLab.Figure 5.4 displays a comparison between the results obtained from the implementation of the RL deconvolution method and our proposed algorithm.Figure 5.4(a) shows the original image, convolving it with a 2-D Gaussian function with h

Figure 5 .
Figure 5.4 represents an image (256 × 256) of clown which is a built-in image in MatLab.Figure 5.4 displays a comparison between the results obtained from the implementation of the RL deconvolution method and our proposed algorithm.Figure 5.4(a) shows the original image, convolving it with a 2-D Gaussian function with h

Figure 5 .
Figure 5.4 represents an image (256 × 256) of clown which is a built-in image in MatLab.Figure 5.4 displays a comparison between the results obtained from the implementation of the RL deconvolution method and our proposed algorithm.Figure 5.4(a) shows the original image, convolving it with a 2-D Gaussian function with h  equal to 2.5 to obtain a blurred image shown in Fig. 5.4(b).Figure 5.4(c)-(e) show the images restored with the RL deconvolution method and our proposed algorithm withh  and the Laplacian filter where

Figure 5 .
4(c)-(e) show the images restored with the RL deconvolution method and our proposed algorithm withh  and the Laplacian filter where , hh

Figure 5 .
6 demonstrates this case where the original and the noisy (σ h = 2.5 and BSNR=30 dB) images are displayed in Fig. 5.6(a), (b), and the restored images are shown in Fig. 5.6(c)-(e) obtained with the use of h   equal to 2, 2.5, and 3, respectively.The MSEs of these results are 163.77,76.82, and 97.97, respectively.Of all the restored images, Fig. 5.6(c) shows a worse image quality than the others, in which noise was intensively produced and hard to be removed although the contrast of the restored image was enhanced.

Fig. 5 . 9 .
Fig.5.9.Reconstruction data through various priors with intensity signals corrupted by Gaussian white noise (SNR=20 dB).Left column: constrained inverse solution with soft prior 1; middle column: constrained inverse solution with soft prior 2; right column: constrained inverse solution with hard prior.

Figure 5 .
Figure 5.9 illustrates the comparisons between constrained solutions using soft priors (Type 1 and 2) and a hard prior, where the left, middle and right columns are the constrained inverse solutions with soft prior 1, soft prior 2, and hard prior [M.-Cheng & M.-Chun Pan, 2010], respectively.Figure5.9(a-f) shows the 2D reconstructions of phantoms with two and three inclusions, where slight discrepancy can be observed.Figure5.9 (g-l) depicts their

Figure
Figure 5.9 illustrates the comparisons between constrained solutions using soft priors (Type 1 and 2) and a hard prior, where the left, middle and right columns are the constrained inverse solutions with soft prior 1, soft prior 2, and hard prior [M.-Cheng & M.-Chun Pan, 2010], respectively.Figure5.9(a-f) shows the 2D reconstructions of phantoms with two and three inclusions, where slight discrepancy can be observed.Figure5.9 (g-l) depicts their

Fig. 5 .
Fig. 5.10.Schematic diagram for the dimensions of two different test cases in simulation.(a) and (b) are Case 1, 2, respectively, where R is radius in the unit of mm.
of p.s.f.) have been introduced.Furthermore, h hp can be designed as a highpass filter such as h lp1h lp2 in general orh lp in the extreme case where h lp1,2 are low-pass filters.part of Eq. (3.21) and use the projection of the right pat in Eq. (3.21) on Δf n for the purpose of true value invariance.Consequently, the new relation function, Eq. (3.23), can be achieved for our novel method and expressed as