Design and fabrication of random phase diffusers for extending the depth of focus

: We present a method for designing non-absorbing optical diffusers that, when illuminated by a converging beam, produce a speciﬁed intensity distribution along the optical axis. To evaluate the performance of the diffusers in imaging systems we calculate the three-dimensional distribution of the mean intensity in the neighborhood of focus. We ﬁnd that the diffusers can be used as depth-of-focus extenders. We also propose and implement a method of fabricating the designed diffusers on photoresist-coated plates and present some experimental results obtained with the fabricated diffusers.


Introduction
In many applications, it is desirable to increase the focal depth of an imaging system without sacrificing its resolution and light-gathering capabilities. This has been a subject of longstanding interest in optics.
It is well-known that the use of annular apertures slows the decay of the axial intensity as one moves away from focus and, thus, that such masks can be used to extend the depth of focus [1]; the core of the Airy pattern also becomes narrower, transferring light to the outer rings of the point spread function. This fact, together with the reduction in the light throughput of the system, limit the usefulness of this approach. Consequently, there have been numerous investigations on the design of phase-only masks for improving the depth of focus (for some recent references see, e.g., [2,3,4,5,6]). Approaches as varied as spatial filtering [7], digital post-processing [8], and polarization coding [9] have been reported.
A common criterion employed in the evaluation of optical elements for improving the focal depth is the rate of decay of the axial intensity of the Point Spread Function (PSF) as one moves away from the focus: the slower the decay the better. Not surprisingly, the design of optical elements that produce a prescribed axial intensity distribution has been the subject of some studies. We can mention reports on computer generated holograms designed by iterative techniques [10,11], and more recent mathematical studies based on the superposition of Bessel beams [12,13].
Although many of the strategies to extend the depth of focus base their operation on inter-ference and diffraction, refractive methods have also been explored. Perhaps the most notable example of an optical element for producing an extended axial intensity distribution is the axicon [14,15]. In particular, generalized axicons [16] can produce a fairly constant axial intensity within some predefined limits. However, the usefulness of these optical elements in imaging systems is compromised because, on the one hand, the transverse distribution is not constant in the region of uniform axial intensity [17] and, on the other, axicons tend to suffer from off-axis aberrations.
In this paper, we consider the design of circularly-symmetric random phase diffusers that, when illuminated by a converging beam, produce a specified distribution of intensity along the optical axis. The designed diffusers have the following characteristics: i) they are random, ii) their operation is based on refraction and, iii) they have good lateral confinement properties. The performance of these diffusers is relatively insensitive to wavelength changes and, thus, can be used in applications that require broad-band illumination, such as in Optical Coherence Tomography (OCT) and white light profilometry.
The approach is motivated by previous work on the design of randomly rough surfaces with prescribed angular scattering properties [18,19,20,21]. Basically, the diffusers proposed here consist of arrays of concentric elements with randomly varying power that focus light with equal probability within the design region. The boundaries of these annular elements are similar in their definition to those of the Fresnel zones. Although the design is based on geometrical optics, calculations based on scalar diffraction theory show that the designed surfaces have the expected properties and can extend the focal depth substantially. We also describe in some detail a method for fabricating these rotationally-symmetric random diffusers on photoresist, and present experimental results with a model imaging system.

