## Abstract

The problem of image restoration of space variant blur is common and important. One of the most useful descriptions of this problem is in its algebraic form I=H∗O, where O is the object represented as a column vector, I is the blur image represented as a column vector and H is the PSF matrix that represents the optical system. When inverting the problem to restore the geometric object from the blurred image and the known system matrix, restoration is limited in speed and quality by the system condition. Current optical design methods focus on image quality, therefore if additional image processing is needed the matrix condition is taken “as is”. In this paper we would like to present a new optical approach which aims to improve the system condition by proper optical design. In this new method we use Singular Value Decomposition (SVD) to define the weak parts of the matrix condition. We design a second optical system based on those weak SVD parts and then we add the second system parallel to the first one. The original and second systems together work as an improved parallel optics system. Following that, we present a method for designing such a “parallel filter” for systems with a spread SVD pattern. Finally we present a study case in which by using our new method we improve a space variant image system with an initial condition number of 8.76e4, down to a condition number of 2.29e3. We use matrix inversion to simulate image restoration. Results show that the new parallel optics immunity to Additive White Gaussian Noise (AWGN) is much better then that of the original simple lens. Comparing the original and the parallel optics systems, the parallel optics system crosses the MSEIF=0 [db] limit in SNR value which is more than 50db lower then the SNR value in the case of the original simple lens. The new parallel optics system performance is also compared to another method based on the MTF approach.

© 2009 Optical Society of America

## 1. Introduction

The development of digital cameras now allows us to design a hybrid optical and image processing imaging system. With a hybrid system we compensate for relatively poor optical performance with sophisticated image processing. The process is done in serial (Fig. 1).

Because image processing is much more flexible and cheaper than optical design, there is a huge advantage in transferring the responsibility for system performance from the optics to image processing. A challenging hybrid system will be one based on a simple lens with a large Numerical Aperture (NA). Such a system will work in weak luminance but suffer from space variant blur due to large aberrations. In a hybrid system, the blurred image is fixed by image processing. We shall see that the main methods for restoring a space variant blurred image, are based on solving an inversion problem. A main difficulty in using inversion is the matrix condition. In this paper we will show a new method for improving this matrix condition by improving the optical design. In Chapter 2 we present a short survey of the basic methods of space variant image processing. It is shown that most of the methods are based on solving a linear system inversion problem in which the optical system is a matrix. It is also shown that in various methods, solution accuracy is limited by the matrix condition whose figure of merit is the matrix condition number. In Chapter 3 we deal with the question of how to improve the system condition. We present a new approach in which we use the original system Singular Value Decomposition (SVD) to identify the needed missing parts of the system, which, if added, would improve the system condition. Then we show that those missing parts can be interpreted as another matrix, which can be added to the original system matrix. The physical expression of this matrix addition, is adding a second imaging system which works parallel to the original one. We suggest to apply the new second system in a ring around the original lens, which we call “Rim-ring” zone. Then we analyze the diffraction model of the assembled system. Further in Chapter 4 we deal with the basic design of such an optical filter and its design factors. In Chapter 5 we present a study case simulation of designing a parallel optics filter to a simple lens. In Chapter 6 we broaden the discussion and give more simulation results and compare them to another method based on the MTF approach. In Chapter 7 we present effective aperture calculations. Finally, Chapter 8 concludes our research.

## 2. Introduction to space variant image processing

A digital image system can be conveniently represented using matrix representation. In this representation the discrete object and image are input and output vector columns x and b respectively, where the optics itself is represented by matrix **A**, the columns of which are the optics Point Spread Function (PSF) for each field point. We get than:

In case of space variant blur the well established methods of de-convolution cannot be used, and other image processing methods are needed. In the following section we present briefly the basic alternatives for restoring a sharp image from the space variant blurred image based on the system matrix.

#### 2.1. Basic restoration methods for space variant blur images

In some cases optical system aberrations are coma-like. For this particular situation, Robbins [1] suggested an inversion method based on Mellin transform and one dimensional Fourier transform. In the general case of space variant blur, under the limitation of boundary condition continuance (to prevent edge effects), one can divide the image into a mosaic of pseudo-space invariant zones [2,3,4]. If the mosaic method is not suitable, other methods should be considered. Probably the most widely used methods for space variant image restoration are those based on matrix representation.

