## Abstract

In this paper the invariance of modulation transfer function (MTF), which describes the insensitivity to perturbation of MTF, is defined to be the evaluating criterion of the wavefront coding system. The rapid optimization of wavefront coding system based on the MTF invariance is proposed by means of introducing the mathematical program Matlab to normal optical design process. The interface called MZDDE between Matlab and Zemax is applied to realize the fast data exchanging and merit function calculating. The genetic algorithm tool (GA) in Matlab is introduced to the optimizing process, which accelerates the converging efficiency considerably. The MTF invariance of optimized system drops to 0.0119 while that of original system is larger than 0.018. If the all the fields of view is taken into consideration, the MTF invariance of optimized system and original system is less than 0.015 and larger than 0.020 respectively. It is proven that the optimization of the unusual optical system with special property can be executed conveniently and rapidly with the help of external program and dynamic data exchange.

©2009 Optical Society of America

## 1. Introduction

The wavefront coding technology (WFC) is known as a computational imaging technique for imaging system which has been widely applied in many fields such as iris recognizion, bar-code scanning, biological microscopy, thermal imaging system, etc. One of the most significant properties of the WFC technology is to control the misfocus and misfocus related aberrations. The depth of focus (DOF) of imaging system can be extended by a factor of ten or more compared to the traditional system with special coding component or surface when the F# keeps constant. Although the direct image obtained by the optical system is degraded (encoded) but nearly remains the same in the extended DOF, the final clear image can be obtained by digital post-processing (decoding) [1–6]. Research on restoration of images by wavefront coded system has also been carried out. Bradbum et al. [7] introduced the traditional least square method. Based on their work, van der Gracht et al. [8] proposed the Modified Residual Norm Steepest Descent (MRNSD method) and excellent result was obtained. In our previous work the restoration based on regularization and simulated point spread function is also proposed [9].

Different from traditional imaging system, the wavefront coding system is not designed to achieve the best image quality but to preserve the image quality as constant as possible in the extended DOF and reduce the difficulty of digital processing. Therefore the optimizing target of wavefront coding system also becomes different from traditional system. The leading factor lies in that proper quality criterion that exactly describes the coincidence of MTF in the extended DOF should be established and applied in the optimizing process. Several criteria have been proposed to optimize and evaluate the wavefront coding system. Fisher information (FI) is the most widely- used evaluating function as a measure of the sensitivity of point spread function (PSF) to defocus of the wavefront coding system [10,11]. Quality factor (QF), which is the integral of the defocused MTF with respect to the spatial frequency has been built as a MTF based criterion to describe the correlation of variation of the MTF with the defocus [12]. Strehl ratio (SR) is also considered as effective evaluation of the wavefront coding system with rotational symmetry [13,14]. Besides, some designers use the combination of SR and encircled energy (EE, which refers to the percentage of the total energy in the PSF contained within a radius measured from the PSF center) as their criterion in the design of triplet objectives [15]. It can be concluded that all the criteria mentioned above are all MTF based directly or indirectly.

Hence, these criteria or merit function can be obtained through simple calculation and enables rapid optimization and convergence.

Theoretically speaking, a perfect wavefront coding system should satisfy the relationship shown in Eq. (1) for given toleranced error budget of defocus aberration${w}_{020}$:

Where m denotes the different position of image plane in the extended DOF corresponding to different defocus value ${w}_{020}$and n denotes the sampled space frequency which ranges from 0 to ${f}_{cut-off}$ and is spaced by proper increment [15].Equation (1) demonstrates the most fundamental and important property of wavefront coded system: MTF values or image quality should be uncorrelated to the defocus aberration within the extended DOF. In another words the MTF values will remain perfectly uniform within the extended DOF if the relationship is fulfilled for a wavefront coded system. Thus, Eq. (1) is actually considered to be the most appropriate criterion or merit function of the optimization since it’s the essential character of the wavefront coding system defined by its fundamental theory. However, the practical application of that criterion means that all the $\partial MTF/\partial {w}_{020}$ values of each sampled space frequency in each position of image plane in the extended DOF will be calculated within each optimizing cycle. It will be of extremely enormous consumption of time and computation for optical design program, which limits the actual application of the merit function in Eq. (1).