Design of the diffusers
In designing a random surface that produces a specified distribution of intensity along the optical axis, we assume that we are dealing with a circularly symmetric aberration-free imaging system. For a point source object the system produces a converging spherical wave. After passing through a diffuser with rotationally-symmetric amplitude transmittance function t(r), the complex amplitude in the neighborhood of focus can be expressed as [22] (see Fig. 1) In this expression A 0 is a constant amplitude and k 0 = 2π/λ = (ω/c) is the wave number in vacuum where, as usual, ω is the frequency, λ is the wavelength, and c the speed of light in vacuum. The length a is the radius of the circular aperture, R is the distance from the pupil plane to the best focus, z 0 is the defocus distance, r 0 = (x 2 0 + y 2 0 ) (1/2) is the transverse radial coordinate, J 0 (z) is a Bessel function of the first kind and order zero, and t(r) is a complex function of the form where H(r) represents the surface profile function of the diffuser. We assume that H(r) is a single-valued function of r that is differentiable, and constitutes a random process, but not necessarily a stationary one. Assuming only small angles of incidence and scattering, we adopt the thin phase screen model and set v 3 = 2(ω/c) for a reflection geometry, and v 3 = (ω/c)Δn for a transmission one [23,24,25], where Δn is the difference between the refractive indices of the material from which the diffuser is made and that of the surrounding medium. For definiteness, we consider here a transmission geometry, but the results can be applied equally well to a reflection geometry. We also assume that the phase excursions introduced by the diffuser are much greater than 2π and, thus, that the coherent component of the transmitted field is negligible.
Let us consider then the complex amplitude along the optical axis, which is given by we rewrite Eq. (2) as where The Fourier transform-like relationship between the radial amplitude transmitance function and the axial complex amplitude given by Eq. (4) is a particular case of a more general result obtained by McCutchen [26]. The similarities with a usual Fraunhofer diffraction problem [see e.g. Eq. (11) in Ref. [19]] are evident. Thus, if one can design one-dimensional diffusers that produce a prescribed angular distribution of the mean intensity [18,19,20], it should also be possible to design diffusers that produce predefined axial distributions. The mean intensity along the optical axis is where the angle brackets denote an average over the ensemble of realizations of the surface profile function H(r).
Our aim is to obtain the surface profile function h(t), and hence H(r), for which the righthand side of Eq. (6) produces a mean intensity I(z 0 , 0) with a prescribed dependence on z 0 . As it stands, Eq. (6) is too complicated for us to invert it to obtain h(t) in terms of I(z 0 , 0) . We therefore pass to the geometrical optics limit of the expression on the right-hand side of Eq.
.., and retaining only terms through the leading nonzero order in (t − t ). The consequences of this approximation will be discussed in the next section. In this way we obtain where h (t) is the derivative of h(t).
We now propose a particular form of h(t) that, as we will see, simplifies the design of the diffusers. We first introduce a characteristic length b through the definition a 2 = Nb 2 , where N is a large integer. We then represent the function h(t) in the form The coefficients {α n } in Eq. (8) are assumed to be independent identically distributed random deviates. The probability density function (pdf) of α n , is therefore independent of n. Its definition indicates that f (γ)dγ is the probability that α n lies in the interval (γ, γ + dγ) in the limit as dγ → 0. In order that the surface be continuous at, say, t = (n + 1)b 2 , the condition must be satisfied. From this recurrence relation the {β n } can be determined from a knowledge of the {α n }, provided that an initial value, e.g. β 0 , is specified. It is convenient to choose β 0 = 0, and we will do so in what follows. The solution of Eq. (10) is then The reasons for our choice of this form for h(t) will become clear below.
With the representation of h(t) given by Eq. (8) the double integral in Eq. (7) becomes where we have used the independence of the {α n } and the fact that they are identically distributed in the last step. Since the integrands in this expression now are independent of n, the sum on n can be carried out immediately, with the result that the double integral becomes This simple result is a consequence of our assumption of a linear dependence of h(t) on t in the interval (nb 2 , (n + 1)b 2 ), so that its derivative h (t) is a constant, h (t) = α n , in this interval. Equations (7) and (12) together give us where sinc x = sin x/x. Note that the condition Nb >> λ implies that a 2 v 3 /(2b) >> 1, and that in the limit as A → ∞, where δ (x) is the Dirac delta function. With the aid of this result Eq. (13) becomes Thus, the mean intensity along the optical axis is given in terms of the pdf of α n . On inverting Eq. (15) we find that From the result given by Eq. (16) a long sequence of {α n } is generated by the rejection method [27], and the surface profile function is constructed on the basis of Eqs. (8), (11), and the fact that t = r 2 . Thus, the surface profile function H(r) can be written as In what follows we consider the design of a diffuser that produces a uniform axial intensity in the interval −z m < z 0 < z m , and zero axial intensity for |z 0 | > z m , and evaluate the possibility of employing it to extend the depth of focus of imaging systems. The mean axial intensity we seek the diffuser to produce is therefore where I 0 is a constant, and θ (z) is the Heaviside unit step function. We obtain from Eqs. (16) and (18) the result The constant I 0 is obtained from the normalization of f (γ), with the result that It therefore follows that and

