Optimisation of Phase-Only Computer-Generated Holograms Jinze Sha Department of Engineering University of Cambridge This dissertation is submitted for the degree of Doctor of Philosophy King’s College January 2025 Declaration This thesis is the result of my own work and includes nothing which is the outcome of work done in collaboration except as declared in the preface and specified in the text. It is not substantially the same as any work that has already been submitted, or is being concurrently submitted, for any degree, diploma or other qualification at the University of Cambridge or any other University or similar institution except as declared in the preface and specified in the text. It does not exceed the prescribed word limit for the relevant Degree Committee. Jinze Sha January 2025 Acknowledgements Completing this thesis has been an incredible journey, and I am deeply grateful to everyone who has been a part of it. I would like to thank my supervisor, Prof. Timothy D. Wilkinson, without whose encouragement I wouldn’t even hit the ‘apply’ button for a PhD study; my parents, for always being proud of me; my wife, Mingyue Liu, for her unwavering love and belief in me; my colleagues at CMMPE, Daoming, Fan, Peter, Ralf, Andrew, Youchao, Yue, Jana, Amr, Miguel, Ben, Oana, Roubing, Zhongling, Dilawer, Antoni, Yutaka and Solomon, for all the useful discussions and fun times together in Cambridge, Tim’s farm, and conference sites around the world; my King’s College friends, Jianghan and Yidong, without whom my cooking, cycling, mahjong and snooker skills wouldn’t have soared; and my high school deskmate, Yiwei Chen, without whom my gaming skills wouldn’t have improved. Q: Why did the photon check into a hotel? 1 1A: Because it needed to travel light! Abstract Currently popular three-dimensional (3D) imaging technologies on the consumer market, including 3D cinema, 3D TV, handheld 3D devices (e.g. Nintendo 3DS, HTC Evo 3D) and Virtual Reality (VR) and Augmented Reality (AR) head sets are stereoscopic displays, which fail to provide a natural experience. These systems rely on artificial tricks to simulate depth, leading to discomfort and limited realism. Holography, by contrast, fully reconstructs the optical wavefront of 3D objects and generates a full 3D light field with true depth cues, which has the potential to provide an unparalleled sense of depth and immersion. However, the practical realisation of holographic displays faces significant challenges due to the limitations of existing Computer-Generated Holography (CGH) methods, which often involve trade-offs between computational efficiency and reconstruction quality. This thesis tackles these challenges through developing novel algorithmic approaches aimed at overcoming the obstacles that prevent the holography technology from dominating the 3D display market. The generation of holograms for practical displays is constrained by the need to compute phase-only holograms compatible with current spatial light modulators (SLMs), which are not yet capable of displaying fully complex holograms but can only modulate either phase or amplitude, among which phase modulation is usually preferred. A method of computing phase-only CGH is called phase retrieval. This thesis proposes several novel phase retrieval methods for holographic displays. The Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) method is proposed to mitigate non-linearities in the existing holographic projection system. The proposal of the Limited-memory Broyden-Fletcher- Goldfarb-Shanno (L-BFGS) optimisation algorithm together with the Sequential Slicing (SS) technique offers an efficient approach for 3D phase-only CGH by evaluating individual depth slices iteratively as opposed to evaluating the entire volume. The Multi-Frame Holograms Batched Optimisation (MFHBO) algorithm is proposed to improve the perceived image quality. Additionally, the information capacity of phase-only CGH is explored from an information-theory perspective to investigate the fundamental limits of phase-only CGH. List of Publications [1] Jana Skirnewskaja, Yunuen Montelongo, Jinze Sha, and Timothy D. Wilkinson. Holo- graphic lidar projections with brightness control. In Imaging and Applied Optics Congress 2022 (3D, AOA, COSI, ISA, pcAOP), page 3F2A.6. Optica Publishing Group, 2022 [2] Jinze Sha, Andrew Kadis, Fan Yang, and Timothy D. Wilkinson. Limited-memory bfgs optimisation of phase-only computer-generated hologram for fraunhofer diffraction. In Digital Holography and 3-D Imaging 2022, page W3A.3. Optica Publishing Group, 2022 [3] Andrew Kadis, Benjamin Wetherfield, Jinze Sha, Fan Yang, Youchao Wang, and Tim- othy D. Wilkinson. Effect of bit-depth in stochastic gradient descent performance for phase-only computer-generated holography displays. London Imaging Meeting, 3:36–40, 7 2022 [4] Jinze Sha, Andrew Kadis, Fan Yang, Youchao Wang, and Timothy D. Wilkinson. Multi- depth phase-only hologram optimization using the l-bfgs algorithm with sequential slicing. J. Opt. Soc. Am. A, 40(4):B25–B32, Apr 2023 [5] Jinze Sha, Adam Goldney, Andrew Kadis, Jana Skirnewskaja, and Timothy D. Wilkinson. Digital pre-distorted one-step phase retrieval algorithm for real-time hologram generation for holographic displays. Journal of Imaging Science and Technology, 67(3):030405–1–030405– 1, 2023 [6] Jana Skirnewskaja, Yunuen Montelongo, Jinze Sha, Phil Wilkes, and Timothy D. Wilkin- son. Accelerated augmented reality holographic 4k video projections based on lidar point clouds for automotive head-up displays. Advanced Optical Materials, 12(12):2301772, 2024 [7] Roubing Meng, Jinze Sha, Zhongling Huang, and Timothy D. Wilkinson. Extending FOV of holographic display with alternating lasers. In Peter Schelkens and Tomasz Kozacki, x editors, Optics, Photonics, and Digital Technologies for Imaging Applications VIII, volume 12998, page 129981J. International Society for Optics and Photonics, SPIE, 2024 [8] Jinze Sha, Andrew Kadis, Benjamin Wetherfield, Roubing Meng, Zhongling Huang, Dilawer Singh, Antoni Wojcik, and Timothy D. Wilkinson. Information capacity of phase- only computer-generated holograms for holographic displays. In Peter Schelkens and Tomasz Kozacki, editors, Optics, Photonics, and Digital Technologies for Imaging Applications VIII, volume 12998, page 129980J. International Society for Optics and Photonics, SPIE, 2024 [9] Jinze Sha, Antoni Wojcik, Benjamin Wetherfield, Jianghan Yu, and Timothy D. Wilkinson. Multi frame holograms batched optimization for binary phase spatial light modulators. Scientific Reports, 14(1):19380, Aug 2024 [10] Jinze Sha, Antoni Wojcik, and Timothy D. Wilkinson. Target image phase optimisation hologram generation method for phase-only spatial light modulators. In Practical Holography XXXIX: Displays, Materials, and Applications, volume 13390. International Society for Optics and Photonics, SPIE, 2025 Table of contents List of figures xv List of tables xxi Nomenclature xxiii 1 Introduction 1 2 Background Theory 5 2.1 The Nature of Light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Wave-Particle Duality . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.2 Wave Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Fundamentals of Holography . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Light Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2 Diffraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.3 Spatial Light Modulator (SLM) . . . . . . . . . . . . . . . . . . . 15 2.3 Computer-Generated Holography (CGH) . . . . . . . . . . . . . . . . . . 19 2.3.1 Direct Phase-Only Hologram Generation . . . . . . . . . . . . . . 20 2.3.2 Direct Binary Search (DBS) Algorithm . . . . . . . . . . . . . . . 25 2.3.3 Simulated Annealing (SA) Algorithm . . . . . . . . . . . . . . . . 28 2.3.4 Gerchberg-Saxton (GS) Algorithm . . . . . . . . . . . . . . . . . . 32 2.3.5 One-Step Phase Retrieval (OSPR) Algorithm . . . . . . . . . . . . 35 2.3.6 Adaptive One-Step Phase Retrieval (AD-OSPR) Algorithm . . . . . 38 2.3.7 3D CGH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm 49 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 xii Table of contents 3.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3 Determining the DPD curve . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Applying the DPD Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4 L-BFGS Optimisation of 2D and 3D CGH 61 4.1 Numerical Optimisation Methods . . . . . . . . . . . . . . . . . . . . . . . 61 4.1.1 Optimisation framework . . . . . . . . . . . . . . . . . . . . . . . 61 4.1.2 Gradient Descent . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.1.3 Newton’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.1.4 Quasi-Newton Method . . . . . . . . . . . . . . . . . . . . . . . . 64 4.1.5 Large Scale Quasi-Newton Method: Limited Memory BFGS (L-BFGS) 64 4.2 Phase-only Hologram Optimisation . . . . . . . . . . . . . . . . . . . . . 66 4.3 Target Image Phase Optimisation (TIPO) . . . . . . . . . . . . . . . . . . . 70 4.4 Multi-Depth Phase-Only Hologram Optimisation . . . . . . . . . . . . . . 74 4.4.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5 Multi-Frame Binary-Phase Holograms Batched Optimisation 91 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.1 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.2 Optical Experiment results . . . . . . . . . . . . . . . . . . . . . . 96 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6 Information Capacity of Phase-Only Computer-Generated Holograms 107 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.2.1 Quantised CGH Algorithm . . . . . . . . . . . . . . . . . . . . . . 108 6.2.2 Measurement of Information . . . . . . . . . . . . . . . . . . . . . 111 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.3.1 Targets in the far field (Fraunhofer region) . . . . . . . . . . . . . . 112 6.3.2 Targets in the near field (Fresnel region) . . . . . . . . . . . . . . . 116 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Table of contents xiii 7 Summary and Future Work 123 7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 References 127 List of figures 1.1 A photo of the holographic portrait of Dennis Gabor [17] . . . . . . . . . . 2 2.1 An illustration of the Young’s double-slit experiment [29] . . . . . . . . . . 6 2.2 Illustration of coherence (a) a point source emitting different wavelengths (temporal coherence) and (b) an extended light source consisting of two iden- tical single wavelength light sources separated spatially (spatial coherence) [36] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Structure of the first laser [37] . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Diffraction geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Huygens-Fresnel wavelet principle [45] . . . . . . . . . . . . . . . . . . . 12 2.6 Fresnel and Fraunhofer region [47] . . . . . . . . . . . . . . . . . . . . . . 13 2.7 A single FLC backplane SLM pixel structure [48] . . . . . . . . . . . . . . 15 2.8 Photo of an FLC SLM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.9 Modulation schemes of the SLMs [49] . . . . . . . . . . . . . . . . . . . . 16 2.10 Rotational symmetry in the projection result using the binary phase SLM . 18 2.11 Sample target image of a mandrill [55] . . . . . . . . . . . . . . . . . . . . 19 2.12 Direct phase-only hologram generation method output . . . . . . . . . . . 21 2.13 Illustration of the intuitive reason for the random phase addition [59]. (a) without random phase (b) with random phase . . . . . . . . . . . . . . . . 22 2.14 Output of the improved phase-only hologram generation method . . . . . . 23 2.15 Output of the improved phase-only hologram generation method with binary- phase quantisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.16 DBS algorithm running on the rotationally symmetrical mandrill target . . . 26 2.17 DBS algorithm running on the low resolution target . . . . . . . . . . . . . 27 2.18 SA algorithm running on the low resolution target . . . . . . . . . . . . . . 30 2.19 SA algorithm running on the rotationally symmetrical mandrill target . . . 31 xvi List of figures 2.20 GS algorithm output on the mandrill target . . . . . . . . . . . . . . . . . . 33 2.21 GS algorithm running on the rotationally symmetrical mandrill target . . . 34 2.22 OSPR algorithm running on the rotationally symmetrical mandrill target . . 37 2.23 AD-OSPR algorithm running on the rotationally symmetrical mandrill target 41 2.24 Multi-slice target consisted of 4 different characters at different distances . 43 2.25 Direct phase-only hologram generation with superposition on the 4-slice target 44 2.26 Direct phase-only hologram generation method method’s result on the 4-slice target after binary quantisation . . . . . . . . . . . . . . . . . . . . . . . . 44 2.27 OSPR algorithm’s result on the 4-slice target . . . . . . . . . . . . . . . . . 45 2.28 GS with superposition method’s result on the 4-slice target . . . . . . . . . 45 2.29 GS with sequential slicing algorithm’s result on the 4-slice target . . . . . . 46 2.30 GS with SS algorithm’s NMSE v.s. iteration number plot . . . . . . . . . . 47 2.31 DCGS algorithm’s result on the 4-slice target . . . . . . . . . . . . . . . . 48 2.32 DCGS algorithm’s NMSE v.s. iteration number plot . . . . . . . . . . . . . 48 3.1 Optical setup of the holographic projection system [71] . . . . . . . . . . . 50 3.2 Mechanical components with part numbers of the holographic projection system [71] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3 Determining the DPD curve. (a) Input linear grey-scale ramp. (b) Corre- sponding CGH of (a) with 24-subframe binary phase encoding. (c) Holo- graphic projection replay field of (b). (d) Plot of non-linearity measurement and according pre-distortion curve. . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Validation of DPD curve on the grey-scale ramp. (a) Pre-distorted ramp. (b) Corresponding CGH of (a) with 24-subframe binary phase encoding. (c) Holographic projection replay field of (b). (d) Non-linearity measurement after DPD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.5 Application of DPD on the 10-step strips. (a) 10 strips with equal step of pixel value. (b) CGH of (a). (c) Holographic projection replay field of (b). (d) After DPD of (a). (e) CGH of (d). (f) Holographic projection replay field of (e). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.6 Application of DPD on two sample real-word images. (a) Sample image 1: City Scene [73]. (b) Sample image 2: Horse. (c) Sample image 1 after DPD. (d) Sample image 2 after DPD. . . . . . . . . . . . . . . . . . . . . . . . . 56 List of figures xvii 3.7 Projection output of the two sample images before and after DPD. (a) Replay field of Sample image 1 before DPD (NMSE=0.06139). (b) Replay field of Sample image 2 before DPD (NMSE=0.04309). (c) Replay field of Sample image 1 after DPD (NMSE=0.04920). (d) Replay field of Sample image 2 after DPD (NMSE=0.03635). . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Flowchart of the optimisation process . . . . . . . . . . . . . . . . . . . . 67 4.2 Convergence plot for comparison between the GD, Adam and L-BFGS optimisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3 Reconstructions at each iteration of the L-BFGS optimisation . . . . . . . . 69 4.4 TIPO flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.5 TIPO iterations on the mandrill target (i is the iteration number, Φi, Hi, Ri are the phase of the target field, the phase hologram, and the reconstruction amplitude at iteration i) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.6 TIPO convergence plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.7 Loss between the multi-depth targets (T1 to Tn) and the reconstructions (R1 to Rn) of hologram H . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.8 Optimisation of CGH with sequential slicing (SS) flowchart . . . . . . . . . 78 4.9 Layout of the 4-slice target (z1 = 1 cm, z2 = 2 cm, z3 = 3 cm, z4 = 4 cm) . . 79 4.10 Final NMSE and run time comparison across the three techniques . . . . . 80 4.11 NMSE of each slice plotted against the iteration number for the optimisation based SS technique implemented with MSE and RE loss functions and (a) GD optimiser, (b) L-BFGS optimiser. . . . . . . . . . . . . . . . . . . . . 81 4.12 NMSE of each slice plotted against the iteration number for the (a) GS with SS algorithm, (b) DCGS algorithm. . . . . . . . . . . . . . . . . . . . . . 82 4.13 (a) Average NMSE among all slices, (b) Maximum difference of NMSE across all slices, plotted against the iteration number for the 6 sequential slicing runs using different techniques . . . . . . . . . . . . . . . . . . . . 84 4.14 Comparison of final holograms and reconstructions on the 4-slice target consisted of letters ‘A’, ‘B’, ‘C’ and ‘D’ . . . . . . . . . . . . . . . . . . . 86 4.15 Layout of the non-binary 4-slice target . . . . . . . . . . . . . . . . . . . . 87 4.16 Comparison of final holograms and reconstructions for non-binary target . . 88 5.1 MFHBO flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.2 An example iteration in the optimisation process . . . . . . . . . . . . . . . 94 xviii List of figures 5.3 Convergence of the MFHBO algorithm on the rotationally symmetrical Mandrill target . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.4 Simulation and optical reconstruction results for different number of frames 97 5.5 Sample target image - ‘holography’ ambigram . . . . . . . . . . . . . . . . 98 5.6 Optical results comparison of the proposed MFHBO method against the existing OSPR and AD-OSPR methods . . . . . . . . . . . . . . . . . . . 99 5.7 4-slice target and according reconstruction results . . . . . . . . . . . . . . 101 5.8 Real-life captured image as target field and their reconstruction results . . . 103 6.1 Gerchberg-Saxton (GS) algorithm flowchart . . . . . . . . . . . . . . . . . 109 6.2 Quantisation of phase holograms . . . . . . . . . . . . . . . . . . . . . . . 110 6.3 Quantised Gerchberg-Saxton (GS) algorithm flowchart . . . . . . . . . . . 110 6.4 Del operation on a sample image . . . . . . . . . . . . . . . . . . . . . . . 111 6.5 Holograms generated at bit depths level from 1 to 8 and their according reconstructions in the far field . . . . . . . . . . . . . . . . . . . . . . . . 112 6.6 The average and standard deviation of the far-field reconstruction errors among the 800 target images plotted against the hologram bit depth . . . . 113 6.7 Scatter plot of the far-field reconstruction errors v.s. entropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.8 Scatter plot of the far-field reconstruction errors v.s. delentropy of target im- age among the 800 target images, with different colours indicating hologram bit depth constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.9 Holograms generated at bit depths level from 1 to 8 and their according reconstructions in the near field . . . . . . . . . . . . . . . . . . . . . . . . 116 6.10 The average and standard deviation of the near-field reconstruction errors among the 800 target images plotted against the hologram bit depth . . . . 117 6.11 Scatter plot of the near-field reconstruction errors v.s. entropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.12 Scatter plot of the near-field reconstruction errors v.s. delentropy of target im- age among the 800 target images, with different colours indicating hologram bit depth constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 List of figures xix 6.13 Hologram entropy (solid lines) and NMSE (dotted lines) plotted against the iteration number in a GS run with the initial hologram phase (∠A) being (a) zeros (b) random . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 List of tables 3.1 Non-linearity results before and after DPD . . . . . . . . . . . . . . . . . . 55 3.2 DPD results for sample images . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1 MFHBO runtime (s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2 Quantitative analysis of the optical results in Fig. 5.6 . . . . . . . . . . . . 100 5.3 Quantitative analysis of the optical results in Fig. 5.7 . . . . . . . . . . . . 102 5.4 Quantitative analysis of the optical results in Fig. 5.8 . . . . . . . . . . . . 104 Nomenclature Roman Symbols B Magnetic flux density D Electric flux density E Electric field intensity H Magnetic field intensity J Current density F Fourier Transform operator F−1 Inverse Fourier Transform operator L Loss function P Propagation function from the hologram plane to the reconstruction plane Q Quantisation function A Hologram aperture (in complex form) E Replay field (in complex form) H Phase-only hologram R Reconstruction intensity of the hologram xxiv Nomenclature T Target image j Imaginary unit, √ −1 k Wavenumber, 2π/λ L Loss value P Point on the reconstruction plane R Distance between a point on the reconstruction plane and the origin of the hologram aperture r Distance between a point on the reconstruction plane and a point on the hologram aperture t Time x Horizontal coordinate in the hologram plane y Vertical coordinate in the hologram plane z Distance between the hologram plane and the reconstruction plane Greek Symbols α Horizontal coordinate in the reconstruction plane β Vertical coordinate in the reconstruction plane ε Electric permittivity λ Wavelength of light µ Magnetic permeability ∇ Gradient operator Nomenclature xxv ω Angular frequency ∂ Partial derivative operator φ Phase angle ρ Volume charge density Acronyms / Abbreviations 3D Three-Dimensional AD-OSPR Adaptive One-Step Phase Retrieval AR Augmented Reality CGH Computer-Generated Holography CRT Cathode Ray Tube DPD Digital Pre-Distortion DPD-OSPR Digital Pre-Distorted One-Step Phase Retrieval FFT Fast Fourier Transform FT Fourier Transform GD Gradient Descent GPU Graphics Processing Unit IFFT Inverse Fast Fourier Transform L-BFGS Limited-memory Broyden-Fletcher-Goldfarb-Shanno LUT Look-Up Table MFHBO Multi-Frame Holograms Batched Optimisation xxvi Nomenclature MSE Mean Squared Error NMSE Normalised Mean Squared Error OSPR One-Step Phase Retrieval PDP-TV Plasma display panel televisions PSNR Peak Signal-to-Noise Ratio SLM Spatial Light Modulator SS Sequential Slicing SSIM Structural Similarity Index Measure TFT LCD Thin Film Transistor Liquid Crystal Display TIPO Target Image Phase Optimisation VR Virtual Reality Chapter 1 Introduction The pursuit of three-dimensional (3D) displays has continued at a pace over the past 20 years. Currently, most commercially-available ‘3D display’ products such as 3D cinema, 3D TV, handheld 3D devices (e.g. Nintendo 3DS, HTC Evo 3D) and Virtual Reality (VR) and Augmented Reality (AR) head sets are in fact stereoscopic displays [11] where two different two-dimensional (2D) images are displayed to the left and right eyes respectively, creating a 3D illusion in the brain. Despite their high image quality, the major issue with stereoscopic displays is that they cannot provide real optical defocusing effect in depth [12]. Modern 3D cinemas are able to provide good comfort because polarisation glasses are as light as regular glasses, and the variable defocusing issue can be avoided by the combinations of good design of point of interest in each scene and the according defocusing effect as captured by the camera, so most audience won’t experience much discomfort for around 2 to 3 hours, which has contributed to the general acceptance of 3D films by the public [13]. However, the content, viewing angle and depth of focus of 3D films are fixed after they are captured. To provide an interactive and real-time rendered immersive experience, VR/AR headsets have frequently been advertised as the ‘gateway to the metaverse’ in recent years [14]. However, personal experiences with VR headsets are far from comfortable, not only because of their heavy weight, but also because the display is physically at a very close distance. This creates an unnatural experience because the brain perceives objects at various distances while keeping them all in focus, but in real life, when focusing on a near object, the background naturally blurs. Moreover, the two displays in the VR headset need to be rendered in real-time based on the location and viewing angle of the user, so delays in rendering often cause dizziness and sickness. Hence, the heavy weight, the lack of real depth 2 Introduction of focus and the delay in rendering collectively cause discomfort, dizziness and sickness to VR headsets users [15, 16]. In comparison, holography techniques can produce a full 3D light field, which does not rely on any head mounted device, has true depth of focus, and does not need to be re-rendered according to change in viewer position and viewing angle. Fig. 1.1 A photo of the holographic portrait of Dennis Gabor [17] Holography, taking its name from the Greek word oλoσ (holos), meaning whole, was first introduced in 1948 by Dennis Gabor [18], originally named as wavefront reconstruction [19]. It is a technology which generates 3D images via the diffraction of light. Similar to 2D photography, the earliest holography used a piece of film to record the diffraction pattern, which was then used to reconstruct the 3D field, as shown in Fig. 1.1 (which is a holographic recording of Dennis Gabor himself). After the invention of digital cameras, digital holography emerged. The limitation of both methods is that they require a physical object prior to recording the hologram. In order to generate hologram for objects that do not physically exist, computer-generated holography (CGH) was developed. In CGH, a hologram can be calculated through various algorithmic approaches and then displayed on a spatial light modulator (SLM) that modulates the wavefront of a coherent light source to produce 3D reconstructions. Currently available SLM’s can only modulate either phase or amplitude, so algorithms are needed to compute amplitude-only or phase-only holograms. The classic phase-retrieval algorithms include direct binary search [20], simulated annealing [21] and Gerchberg-Saxton [22]. With the developments in modern numerical optimisation methods and increases in computational power, phase retrieval using numerical optimisation methods has also been found in the literature, such as the gradient descent [23, 24] and its 3 stochastic variations [25, 26, 3]. However, all the existing phase retrieval methods still have some fundamental issues, primarily related to their poor image quality and/or the heavy computation required, the solutions of which are the ultimate goals of this research. This thesis therefore explores the development and optimisation of phase-only CGHs for holographic displays. The research investigates various phase retrieval algorithms and proposes novel methods to enhance the reconstruction quality and computational efficiency of CGHs, and then investigates the fundamental limits of discretised phase holograms from an information theory point of view. The thesis is structured as following. Chapter 2 provides a comprehensive literature review, covering the fundamental theories of light, the principles of holography and the evolution of CGH methods. The chapter reviews various phase retrieval algorithms, including the Direct Binary Search (DBS), Simulated Annealing (SA), Gerchberg-Saxton (GS), One-Step Phase Retrieval (OSPR) and Adaptive One-Step Phase Retrieval (AD-OSPR) algorithms and their adaptations for 3D CGHs, em- phasising the limitations of current algorithms and introducing the motivation for pursuing more advanced CGH techniques. Chapter 3 proposes the Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) method. By experimentally evaluating the non-linearities in the holographic projection system and applying a digital pre-distortion (DPD) curve, this chapter demonstrates significant improve- ments in reconstruction quality and reductions in mean squared error. Chapter 4 introduces the optimisation of phase-only holograms using the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimisation algorithm, and then proposes the novel Target Image Phase Optimisation (TIPO) technique, which optimises the phase of the target image instead of the phase of the hologram. Then the multi-depth phase-only hologram optimisation for 3D targets is investigated, and a novel technique called Sequential Slicing (SS) is proposed, which evaluates the loss for a single slice of the 3D target at each iteration instead of a full evaluation on all slices, reducing computational time while maintaining overall quality and minimising quality imbalances across all slices. Chapter 5 extends the optimisation method on generating time-multiplexing binary-phase holograms and proposes the novel Multi-Frame Holograms Batched Optimisation (MFHBO) technique. By optimising a batch of holograms for time multiplexing, the finite response time of human vision is leveraged to average out noise, resulting in improved visual quality. 4 Introduction The MFHBO algorithm demonstrates much better quality compared to the existing time- multiplexing hologram generation methods, such as OSPR and AD-OSPR. Chapter 6 investigates the information capacity of phase-only CGH. This chapter examines the effects of quantisation on hologram bit depth and their impact on reconstruction quality, and looks for the correlation between the entropy of the target image and the reconstruction error, providing insights into the fundamental limits of discretised CGH. Chapter 7 marks the conclusion and lists potential further work that could be done. Chapter 2 Background Theory This chapter lays down the fundamental theories of optoelectronics and search algorithms, which are essential to the research outlined in the later chapters of this thesis. 2.1 The Nature of Light 2.1.1 Wave-Particle Duality The problem of how light propagates has been troubling scientists for centuries. In the 17th century, Sir Isaac Newton made a significant step. He proposed that light consists of particles, or ‘corpuscles’, with a mass varying with colour, which explained phenomena such as reflection and refraction [27]. In contrast, Christiaan Huygens, a contemporary of Newton, demonstrated that light behaves as a wave, as it is capable of diffraction, which is the bending of light around the edge of an object, leading to the non-sharp edges of shadows [28]. 6 Background Theory Fig. 2.1 An illustration of the Young’s double-slit experiment [29] The wave theory gained significant support in the early 19th century through the experiments by Thomas Young. Young’s double-slit experiment in 1801 provided clear evidence of the wave nature of light by showing that light passing through two slits creates an interference pattern on a screen [30], as illustrated in Fig. 2.1. Augustin-Jean Fresnel further advanced the wave theory by developing a comprehensive mathematical framework to describe light as a wave, explaining phenomena such as polarization and the diffraction of light [31]. This understanding was further reinforced by James Clerk Maxwell in the mid-19th century. In 1864, James Clerk Maxwell organised a set of four equations describing the space and time dependence of the electromagnetic field, which are: ∇×E=−∂B ∂ t (2.1) ∇×H= J+ ∂D ∂ t (2.2) ∇ ·D= ρ (2.3) ∇ ·B= 0 (2.4) 2.1 The Nature of Light 7 where D is the electric flux density, E is the electric field intensity, B is the magnetic flux density, H is the magnetic field intensity, ρ is the volume charge density, and J is the current density [32]. The relation between D and E and between B and H for linear materials are: B= µH (2.5) D= εE (2.6) where µ is the magnetic permeability and ε is the dielectric permittivity of the material [33]. Maxwell’s equations unified electricity and magnetism into a single theory of electromag- netism, predicting that light is an electromagnetic wave that propagates through space [34]. Despite its success, the wave theory could not explain all light-related phenomena. The early 20th century brought a pivotal development with Albert Einstein’s explanation of the photoelectric effect. In 1905, Einstein proposed that light also behaves as particles, or ‘quanta’ (later called photons), which could eject electrons from a metal surface when light is shone upon it [35]. This particle nature of light was critical in explaining observations that wave theory alone could not address and earned Einstein the Nobel Prize in Physics in 1921. These discoveries collectively revealed that light exhibits both wave and particle properties, depending on the experimental context. This wave-particle duality became a cornerstone of quantum mechanics, fundamentally altering human’s understanding of the nature of light. Although the exact nature of light remains unknown, its behaviour is now well understood. 2.1.2 Wave Equation To mathematically describe the propagation of light in free space (i.e. in absence of free charge), the Maxwell equations in Eq. (2.1) - Eq. (2.4) can be combined to derive the wave equation, relating the space and time domain relation of electromagnetic waves propagating in free space: ∇ 2E= µε ∂ 2E ∂ t2 (2.7) A valid solution to Eq. (2.7) is: E= E0e j(ωt−kr) (2.8) 8 Background Theory where ω is the angular velocity of the wave, t is time, r is the propagation distance and k is called the wave number (k = 2π λ , where λ is the wavelength). From Eq. (2.8) we can see that the propagation of light in free space is essentially a phase shift. This suggests that, if we have a coherent light source and a device to manipulate light (called SLM, further explained in Section 2.2.3), we can produce an interference pattern reconstructing the target field we desire, and such method is called holographic projection. 2.2 Fundamentals of Holography 9 2.2 Fundamentals of Holography Holography is a technology that can fully reconstruct the wavefront of 3D objects, which is usually achieved by modulating a coherent light source. This section explains what a coherent light source is and how it is modulated and diffracted. 2.2.1 Light Source The mechanism of holographic projection is to control the propagation of light in a way that, after diffraction, it reconstructs a wavefront that matches the target field. We usually prefer to start from a coherent light source rather than a random one, as it will be a lot more difficult or sometimes impossible to analyse and predict the interference pattern of a random light source. Fig. 2.2 Illustration of coherence (a) a point source emitting different wavelengths (temporal coherence) and (b) an extended light source consisting of two identical single wavelength light sources separated spatially (spatial coherence) [36] The coherence of light refers to the property of light waves where the phase relationship between the waves is consistent over time and space, corresponding to temporal and spatial coherence: • Temporal coherence: Temporal coherence describes the correlation between the phases of a light wave at different points along its propagation direction. It indicates how monochromatic (i.e. single-frequency) a light source is. • Spatial coherence: Spatial coherence describes the correlation between the phases of a light wave at different points across the wavefront, perpendicular to the direction of propagation. It indicates the uniformity of the phase across the wavefront, as illustrated 10 Background Theory in Fig. 2.2. High spatial coherence means that the light waves across different points on the wavefront are in phase. Fig. 2.3 Structure of the first laser [37] The most common coherent visible light source is Laser, which stands for Light Amplification by the Stimulated Emission of Radiation. It was first invented by Theodore Maiman in 1959, with the structure shown in Fig. 2.3 [37–39]. It differs from other sources of light in that it emits coherent light, which is suitable for holographic projection. However, the coherent and monochromatic property of laser also has a side effect of speckle noise in the reconstructed image [40], which is one of the major problems affecting the image quality of holographic projections and has seen lots of efforts to cope with it in the literature [41–44]. 2.2.2 Diffraction Diffraction is the deviation of waves from straight-line propagation without any change in their energy due to an obstacle or through an aperture [19]. The principles of diffraction and interference underpin the essential process of holographic projection, making it possible to accurately recreate complex wavefronts and achieve true 3D visualization. This section delves into how light interacts with apertures, leading to diffractions. Understanding diffraction is essential for holography, as it explains how light can be manipulated to reconstruct three- dimensional light fields. 2.2 Fundamentals of Holography 11 Fig. 2.4 Diffraction geometry To model how light diffracts through a 2D aperture, we first set up a coordinate system as shown in Fig. 2.4, where the aperture is denoted by A(x,y) with x,y denoting the coordinates in the hologram plane, and the diffracted field is denoted by E(α,β ,z) with α,β denoting the coordinates in the reconstruction plane and z denoting the distance from the hologram aperture A(x,y). The R in Fig. 2.4 defines the distance between point P and the origin of the aperture ((x,y) = (0,0)), r defines the distance between point P and a point on the aperture, and θ defines the angle r from the z-axis. Then by trigonometry we can have the following identities: cos(θ) = z r (2.9) R2 = α 2 +β 2 + z2 (2.10) r2 = (α− x)2 +(β − y)2 + z2 (2.11) 12 Background Theory Fig. 2.5 Huygens-Fresnel wavelet principle [45] The Huygens-Fresnel principle states that every point on a wavefront is itself the source of outgoing secondary spherical wavelets, which can be expressed mathematically as follows when r≫ λ [46]: E(α,β ,z) = 1 jλ ∫∫ A(x,y) e jkr r cos(θ)dxdy (2.12) Applying the identities in Eq. (2.9) - Eq. (2.11), Eq. (2.12) becomes: E(α,β ,z) = z jλ ∫∫ A(x,y) e jkr r2 dxdy (2.13) = z jλ ∫∫ A(x,y) e jk √ (α−x)2+(β−y)2+z2 (α− x)2 +(β − y)2 + z2 dxdy (2.14) Unfortunately, Eq. (2.14) cannot be solved analytically except for few specific aperture functions A(x,y), so we have to make some approximations in order to solve for arbitrary A(x,y). The common methods are Fresnel and Fraunhofer approximations for regions depicted in Fig. 2.6. 2.2 Fundamentals of Holography 13 Fig. 2.6 Fresnel and Fraunhofer region [47] Fresnel Approximation √ 1+d = 1+ 1 2 d− 1 8 d2 + · · · (2.15) Fresnel approximation replaces expressions for spherical waves by quadratic-phase exponen- tials, using the binomial expansion of the square root (given in Eq. (2.15)) to approximate r in Eq. (2.13) [46]. Retaining only the first two terms of the expansion gives: r = √ (α− x)2 +(β − y)2 + z2 (2.16) = z √ 1+ ( α− x z )2 + ( β − y z )2 (2.17) ≈ z [ 1+ 1 2 ( α− x z )2 + 1 2 ( β − y z )2 ] (2.18) For the r2 in the denominator of Eq. (2.13), the error introduced by dropping all terms but z is generally acceptably small (i.e. r2 ≈ z2), and for the r appearing in the exponent in the numerator of Eq. (2.13), errors are much more critical [46]. So, by substituting Eq. (2.18) for 14 Background Theory the r in the numerator of Eq. (2.13) and substituting z for the r in the denominator, we have: E(α,β ,z)≈ z jλ ∫∫ A(x,y) e jkz [ 1+ 1 2( α−x z ) 2 + 1 2 ( β−y z )2 ] z2 dxdy (2.19) = e jkz jλ z e j k 2z (α 2+β 2) ∫∫ { A(x,y)e j k 2z (x 2+y2) } e− j 2π λ z (αx+βy)dxdy (2.20) = e jkz jλ z e j k 2z (α 2+β 2)F { A(x,y)e j k 2z (x 2+y2) } (2.21) where F denotes the Fourier Transform, implemented on computers using the Fast Fourier Transform (FFT) function. Such method of including Fourier Transform (FT) in the study of optics is also named ‘Fourier Optics’. Now we have a more simple and solvable expression than Eq. (2.14). And also, as we are only interested in the scaling of relative points at P with respect to each other, it is safe to normalise the multiplier term before the Fourier Transform to 1 [47]. So we can express the diffraction pattern in Fresnel region as: EFresnel region(α,β ,z) = F { A(x,y)e j k 2z (x 2+y2) } (2.22) Fraunhofer Approximation Fraunhofer approximation is very stringent, it assumes that the distance between the light source and the receiving screen are in effect at infinite, with z satisfying the condition of z≫ k(x2 + y2)max 2 (2.23) so that the wave fronts can be treated as planar rather than spherical [32], then the e j k 2z (x 2+y2) term tends to 1, and Eq. (2.22) becomes: EFraunho f er region(α,β ) = F {A(x,y)} (2.24) which suggests that the far field pattern is simply the Fourier Transform of the aperture function. 2.2 Fundamentals of Holography 15 2.2.3 Spatial Light Modulator (SLM) SLMs are critical components in computer-generated holography (CGH). An SLM is a device used to control the amplitude or phase of light waves in a spatially varying manner. SLMs typically consist of a two-dimensional (2D) array of pixels, each of which can modulate the light either passing through or reflected from it. These pixels are usually addressed by electronic signals, allowing precise manipulation of the light wavefront. The modulation can be achieved through various mechanisms, such as liquid crystal (LC) SLMs, magneto-optic SLMs, deformable mirror SLMs, multiple-quantum-well SLMs, or acousto-optic Bragg cells [46]. Fig. 2.7 A single FLC backplane SLM pixel structure [48] The typical structure of a Ferroelectric Liquid Crystal (FLC) SLM pixel is shown in Fig. 2.7. The horizontal dimension of a single pixel, which is called the pixel pitch, is on scale of a few µm. LCs are anisotropic materials, which have different refraction index in different axis. And modulation is achieved by controlling the tilt of the angle of the LCs. For the FLCs, the tilt angle is controlled by the electric field across them. Therefore, the FLCs are sandwiched between the Indium Tin Oxide (ITO) conductor and the aluminium (Al) mirror, aside which is a Very Large Scale Integration (VLSI) circuitry to control the voltage and the electric field across the FLCs, hence the phase of the incident light is modulated. Aligning millions of pixels in a grid produces an SLM which can spatially modulate the wavefront of the incident light. 16 Background Theory Fig. 2.8 Photo of an FLC SLM A photo of an FLC SLM is shown in Fig. 2.8. The SLM used in this research is a binary phase SXGA-R2 ForthDD FLC SLM with a refresh rate of 1440Hz, a pixel pitch of 13.6 µm and a resolution of 1280px×1024px. Fig. 2.9 Modulation schemes of the SLMs [49] Currently available SLMs have one of the four modulation schemes, as illustrated in Fig. 2.9 [49], which are: 2.2 Fundamentals of Holography 17 • Multi-level Amplitude modulators can modulate each pixel from zero transmission (0) to full transmission (1), either continuously or in discrete steps (e.g. nematic liquid crystal display [50], found for example in laptops and many conventional video projectors). • Binary Amplitude modulators can switch each pixel to zero transmission (0) or full transmission (1), but nothing in between (e.g. deformable mirror device [51], ferroelectric liquid crystal display [52], both used in high-end video projectors). • Multi-level Phase modulators can modulate the phase shift imparted by each pixel from 0 to 2π radians, either continuously or in discrete steps (e.g. Nematic liquid crystal devices [53]). • Binary Phase modulators can switch each pixel for a phase shift of either 0 or π radians (e.g. Ferroelectric liquid crystal displays [54]). Among the four modulation schemes, phase modulations are of higher interest for the purpose of holography, because amplitude modulators, either multi-level or binary, block light at the SLM, causing waste of energy and leading to poorer energy efficiency. Moreover, amplitude modulation always has a zero-order (forming a central bright spot), because the average amplitude is always between 0 and 1; on the contrary, phase modulation can suppress the zero-order by designing the hologram to have zero average. As there is no fully complex modulator available yet, we need algorithms to generate phase-only holograms, and such process is called phase retrieval. There are currently many algorithms for such purpose, which will be discussed in Section 2.3. Rotational symmetries in the binary phase modulation The spatial light modulator (SLM) used in this thesis is a binary phase modulator. As the binary phase modulation is purely real (i.e. it’s only switching between 0◦ and 180◦, corresponding to 1 and -1 values of A∗(x,y)), the complex conjugate A∗(x,y) is the same as A(x,y): A∗(x,y) = A(x,y) (2.25) because the Fourier transform of A∗(x,y) is the same as the Fourier transform of A(x,y) E(−α,−β ) = F [A∗(x,y)] = F [A(x,y)] = E(α,β ) (2.26) 18 Background Theory So, at the Fraunhofer region, there is no distinction between the desired image and its 180◦ rotation in the replay field, resulting in a symmetrical conjugate image rotated 180◦ from the target image. (a) Target image (b) Projection result Fig. 2.10 Rotational symmetry in the projection result using the binary phase SLM To demonstrate this phenomenon, the example target image shown in Fig. 2.10a was used to generate the binary-phase hologram using the algorithm known as one-step phase-retrieval (OSPR), which will be explained in detail in Section 2.3.5, and the projection output is shown in Fig. 2.10b. The simplest workaround for this issue is to use only half of the reconstruction field, like the target image in Fig. 2.10a. However, the side effect of this is that half of the energy will be wasted, leading to higher power consumption and heat dissipation. 2.3 Computer-Generated Holography (CGH) 19 2.3 Computer-Generated Holography (CGH) Different from traditional analogue optical holography and digital holography, which both need a physically existing object to record the hologram, the CGH is a method to use computer algorithms to generate holograms without the need for a physical target object. Ideally, holograms can be computed using the inverse Fourier Transform, using the inverse functions of Eq. (2.22) and Eq. (2.24) derived in Section 2.2.2. The inverse Fourier Transforms on most targets will end up with results in complex numbers; however, as mentioned in Section 2.2.3, currently available SLMs cannot achieve fully complex modulation yet. Therefore, computer algorithms are needed to compute either phase-only or amplitude-only holograms, among which the former is preferred, as explained in Section 2.2.3. This section reviews the existing methods in the literature for calculating phase-only holo- grams H, which are the angles of the complex-valued hologram apertures A (i.e. H = ∠A). And the propagation function is unified as P , which can be either the Fraunhofer propa- gation equation in Eq. (2.22) or the Fresnel propagation equation in Eq. (2.24) depending on the distance from the hologram plane to the target field plane. Lastly, R denotes the reconstruction intensity of the hologram (i.e. R = |P[e jH]|2). Fig. 2.11 Sample target image of a mandrill [55] A sample image of a mandrill (shown in Fig. 2.11) is chosen from the University of Southern California Signal and Image Processing Institute (USC-SIPI) Image Database [55] to test and compare the classical phase retrieval algorithms in the literature. As the square of the amplitude, which is the intensity of light, is the only visible component with the human eye [56], the diffracted electric field amplitude (|E|) will be set to match the square root of the target image (T). 20 Background Theory To quantitatively analyse the reconstruction quality, two metrics are used in this section, which are the normalised mean squared error (NMSE) [57] and the structural similarity index (SSIM) [58]. The NMSE is calculated using Eq. (2.27): NMSE(x,y) = 1 n ∗∑(x− y)2 ∑(x)2 (2.27) where x is the target, y is the measured output and n is the number of pixels in x and y. As the NMSE quantifies the total error between the measured output and the target, a lower NMSE value indicates better quality. Different from the NMSE which calculates absolute errors, the SSIM is a perception-based metric that considers the structural similarity as perceived, as defined in Eq. (2.28): SSIM(x,y) = (2µxµy +C1)+(2σxy +C2) (µ2 x +µ2 y +C1)(σ2 x +σ2 y +C2) (2.28) where µx and µy are the mean values of the pixels in x and y respectively, σ2 x and σ2 y are the variances, σxy is the covariance of x and y, C1 = (k1L)2 and C2 = (k2L)2, where L is the dynamic range of the pixel values and k1 = 0.01 and k2 = 0.03 by default [58]. A higher SSIM value indicates better reconstruction quality. 2.3.1 Direct Phase-Only Hologram Generation The simplest method to get a phase-only hologram H is by directly extracting the phase of the reverse propagation from the target field, discarding the amplitude component (e.g. for Fraunhofer propagation, H will simply be the phase of the inverse Fourier transform F−1 of the target field). The pseudocode of this method is shown in Algorithm 1 below: Algorithm 1 Direct phase-only hologram generation method Input: Target image T, Propagation function P (e.g. Fresnel or Fraunhofer propagation) Output: Phase hologram H and its reconstruction intensity R E← √ T A←P−1[E] H← ∠A R← |P[e jH]|2 2.3 Computer-Generated Holography (CGH) 21 where j = √ −1, the ∠ sign takes the arguments of the complex numbers element-wise in the matrix, and e jH converts the angles back to complex numbers. All exponentials, modulus and square-root operators were carried out in an element-wise manner, so that the dimensions of T, A, H, R and E are all the same. The method described in Algorithm 1 was then implemented in MATLAB and run on the sample target image shown in Fig. 2.11, with the distance set to infinity (i.e. in the Fraunhofer region using the propagation formula Eq. (2.24)). The results are shown in Fig. 2.12 below: (a) Discarded amplitude (|A|) (b) Phase Hologram (H) (c) Reconstruction (R) Fig. 2.12 Direct phase-only hologram generation method output After taking the inverse Fourier Transform, the amplitude and the phase of A are shown in Fig. 2.12a and Fig. 2.12b respectively. The next step then discarded the amplitude component (Fig. 2.12a) and used the phase component (in Fig. 2.12b) as the phase hologram H. Then the phase hologram went through a forward propagation and the resulting reconstruction intensity is shown in Fig. 2.12. The reconstruction intensity shown in Fig. 2.12c is very different from the desired target image in Fig. 2.11. This shows that, discarding the amplitude has introduced a significant loss of information. From a signal processing point of view, the peak around the centre in Fig. 2.12a corresponds to low spatial frequency signals, and discarding them causes the reconstruction in Fig. 2.12c to lose low frequency components and effectively becomes an edge detector. The reason for this phenomenon is shown in Fig. 2.13 below. 22 Background Theory Fig. 2.13 Illustration of the intuitive reason for the random phase addition [59]. (a) without random phase (b) with random phase In the case without random phase, shown in Fig. 2.13(a), the low frequency object parts can- not spread light widely over the CGH because the spread angle is proportional to sin−1(λν) where λ and ν are the wavelength and the spatial frequency, respectively. Therefore, the object information of parts (2) and (3) can be recorded on the CGH, whereas that of part (1) cannot be recorded. Conversely, as shown in Fig. 2.13(b), random phase acts as an equivalent of a physical diffuser so that light from all parts of the object can spread over the CGH owing to high spatial frequency provided by the random phase. [59] The improvement method of adding a random phase to the target is demonstrated in the pseudocode below: Algorithm 2 Improved direct phase-only hologram generation method with random phase added to the target field Input: Target image T, Propagation function P (e.g. Fresnel or Fraunhofer propagation) Output: Phase hologram H and its reconstruction intensity R E← √ T×RandomPhase() A←P−1[E] H← ∠A R← |P[e jH]|2 The improved direct phase-only hologram generation method was implemented in MATLAB and produced the results shown in Fig. 2.14 below: 2.3 Computer-Generated Holography (CGH) 23 (a) Discarded amplitude (|A|) (b) Multi-level Phase Holo- gram (H) (c) Reconstruction (R) Fig. 2.14 Output of the improved phase-only hologram generation method As shown in Fig. 2.14c, the reconstruction quality has been greatly improved, although it is still quite noisy. The discarded amplitude of the hologram (shown in Fig. 2.14a) is a lot more uniformly distributed than the one in Fig. 2.12a, so that the loss of information is more evenly spread across all spatial frequencies, leading to the much better reconstruction quality. However, the reconstruction quality is still quite noisy. The reconstruction in Fig. 2.14c has an NMSE of 1.0228×10−6 and an SSIM of 0.1603. Moreover, additional error will be introduced during the quantisation step, which is necessary for the phase hologram to be displayed on SLMs with limited bit depth. As the SLM used in this thesis is a binary phase SLM, which has a rotational symmetry property as explained in Section 2.2.3, a rotationally symmetrical target image is specifically designed, as shown in Fig. 2.15a. 24 Background Theory (a) Target image for the binary- phase SLM (b) Binary-phase Hologram (c) Reconstruction Fig. 2.15 Output of the improved phase-only hologram generation method with binary-phase quantisation The binary phase hologram in Fig. 2.15b is generated by adding an additional binary quanti- sation step (Q) on the phase hologram computed using Algorithm 2, which is achieved by rounding all phases to 0 and π radians, using Eq. (2.29): θbinary quantised { 0 −1 2π ⩽ θ < 1 2π π θ <−1 2π or θ ⩾ 1 2π (2.29) where θ is the input phase value in range [−π,π] and θbinary quantised is the output quantised binary phase value. The reconstruction of the binary-phase hologram is shown in Fig. 2.15c, which has an NMSE of 4.5452×10−7 and an SSIM of 0.0603. To improve the reconstruction quality, better algorithms are needed. The following sections explore predecessors efforts in quality improvement. 2.3 Computer-Generated Holography (CGH) 25 2.3.2 Direct Binary Search (DBS) Algorithm Direct Binary Search (DBS) algorithm [20] is an algorithm that generates the hologram by randomly flipping each pixel in the SLM between binary states (0 and π) one by one for many times in order to minimise the difference between the hologram’s reconstruction intensity R and the target image T . The detailed algorithm is described in Algorithm 3 below: Algorithm 3 Direct Binary Search (DBS) algorithm Input: Target image T, Propagation function P , Loss function L (e.g. mean-squared error), Number of iterations N Output: Phase hologram H and its reconstruction intensity R // Start with a random hologram with a size matching T H← Rand(Size(T)) R← |P[e jH]|2 L←L [R,T] for n = 1 to N do // Flip a random pixel in the hologram Hn← FlipRandomPixel(H) // Calculate the loss function for the new hologram Rn← |P[e jHn]|2 Ln←L [Rn,T] // Compare the new loss with the old one if Ln < L then // Accept the new hologram if loss is lower H←Hn R← Rn L← Ln end if end for Although the DBS algorithm is specifically suited for generating binary phase holograms, it can also be adapted for generating multi-level phase holograms, by representing each level as binary numbers at the cost of more computation. 26 Background Theory (a) Target image (1024px× 1024px) (b) Binary-phase Hologram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.16 DBS algorithm running on the rotationally symmetrical mandrill target DBS algorithm can sometimes find very accurate holograms if the run is lucky; however, it is extremely slow due to the large number of iteration required (as shown in Fig. 2.16d, even 105 iterations have not reached a good convergence). The example run on the target image with resolution 1024px× 1024px in Fig. 2.16a took more than one hour to finish the 105 iterations, which is still far from reaching a good convergence. The binary-phase hologram produced is shown in Fig. 2.16b, and its corresponding reconstruction in Fig. 2.16c has an NMSE of 8.0717×10−7, which is 78% larger than the NMSE of the direct phase-only hologram generation method(4.5452×10−7). The reconstruction in Fig. 2.16c has an SSIM of only 0.0076, indicating a very low structural similarity to the target image, and is much 2.3 Computer-Generated Holography (CGH) 27 worse than the direct phase-only hologram generation method, which has an SSIM of 0.0603, as discussed in Section 2.3.1. (a) Target image (128px × 128px) (b) Binary-phase Hologram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.17 DBS algorithm running on the low resolution target As the DBS algorithm only flips one pixel per iteration, it naturally takes significantly longer to generate holograms with higher resolution. To test the programme on a smaller image for much quicker convergence, another target image has been designed, as shown in Fig. 2.17a, which is also rotationally symmetrical as it is used for the binary-phase hologram generation. After 106 iterations, which took 10 minutes, the hologram generated is shown in Fig. 2.17b and its resulting reconstruction is shown in Fig. 2.17c, which has an NMSE of 9.6881×10−6 and an SSIM of 0.2871. The NMSE v.s. iteration plot in Fig. 2.17d shows that it reaches a 28 Background Theory good convergence at around 2×105 iterations, which needs around 2 minutes. Moreover, the curve of NMSE is monotonically decreasing with iteration number, as only holograms with better results are accepted during the iterations. In summary, the DBS algorithm is a slow but working algorithm for binary phase hologram generation. The programme running time scales up significantly when the target image’s resolution gets higher. Furthermore, as the DBS algorithm only cares about local optimality at each iteration, it is a greedy algorithm that only follows the steepest descent route, which can easily get trapped in a local minimum where flipping any bit will not result in better reconstructions. Another consequence of the random nature is that the generated hologram will be different at each run, so the quality of the resulting reconstruction (R) will depend on how ‘lucky’ each run is. The Simulated Annealing (SA) algorithm [21] in the next session aims to resolve this issue. 2.3.3 Simulated Annealing (SA) Algorithm The Simulated Annealing (SA) algorithm [21] is a variant of the DBS algorithm. It adopts a probabilistic approach to avoid the steepest gradient descent. Its name derives from the fact that it approximates the recrystallisation process during metal annealing and is particularly well-suited to avoiding the trap of local minima [60]. To implement this idea, we need a function (Z ) to calculate the probability of the hologram (H), and a threshold pt to decide whether the probability is high enough for the according hologram to be accepted. In this thesis, the probability is selected to be a random function and the threshold is chosen to be 0.9. The pseudocode for this algorithm is listed in Algorithm 4. 2.3 Computer-Generated Holography (CGH) 29 Algorithm 4 Simulated Annealing (SA) algorithm Input: Target image T, Propagation function P , Loss function L , Number of iterations N, Probability function Z , Probability threshold pt Output: Phase hologram H and its reconstruction intensity R // Start with a random hologram with a size matching T H← Rand(Size(T)) R← |P[e jH]|2 L←L [R,T] for n = 1 to N do // Flip a random pixel in the hologram Hn← FlipRandomPixel(H) // Calculate the loss function for the new hologram Rn← |P[e jHn]|2 Ln←L [Rn,T] // Compare the new loss with the old one if Ln < L then // Accept the new hologram if loss is lower H←Hn R← Rn L← Ln else // Calculate the probability of the hologram pn←Z [Hn] if pn > pt then // Accept the new hologram if the probability exceeds the threshold H←Hn R← Rn L← Ln end if end if end for 30 Background Theory (a) Target image (128px × 128px) (b) Binary-phase Hologram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.18 SA algorithm running on the low resolution target An implementation of SA algorithm with pt = 0.9 was carried out on the low resolution target in Fig. 2.18a, and the resulting binary-phase hologram and its reconstruction are shown in shown in Fig. 2.18b and Fig. 2.18c respectively. From the NMSE v.s. iteration plot in Fig. 2.18d, it can be seen that, instead of the monotonic decrease observed in Fig. 2.17d for DBS algorithm, the SA algorithm has occasional rises in NMSE, which happens when the probability pn exceeds the threshold pt . The final NMSE was recorded to be 1.0542×10−5 and the final SSIM was 0.2750, which are both slightly worse than the DBS algorithm in this case. Due to the probabilistic nature of the SA algorithm, although it can avoid being trapped in local optimal points, the ‘jump backs’ can also cause delays in convergence. 2.3 Computer-Generated Holography (CGH) 31 (a) Target image (1024px× 1024px) (b) Binary-phase Hologram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.19 SA algorithm running on the rotationally symmetrical mandrill target Then the SA algorithm was run for the mandrill target in Fig. 2.19a. The convergence plot in Fig. 2.19d shows that it did not converge within 105 iterations, which took around one hour. The binary phase hologram generated is shown in Fig. 2.19b and its corresponding reconstruction intensity is shown in Fig. 2.19c, which has an NMSE of 8.2054×10−7 and an SSIM of 0.0073, which are slightly worse than the DBS algorithm results in Fig. 2.16c with an NMSE of 8.0717×10−7 and an SSIM of 0.0076. Both the DBS and the SA algorithms rely on flipping only a single pixel per iteration, which is very inefficient. A better algorithm should change the values of most pixels at every iteration for better efficiency, and to converge within much fewer iterations for lower computational 32 Background Theory time. The Gerchberg-Saxton (GS) algorithm [22] is a classical example, which will be further explained in Section 2.3.4. 2.3.4 Gerchberg-Saxton (GS) Algorithm The Gerchberg-Saxton (GS) algorithm [22] is a revolutionary algorithm and is much better and more robust than the algorithms introduced in the previous sections. Although being more than 50 years old, the GS algorithm is still frequently used and has lots of variants [61–63]. It functions by iteratively determining the phase profile of the hologram required to reconstruct a target image, looping between the hologram and the reconstruction plane, and applying constraints to each plane accordingly during each iteration. GS algorithm is very easy to implement, its pseudocode is shown in Algorithm 5. Algorithm 5 Gerchberg-Saxton (GS) Algorithm Input: Target image T, Propagation function P , Number of iterations N, Initial phase Φ (e.g. random, zeros, or other patterns) Output: Phase hologram H and its reconstruction intensity R // Initiate E with amplitude √ T and initial phase Φ E← √ T× e jΦ for n = 1 to N do // Compute the hologram plane A←P−1[E] // Apply the phase-only constraint at the hologram plane A← e j∠A // Compute the propagation for the new hologram E←P[A] // Apply the target field amplitude constraint at the reconstruction plane E← √ T× e j∠E end for H← ∠A R← |P[A]|2 2.3 Computer-Generated Holography (CGH) 33 (a) Target image (512px × 512px) (b) Multi-level Phase Holo- gram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.20 GS algorithm output on the mandrill target The GS algorithm described in Algorithm 5 was implemented in MATLAB and was first run on the mandrill target image in Fig. 2.11, and the output results are shown in Fig. 2.20. Fig. 2.20b is the multi-level hologram generated after 30 iterations, and its reconstruction is shown in Fig. 2.20c. The reconstruction in Fig. 2.20c has an NMSE of 2.6612×10−8 and an SSIM of 0.7940, which are both much better than the single-iteration direct phase-only hologram generation method’s result in Fig. 2.14c (with an NMSE of 1.0228×10−6 and an SSIM of 0.1603). The NMSE v.s. iteration plot in Fig. 2.20d shows that the GS algorithm converged quickly, providing very good result in tens of iterations, which is much fewer than the DBS and SA algorithms. Although the GS algorithm is more computationally expensive 34 Background Theory at each iteration, as it needs to compute both a forward and an backward propagation, leading to two Fourier transforms every iteration, the GS algorithm is still much faster and provides much better reconstruction quality than the DBS and SA algorithms. (a) Target image (1024px× 1024px) (b) Binary-phase Hologram (c) Reconstruction (d) NMSE v.s. iteration plot Fig. 2.21 GS algorithm running on the rotationally symmetrical mandrill target Then the GS algorithm was adapted to generate binary-phase holograms, for use on the binary- phase SLM in this thesis. The change was simply implemented by adding a quantisation function (Q) when applying the phase-only constraint at the hologram plane (i.e. the line ‘A← e j∠A’ in Algorithm 5 is changed to ‘A← e jQ[∠A]’). The results are shown in Fig. 2.21. Fig. 2.21d shows good convergence within 20 iterations. The resulting binary-phase hologram is shown in Fig. 2.21b and its corresponding reconstruction in Fig. 2.21c has an NMSE 2.3 Computer-Generated Holography (CGH) 35 of 1.6968e-07 and an SSIM of 0.0619, which are both better than the direct phase-only hologram generation method’s result in Fig. 2.15 (with an NMSE of 4.5452×10−7 and an SSIM of 0.0603). In summary, the GS algorithm is quick and robust. On my laptop computer of model ASUS ROG Zephyrus M16, which has a CPU of model i7-11800H and a GPU of model RTX3060, the 30 iterations took 1.5 seconds to complete. It reached convergence in tens of iterations. However, as it is still iterative, generating holograms in real-time is still a challenge, and the reconstruction still suffers from noise. 2.3.5 One-Step Phase Retrieval (OSPR) Algorithm The OSPR algorithm was first demonstrated by Cable and Buckley [41]. It is a solution to high-quality hologram reconstruction that relies on time multiplexing of holograms, exploiting the response time of eye in order to reduce noise in the replay field [49]. The random noise is averaged by the eye, while the target image stays, so that the average noise can be reduced. The perceived noise is lessened by the temporal average detected by the eye, rather than computational optimisation of the hologram [49]. The pseudocode for OSPR is shown in Algorithm 6 below. 36 Background Theory Algorithm 6 One-Step Phase Retrieval (OSPR) algorithm Input: Target image T, Propagation function P , Number of sub-frames S, Quantisation function Q Output: List of phase holograms H[1 . . .S] // Compute a list of hologram sub-frames for i = 1 to S do E← √ T· RandomPhase() A←P−1[E] H[i]←Q[∠A] end for // Then display the sub-frames on the phase modulator sequentially i← 1 while True do Display(H[i]) i← i+1 if i > S then i← 1 end if end while When generating the list of holograms, it repetitively computes the inverse propagation of the target amplitude (which is the square-root of the target intensity T) multiplied by different random phases, for a total of S times to generate S hologram sub-frames (H[1 . . .S]). The computation of each hologram sub-frame is the same as the direct phase-only hologram generation method in Algorithm 2 discussed in Section 2.3.1. Then the S hologram sub- frames are displayed sequentially on a SLM having a refresh rate being so fast that the average reconstruction intensity is perceived by the human eyes. As currently available fast SLMs are binary-phase modulators, an example run on the rotationally symmetrical target previously used in Fig. 2.16a was carried out for 24 sub-frames (S = 24). The results are summarised in Fig. 2.22. 2.3 Computer-Generated Holography (CGH) 37 H1 H2 H24 Hencoded R1 R2 R24 Ravg … , , , , … { { } } T Fig. 2.22 OSPR algorithm running on the rotationally symmetrical mandrill target In Fig. 2.22, a total of 24 binary-phase hologram sub-frames (H1,H2, . . . ,H24) were generated for the target image T. For easier data transfer and to fit with the common display signal formats, the 24 binary-phase hologram sub-frames are encoded into a single file with 8 bit depth and RGB (red-green-blue) channels, so that each of the (8× 3 =)24 bit-planes corresponds to a single binary-phase hologram sub-frame. For a quantitative analysis, the reconstruction intensities of the hologram sub-frames are computed in R1,R2, . . . ,R24 respectively, whose average is Ravg in Fig. 2.22. The average reconstruction intensity Ravg has an NMSE of 9.8632×10−8 and an SSIM of 0.1321, which are both significantly better than the GS algorithm (NMSE=1.6968× 10−7, SSIM=0.0619) and the direct phase-only hologram generation method (NMSE=4.5452×10−7, SSIM=0.0603). The major advantage of the OSPR algorithm is that it is fast. It is non-iterative and requires only one Fourier Transform per frame, taking less than a second to generate a 24-frame hologram set. Its non-iterative nature also allows it to be parallelised to further improve computation speed, which is crucial for Light Blue Optics who made a real-time holographic laser projector commercially available in 2010 [64], although the product was later dis- continued for financial reasons. The downside of this algorithm is that the sub-frames are 38 Background Theory independent from each other, and the final reconstruction output is still subject to some noise, as they are only more randomly distributed instead of being reduced. There is a variant of improvement on the OSPR algorithm, called Adaptive OSPR (AD-OSPR), to be introduced in Section 2.3.6. 2.3.6 Adaptive One-Step Phase Retrieval (AD-OSPR) Algorithm The AD-OSPR algorithm [65] is a variant of the OSPR algorithm. It aims to improve the reconstruction quality without introducing a significant amount of additional computational cost. It functions in such a way that, when computing the second sub-frame onwards, it subtracts the average reconstruction from the target image to get the error, so that it can compensated in the next iteration. To help explain the process in detail, a pseudocode is written for the AD-OSPR algorithm, as shown in Algorithm 7. 2.3 Computer-Generated Holography (CGH) 39 Algorithm 7 Adaptive One-Step Phase Retrieval (AD-OSPR) algorithm Input: Target image T, Propagation function P , Number of sub-frames S, Quantisation function Q Output: List of phase holograms H[1 . . .S] // Compute a list of hologram sub-frames T[1]← T Rtotal ← 0 for i = 1 to S do E← √ T[i]· RandomPhase() A←P−1[E] H[i]←Q[∠A] // Compute the reconstruction intensity R← |P[e jH[i]]|2 // Carry out energy conservation to match the total energy of the target image R← R · √ sum(T2) sum(R2) // squares and sums are taken element wise // Compute the total reconstruction intensity so far Rtotal ← Rtotal +R // Compute the new target for the next iteration for [x,y] = [1,1] to size(T) do // loop among all the pixels if (i+1) ·T[x,y] > Rtotal[x,y] then T[i+1][x,y]← (i+1) ·T[x,y]−Rtotal[x,y] else T[i+1][x,y]← 0 end if end for end for // Then display the sub-frames on the phase modulator sequentially i← 1 while True do Display(H[i]) i← i+1 if i > S then i← 1 end if end while 40 Background Theory In the pseudocode in Algorithm 7, the target intensity for the first iteration (T[1]) is initialised as the input target image T, and the total reconstruction intensity for the holograms generated up to each iteration (Rtotal) is initialised as 0. Then the for loop starts from the same routine as the OSPR algorithm, multiplying a random phase to the square root of target intensity, taking an inverse propagation, and applying the binary phase quantisation. Then, the total reconstruction is computed by forward propagating the quantised hologram, conserving the energy to match the target image, and adding to the total reconstruction of the last iteration. The total reconstruction is then used to compute the target intensity for the next iteration, by subtracting the total reconstruction from (i+1) times the target image, with an if-else statement to avoid negative intensity values. To provide a visual illustration, an example run of the AD-OSPR algorithm was carried out on the mandrill target as shown in Fig. 2.23. 2.3 Computer-Generated Holography (CGH) 41 Ti Hi Ri Ravg, ii 1 -2× + -2× + 2 + 3 ×(i-1) + ÷i + … 24 … … … … ×(i-1) ÷i ×(i-1) ÷i … Fig. 2.23 AD-OSPR algorithm running on the rotationally symmetrical mandrill target In Fig. 2.23, each row of images correspond to one iteration that is indexed in the ‘i’ column, where only iterations 1, 2, 3, 24 are shown here, while iterations 4-23 are omitted to save space. At each iteration, the binary-phase hologram Hi is computed from its target Ti, and the corresponding reconstruction is shown in Ri. Instead of showing the Rtotal,i, the Ravg,i = Rtotal,i i is shown here as Rtotal,i has increasing dynamic range at each iteration. 42 Background Theory For the first iteration, the target image of the same one as in Fig. 2.22 is used, and the first average reconstruction Ravg,1 is just the first reconstruction R1. Then, from the second iteration onwards, the target intensity Ti is updated following the pseudocode in Algorithm 7. After 24 iterations, the 24 hologram sub-frames are computed and their average reconstruction (Ravg,24) had an NMSE of 9.1161×10−8 and a SSIM of 0.1992, which are both better than the OSPR algorithm in Section 2.3.5 who achieved an NMSE of 9.8632×10−8 and an SSIM of 0.1321. The programme running time of the AD-OSPR algorithm is nearly the same as the OSPR algorithm, both being around 0.7 second. Therefore the AD-OSPR algorithm is quite an effective improvement on the original OSPR algorithm, achieving a 7.6% reduction in NMSE and 50.8% improvement in SSIM without adding much computation. The disadvantage is however, that it cannot be parallelised, unlike the OSPR algorithm, as it needs to have the total reconstruction intensity result from the previous iteration for use in the next iteration. 2.3 Computer-Generated Holography (CGH) 43 2.3.7 3D CGH Section 2.3.1 - Section 2.3.6 described several algorithms to generate a phase hologram for a single target image. However, the major benefit of holography is that it can produce a true 3D light field, much more than just a single slice 2D image. Then the problem arises as how to generate a hologram for 3D targets to make full use of the major benefit of holography. The simplest method is to slice the 3D target into a set of 2D layers, like computed tomography (CT) scanning, and then generate a hologram so that its Fresnel propagation (in Eq. (2.22)) at each depth (z) matches the according layer. This subsection therefore reviews how the current methods in Section 2.3.1 - Section 2.3.6 can be adapted to multi-depth hologram generation. zz 4z3 z2z1 T4T3T2T1 Fig. 2.24 Multi-slice target consisted of 4 different characters at different distances For fair tests on algorithms, an example 3D target field consisted of 4 slices has been created using alphabets, as shown in Fig. 2.24. T1 to T4 have the same resolution of 512px×512px, and the distances z1 to z4 are set to 1,2,3,4cm respectively. 44 Background Theory Direct Phase-Only Hologram Generation with Superposition To adapt the direct phase-only hologram generation method in Section 2.3.1 to compute multi-depth hologram, firstly, a set of holograms are generated corresponding to each layer of the target field respectively, using the inverse of the Fresnel propagation function in Eq. (2.22). Then, based on the principle of superposition, the set of holograms are added up to form the final hologram, whose phase is directly extracted to be the phase-only hologram while the amplitude component is discarded. R4R3R2R1H Fig. 2.25 Direct phase-only hologram generation with superposition on the 4-slice target The simulation results shown in Fig. 2.25 demonstrates the effectiveness of the direct phase- only hologram generation method. The reconstructions at each depth are legible, despite the presence of some noise. However, after applying a binary-phase quantisation on the hologram, the reconstruction quality deteriorates significantly, as shown in Fig. 2.26. R4R3R2R1H Fig. 2.26 Direct phase-only hologram generation method method’s result on the 4-slice target after binary quantisation To improve the reconstruction quality, the time-multiplexed OSPR algorithm (in Sec- tion 2.3.5) is adapted to produce multi-frame holograms for 3D targets. 2.3 Computer-Generated Holography (CGH) 45 OSPR algorithm adaptation R4R3R2R1H Fig. 2.27 OSPR algorithm’s result on the 4-slice target The OSPR algorithm in Section 2.3.5 was adapted to generate multi-depth targets by doing the direct phase-only hologram generation method in the previous paragraph for a total of S times to form S sub-frames. The results of the example run with S = 24 is shown in Fig. 2.27. The average reconstruction has much less noise than the single frame one in Fig. 2.26. However, the contrast is not high enough. Therefore, one of the major objective of this thesis is to search for time-multiplexed binary-phase holograms with better reconstruction quality than the existing methods, which will be further explored in Chapter 5. GS algorithm adaptations The GS algorithm in Section 2.3.4 can also be adapted to generate phase-only holograms for multi-depth 3D targets. Three different adaptations are reviewed in the following paragraphs. GS algorithm adaptation 1 - Superposition: Similar to the OSPR adaptation, the first method of adapting the GS algorithm for 3D CGH is by superposition. Firstly, a set of holograms are generated individually using the GS algorithm on each slice of the 3D target, and then those holograms are superposed into a total hologram whose phase is then used as the final phase hologram. R4R3R2R1H Fig. 2.28 GS with superposition method’s result on the 4-slice target 46 Background Theory The GS with superposition method was run on the target field in Fig. 2.24 and produced the result in Fig. 2.28. Similar to the observations with the OSPR adaptation, the superposition leads to defocusing between each slice, giving rise to the background noise. The advantage of this method is its simplicity of implementation, while the disadvantages of this method are its poor reconstruction quality and its high number of iterations, which grows linearly with the number of slices. GS algorithm adaptation 2 - GS with Sequential Slicing (GS-SS): The second adaptation is called GS-SS, Instead of propagating to a fixed target image at each iteration, the GS algorithm is modified to propagate to a different distance at each iteration, and the target amplitude constraint of the according distance is applied (i.e. line E← √ T× e j∠E in Algorithm 5 becomes E← √ Ti%n× e j∠E, where i is the iteration number, n is the total number of slices and the % sign takes the remainder of the division). Such method is named sequential slicing, where the algorithm sweeps through the slices one by one sequentially during its iterations. R4R3R2R1H Fig. 2.29 GS with sequential slicing algorithm’s result on the 4-slice target The GS with sequential slicing (SS) algorithm was run on the target field in Fig. 2.24 for 100 iterations, and the result is shown in Fig. 2.29. The interesting phenomena is that, as the 100 iterations terminated at the 4th layer, the 4th reconstruction (R4) has the best quality among all slices, and the reconstruction quality deteriorates as the slice index reduces. 2.3 Computer-Generated Holography (CGH) 47 0 20 40 60 80 100 iteration number 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 NM SE 1e 6 GS with SS (Slice 1) GS with SS (Slice 2) GS with SS (Slice 3) GS with SS (Slice 4) Fig. 2.30 GS with SS algorithm’s NMSE v.s. iteration number plot As a further investigation, the NMSE of each slice is plotted against iteration number in Fig. 2.30. It shows that the iterations did not converge. The NMSE curves are oscillating severely. When the NMSE of one slice decreases, the NMSE of all other slices increases. Applying the amplitude constraint at a single depth worsen all the other slices. There is a solution to this issue in the literature, called Dynamic Compensatory GS (DCGS) [63], as introduced in the following paragraph. GS algorithm adaptation 3 - Dynamic Compensatory GS (DCGS): The DCGS algo- rithm ‘softens’ the target field amplitude constraint [63]. Instead of forcing the amplitude of the reconstructed field to the target amplitude directly, it allows a fraction α of the orig- inal amplitude to be retained (i.e. E← √ Ti%n× e j∠E becomes E← [α×|E|+(1−α)× √ Ti%n]× e j∠E), where α is adjusted dynamically at each iteration. The DCGS algorithm was implemented and run on the 4-slice target in Fig. 2.24. The result is shown in Fig. 2.31. 48 Background Theory R4R3R2R1H Fig. 2.31 DCGS algorithm’s result on the 4-slice target Fig. 2.31 shows the effectiveness of the DCGS algorithm, where all slices have good quality instead of the huge inter-slice quality imbalance observed in Fig. 2.29. The NMSE of each slice is plotted against the iteration number in Fig. 2.32. 0 20 40 60 80 100 iteration number 2.0 2.5 3.0 3.5 4.0 4.5 NM SE 1e 6 DCGS (Slice 1) DCGS (Slice 2) DCGS (Slice 3) DCGS (Slice 4) Fig. 2.32 DCGS algorithm’s NMSE v.s. iteration number plot Fig. 2.32 shows a much better convergence than Fig. 2.30. Although the bounce-backs still present (i.e. when the NMSE of one slice decreases, those of other slices increase), the overall trend of NMSE is decreasing as the iterations continue. Chapter 3 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm Note: The work in this chapter has been published in Ref. [5] In a computer-generated holographic projection system, the image is reconstructed via the diffraction of light from a spatial light modulator. In this process, several factors can contribute to non-linearities between the reconstruction and the target image. This chapter first evaluates the non-linearity of the overall holographic projection system experimentally, using binary phase holograms computed by the OSPR algorithm introduced in Section 2.3.5. It then applies a digital pre-distortion (DPD) method to correct for the non-linearity. Both a notable increase in reconstruction quality and a significant reduction in mean squared error were observed, which proves the effectiveness of the proposed DPD-OSPR algorithm. 3.1 Introduction Chapter 2 reviewed the real-time CGH method called OSPR (in Section 2.3.5), which was fast enough for real-time holography, but its reconstruction quality can still be potentially improved. Hence, a computationally inexpensive method is needed to improve the recon- struction quality whilst maintaining the real-time property of the OSPR algorithm. 50 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm Several factors contribute to the non-linearities between the reconstruction of hologram and the target image, including the calculation and quantisation of the hologram, the modulation of the light and the imperfections in the optical setup. This chapter proposes the DPD-OSPR algorithm. DPD is a technique originally used to compensate for the non-linearities of power amplifiers in wireless communication systems. Here, the DPD is carried out on the holographic projection system to compensate for non-linearities between input and output intensities. By measuring the non-linearity experimentally and applying the according pre- distortion curve on target images, DPD can be implemented via a one-to-one correction curve or through a look-up table (LUT) which allows the relationship between the input and output to be adjusted without requiring any heavy computation. The intuition of the proposed DPD-OSPR algorithm for CGH originates from the gamma correction method for conventional displays, such as cathode-ray tube (CRT) monitors [66], plasma display panel televisions (PDP-TV) [67] and thin-film-transistor liquid-crystal displays (TFT LCD) [68, 69]. Gamma correction for conventional displays were originally developed to mimic the perceptual response of human vision [70]. The work presented here builds logically on this approach and extends its application to holographic displays. 3.2 Experimental setup Fig. 3.1 Optical setup of the holographic projection system [71] The holographic projector used in this experiment is a Fourier projection system developed by Freeman [71], as shown in Fig. 3.1. The design is consisted of a diode-pumped solid-state (DPSS) 532 nm 50mW laser source, focussed down by an aspheric singlet (A), the focus of which becomes the diffraction limited point source (B) for the projector. The beam then passes through a polarizing beam splitter cube (C) to a collimating lens (D), which illuminates the SLM (E). The SLM is a binary phase SXGA-R2 ForthDD ferroelectric Liquid crystal on 3.2 Experimental setup 51 silicon (LCOS) micro-display with a refresh rate of 1440Hz, a pixel pitch of 13.6 µm and a resolution of 1280px×1024px. An aperture at point (F) spatially filters out the other orders, leaving only one first order, which is then magnified up by a finite conjugate lens group (G) to produce an image, of the required size, on the screen (H). [71] Fig. 3.2 Mechanical components with part numbers of the holographic projection system [71] The mechanical components are listed in Fig. 3.2 with part numbers. The holograms displayed on the SLM are generated using the OSPR algorithm [41], and as the SLM is a binary-phase modulator, each hologram sub-frame needs to be binary quantised. Then each group of the 24 binary-phase hologram sub-frames are encoded as the 8-bit red, green, blue (RGB) channels of a 24-bit image to interface with the SLM driver electronics. The SLM displays each bit plane sequentially, with ones and zeros mapping to opposing phase modulations at each pixel. The reconstructions were captured using a Canon 550D camera with an EFS 18-55 mm lens. To ensure fair comparison, the camera was set to the same manual setting when comparing each pair of replay fields before and after DPD. It takes 24/1440 = 1/60s to display all 24 sub-frames on a 1440Hz SLM, so the camera shutter speed was set to 1/30s to capture all frames twice. The images captured are raw format in RGB colour, which are subsequently converted to grey-scale (using the rgb2gray function in Matlab[72]) when calculating NMSE 52 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm against the target images, using the equation NMSE = 1 n∗∑(x−x̂)2 ∑(x)2 , where x is the target, x̂ is the measured output and n is the dimension of x and x̂. 3.3 Determining the DPD curve Fig. 3.3 Determining the DPD curve. (a) Input linear grey-scale ramp. (b) Corresponding CGH of (a) with 24-subframe binary phase encoding. (c) Holographic projection replay field of (b). (d) Plot of non-linearity measurement and according pre-distortion curve. 3.3 Determining the DPD curve 53 To determine the DPD curve of the holographic projection system, the non-linearity needs to be measured first. The hologram in Fig. 3.3(b) was first generated using OSPR algorithm for the linear grey-scale ramp of pixel value increasing linearly from 0 to 255, as shown in Fig. 3.3(a), along with a single pixel white (255) strip at the left end as a fiducial marker to demonstrate the beginning of the grey-scale region [49]. The projection output of the linear grey-scale ramp was then captured and cropped as shown in Fig. 3.3(c), from which the non-linearity curve was determined, by averaging each column of pixels in the image and discarding the fiducial marker, forming the blue line in Fig. 3.3(d). A third-order polynomial fit was applied, generating a smoothed non-linearity curve (yellow line in Fig. 3.3(d)). There exhibits a high degree of non-linearity. By taking the mean of the square of the error between the measured output (blue line) and the linear reference (green dashed line), the normalized mean squared error (NMSE) of the measured output was calculated to be 0.0858. To correct for the non-linearity, the DPD curve (red line) was formed by inverting the smoothed non-linearity curve (yellow line) in Fig. 3.3(d). 54 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm Fig. 3.4 Validation of DPD curve on the grey-scale ramp. (a) Pre-distorted ramp. (b) Corresponding CGH of (a) with 24-subframe binary phase encoding. (c) Holographic projection replay field of (b). (d) Non-linearity measurement after DPD. Subsequently, the DPD curve (red line in Fig. 3.3(d)) was used to adjust the grey-scale ramp, achieving the pre-distorted grey-scale ramp as shown in Fig. 3.4(a). The according projection output was then captured as shown in Fig. 3.4(c). By using the same method of averaging columns of pixels, the measured output was plotted in Fig. 3.4(d). It can be seen that the 3.4 Applying the DPD Curve 55 corrected non-linearity was much closer to linear comparing to the original non-linearity, and the NMSE was calculated to be 0.0039. NMSE Percentage Before DPD 0.0858 100% After DPD 0.0039 4.55% Table 3.1 Non-linearity results before and after DPD Hence, as demonstrated in Table 3.1, DPD achieved a 95.45% reduction in NMSE, which was a significant improvement in non-linearity, therefore the DPD curve measured is validated. 3.4 Applying the DPD Curve Fig. 3.5 Application of DPD on the 10-step strips. (a) 10 strips with equal step of pixel value. (b) CGH of (a). (c) Holographic projection replay field of (b). (d) After DPD of (a). (e) CGH of (d). (f) Holographic projection replay field of (e). To qualitatively demonstrate the effectiveness of our approach, we project a simple test pattern of a graduated ramp test pattern consisted of 10-step strips in Fig. 3.5(a), which is commonly employed in gamma-correction calibration of many display systems. As shown 56 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm in the projection replay field captured in Fig. 3.5(c), before DPD, the right few strips are barely distinguishable. In comparison, after carrying out DPD, it can be seen that each pair of adjacent strips in Fig. 3.5(f) are much more distinguishable, qualitatively showing the effectiveness of the DPD method. Fig. 3.6 Application of DPD on two sample real-word images. (a) Sample image 1: City Scene [73]. (b) Sample image 2: Horse. (c) Sample image 1 after DPD. (d) Sample image 2 after DPD. Then the DPD curve was applied to the two sample images as shown in Fig. 3.6 (a) and (b), producing pre-distorted images in Fig. 3.6 (c) and (d). Holograms were generated for each image using the OSPR algorithm and loaded onto the SLM respectively. The according replay fields were captured as shown in Fig. 3.7. 3.4 Applying the DPD Curve 57 Fig. 3.7 Projection output of the two sample images before and after DPD. (a) Replay field of Sample image 1 before DPD (NMSE=0.06139). (b) Replay field of Sample image 2 before DPD (NMSE=0.04309). (c) Replay field of Sample image 1 after DPD (NMSE=0.04920). (d) Replay field of Sample image 2 after DPD (NMSE=0.03635). The replay fields of the holographic projection of original images are shown in Fig. 3.7 (a) and (b), and the replay fields of the holographic projection of images after DPD are shown in Fig. 3.7 (c) and (d). As shown in Fig. 3.7(a), it can be seen that, before DPD, the edges between the buildings and the sky were quite ambiguous, with most detail of the sky being lost. In comparison, after DPD, the replay field in Fig. 3.7(c) provided not only sharper edges between buildings and 58 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm the sky, but also more detail of clouds in the sky. The NMSE of the replay field for sample image 1 decreased from 0.06139 to 0.04920, which was a 19.86% reduction. In Fig. 3.7(b), before DPD, the horse was difficult to distinguish from the background, especially around the horse’s back area. But after DPD, as shown in Fig. 3.7(d), contrast has been significantly boosted and the fine detail around this part of the horse is more evident. The NMSE of the replay field for sample image 2 decreased from 0.04309 to 0.03635, which was a 15.64% reduction. Sample image 1 NMSE Percentage Before DPD 0.06139 100% After DPD 0.04920 80.15% Sample image 2 NMSE Percentage Before DPD 0.04309 100% After DPD 0.03635 84.36% Table 3.2 DPD results for sample images Hence, as summarised in Table 3.2, DPD achieved a 19.86% reduction in NMSE for sample image 1 and a 15.64% reduction in NMSE for sample image 2, quantitatively proving the effectiveness of DPD method for CGH of real-world test images using OSPR algorithm. Lastly, as the DPD is a one-to-one mapping, the computation time is negligible. In practice, the computational overhead is too small to be measured against randomness between subse- quent runs. DPD can also be further accelerated in hardware using a hardware LUT, so that the DPD can be carried out instantly. This approach is widely adopted in gamma correction for displays. 3.5 Summary The non-linearity between the target image and the reconstructed image was measured for the overall holographic projection system by projecting a linear grey-scale ramp. The DPD curve was then applied to the grey-scale ramp, successfully reducing the NMSE by 95.45%. To evaluate its effectiveness on real world images, the DPD method was applied on two sample images. It was observed that more details were shown in the replay field after DPD, and the NMSE’s of the two example images were reduced by 19.86% and 15.64%, respectively. 3.5 Summary 59 As the DPD is a one-to-one mapping, the extra computation required is negligible. Hence, we have demonstrated the effectiveness of the proposed DPD-OSPR method in improving reconstruction quality over the original OSPR algorithm while still maintaining its capability for real-time holography. The proposed DPD method needs to be applied only once for every holographic projector; however, compared to the original OSPR method, the proposed DPD-OSPR method requires an additional device of a camera, whose accuracy and linearity can adversely affect the effectiveness of the DPD process. Chapter 4 L-BFGS Optimisation of 2D and 3D CGH Note: The work in this chapter have been published in Ref. [2, 4, 10] As previously introduced in Chapter 2, currently available SLMs can only modulate either phase or amplitude, so algorithms are needed to compute amplitude-only or phase-only holograms. Phase-only holograms are usually preferred among the two due to their higher energy efficiency, leading to the emergence of the classical phase retrieval algorithms re- viewed in Section 2.3. With the developments in modern numerical optimisation methods and computational power, advances in CGH algorithms can be made. This chapter therefore implements numerical optimisation methods for CGHs and proposes a novel target image phase optimisation (TIPO) algorithm. It also explores the use of sequential slicing (SS) techniques in optimisation algorithms. 4.1 Numerical Optimisation Methods 4.1.1 Optimisation framework Numerical optimisation methods aim to find an optimal solution which minimise an objective function numerically. They begin with an initial guess of the optimal solution (x0) and then, after iterations, generate a sequence of gradually improved estimates until they reach a solution [74]. If we have x as the vector of variables, and denote f (x) as the objective function, which is a function of x we want to minimise, any unconstrained optimisation 62 L-BFGS Optimisation of 2D and 3D CGH problem can be written as: minimise x∈Rn f (x) (4.1) Numerical optimisation then calculate the optimal solution x∗ iteratively, the iteration is given by: xk+1 = xk +αkpk (4.2) where the positive scalar αk is called step length, or sometimes may be referred as ‘learning rate’ in some contexts, and the vector pk is the search direction, which usually takes the form of pk =−B−1 k ∇ fk (4.3) where Bk is a nonsingular matrix that varies for different optimisation methods, and ∇ fk is the gradient, which, if unable to evaluate directly, can be approximated by: ∇ fk ≈ fk+1− fk xk+1−xk where fk denotes f (xk) (4.4) The strategy used to determine pk distinguishes one algorithm from another. Most methods make use of the values of f , ∇ f and ∇2 f , and some methods even make use of the accu- mulated historical values of those derivatives, which are further discussed in Section 4.1.2 - Section 4.1.5. 4.1.2 Gradient Descent As reviewed in Ref. [75], Gradient descent (GD) is a numerical optimisation method that is frequently seen in the literature for CGH. GD is a first-order optimisation method, it finds a local minimum by following the negative of the gradient (i.e. the steepest descent direction). The Bk (in Eq. (4.3)) for gradient descent simply takes the value of I, which is the identity matrix. Then the search direction becomes: pk =−∇ fk (4.5) The steepest descent method is very intuitive: among all possible directions to move away from xk, the steepest gradient direction is the one which f decreases most rapidly. The advantage of this method is that it requires few computations and memory resources, because 4.1 Numerical Optimisation Methods 63 it only requires a computation of the first derivative, and it does not require any accumulation of historical gradients. However, it is a greedy method that only considers the current iteration without any global consideration, so it can be extremely slow on complicated problems. [74] To work around the disadvantage, a few variants have emerged, such as AdaGrad [76], RMSProp [77] and Adam [78] which combines the advantages of AdaGrad and RMSProp. The Adam method is an iconic variant of GD, often referred to as ‘gradient descent with momentum’, as the name Adam is derived from ‘adaptive moment estimation’. Adam algorithm is based on adaptive estimates of lower-order moments [78]. It computes individual adaptive learning rates for different parameters from estimates of first and second moments of the gradients [78]. 4.1.3 Newton’s Method Newton’s method is a second-order optimisation method. Its search direction is derived from the second-order Taylor series approximation to f (xk +p), which is f (xk +p)≈ fk +pT ∇ fk + 1 2 pT ∇ 2 fkp def = mk(p) (4.6) The Newton direction can then be obtained by finding the vector p that minimises mk(p). By setting the derivative of mk(p) to zero, p can be obtained as: pk =−∇ 2 f−1 k ∇ fk (4.7) By comparing Eq. (4.7) to Eq. (4.3), it can be seen that the Newton’s method has a Bk of ∇2 fk. Unlike the gradient descent method, there is a "natural" step length of 1 associated with the Newton direction, so αk = 1 by default and is only adjusted when it does not produce a satisfactory reduction in the value of f . The Newton direction is reliable when the difference between the true function f (xk +p) and its quadratic model mk(p) is not too large. Methods that use the Newton direction have a fast rate of local convergence, typically quadratic. After a neighbourhood of the solution is reached, convergence to high accuracy often occurs in just a few iterations. The main drawback of the Newton direction is the need for the Hessian ∇2 fk. Explicit computation of this matrix of second derivatives can sometimes be a cumbersome, error-prone, and expensive process. [74] 64 L-BFGS Optimisation of 2D and 3D CGH 4.1.4 Quasi-Newton Method Quasi-Newton method provides an attractive alternative to Newton’s method, in that it does not require computation of the Hessian and yet still attains a linear rate of convergence. In place of the true Hessian ∇2 fk, they use an approximation Hk def = B−1 k , which is updated after each step to take account of the additional knowledge gained during the step. The updates make use of the fact that changes in the gradient provide information about the second derivative of f along the search direction. The most popular quasi-Newton algorithm is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, named for its discoverers Broyden, Fletcher, Goldfarb, and Shanno. [74] The process of the BFGS method is shown below: denote { Hk = B−1 k pk =−Hk∇ fk (4.8) Initiate H0← yT k sk yT k yk I (4.9) update Hk+1 = (I−ρkskyT k )Hk(I−ρkyksT k )+ρksksT k (4.10) where  sk = xk+1−xk yk = ∇ fk+1−∇ fk ρk = 1 yT k sk (4.11) The algorithm is robust, and its rate of convergence is linear, which is fast enough for most practical purposes. Even though Newton’s method converges more rapidly (that is, quadratically), its cost per iteration usually is higher, because of its need for second derivatives and solution of a linear system. The drawback is that, it is not directly applicable to large optimisation problems because Hk’s are usually dense, requiring large storage and computational requirements. [74] 4.1.5 Large Scale Quasi-Newton Method: Limited Memory BFGS (L- BFGS) L-BFGS algorithm [79] modifies the technique described in Section 4.1.4 to obtain Hessian approximations that can be stored compactly in just a few vectors of length n, where n is the number of unknowns in the problem. The main idea of this method is to use the curvature information from only the most recent iterations to construct the Hessian approximation. 4.1 Numerical Optimisation Methods 65 Curvature information from earlier iterations, which is less likely to be relevant to the actual behaviour of the Hessian at the current iteration, is discarded in the interest of saving storage. [74] Denoting Vk = I−ρkyksT k , Eq. (4.10) can be written as: Hk+1 = VT k HkVk +ρksksT k (4.12) The inverse Hessian approximation Hk will generally be dense, so that the cost of storing and manipulating it is prohibitive when the number of variables is large. To circumvent this problem, we store a modified version of Hk implicitly, by storing a certain number (say, m) of the vector pairs {si,yi} used in the Eq. (4.10) and Eq. (4.11). The product Hk∇ fk can be obtained by performing a sequence of inner products and vector summations involving ∇ fk and the pairs {si,yi}. After the new iterate is computed, the oldest vector pair in the set of pairs {si,yi} is replaced by the new pair {sk,yk} obtained from the current step (Eq. (4.11)). In this way, the set of vector pairs includes curvature information from the m most recent iterations. Practical experience has shown that modest values of m (between 3 and 20, say) often produce satisfactory results. We now describe the updating process in a little more detail. At iteration k, the current iterate is xk and the set of vector pairs is given by {si,yi} for i = k−m, . . . ,k−1. We first choose some initial Hessian approximation H 0 k (in contrast to the standard BFGS iteration, this initial approximation is allowed to vary from iteration to iteration) and find by repeated application of Eq. (4.10) that the L-BFGS approximation Hk satisfies the following formula: [74] Hk = (VT k−1 · · ·VT k−m)H 0 k (Vk−m · · ·Vk−1) +ρk−m(VT k−1 · · ·VT k−m+1)sk−msT k−m(Vk−m+1 · · ·Vk−1) +ρk−m+1(VT k−1 · · ·VT k−m+2)sk−m+1sT k−m+1(Vk−m+2 · · ·Vk−1) + · · · +ρk−1sk−1sT k−1 (4.13) From this expression we can derive a recursive procedure (Algorithm 8) to compute the product Hk∇ fk efficiently. 66 L-BFGS Optimisation of 2D and 3D CGH Algorithm 8 L-BFGS two-loop recursion [74] Input: Iteration number k, Gradient at current iteration ∇ fk, Historic gradient window length m, Initial hessian matrix H 0 k , ρ , s and y are defined in Eq. (4.9) Output: Search direction −r q← ∇ fk for i = k−1,k−2, ...,k−m do αi← ρisT i q q← q−αiyi end for r←H 0 k q for i = k−m,k−m+1, ...,k−1 do β ← ρiyT i r r← r+ si(αi−β ) end for Step with pk←−Hk∇ fk =−r Apart from being inexpensive, L-BFGS has the advantage that the multiplication by the initial matrix H 0 k is isolated from the rest of the computations, allowing this matrix to be chosen freely and to vary between iterations. A method for choosing H 0 k that has proved effective in practice is to use the same as BFGS as stated in Eq. (4.9). [74] 4.2 Phase-only Hologram Optimisation To implement the optimisation methods listed in Section 4.1 on the CGH, firstly the optimi- sation framework needs to be adapted. The objective of CGH is to find the phase hologram (H) that has the optimal reconstruction (R) matching the target image (T), which can be formulated as minimising the difference between R and T by varying H, leading to the mathematical expression below: argmin H Loss(T,R) (4.14) where ‘Loss’ denotes an error function quantifying the difference between T and R, and ‘arg’ returns the argument (H) instead of the error value. For the Loss function, the classic error 4.2 Phase-only Hologram Optimisation 67 function mean-squared error (MSE) [57] is selected, with its definition as shown below: MSE(T,R) = 1 X×Y X ∑ x=1 Y ∑ y=1 (Tx,y−Rx,y) 2 (4.15) Compute Loss (e.g. MSE(T, R)) Update H Compute Reconstruction (R) Step optimiser (e.g. L-BFGS) Yes No Terminate? Load target image (T) Initialise random phase hologram (H) Output H Fig. 4.1 Flowchart of the optimisation process The optimisation process for generating a phase hologram to reconstruct a target image is depicted in Fig. 4.1. The process begins by loading the target image (T) and initializing a random phase hologram (H). Using this initial phase hologram, a reconstruction (R) of the target image is computed via either the Fraunhofer or Fresnel propagation equation introduced in Section 2.2.2 based on the distance needed. The difference between the target image (T) and the reconstruction (R) is quantified by computing the loss, such as the MSE in Eq. (4.15). An optimiser (such as the GD, Adam or L-BFGS algorithms mentioned in Section 4.1) then updates the phase hologram based on the computed loss. This iterative process of reconstruction, loss calculation, and hologram update continues until a termination condition is met, usually a fixed total number of iterations, or when optimisation has reached convergence. The final optimised phase hologram (H) is then outputted. 68 L-BFGS Optimisation of 2D and 3D CGH iteration number Fig. 4.2 Convergence plot for comparison between the GD, Adam and L-BFGS optimisations The optimisation flowchart in Fig. 4.1 is run on the example target image in Fig. 2.11 for a total of 50 iterations, using the GD, Adam and L-BFGS algorithms. The NMSE (defined in Eq. (2.27)) is plotted against the iteration number for each optimiser, in Fig. 4.2. Fig. 4.2 shows that, the proposed method using L-BFGS algorithm (the orange line in Fig. 4.2) stagnates for the first two iterations as it needs to estimate the Hessian as explained in Section 4.1.5. After the brief stagnation, the L-BFGS converges quickly, surpassing the existing methods in the literature using GD and Adam optimiser, corresponding to the blue line and the green line in Fig. 4.2 respectively. The final reconstructions are shown aside the NMSE plots, from which it can be seen that the L-BFGS algorithm produces a reconstruction of comparable quality to the Adam algorithm, and both of them are much better than the GD algorithm. Such observation matches the final NMSE values. The NMSE value of the L-BFGS algorithm (the orange line in Fig. 4.2) stays lower than the Adam optimiser (the blue line in Fig. 4.2) after the 3rd iteration, and reaches convergence earlier at the 30th iteration. In summary, the L-BFGS algorithm was proven to be able to optimise a phase-only hologram for a target image, and is much better than the GD optimiser and converges quicker than the Adam optimisation algorithm. However, the L-BFGS algorithm is not without its drawbacks, as it requires more memory and computational resources than the GD and Adam algorithms. 4.2 Phase-only Hologram Optimisation 69 On the laptop of model Mechrevo 16 Ultra, which has a CPU of model Intel Ultra 9 275HX and a GPU of model RTX5070Ti, using python with PyTorch and CUDA, the example run of the GD algorithm takes up 1311MB of memory while the L-BFGS algorithm takes up 1867MB of memory, which is a 556MB ( 42.4%) increase in memory usage. R₁ R2 R3 R4 R5 R6 R7 R8 R9 R10 R1₁ R12 R13 R14 R15 R46 R47 R48 R49 R50 … Fig. 4.3 Reconstructions at each iteration of the L-BFGS optimisation The reconstruction (R) at each iteration of the L-BFGS optimisation is listed in Fig. 4.3, with iterations 16 to 45 omitted to save space. The list of reconstructions demonstrate visually how the optimisation converges to the resulting reconstruction (R50) which has a very good quality. 70 L-BFGS Optimisation of 2D and 3D CGH 4.3 Target Image Phase Optimisation (TIPO) This section proposes a novel method of optimising the phase of the target image instead of optimising the phase of the hologram as previously introduced in Section 4.2, under the same ‘phase only’ constraint of the hologram. The objective is to find an optimal phase profile to be attached to the target image, so that its inverse propagation to the hologram plane produces a hologram whose phase has a reconstruction that best matches the target image. As it optimises the phase of the target image instead of the traditional way optimising the phase of the hologram, this method is named as Target Image Phase Optimisation (TIPO). A flowchart is drawn in Fig. 4.4 to clarify the detailed operation of the proposed TIPO algorithm. Compute Loss (e.g. MSE(T, R)) Compute Reconstruction R = |P (ejH)| Update Φ (e.g. using L-BFGS) Yes No Terminate? Load target (T) Initialise phase (Φ) Output H Compute target field E = T·ejΦ Compute hologram plane A = P -1(E) Compute phase hologram H = ∠A Fig. 4.4 TIPO flowchart The flowchart in Fig. 4.4 outlines the TIPO algorithm process for generating a phase-only hologram (H) to reconstruct a target image (T). The procedure starts with loading the target image and initializing the phase (Φ). The complex target field (E) is then computed from the target image and the target image phase (E = T · e jΦ). Next, the hologram aperture (A) is computed by applying the inverse propagation (P−1) to the target field, where P is chosen to be the Fraunhofer diffraction equations in Section 2.2.2. The phase-only hologram (H) is then derived from the angle of the complex hologram aperture (H = ∠A). Subsequently, a reconstruction (R) is computed using the forward propagation (R = |P(e jH)|). The loss, 4.3 Target Image Phase Optimisation (TIPO) 71 such as the Mean Squared Error (MSE), is calculated between the target image (T) and the reconstruction (R). An optimiser (e.g. SGD or L-BFGS) then updates the phase (Φ) based on the computed loss. This iterative process of computing the target field, hologram plane, phase hologram, reconstruction, and loss calculation continues, followed by phase updates, until a termination condition is met, which is usually set for a fixed number of iterations. Upon convergence, the final optimised phase hologram (H) is produced. 72 L-BFGS Optimisation of 2D and 3D CGH Φi Hi Rii 1 50 100 30 … … … Fig. 4.5 TIPO iterations on the mandrill target (i is the iteration number, Φi, Hi, Ri are the phase of the target field, the phase hologram, and the reconstruction amplitude at iteration i) The results for an example run of the TIPO algorithm on the mandrill target (in Fig. 2.11) with the total number of iterations set to 100 are shown in Fig. 4.5 (only iterations number 1, 30, 50 and 100 are shown for space saving). As the first hologram is computed by the inverse propagation (P−1) of the target image with an random phase, the reconstruction (R1) is 4.3 Target Image Phase Optimisation (TIPO) 73 already showing the target image. Then as the iterations continue, the reconstruction quality gets better, proving the effectiveness of the proposed TIPO algorithm. 0 20 40 60 80 100 iteration number 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NM SE 1e 6 Phase hologram optimisation using SGD Phase hologram optimisation using L-BFGS Target image phase optimisation using SGD Target image phase optimisation using L-BFGS Fig. 4.6 TIPO convergence plot To quantitatively analyse the results, the convergence of the TIPO algorithm is plotted in Fig. 4.6, where the NMSE between the reconstruction and the target image are plotted against the number of iterations. Both SGD and L-BFGS optimisers are run for the TIPO algorithm (corresponding to the green and red line in Fig. 4.6 accordingly), and are compared against the regular phase hologram optimisation algorithm in Section 4.2 (corresponding to the dotted blue line and the dotted orange line in Fig. 4.6 respectively). The TIPO algorithms using both SGD and L-BFGS optimisers start with a lower NMSE of the reconstruction as their holograms in the first iteration are extracted from the inverse propagation of the target image, instead of pure random holograms as done in the regular phase hologram optimisation. However, the regular phase hologram optimisation using L-BFGS optimiser quickly surpassed the TIPO algorithm within 10 iterations, and reached the lowest reconstruction error. For the SGD optimiser, the TIPO method (green line) shows a significant improvement over 74 L-BFGS Optimisation of 2D and 3D CGH the regular phase hologram optimisation method (dotted blue line). The two TIPO methods (in solid lines) lie between the two regular phase hologram optimisation method (in dotted lines), showing that the conventional L-BFGS phase hologram optimisation method has the fastest convergence, while the two TIPO methods are better than the SGD phase hologram optimisation method, proving the proposed TIPO method to be an effective method for phase hologram generation. 4.4 Multi-Depth Phase-Only Hologram Optimisation In addition to the 3D CGH methods reviewed in Section 2.3.7, further search in the literature has found several recent advances in 3D CGH. For example, Chen et al. [80] proposed an angular-tiling layer-based approach that treats a 3D object as stacked depth slices and uses a graphics rendering pipeline to incorporate depth cues into the hologram. By combining these techniques with parallel GPU computation, they reported generating a single-view XGA-resolution hologram in around 176 ms, enabling realistic 3D image reconstruction. In a follow-up study, Chen and Chu [81] introduced several refinements to the layer-based method: a more efficient coding scheme, a multi-layer depth-fusion technique, and a fractional- depth encoding. These innovations yielded over a 4× reduction in computation time and improved depth fidelity (including accommodation cues), allowing real-time interactive holography (e.g. for head-mounted displays). Jia et al. [82] similarly exploited layer sparsity by developing a two-step algorithm that computes only the nonzero (visible) regions of each depth layer, together with a “sub-sparse” 2D FFT that performs two 1D FFTs while skipping rows/columns of zeros. This significantly cuts the active workload, resulting in roughly 5−10× speedups on 3D scenes over conventional layer-based hologram methods. In summary, these works demonstrate that by combining layer-based scene slicing with optimized FFT strategies and parallelism, one can dramatically accelerate CGH computation while preserving or enhancing image realism. There are also some other work in the literature that compute CGH for 3D targets using numerical optimisation methods [23, 24, 26, 25], but speed and quality are still the major challenges. They either evaluate the error of reconstructions against the entire 3D target, which is time-consuming, or evaluate the hologram for each plane and then sum the holograms, which introduces quality degradation. This section reviews the traditional multi-depth hologram optimisation methods, and proposes the novel use of sequential slicing (SS) techniques during the optimisation process, which only evaluates the loss for a single slice at each iteration. The proposed technique aims for a 4.4 Multi-Depth Phase-Only Hologram Optimisation 75 quicker hologram generation with proper overall quality and low quality imbalance across the multiple depths. 4.4.1 Methods The optimiser used for CGH in this section is the L-BFGS optimiser as introduced in Section 4.1.5, with the GD optimiser in Section 4.1.2 also implemented as a reference. The phase-only constraint of CGH is applied by fixing a constant amplitude of the hologram, while keeping its phase varying and being the argument of optimisation (x in Eq. (4.1)). Fig. 4.7 Loss between the multi-depth targets (T1 to Tn) and the reconstructions (R1 to Rn) of hologram H As shown in Fig. 4.7, the multi-depth target is set up as a collection of n slices (T1 to Tn), each slice Ti is at a distance zi to the hologram plane. And for the hologram H, its reconstruction at each distance zi is computed using Fresnel diffraction formula in Eq. (2.22), which is labelled as the propagation function P(H,zi). 76 L-BFGS Optimisation of 2D and 3D CGH Loss Functions To formulate an objective function f (x) in Eq. (4.1), we need to quantify the difference between each target slice (Ti) and the respective reconstruction (Ri) numerically, which is denoted as Loss(Ti,Ri) in Fig. 4.7. In addition to the mean squared error (MSE) [57] previously used in Section 4.2, the cross entropy (CE) [83] and relative entropy (RE) [84] are also implemented. To adapt the loss functions for two-dimensional (2D) target image Ti and reconstructed image Ri of dimension X×Y , the loss functions are adapted as shown in Eq. (4.16) to Eq. (4.18). MSE(Ti,Ri) = 1 X×Y X ∑ x=1 Y ∑ y=1 (Ti;x,y−Ri;x,y) 2 (4.16) CE(Ti,Ri) =− X ∑ x=1 Y ∑ y=1 Ti;x,y log(Ri;x,y) (4.17) RE(Ti,Ri) =− X ∑ x=1 Y ∑ y=1 Ti;x,y log ( Ri;x,y Ti;x,y ) (4.18) MSE is a traditional metric averaging the squared error between the target and observed values. CE, adapted as shown in Eq. (4.17), is usually used in classification problems, such as language modelling [85]. RE, also called Kullback-Leibler divergence (usually denoted as DKL(P||Q), but is denoted as RE here for uniformity), is adapted as shown in Eq. (4.18). RE is usually used to measure how much a probability distribution P is different from another probability distribution Q. Both CE and RE are usually computed between the true probabilistic distribution and the predicted probabilistic distribution. While the images are not probability distributions, the pixel values can be normalized to decimal numbers in the range from 0 to 1 so that CE and RE can be applied. Sum-of-Loss (SoL) technique The conventional technique to compute 3D CGH is to sum the losses for each slice at each iteration during optimisation, which is called the Sum-of-Loss (SoL) method here. At every iteration, it computes the full 3D reconstructions (R1, · · · ,Rn) of the hologram H at every distance zi, and then sum the losses between each pair of reconstruction Ri and target image Ti. The mathematical expression of the SoL technique’s optimisation objective is shown 4.4 Multi-Depth Phase-Only Hologram Optimisation 77 below: SoL: argmin H n ∑ i=1 Loss(Ti,Ri) (4.19) The SoL method requires a total of n Fourier Transforms to fully evaluate the hologram at each step, making it computationally heavy. Sum-of-Hologram (SoH) technique Another technique for 3D CGH is to compute the complex sum of the sub-holograms Hi generated for each target slice Ti to form a total hologram based on the principle of super- position, which is called the Sum-of-Hologram (SoH) method here. The SoH technique’s mathematical expression is shown in Eq. (4.20) below: SoH: ∠ n ∑ i=1 e j·argmin Hi Loss(Ti,Ri) (4.20) The sub-hologram Hi for each slice number i is computed independently with the optimisation objective of min Hi Loss(Ti,Ri), the arguments of which are then summed in a complex manner, as each phase hologram is the angle, the exponential operator is needed to turn angles into complex numbers. Lastly, the angle of the complex sum is taken so that it can meet the ‘phase-only’ constraint. If using a fixed number of iterations per sub-hologram, the total computation scales up linearly with the number of slices n. The SoH method’s advantage is its ease of implementation, that any existing single slice CGH algorithm can be quickly converted to multi-depth 3D CGH algorithm. Its major disadvantage for phase-only hologram generation is that, the final summing of sub-holograms will result in a non-uniform amplitude hologram, and taking the phase of which will result in discarding the amplitude information of the summed hologram, leading to deprecations in reconstructions quality. And also, the SoH method suffers from the defocusing effect between the each slice to another, causing additional noise. 78 L-BFGS Optimisation of 2D and 3D CGH Sequential Slicing (SS) technique Fig. 4.8 Optimisation of CGH with sequential slicing (SS) flowchart This section proposes the novel CGH optimisation with sequential slicing (SS) technique, as shown in the flowchart in Fig. 4.8, that only computes the loss for a single slice at each iteration (between a reconstruction Ri at a single distance zi and its according target slice Ti), where i sweeps through the multi-layer 3D target sequentially when the algorithm iterates. When the final slice is reached (i = n), it goes back to the first slice (i← 1). The proposed method only needs to carry out one Fourier Transform at each iteration, and the number of iterations does not need to scale up with n. So it is expected to be much quicker than SoL and SoH techniques while producing a proper resulting hologram. 4.4 Multi-Depth Phase-Only Hologram Optimisation 79 4.4.2 Results CGH for a 4-slice target consisted of letters ‘A, B, C, D’ Fig. 4.9 Layout of the 4-slice target (z1 = 1 cm, z2 = 2 cm, z3 = 3 cm, z4 = 4 cm) The first example 3D target used consisted of 4 slices made from letters ‘A’, ‘B’, ‘C’, ‘D’, each having 512×512 pixels. The positions of the four slices range from 1 cm to 4 cm with 1 cm gap between each other (i.e. zi = i cm). The overall layout is shown in Fig. 4.9. As there are two optimisers (GD and L-BFGS), three techniques (SoL, SoH and SS) and three loss functions (MSE, CE, RE) in consideration, they can form a total of 18 combinations. In order to control the number of variables, all 18 combinations were set to start from the same initial random hologram and run for the same amount of 100 iterations for the optimisation of each hologram, on the same laptop of model ASUS ROG Zephyrus M16, which has a CPU of model i7-11800H and a GPU of model RTX3060. For the L-BFGS algorithm, the gradient history of size 10 (m = 10 in Algorithm 8) was used for all techniques and loss functions. And to ensure a sensible comparison, although three different loss functions are used for optimisation of CGH, the metric used to assess the final quality of multi-depth reconstructions of the hologram is the normalized mean squared error (NMSE). As there are a total of 4 slices in this example, the final NMSE of each are computed separately and the total optimisation run time is recorded. The final results are gathered in Fig. 4.10. As each slice has a different final error, their mean and standard deviation (SD) are also computed for 80 L-BFGS Optimisation of 2D and 3D CGH investigations. Three columns are colour coded where green indicates better a result while red indicates a worse result. Fig. 4.10 Final NMSE and run time comparison across the three techniques Comparing the mean of final NMSE and the run time of the proposed sequential slicing (SS) technique to those of the sum-of-loss (SoL) and sum-of-hologram (SoH) techniques in Fig. 4.10, it can be concluded that, for all combinations of optimisers and loss functions attempted, the proposed SS technique runs much quicker than the existing SoL and SoH techniques, while still providing a proper result, sitting between the SoL and SoH techniques. The SS technique is both quicker and has better reconstruction quality than the SoH technique, demonstrating an absolute advantage. Meanwhile, the runs with CE as loss function are much slower and has not demonstrated any advantage in the final NMSE, demonstrating an absolute disadvantage, so the results of runs using CE loss or SoH method will not be shown in the per-iteration plots later on in this section. 4.4 Multi-Depth Phase-Only Hologram Optimisation 81 (a) (b) Fig. 4.11 NMSE of each slice plotted against the iteration number for the optimisation based SS technique implemented with MSE and RE loss functions and (a) GD optimiser, (b) L-BFGS optimiser. 82 L-BFGS Optimisation of 2D and 3D CGH (a) (b) Fig. 4.12 NMSE of each slice plotted against the iteration number for the (a) GS with SS algorithm, (b) DCGS algorithm. 4.4 Multi-Depth Phase-Only Hologram Optimisation 83 For comparison among the runs using the sequential slicing (SS) techniques, the NMSE for each slice are plotted for the GD and LBFGS optimisers with MSE and RE loss functions in Fig. 4.11, and the GS-based sequential GS and DCGS [63] methods (mentioned in Section 2.3.7) are also plotted in Fig. 4.12 for reference. Looking at the plots of NMSE of each slice against iteration numbers (in Fig. 4.11 and Fig. 4.12), apart from the L-BFGS algorithm (Fig. 4.11b), all the other algorithms (in Fig. 4.11a Fig. 4.12a Fig. 4.12b) are showing a staircase-like trend, where a decrease in error on one slice results in an increase in error on all the other slices, and the final NMSE of each slice distinguishes a lot from another. The sequential GS algorithm suffers the most from the quality imbalance between each slice, and the sequential GD algorithm follows. The DCGS algorithm benefits from its modification of the inclusion of a weighting factor consisting historical amplitude, therefore managed to converge. The average and maximum difference of the NMSE values across the 4 slices are then computed and summarised in Fig. 4.13. 84 L-BFGS Optimisation of 2D and 3D CGH (a) (b) Fig. 4.13 (a) Average NMSE among all slices, (b) Maximum difference of NMSE across all slices, plotted against the iteration number for the 6 sequential slicing runs using different techniques 4.4 Multi-Depth Phase-Only Hologram Optimisation 85 Fig. 4.13 plots the average and maximum difference values of NMSE across the 4 slices against the iteration number, for all the six runs (sequential GS, DCGS, GD optimiser with MSE as loss function, GD optimiser with RE as loss function, L-BFGS optimiser with MSE as loss function, and L-BFGS optimiser with RE as loss function). From the average NMSE plot (Fig. 4.13a), the proposed sequential L-BFGS method does not have the lowest average NMSE, but it has the lowest quality imbalance across the slices as shown in the maximum difference plot (Fig. 4.13b). The L-BFGS algorithm mainly benefits from its inclusion of curvature information during optimisation, so that the update of hologram H at each iteration takes into account not only the loss for the current slice, but also all the historical iterations up to the set history size (m in Algorithm 8). So for the L-BFGS algorithm, at each iteration, the NMSE of all slices behave in the same way, ensuring each slice to have the similar quality (i.e. lower difference of NMSE across slices). Hence, the proposed SS technique with L-BFGS optimiser is shown to have the lowest quality imbalance across all slices, although the average quality across the 4 slices is worth than the GS based algorithms. 86 L-BFGS Optimisation of 2D and 3D CGH Fig. 4.14 Comparison of final holograms and reconstructions on the 4-slice target consisted of letters ‘A’, ‘B’, ‘C’ and ‘D’ The final holograms and their reconstructions for L-BFGS algorithm with SoL, SoH and SS techniques are shown in Fig. 4.14, with GS with SS technique and DCGS also shown as reference. The reconstructed images confirm the SS technique having a quality between SoH and SoL method (for the same amount of iterations), and has a much better quality imbalance than sequential GS, which has a very clear reconstruction at the fourth slice (letter ‘D’) because the iteration stopped at the fourth slice but much worse reconstruction at other slices. Admittedly the proposed L-BFGS with SS method cannot surpass the GS-based DCGS algorithm yet, but among the optimisation-based methods, the proposed L-BFGS with SS 4.4 Multi-Depth Phase-Only Hologram Optimisation 87 combination has shown better reconstruction quality than the SoH method inFig. 4.14, and takes fewer time to run than the SoH and the SoL methods as previously shown in Fig. 4.10. CGH for a 4-slice target consisted of letters and real-life images Fig. 4.15 Layout of the non-binary 4-slice target To prove that the proposed method also works for non-binary valued target images, another example of a 4-slice 3D target is attempted, as shown in Fig. 4.15, where two of the slices are replaced by an image of the mandrill [55] and an image of the city scene [23] respectively. 88 L-BFGS Optimisation of 2D and 3D CGH Fig. 4.16 Comparison of final holograms and reconstructions for non-binary target As shown in the final holograms and reconstructions in Fig. 4.16, the proposed LBFGS with SS technique still managed to converge, with final reconstruction quality sitting between the SoL and the SoH method, and also having a good quality balance across all slices. As the SS technique is faster than both the SoL and the SoH method, its quality performance is impressive. 4.5 Summary 89 4.5 Summary This chapter began by covering the fundamentals of numerical optimisation, including the L- BFGS optimiser, and proceeded to introduce and perform the optimisation of phase hologram. The novel TIPO algorithm was then proposed, which optimises the phase of the target image instead of the phase of the hologram. The TIPO showed its effectiveness for CGH but did not exhibit any significant advantage over conventional phase hologram optimisation methods. Finally, a method using the L-BFGS optimiser in conjunction with the SS technique was proposed for multi-depth CGH. This method effectively suppressed quality imbalances across multi-depth slices, leveraging the second order nature of the L-BFGS optimiser, which implicitly records the historical gradients from other slices to determine the descent direction. For both GD and L-BFGS optimisation algorithms, the SS technique runs faster and produces better reconstruction quality compared to the simple SoH technique and is much quicker than the SoL technique. Therefore, the proposed L-BFGS optimisation with SS method has demonstrated strong capability for time-constrained optimisation of multi-depth 3D CGH. However, admittedly, the L-BFGS optimiser introduces more memory usage than the GD optimiser, which can lead to out-of-memory issues for memory resource limited applications, and the L-BFGS optimiser is an iterative method that is intrinsically slower than non-iterative CGH methods such as the double-phase method. However, all the CGH methods investigated in this chapter produce single-frame multi-level phase hologram, which cannot be used directly on the binary-phase SLM in the optical setup described in Section 3.2. Therefore, the next chapter will propose a multi-frame binary-phase holograms batched optimisation (MFHBO) method to optimise multi-frame binary-phase holograms which are the most suitable for the optical setup in Section 3.2. Chapter 5 Multi-Frame Binary-Phase Holograms Batched Optimisation Note: The work in this chapter has been published in Ref. [9] This chapter builds on the idea of time-averaging multiple hologram frames and proposes a novel technique called Multi-Frame Holograms Batched Optimisation (MFHBO). This technique employs the L-BFGS optimisation algorithm to simultaneously generate a batch of phase-only holograms, which results in an average reconstructed image of improved fidelity and rapid algorithmic convergence, in both the Fraunhofer and the Fresnel regions. 5.1 Introduction Chapter 4 implemented optimisation algorithms to generate single-frame holograms. This chapter proposes the novel MFHBO method which optimise a batch of multi-frame holograms for time multiplexing. Time multiplexing aims to improve a time-averaged response by displaying different hologram sub-frames at a high refresh rate [86]. This approach can exploit the finite response time of human vision, wherein human eyes average out the unwanted noise while preserving the desired signal. A few time-multiplexed multi-frame holograms generation methods have been explored in the literature, including the OSPR algorithm in Section 2.3.5 and the AD-OSPR algorithm in Section 2.3.6; however, both 92 Multi-Frame Binary-Phase Holograms Batched Optimisation OSPR and AD-OSPR are still subject to noise and defects in reconstruction quality. The objective of the proposed MFHBO algorithm is therefore to produce a batch of holograms with better reconstruction quality compared to the existing OSPR and AD-OSPR methods. 5.2 Methods The use of the L-BFGS algorithm for single-frame hologram optimisation has been introduced in Chapter 4. To extend it to multi-frame hologram generation, the optimisation variable becomes the set of n hologram sub-frames ({H1,H2, ...,Hn}), each with a resolution of X×Y pixels matching that of the target image. The objective function to minimise is the difference between the average reconstruction amplitude (Ravg = 1 n ∑ n i=1 Ri) and the target image (T), denoted as Loss(T,Ravg). Here, Ri denotes the reconstruction from the individual sub-frame Hi, for i ∈ [1,n]. Each Ri is computed from Hi using the Fresnel diffraction formula in Eq. (2.22). Since phase-only SLMs are used, the hologram aperture A in Eq. (2.22) is given by a uniform amplitude with phase H, i.e., A = e jH, with the exponential applied element-wise. The amplitude of the replay field E gives the reconstruction, i.e., R = |E|. Quantize [H1, H2, ..., Hn] Update [H1, H2, ..., Hn] using L-BFGS optimizer Compute each reconstruction sub-frame [R1, R2, ..., Rn] Terminate? No Yes Load target image T Generate initial holograms [H1, H2, ..., Hn] Output [H1, H2, ..., Hn] Compute average Ravg = (�Ri)/n i=1 n Compute Loss(T, Ravg) Fig. 5.1 MFHBO flowchart To illustrate the optimisation process, a flowchart is provided in Fig. 5.1. Initially, the target image T is loaded, and a set of n hologram sub-frames {H1,H2, ...,Hn} is randomly generated. In each iteration, the sub-frames Hi are quantised to the SLM’s bit-depth constraint. To improve the optimisation process, a sigmoid function is used to perform soft quantisation. This allows for a differentiable approximation to discrete phase levels, maintaining gradient flow and stabilising convergence. After quantisation, each Hi is propagated to generate its 5.3 Results 93 reconstruction Ri, and the average amplitude Ravg is computed. The loss Loss(T,Ravg) is then evaluated using the relative entropy metric [84], as defined in Eq. (4.18). Based on this loss, the L-BFGS optimiser computes the search direction and updates all hologram sub-frames accordingly. The optimisation process is repeated for a fixed number of iterations, after which it is terminated. Since fast SLM’s available in the lab are binary-phase devices, the quantisation step in the flowchart in Fig. 5.1 is carried out with bit-depth limit of 1, hence producing binary-phase holograms. However, the optimisation algorithm does not converge with a straight binary quantisation as integers are discrete, therefore a Sigmoid function [87] is used for a smoother and differentiable quantisation, as defined in Eq. (5.1). The output of the Sigmoid function is then scaled by π so that the binary phase levels are 0 and π . sigmoid(x) = 1 1+ e−x (5.1) Finally, when displaying the multi-frame holograms, each of the n frames generated are then rounded to binary phase values and displayed on the binary phase SLM sequentially. When the first round finished, the second round starts with the first frame again (i.e. after frame n, the next frame displayed is frame 1), and such infinite loop doesn’t stop until another set of holograms are uploaded. 5.3 Results 5.3.1 Simulation results To test the proposed MFHBO method, a target image T as shown in Fig. 5.2 was used. It was designed from the widely used mandrill image [55] in Fig. 2.11. A rotational symmetry was introduced to match the rotational symmetry of the far field projection from binary phase holograms, as explained in Section 2.2.3. It was then zero padded to a resolution of 1024px× 1024px and subsequently interpolated (using the ‘torch.nn. f unctional.interpolate’ function in the PyTorch module[88]) to a resolution of 1280px×1024px to match the resolution of the SLM in our lab. Note that the target image was zero padded to a square aspect ratio and then stretched to the non-square aspect ratio because more pixels in the horizontal axis only means higher sampling rate as part of the features of the FFT, the replay field is continuous 94 Multi-Frame Binary-Phase Holograms Batched Optimisation and is not pixelated and the simulated reconstruction of 1280px×1024px resolution is the sampled result, which will be illustrated visually in Fig. 5.4 later. H1 H2 H3 H24 .... R1 R2 R3 R24 .... Ravg T Loss Forward Propagation Averaging Ravg Update Fig. 5.2 An example iteration in the optimisation process To further explain the optimisation process described in Fig. 5.1, an example iteration with n = 24 is shown in Fig. 5.2. At each iteration, every hologram is quantised and propagated to the reconstruction plane, forming {R1,R2, ...,R24}. The average reconstruction amplitude Ravg is then compared against the target image T, using the loss function in Eq. (4.18). The holograms {H1,H2, ...,H24} are then updated according to the search direction calculated using the L-BFGS optimiser. After setting the optimisation to terminate when the number of iterations reach 1000, the same algorithm was run on the same target for different number of frames (n), the normalised mean squared error (NMSE) and the peak signal-to-noise ratio (PSNR) between the average reconstruction Ravg and the target image T were calculated at every iteration and plotted in Fig. 5.3a and Fig. 5.3b respectively. 5.3 Results 95 (a) NMSE v.s. Iteration Number (b) PSNR v.s. Iteration Number Fig. 5.3 Convergence of the MFHBO algorithm on the rotationally symmetrical Mandrill target The plots in Fig. 5.3 show that the proposed MFHBO method has achieved good convergence within 400 iterations, for the various number of frame settings n in {1,2,3,4,6,8,12,24}. The final NMSE values in Fig. 5.3a are difficult to distinguish in the plot, therefore it will be further compared in the bar chart in Fig. 5.4. The number of frames are chosen to be integer 96 Multi-Frame Binary-Phase Holograms Batched Optimisation factors of 24, which is determined by our experimental setup, further explained in the next subsection. n 200 iterations 400 iterations 600 iterations 800 iterations 1000 iterations 1 1.79 3.59 5.31 7.00 8.72 2 2.82 5.59 8.34 11.11 13.88 3 3.84 7.67 11.45 15.21 19.00 4 4.93 9.83 14.70 19.58 24.47 6 6.95 13.87 20.76 27.58 34.50 8 8.87 17.67 26.54 35.60 44.56 12 12.95 25.79 38.63 51.47 64.30 24 51.81 101.43 151.09 201.08 251.15 Table 5.1 MFHBO runtime (s) The programme runtime of the proposed MFHBO method has been measured on a laptop computer of model ASUS ROG Zephyrus M16 (GU603H) with a CPU of model i7-11800H and a GPU of model NVIDIA RTX3060 and the results for different combinations of number of frames and number of iterations are listed in Table 5.1. It can be concluded that the application of the proposed method is for pre-computed high-quality holograms, instead of real-time holographic projection. 5.3.2 Optical Experiment results The holographic projection system used in this experiment is the same as the one described in Section 3.2, with the optical setup shown in Fig. 3.2. Since the SLM has a refresh rate of 1440Hz and modern computer monitors have refresh rate of at least 60Hz, the maximum number of frames was chosen to be 1440/60 = 24, so that each set of 24 frames will take a total of 1/60 seconds to display, therefore giving an equivalent refresh rate of 60Hz. Then the integer factors of 24 were chosen so that the equivalent refresh rate becomes integer multiples of 60Hz. The number of frames starts from 1 to help illustrate how the increase in number of frames positively affect the reconstruction quality. 5.3 Results 97 1 2 3 4 6 8 12 24 12.0 11.5 11.0 10.5 10.0 9.5 9.0 8.5 Si m ul at io n re su lts O pt ic al r es ul ts N M SE (× 10 -8 ) Number of frames 0.14 0.12 0.20 0.18 0.26 0.24 0.30 0.28 0.22 0.16 0.10 SSIM Fig. 5.4 Simulation and optical reconstruction results for different number of frames The results in Fig. 5.4 further compares the final results for different number of frames. The histogram in Fig. 5.4 shows that, as the number of frames increases, the NMSE between the average reconstructions Ravg and the target image T decreases and the structural similarity index (SSIM)[58] increases, showing a trend of better reconstruction quality with higher number of frames. Such a trend is expected as more frames provide higher information capacity, which agrees with the information capacity research in Chapter 6 where holograms with higher bit depth were found to achieve better reconstruction quality. The trend is also shown visually via the simulation results and their detail enlargements. The corresponding multi-frame holograms are then loaded onto the SLM, and the reconstructed field is captured using a camera of model Canon EOS 1000D. Only the bottom halves of the reconstructed field were captured as the symmetrical conjugates were unwanted feature of far field projections from binary-phase SLM’s. The raw data including multi-frame binary-phase holograms, simulated reconstruction and optical results captured are accessible in the database [89]. 98 Multi-Frame Binary-Phase Holograms Batched Optimisation Fig. 5.5 Sample target image - ‘holography’ ambigram Then another target image was tested, which is the holography ambigram shown as shown in Fig. 5.5.1 The term ambigram is used to refer to (often typographical) designs that are invariant under a reflection, rotation or other symmetry. The ‘holography’ design contains 180-degree rotational symmetry, which makes it especially well suited to binary Fourier- holographic projection, where this symmetry is unavoidable. Multi-frame holograms were then generated using the proposed MFHBO method and the existing OSPR and AD-OSPR methods, for the same number of frames n = 24. The optical results are shown in Fig. 5.6. 1Adapted, with colours reversed, from holography - Benjamin Wetherfield, 2022. 5.3 Results 99 MFHBO OSPR AD-OSPR M an dr ill A m bi gr am Fig. 5.6 Optical results comparison of the proposed MFHBO method against the existing OSPR and AD-OSPR methods As shown in Fig. 5.6, for the Mandrill target image, it can be seen that the proposed MFHBO method achieved a much better optical reconstruction quality than the existing OSPR and AD-OSPR methods, with clearer details and better contrasts; for the ‘holography’ ambigram target image, the proposed MFHBO method is shown to have a much lower background noise around the centre, than the existing OSPR and AD-OSPR methods. The intended black regions are represented much more cleanly, with an elimination of speckle-like artefacts in the zero-valued space around the lettering, and an overall increase in discernible contrast. 100 Multi-Frame Binary-Phase Holograms Batched Optimisation Image Metric MFHBO OSPR AD-OSPR Mandrill NMSE (×10−4) 0.84 1.00 1.12 SSIM 0.124 0.076 0.078 Ambigram NMSE (×10−5) 2.29 3.31 3.23 SSIM 0.795 0.826 0.827 Table 5.2 Quantitative analysis of the optical results in Fig. 5.6 A quantitative analysis was then conducted on the optical results in Fig. 5.6, the NMSE and SSIM between the captured reconstructions and there corresponding targets are computed and listed in Table 5.2. The NMSE results of the proposed MFHBO method are lower than those of the existing OSPR and AD-OSPR methods, with a 25% reduction on average among both target images. On the other hand, the SSIM results have shown a 62% increase using MFHBO than OSPR and AD-OSPR for the mandrill target image, but a slight decrease of 3.7% for the ‘holography’ ambigram target image, which is negligible as it is less than 5% and the SSIM metric is not originally designed for binary-valued non-greyscale images. 5.3 Results 101 3D Holography O SP R M FH B O zzzzz1 2 3 4 TTTT1 2 3 4 Fig. 5.7 4-slice target and according reconstruction results The proposed MFHBO method was extended to multi-slice targets, by computing the loss between all 4 slices of reconstructions and target images (the Sum-of-Loss method in [4]). An example 4-slice target made from letters ‘A, B, C, D’ is shown in Fig. 5.7. The z values, corresponded to the z variable in Eq. (2.22), were chosen to be 1.1m,1.9m,3.5m,7.7m for the 4 slices respectively (as there’s no correlation between each slice, larger separation was chosen for fewer cross-talk across different planes). It can be seen that the proposed MFHBO method has produced sharper edges in reconstruction than the existing OSPR method. (The AD-OSPR method was not attempted here as its application to multi-slice targets was not defined). 102 Multi-Frame Binary-Phase Holograms Batched Optimisation Method Metric Slice 1 Slice 2 Slice 3 Slice 4 Average OSPR NMSE(×10−4) 4.980 4.484 5.644 4.846 4.988 SSIM 0.072 0.061 0.048 0.060 0.060 MFHO NMSE(×10−4) 4.484 4.230 4.990 4.289 4.498 SSIM 0.063 0.067 0.058 0.072 0.065 Table 5.3 Quantitative analysis of the optical results in Fig. 5.7 Then a quantitative analysis was carried out, with NMSE and SSIM values measured and shown in Table 5.3. The proposed MFHBO method has shown a 10% reduction in NMSE and a 8% improvement in SSIM on average than the existing OSPR method, demonstrating the effectiveness of the proposed method. On the laptop of model Mechrevo 16 Ultra, which has a CPU of model Intel Ultra 9 275HX and a GPU of model RTX5070Ti, using python with PyTorch and CUDA enabled, the optimisation took 31.249 seconds to complete for 1000 iterations with 24 frames, and the memory usage was 2GB. The runtime is much longer than the existing OSPR method, which takes < 0.5 seconds to complete for 1000 iterations with 24 frames, but the proposed MFHBO method achieves much better reconstruction quality, therefore it is suitable for pre-computed high-quality hologram applications. 5.3 Results 103 O SP R M FH B O T1 T2 T3 Fig. 5.8 Real-life captured image as target field and their reconstruction results Lastly, a set of real-life scene was captured in the lab using near, middle and far focus, as shown in T1,T2,T3 in Fig. 5.8 respectively. The z values were set to 1.1m,1.2m,1.3m for hologram generation, and the reconstruction results of the existing OSPR and the proposed MFHBO methods are compared in Fig. 5.8. The proposed MFHBO method is shown to have achieved better reconstruction quality than the existing OSPR method. 104 Multi-Frame Binary-Phase Holograms Batched Optimisation Method Metric Slice 1 Slice 2 Slice 3 Average OSPR NMSE(×10−6) 3.70 3.69 3.47 3.62 SSIM 0.37 0.28 0.34 0.33 MFHO NMSE(×10−6) 3.20 2.78 3.06 3.01 SSIM 0.42 0.32 0.32 0.35 Table 5.4 Quantitative analysis of the optical results in Fig. 5.8 A quantitative analysis was conducted again, with NMSE and SSIM values measured and listed in Table 5.4. The proposed MFHBO method has shown a 17% reduction in NMSE and a 7% improvement in SSIM on average than the existing OSPR method, proving the effectiveness of the proposed method. On the laptop of model Mechrevo 16 Ultra, which has a CPU of model Intel Ultra 9 275HX and a GPU of model RTX5070Ti, using python with PyTorch and CUDA enabled, the optimisation took around half an hour to complete for 1000 iterations with 24 frames, and the memory usage was 20GB. Therefore the proposed method is not suitable for real-time holographic projection or computation resource limited situations, but it is suitable for pre-computed high-quality hologram applications. 5.4 Summary This chapter proposed the MFHBO method, which generates multi-frame binary-phase holograms that are displayed on a binary-phase SLM with a high refresh rate. The proposed MFHBO method was shown to achieve much better reconstruction quality and higher contrast compared to the existing multi-frame binary-phase hologram generation methods OSPR [41] and AD-OSPR [65] on the holographic projector with a binary-phase SLM, for all the single-slice far-field targets and the multi-slice near-field targets tested. Although the proposed MFHBO method takes much longer time to compute than the existing OSPR and AD-OSPR methods, its much better reconstruction quality makes it suitable for pre-computed high-quality hologram applications. Its strong advantage for high contrast targets (such as the ‘holography’ ambigram in Fig. 5.5, which has a much reduced background noise) makes the MFHBO method particularly well-suited for photo-lithography applications. The proposed method can also be adapted for multi-level SLMs by simply removing the quantisation step in Fig. 5.1. This approach can be the suitable for applications such as photo-lithography, where 5.4 Summary 105 the system’s time response is much longer than that of human vision, and high refresh rates of the SLM are not required. Compared to the Lee’s method [90], the proposed MFHBO method is much more computationally heavy but it requires less optical components in the optical setup, the multi-frame technique can also exploit the high refresh rate of binary phase SLMs to average out the perceived speckle-like noise in the reconstruction. Chapter 6 Information Capacity of Phase-Only Computer-Generated Holograms Note: The work in this chapter has been published in Ref. [8] Despite many years of development in computer-generated holography, perfect phase-only holograms, whose reconstructions exactly match the targets with 0 error, are still not yet possible to compute. All computational phase retrieval algorithms end up with some error between the target image and the reconstruction of the computer-generated hologram (CGH), except for certain specific targets. This chapter focuses on the fundamental limits of phase- only CGH quantised to limited bit-depth levels, from an information theory perspective. It explores the information capacity of a CGH and its impact on reconstruction quality, and aims to quantify the difficulty of computing a phase-only CGH for a target image. 6.1 Introduction As introduced in the previous chapters, Holography is a technology that can fully reconstruct the wavefront of 3D objects. CGH is a technique for converting a 3D object scene into a 2D complex-valued hologram [91], without the need for the 3D object to physically exist. However, despite many advancements in liquid crystals and micro mirrors technologies, complex-valued SLM’s are still not available yet, and the currently available SLM’s can 108 Information Capacity of Phase-Only Computer-Generated Holograms only achieve either amplitude or phase modulation, among which the phase modulation is usually preferred in holographic projections for its lower zero order and higher energy efficiency due to less blockage of light. There are many algorithms available to compute good quality phase-only holograms as shown in the previous chapters; however, none of them can guarantee to compute a perfect phase hologram having zero error between the reconstructions and the according 3D or 2D targets, where they always end up with some error between the reconstruction of the hologram and the target scene, especially when the phase holograms are quantised to be able to display on SLM with limited bit depth. An intuition therefore arose that the bit depth of the phase hologram is limiting the reconstruction quality, and that the target scene’s entropy seems to denote how difficult it is for phase hologram generation. It is obvious that the entropy of the target can certainly never exceed the bit depth limit of its corresponded perfect hologram, otherwise a lossless compressor breaking the Shannon’s information theory [92] would be invented; however, such statement is not quite useful in practice as perfect holograms are generally impossible to find. The literature review had found no related work on the information entropy of computer- generated holograms, with the closest match being the research on hologram compression using optical method to achieve lossy compression [93], which cannot be integrated in CGH processes. Therefore, this chapter aims to investigate how much information content can a bit-depth-limited phase hologram contain, taking from an information theory point of view. Previously, work had been done to investigate the effect of bit-depth in Stochastic Gradient Descent (SGD) performance for phase-only computer-generated holography displays [3]. This chapter extends on previous research onto the Gerchberg-Saxton (GS) [22] algorithms for hologram generation, and investigates the correlation between the quality of the reconstructed image from the hologram and the information entropy of the target image, with an aim to reduce the hologram’s entropy during the CGH process, for smaller size holograms. 6.2 Methods 6.2.1 Quantised CGH Algorithm To compute phase-only holograms, the Gerchberg-Saxton (GS) [22] algorithm, which is the classic and robust phase retrieval algorithm introduced in Section 2.3.4, is selected. Taking the Fraunhofer propagation in Eq. (2.24) as an example, where the reconstructed field is simply the Fourier Transform (FFT when discretised) of the hologram field, the operation flowchart of GS algorithm is illustrated in Fig. 6.1. GS algorithm functions that it iteratively 6.2 Methods 109 determines the phase profile (∠A) of the hologram (A) required to reconstruct a target image (T ); it loops between the hologram (A) and the diffracted plane (E), and applying constraints to each plane accordingly during each iteration [22], with the total number of iterations denoted by N. Load T, N E=Te j∠E A=IFFT(E) E=FFT(A) A=e j∠A Output ∠A and |E| if n < N (∠A is random for n=0) n=n+1 Initiate n=0 if n = N Fig. 6.1 Gerchberg-Saxton (GS) algorithm flowchart Then, as illustrated in Fig. 6.2, a quantisation function (Q) can be defined by finding the closest point from one of the 2d quantisation levels, given the phase (∠A) and the quantisation bit depth d as input. 110 Information Capacity of Phase-Only Computer-Generated Holograms (a) bit depth = 1 (b) bit depth = 2 (c) bit depth = 4 (d) bit depth = 8 Fig. 6.2 Quantisation of phase holograms Load T, N E=Te j∠E A=IFFT(E) E=FFT(A) A=e jQ(∠A,d) Output ∠A and |E| if n < N (∠A is random for n=0) n=n+1 Initiate n=0 if n = N Fig. 6.3 Quantised Gerchberg-Saxton (GS) algorithm flowchart 6.2 Methods 111 To compute phase holograms quantised to certain bit depth (d), the GS algorithm is modified to include an additional quantisation operation (Q) when applying the ’phase-only’ constraint as shown in Fig. 6.3. Such a method is better than applying the quantisation at the end of the loop, as it includes the quantisation constraint throughout the iterations, instead of introducing significant quantisation noise in the end. 6.2.2 Measurement of Information Shannon entropy To quantify the information content, the classical one-dimensional (1D) Shannon entropy [92] with equation shown in Eq. (6.1) is selected. H(X) =−∑ x∈X p(x) log2 p(x) (6.1) Although the Shannon entropy was designed for 1D data, it can also be implemented to two- dimensional (2D) data by ignoring the 2D spatial correlations and summing p(x) log2 p(x) over the histogram of the 2D data. As only discrete data can have a meaningful Shannon entropy, the entropy can only be calculated for quantised holograms and target images, which are usually quantised to less than 8 bit depth. Delentropy To account for 2D spatial correlation, the delentropy [94] is also used. Delentropy is an extension of the 1D Shannon entropy that it first computes the gradient (del) vector field image, whose entropy is then named as the delentropy, so that the spatial image structure and pixel co-occurrence can be captured [94]. (a) Sample image (b) x-derivative (c) y-derivative Fig. 6.4 Del operation on a sample image 112 Information Capacity of Phase-Only Computer-Generated Holograms As an example, the sample image in Fig. 6.4 (a) is the file of filename ‘0500.png’ under the ‘DIV2K_train_HR’ folder sourced from the DIV2K dataset [95]. The sample image is calculated to have a Shannon entropy of 7.502 bits/pixel. HPGS( fx, fy)≤ H( fx)+H( fy) 2 (6.2) By taking the x-derivative ( fx) and y-derivative ( fy) as shown in Fig. 6.4 (b) and (c), and using the Papoulis Generalized Sampling (PGS) [96] theory, the delentropy is calculated using Eq. (6.2)[94] to be 5.867 bits/pixel. 6.3 Results 6.3.1 Targets in the far field (Fraunhofer region) In this subsection, the target images are set in the far field, so the Fraunhofer diffraction formula in Eq. (2.24) is used. Bit depth = 1 Bit depth = 2 Bit depth = 3 Bit depth = 4 Bit depth = 5 Bit depth = 6 Bit depth = 7 Bit depth = 8 NMSE=7.09×10-8 NMSE=3.67×10-8 NMSE=1.76×10-8 NMSE=9.18×10-9 NMSE=5.47×10-9 NMSE=4.23×10-9 NMSE=4.00×10-9 NMSE=3.97×10-9 H o lo g ra m R ec o n st ru ct io n Fig. 6.5 Holograms generated at bit depths level from 1 to 8 and their according reconstruc- tions in the far field An example run of the quantised GS algorithm (in Fig. 6.3) for the sample image in Fig. 6.4 (a) is shown in Fig. 6.5, which demonstrates qualitatively how the reconstruction quality im- proves with the increase in the bit depth of the hologram, and also quantitatively as the NMSE between the reconstruction and target image has shown a decreasing trend as the bit depth of hologram increases. The rotational symmetry in the reconstruction of the hologram with bit depth 1 is explained by the conjugate properties of Fourier Transforms in Section 2.2.3. 6.3 Results 113 Then the same quantised GS algorithm is run on the 800 images in the ‘DIV2K_train_HR’ folder in the DIV2K dataset [95], for hologram bit depth set to integers ranging from 1 to 8, with the total number of iterations (N) set to 100. 1 2 3 4 5 6 7 8 Hologram bit depth 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 AVG and STD for the 800 images Fig. 6.6 The average and standard deviation of the far-field reconstruction errors among the 800 target images plotted against the hologram bit depth For the 800 target images, the average (AVG) and standard deviation (STD) of the NMSE values between the far-field reconstructions of the holograms and their corresponding target images are plotted against the hologram bit depth in Fig. 6.6 (the raw data is accessible from the published research dataset [97]). It can be concluded that, the NMSE between the resulting reconstructions and their according target images decreases as the hologram bit depth increases. As holograms with higher bit depths can carry more information according to the information theory, they are more capable of storing the information needed to better reconstruct the target images. 114 Information Capacity of Phase-Only Computer-Generated Holograms 4 5 6 7 8 Target image entropy 0.0 0.5 1.0 1.5 2.0 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 Hologram bit depth = 1: Results of CGH on 800 images Hologram bit depth = 1: linear regression line (R2=0.01377) Hologram bit depth = 2: Results of CGH on 800 images Hologram bit depth = 2: linear regression line (R2=0.01123) Hologram bit depth = 3: Results of CGH on 800 images Hologram bit depth = 3: linear regression line (R2=0.01899) Hologram bit depth = 4: Results of CGH on 800 images Hologram bit depth = 4: linear regression line (R2=0.05290) Hologram bit depth = 5: Results of CGH on 800 images Hologram bit depth = 5: linear regression line (R2=0.10571) Hologram bit depth = 6: Results of CGH on 800 images Hologram bit depth = 6: linear regression line (R2=0.15517) Hologram bit depth = 7: Results of CGH on 800 images Hologram bit depth = 7: linear regression line (R2=0.19216) Hologram bit depth = 8: Results of CGH on 800 images Hologram bit depth = 8: linear regression line (R2=0.21355) Fig. 6.7 Scatter plot of the far-field reconstruction errors v.s. entropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints Then the entropy of each of the 800 target images are calculated using Eq. (6.1), and a scatter plot of the NMSE between the far-field reconstruction and target image against the target image entropy is plotted for all 800 target images run with each of the 8 hologram bit depth levels, as shown in Fig. 6.7. To avoid the effect of the different initial random phases on the final result, 5 different randomly generated initial phases are used for each run; however, the result [97] has shown negligible difference between the runs using different random initial phase holograms. Unfortunately, as the linear regressions between NMSE and entropies of target image have coefficients of determination (R2) much less than 0.5, no correlation has 6.3 Results 115 been found between the NMSE and the target image entropy. Therefore, the Shannon entropy cannot be used to quantify the difficulty of CGH for a given target image. 3 4 5 6 7 Target image delentropy 0.0 0.5 1.0 1.5 2.0 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 Hologram bit depth = 1: Results of CGH on 800 images Hologram bit depth = 1: linear regression line (R2=0.00003) Hologram bit depth = 2: Results of CGH on 800 images Hologram bit depth = 2: linear regression line (R2=0.00019) Hologram bit depth = 3: Results of CGH on 800 images Hologram bit depth = 3: linear regression line (R2=0.00001) Hologram bit depth = 4: Results of CGH on 800 images Hologram bit depth = 4: linear regression line (R2=0.00139) Hologram bit depth = 5: Results of CGH on 800 images Hologram bit depth = 5: linear regression line (R2=0.00660) Hologram bit depth = 6: Results of CGH on 800 images Hologram bit depth = 6: linear regression line (R2=0.01410) Hologram bit depth = 7: Results of CGH on 800 images Hologram bit depth = 7: linear regression line (R2=0.02149) Hologram bit depth = 8: Results of CGH on 800 images Hologram bit depth = 8: linear regression line (R2=0.02679) Fig. 6.8 Scatter plot of the far-field reconstruction errors v.s. delentropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints The delentropy of each of the 800 target images are then calculated using the method in Section 6.2.2, and a scatter plot of the NMSE between the far-field reconstruction and target image against the target image delentropy is plotted in Fig. 6.8 for the 800 target images. As all the R2 are much less than 0.5, no correlation is shown between the NMSE and the target image delentropy, inferring that the 2D delentropy also fails to predict the error of CGH for a given target image. 116 Information Capacity of Phase-Only Computer-Generated Holograms 6.3.2 Targets in the near field (Fresnel region) The target images are now set to near field, where the Fresnel diffraction formula in Eq. (2.22) applies; therefore the FFT and IFFT stages in Fig. 6.3 are modified to include the phase term in Eq. (2.22), and for experimental purpose, the distance (z) is set at 10cm, the hologram’s pixel pitch (sampling resolution of x and y) has a size of 13.62µm and the incident light’s wavelength (λ ) is 532nm. Bit depth = 1 Bit depth = 2 Bit depth = 6 Bit depth = 7Bit depth = 3 Bit depth = 4 Bit depth = 5 Bit depth = 8 NMSE=7.19×10-8 NMSE=3.6×10-8 NMSE=1.74×10-8 NMSE=8.69×10-9 NMSE=4.48×10-9 NMSE=2.39×10-9 NMSE=1.32×10-9 NMSE=7.78×10-10 H o lo g ra m R ec o n st ru ct io n Fig. 6.9 Holograms generated at bit depths level from 1 to 8 and their according reconstruc- tions in the near field Fig. 6.9 shows both qualitatively and quantitatively how the reconstruction quality improves (i.e. NMSE decreases) with the increase in the bit depth of the hologram. Such a trend is the same for target images placed in the near field as those placed in the far field in Fig. 6.5. The rotational symmetry is gone for the binary phase hologram (bit depth = 1) due to the extra phase term in the Fresnel diffraction formula making the product of binary phase hologram and the phase term to be complex-valued whose complex conjugate does not equal to itself; however, the complex conjugate wouldn’t disappear, but it will appear at a different distance to where the target image is set at, leading to extra defocused noise onto the reconstruction plane. Nevertheless, the trend of NMSE shows that holograms with higher bit depth produces better quality in the reconstruction plane. 6.3 Results 117 1 2 3 4 5 6 7 8 Hologram bit depth 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 AVG and STD for the 800 images Fig. 6.10 The average and standard deviation of the near-field reconstruction errors among the 800 target images plotted against the hologram bit depth Fig. 6.10 plots the average (AVG) and standard deviation (STD) of the NMSE values between the near-field reconstructions of the holograms and their corresponding target images against the hologram bit depth (the raw data is accessible from the published research dataset [97]). In Fig. 6.10, the same trend as the one for far field in Fig. 6.6 can be observed, where NMSE decreases as hologram bit depth increases, for every single target image. Therefore, the previous conclusion of higher hologram bit depth being able to produce better reconstruction quality is still valid. 118 Information Capacity of Phase-Only Computer-Generated Holograms 4 5 6 7 8 Target image entropy 0.0 0.5 1.0 1.5 2.0 2.5 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 Hologram bit depth = 1: Results of CGH on 800 images Hologram bit depth = 1: linear regression line (R2=0.02802) Hologram bit depth = 2: Results of CGH on 800 images Hologram bit depth = 2: linear regression line (R2=0.01169) Hologram bit depth = 3: Results of CGH on 800 images Hologram bit depth = 3: linear regression line (R2=0.01890) Hologram bit depth = 4: Results of CGH on 800 images Hologram bit depth = 4: linear regression line (R2=0.05302) Hologram bit depth = 5: Results of CGH on 800 images Hologram bit depth = 5: linear regression line (R2=0.10570) Hologram bit depth = 6: Results of CGH on 800 images Hologram bit depth = 6: linear regression line (R2=0.15503) Hologram bit depth = 7: Results of CGH on 800 images Hologram bit depth = 7: linear regression line (R2=0.19222) Hologram bit depth = 8: Results of CGH on 800 images Hologram bit depth = 8: linear regression line (R2=0.21349) Fig. 6.11 Scatter plot of the near-field reconstruction errors v.s. entropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints A scatter plot of the NMSE between the near-field reconstruction and target image against the target image entropy is plotted for all 800 target images run with each of the 8 hologram bit depth levels as shown in Fig. 6.11, with the hologram bit depth distinguished by different colours. Unfortunately, as the linear regressions between NMSE and entropies of target image have coefficients of determination (R2) much less than 0.5, no correlation has been found between the NMSE and the target image entropy. Therefore, the Shannon entropy cannot be used to quantify the difficulty of CGH for a given target image. 6.3 Results 119 3 4 5 6 7 Target image delentropy 0.0 0.5 1.0 1.5 2.0 2.5 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e 1e 7 Hologram bit depth = 1: Results of CGH on 800 images Hologram bit depth = 1: linear regression line (R2=0.00003) Hologram bit depth = 2: Results of CGH on 800 images Hologram bit depth = 2: linear regression line (R2=0.00016) Hologram bit depth = 3: Results of CGH on 800 images Hologram bit depth = 3: linear regression line (R2=0.00001) Hologram bit depth = 4: Results of CGH on 800 images Hologram bit depth = 4: linear regression line (R2=0.00145) Hologram bit depth = 5: Results of CGH on 800 images Hologram bit depth = 5: linear regression line (R2=0.00661) Hologram bit depth = 6: Results of CGH on 800 images Hologram bit depth = 6: linear regression line (R2=0.01410) Hologram bit depth = 7: Results of CGH on 800 images Hologram bit depth = 7: linear regression line (R2=0.02149) Hologram bit depth = 8: Results of CGH on 800 images Hologram bit depth = 8: linear regression line (R2=0.02680) Fig. 6.12 Scatter plot of the near-field reconstruction errors v.s. delentropy of target image among the 800 target images, with different colours indicating hologram bit depth constraints The scatter plot ofthe NMSE between the near-field reconstruction and target image against the target image delentropy in Fig. 6.12 again shows no correlation between the NMSE and the delentropy of target image. Such ‘no correlation’ result is the same as the results for far field investigated in Section 6.3.1, confirming that neither entropy nor delentropy is suitable for quantifying how difficult a target image is for CGH. 120 Information Capacity of Phase-Only Computer-Generated Holograms 0 1 2 3 4 5 6 7 8 9 10 iteration number 0 1 2 3 4 5 6 7 8 9 En tro py o f h ol og ra m (s ol id li ne s) holo bit depth = 8 holo bit depth = 7 holo bit depth = 6 holo bit depth = 5 holo bit depth = 4 holo bit depth = 3 holo bit depth = 2 holo bit depth = 1 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e (d ot te d lin es ) 1e 7 (a) 0 1 2 3 4 5 6 7 8 9 10 iteration number 0 1 2 3 4 5 6 7 8 9 En tro py o f h ol og ra m (s ol id li ne s) holo bit depth = 8 holo bit depth = 7 holo bit depth = 6 holo bit depth = 5 holo bit depth = 4 holo bit depth = 3 holo bit depth = 2 holo bit depth = 1 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 NM SE b et we en re co ns tru ct io n an d ta rg et im ag e (d ot te d lin es ) 1e 7 (b) Fig. 6.13 Hologram entropy (solid lines) and NMSE (dotted lines) plotted against the iteration number in a GS run with the initial hologram phase (∠A) being (a) zeros (b) random 6.3 Results 121 The research moves on to investigate the entropy of holograms. In an example run for a quantised hologram generation in the near field, both the hologram entropy and the NMSE between reconstruction and target image are recorded and plotted in Fig. 6.13 (a) and (b) for the first 20 iterations. As randomly generated holograms will by definition create holograms with high entropy, additional experiment has been carried out with initial phase set to zeros (i.e. ∠A is set to zeros at n = 0 in Fig. 6.3). Therefore, two diagrams are plotted in Fig. 6.13, with Fig. 6.13 (a) having initial phase of zeros and Fig. 6.13 (b) having random initial phase. Both Fig. 6.13 (a) and (b) have the horizontal axis to be the iteration number n, and the vertical axis in the left is the entropy of hologram (corresponding to solid lines in the plot), while the vertical axis in the right is the NMSE (corresponding to dotted lines in the plot). Colour coding is used to distinguish between the 8 different runs where hologram bit depth is set from 1 to 8. The solid lines in Fig. 6.13 (a) show that the entropy of hologram keeps increasing towards a value lower than the bit depth, with their corresponding NMSE between reconstruction and target image (dotted lines) decreasing. Such a trend can be explained qualitatively that, as the iteration goes on, the hologram is attempting to contain more information to sustain a better reconstruction, while the entropy cannot exceed or even reach the bit depth level. On the other hand, the random initial phase plotted in Fig. 6.13 (b) has a constant entropy approximately equal to the bit depth, which infers that the iterations are improving the reconstruction quality without reducing the information entropy of the hologram. In both cases, the entropy of the hologram does not decrease at any iteration, and the final NMSE does not have significant reduction when the hologram’s bit depth exceeds 5. The entropy of final hologram in Fig. 6.13 (a) is lower than that in Fig. 6.13 (b), although the final NMSE are similar. If the computer-generated holograms were to undergo a lossless compression, the hologram in Fig. 6.13 (a) would be better compressed than the hologram in Fig. 6.13 (b) as the entropy denotes the compression limit, assuming that the holograms are treated as 1D arrays when compressing. Therefore if hologram compression is of a concern, then starting with a low entropy initial hologram in the CGH process is recommended. In case if the reconstruction quality is not as high of a priority than making holograms to occupy less storage space, it is also recommended to use quantised GS algorithm with 5 bit depth instead of 6-7 bit depth as the final reconstruction quality will not degrade significantly. 122 Information Capacity of Phase-Only Computer-Generated Holograms 6.4 Summary By carrying out the quantised GS algorithm on 800 sample target images placed at both far field (Fraunhofer diffraction) and near field (Fresnel diffraction), and computing the entropy and delentropy of the target images and holograms, this chapter concluded that holograms with higher bit depth can sustain more information, thereby producing reconstructions with better quality. However, the quality of the reconstruction is not correlated to either the entropy or the delentropy of the target image, so neither entropy nor delentropy can quantify how difficult an image is for phase-only hologram generation. Moreover, the entropy of the hologram generated using the quantised GS algorithm is not only bounded by the hologram bit depth, but also affected by the entropy of the initial phase. For applications where hologram file size is a high priority, it is advised to start with a low entropy initial phase (e.g. all zeros) rather than a random initial phase. It is also recommended to reduce the hologram bit depth limit, so that holograms with lower entropy will be generated. Chapter 7 Summary and Future Work 7.1 Summary The research presented in this thesis has focused on advancing the field of computer-generated holography (CGH) through the development and optimisation of phase-only holograms. By exploring and implementing various algorithms, this thesis has identified the strengths and weaknesses of each approach in terms of computational efficiency and reconstruction quality. Chapter 3 proposed the Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) al- gorithm, described the experimental setup, explained the method for determining the DPD curve and then applied the DPD curve to improve holographic projection quality. The results demonstrated that DPD can significantly enhance the quality of the reconstructed images. On the grey-scale ramp target, the DPD-OSPR method reduced the NMSE by a stunning 95.45%. Then the DPD-OSPR was applied on two sample images. It was observed that more details were shown in the replay field, and the NMSE’s of the two example images were reduced by 19.86% and 15.64% respectively. As the DPD is a one-to-one mapping, the extra computation required is negligible. The effectiveness of the proposed DPD-OSPR method to improve reconstruction quality on the existing OSPR algorithm while still keeping its ability for real-time holography was demonstrated. Chapter 4 focused on the optimisation of phase-only holograms, and proved the effectiveness of using the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm for CGH. Then the novel Target Image Phase Optimisation (TIPO) technique was proposed, which optimises the phase of the target image instead of the phase of the hologram. The 124 Summary and Future Work L-BFGS algorithm was shown to be able to offer efficient convergence, improving the quality of the generated holograms. Then the two existing 3D CGH optimisation methods - the sum-of-hologram (SoH) and sum-of-loss(SoL) techniques - were investigated. The novel method using L-BFGS optimiser with Sequential Slicing (SS) technique was proposed to generate phase-only hologram for multi-depth target, which was shown to be faster than the SoL technique and have better quality than the SoH technique. The L-BFGS with SS technique has demonstrated a good suppression on the quality imbalance across the multi- depth slices, benefiting from the nature of L-BFGS being a second order optimiser, which implicitly records the historical gradients by other slices for the determination of the descent direction. For both GD and L-BFGS optimisation algorithms, the SS technique runs faster and produces better reconstruction quality than the simple SoH technique. The proposed SS method also proved effective for complicated 3D targets and demonstrated great ability of time-limited applications. Chapter 5 proposed the Multi-Frame Holograms Batched Optimisation (MFHBO) algorithm to generate multi-frame binary-phase holograms to be displayed on the high-refresh-rate binary-phase SLM in the lab. Comparisons of simulation and optical experiment results showed that the proposed MFHBO method had superior performance in reducing noise and improving reconstruction quality for multi-frame holograms than did the existing multi- frame binary-phase holograms generation methods OSPR and AD-OSPR on the holographic projector with binary-phase SLM, for all the single-slice far-field targets and the multi-slice near-field targets tested. Although the proposed MFHBO method is slower than the existing OSPR and AD-OSPR methods, its much better reconstruction quality makes it especially suitable for pre-computed high-quality hologram applications. Its strong advantage for high contrast target also makes it well-suited for photo-lithography applications. The proposed method can also be adapted for multi-level SLM’s in the future once high-refresh-rate multi-level phase SLM’s are available. Finally, Chapter 6 explored the information capacity of phase-only CGH. This chapter examined quantisation effects on hologram bit depth and their impact on reconstruction quality, and reached the conclusion that holograms with higher bit depth can sustain more information and can therefore produce reconstructions with better quality. However, the quality of the reconstruction is not correlated to either the entropy or the delentropy of the target image, so neither entropy nor delentropy can quantify how difficult an image is for phase-only hologram generation. Moreover, the entropy of the hologram generated using quantised GS algorithm is not only bounded by the hologram bit depth, but also affected by 7.2 Future Work 125 the entropy of the initial phase. For applications where holograms file size is a high priority, it is advised to start with a low entropy initial phase rather than a random initial phase, and it is recommended to reduce the hologram bit depth limit. In conclusion, this thesis has contributed a series of novel phase retrieval techniques, with improvements achieved in the computational speed and reconstruction quality of phase-only CGH. While current methods provide a foundation for practical implementations, there still remains substantial room for improvement, particularly in achieving real-time, high-quality holography. Future research should continue to build on these findings and explore novel algorithms and advancing technologies to realise the full potential of CGH. 7.2 Future Work The experimental work in this thesis in Chapter 3 and Chapter 5 are based on a high-refresh- rate binary-phase SLM (as introduced in Section 3.2). In the future, if high-refresh-rate multi-level-phase modulators become available, it will be recommended to adapt the MFHBO method proposed in Chapter 5 to generate multi-frame multi-level phase holograms by simply setting the quantisation step in Fig. 5.1 to match the bit-depth level of the SLM. With a higher bit-depth level of the SLM, the reconstruction quality can be further improved. Currently available SLMs on the market not only lack in bit-depth levels and refresh rates, but also lack in the size and the number of pixels. If an SLM of the size of a regular computer monitor (e.g. more than 13 inch) having trillions of pixels were to be available in the future, it should have similar or better visual effects compared to the analogue film-based holograms such as the one in Fig. 1.1, but with computer-generated contents. In Chapter 5, the MFHBO method was proposed to generate multi-frame binary-phase holograms for high-refresh-rate binary-phase SLMs. It provided good reconstruction quality but the computational speed was much slower than the existing binary phase hologram generation methods. In the future, it would be ideal to further improve the computational speed of the MFHBO method, so that it can be used for real-time applications. The proposed MFHBO method also consumes a lot of memory, as it requires storing and optimising the holograms for each frame in the batch, which is not suitable for applications with limited memory resources. Future work could focus on reducing the memory footprint of the MFHBO method, for example by optimising the holograms in a pipeline manner instead of storing all frames in memory at once. 126 Summary and Future Work In Chapter 6, attempts were made to find a metric to quantify how difficult any image is for phase hologram generation, with entropy and delentropy being ruled out. In the future, more research can be done to find a suitable metric to predict phase hologram generation difficulty. If found, it will set a general standard on how CGH algorithms perform and how far away they are from the fundamental limit. The analogy is information entropy in data compression, where entropy denotes the lower bound of compression rate for any compression algorithm (i.e. a metric is needed to indicate the limit of reconstruction quality of the phase hologram generated for any target image). However, if a fully complex modulator could be built in the future, all the issues with phase retrieval algorithms can be avoided. Before the invention of a fully complex SLM, the current phase-only SLMs still has rooms of improvements as they are currently pixelated and arranged in a grid, leading to errors and noise in quantisation and discretisation. The gaps between the adjacent pixels reduce the active area of the SLM and cause mis-matches between simulations and optical experiments. In the future, it would be ideal to develop a new type of non-pixelated SLM that can modulate the phase continuously. It could be addressed acoustically, for example a thin layer of liquid that is modulated by ultrasonic waves. The continuous phase modulation will eliminate the quantisation and discretisation errors, and the continuous phase profile will also eliminate the gaps between pixels, leading to 100% active area of the SLM and achieve better matches between simulations and optical experiments. References [1] Jana Skirnewskaja, Yunuen Montelongo, Jinze Sha, and Timothy D. Wilkinson. Holo- graphic lidar projections with brightness control. In Imaging and Applied Optics Congress 2022 (3D, AOA, COSI, ISA, pcAOP), page 3F2A.6. Optica Publishing Group, 2022. [2] Jinze Sha, Andrew Kadis, Fan Yang, and Timothy D. Wilkinson. Limited-memory bfgs optimisation of phase-only computer-generated hologram for fraunhofer diffraction. In Digital Holography and 3-D Imaging 2022, page W3A.3. Optica Publishing Group, 2022. [3] Andrew Kadis, Benjamin Wetherfield, Jinze Sha, Fan Yang, Youchao Wang, and Timothy D. Wilkinson. Effect of bit-depth in stochastic gradient descent performance for phase-only computer-generated holography displays. London Imaging Meeting, 3:36–40, 7 2022. [4] Jinze Sha, Andrew Kadis, Fan Yang, Youchao Wang, and Timothy D. Wilkinson. Multi- depth phase-only hologram optimization using the l-bfgs algorithm with sequential slicing. J. Opt. Soc. Am. A, 40(4):B25–B32, Apr 2023. [5] Jinze Sha, Adam Goldney, Andrew Kadis, Jana Skirnewskaja, and Timothy D. Wilkin- son. Digital pre-distorted one-step phase retrieval algorithm for real-time hologram generation for holographic displays. Journal of Imaging Science and Technology, 67(3):030405–1–030405–1, 2023. [6] Jana Skirnewskaja, Yunuen Montelongo, Jinze Sha, Phil Wilkes, and Timothy D. Wilkinson. Accelerated augmented reality holographic 4k video projections based on lidar point clouds for automotive head-up displays. Advanced Optical Materials, 12(12):2301772, 2024. [7] Roubing Meng, Jinze Sha, Zhongling Huang, and Timothy D. Wilkinson. Extending FOV of holographic display with alternating lasers. In Peter Schelkens and Tomasz Kozacki, editors, Optics, Photonics, and Digital Technologies for Imaging Applications VIII, volume 12998, page 129981J. International Society for Optics and Photonics, SPIE, 2024. [8] Jinze Sha, Andrew Kadis, Benjamin Wetherfield, Roubing Meng, Zhongling Huang, Dilawer Singh, Antoni Wojcik, and Timothy D. Wilkinson. Information capacity of phase-only computer-generated holograms for holographic displays. In Peter Schelkens and Tomasz Kozacki, editors, Optics, Photonics, and Digital Technologies for Imaging 128 References Applications VIII, volume 12998, page 129980J. International Society for Optics and Photonics, SPIE, 2024. [9] Jinze Sha, Antoni Wojcik, Benjamin Wetherfield, Jianghan Yu, and Timothy D. Wilkin- son. Multi frame holograms batched optimization for binary phase spatial light modula- tors. Scientific Reports, 14(1):19380, Aug 2024. [10] Jinze Sha, Antoni Wojcik, and Timothy D. Wilkinson. Target image phase optimisation hologram generation method for phase-only spatial light modulators. In Practical Holography XXXIX: Displays, Materials, and Applications, volume 13390. Interna- tional Society for Optics and Photonics, SPIE, 2025. [11] John P. McIntire, Paul R. Havig, and Eric E. Geiselman. Stereoscopic 3d displays and human performance: A comprehensive review. Displays, 35(1):18–26, 2014. [12] Simon J. Watt, Kurt Akeley, Marc O. Ernst, and Martin S. Banks. Focus cues affect perceived depth. Journal of Vision, 5(10):7–7, 12 2005. [13] Barbara Klinger. Three-dimensional cinema: The new normal. Convergence, 19(4):423– 431, 2013. [14] Catherine Allen and Verity McIntosh. Safeguarding the metaverse. The Institution of Engineering and Technology, 2022. [15] Michael E. McCauley and Thomas J. Sharkey. Cybersickness: Perception of self- motion in virtual environments. Presence: Teleoperators and Virtual Environments, 1(3):311–318, 08 1992. [16] Hyun Taek Kim Eunhee Chang and Byounghyun Yoo. Virtual reality sickness: A review of causes and measurements. International Journal of Human-Computer Interaction, 36(17):1658–1682, 2020. [17] Ivan Y. Lo. A photo of the holographic portrait of dennis gabor, 2018. [18] D. Gabor. A new microscopic principle. Nature, 161:777–778, 1948. [19] Eugene Hecht. Optics. Pearson Education Limited, 5 edition, 2017. [20] Michael A. Seldowitz, Jan P. Allebach, and Donald W. Sweeney. Synthesis of digital holograms by direct binary search. Applied Optics, 26, 1987. [21] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 220, 1983. [22] R W Gerchberg. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik, 35:237–246, 1972. [23] Jingzhao Zhang, Nicolas Pégard, Jingshan Zhong, Hillel Adesnik, and Laura Waller. 3d computer-generated holography by non-convex optimization. Optica, 4:1306, 10 2017. [24] Shujian Liu and Yasuhiro Takaki. Optimization of phase-only computer-generated holograms based on the gradient descent method. Applied Sciences (Switzerland), 10, 2020. References 129 [25] Chun Chen, Byounghyo Lee, Nan-Nan Li, Minseok Chae, Di Wang, Qiong-Hua Wang, and Byoungho Lee. Multi-depth hologram generation using stochastic gradient descent algorithm with complex loss function. Optics Express, 29:15089, 5 2021. [26] Suyeon Choi, Jonghyun Kim, Yifan Peng, and Gordon Wetzstein. Optimizing image quality for holographic near-eye displays with michelson holography. Optica, 8:143, 2 2021. [27] Isaac Newton. Opticks. Dover Press, 1704. [28] Christiaan Huygens. Traite de la lumiere. Où sont expliquées les causes de ce qui luy arrive dans la reflexion, & dans la refraction. Et particulierment dans l’etrange refraction du cristal d’Islande, par C.H.D.Z. Avec un Discours de la cause de la pesanteur. chez Pierre Vander Aa marchand libraire, 1690. [29] Alasdair Beal. Who invented young’s modulus? 78, 07 2000. [30] Thomas Young. Ii. the bakerian lecture. on the theory of light and colours. Philosophical Transactions of the Royal Society of London, 92:12–48, 1802. [31] Augustin Jean Fresnel. Memoir on the Diffraction of Light. 1826. [32] John Daintith. A Dictionary of Physics. Oxford University Press, 2009. [33] Timothy D. Wilkinson. Electrical Data Book. Cambridge University Engineering Department, 2017. [34] James Clerk Maxwell. Viii. a dynamical theory of the electromagnetic field. Philosoph- ical Transactions of the Royal Society of London, 155:459–512, 1865. [35] A. Einstein. On a heuristic point of view about the creation and conversion of light. Annalen der Physik, 322(6):132–148, 1905. [36] Yuanbo Deng and Daping Chu. Coherence properties of different light sources and their effect on the image sharpness and speckle of holographic displays. Scientific Reports, 7(1):5893, Jul 2017. [37] Jeff Hecht. Solid-State and Fiber Lasers, chapter 8, pages 223–263. John Wiley & Sons, Ltd, 2008. [38] Gould R. Gordon. The laser, light amplification by stimulated emission of radiation. page 128. Ann Arbor, 6 1959. [39] Edwin Cartlidge. Theodore maiman 1927-2007. Physics World, 20, 2007. [40] John B Develis, Merrimack College, North Andover, and George O Reynolds. Three dimensional hologram reconstruction and image speckle, 1966. [41] A. J. Cable, E. Buckley, P. Mash, N. A. Lawrence, T. D. Wilkinson, and W. A. Cross- land. 53.1: Real-time binary hologram generation for high-quality video projection applications. SID Symposium Digest of Technical Papers, 35:1431, 2004. 130 References [42] Tim Stangner, Hanqing Zhang, Tobias Dahlberg, Krister Wiklund, and Magnus An- dersson. Step-by-step guide to reduce spatial coherence of laser light using a rotating ground glass diffuser. Applied Optics, 56:5427, 7 2017. [43] Linxiao Deng, Tianhao Dong, Yuwei Fang, Yuhua Yang, Chun Gu, Hai Ming, and Lixin Xu. Speckle reduction in laser projection based on a rotating ball lens. Optics and Laser Technology, 135, 3 2021. [44] Philip J W Hands, Calum M Brown, Daisy K E Dickinson, Stephen M Morris, and Jia-De Lin. Liquid-crystal lasers: Recent advances and future opportunities, 2022. [45] Arne Nordmann. Wave diffraction in the manner of huygens and fresnel, 2007. [46] Joseph W. Goodman. Introduction to Fourier Optics, Fourth Edition. W. H. Freeman, 2017. [47] Timothy D. Wilkinson. Lecture notes of 4b11 photonics systems course, 2019. Univer- sity of Cambridge. [48] T D Wilkinson, W A Crossland, S T Warr, T C B Yu, A B Davey, and R J Mears. New applications for ferroelectric liquid crystals. Liquid Crystals Today, 4(3):1–6, 1994. [49] A. J. Cable. Real-time high-quality two and three-dimensional holographic video projec- tion using the one-step phase retrieval (ospr) approach, 2006. PhD thesis, Department of Engineering, University of Cambridge, United Kingdom. [50] M. Schadt and W. Helfrich. Voltage-dependent optical activity of a twisted nematic liquid crystal. Applied Physics Letters, 18, 1971. [51] Dennis R. Pape and Larry J. Hornbeck. Characteristics of the deformable mirror device for optical information processing. Optical Engineering, 22, 1983. [52] Kristina M. Johnson, Douglas J. McKnight, and Ian Underwood. Smart spatial light modulators using liquid crystals on silicon. IEEE Journal of Quantum Electronics, 29, 1993. [53] Yongmin Lee, James Gourlay, William J. Hossack, Ian Underwood, and Anthony J. Walton. Multi-phase modulation for nematic liquid crystal on silicon backplane spatial light modulators using pulse-width modulation driving scheme. Optics Communications, 236, 2004. [54] S. E. Broomfield, M. A.A. Neil, E. G.S. Paige, and G. G. Yang. Programmable binary phase-only optical device based on ferroelectric liquid crystal slm. Electronics Letters, 28, 1992. [55] Allan Weber. Sipi image database - misc, 2022. [retrieved 17 Nov 2022]. [56] Zhengzhong Huang and Liangcai Cao. Quantitative phase imaging based on holography: trends and new perspectives. Light: Science & Applications, 13(1):145, Jun 2024. [57] Claude Sammut and Geoffrey I. Webb, editors. Mean Squared Error, pages 653–653. Springer US, Boston, MA, 2010. References 131 [58] Zhou Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli. Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 13(4):600–612, 2004. [59] Tomoyoshi Shimobaba and Tomoyoshi Ito. Random phase-free computer-generated hologram. Opt. Express, 23(7):9549–9554, Apr 2015. [60] Han Jin Yang, Jeong Sik Cho, and Yong Hyub Won. Reduction of reconstruction errors in kinoform cghs by modified simulated annealing algorithm. Journal of the Optical Society of Korea, 13, 2009. [61] Guo zhen Yang, Bi zhen Dong, Ben yuan Gu, Jie yao Zhuang, and Okan K. Ersoy. Gerchberg–saxton and yang–gu algorithms for phase retrieval in a nonunitary transform system: a comparison. Appl. Opt., 33(2):209–218, Jan 1994. [62] Haichao Wang, Weirui Yue, Qiang Song, Jingdan Liu, and Guohai Situ. A hybrid gerchberg-saxton-like algorithm for doe and cgh calculation. Optics and Lasers in Engineering, 89:109–115, 2017. 3DIM-DS 2015: Optical Image Processing in the context of 3D Imaging, Metrology, and Data Security. [63] Pengcheng Zhou, Yan Li, Shuxin Liu, and Yikai Su. Dynamic compensatory gerchberg- saxton algorithm for multiple-plane reconstruction in holographic displays. Optics Express, 27:8958, 3 2019. [64] Edward Buckley. 70.2: Invited paper: Holographic laser projection technology. SID Symposium Digest of Technical Papers, 39(1):1074–1079, 2008. [65] Andrzej Kaczorowski, George S. D. Gordon, and Timothy D. Wilkinson. Adaptive, spatially-varying aberration correction for real-time holographic projectors. Opt. Ex- press, 24(14):15742–15756, Jul 2016. [66] Xu Guan, Su Jian, Pan Hongda, Zhang Zhiguo, and Gong Haibin. An image enhance- ment method based on gamma correction. In 2009 Second International Symposium on Computational Intelligence and Design, volume 1, pages 60–63, 2009. [67] Sung Jin Kang and Sung Il Chien. Apl-adaptive inverse gamma correction for improving gray-level linearity of pdp-tv. Molecular Crystals and Liquid Crystals, 499(1):185/[507]–192/[514], 2009. [68] Po-Ming Lee and Hung-Yi Chen. Adjustable gamma correction circuit for tft lcd. In 2005 IEEE International Symposium on Circuits and Systems (ISCAS), pages 780–783 Vol. 1, 2005. [69] C. Alejandro Párraga, Jordi Roca-Vila, Dimosthenis Karatzas, and Sophie M. Wuerger. Limitations of visual gamma corrections in lcd displays. Displays, 35:227–239, 2014. [70] Charles Poynton. Digital Video and HD: Algorithms and Interfaces. Morgan Kaufmann, 2 edition, 2012. [71] J. Freeman. Visor projected helmet mounted display for fast jet aviators using a fourier video projector, 2009. PhD thesis, Department of Engineering, University of Cambridge, United Kingdom. 132 References [72] The MathWorks Inc. Matlab version: 9.13.0 (r2022b), 2022. [73] Xuetun Zhao. Suzhou center mall, 2017. Suzhou, Jiangsu, China. [74] Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, 2006. [75] Xiaomeng Sui, Zehao He, Daping Chu, and Liangcai Cao. Non-convex optimization for inverse problem solving in computer-generated holography. Light: Science & Applications, 13(1):158, Jul 2024. [76] John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12:2121– 2159, 2011. [77] Tijmen Tieleman, Geoffrey Hinton, et al. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural networks for machine learning, 4:26–31, 2012. [78] Diederik P. Kingma and Jimmy Lei Ba. Adam: A method for stochastic optimization. 2015. [79] Dong C. Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical Programming, 45, 1989. [80] Jhen-Si Chen, Daping Chu, and Quinn Y. Smithwick. Rapid hologram generation utilizing layer-based approach and graphic rendering for realistic three-dimensional image reconstruction by angular tiling. Journal of Electronic Imaging, 23(2):023016, 2014. [81] J.-S Chen and Daping Chu. Improved layer-based method for rapid hologram generation and real-time interactive holographic display applications. Optics Express, 23:18143– 18155, 07 2015. [82] Jia Jia, Jhen Si, and Daping Chu. Fast two-step layer-based method for computer gener- ated hologram using sub-sparse 2d fast fourier transform. Opt. Express, 26(13):17487– 17497, Jun 2018. [83] G. Cybenko, D.P. O’Leary, and J. Rissanen. The Mathematics of Information Coding, Extraction and Distribution. The IMA Volumes in Mathematics and its Applications. Springer New York, 1998. [84] S. Kullback and R. A. Leibler. On Information and Sufficiency. The Annals of Mathematical Statistics, 22(1):79 – 86, 1951. [85] Xunying Liu, Shansong Liu, Jinze Sha, Jianwei Yu, Zhiyuan Xu, Xie Chen, and Helen Meng. Limited-memory bfgs optimization of recurrent neural network language models for speech recognition. ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings, 2018-April:6114–6118, 9 2018. [86] Jun Amako, Hirotsuna Miura, and Tomio Sonehara. Speckle-noise reduction on kino- form reconstruction using a phase-only spatial light modulator. Appl. Opt., 34(17):3165– 3171, Jun 1995. References 133 [87] Nicolas Bacaër. Verhulst and the logistic equation (1838), pages 35–39. Springer London, London, 2011. [88] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmai- son, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. volume 32, 2019. [89] Jinze Sha, Antoni Wojcik, Benjamin Wetherfield, Jianghan Yu, and Timothy D. Wilkin- son. Research data supporting “multi-frame holograms batched optimization”. Apollo - University of Cambridge Repository, 2024. [90] Wai-Hon Lee. Iii computer-generated holograms: Techniques and applications. vol- ume 16 of Progress in Optics, pages 119–232. Elsevier, 1978. [91] P. W. M. Tsang, T.-C. Poon, and Y. M. Wu. Review of fast methods for point-based computer-generated holography. Photon. Res., 6(9):837–846, Sep 2018. [92] C. E. Shannon. A mathematical theory of communication. Bell System Technical Journal, 27, 1948. [93] Joel S Kollin. Design and information considerations for holographic television, 1988. [94] Kieran G. Larkin. Reflections on shannon information: In search of a natural information-entropy for images, 2016. [95] Eirikur Agustsson and Radu Timofte. Ntire 2017 challenge on single image super- resolution: Dataset and study. July 2017. [96] A. Papoulis. Generalized sampling expansion. IEEE Transactions on Circuits and Systems, 24(11):652–654, 1977. [97] Jinze Sha, Andrew Kadis, Benjamin Wetherfield, Roubing Meng, Zhongling Huang, Dilawer Singh, Antoni Wojcik, and Timothy D. Wilkinson. Research data supporting ‘information capacity of phase-only computer-generated holograms for holographic displays’, 2024. Table of contents List of figures List of tables Nomenclature 1 Introduction 2 Background Theory 2.1 The Nature of Light 2.1.1 Wave-Particle Duality 2.1.2 Wave Equation 2.2 Fundamentals of Holography 2.2.1 Light Source 2.2.2 Diffraction 2.2.3 Spatial Light Modulator (SLM) 2.3 Computer-Generated Holography (CGH) 2.3.1 Direct Phase-Only Hologram Generation 2.3.2 Direct Binary Search (DBS) Algorithm 2.3.3 Simulated Annealing (SA) Algorithm 2.3.4 Gerchberg-Saxton (GS) Algorithm 2.3.5 One-Step Phase Retrieval (OSPR) Algorithm 2.3.6 Adaptive One-Step Phase Retrieval (AD-OSPR) Algorithm 2.3.7 3D CGH 3 Digital Pre-Distorted One-Step Phase Retrieval (DPD-OSPR) Algorithm 3.1 Introduction 3.2 Experimental setup 3.3 Determining the DPD curve 3.4 Applying the DPD Curve 3.5 Summary 4 L-BFGS Optimisation of 2D and 3D CGH 4.1 Numerical Optimisation Methods 4.1.1 Optimisation framework 4.1.2 Gradient Descent 4.1.3 Newton's Method 4.1.4 Quasi-Newton Method 4.1.5 Large Scale Quasi-Newton Method: Limited Memory BFGS (L-BFGS) 4.2 Phase-only Hologram Optimisation 4.3 Target Image Phase Optimisation (TIPO) 4.4 Multi-Depth Phase-Only Hologram Optimisation 4.4.1 Methods 4.4.2 Results 4.5 Summary 5 Multi-Frame Binary-Phase Holograms Batched Optimisation 5.1 Introduction 5.2 Methods 5.3 Results 5.3.1 Simulation results 5.3.2 Optical Experiment results 5.4 Summary 6 Information Capacity of Phase-Only Computer-Generated Holograms 6.1 Introduction 6.2 Methods 6.2.1 Quantised CGH Algorithm 6.2.2 Measurement of Information 6.3 Results 6.3.1 Targets in the far field (Fraunhofer region) 6.3.2 Targets in the near field (Fresnel region) 6.4 Summary 7 Summary and Future Work 7.1 Summary 7.2 Future Work References