In this paper an innovative optimizing method is proposed based on the MTF invariance deduced from the merit function in Eq. (1) to improve the coincidence of MTF within the extended DOF and accelerate the optimizing process with less time and computational resource by introducing Matlab and effective genetic algorithm. The widely-used mathematical program Matlab is introduced to the optical design software program, which provides with friendly interface to other Windows program and enables data link and exchange conveniently. The optimization can be executed in Matlab environment with the MTF value exported from Zemax, which just serves as a MTF calculator. Since the Matlab programs are compiled, they run significantly faster, from which the complex calculations of the merit function benefits dramatically. The total optimizing time is further shortened by means of the application of GA tool in Matlab. Although the accurate analytical expressions of the MTF of wavefront coded systems has been developed and applied as criteria in the optimization of such systems [3], the numerical approach is applied in our design process rather than analytical form of MTF. Firstly, the MTF value obtained by Zemax is digital not continuous corresponding to each sampled spatial frequency, which enables the numerical calculation of MTF invariance more conveniently and directly. Secondly, the analytical expression of MTF may be more complicated if extra parameters are introduced into the optimizing process, besides the analytical form of MTF is not concerned in the optimization (even in the digital image restoration). Finally, the numerical approach is easier to execute. Therefore, the numerical approach is adopted.

## 2 Work foundations

The research is based on a typical Cook’s three-mirror anastigmatic (TMA) system (invented by Wetherell and Womble of ITEK Corp., US Patent 4240707 in Dec. 23, 1980) for compact space-borne camera. The layout of the system and its main parameters are shown in Fig. 1 . All the mirrors have been finished and we have assembled the system since May 2009. Because the WFC technology can control the misfocus and misfocus related aberration and loosen the tolerance error budget of fabrication and alignment, we attempt to apply the technology on the off-axis TMA system shown in Fig. (1) as a new method to maintain its focal plane in the DOF and keep the image quality nearly invariable when there is mechanical deformation and thermal gradient [16]. In our previous work presented in ref. 13 and 14 the basic designing of TMA system with WFC technology and the image restoration of encoded system has been researched and in this paper an optimizing method of this kind of system is proposed based on the basic design to improve its insensitivity to defocus aberration.

The secondary mirror is selected to be the wavefront coding element for it serves as the stop aperture and has relatively small geometrical size. The surface sag of the secondary mirror before optimization is shown in Eq. (2) [16]

where x, y is mirror coordinate, z is surface sag, c is the curvature of the vertex, k is the conic constant and *β* denotes the magnitude of cubic terms.

where x, y is normalized coordinate and *α* indicates the magnitude of cubic phase term from which *β* in Eq. (2) is deduced.

The actual cut-off space frequency of the system is about 57 lp/mm, which is determined by the Nyquist frequency of digital image detector (the size of CCD pixel is$8.75\mu m\times 8.75\mu m$). Considering authentic restoration of mid images, the dynamic MTF must be higher than 0.05 to decrease the difficulty of image restoration, which demand the design MTF should be 0.18 at least because of the push-broom working mode (In a push broom sensor, a line of sensors arranged perpendicular to the flight direction of the spacecraft is used, thus different areas of the surface are imaged as the spacecraft flies forward), residual manufacturing and alignment error and other factor that may decrease the MTF value. Thus, the value of *α* is set about 24 to satisfy the requirement although the value of *α* is suggest to be set much larger than 20 to ensure adequate insensitivity to misfocus of the MTF according to the theory of WFC technology. Therefore, the MTF curves differ from each other a little especially in the low-frequency domain, and then further compensation and optimization can be carried out to improve the performance.

## 3. Definition of MTF invariance

The most distinct character of wavefront coding system is the insensitivity of the MTF to various defocus aberration in its given range. The MTF value of arbitrary spatial frequency should maintain nearly the same in the extended DOF. MTF invariance is defined in Eq. (4) to describe the degree of consistency of MTF in the extended DOF.

The $P{V}^{MTF}$ gives the largest departure of the $MT{F}^{test}$from$MT{F}^{ref}$, while the $RM{S}^{MTF}$ describes the average deviation between the $MT{F}^{test}$from $MT{F}^{ref}$during the whole region of space frequency, which is a more important parameter to characterize the invariance. The $RM{S}^{MTF}$is applied as the optimizing criterion in the next section and the $P{V}^{MTF}$is used as verification of the optimization result.

## 4. Definition of optimizing criterion base on MTF invariance

The optimizing criterion is established in Eq. (5):

