Open Access
Add to BibSonomyAbstract
We propose a numerical analysis method for evaluating GRIN lenses using the Monte Carlo method. Actual measurements of the modulation transfer function (MTF) of a GRIN lens using this method closely match those made by conventional methods. Experimentally, the MTF is measured using a square wave chart, and is then calculated based on the distribution of output strength on the chart. In contrast, the general method using computers evaluates the MTF based on a spot diagram made by an incident point light source. However the results differ greatly from those from experiments. We therefore developed an evaluation method similar to the experimental system based on the Monte Carlo method and verified that it more closely matches the experimental results than the conventional method.
©2010 Optical Society of America
1. Introduction
Gradient index lenses (GRIN lenses) have excellent optical design characteristics compared with homogeneous lenses and are widely used in optical information processing and communication applications as new optical devices. Historically, Uchida, et al. [1] from Nippon Sheet Glass Co., Ltd. reported the development of a glass GRIN lens in 1970. In 1978, a lens array with several aligned GRIN lenses was developed and commercialized as the Selfoc® lens array. This has been used in various applications, such as for the optical information reading element in scanners, fax machines, and copy machines.
Cylindrical GRIN lenses are classified into two types depending on the direction of refractive index distribution: axial-GRIN lens and radial-GRIN lens. The axial-GRIN lens has a refractive index that changes one-dimensionally in the direction of cylinder length, and may be useful as an aplanar lens [2]. In contrast, in the radial-GRIN lens the refractive index is distributed two-dimensionally, gradually decreasing from the central axis to the circumference. GRIN lenses are characterized by bending rays of light toward higher refractive index within the medium and thus produce lens effects without curvature on either end (even when the ends are flat). In this paper, we studied the radial-GRIN lens (hereafter, “GRIN lens”). The cylindrical GRIN lens, which is used as the optical information reading element in scanners and copy machines, is also called a rod lens, and its refractive index distribution decreases in the shape of a parabolic curve from the central axis to the circumference along the radius [3,4]. As this lens is capable of forming inverted real images or erect magnified images, it requires only a short working distance, which makes it suitable for miniature, high-performance devices.
Whether a rod lens satisfies the specified performance is usually evaluated by measuring the modulation transfer function (MTF). The most widely used method for doing this experimentally is to use a square wave chart: incident light is applied to the rod lens through the square wave chart, and the MTF is calculated based on the maximum and minimum values of the output intensity distribution. For their part, the MTF of rod lenses is often evaluated using computers running commercial optical design software programs such as CODE V® and ZEMAX®. In this method, an incident point light source is applied to the rod lens instead of the square wave chart to calculate the MTF based on the spot diagram. However, these programs use a different MTF calculation method from the theoretically true evaluation method and so their results may differ greatly from experimental results, as will be shown later. Thus, a computer-based method for accurately evaluating the performance of rod lenses is needed.
We propose a rod lens evaluation method using a square wave chart based on the Monte Carlo method. Incident light applied through a diffuser panel is filtered through a square wave chart into a rod lens. The incident image is reconstructed by integrating the light ray passing through the rod lens onto the imaging plane to evaluate the MTF. We show that the method precisely reproduces the experimental results.
2. Principle of imaging in a gradient index lens and MTF evaluation method by experiment
As shown in Fig. 1 , a light ray that passes through a gradient index lens is bent toward the center with higher refractive index as a result of the refractive index difference formed within the lens, and the light is propagated as a sine wave.
Nishizawa [5] showed that the lens produces acceptable practical imaging performance if the refractive index in the radial direction is expressed by the following equation based on theoretical analyses by Fletcher, et al. [3] and Rawson, et al. [4]:
where, r is the distance from the central axis, n(r) the refractive index at distance r from the central axis, n 0 the refractive index on the central axis, and g the refractive index distribution constant (g value). The cycle of light P propagating through a GRIN lens with refractive index distribution expressed by Eq. (1) is expressed by [5]:A GRIN lens has the imaging characteristics shown in Fig. 2 , and forms an erect real image of the same size when it is cut to the length of 3/4 P [6]. An alignment of many such GRIN lenses of length forming an erect real image of the same size is called a rod lens array, and such arrays are used in scanners, fax machines, copy machines and so forth.
Figure 3 shows the measurement system of the MTF. The principle for measuring the MTF of a rod lens is as follows. The incident light enters the rod lens array through a square wave chart and forms a pattern behind the lens array. Here, a diffuser plate and a wavelength filter are also used to generate the light with a perfect diffuse surface and specific wavelengths. A formed pattern on the imaging plane is measured with a one-dimensional CCD sensor. The MTF value is calculated from an intensity distribution of the imaged pattern.
The MTF measurement method using a square wave chart is shown in Fig. 4 .The image for precisely-printed square wave grid patterns is received by a CCD sensor via the rod lens array and measured based on the quantity of light. An image closer to the original image is transmitted as the MTF value approaches 100%. Figure 5 shows how to set up the working distance in MTF measurement using a square wave chart.
As shown, l is the working distance, Z the lens length and TC the conjugation length. To measure the MTF, the system is arranged so that an image of the same size is formed by setting the distance from the rod lens to the square wave chart and the distance to the sensor both equal to l.
3. MTF evaluation method by computer simulation
In evaluating the MTF of a rod lens using a computer, first an array of GRIN lenses is created so that the phenomenon by which light from outside the central lens affects the MTF is taken into consideration; this enables imaging performance to be analyzed using actual conditions. In the simulation, multiple lenses are aligned as shown in Fig. 6 to evaluate the MTF.
Based on the conditions described above, the MTF was evaluated using the ray tracing method.
3.1 MTF calculation method in the conventional method
The MTF evaluation method used in commercial optical design software programs such as CODE V® and ZEMAX® is as follows. An incident point light source is applied to the rod lens as shown in Fig. 6 and the light ray distribution on the image plane is calculated using the ray tracing method. In this case, the point light source is assumed to be a perfectly diffusing surface. The MTF is calculated from the distribution of light rays on the imaging plane or the spot diagram.
When normalization of the Fourier-transformed distribution of point image intensity is defined as the optical transfer function (OTF) and I(x,y) as the point image intensity distribution on the imaging plane (PSF: point spread function), s as the spatial frequency in the sagittal direction and t as the spatial frequency in the tangential direction, we obtain:
Equation (3) generally expresses a complex number and its absolute value is called the MTF. Based on the spot diagram, the coordinate system for calculating OTFG(s,t), which is geometrical-optical OTF, is as shown in Fig. 7 .
Fig. 7 Coordinate system for calculating the geometrical-optical intensity distribution from the spot diagram.
Supposing the area of the pupil for the optical system is A and the energy passing through the unit area of the pupil is uniform at 1, then the total light quantity passing through the optical system is A. When the intensity distribution obtained from the spot diagram is Ip(x,y), from Eq. (3) we obtain:
Dividing the pupil plane into N rectangles of equal area, letting the coordinates at which the rays are incident on the imaging plane having passed through the vertex of the jth grid be (xj,yj), and supposing A = Ndudv, we obtain:
The MTF is obtained by taking the absolute value of Eq. (5).The results obtained by this method are a sine wave response function, whereas a square wave response function is measured using a square wave chart by experiment. To compare the experimental results with the simulation results, one of them must be converted to match the other; this conversion method is called Coltman's correction [7].
A square wave response function is converted into a sine wave response function by:
where, ω denotes the spatial frequency, the sine wave MTF, and the square wave MTF. If the square wave MTF for spatial frequency is measured for a wide range, the sine wave MTF can be calculated using Eq. (6).In this paper, the MTF values measured in the experimental system are all for the spatial frequency of 6 LP/mm, and thus the sine wave MTF obtained by simulation is converted into a square wave MTF by reverse operation using the following equation to compare the experimental and simulation results:
3.2 MTF calculation method using the Monte Carlo method
We propose the following method of evaluating the MTF by the Monte Carlo method. The Monte Carlo method is a numerical analysis using pseudorandom numbers, and has been applied in optical fields such as illumination optics and computer graphics [8,9].
The square wave charts used in experiments have several light and dark parts drawn in each millimeter depending on the spatial frequency, and the light passes through only the light parts. These parts are a type of surface light source, and are located at the positions shown in Fig. 6. In this case, light has passed through the diffuser plate in front of the square wave chart and thus the incident rays coming into the rod lens from the chart are considered to advance in random directions as shown in Fig. 8 . Therefore, the direction in which each ray advances is determined using pseudorandom numbers, which we generate using a Mersenne Twister [10].
Rays are generated from an incoherent light source from the white parts of the chart in Fig. 8, and the geometrical-optical irradiance distribution on the illuminated plane is calculated. The flow chart for the simulation procedure is shown in Fig. 9 .
Pseudorandom numbers are generated, and the random values are used to determine the positional coordinates and the direction of emission of the point light source. As the coordinates and generalized momentum are given as initial values, the ray distribution on the square wave chart can be obtained by implementing ray tracing for the point light source and repeating the procedures for the desired number of incident rays. However, a very large number of incident rays is needed as the Monte Carlo method is stochastic.
While the ray distribution on the imaging plane can be obtained by implementing the Monte Carlo ray tracing method, the output image should be reconstructed based on the ray distribution in order to calculate the MTF. While methods using window functions and averaging filters can be used for reconstructing the image, we need to obtain the precise light quantity obtained from the image, instead of obtaining the image itself. The light quantity ratio is changed when folding by a window function or an averaging filter is applied depending on the parameter settings, and so the resulting MTF value may no longer be absolute. Therefore, we reconstructed the image by dividing a square of side 1 mm into grids and integrating the rays arriving at a specific pixel. A conceptual drawing is shown in Fig. 10 .
4. Simulation and experimental results
The results of the simulation and experiment were as follows. First, we verified whether simulation using the Monte Carlo method as described in the previous section is effective or not. The distribution of the light of single GRIN lens when scanned in the y direction is expressed analytically as the following Eq. (8) [11]
where, k is the product of the transmittance of the lens and luminance of the light source, θ the maximum semifield of the lens, l the image surface distance, D the diameter of the lens. Simulation results showing a comparison with the analytic formula are shown in Fig. 11 .Where, Z = 4.00 mm, l = 2.80 mm, D = 200 μm, θ = 0.134 rad, and the effective number of light rays is 5.0 × 107. We defined a ray entering the lens as an effective ray. These results show the proposed method reproduces analytic formula nicely, and we could confirm this method is effective. Next, we evaluated the minimum necessary number of rays for this method. We studied the number of rays where the root mean square error (RMSE) between the simulation result and the analytical solution became stable by increasing the number of rays evaluated for a rod lens of ideal refractive index distribution and plotting the changes in RMSE. The resulting RMSE values for the number of rays evaluated are shown in Fig. 12 .
Figure 12 shows that the RMSE value stabilizes when the number of rays evaluated reaches about 50 million. As Fig. 12 shows, there is a trade off relation between the accuracy and computational costs. Since image reconstruction on the imaging plane is required in this method, about 104 times the number of rays is needed compared with the conventional method. Fig. 13 shows the resulting image produced by inputting an image into a single GRIN lens (φ265 µm) with ideal refractive index distribution (g = 836.6 1/m, n 0 = 1.544) and reconstructing it by integration. Here, the number of pixels was set to 256 × 256 and the length of one side to 1 mm. The length of the lens was 4.4 mm and its working distance was 2.8 mm.
As the results clearly show, the image is reconstructed with good precision using this method.
4.1 Comparison between experimental results and simulation results
The refractive index distributions used in comparing the MTF values measured by experiment and those by simulation analysis are shown in Fig. 14 . These figures show the relation between radial distance from the center of the GRIN lens and the refractive index. The distribution has uniformity in the azimuthal direction.
The refractive index distribution data is measured using a transmission interference microscope with dual beam, and GRIN lenses are manufactured by the continuous lamination process [12,13].
The results of evaluating the MTF values by experiment and simulation are shown in Table 1 . The experimental results were obtained using a square wave chart of 6 LP/mm, and the evaluation results by the conventional simulation were obtained by implementing Coltman's correction considering up to the second term in Eq. (7) based on the MTF values of spatial frequencies 6 LP/mm and 18 LP/mm. In the Monte Carlo method, the pseudorandom number distribution was set as x (mm) is [−0.3,0.3] and y (mm) is [−0.5,0.5], px is [−0.1,0.1] and py is [−0.1,0.1], and the number of incident rays was set to 500 million. The width of the calculation was Δz = 5 µm, Δx = Δy = Δz/10000, respectively. The chart position was set up so that the dark part at the center of 6 LP/mm would come to the center of the array of 7 lenses.
Table 1. MTF evaluation results by experiment and simulation
In an ideal condition (g = 836.6 1/m, n 0 = 1.544), there is no difference between the conventional method and proposed method, and the MTF value reaches 100% in both methods. Differences appear in the evaluation of practical rod lenses. The results of the Monte Carlo method clearly match the experimental results more closely than those of the conventional method. This is because that Eq. (3) is based on linearity of the optical system. The positional dependence of the PSF decreases the evaluation accuracy of the MTF with Eq. (5). Figure 15 shows spot diagrams when the position of the light source is changed.
Fig. 15 Spot diagrams of rod lens arrays with changing the incident position. (a) Incident position is x = 0 mm (ideal condition). (b) Incident position is x = 0.04 mm (ideal condition). (c) Incident position is x = 0 mm (data no. 1). (d) Incident position is x = 0.4 mm (data no. 1).
In the ideal rod lens, spot diagrams are localized to a small region, and these distributions are nearly independent of the position of light source. On the other hand, spot diagrams of the practical lens are distributed over a wide range, and vary with location. This means that the PSF depends on the location. The output images and their intensity distributions calculated using the refractive index distribution for data no. 1, 2 and 4 are shown in Fig. 16 . The light has diffused widely in a rod lens with a low MTF. These results match the experimental results.
Fig. 16 Reconstructed images and intensity distribution obtained by the Monte Carlo method. (a) Data no. 1 (Z = 5.01 mm, l = 1.40 mm). (b) Data no. 2 (Z = 4.71 mm, l = 1.20 mm). (c) Data no. 4 (Z = 4.50 mm, l = 2.20 mm).
Figure 17 shows similar simulation results with changing spatial frequencies. In this simulation, spatial frequencies of the square wave chart are 4 and 8 LP/mm, and the refractive index distribution for data no. 2 is used. MTF values of these results are 67% at 4 LP/mm and 32% at 8 LP/mm. As results show, the MTF decreases with increasing the spatial frequency, and the reconstructed image and the intensity distribution of 8 LP/mm chart are distorted.
Fig. 17 Reconstructed images and intensity distribution with changing spatial frequencies (Z = 4.71 mm, l = 1.20 mm). (a) 4 LP/mm. (b) 8 LP/mm.
5. Conclusion
We have developed a numerical analysis method for evaluating rod lenses based on the Monte Carlo method. We compared the actual MTF measurement results and the evaluation results using this method, and confirmed that our method reproduces the experimental results with good precision. We also showed that this method is superior to the conventional method under nearly all of the conditions evaluated in this study.
We believe this method will be useful for designing optical systems using rod lenses, and expect widespread application in the field of optical design.
References and links
1. T. Uchida, M. Furukawa, I. Kitano, K. Koizumi, and H. Matsumura, “Optical characteristics of a light-focusing fiber guide and its applications,” IEEE J. Quantum Electron. 6(10), 606–612 (1970). [CrossRef]
2. Y. Koike, Y. Takezawa, and Y. Ohtsuka, “New interfacial-gel copolymerization technique for steric GRIN polymer optical waveguides and lens arrays,” Appl. Opt. 27(3), 486–491 (1988). [CrossRef] [PubMed]
3. A. Fletcher, T. Murphy, and A. Young, “Solutions of Two Optical Problems,” Proc. R. Soc. Lond. A Math. Phys. Sci. 223(1153), 216–225 (1954). [CrossRef]
4. E. G. Rawson, D. R. Herriott, and J. McKenna, “Analysis of refractive index distributions in cylindrical graded-index glass rods (GRIN Rods) used as image relays,” Appl. Opt. 9(3), 753–759 (1970). [CrossRef] [PubMed]
5. K. Nishizawa, “Chromatic aberration of the Selfoc lens as an imaging system,” Appl. Opt. 19(7), 1052–1055 (1980). [CrossRef] [PubMed]
6. N. Yamamoto and K. Iga, “Evaluation of gradient-index rod lenses by imaging,” Appl. Opt. 19(7), 1101–1104 (1980). [CrossRef] [PubMed]
7. J. W. Coltman, “The Specification of Imaging Properties by Response to a Sine Wave Input,” J. Opt. Soc. Am. 44(6), 468–469 (1954). [CrossRef]
8. H. W. Jensen, “Global Illumination using Photon Maps,” in Proceedings of the 7th Eurographics Workshop on Rendering, P. Schröder, ed. (Porto, Portugal, 1996), pp. 21–30.
9. H. W. Jensen, S. R. Marschner, M. Levoy, and P. Hanrahan, “A Practical Model for Subsurface Light Transport,” in Proceedings of the 28th annual conference on Computer graphics and interactive techniques, L. Pocock, ed. (Los Angeles, Calif., 2001), pp. 511–518.
10. M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul. 8(1), 3–30 (1998). [CrossRef]
11. M. Kawazu and Y. Ogura, “Application of gradient-index fiber arrays to copying machines,” Appl. Opt. 19(7), 1105–1112 (1980). [CrossRef] [PubMed]
12. M. Oda, S. Suga, H. Yoshii, and T. Furuta, “Multilayer coating by drawing a thin plastic fiber through a polymer solution,” Asia Pac. J. Chem. 3(1), 63–69 (2008). [CrossRef]
13. M. Oda, S. Suga, H. Yoshii, and T. Furuta, “Multi-layer coating by continuous withdrawal of a thin plastic fiber through polymer solution,” Kagaku Kogaku Ronbunshu 34(1), 187–193 (2008). [CrossRef]
