Adaptive, spatially-varying aberration correction for real-time holographic projectors.

A method of generating an aberration- and distortion-free wide-angle holographically projected image in real time is presented. The target projector is first calibrated using an automated adaptive-optical mechanism. The calibration parameters are then fed into the hologram generation program, which applies a novel piece-wise aberration correction algorithm. The method is found to offer hologram generation times up to three orders of magnitude faster than the standard method. A projection of an aberration- and distortion-free image with a field of view of 90x45 degrees is demonstrated. The implementation on a mid-range GPU achieves high resolution at a frame rate up to 12fps. The presented methods are automated and can be performed on any holographic projector.


Introduction
2D holographic projection presents numerous advantages over traditional projection, namely: the potential for miniaturization, high efficiency, shock resistance due to the absence of moving parts [1][2][3][4][5][6][7] and aberration correction [7][8][9], the subject of this paper. To achieve optimal image reproduction using conventional methods, high-quality optical elements need to be used. Lenses of such quality have complex optical designs in order to minimize optical aberrations, but this makes them bulky and expensive. Holographic projectors rely on the use of Spatial Light Modulators (SLMs), but SLMs that exhibit the required flatness are difficult to manufacture and are therefore also expensive [1,7]. Another source of aberrations is the misalignment of optical components, which is always introduced to some degree in the manufacturing process.
It is therefore substantially more cost-effective to correct aberrations in situ using a holographic aberration correction algorithm rather than optimizing all of the factors mentioned above. Aberrations are disturbances in the phase of the wave [10]. In contrast to the ideal, spherical wavefront shown in Fig. 1(a), which focuses to a diffraction-limited spot on a screen as in Fig. 1(b), an aberrated wavefront like that in Fig. 1(c) will depart from the ideal shape and, after focusing, will lead to a smeared-out point with low contrast, like that of Fig. 1(d). Holography allows a high degree of control over the wavefront, and hence enables aberration correction [1,[7][8][9]. When the aberration profiles of the lens and the SLM are known, the virtual counter-aberration phase mask can be added to the wavefront. When such a wavefront passes through a lens, the virtual aberration cancels out the actual aberrations of the system, leading to an aberration-free, diffraction limited image, as illustrated in Fig. 2. To demonstrate the principle of aberration correction, in this work we deliberately use a low-quality sapphire ball lens [7][8][9]. Because of its wide field-of-view, such a lens produces a substantial amount of aberrations, which varies spatially across the replay field (explained in detail in section 2.2).
The method most commonly used to compensate the spatially-varying aberrations of a ball lens is the Pixel to Wrapped Phase Summation (PWPS). Instead of employing a Fourier Transform operation, the PWPS method splits the image into single-pixel contributions (gratings) on which the correction is then superimposed [8,9]. All existing implementations require an initial aberration characterization of every point in the replay field, which can only realistically be achieved by simulating the system. In this case the corrections are of limited quality due to discrepancies between simulations and the real physical system. In this work, an experimental approach to spatially varying aberration correction is therefore presented. Traditional PWPS is also very slow as a different correction hologram must be applied to every point. A MATLAB implementation can take up to 15 hours to generate a single corrected hologram. A parallelized GPU implementation of the same algorithm can speed up the process significantly reducing running time to a few minutes. Nonetheless, this method is still too slow to provide real-time hologram generation.
Currently, the only algorithm that can achieve real-time hologram generation is the One-Step Phase Retrieval (OSPR) Algorithm [1][2][3][4]11], but it is not able to achieve spatiallyvarying aberration correction such as aberrations of the ball lens discussed here. Another benefit of using an OSPR-type approach is the fact that there exist a number of variations of OSPR, such as an Adaptive OSPR or Adaptive OSPR with Liu-Taghizadeh optimization, that can be used to further improve the image quality [1,11].
In this paper we present a novel approach that combines these two algorithms by dividing a replay field into a set of regions and assuming that within each region, the aberrations are approximately spatially invariant. The regions chosen approximate concentric rings, a selection that reflects the spherical symmetry of the system. This allows an approximate version of PWPS to be applied that relies on experimentally characterized aberrations. Further, because only a small number of different corrections need to be applied, rather than one for every single point in the replay field, the algorithm can be applied in real-time. For the projector used in this research, 6 regions proved sufficient to correct a replay field of 90 by 45 degrees. It is possible to use an even wider field of view but further from the center, the size of the point spread function grows significantly because of the diffraction limit of the system. It is shown that this approach is substantially faster than existing algorithms and allows real-time correction for video-rate projection.