#### 2.2. Important matrix properties for numerical solution

Knowing the PSF matrix, one will naturally try to restore the object quality by multiplying the image with the inverse matrix. Matrix inversion depends on the matrix rank and matrix condition. The matrix condition and its influence on restoration accuracy depends on the matrix singular values. The figure of merit for the matrix condition is the condition number:

where σ_{1},σ_{n} are the highest and lowest singular values of matrix **A** SVD. A brief explanation for the SVD method is presented below. Even in cases in which **A** is an invertible matrix det(A)≠0, matrix inversion is a sensitive procedure. We define relative accuracy as the following ratio:

‖ ‖_{2} is the vector norm. When applying Gaussian Elimination with pivoting [7] the relative accuracy is proportional to the system conditional number. Wilkinson [6] Shows that when adding an error to the blurred image before restoration the relative accuracy is proportional to the noise to signal ratio multiplied by the condition number.

A wide and extensive field of solutions to linear systems is using least square methods [7, 8]. When solving a linear system one might find that there is no exact solution. The source of this ambiguity might be perturbations in matrix **A**, perturbations in vector **b** or due to a rank deficiency or ill condition problem. In all of those cases we are forced to abandon the exact solution and replace it, with a minimal error in some sense. In the new approach of minimal error we exchange Eq. (1) in the moderate problem

We define the error as:

We define the norm2 over vector r as cost function. To find the least square solution of system Eq. (5) we look for the minimal cost function:

Every matrix has a SVD [7] (we use this decomposition later on in the paper Eq. (24)). This method is very comfortable to use for least square solution. Assuming a mxn matrix, m>n and with rank r=n, the least square solution will be [7]:

The error of this least square solution (LS) will be:

For a full column rank system with some perturbation in matrix value and additive noise in the blur image, the error in the relative accuracy Eq. (3) will rise:

A development of the SVD is the Truncated SVD (TSVD). The TSVD contains a variety of algorithms all based on truncating those part of the SVD which amplify noise and rounding error more than they will add to the accuracy of the solution. For example, if the signal is an eigen vector, only the singular value which associates with this eigen vector is important to the solution, where all the other can be truncated [5]. The linear system in Eq. (1) after truncating is reduced to:

where (k) stands for the order of the eigen component (σ_{k}) involved in the solution. Higher eigen components between k<i<n are truncated. The description of **b** will then be:

Varah [5] analyses a case of perturbation in matrixes **S**,**V**,**D**. He also assumes that most of the information is stored in the low component of **b** decomposition Eq. (14). In this case the relative error between the least square solution to the k order solution is:

η1 is the roundoff error, K_{2} is the arithmetic precision factor. The solution error rises with inverse proportion to the system’s weakest singular values. Since σ_{1}=1 the solution error is roughly proportional to the system condition number Eq. (1). Nevertheless Gorodnitsky and Rao [9] show that when decomposition of **b**
Eq. (14) includes components in high values of i, (β_{i}/σ_{i}) increases with index (i). In that case the error is very much dominated by the image **b** an less by the matrix condition [9]. Another closely related method for solving an inverse problem is regularization. Regularization is based on the Tikhonov method which enlarges the least square cost function in a way which adds a damping component to the TSVD solution. The cost function in this case is [9]:

(b.f) in this case acts as the regularization operator and is represented elsewhere as **C** [10]. For low λ value the solution tends to the least square solution as in Eq. (8) and therefore to accuracy. For high values of λ the dumping component is more significant so the solution tends to smoothness [9], the solution then is:

The regularization error is:

Here again the error is sensitive to the inverse of the system singular values, hence the sensitivity is proportional to a function in the form of the condition number. Iterative methods are in use for solving the inverse problem as well. Following [11] conjugated gradient convergence rate is roughly proportional to:

Apart of the condition number, an additional criterion for solution accuracy for Eq. (5) is the so called Average Relative Error Magnification (AREMF) [12]. Following [13] the AREMF is more tolerant and a better estimation for image restoration than the condition matrix. In the singular value related form:

In this chapter we presented a view in which the major method for solving a space variant imaging system inverse problem is by the mean of its algebraic representation. The quality of this solution depends on the matrix condition and proportional to the matrix condition number or to other close figures.