Where$RM{S}_{j}^{MTF}={[(\underset{1}{\overset{n}{\Sigma}}\Delta MT{F}_{i}^{2})/n-{(\underset{1}{\overset{n}{\Sigma}}\Delta MT{F}_{i}/n)}^{2}]}^{1/2}{|}_{i=1\cdots n,j=1\cdots m}$ $\Delta MT{F}_{i}{|}_{j=1\cdots m}=MT{F}^{test}{}_{i}-MT{F}^{ref}{}_{i}{|}_{i=1\cdots n,j=1\cdots m}$ $j=1,\mathrm{2......}m$, m is the number of sampling positions where the $MT{F}^{test}$is calculated corresponding to different defocus value. The value of m should be assigned properly to balance the tradeoff between calculation speed and sampling intensity. When the value of Eq. (5) is minimized, the optimal design is considered to be achieved. It can be observed in Fig. 4 that the $RM{S}^{MTF}$curve of the original system ($\beta =8\times {10}^{-4}$)with different defocus aberration in the extended DOF is not continuous or monotonous but do not oscillate fiercely or vary sharply, it is enough to set m equal 11 to get sufficient sampling intensity with the increment of $0.25\lambda $ defocus aberration. It also can be seen from Fig. 4 that the maximum $RM{S}^{MTF}$value of the original system is about 0.0187. As the assignment of the value of m, the same consideration should be taken into in setting the value of n. Finally, the interval of sampling space frequency is 0.5 lp/mm so that n is set 115.This criterion is based on the consideration that if the maximum value of$RM{S}^{MTF}$in the extended DOF is minimized, all the $RM{S}^{MTF}$value will also achieve its minimum. Thus, the optimizing cycles are concentrated on minimizing the maximum value of the 11 sampling$RM{S}^{MTF}$.

The surface sag of the secondary mirror is transformed to the optimizing form as shown in Eq. (6):

It is presented above that the DOF of the wavefront coding system under research is extended as ten times as traditional system and the image quality will remain the same in the extended DOF. However, the concept of Gaussian image plane of traditional system is not suitable for the wavefront coding system or in another words the extended DOF of wavefront coding system may not be symmetric about the original “Gaussian image plane” for optimum image quality coincidence mainly because of the residual off-axis aberration. A nominal Gaussian image plane can be defined for describe the location of the extended DOF for convenience as shown in Fig. 5 .

Furthermore, in the unoptimizable system the location of the extended DOF is fixed and just symmetric about the original Gaussian image plane. Nevertheless for the system to be optimized the location of the extended DOF is made another variable ${a}_{5}$to gain the most appropriate position, which means the range of the extended DOF can be moved as a whole to obtain the global minimum value of MTF similarit. In the unoptimizable system the position of $MT{F}^{ref}$is just set in the middle of the extended DOF, which is namely the Gaussian image. But for the system to be optimized, the position of Gaussian image may not the optimal choice to obtain the minimum of$RM{S}_{MTF}$. Hence the position where the$MT{F}^{ref}$is calculated is also made variable ${a}_{6}$to obtain the most suitable $MT{F}^{ref}$for minimizing the $RM{S}_{MTF}$, which means that the position of$MT{F}^{ref}$can be select within the whole extended DOF in optimization. Thus, the optimizing variable vector is established as$A={[{a}_{i}]}_{i=6}$. Several restrictions should be imposed on the variable to avoid invalid results. Firstly, ${a}_{1}$and${a}_{4}$should be limited a little smaller than the initial value of *β* to maintain the fundamental property of the system. Secondly, the maximum value of${a}_{2}$and ${a}_{3}$should be limited less than one tenth of *β* for the same purpose. Thirdly, the absolute values of ${a}_{5}$ should be restricted less than 3 to limit the deviation of the referenced image plane from original Gaussian image plane not too far away. Finally, ${a}_{6}$should be limited properly so as to maintain the reference position where the $MT{F}^{ref}$is calculated within the bound of the extended DOF. All the restrictions are appended to prevent the nonsensical result which achieves the small $RM{S}_{MTF}$ at the cost of excessively low MTF value.

## 5. Optimization of the wavefront coding system with MZDDE and GA tool