Aberrations in holography
The replay field of a hologram in the presence of aberrations can be represented as [1]: Transform operation. In this work, it is assumed that the illumination is uniform. In order to correct the distorted replay field, a corrective phase profile must be determined that compensates ( , ) x y φ . This corrective profile can be approximated as a linear combination of Zernike Polynomials.
Zernike Polynomials are frequently used in optics [7][8][9][10], because they correspond to common aberrations, such as defocus, astigmatism, coma, spherical aberration and higherorder aberrations [10]. They also form a complete, orthogonal set, which means that any smooth function defined on the unit circle can be approximated up to an arbitrary precision given a large number of terms. We will use this property and rewrite the corrective phase profile as a sum of Zernike Polynomials: where i a is the i -th Zernike Coefficient and i Z is the i -th Zernike Polynomial (using a single-numbering scheme used in ZEMAX as Zernike Fringe Phase [12] and described by Wyant and Creath [10]) and N is the index of the last term used in the expansion. We have ignored the first three terms: piston and tip and tilt, hence, the summation starts at 4 a (defocus coefficient).
This corrective phase profile can be thought of as being like a virtual Fresnel lens typically used to correct for defocus [1]. However, in this case instead of being a quadratic surface, the phase profile has a generalized high-order polynomial surface. If this virtual aberration component is chosen so as to exactly cancel out the system's aberrations, diffraction-limited performance can be achieved.

Spatially-varying aberrations
In the case of spatially-invariant aberrations, i.e. dependent only on the hologram coordinates but not on the image coordinates, one phase mask ( , ) x y φ is sufficient to correct the entire replay field: However, for spatially-varying aberrations in their most generic form, the aberration phase profile depends not only on the hologram coordinates, but also on the spatial coordinates: ( , , , ) x y u v φ [10,13]. The problem of spatially varying aberrations can alternatively be viewed as considering the spatially sampled optical system to have a transmission matrix that contains non-zero off-diagonal elements. As such, it represents a subset of a broader class of problems that seeks to shape wavefronts to correct aberrations of generalized scattering media [14]. A system with spatially invariant aberrations can be considered to have a perfect memory effect, while this effect is reduced in the case of spatially varying aberrations [15]. Now, assume that we perform aberration correction at some spatial coordinate 0 0 ( , ) u v and then attempt to apply the correction to the entire hologram: It is noted that this phase correction will cancel out exactly only at the point 0 0 ( , ) u v , but some degree of aberration will remain as the spatial coordinates deviate from this, due to a memory effect. In this work, we present a novel variant of One-Step Phase Retrieval (OSPR) suited for correcting such wavefront errors termed Piecewise Corrected OSPR and compare it with the previous Pixel-To-Wrapped Phase Summation method developed in [8,9].

Aberration correction algorithms
In the following discussion, it is assumed that the aberration coefficients for the different regions of the reply field have been previously determined by an algorithm such as that presented in [7] or [8,9]. Consequently, we are only concerned with the generation of appropriate holograms that correct these aberrations in multiple spatial regions to form high quality images.