## 3. Optical design for improvement of the optical system Algebraic representation

In the previous chapter we showed that low condition is associated with low singular value of the matrix SVD. The missing singular values and their eigen-images have optical interpretation. In this chapter we show this interpretation, and that one can improve those singular values by adding another imaging system which works parallel to the original system.

#### 3.1. Matrix representation of a system

Optical space variant imaging systems tend to be ill conditioned [14]. The imaging matrix can be described in matrix formalism as follows:

H is an LxL matrix, it’s columns are the system PSFs for each field point in vector representation Eq. (22). The 2D object and image matrixes transfer into Lx1 column vectors where matrix cell (i,j) is transferred into the column cell as follows:

#### 3.2. Improving the matrix condition using optical filtering

The basic tool in this analysis is the singular value decomposition (SVD). By using SVD one can find the singular values and eigen images of the PSF matrix (H) [15] so that:

Were U and V are found by solving the following eigen problems [15]:

We define a singular value matrix S from matrix Δ, S=Δ^{½}. σ_{i} are positive square roots of Δ. The order of σ_{i} values is high to low [7]:

For each singular value σ_{i}, there exists an eigen matrix M [15]:

If σ_{i}≠0 for all i, one can calculate H^{-1} as follows [7]:

In this representation the inversion of S is solely responsible for the inversion of H. Now, if some of the singular values are equal 0 or numerically approaching zero (assuming the highest value of σ_{max}=1, that might be 10^{-17}), the matrix is un invertible. An intermediate stage might be when σ_{N}>0 but σ_{max}≫σ_{N}>0. We then get a low condition number. To force the inversion or improve the matrix condition one can define a new diagonal matrix S1 which has sufficient values on its diagonal. One can define S1 as follows:

Where ΔS is a diagonal matrix that contains real values at least where they are missing in S. We get a new invertible matrix H1 with a higher condition:

On the right hand of Eq. (32) we have two expressions. The first is equal to H Eq. (23) and the second is an additional new optical system we designate as O:

O is the “optical filter” that works parallel to the original system H. Both optical systems H and O see the same object and bore-sight to the same image plane. A parallel connection of this kind can be done by dividing the lens area into two parts (or more in general). The division can be done simply by dividing an extended lens to a “lens zone” in the center and an “optical filter” in the lens rim zone, named “Rim-ring zone” (Fig. 2). The linear nature of diffraction allows such a separation, which is a division of the integration into two zones. The integration boundaries for the “lens zone” are like those of inner circular pupil function, while the integration boundaries for the “Rim-ring zone” are like those of annular pupil function. In simulation those are continuous pupil functions with two different zones for two different phases. One should note that under this approach the lens area can be divided into as many zones as needed. The zones can be in arbitrary shapes and transmission values.

#### 3.3. Constraints

O design is very flexible, however, as an optical component, matrix O columns are PSFs and therefore must be real.

#### 3.4. Point spread function of parallel optics

Translating the above to scalar diffraction theory, following [16], the point source response for monochrome coherence luminance will be:

Where (x_{img},y_{img}) are the image point coordinate, S_{img} is the image distance. The explicit form of the pupil function is:

P() is the amplitude which is effected by local transmittance and KW() is the phase effected by both lens aberration and phased elements (the latter affecting the “rim-ring zone” only). The system impulse response is super-position of the two “optics” responses:

However, in regular photography, system Eq. (1) describes imaging in incoherent illumination. Hence the matrix columns are the system PSF. Following [16]:

In that case there are potential cross-products between the “lens” and the “rim-ring” and therefore their contributions are not entirely parallel as in field response. Theoretically, the cross products damage the parallelism assumption. However, the nature of the problem is that the “rim-ring” PSF and “lens” PSF work in different zones of the FOV. Thus, the cross-products’ power is relatively low so the parallelism assumption is generally reasonable. To define a limit under which the cross-products can be neglected for filter design, we note that in many cases in engineering a 10% error can be assumed “engineering neglected”. Following that we can define more a rigorous limit to the parallelism assumption in Eq. (39) for all points in the FOV:

For example we take the case in which we use a rim ring with wide PSF that overlaps the lens PSF. Assuming the lens blur is Gaussian with a standard deviation of 0.7 pixel, and the rim-ring response is wide spread Gaussian with a standard deviation of 3 pixel, in Eq. (39) we get 0.088, while for an equal spreading blur twice as wide as the FOV Eq. (39) gives 0.027.

#### 3.5. Lens aberrations in space variant optics

Simple lenses suffer from high aberration and therefore yield blur images. Using the Seidel sums to establish the lens wave front aberration [17] we get:

$$+\genfrac{}{}{0.1ex}{}{1}{2}.S3.\genfrac{}{}{0.1ex}{}{{\mathrm{yp}}^{2}}{{\mathrm{hp}}^{2}}.\genfrac{}{}{0.1ex}{}{{\eta}^{2}}{\eta {max}^{2}}+\genfrac{}{}{0.1ex}{}{1}{4}.\left(S3+S4\right).\genfrac{}{}{0.1ex}{}{{\mathrm{xp}}^{2}+{\mathrm{yp}}^{2}}{{\mathrm{hp}}^{2}}.\genfrac{}{}{0.1ex}{}{{\eta}^{2}}{\eta {max}^{2}}+\genfrac{}{}{0.1ex}{}{1}{2}.S5.\genfrac{}{}{0.1ex}{}{\mathrm{yp}}{\mathrm{hp}}.\genfrac{}{}{0.1ex}{}{{\eta}^{3}}{{\eta max}^{3}}$$

(xp,yp) are the exit pupil coordinates, (hp) is the pupil radius, (η) is the field point distance from the optical axis. The wave front aberration function is common to both “lens” zone and “rim-ring” zone. It yields a space variant point spread function which is different for each field point. The “rim-ring” phase is a summation of the local lens aberrations and the filter phase. The pupil function of the “rim-ring” zone is:

In general, aberrations rise in large angles. On the other hand Eq. (34) assumed paraxial imaging. This gap suggests that in case of a large angle or high accuracy Eq. (34) should be replaced with a more exact model such as presented in [18].

## 4. Determination of the optical filter

As explained above the optical filter **O** is the complementary design of the original system **H**. The two systems together act as a new system **H1** with a better condition. Bellow we will show how to define such a correction **O** filter based on the optics SVD information.

#### 4.1. First design assumption

In general, according to Eq. (37) the system impulse response of **H** and **O** must have cross-products before it becomes PSF in Eq. (38) [19]. As mentioned, one might doubt the parallelism assumption. This assumption might be straight forward in field, but is not obvious in incoherent light intensity. Fortunately, the nature of our problem helps us, because the missing parts of SVD which we add to the original PSF are mainly in regions that do no exist in the original PSF. If the two systems’ impulse responses are spatially separate we can assume that no significant cross-product exists. Hence, our first assumption for **O** filter design is that for each field point the **O** filter PSF is spatially separated from **H** PSF.

#### 4.2. Second design assumption - BMSD concept

To improve the condition of system **H**, **O** should fulfill a few conditions:

1) **O** must be a linear combination of the missing eigen matrices of **H**, the coefficients of the missing eigen matrix are the missing singular values.

2) **O** must be a physical solution, mathematically there are enormous available combinations from which we must find a valid physical solution which yields a positive PSF for each field point.

3) **H1** should be invertible, and hence must include all eigen matrices.

There might be many possible combinations which fulfill the above conditions. All those combinations share the fact that the light reaches certain points and doesn’t reach others. If we put aside the issue of how much light reaches each point, we are left with a common binary map for all valid solutions. The binary map’s relation to **O** design is similar to the relation between a spot diagram and PSF. We call this concept binary matrix spot diagram (BMSD). We suggest to use the BMSD as an initial condition for **O** design.

#### 4.3. Power ratio and lens transmission as a design factor

Extending the original lens with another ring around the original rim, adds more power to the system. The “rim ring” power is proportional to the “rim-ring” area, the “lens” power is proportional to the “lens” area. The ratio between the “lens” area (A_{lens}) and the “rim-ring” area (A_{Rim-ring}) influences the total superposition impulse response in Eq. (37) and Eq. (38). That makes that ratio an important design factor. Another factor for the superposition is the more flexible factor of the apertures transmission for the field. Once setting the geometry one can still control the system transmission. We use Tlens and T_{Rim_ring}. By reducing Tlens we reduce the “lens” zone contribution to the overall PSF shape. The same applies to T_{Rim_ring}. Assuming T_{Rim_ring}=1:

The overall contribution to the total impulse response will be:

Power Rim-ring is the optical power that passes through the “rim-ring” zone (**O**). Power lens is the optical power that passes through the “lens” zone (**H**). Rlens is the “lens” zone radius, ΔRrim is the “rim-ring” radial width.

#### 4.4. The shape of the optical filter phase

In Fig. 3 below we present the SVD information of our study case. Figure 3(a) is a typical singular value graph in which the value decays as the singular value order rises. From the 90th singular value approximately, the numerical values are very low (Figs. 3(a) and 3(b)). If we want to improve the system condition we should focus on improving the low parts of the SVD of **H**, i.e., to include some of the associated eigen images in the BMSD. In our case we chose the last 6 eigen images to compose the BMSD. By defining the BMSD we define the initial condition PSF of **O** for each field point in the system’s FOV. Each BMSD column stands for a different point in the FOV. In Fig. 3(c) we present the PSF shape of the BMSD after translating it back into a 2D image (for convenience, we show only every fourth column). Since there can be no negative value in the PSF, negative values were trimmed. From Fig. 3(c) we see that the typical PSF blur diameter is expected to be in the size of the FOV (10×10 pixels). By calculating second order statistics over the absolute value of each PSF of the FOV, the mean standard deviations of the PSFs’ radius is std=3.1pixels in x direction, std=2 pixels in y direction. Hence in 3 sigma range the PSF radius is of the size of the FOV. The mean distance of the PSF center of mass from its paraxial image in 4.5 pixel with std of 2 pixel, hence to cover all possibilities a single space invariant estimator width should be double of the FOV. For such a wide spread we suggest to use a quadrate filter as a gross estimator of the BMSD. We will calculate the quadrate filter using aberration theory. The fact that the singular values drop sharply is an inherent advantage for the gross estimator, because it insures that even if the parallel system (**O**) effects other zones except the BMSD, the contribution will be weak enough compared to the lower modes power. Hence, in first approximation high spatial precision in not needed. In Figs. 4(a)–4(b) we illustrate the optical system data including the system cross-section, aberration coefficients and the blur geometry for on axis imaging. Assuming on axis imaging, β is the angle of intersection between a ray from the middle of the rim-ring zone and the optical axis in the image plane:

Dlens is the lens diameter. From geometry the relation between Si and R will be:

Rim is the ray aberration which acts as the blur radius. Following [17] the ray aberration is related to the wave front derivative. For the maximal ray aberration we get:

(xp_{rim},yp_{rim}) are the Rim-ring coordinates, n is the refraction index of the surrounding, we get:

Assuming quadratic shape for the “rim-ring” zone:

Deriving the ray aberration for this expression we get:

Applying the maximum value:

Combining Eq. (49) and Eq. (52) we get:

## 5. Study case simulation and results.

Our study case is a simple, bare 0.4 mm pupil diameter plano-convex lens which suffers from space variant aberrations.

#### 5.1. OSLO Lens simulation for Seidel sums

As a study case we took a 0.4 mm diameter micro-lens, ΔRrim=0.0324 mm, hence the overall diameter equals 0.465 mm. The image distance Si=0.69 mm. The Field of view is ±4.673 deg (Fig. 4(c)). The bare lens was first simulated by OSLO EDU version, from that we plotted the Seidel sums which served us in our simulation. In order to have significant space variant blur the original S1,S2,S3,S4 were enlarged by a factor of ten and more. The Seidel sums values are listed in Fig. 4.

#### 5.2. Initial conditions values for the simulation

FOV=0.113 mm, Tlens=1, T_{Rim_ring}=1, For Sshape calculation we use ΔRrim=0.03 mm, Rlens=0.2mm,Si=0.69 mm, from Eq. (42) PR=0.324, from Eq. (46) β=17.3°, from Eq. (54) Sshape1=2.235µm, pixel size 11.3µm, image size 10×10, n=1, Rim=0.113 mm. To present our correction method capabilities we first make the system matrix more ill. We do that by de-centering the imager field of view relative to the optical axis by -0.6 pixel in x direction and +1.4 pixel in y direction. Compromise between sampling consideration and computing resources brought the simulation FFT stage to be 512×512.