Three-dimensional distribution in the neighborhood of focus
If we intend to use the diffuser in an imaging system, the axial intensity distribution is not the only function of interest: it is also important to know the intensity distribution in the radial direction away from the optical axis. Consequently, in this section, we calculate the mean intensity in the focal region. We write the expression for the field in the form With the form of H(r) given by Eq. (17), this expression can be written in the form where The field ψ n (z 0 , r 0 ; α n ) represents the diffraction pattern of an annular pupil function with defocus z 0 + v 3 α n /(bκ), where α n is a random quantity. Diffraction integrals like the one represented by Eq. (25) have been well-studied in the past [1,22]. Here, we evaluate them using the Nijboer expansion [22,28].
The mean intensity obtained from Eq. (24) is Recall that we are assuming that the random numbers {α n } are statistically independent. If we further assume that the coherent component is negligible (i.e. that ψ n (z 0 , r 0 ; α n ) = 0), the expression for the mean intensity simplifies to The results of calculations not presented here show that the relative error made in using Eq. (27) instead of (26) is of the order of one part in 10 4 for all the parameters and values of z 0 and r 0 assumed in this work. Consequently, in what follows we will use the simpler expression (27) in calculating the mean intensity I(z 0 , r 0 ) . Let us consider first the mean axial intensity obtained from this expression. The axial amplitude produced by the n-th ring is given by After integration we obtain which yields for the mean axial intensity where we have used the relation a 2 = Nb 2 and the fact that the width of the diffraction structure obtained in Eq. (29) does not depend on the order of the ring. We note that the first term in the argument of the sinc function is a deterministic defocus term, while the second one contains the random numbers α n . Apart from this random term, the result agrees with the well-known expression for the axial intensity produced by an annular aperture [1,29]. So, apart from the random defocus, all the annular subapertures produce exactly the same axial intensity distribution. It is also interesting to compare Eq. (30) with the expression for the mean intensity obtained in Sect. 2. From Eq. (13), we can write Forgetting for a moment the random defocus term, we immediately recognize in this expression the axial response of an aberration-free optical system [22]. We observe that the parameter N that appears in Eq. (30) is absent in Eq. (31). This is due to the fact that with the approximation employed in its derivation (geometrical optics limit) we have neglected diffraction effects by the small scale of the diffuser, retaining only the diffraction effects of the whole pupil. Consequently, the diffraction pattern obtained is related to the complete pupil, rather than the individual annuli, as it should be. The use of relation (14), removes the remaining diffraction effects, producing a final result that is consistent with geometrical optics [Eq. (15)]. So, the results obtained in Sect. 2 are useful in the design of the diffuser but cannot be used for diffraction calculations. The result given by Eq. (27), expresses the mean intensity as a sum of the intensity diffraction patterns |ψ n (z 0 , 0; α n )| 2 produced by different zones of the diffuser. To help in the visualization of this result, in Fig. 2 we show two of these contributions for the particular case α n = 0. The optical system has a = 2 cm, R = 15 cm, and λ = 0.633 μm. The circular aperture was divided into N = 100 annular apertures, resulting in a value of b = a/ √ N = 0.2 cm. The diffraction pattern of Fig. 2(a) corresponds to the zone with n = 0, and is just the diffraction pattern produced by a circular aperture of radius b. Figure 2(b) corresponds to the annular pupil with n = 25. The figures represent meridional sections of the rotationally symmetric three-dimensional diffraction patterns. These figures confirm our earlier remark that the distribution of the axial intensity is independent of the order of the ring. In this case, the first axial zeroes occur around z 0 = ±0.71 cm. On the other hand, the transverse distribution decreases in size with the order of the ring. We consider now, as an example, the design of a diffuser that produces a uniform distribution of intensity in the range −2cm ≤ z 0 ≤ 2 cm of the optical axis. The pdf of α n for this case is given by Eq. (22), so that from Eq. (27) we have As before, we assume that a = 2 cm, R = 15 cm, λ = 0.633 μm, and N = 100. The parameter Δn = 0.6.
The result for the rotationally-symmetric mean intensity I(z 0 , r 0 ) is illustrated in Fig. 3. From this meridional intensity map we see that the intensity is fairly constant in the design region of the optical axis, and that it decreases rapidly outside it and for off-axis points. These same data were used to generate Fig. 4, where we show the mean intensity distribution along the optical axis (a) and the mean intensity in the transverse direction for the plane z 0 = 0 (b). From Eq. (30), we observe that the mean axial intensity can be expressed as the convolution of a sinc 2 function with the desired distribution. Consequently, the axial intensity distribution is not perfectly rectangular, but smoothed by diffraction effects whose extent depend on the parameter 2λ (R/b) 2 . For reference, we mention that the first zeroes of the axial response of the system without diffuser occur at z 0 = ±71 μm. Also, in Fig. 4(b), we show the Airy pattern that corresponds to the transverse response of the system with a clear pupil function. We see that the width of the central core is similar in the two cases and that the main difference is the relatively slow decay of the tails in the response of the system with the diffuser. Summarizing, the mean PSF consists of a bright central core with a broad halo, and its structure is fairly invariant within the design region. A predictable consequence of the halo is that the high frequency components will be transfered to the image with low contrast. Despite this, the sharpness of the central core indicates that the diffusers can be useful in imaging applications. It should be mentioned, however, that in normal circumstances, the estimation of the ensemble average is not a practical proposition and one has to work with a single realization of the diffuser. For a particular realization, different diffuser rings focus light on different axial points, with equal probability within the design region. This is illustrated in Fig. 5. For monochromatic illumination, the interference between all these randomly phased contributions produces speckle, which is manifested by random intensity variations in the axial intensity. At the same time, due to the rotational symmetry of the system, the transverse intensity pattern contains rings that change rapidly as one moves along z 0 . As an example, we show in Fig. 6 (a) and (b) the calculated transverse intensity images at z 0 = −1 and z 0 = 1 cm.  These fluctuations can be smoothed through the use of broadband illumination. In Figs. 6 (c) and (d) we show the PSFs obtained at z 0 = −1 cm and z 0 = 1 cm by averaging over 400 different wavelengths in the visible region of the spectrum. One can see that the ring structure has practically disappeared and that the shape of the PSFs approaches that of the ensembleaveraged PSF. This resemblance improves with the number of independent patterns employed and with the number of rings N.