Pixel-to-Wrapped Phase Summation (PWPS)
The PWPS algorithm is a standard aberration correction algorithm for spatially varying aberrations [8,9]. In order to introduce the development of this algorithm, it is noted that a hologram of a single point in the replay field is a continuous phase surface -an unwrapped grating: where ( , ) H x y is the complex hologram, A is the amplitude of the pixel, and uv pixel Φ is the continuous phase surface corresponding to a pixel at a position ( , ) where ( , ) u v ϑ is the phase of the pixel at a position ( , ) u v and max max ( , ) u v is the size of the replay field in pixels. In the case of this research, it is 1280x640. However, in general such a hologram will not produce a clean point in the replay field due to optical aberrations in the system. As discussed in section 2.2, in order to correct an arbitrary optical aberration for a point, a corrective phase mask has to be applied: is a phase mask defined in Eq. (2). However, this model does not place any constraints on the corrective phase profile. Therefore, we can instead use the spatially-varying representation of the aberration phase mask: ( , , , ) x y u v φ . The complex hologram that corrects the entire replay field containing spatially-varying aberrations can then be constructed: If the phase mask was not present, the preceding equation would be equivalent to an ordinary Discrete Fourier Transform. With the introduced phase component, however, it can be thought of as a summation of virtual aberrating Fresnel lenses. Each lens is designed to exactly cancel out the aberrations of the optical system at the particular replay field position.
In the final step, the hologram is quantized to binary phase due to the ferroelectric binaryphase SLM used in this research: It would be memory-inefficient to store all the values of ( , , , ) x y u v φ , therefore we represent it as: where ( , ) i zmap u v are arrays that store the Zernike Coefficients for the i -th Zernike Polynomial i Z , and are termed Zernike Maps [8,9]. This procedure reduces the memory requirement significantly (from 11 8.5 10 × entries to 7 3 10 × in our case). This procedure again produces a complex-valued hologram, which is quantized to become a binary hologram. In order to reduce the quantization noise, many binary holograms with different random phases need to be produced. Therefore, the same set of operations with different input ( , ) u v ϑ has to be executed. Constructing such a set of holograms significantly increases computation time.
The key shortcoming of all implementations of PWPS to date is the requirement to know the spatially varying aberration Zernike Maps for every point in the replay field. This can only feasibly be achieved by modelling the optical system as a point-by-point characterization of the replay field. Furthermore, applying this correction to each hologram is very slow. A MATLAB implementation of this algorithm can take up to 15 hours to generate a single corrected hologram, excluding the time it takes to obtain correction parameters initially. Using General-Purpose Graphical Processing Unit (GPGPU) programming, this time can be reduced to less than 4 minutes, but that is still orders of magnitude too long for real-time video projection.

One-Step Phase Retrieval (OSPR) with field-independent correction
Over a small area of the replay field, the aberrations can be assumed to be spatially invariant and so a single set of Zernike phase polynomials can be used to describe a correction valid over the whole region. This can be determined using a point-spread function in an optimization routine, as described in [7].
However, even after this correction is applied, there is still noise in the produced replay field caused by the phase quantization present when using a binary-phase hologram. One-step phase retrieval (OSPR) is a technique to overcome this, commonly used in projectors. OSPR evolved from the realization that the human perception of noise is more sensitive to the noise variance rather than noise mean [1][2][3][4]11]. Therefore, for the purpose of video projection, it is better to display a series of noisy images, which are straight-forward to calculate, rather than trying to optimize a single frame trying to reduce noise. This is achieved by calculating multiple Inverse Fourier Transforms of the image with different random phases. To display these different frames in quick succession, binary-phase ferroelectric devices are used due to their high switching speeds. The algorithm in its simplest form is: Step 1: Add a random phase to a target image: where ( , ) A u v is the input grayscale image and ( , ) u v ϑ is a uniformly distributed random variable.
Step 2: Perform an inverse Fourier Transform: Step 3: Apply a phase mask determined using an optimization algorithm in order to correct aberrations: Step 4: Quantize the hologram to binary phase modulation as in Eq. (10).
Step 5: Repeat steps 1-4 multiple times with different random phases ( , ) u v ϑ and display the results sequentially on the Spatial Light Modulator with a high frame-rate. This approach is extremely efficient at reducing noise, but is not suited for spatiallyvarying aberration correction, because all of the replay field points are processed at the same time in a single Fourier Transform operation.

Piecewise-Corrected OSPR (PC-OSPR)
Since all of the optical components are continuous in their shape, one can expect the optical aberration profile to be a continuous function with respect to spatial coordinates ( , ) u v .
A first inspection of the correction using an array of spots suggested that it changes relatively slowly with distance, shown in Fig. 3. In this case, the correction was designed for the point indicated by the circle (inset figure), but because of the approximate cylindrical symmetry of the system, is also valid within the ring of pixels enclosing approximately 10% of the replay field (the area enclosed between red circles). This is an indication that it might be possible to correct an entire replay field by splitting it up into a small number of regions. Mathematically, if a replay field can be divided into a set of n regions 1 R Rn  for which the aberration-correcting phase mask can be approximated as constant: Then, Eq. (9) can be rewritten: The term inside the curly braces is equivalent to a Discrete Fourier Transform, so we can write: where the mask ( , ) q M u v has been introduced in order to filter only the points that belong to the region Rq , defined as: This important result demonstrates that PWPS and OSPR can be applied over a large replay field using a piecewise approximation consisting of fixed phase masks ( , ) q x y φ . To correct an entire replay field, it is sufficient to correct each region with a separate phase mask and sum the resulting complex holograms. This novel algorithm is termed Piecewise-Corrected OSPR (PC-OSPR).
After the optimization procedure is carried out for a number of different field positions, the experimental points are fit into polynomials to get a smooth variation of parameters. One of these fits corresponding to the variation of the fourth Zernike Polynomial is shown in Fig.  5.
To get the aberration-correcting phase mask for the PWPS method, one then needs to calculate the value of each Zernike coefficient for the particular field position and then perform the summation as indicated in Eq. (2). It is also possible to adapt the Zernike correction for PC-OSPR by sampling this function at particular points. However, this correction proved substantially inferior to the Adaptive-Optical method.  Fig. 6. A schematic view of the feedback loop mechanism: an image of a single point is displayed on a projector, the replay field is then fed back through the webcam to the computer, which performs the correction.