#### 5.3. Simulation log

After determining the system geometry and the filter phase we start the simulation. The below log in Table 1 presents the progress towards improvement in the condition number. In line 1 we see the system properties before adding the “rim-ring”. In line 2, after adding the “rim-ring” and applying initial conditions we get an instant improvement by factor of 35.7. From line 3 to line 8 we try to maximize the improvement in the condition number by influencing the PSF superposition according to Eq. (38), Eq. (45). Since the system condition rises sharply when reducing Tlens, we fix Tlens=1 and start tuning T_{Rim_ring}. We get the best results using T_{Rim_ring}=0.9. In lines 6 to 9 we scan some values of Sshape. Results show that Eq. (54) predicts the best value in this neighborhood of solutions space. Finally we see that line 4 values yield a clear minimum of the condition number. The overall improvement of the system condition number is by a factor of 38.3. Additional information in Table 1 refer to the “matrix sum”. We use PSF matrix summation in order to calculate the amount of transferred illumination to the image plane. The power transmittance and effective NA comparison between the systems with and without the “rim ring” (the original simple lens and the parallel optics system) is done in chapter 6. In Fig. 5 below we compare the immunity to noise of the original simple lens with that of the new parallel optics system. The systems were subject to additive white Gaussian noise, the image restoration was done based on simple matrix inversion Eq. (55). The performance rate was the value of the Mean Square Error Improvement Factor (MSEIF) Eq. (56) [13]. In this function the restoration error and the blur error are both relative to the ideal object (I^{object}). When MSIEF<0 db the restored image (I^{res}) is worse then the blur image (I^{img}) so there is no use in restoration. For each SNR we perform a 10,000 measurements ensemble, we calculate the MSIEF for each measurement and present the average value over the ensemble.

Results show that the new system immunity to noise in this case is significantly better with a gap of 50db in SNR before reaching the 0 db equilibrium. In Fig. 6 below we present the imaging results in two SNR levels. In SNR=133db the results MSEIF are 0.236db and 40.7db for the simple lens and the parallel optics respectively, and in SNR=100 db are (-)8.15db and 12.3db respectively. These results suggest significant improvement in noise immunity of the parallel optics design relative to the original simple lens.

## 6. Ensemble of results and comparison with the MTF approach

In this section we provide some simulation results for images created using our proposed system and a comparison to a typical MTF approach. In previous hybrid imaging systems (optical and image processing) the strategy was to find a way to design a phase mask that makes the system’s OTF independent of the aberration parameter. The typical result is an invariant OTF which has a somewhat narrow bandwidth. Next, a sharp image was achieved by image post processing. Authors used several approaches to define the needed phase mask. In [20] a phase mask was defined for extending the depth of field. That was done by solving the stationary phase approximation of the defocus related ambiguity function. Another approach to define a phase mask for extending the depth of field was by using the Fermat principle [21]. In [22] a phase mask was made for reducing defocus and spherical aberration dependency, by solving the on-axis intensity stationary phase approximation. In a recent work [23] an effective quartic phase retardation function (QF) of 2^{nd} and 4^{th} order was suggested to be used for reducing the variation of intensity response in the center of various field points. In result, the Sterhl ratio increases and the modulation transfer function (MTF) becomes practically invariant. The work shows remarkable results for a moderate amount of space variant coma or astigmatism aberrations (up to 3λ) separately. The principle is to add the system QF that is optimized to cancel the dependency of the PSF intensity (stationary phase approximation) on the space variant aberration term. Before Applying this model in our case we must emphasize two major differences between our case and the case in [23]. First, our case is a 5 term aberration polynomial, while in [23] it is a single kind residual aberration. Second, we focus on high aberration while [23] solves a moderate aberration. In order to implement the QF phase in our case we use the fact that in some conditions the filter parameter for both coma and astigmatism are the same [23]. Following that we define a QF in Eq. (57). r is the normalized pupil radius, coefficients are in [mm].