Experimental techniques and results
Some surfaces of the type considered here have been fabricated by exposing photoresist-coated plates to blue light (λ = 442 nm) from a He-Cd laser transmitted through a rotating ground glass (to reduce its coherence). A schematic diagram of the experimental setup employed in the fabrication is shown in Fig. 7. An incoherent image of a disk-shaped mask is formed on the rotating photoresist-coated plate by a well-corrected imaging system with magnification m. The plate is exposed during a time T e , during which it executes a large number of revolutions. As explained below, the arrangement is such that the total exposure of the plate is a scaled version of the profile function employed in the generation of the mask. Before proceeding further with the description of the experimental it is worth discussing the consequences of scaling the profiles H(r) generated according to Eqs. (17), (22), and (11) in both, the vertical and horizontal directions. That is, we are interested on the mean intensity produced by a diffuser with a surface profile function AH(mr), where A and m are dimensionless constants, with reference to the mean intensity obtained with the original function H(r). From Eq. (17) we see that the transformation is equivalent to chosing new random deviates α n = Am 2 α n . In consequence, the parameter that defines the new region of constant intensity is z m = Am 2 z m . In other words, by scaling the function H(r) in the vertical or horizontal direction one changes the length of the region of constant intensity. These transformation are almost inevitable in our fabrication scheme and are defined by the exposure of the plate and by the magnification of the optical system (see Fig. 7).
To produce a suitable mask, one first needs to generate a realization of a profile H(r). An example is shown in Fig. 8(a). Then, we define the function Δ H (r) = KH(r) that, with an appropriate choice of the units of the constant K, can be interpreted as an angle. For a given radius, the angles θ that fall in the transparent section of the mask are defined by the condition Δ H (r) > θ > Δ 0 , where Δ 0 is a constant smaller than the minimum value of Δ H (r) [see Fig.   8(b)].
An incoherent image of the mask is formed on the surface of the rotating photoresist plate, producing an exposure that is circularly symmetric with a radial dependence of the form E(r) = I e T e [Δ H (mr) − Δ 0 ]/2π, where I e is a constant related to the intensity of the illumination. The exposure of the plate can then be put in the final form E(r) = E 0 + AH(mr), where E 0 and A are constants that one can adjust by varying the intensity of the light reaching the plate, the aperture of the mask and the exposure time. Assuming that the relation between exposure and resulting height on the plate is linear, the developed surface will have the desired properties of being proportional to AH(mr). A number of samples have been fabricated with the procedure described above. We present some experimental results that illustrate the effects of using these diffusers on the formation of extended images. Photographs of a wheel-like object taken under white light illumination are shown in Fig. 9. The photographs on the top row correspond to the best-focus position. The image on the left (a) was formed by an aberration-free optical system with a clear pupil function. The one on the right (b) was formed by the same system, but with the inclusion of a diffuser in the pupil plane. One can see that, due to the halo of the mean PSF, the edges of the image show a reduction in contrast. Images corresponding to an out-of-focus situation are shown on the bottom row of Fig. 9. The image on the left (c) is noticeably blurred and its quality has been degraded with respect to the one on top. On the other hand, the image on the right (d), obtained with the system that includes the diffuser, has changed only slightly. The resolution of the images can be assessed by inspection of the fine details in the circled regions. These results show that, as expected, the optical system with the diffuser has a greater depth of focus.

Summary and concluding remarks
In this paper, we have studied the problem of designing and fabricating optical diffusers that, when illuminated by a converging beam, produce a uniform intensity along the optical axis within a specified range.
The technique is also suited for the production of such surfaces in photoresist, and we have implemented a procedure for their optical fabrication. The results indicate that good approximations to the desired uniform axial intensity distribution can be obtained with surfaces fabricated with the proposed method. Many of the techniques employed for extending the depth of focus have a preferred wavelength of operation. In contrast, the performance of the optical elements studied here improves with polychromatic light. The method is verified through calculations based on scalar diffraction theory. The results show that the mean intensity produced by these diffusers is indeed fairly constant over the specified region and decays as the point of observation moves away from the optical axis. The diffusers are potentially useful as focal depth extenders in imaging systems that use broad band illumination.