Aberration measurement using an Adaptive-Optical Feedback Loop Mechanism
Alternatively, aberrations can be experimentally characterized using an adaptive optical feedback loop mechanism, such as that illustrated in Fig. 6. We have implemented such a system here which, based only on the feedback from the projector, is able to measure and correct arbitrary spatially varying aberrations. This is essentialy achieved by applying the method developed in [7] repeatedly to contiguous spatial regions. Here, we will only summarize the method.
The projector under test displays an image of a single point at a desired location with some correction parameters 4 1 6 ( ... ) a a . The image of that point is taken with a webcam and the fit function ( ff ) value based on the peak intensity of the point and its physical size is calculated. That procedure can be thought of as a mapping in 13-dimensional space: 4 1 6 ( ... ) a a ff → .
In order to find a perfect correction, one has to seek the global minimum of this function. This is an optimization problem and the solution developed [7] uses a hybrid between Genetic and Steepest Descent Algorithms. Such a correction is performed once for the lifetime of a projector and completes within 5 hours. For the spatially varying case, this procedure is performed multiple times for different field positions and the correction is found in each case. When this system is implemented on a commercial scale, it is possible to substantially speed up the process by correcting multiple points at once and employing state-of-the-art hardware and software.

Distortion measurement and correction
In the general case an optical system will exhibit both aberrations and distortions that need to be corrected. Distortion of a cylindrically-symmetric optical system is always radial, i.e. in order to correct for distortion it is sufficient to know a curve of distorted radius vs. paraxial radius. That measurement was made by displaying an image of concentric circles of known radii, shown in Fig. 7, and measuring the actual radius of each circle when projected. The resulting output curve is presented in Fig. 8. Fig. 7. Image of concentric circles used as an input to the distortion correction algorithm. In order to improve the performance of the algorithm, many images of the same target with different corrections were averaged to get a sharp image. Fig. 8. Measured distortion curve. The line is tangential to y = x at x = 0, indicating that no distortion is present close to the centre, but curves upwards as y rises, indicating pincushion distortion, which drags distorted points outwards.
In order to correct for distortion, it is sufficient to predistort the image before passing it to the hologram generation routine. Figure 9(a) shows the grid of pixels and Fig. 9(b) shows the same grid pre-distorted so as to exactly counteract the distortion of the system.

Region boundary assignment
Having determined optimal aberration-correcting phase masks at 6 radially separated points, the region boundaries can then be assigned. This is done by displaying a grid of single points and algorithmically recognizing which of the 6 correction masks minimizes the aberrations at each point. From this, groups of contiguous points using the same correction mask emerge. Such a grid with one of the corrections applied can be seen in Fig. 10(a). The chosen group of points can be seen in Fig. 10(b) and the assigned region can be seen in Fig. 10(c).  Fig. 11. A schematic view of the holographic projector used. Lenses L1-L3 have the same parameters as indicated in Fig. 4, the laser wavelength is 532 nm

Experimental setup
The holographic projector used in this research is shown in Fig. 11. A 532nm laser beam is expanded and collimated by lenses L1 and L2, modulated by an SLM in reflection mode. Then, the combination of lenses L2 and L3 is used to demagnify the hologram (i.e. expand the image size). Most of the aberrations are introduced by the ball lens L3. Then, the demagnified hologram, formed at the focal plane P2, continues to diffract as it propagates to form a far field image at the webcam. By displaying images at different positions in the replay field and moving the webcam appropriately, the multiple points in the replay field can be measured and then corrected.