The optimizing criterion is defined in Eq. (4) as the merit function. Therefore the calculation of this merit function becomes the key problem. It can be seen that the MTF value of totally 11 particular positions with the give sampling space frequency is computed and the compound calculation of summing, subtracting and root-mean-square computation is performed to obtain the value of merit function in each optimizing cycle. However, Zemax is not good at dealing with this kind of very complicated computation. Of course the merit function can be founded by applying the user defined operand (UDO) programmed with Zemax programming Language (ZPL) and calculated by Zemax directly, but it will be of extreme consumption of time and computing resource, which makes the optimizing process badly insufficient and completely unacceptable. Therefore, a different effective method must be sought to execute the optimization. Fortunately, Zemax has a very powerful feature which allows another Windows program to establish a communication link with Zemax and to extract data from Zemax through DDE (Dynamic Data Exchange) for further analysis or computation, where Zemax just serves as a ray-tracing tool. This feature indicates that the optimization can be performed with the help of expert program which is more effective in treat with complex computation. Thus the widely used mathematical program Matlab is selected to perform the optimization to accelerate the converging speed. The interface between Zemax and Matlab called MZDDE (Zemax DDE Toolbox for Matlab), which is published by Gerek Griffith at CSIR, provides with the convenient and fast data exchange between the two programs and enables the operation of Zemax in Matlab environment [17]. Hence, the compound computation and optimizing cycles is executed in Matlab program and the MTF of 12 particular positions is calculated in Zemax taking advantage of its mature and rapid algorithm. The operating step in one optimizing cycle is shown as follows:

- 1 Create an optimizing vector $A={[{a}_{i}]}_{i=6}$ in Matlab.
- 2 Push these data to Zemax and make them take effect through MZDDE.
- 3 Calculate the MTF of the 12 particular positions with the given sampled space frequency in Zemax.
- 4 Extract the MTF value from Zemax through MZDDE and create the variable vector of Matlab according to the MTF value.
- 5 Calculate the value of opt_cri for comparison and filtration.

It should be noticed that the Zemax file to be optimized should be opened and the DDE link must be established before the optimization is executed.

Of course only one cycle is not enough for the minimum-seeking problem; thereby a proper algorithm is needed to select the optimum solution in thousands of cycles. Zemax uses an actively damped least squares method in optimization. For the problem presented in this paper, the GA tool is utilized as the optimizing algorithm because the optimization is performed in Matlab environment and the embedded tool of Matlab is convenient to invoke, which can further accelerate the optimizing process. The genetic algorithm (GA) is a method for solving optimization problems that is based on natural selection, the process that drives biological evolution. The genetic algorithm repeatedly modifies a population of individual solutions. At each step, the genetic algorithm selects individuals at random from the current population to be parents and uses them produce the children for the next generation. Over successive generations, the population “evolves” toward an optimal solution [18,19]. The GA can be applied to solve a variety of optimization problems that are not well suited for standard optimization algorithms. Here GA tool in Matlab is used directly to handle the problem under research. The fitness function or objective function is founded by converting the optimizing cycle mentioned above to M-file properly. The number of independent variables is 6. The total generation is 200, the population of which is 50. The final optimized results are: a_{1} = 7.99958 × 10^{−4}, a_{2} = −5.45392 × 10^{−4}, a_{3} = −7.56767 × 10^{−4} and a_{4} = 7.95896 × 10^{−4}. It can be observed from Fig. 6
that the maximum $RM{S}_{MTF}$value drops to less than 0.012 compared with the original value of larger than 0.018 and all the restrictions added ahead are all satisfied.

It can be seen from Fig. 7 that these MTF curves seem to coincide with each other not so closer as the original system shown in Fig. (2) left. But in fact t the $P{V}^{MTF}$of optimized system drops to 0.046 from 0.057 of original system since the selection of$MT{F}^{ref}$is also optimized. The MTF value has nearly no reduction compared with the original system. The most distinguished advantage of this optimizing method is much more time-saving than the optimization directly operated in Zemax. The total execution time of the optimizing process is about 30 minutes. Nevertheless the direct optimization with the same criterion merit function with Zemax is very difficult to operate, which always terminate abnormally without any significative result mainly because the system resource is exhausted by the massive computation. And it can be seen from the Fig. 8 that the optimization result is also valid for all the fields of view (FOV).

The$RM{S}_{MTF}$of the optimized system of all the other FOV is less than 0.015 while the corresponding$RM{S}_{MTF}$of the original system is larger than 0.020. It can be observed that the $RM{S}_{MTF}$of the other FOV is not as fine as the one of the central FOV used in optimization, which is obvious because the $RM{S}_{MTF}$of the other FOV is not considered adequately in optimization but just supposed to be minimized so long as the $RM{S}_{MTF}$of central FOV approaches its optimum instead. The assumption is applied practically although it is not absolutely exact. After all the optimizing result can be accepted without insufferable deviation and the dominating excuse is that the cost time will be extremely prolonged based on the existing criterion if the FOV ingredient is involved. The improved optimizing criterion including the FOV and requiring less computational time and resource is under research

## 6. Conclusions