Figure 7 presents the average MSEIF Eq. (56) results for an ensemble of 11 different objects Fig. 8, (column 3). Results include the “lens”, the “rim-ring”, and the “QF” restorations. Following [22] we assume that the QF achieved a space invariant response, and so for the QF image restoration we can use the on axis OTF for Wiener filtering [24]. In order that Eq. (56) results reflect the improvement relative to the “lens” imaging we use the “lens” image to calculate the MSEIF nominator for the three restorations. In Fig. 8 we present representative ensemble imaging simulation in SNR=130 [db]. The results show that the ensemble MSEIF behaves very much like in the MSEIF in Fig. 5. In Fig. 7 we also see that the QF image restoration MSEIF has a low value and is nearly constant as a function of the SNR, i.e. the Wiener restoration is less sensitive to noise than the matrix inversion but yields low restoration due to the space variant nature of the problem which is still unsolved by the QF filter for our model. The space variant behavior which is still present after the QF can be seen very clearly in the “Checkerboard” QF image Fig. 8 in line 10 column 7. In general, the results in Fig. 8 show that the QF restoration appearance (column 6) is not significantly different from the simple lens image (column 1). In fact its dynamic range which is not seen is broaden to hundred below and above the 8 bit gray scale restoration kept in the “lens” and “rim-ring” restorations. As a result of the broadened dynamic range the norm calculation in the MSEIF denominator is very large compared to the nominator and the final result is a very low QF MSEIF. On the other hand, comparing the images before restoration, Fig. 8, (column 1, column 5, column 7) it seems that the QF image, far as it is from the object, is the most space invariant. This result is expected since in the QF approach the optical filter comes to reduce the PSF broadening.

## 7. Power transmittance and effective numerical aperture

Illumination transmission to the image is an important property of the system. Several factors influence the image plane illumination. By adding the “rim-ring” aperture to the original system we add more power to the image. By reducing Tlens and T_{Rim_ring} we reduce the power transmission to the image. Bellow we investigate this important issue. By summation over the PSF matrix we can get a basis for comparing the amount of power which each of the systems transfer to the FOV. The three important references for power transmittance are P_{lens} (power of the original system) P2 (power of the system with the “rim-ring” when the lens are in full transmission), P_{parallel_opt} (power of the final parallel optics system). Aberrations in a single lens optical system might be prevented by reducing the system NA, therefore, one of the main advantages of the above method is that it allows producing simple optics with relative high NA. We can calculate the ratio between the new system power and the original system power as follows:

Power is related to the square of the aperture, so the effective parallel optics aperture diameter can be calculate as:

The effective NA of the optical system is:

This result is very encouraging because it shows that our suggested method allows an effectively high NA optical system. As part of our design we assume that the cross products values are low, so the “lens” and “rim-ring” are effectively parallel. In a somewhat different approach than Eq. (39) we can perform a global check of this assumption with the simulation data in Table 1. Knowing the power of the “lens” zone and the “lens” with the “rim-ring”, all in full transmission, the transmitted power of the “rim-ring” into the FOV equals P2-Plens. Assuming it is constant we expect the power which does transmit to the FOV to be:

If now we compare Eq. (60) to the sum(matrix) in Table 1, the difference will be the cross-product:

Calculating P_{cross_product} for each line of Table 1, the maximal value equals 3.18% of the transmitted power, hence the parallel assumption in our case holds.

## 8. Summary and conclusions

Matrix condition number is one of the must powerful figure to measure system immunity to noise. We suggested a new approach in which the optical design’s aim is to improve the matrix condition. In the case studied we showed a significant improvement in the matrix condition number. We showed that the SVD of the system’s matrix can direct us to the missing parts needed for improving system condition. We showed that those missing parts can be interpreted as a parallel optical system. Further we suggested an approach for such filter design for a case where a wide spread parallel filter PSF is needed. Finally we simulated a study case and compared the condition number and noise immunity in both original and parallel optics systems. The simulation results showed a significant improvement in noise immunity in the process of image restoration by matrix inversion. In this study case the new proposed design crosses the MSEIF=0 db threshold in a SNR level which is more than 50 db lower than the SNR level in which the original simple lens system crosses that threshold. Comparing the new model study case simulation results with those of the QF model, results are to the favor of the new model. Another improvement is in the system NA. The new system NA equals 0.35 and is higher than the original system NA witch equal to 0.28. Based on these results it seems possible to design a high NA hybrid optical system based on simple optics and image processing.

## Acknowledgment