Results
The newly developed Piecewise Corrected One-Step Phase Retrieval method is presented here and finally compared with the Pixel-To-Wrapped Phase Summation method. There are various corrections, which are compared against each other. Distortion correction, which is independent of the algorithm, is presented first. Then, there is Piecewise aberration correction done using PWPS and OSPR and finally, a Zemax-PWPS correction.

Distortion correction
We projected a regular array of pixels, which, due to a severe distortion, resulted in a hyperbolic projection seen in Fig. 12(a). In Fig. 12(b), it can be seen that this error is eliminated.

Aberration correction
To visually assess the quality of the correction, a box with 5 pixels forming a cross in the middle (1 pixel apart) is displayed within each different corrected region. This test shape allows visual assessment of various properties of the image, such as quality of vertical and horizontal lines, and the interaction between individual pixels. Figure 13(a) shows the points selected to analyze the correction. Figure 13(b) shows the uncorrected box test image, Fig.  13(c) shows the corrected replay field using Zemax and Fig. 13(d) shows the corrected replay field using the adaptive-optical feedback mechanism. It is observed that all of the images are significantly improved by this correction as seen in the insets 1-6. The correction at points 1-3 is sub-optimal due to severe distortion of the image in those areas. However, as will be demonstrated later this does not noticeably affect the quality of the resultant projected image.

Boundary setting
For each of the 6 correction masks, the methods described in section 4.5 were applied and the resulting set of regions is shown in Fig. 14. It can be seen that there is some degree of cylindrical symmetry that disappears further from the center. This set of masks can be then used to isolate parts of the image before the Fourier Transform operation is performed. In Fig.  15 one of these regions is presented before and after the aberration correction.

Speed comparison
In Fig. 16 we can see the comparison of execution times between different algorithms. As we can see, the novel PC-OSPR algorithm is independent of the number of pixels and executes within 60s on MATLAB and 70ms on CUDA. The PWPS algorithm, however, depends linearly on the number of non-zero image pixels. If every pixel in the image is present, a worst-case scenario, the CUDA implementation takes 230 seconds to execute. MATLAB would need 15 hours to compute such holograms.  The green line at the bottom center of (b)-(e) is zero-order diffraction caused by light specularly reflecting off the glass SLM cover. This could be reduced with improved anti-reflection coatings or by using off-axis holographic projection.

Final correction
The final corrections as applied to a real image are shown in Fig. 17. It can be seen that both Zemax correction, shown in Fig. 17(c) and adaptive-optical corrections, shown in Figs. 17(d) and 17(e) improve the image, with the adaptive-optical method performing better in the outer edges. Moreover, it can be noticed that the OSPR hologram generation method, shown in Fig.  17(e) is significantly better at reducing noise in the image than PWPS method shown in Fig.  17(d).

Real-time operation
Due to its fast execution, the algorithm allowed a real-time operation at up to 12 frames per second. This process is pictured in Visualization 1 (see Fig. 18). The program takes input from one of the monitors, applies the corrections and displays the resultant hologram on the SLM. The program was executed on a mid-range nVidia GTX 760 GPU. Because of various hardware and software limitations, the algorithm had to be executed with only 8 OSPR frames resulting in a slight reduction of image quality. However, when state-of-the-art hardware is employed, it is possible to achieve both exceptional image quality and video frame rate of 30 frames per second.

Conclusions
A novel piecewise variant of the Pixel-to-Wrapped Phase algorithm for correction spatially varying aberrations is presented. When combined with OSPR technique, the novel algorithm offers substantial speed improvement over pre-existing PWPS methods while improving image quality.
A working Adaptive-Optical Feedback Loop mechanism is demonstrated that can characterize spatially-varying aberrations of the holographic projector. This method is contrasted with traditional PWPS aberration correction using ZEMAX ray-tracing software. It is shown that while ZEMAX simulations can reduce these errors, the piece-wise adaptiveoptical approach provides better experimental image quality.
Projection of distortion-and aberration-corrected images over a field of view of 90x45 degrees is demonstrated. Using highly-parallel GPU programing, a continuous operation of the algorithm at a frame rate up to 12 frames per second on a mid-range graphics card can be achieved, which is sufficient for real-time projection.