In this paper the MTF invariance is defined as a new criterion for the optimization of wavefront coding system. The new optimizing method based on MZDDE and GA tool is proposed to deal with the complex computation of MTF invariance. The optimizing result and the cost time are satisfied. In the future work the factor of FOV is planned to be introduced to the optimization and the more suitable result for all FOV may be obtained. Another important work calls for great attention and effort lies in the manufacturing and testing of the optimized secondary mirror with adequate accuracy, which is actually a free-form surface.

## Acknowledgement

This research is supported by a grant from the National High Technology Research and Development Program of China (863 Program) (No. 2009AA12Z105).

## Reference and Links:

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

**2. **J. Ojeda-Castañeda, J. E. Landgrave, and C. M. Gómez-Sarabia, “Conjugate phase plate use in analysis of the frequency response of imaging systems designed for extended depth of field,” Appl. Opt. **47**(22), E99–E105 (2008). [CrossRef] [PubMed]

**3. **S. Bagheri, P. E. X. Silveira, and D. Pucci de Farias, “Analytical optimal solution of the extension of the depth of field using cubic phase Wavefront Coding,” J. Opt. Soc. Am. **25**(5), 1051–1063 (2008). [CrossRef]

**4. **W. Chi, and N. George, “Smart Camera with Extended Depth of Field,” Proc. SPIE **6024**, 602424–1—602424–6(2005).

**5. **M. Somayaji and M. P. Christensen, “Frequency analysis of the wavefront-coding odd-symmetric quadratic phase mask,” Appl. Opt. **46**(2), 216–226 (2007). [CrossRef] [PubMed]

**6. **K. H. Brenner, A. W. Lohmann, and J. Ojeda-Castaneda, “The ambiguity function as a polar display of the OTF,” Opt. Commun. **44**(5), 323–326 (1983). [CrossRef]

**7. **S. Bradburn, W. T. Cathey, and E. R. Dowski, “Realizations of focus invariance in optical-digital systems with wave-front coding,” Appl. Opt. **36**(35), 9157–9166 (1997). [CrossRef]

**8. **J. van der Gracht, J. G. Nagy, V. Pauca, and R. J. Plemmons, “Iterative restoration of wavefront coded imagery for focus invariance,” in Integrated Computational Imaging Systems, OSA Technical Digest Series (Optical Society of America, 2001), paper ITuA1.

**9. **F. Yan, Z. Li-gong, and Z. Xue-jun, “Image Restoration of an Off-axis Three mirror Anastigmatic Optical System with Wavefront Coding Technology,” Opt. Eng. **47**(1), 0170081–0170088 (2008).

**10. **Z. Ting-yu, W.-Z. Zhang, Z. Ye, and F.-H. Yu, “Design of wavefront coding system based on evaluation function of fisher information,” Acta Opt. Sin. **27**(6), 1096–1101 (2007) (in Chinese).

**11. **V. Sudhakar Prasad, P. Panca, R. J. Plemmons, T. C. Torgersen, and J. van der Gracht, “Pupil-phase optimization for extended-focus, aberration-corrected imaging systems,” Proc. SPIE **5559**, 335–345 (2004). [CrossRef]

**12. **S. Mezoouari and A. R. Harvey, “Wavefront coding for aberration compensation in thermal imaging systems,” Proc. SPIE **4442**, 34–42 (2001). [CrossRef]

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

**14. **S. Mezouari, G. Muyo, and A. R. Harvey, “Circular symmetric phase filters for control of primary third-order aberrations:coma and astigmatism,” J. Opt.Soc.Am.A **23**, 1058–1062 (2006).

**15. **L. Guang-zhi, X. Zhang, Z. Jian-pin, Y. Hao-ming, H. Feng-yun, and X. Zhang, “Novel optimization method for wavefront coding system,” Opt.Precision Eng. **16**(7), 1171–1176 (2008) (in Chinese).

**16. **F. Yan, Z. Li-gong, and Z. Xue-jun, “A Design of Off-axis Three Mirror Anastigmatic Optical System with Wavefront Coding Technology,” Opt. Eng. **47**(6), 063001 (2008). [CrossRef]

**17. **http://www.mathworks.com/matlabcentral/fileexchange/7507.

**18. **D. C. van Leijenhorst, C. B. Lucasius, and J. M. Thijssen, “Optical design with the aid of a genetic algorithm,” Biosystems **37**(3), 177–187 (1996). [CrossRef] [PubMed]

**19. **P. Siarry, A. Petrowski, and M. Bessaou, “A multipopulation genetic algorithm aimed at multimodal optimization,” Adv. Eng. Software **33**(4), 207–213 (2002). [CrossRef]