We would like to thank Professor Uri Shaked of Tel Aviv University Electrical Engineering School, for helpful discussions on matrix inversion and least square problems. In Fig. 4(c) we used an OSLO EDU version for calculating Seidel sums and the optical system cross-section respectively.

## References and links

**1. **G.M. Robbins, “Inverse Filtering for linear shift variant imaging system”, Proc. IEEE , **60**,862–872 (1972).
[CrossRef]

**2. **H. Trussell and B. Hunt, “Image restoration of space-variant blurs by sectioned methods,” IEEE. Trans. ASSP. **26**608–609 (1978).
[CrossRef]

**3. **H. Trussell and B. Hunt, “Image restoration of space-variant blurs by sectioned methods,” Proc. IEEE. ASSP ,**3**,196–198 (1978).

**4. **T.P. Costello and W.B. Mikhael, “Efficient restoration of space variant blurs from physical optics by sectioning with modified wiener filtering,” Digital and image processing ,**13**,1–22 (2003).
[CrossRef]

**5. **J.M. Varah, “On the numerical solution of ill-conditioned linear system with applications to ill posed problems,” SIAM J. Numer. Anal ,**10**, 257–265 (1972).
[CrossRef]

**6. **J.H. Wilkinson, “Rounding Errors in Algebraic Processes,” 91–93 (Her Majesty’s stationery office, 1963).

**7. **G.H. Golub and C.F. Van-loan, “Matrix Computation,” 17–185 (North oxford academic,1983).

**8. **C.L. Lawson and R.J. Hanson, “Solving Least Squares Problems,” SIAM J. Numer. Anal, 185–8 (1995)

**9. **I.F. Gorodnitsky and D. Rao, “Analysis of Error Produced by Truncated SVD and TIkhonov Regularization Methods, Conference Record of the Twenty-Eighth Asilomar Conference on Signals, Systems and Computers,” **1**, 25–29 (1994).

**10. **C.M. Leung and W.S. Lu, “An L curve approach to optimal determination of regularization parameter in image restoration,” Proc. IEEE,1021–1024 (1993).

**11. **X. Wang, “Effect of small rank modification on the condition number of matrix,” Computers and mathematics with applications **54**,819–825 (2007).
[CrossRef]

**12. **S. Twomay, “Information content in remote sensing,” Appl. Opt , **13**, 942–945 (1974).

**13. **M. Bertero and P. Boccacci, “Introduction to inverse problems in imaging,” **86**,252 (IOP,1998).

**14. **A. A. Sawchuk and M. J. Peyrovian, “Restoration of astigmatism and curvature of field,” J. Opt. Soc. Am. **65**, 712–715 (1975).
[CrossRef]

**15. **H.C. Andrews and C.L. Paterson, “Singular value Decomposition and digital image processing”, IEEE. Trans. ASSP. **24**, 26–53 (1976).
[CrossRef]

**16. **J.W. Goodman, “Introduction to Fourier Optics,” (Mcgraw-Hill, 1996).

**17. **W.T. Welford, “Aberrations of Optical Systems,” (Adam-Hilger, 1991).

**18. **V. Shaoulov et al, “Model of wide angle optical field propagation using scalar diffraction theory,” Opt. Eng , **43**, 1561–1567 (2004).
[CrossRef]

**19. **A.W. Lohmann and W.T. Rhodes, “Two pupil synthesis of optical transfer function,” Appl. Opt. **17**, 1141–1146 (1978).
[CrossRef] [PubMed]

**20. **E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. **34**, 1859–1866 (1995).
[CrossRef] [PubMed]

**21. **W. Chi and N. George, “Electronic imaging using a logarithmic asphere,” Opt. Lett. **26**, 875–877, (2001).
[CrossRef]

**22. **S. Mezouari and A.R. Harvey, “Phase pupil function for reduction of defocus and spherical aberration,” Opt. Lett. **28**, 771–773, (2003).
[CrossRef] [PubMed]

**23. **S. Mezouari et al, “Circularly symmetric phase filter for control of primary third order aberration: coma and astigmatism,” J, Opt,.Soc.Am.A. **23**,1058–1062 (2006).
[CrossRef]

**24. **N.S. Kopeika, “A system Engineering approach to imaging,” 517–520 (SPIE, 1998).