## Improvement of matrix condition of Hybrid, space variant optics by the means of Parallel Optics design

Optics Express, Vol. 17, Issue 14, pp. 11673-11689 (2009)

http://dx.doi.org/10.1364/OE.17.011673

Acrobat PDF (673 KB)

### 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

## 2. Introduction to space variant image processing

**A**, the columns of which are the optics Point Spread Function (PSF) for each field point. We get than:

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

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]

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]

### 2.2. Important matrix properties for numerical solution

_{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**, 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

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]

_{k}) involved in the solution. Higher eigen components between k<i<n are truncated. The description of

**b**will then be:

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]

**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:

**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:

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

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

### 3.1. Matrix representation of a system

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

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

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

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

^{½}. σ

_{i}are positive square roots of Δ. The order of σ

_{i}values is high to low [7]:

_{i}, there exists an eigen matrix M [15

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

_{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:

### 3.3. Constraints

### 3.4. Point spread function of parallel optics

_{img},y

_{img}) are the image point coordinate, S

_{img}is the image distance. The explicit form of the pupil function is:

### 3.5. Lens aberrations in space variant optics

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

## 4. Determination of 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

**H**and

**O**must have cross-products before it becomes PSF in Eq. (38) [19

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

**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

**H**,

**O**should fulfill a few conditions:

**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.

**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.

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

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

_{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:

**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

**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:

_{rim},yp

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

## 5. Study case simulation and results.

### 5.1. OSLO Lens simulation for Seidel sums

### 5.2. Initial conditions values for the simulation

_{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

_{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.

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

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]

## 7. Power transmittance and effective numerical aperture

_{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:

_{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

## Acknowledgment

## References and links

1. | G.M. Robbins, “Inverse Filtering for linear shift variant imaging system”, Proc. IEEE , |

2. | H. Trussell and B. Hunt, “Image restoration of space-variant blurs by sectioned methods,” IEEE. Trans. ASSP. |

3. | H. Trussell and B. Hunt, “Image restoration of space-variant blurs by sectioned methods,” Proc. IEEE. ASSP , |

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 , |

5. | J.M. Varah, “On the numerical solution of ill-conditioned linear system with applications to ill posed problems,” SIAM J. Numer. Anal , |

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,” |

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 |

12. | S. Twomay, “Information content in remote sensing,” Appl. Opt , |

13. | M. Bertero and P. Boccacci, “Introduction to inverse problems in imaging,” |

14. | A. A. Sawchuk and M. J. Peyrovian, “Restoration of astigmatism and curvature of field,” J. Opt. Soc. Am. |

15. | H.C. Andrews and C.L. Paterson, “Singular value Decomposition and digital image processing”, IEEE. Trans. ASSP. |

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 , |

19. | A.W. Lohmann and W.T. Rhodes, “Two pupil synthesis of optical transfer function,” Appl. Opt. |

20. | E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. |

21. | W. Chi and N. George, “Electronic imaging using a logarithmic asphere,” Opt. Lett. |

22. | S. Mezouari and A.R. Harvey, “Phase pupil function for reduction of defocus and spherical aberration,” Opt. Lett. |

23. | S. Mezouari et al, “Circularly symmetric phase filter for control of primary third order aberration: coma and astigmatism,” J, Opt,.Soc.Am.A. |

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

**OCIS Codes**

(070.6110) Fourier optics and signal processing : Spatial filtering

(080.1010) Geometric optics : Aberrations (global)

(100.3190) Image processing : Inverse problems

(110.3010) Imaging systems : Image reconstruction techniques

**ToC Category:**

Image Processing

**History**

Original Manuscript: October 2, 2008

Revised Manuscript: April 1, 2009

Manuscript Accepted: April 12, 2009

Published: June 26, 2009

**Citation**

Iftach Klapp and David Mendlovic, "Improvement of matrix condition of Hybrid, space variant optics by the means of Parallel Optics design," Opt. Express **17**, 11673-11689 (2009)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-14-11673

Sort: Year | Journal | Reset

### References

- G. M. Robbins, "Inverse Filtering for linear shift variant imaging system," Proc. IEEE, 60, 862-872 (1972). [CrossRef]
- H. Trussell and B. Hunt, "Image restoration of space-variant blurs by sectioned methods," IEEE.Trans. ASSP 26, 608-609 (1978). [CrossRef]
- H. Trussell and B. Hunt, "Image restoration of space-variant blurs by sectioned methods," Proc. IEEE. ASSP 3, 196-198 (1978).
- 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 Processing13, 1-22 (2003). [CrossRef]
- 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]
- J. H. Wilkinson, " Rounding Errors in Algebraic Processes," 91-93 (Her Majesty’s stationery office, 1963).
- G. H. Golub and C. F. Van-loan, "Matrix Computation," 17-185 (North oxford academic, 1983).
- C. L. Lawson and R. J. Hanson, "Solving Least Squares Problems," SIAM J. Numer. Anal, 185-8 (1995)
- 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).
- 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).
- X. Wang, "Effect of small rank modification on the condition number of matrix," Computers and mathematics with applications 54, 819-825 (2007). [CrossRef]
- S. Twomay, "Information content in remote sensing, " Appl. Opt, 13, 942-945 (1974).
- M. Bertero and P. Boccacci, "Introduction to inverse problems in imaging," 86, 252 (IOP,1998).
- A. A. Sawchuk and M. J. Peyrovian, "Restoration of astigmatism and curvature of field, " J. Opt. Soc. Am. 65, 712-715 (1975). [CrossRef]
- H. C. Andrews and C.L. Paterson, "Singular value Decomposition and digital image processing," IEEE. Trans. ASSP. 24, 26-53 (1976). [CrossRef]
- J. W. Goodman, "Introduction to Fourier Optics," (Mcgraw-Hill, 1996).
- W. T. Welford, "Aberrations of Optical Systems," (Adam-Hilger, 1991).
- V. Shaoulov et al, " Model of wide angle optical field propagation using scalar diffraction theory," Opt. Eng, 43, 1561-1567 (2004). [CrossRef]
- A. W. Lohmann and W. T. Rhodes, "Two pupil synthesis of optical transfer function," Appl. Opt. 17, 1141-1146 (1978). [CrossRef] [PubMed]
- E. R. Dowski and W. T. Cathey, "Extended depth of field through wave-front coding," Appl. Opt. 34, 1859-1866 (1995). [CrossRef] [PubMed]
- W. Chi and N. George, "Electronic imaging using a logarithmic asphere," Opt. Lett. 26, 875-877, (2001). [CrossRef]
- S. Mezouari and A. R. Harvey, "Phase pupil function for reduction of defocus and spherical aberration," Opt. Lett. 28, 771-773, (2003). [CrossRef] [PubMed]
- 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]
- N. S. Kopeika, "A system Engineering approach to imaging," 517-520 (SPIE, 1998).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.