Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Peng, Yifan
Doctor of Philosophy - PhD
University of British Columbia
Diffractive optical elements (DOEs) have been studied extensively in optics for decades, but have recently received a lot of renewed interests in the context of computational imaging and computational display because they can drastically reduce the size and weight of devices. However, the inherent strong dispersion is an obstacle that limits the use of DOEs in full spectrum imaging, causing unacceptable color fidelity loss in the captured or reconstructed images. Despite the benefits of facilitating compact form factors, DOEs have sufficient degrees of freedom that one can manipulate to encode the desirable light modulation functionality. In this dissertation we theoretically and experimentally investigate the practicability of introducing numerical optimization into the design procedure of optics, to enable a variety of diffractive optics subject to different application scenarios. Regarding imaging applications, we first validate the practicality of introducing diffractive-refractive hybrid elements as simplified optics in conventional computational imaging. We re-implement a cross-channel based deconvolution to correct the chromatic aberration. The full fabrication cycle of photolithography technique is developed, serving as the basis for all designs. Then, the desirable focal powers in both spectral and spatial domain are encoded onto the DOEs. Precisely, we develop the diffractive achromat that balances the focusing contributions in full visible spectrum. The color fidelity can thus be well preserved. Meanwhile, we develop the encoded lenses that provide focus tunable imaging ability and multi-focal sweep imaging compromise. The trade-off light loss and residual aberrations are tackled by a deconvolution step. In addition, we extend the design paradigm from image capture to image display, where the holograms of visualization of multiple narrowband spectra are encoded onto a pair of ultrathin diffractive phase plates. The mix-and-match scheme of designing high degree-of-freedom diffractive optics is exploited via the complex matrix factorization. The essential random appearance of holograms in fact benefits to enabling a wide range of visualization applications. We envision the work on computational imaging, involving both capture and display end, provide new insights on incorporating optics and computation algorithms to better record, understand and deliver desirable visual information under the constraint of current data bandwidth of hardware.
Computational Imaging with Diffractive OpticsbyYifan PengMSc, Zhejiang University, 2013BE, Zhejiang University, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)April 2018c© Yifan Peng, 2018AbstractDiffractive optical elements (DOEs) have been studied extensively in optics fordecades, but have recently received a lot of renewed interests in the context ofcomputational imaging and computational display because they can drastically re-duce the size and weight of devices. However, the inherent strong dispersion is anobstacle that limits the use of DOEs in full spectrum imaging, causing unaccept-able color fidelity loss in the captured or reconstructed images. Despite the benefitsof facilitating compact form factors, DOEs have sufficient degrees of freedom thatone can manipulate to encode the desirable light modulation functionality.In this dissertation we theoretically and experimentally investigate the practi-cability of introducing numerical optimization into the design procedure of optics,to enable a variety of diffractive optics subject to different application scenarios.Regarding imaging applications, we first validate the practicality of introducingdiffractive-refractive hybrid elements as simplified optics in conventional compu-tational imaging. We re-implement a cross-channel based deconvolution to correctthe chromatic aberration. The full fabrication cycle of photolithography techniqueis developed, serving as the basis for all designs.Then, the desirable focal powers in both spectral and spatial domain are en-coded onto the DOEs. Precisely, we develop the diffractive achromat that balancesthe focusing contributions in full visible spectrum. The color fidelity can thus bewell preserved. Meanwhile, we develop the encoded lenses that provide focus tun-able imaging ability and multi-focal sweep imaging compromise. The trade-offlight loss and residual aberrations are tackled by a deconvolution step.In addition, we extend the design paradigm from image capture to image dis-play, where the holograms of visualization of multiple narrowband spectra are en-iicoded onto a pair of ultrathin diffractive phase plates. The mix-and-match schemeof designing high degree-of-freedom diffractive optics is exploited via the com-plex matrix factorization. The essential random appearance of holograms in factbenefits to enabling a wide range of visualization applications.We envision the work on computational imaging, involving both capture anddisplay end, provide new insights on incorporating optics and computation algo-rithms to better record, understand and deliver desirable visual information underthe constraint of current data bandwidth of hardware.iiiLay SummaryDespite the benefits of facilitating compact form factors, diffractive optical ele-ments (DOEs) have sufficient degrees of freedom that one can manipulate to en-code the desirable light modulation functionality. We theoretically and experimen-tally investigate the practicability of introducing numerical optimization into thedesign procedure of ultrathin diffractive optics.With respect to image capture end, we anticipate enabling more powerful com-putational cameras with more compact form factors. This imaging modality canbe very promising in the consumer-level camera market as well as many scien-tific imaging scenarios. With respect to image display end, we anticipate enablingmix-and-match holographic designs that can drive novel encoding and decodingmodalities in the visualization, education, and security markets.Computational imaging with diffractive optics, involving both capture and dis-play end, provide new insights on incorporating optics and computation algorithmsto better record, understand and deliver desirable visual information under the con-straint of current data bandwidth of hardware.ivPrefaceThe content of this thesis dissertation is based on the following publications (sub-missions) under the framework of Computational Imaging with Diffractive Optics.For each publication (submission), the contributions of each author are listed.- Computational Imaging using Lightweight Diffractive-refractive Optics (15p.)Y Peng, Q Fu, H Amata, S Su, F Heide, W Heidrich - Optics Express (OSA), 2015.W. Heidrich conceived the initial idea of reducing the complexity of optics incamera designs by replacing pure refractive optics with diffractive-refractive op-tics. The implementation of lens design and deconvolution was raised up duringdiscussions between Y. Peng, Q. Fu, and F. Heide. Y. Peng derived the lens de-sign and the optimization algorithm, and wrote an initial paper manuscript. Q. Fufabricated the prototype lenses. Y. Peng, S. Su and Q. Fu performed all experi-ments. F. Heide contributed to evaluating the deconvolution algorithm. All authorscontributed to writing the manuscript.- The Diffractive Achromat: Full Spectrum Computational Imaging with Diffrac-tive Optics (11p.) Y Peng, Q Fu, F Heide, W Heidrich - ACM SIGGRAPH’16Technical Paper (TOG), 2016.W. Heidrich and Y. Peng conceived the initial idea of designing a pure diffrac-tive lens that tackle both the form factor issue as well as the metamarism issue. Themodel derivation for the optical design was raised up during discussions betweenY. Peng, Q. Fu, and W. Heidrich. Y. Peng developed the optimization methodsof both optical design and image deconvolution, and wrote the initial manuscript.Q. Fu fabricated the prototype lenses. Y. Peng and Q. Fu performed all experi-ments. F. Heide contributed to evaluating the deconvolution algorithm. All authorsvcontributed to writing the manuscript.- Encoded Diffractive Optics for Full-spectrum Computational Imaging (11p.)F Heide, Q Fu, Y Peng, W Heidrich - Scientific Reports (NPG), 2016.W. Heidrich and F. Hedide conceived the initial idea of designing a pair ofdiffractive optics to encode different optical functionality subject to different geo-metric alignment. The implementation of lens design and deconvolution was raisedup during discussions between Q. Fu, and F. Heide. F. Hedide developed the op-timization methods of both matrix factorization and image deconvolution. Q. Fufabricated the prototype lenses. F. Heide, Q. Fu, and Y. Peng performed all exper-iments. Y. Peng implemented the simulation and analysis. All authors contributedto writing the manuscript.- Revisiting Cross-channel Information Transfer for Chromatic Aberration Cor-rection (9p.) T Sun, Y Peng, W Heidrich - IEEE Conference on Computer Vision(ICCV) Proceeding, 2017.W. Heidrich and T. Sun conceived the initial idea of tackling the chromaticaberration issue in a form of strong pixel-wise fitting in a blind deconvolutionmanner. The model derivation for the image formation was raised up during discus-sions between T. Sun, and W. Heidrich. T. Sun developed the optimization methodimage deconvolution and wrote the initial manuscript. Y. Peng evaluated the de-convolution algorithm and implemented the simulation comparisons. All authorscontributed to writing the manuscript.- Focal Sweep Imaging with Multi-focal Diffractive Optics (8p.) Y Peng, XDun, Q Sun, F Heide, W Heidrich - To appear at IEEE International Conferenceon Computational Photography (ICCP), 2018.Y. Peng conceived the initial idea of using a focal sweep scheme together withthe assistance of multi-focal diffractive optics to resolve an all-in-focus image com-putationally. The encoded focal range for the optical design was raised up duringdiscussions between Y. Peng and X. Dun. Y. Peng designed the multi-focal lenses,developed the optimization method of image deconvolution, and wrote the initialmanuscript. X. Dun fabricated the prototype lenses. X. Dun and Q. Sun performedall experiments. F. Heide contributed to evaluating the deconvolution algorithm.All authors contributed to writing the manuscript.vi- Mix-and-match Holography (11p.) Y Peng, X Dun, Q Sun, W Heidrich - ACMSIGGRAPH Asia’17 Technical Paper (TOG), 2017.Y. Peng conceived the initial idea of using a pair of diffractive optics to multi-plex multiple holograms subject to different geometric alignment. The multiplex-ing scheme for encoding and decoding the holograms was raised up during discus-sions between Y. Peng and W. Heidrich. Y. Peng developed the modified phaseretrieval algorithm and complex matrix factorization algorithm, and wrote the ini-tial manuscript. X. Dun fabricated the prototype lenses. Y. Peng, X. Dun and Q.Sun performed all experiments. Q. Sun contributed to evaluating the factorizationalgorithm. All authors contributed to writing the manuscript.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Research subject overview . . . . . . . . . . . . . . . . . . . . . 31.2 Noted technical contributions . . . . . . . . . . . . . . . . . . . . 41.3 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Background and Related Work . . . . . . . . . . . . . . . . . . . . . 92.1 Optical aberration and lens design space . . . . . . . . . . . . . . 102.2 Diffractive imaging . . . . . . . . . . . . . . . . . . . . . . . . . 112.2.1 Fresnel diffractive propagation . . . . . . . . . . . . . . . 12viii2.2.2 Broadband diffractive imaging . . . . . . . . . . . . . . . 132.2.3 Computational imaging with flat optics . . . . . . . . . . 142.3 General image restoration model . . . . . . . . . . . . . . . . . . 152.3.1 Iterative optimization schemes . . . . . . . . . . . . . . . 162.4 Natural image deblurring . . . . . . . . . . . . . . . . . . . . . . 192.4.1 Image deconvolution . . . . . . . . . . . . . . . . . . . . 192.4.2 Color fidelity preservation . . . . . . . . . . . . . . . . . 202.5 Matrix factorization . . . . . . . . . . . . . . . . . . . . . . . . . 212.6 Computational image representation . . . . . . . . . . . . . . . . 212.6.1 Refractive caustic designs . . . . . . . . . . . . . . . . . 222.6.2 Diffractive displays . . . . . . . . . . . . . . . . . . . . . 232.6.3 Computer-generated holography . . . . . . . . . . . . . . 232.6.4 Compressive computational displays . . . . . . . . . . . . 232.6.5 Multiplexed holograms . . . . . . . . . . . . . . . . . . . 242.7 Fabrication techniques of flat optics . . . . . . . . . . . . . . . . 253 Lightweight Diffractive-refractive Hybrid Lenses . . . . . . . . . . . 263.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2 Diffractive-refractive optics imaging . . . . . . . . . . . . . . . . 283.2.1 Optical characteristics of Fresnel lenses . . . . . . . . . . 293.2.2 Imaging pipeline . . . . . . . . . . . . . . . . . . . . . . 303.3 Optical designs . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.4 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 333.4.1 Imaging model . . . . . . . . . . . . . . . . . . . . . . . 333.4.2 Optimization method . . . . . . . . . . . . . . . . . . . . 343.4.3 Analysis on efficient optimization . . . . . . . . . . . . . 353.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.5.1 Prototypes . . . . . . . . . . . . . . . . . . . . . . . . . 363.5.2 Experimental results . . . . . . . . . . . . . . . . . . . . 383.5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 403.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44ix4 The Diffractive Achromat . . . . . . . . . . . . . . . . . . . . . . . . 454.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.1.1 Diffractive imaging model . . . . . . . . . . . . . . . . . 474.2 Optical design optimization . . . . . . . . . . . . . . . . . . . . . 504.2.1 Optimization model . . . . . . . . . . . . . . . . . . . . 504.2.2 Stochastic optimization algorithm . . . . . . . . . . . . . 534.3 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 544.3.1 Optimization method . . . . . . . . . . . . . . . . . . . . 554.3.2 Efficiency and robustness analysis . . . . . . . . . . . . . 584.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.5 Evaluation and discussion . . . . . . . . . . . . . . . . . . . . . . 624.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675 Multi-focal Sweep Imaging . . . . . . . . . . . . . . . . . . . . . . . 685.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.2 Optical design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.2.1 Multi-focal diffractive lens . . . . . . . . . . . . . . . . . 705.2.2 PSF behaviors . . . . . . . . . . . . . . . . . . . . . . . 725.3 Image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 735.3.1 Sweeping image formation . . . . . . . . . . . . . . . . . 735.3.2 Optimization method . . . . . . . . . . . . . . . . . . . . 745.3.3 Efficiency and robustness analysis . . . . . . . . . . . . . 755.4 Implementation and discussion . . . . . . . . . . . . . . . . . . . 765.4.1 Prototype parameters . . . . . . . . . . . . . . . . . . . . 765.4.2 Simulation results . . . . . . . . . . . . . . . . . . . . . 775.4.3 Real world results . . . . . . . . . . . . . . . . . . . . . 775.4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 805.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816 Encoded Diffractive Lenses . . . . . . . . . . . . . . . . . . . . . . . 826.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 826.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 856.2.1 Encoded diffractive lens design . . . . . . . . . . . . . . 85x6.2.2 Reconstruction algorithm . . . . . . . . . . . . . . . . . . 876.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 896.3.1 Encoded diffractive lens . . . . . . . . . . . . . . . . . . 906.3.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . 916.3.3 Static colour imaging . . . . . . . . . . . . . . . . . . . . 936.3.4 Rotational refocusing . . . . . . . . . . . . . . . . . . . . 946.3.5 Rotation-only zooming . . . . . . . . . . . . . . . . . . . 946.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 987 Mix-and-match Holography . . . . . . . . . . . . . . . . . . . . . . 1007.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1017.2 Computational design . . . . . . . . . . . . . . . . . . . . . . . . 1027.2.1 Diffractive image formation . . . . . . . . . . . . . . . . 1027.2.2 Hologram optimization . . . . . . . . . . . . . . . . . . . 1047.2.3 Mix-and-match multiplexing . . . . . . . . . . . . . . . . 1077.3 Applications and results . . . . . . . . . . . . . . . . . . . . . . . 1127.3.1 Application scenarios . . . . . . . . . . . . . . . . . . . . 1137.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1197.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1248 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129A Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146A.1 Supplementary details of fabrication . . . . . . . . . . . . . . . . 146A.1.1 Specifying a Fresnel lens. . . . . . . . . . . . . . . . . . 146A.1.2 Fabrication procedure illustration. . . . . . . . . . . . . . 147A.2 Supplementary details of image reconstruction algorithms . . . . . 149A.2.1 Blind deconvolution. . . . . . . . . . . . . . . . . . . . . 150A.2.2 Non-blind deconvolution. . . . . . . . . . . . . . . . . . 152xiList of Tables4.1 Quantitative evaluation of state-of-the-art algorithms. . . . . . 595.1 Quantitative evaluation of different image reconstructions. . . . 767.1 Quantitative evaluation under different experimental conditions. 120xiiList of Figures1.1 Overview of Computational Imaging with Diffractive Optics. . 52.1 Illustration of first-order lens aberrations. . . . . . . . . . . . 113.1 Imaging pipeline of diffractive-refractive computational imaging. 283.2 Illustration of multi-order diffraction. . . . . . . . . . . . . . 303.3 Comparison of diffractive optics and refractive optics. . . . . 313.4 Illustration of prototyped lenses and their PSFs. . . . . . . . . 333.5 Imaging quality comparison of five lenses. . . . . . . . . . . . 373.6 Experimental results of the stacked diffractive optics. . . . . . 393.7 Experimental results of the hybrid lens. . . . . . . . . . . . . 403.8 Experimental results of the hybrid phase-coded-aperture lens. 413.9 Comparison results of capturing a resolution chart of three lenses. 423.10 Comparison results of the off-axis performance of three lenses. 433.11 Experimental results of an outdoor and large depth range scene. 444.1 Illustration of initializing the height profile for optimization. . 534.2 Cross section comparison of normalized PSFs of three lenses. 554.3 Deconvolution comparison of state-of-the-art algorithms. . . . 584.4 Synthetic results of our diffractive achromat. . . . . . . . . . 614.5 Experimental results of our diffractive achromat (Scenes I). . . 614.6 Experimental results of our diffractive achromat (Scenes II). . 624.7 Experimental results capturing a standard resolution chart. . . 634.8 Comparison of on-axis and off-axis behavior of two lenses. . . 644.9 Experimental comparison of color performance of two lenses. 65xiii4.10 Measured spectral PSFs of our diffractive achromat. . . . . . 665.1 Focused distances subject to object distance and focal length. . 715.2 Comparison of synthetic PSF behavior of two sweeping lenses. 725.3 Graph fusion and microscope images of prototyped lenses. . . 735.4 Photograph of the experimental setup. . . . . . . . . . . . . . 775.5 Illustration of selected simulation results. . . . . . . . . . . . 785.6 Experimental results of our multi-focal sweep imaging (Lens I). 795.7 Experimental results of our multi-focal sweep imaging (Lens II). 795.8 DOF comparison of a multi-focal lens and a standard lens. . . 806.1 Overview of the encoded lenses computational imaging. . . . 846.2 Encoding optimization examples for different target functions. 916.3 Illustration of experimental set-up. . . . . . . . . . . . . . . . 926.4 Results of static color imaging with a encodable lens. . . . . . 956.5 Results of capturing a resolution chart. . . . . . . . . . . . . . 966.6 Results of refocusing with relative rotation between DOEs. . . 976.7 Results of rotation-only zooming with encoded lenses. . . . . 987.1 Processing pipeline of mix-and-match holography. . . . . . . 1037.2 Flowchart of our modified double-sampling GS algorithm. . . 1057.3 Configuration of diffractive wave propagation. . . . . . . . . 1067.4 Illustration of the matrix structure for our multiplexing scheme. 1077.5 Synthetic results of the sub-cell offset design. . . . . . . . . . 1147.6 Results of projecting an image through an offset design. . . . 1157.7 Results of projecting images through multiple pairing designs. 1167.8 Experimental demonstration of the fast alignment manipulation. 1177.9 Results of projecting an animation through translation pairings. 1187.10 Diagram and results of projecting multiview images. . . . . . 1197.11 Decoded image visualizations of simulating misalignments. . 1217.12 Comparison of convergence for two factorization methods. . . 1237.13 Results of broadband imaging of an offset pairing example. . . 124A.1 Illustration of photolithography fabrication procedure. . . . . 148xivA.2 Illustration of microscope images of fabricated lenses (I). . . . 148A.3 Illustration of microscope images of fabricated lenses (II). . . 149xvGlossaryBM3D block matching and 3D filteringBRDF bidirectional reflectance distribution functionCGH computer generated hologramCNC computer numeric controlDOE diffractive optical elementDOF depth of fieldDSLR digital single-lens reflexFOV field of viewFT Fourier transformHDR high dynamic rangePSF point spread functionPSNR peak signal-to-noise ratioRIE reactive ion etchSSIM structural similarity indexSLM spatial light modulatorTV total variationUV ultravioletxviAcknowledgementsFirst and foremost, I would like to express the thanks to my supervisor WolfgangHeidrich for his guidance and support throughout my time at The University ofBritish Columbia. Thanks to my thesis and RPE committee members James Little,Mark Schmidt and Robert Woodham for their timely feedback and advice over theyears.Thanks to UBC fellows Lei Xiao, Gerwin Damberg, James Gregson, and RobinSwanson. Thanks to all my research collaborators: Felix Heide, Gordon Wetzstein,Shuochen Su, Qiang Fu, Xiong Dun, Qilin Sun, Tiancheng Sun, Hadi Amata, RuiWang, Chao Yu, Chen Su, Qing Zhong, Haifeng Li, Liang Xu, Xu Liu, for all theconstructive discussions. Thanks to Tina Smith and Laura Selander for administra-tive help.I would also like to thank the following institutes, Imager Lab at UBC, Compu-tational Imaging Group at Stanford University, Visual Computing Center at KAUST,Advanced Nanofabrication Imaging and Characterization Core Lab Lab at KAUST,State Key Lab of Modern Optical Instrumentation at Zhejiang University, State KeyLab of CAD & CG at Zhejiang University, for providing fruitful support.Special thanks are owed to my parents, whose have supported me throughoutmy years of education, both morally and financially.xviiDedicationTo my parents, for their unsparing love and support.xviiiChapter 1IntroductionLet there be light to illuminate the world. Human beings and many mammalianspecies acquire most of external natural information through the vision systems.With decades of evolution, modern imaging technology has evolved to equip hu-man beings with remarkable capabilities to better record, understand, visualize,and analyze the world. The imaging field actually represents a procedure of infor-mation recording and reconstruction from capture to display, whose rise is mostlyattributed to the advances in optics and semiconductor electronics research. No-tably, more precise optics and sensors contribute to an increasing supply for highquality imaging solutions.In conventional imaging systems, for both capture and display, increasing opti-cal complexity is inevitable for either correcting a wide variety of aberrations [143]or enabling additional optical functionalities [87]. Simultaneously, high qualityimaging with reduced optical complexity is an ongoing target for investigationfrom both academic and industrial research.The power of computing nowadays refreshes our perspective on many conven-tional imaging problems. Notably, recent advances in computational imaging haveintroduced computation as a virtual component in optical imaging procedures todrastically shift the burden from optics to algorithms. This is achieved throughmathematical formulations of image restoration as inverse problems that can beefficiently and robustly solved by manipulating optimization algorithms of somesense. The trends in imaging development allow for significantly reduced opti-1cal complexity while maintaining high image fidelity and realistic light efficiency(e.g., [57, 131]). This is a core motivation for this thesis work.In the visual computing community, point spread function (PSF) engineering isa well-known technique that has proven useful for applications such as extendeddepth-of-field [30], motion [125] and defocus [180] deblurring, and 3D or ligthfield imaging [88, 97, 156]. Recent advances in computational imaging hold lim-ited pre-conditioning in the design space, especially very limited with consideringthe mutual effect of the joint design between optics and computation.For decades, diffractive optical element (DOE)s have drawn considerable atten-tion in the optics community because of their ultrathin and lightweight physicalstructure, unprecedented large and flexible design space, availability of state-of-the-art fabrication techniques, as well as relatively spatial invariant focusing be-havior. During the past few years, integrating diffractive imaging elements andcomputational methods in a single imaging system has resulted in several newcomputational imaging devices with ultra compactness (e.g., [47, 117, 147, 159]).In this thesis work, we leverage phase modulation DOE as a lens platform, aroundwhich we build novel computational imaging modalities. We acknowledge that ourdesign intentions differ from lensless imaging in that we are not aiming for entirelyremoving imaging optics to create flat cameras, an approach that will be analyzedin the next chapter.It is worth noting two main problems are brought to attention for diffractiveimaging applied in the full visible spectrum. First, the wavelength-dependencyof diffraction leads to strong chromatic aberrations that degrade the image qualitywith blurs of a very large diameter. Second, the strong chromatic aberrations,in addition, cause significant wavelength dependence of the point spread func-tion (PSF) even within a single color channel. In particular, this wavelength de-pendency means that objects with the same RGB color are blurred differently if theunderlying spectral distributions differ [165].This metamerism problem means that image restoration algorithms in prac-tice use approximate PSFs based on some fixed spectral distribution that is im-plicitly derived from the lighting conditions during calibration. We note that thisapproximation can result in a significant loss of color fidelity even after imagereconstruction [117]. While metamerism affects all imaging systems exhibiting2chromatic aberration including regular refractive optics, the problem is particu-larly pronounced for the strong dispersion present in diffractive optics. That said,when designing an imaging approach with diffractive optics, we need to carefullytake consideration of both spatial and spectral focal power [118]. Effectively andefficiently compromising these two kinds of focal powers subject to desirable ap-plication scenarios is one of the core goals throughout this thesis work.On the other hand, when applied to display end, this chromatic variance inmodulating light rays limits general holographic solutions to only be viable in sin-gle color. A conventional approach to deliver color visualization relies on spatialor temporal multiplexing [11], that practically introduces extra hardware cost andengineering effort. Despite the spectral variance, the great design freedom of DOEis beneficial for enabling a wide range of holographic visualizations. This insightinspires us to exploit a mix-and-match holographic visualization modality to bridgethe gap between conventional goal-driven caustics or computational light steeringand holographic projections with small form factor phase holograms.From a broad perspective of the other key aspect in computational imaging,image reconstruction algorithms, significant progress has been made recently insolving inverse problems, as is reviewed in a comprehensive survey report [102].As the rise of mobile imaging systems such as smartphone modules, developingmore powerful and more efficient optimization techniques are highly desirable.This is another point of investigation for this thesis dissertation. State-of-the-artmethods can be grouped into two classes: prior-based optimization and machinelearning (data driven) approaches [56]. Both have guided successful applications inrecent years. At the current stage of the proof-of-concept, in this thesis dissertationwe mainly exploit prior-based methods for image reconstruction.1.1 Research subject overviewThis thesis dissertation presents our research on the above topics by introducingnot only improvements in the image reconstruction (post-processing), but moreimportantly, optimizing the diffractive optics itself subject to application intentions.Therefore, a principal argument that covers almost all aspects throughout this workis joint design, including the optimization in optical design and the optimization in3reconstruction process.Figure 1.1 presents the subjects and the methods investigated in this thesis dis-sertation. Falling under the scope of Computational Imaging with Diffractive Op-tics, the research subjects can be cast into two directions of imaging — captureand display. From any perspective, the technical methodologies exploited includediffractive optical design, numerical optimization algorithm, and prototyping ex-perimentation.We imagine a wide range of applications for the proposed computational imag-ing modality. Under each branch, we present research projects that are built uponthese three ideas. With respect to the image capture end, we anticipate enablingmore powerful computational cameras with more compact form factors. This imag-ing modality can be very promising in the consumer-level camera market as wellas many scientific imaging scenarios. With respect to the image display end, weanticipate enabling mix-and-match holographic designs that can drive novel encod-ing and decoding modalities in the visualization, education, and security markets.Refer to individual chapter and corresponding publication for details.1.2 Noted technical contributionsThe technical contributions of this thesis dissertation are as follows:We investigate the principle and the implementation practicability of applyingnumerical optimization in the joint design of optics and computation, to enablea variety of imaging applications with compact form factor optics. To our bestknowledge, we are one of the very first attempts in this community that has inves-tigated this interdisciplinary direction.We derive both the diffractive image formation model and the inverse prob-lem model that can be generalized to alternative broadband diffractive imagingsolutions. We start by validating the direct use of diffractive-refractive optics inbroadband computational imaging to lay a foundation of prototyping. We nextmove to design the diffractive achromat, a pure diffractive imaging system for fullspectrum visible light that combines optimization in both diffractive optical designand post-capture image reconstruction.By tuning the functions and the terms in the model subject to a variety of light4Figure 1.1: Overview of research subjects and methods investigated for Com-putational Imaging with Diffractive Optics. Our work is motivated bythe essential benefits of diffractive optics and modern high computingpower (top-left). The joint design research scope can be cast into twobranches — capture and display, respectively. Each end includes a seriesof research projects and applications that correspond to the Chapters (3-7). The methodologies of diffractive optical design, optimization-basedimage reconstruction, and engineering implementation bridge the gapamong different application scenarios throughout this thesis disserta-tion. Specially, the fabrication serves as a platform for all investigatingapplications (bottom-left).5bending requirements, we further exploit the practicality of encoding in spatialdomain. Precisely, we enable a multi-focal sweep imaging modality with the ex-tended depth-of-field and the focus tunable lenses with manipulating different ge-ometrical alignment of a pair of phase plates.This work covers the computational challenges of computational imaging, thatare represented as exploiting formal optimizations in both optical design and im-age reconstruction. Precisely, we exploit the following optimization algorithms:deconvolution with different image priors, matrix factorization in complex valuedomain, modified iterative phase retrieval, and heuristic random search with phys-ical constraints.On the one hand, we employ a non-convex but effective optimization for de-signing achromatic DOEs subject to diverse spectra requirements, which re-arrangesthe spatial and spectral distributions of PSFs so that chromatic aberrations and colorcorruption are largely eliminated in hardware.In addition, we expand the regular non-negative matrix factorization in opti-mization to the complex value domain, and apply it to extract the phase distribu-tions of individual layers of both encoded lenses and multiplexed holograms. Amodified projection rule is specially exploited to enforce better and faster conver-gence under our application scenarios.On the other hand, the inverse problems to solve for most of the applications arenon-convex or convex optimization with sparse representation. We re-implementthe cross-channel prior in the deconvolution with modern optimization trick andsolver. This can be directly applied in most broadband diffractive imaging solu-tions, to remove chromatic aberration digitally. We also exploit a cross-scale priorto further mitigate residual aberrations introduced by the designed diffractive op-tics. With the optimized PSFs that exhibiting spectral invariance and preservingcentral sharp peak, this prior exhibits more robustness and higher efficiency thanstate-of-the-arts.1.3 OrganizationThe remaining part of this thesis dissertation is structured as follows.Chapter 2: Background and Related Work. This chapter reviews the working6principles of diffractive optics and computational imaging. Precisely, we surveystate-of-the-art work that is closely relevant to the investigated framework fromboth the perspective of optics and that of computation. Besides, the available fab-rication techniques of are overviewed. Thus, this chapter lays the foundation ofderiving and understanding the techniques described in the chapters following.Chapter 3: Computational Imaging with Lightweight Diffractive-refractive Op-tics. This chapter summarizes the first attempt of introducing diffractive optics intoconsumer level photography solutions, together with the validation of photolithog-raphy fabrication technique, serving as a proof-of-concept project under the high-level computational imaging framework.Chapter 4: The Diffractive Achromat. This chapter details a pure diffractiveoptics camera that incorporates optimization both in optical design and in image re-construction. This is realized by engineering desirable point spread function (PSF)behaviours along full visible spectrum, further leading to obtaining high-qualitycomputational imaging in full color. The optical design step is formed as an op-timization problem of updating the phase profile of an ahcromat, meanwhile theimage reconstruction step is formed as a deconvolution problem of removing theresidual achromatitc blur.Chapter 5: Multi-focal Sweep Imaging with Extended DOF. This chapter ex-tends from the projects presented in the previous two chapters, and encodes multi-ple focal powers subject to depths and spectra onto diffractive optics. We developan extended depth of field (DOF) imaging with multi-focal sweep imaging modal-ity. This is realized by fusing multiple focal powers onto one single element soas to obtain approximated depth-invariant kernel behavior after sweeping. The re-maining residual blur is tackled by the deconvolution algorithms.Chapter 6: Encoded Diffractive Optics for Color Imaging. This chapter ad-dresses another encoding scheme in spatial domain, and validates an efficient de-sign manner to encode multiple optical functionalities onto diffractive optics, in-cluding refocusing and zooming effect. In particular, we develop focus-tunablelenses with a pair of diffractive optics. Despite the optimization (i.e. deconvo-lution) in the image reconstruction step to remove aberrations, a complex matrixfactorization scheme is specially introduced in deriving the pattern of each phaselayer.7Chapter 7: Mix-and-match Holography. This chapter extends from the par-tial optimization scheme described in the previous chapter but fits its practica-bility on display end. Precisely, the combination of controlled diffractive prop-agation and geometrical aligned holograms enables a mix-and-match holographytechnique that can deliver multiple different target images. The modified phaseretrieval methods and the complex matrix factorization are exploited.Chapter 8: Conclusion. This chapter concludes the technical concepts andimplications that are implemented in this thesis dissertation. The features of sucha diffractive imaging system are outlined, major challenges on the way toward thisvision are discussed, and the next immediate steps explored.Chapter A: Appendix. This chapter includes the supporting materials includingfabrication details of photolithography technique and image deconvolution deriva-tions.8Chapter 2Background and Related WorkIn this chapter we review the background techniques and the recent work related tousing diffractive phase modulation to enable desirable imaging and display func-tionalities. Until our recent work, diffractive optics has almost exclusively beenused for single wavelength or narrow band imaging in scientific application scenar-ios. In addition, almost all existing optical designs are based on analytic solutionsand heuristic design manners.We here firstly explain the fundamental image degradation sources and thecommon technical solutions in optics (Sec. 2.1). In this dissertation we describeimaging systems relying on wave behavior of light, more specifically the scalardiffraction theory (Sec. 2.2). We then review the advantages of introducing com-putation schemes into the joint design of optics and image reconstruction (Sec. 2.3and Sec. 2.4), in order to resolve high quality imaging with broadband illumination.In these two parts, an amount of relevant work in image deblurring is analyzed. Inaddition, the current research patterns and the development tendency lead to thepracticality of extending our design framework from capture end to display end(Sec. 2.6). This includes intensive discussion on both refractive caustic designsand holographic approaches. Last, we briefly discuss state-of-the-art fabricationtechniques of flat optics (Sec. 2.7). More detailed discussion on related techniquesare to be presented in the corresponding chapters.92.1 Optical aberration and lens design spaceOptical systems suffer from both monochromatic and chromatic aberrations [38],which are derived from the optical path deviations when light travels through dif-ferent regions of a real lens (i.e. having a certain of aperture size and thickness)from different incident angles. That said, optical aberrations accurately describethe errors introduced by using the paraxial and thin-lens approximations in the de-sign space.Monochromatic aberrations describe the deviations from the idealized linearGaussian optics [40] that occur even for a single wavelength (or narrow spectra).Five primary aberrations are classified into this type, including spherical aberra-tion, coma, astigmatism, field curvature, and distortion. The aberrations of a reallens are illustrated in Fig. 2.1, interpreted as a fact that the beam emitted from apoint source may be focused at many points that vary in the spatial domain. Chro-matic aberrations result from the wavelength-dependency of the refractive indexof the material. Consider the scenario of imaging a white point light source, theshift or spread of the imaged spot along the optical axis is denoted as axial chro-matic aberration, while the remaining component in the image plane is denoted aslateral chromatic aberration [76]. Chromatic aberrations will be discussed in moredetail in later chapters, which present one of our key contributions for preserv-ing color fidelity in post-capture processing. All these aberrations, including bothmonochromatic and chromatic aberrations, are visualized as unwanted spatial andspectral variant blur of the image, that may become more severe for a lens havinga relatively large numerical aperture (aka. a small f -number).The well-known aberration theory was firstly introduced by Seidel in 1857 [135].Ever since then lens design in optics attempts to eliminate aberrations of all kindsby designing increasingly complex lens structures that consist of a large numberof elements [143]. This involves designing aspherical surfaces as well as findingmaterials with better optical properties. Therefore, lens design is actually a heuris-tic compromise among a variety of optical evaluation criteria [41, 94]. Accord-ingly, state-of-the-art commercial optical design softwares, like Zemax, Code-V,etc., have been intensively used in designing refractive lens profiles for a widerange of applications. Intuitively, a high level insight in the design space is to en-10Figure 2.1: Pinhole model (top-left), monochromatic aberrations (centerfive), and chromatic aberrations (bottom two) of a real lens.force the focusing contribution to be as sharp as possible subject to field-of-viewsand wavelengths. Only until the recent decade, researchers have recognized thatoptical aberrations can not only be deliberately minimized but designed to exhibitadditional behavior for obtaining extra information of scene. Computational imag-ing then emerged, which is to be reviewed in Sec. 2.4.2.2 Diffractive imagingDespite using refractive optics as the imaging or projection lens, diffractive opticalelement (DOE) is very flexible in modulating light, thereby delivering greater free-dom in design parameters [117, 158]. DOEs have been investigated intensively inthe optics community, for applications such as extending depth-of-field in micro-scopes [155], or realizing an optimal phase-only filter attached with lenses to gainexpected optical transfer function [33]. DOEs are operated by means of interfer-ence and diffraction to produce arbitrary distributions of light [9, 112]. We addressthe following inherent advantages. First, DOE can be fabricated on a thin sheet;11second, a single DOE can perform multiple optical operations simultaneously toact as an efficient light modulation platform; third, it suffers from less off-axisaberrations due to the flat lens profile.2.2.1 Fresnel diffractive propagationWe start by reviewing the propagation of light from the DOE plane to the imageplane. To this end, we apply scalar diffraction theory [52] under the assumption ofparaxial incident light. Using the Fresnel approximation, the light propagation isgiven in Eq. 2.1U˜(x, y) =exp(jkz1)jλ z1∫ ∫∑U˜(x1, y1)exp{jk2z1[(x− x1)2+(y− y1)2]}dx1dy1,(2.1)where (x,y) and (x1,y1) refer to coordinates on the image plane and the DOE plane,respectively. j is the imaginary unit, k is the wave number, λ is the wavelength, andz1 stands for the axial distance from the DOE plane to the image plane. Accord-ing to Eq. 2.1, one can calculate the integral over the whole DOE transmittanceaperture ∑ to compute the complex amplitude distribution U˜(x,y), including bothamplitude and phase information.We observe that the integral range ∑ in Eq. 2.1 can be released to the infiniteplane of (x1,y1), considering a fact that the complex amplitude U˜(x1,y1) outsidethe aperture is actually zero. We next mathematically expand the terms of theexponential part inside the integral, and extract those terms that do not include(x1,y1) outside the integral. Further by partially applying the definition k = 2piλ , thediffractive image formation is rewritten as in Eq. 2.2:U˜(x, y) =exp(jkz1)jλ z1exp[jk2z1(x2+ y2)]∫ ∫ ∞−∞U˜(x1, y1) · exp[jk2z1(x12+ y12)]·exp[−j2pi(x1 xλ z1 + y1yλ z1)]dx1dy1.(2.2)Specifying transmittance function. We next revisit the first half of propagation,that the wave is propagated from the light source and then modulated by the DOE.12Particularly, U˜(x1,y1) is the wave exiting the DOE plane, i.e. after diffraction inthe DOE has already occurred. It is related to the incident wave U˜0 at the DOE viathe complex transmittance function t(x1,y1):U˜(x1,y1) = U˜0(x1,y1) · t(x1,y1) (2.3)witht(x1,y1) = A(x1,y1) · exp(jΦ(x1,y1)), (2.4)where A represents the amplitude and Φ represents the phase of each point on theDOE plane.Introducing the lens phase. We next consider imaging systems where we repre-sent the complex field of a spherical wave (convex) as follows:U˜s(x1, y1) =Asz1exp(−jkz1)exp[− jk2z1(x12+ y12)]. (2.5)We recognize Eq. 2.2 to the Fourier transform (FT) of the product of the com-plex field just to the right of the aperture and a quadratic-phase exponential. Wecan define u = xλ z1 and v =yλ z1 , which is practically implementable. Note that atthe focal plane of a lens, the diffractive propagation in fact becomes Fraunhoferdiffraction. By incorporating Eq. 2.5, the diffraction is described as:U˜(x, y) =1jλ z12exp[jk2z1(x2+ y2)]F[U˜0(x1,y1) · t(x1,y1)], (2.6)whereF [·] denotes the Fourier transformation (FT). Our diffractive image forma-tions in all the following projects use the above model, with the details explainedin each individual chapter.2.2.2 Broadband diffractive imagingDespite the advantages of ultrathin and lightweight physical appearance, as wellas large design freedom, the severe color dispersion in diffraction has limited theirapplications in imaging under broadband illumination. This can be noted from13Eq. 2.2. In practice, a limited amount of work has applied diffractive optical ele-ment (DOE) in consumer level imaging products, but only in collaboration with re-fractive lenses. For example, a multilayer DOE has been used to correct chromaticaberrations in refractive lens systems [107], although it still relies on a multitudeof refractive lens elements to contribute to most of the focal power.Recent advances on achromatic metasurfaces has revealed the potential foruse in lightweight collimators or in microscopy imaging [3, 74]. Unlike con-ventional diffractive imaging that has been investigated for years, metasurfacesmodulate the behaviors of electromagnetic waves through the specific boundaryconditions, which can also be recognized as an extreme case of diffraction. Thevery recent emerged metalenses indeed investigated a field of fabricating elementsin the feature scale of sub-wavelength. They are promising in shaping wave-front while it remains challenging to create realistic lenses for now. As such, achromatic-aberration-corrected diffractive lens for broadband focusing has beendesigned [159, 160]. We note that these reported works don’t address the recon-struction problem and they are only designed at three wavelengths rather than themore practical broadband usage.A huge amount of research has been done in designing DOEs for multiplewavelengths, either by using harmonic diffraction and multiple layers to createhigh efficiencies in multiple wavebands [140], or by introducing computation tothe design footprint to redefine the light transmittance function. The latter case isof significance not only to mitigate the dispersion, but also to encode the engineeredphase distribution on the DOE profile for the tasks of image recovery or enablingdesirable optical functionalities [124].2.2.3 Computational imaging with flat opticsIn computational imaging, two types of devices based on diffractive optical element(DOE) have recently been raised — lensless computational sensors [47, 105] andFresnel lens imaging with post processing [110, 117]. The former two integrateDOEs into the sensor, aiming for an ultra-miniature structure with medium imagequality. The latter two share a similar concept to apply Fresnel lens to replace bulkyand heavy refractive lenses in a camera design. Although the chromatic aberrations14can be partially mitigated by optimizing for three discrete wavelengths or digitallyremoved in a post-capture step, the color fidelity may be significantly reduced infinal images due to metamerism (recall the discussion in Chapter 1). That said,naively relying on post-processing to remove chromatic aberration digitally couldbe ill-conditioned. We will expand the analysis on the color fidelity issue in laterchapters.Another promising solution of enabling flat optics is to place a coded apertureor a diffuser in front of the sensor. Very recently, Antipa et al. [4] propose a single-shot lensless solution that uses a diffuser to enable compressive sensing ability andthen use the designed algorithms to resolve the 3D imaging from the compressedcaptured caustic patterns. Regarding this kind of lensless schemes, the limited fieldof view (FOV) and resolved image resolution remain problematic. In this work, weare not intended to completely get rid of lenses, but to provide flexible modulationon the diffractive lenses such that a realistic image resolution with reasonable FOVcan be achieved.2.3 General image restoration modelBroadly speaking, state-of-the-art image restoration methods involve removingblur, noise, over-exposure pixels, and etc. While these different kinds of imagedegradation sources are usually corrupted across the image thereby meaningless tobe separately tackled, in this thesis work blurry degradation is the major problemwe aim to resolve when investigating diffractive imaging.In information theory, the defocus effect can be formulated as a latent imageconvolved with a blur kernel. Then, we can write the recorded image in channel cin vector-matrix format as in Eq. 2.7:bc = Kcic+n, (2.7)where for a channel c, b, K, i, n are the captured image, convolution matrix, sharplatent image and additive noise in the capture, respectively.Regarding an imaging system, Kc can be derived from the PSF Pc integrated15over the spectrum Λ as in Eq. 2.8:Pc(x,y,z) =∫ΛQc(λ ) ·P(x,y,z;λ )dλ , (2.8)where P(x,y,z;λ ) is the spatial spectral variant PSF describing the aberrations oflens, which is a function of both spatial position (x,y,z) and spectral component λ .Qc represents the sensor response.Image restoration tasks are typically formulated as inverse problems. As theoriginal problems are typically ill-posed, prior knowledge of the latent images areusually added as regularization. Thus, the image restoration seeks to solve a mini-mization problem, with its general formation given in Eq. 2.9:ic = argminicµc2‖bc−Kcic‖22+Γ(ic), (2.9)where the first term in the objective is a standard least-squares data fitting termthat is derived from likelihood, µc is the weight of data fitting term for channelc = 1,2,3. The second term Γ(ic) is a regularization term that enforces naturalimage priors in the solution.From the perspective of statistics, solving restoration tasks as maximum a pos-teriori probability (MAP) estimation problems not only requires physically moti-vated models for the data likelihood, but relies on effective prior knowledge onthe latent images. For instance, the well-known local statistical prior, total varia-tion (TV) enforces sparsity and preserves high frequency edges [50, 128]. In ad-dition, the highly successful block matching and 3D filtering (BM3D) method [26,83] searches for similar patches within the image and combines the informationthrough a collaborative filtering step. Both deblurring and denoising can then betackled by exploiting good image priors. This is a key research topic that we willspecially review in the next section. Now, we may engineer different forms of theregularization term subject to different image restoration tasks.2.3.1 Iterative optimization schemesNumerical optimization has become a standard tool for efficiently solving a num-ber of classical image restoration problems with the priors [13]. Recent years the16graphics and computer vision community have witnessed the widespread usage ofmany representative iterative optimization algorithms such as Half Quadratic Split-ting (HQS) [42], Split Bregman [51], Alternating Direction Method of Multipliers(ADMM) [12], First Order Primal Dual method (PD) [18], etc. These formal opti-mization methods significantly improve the performance of image restoration bothin time and quality.Note that most of the above well-known algorithms are recognized as proxi-mal algorithms. The concept behind this is to evaluate the proximal operator ofa function and involves a small optimization problem [115]. That said, they mayhave a closed-form solution under some certain assumption so as to be well-suitedto large-scale problems. In this thesis work, we widely investigate the two well-defined proximal algorithms (aka. HQS, ADMM) in deconvolving latent imagesfrom the degraded ones.Proximal operators. The basic operation is evaluating the proximal operator of afunction. Thus, given a function f , we can define its proximal operator as:proxf/λ(v) = argminx f (x)+λ2‖v−x‖22. (2.10)proxf/λ(v) can be interpreted as a point that compromises between minimizingthe function f and being close to point x. The positive scalar λ is the trade-offparameter. More details of proximal operators can be found in [115]. Below wepresent several representative functions and corresponding proximal operators.When f is the l0 norm, its proximal operator is called hard-shrinkage:proxf/λ(v) ={0,−1/λ ≤ v≤ 1/λv, otherwise(2.11)When f is the l1 norm, its proximal operator is called soft-shrinkage:proxf/λ(v) =v+1/λ , v <−1/λ0,−1/λ ≤ v≤ 1/λv−1/λ , v > 1/λ(2.12)17Alternating direction method of multipliers (ADMM) solver. Without loss of gen-erality, we here take the widely used TV prior as an instance. As discussed in theprevious chapter, the objective function to minimize can be formulated as:x = argminxλ2‖b−Ax‖22+‖∇x‖1. (2.13)Note that the above objective is non-linear and can be hard to solve with con-ventional optimization tricks. The concept is to re-formulate the objective byadding different slack variables.Alternating direction method of multipliers (ADMM) [12, 157] relaxes theoriginal objective in Eq. 2.13 to be:x = argminx,zλ2‖b−Ax‖22+‖z‖1+ρ‖z−∇x+u‖22, (2.14)where slack variables z and u are introduced. Eq. 2.14 is iteratively minimized bysolving for the latent image x, slack variable z and u alternately, that actually opti-mizes in the dual domain with a gradient descent scheme, as given in the followingequations:xt = argminx,z,uλ‖b−Ax‖22+ρ‖zt−1−∇x+ut−1‖22,zt = argminzλ‖z‖1+ρ‖z−∇xt +ut−1‖22,ut = ut−1+ zt −∇xt .(2.15)Half quadratic splitting (HQS) solver. Half quadratic splitting (HQS) [42] relaxesthe original objective in Eq. 2.13 to be:x = argminx,zλ2‖b−Ax‖22+‖z‖1+ρ‖z−∇x‖22, (2.16)where a slack variable z is introduced to approximate x, and ρ is a positive scalar.Eq. 2.16 is iteratively minimized by solving the latent image x and the slack vari-18able z alternately, as given in the following equations:xt = argminxλ‖b−Ax‖22+ρ t‖zt−1−∇x‖22,zt = argminzλ‖z‖1+ρ t‖z−∇xt‖22.(2.17)Importantly, ρ t shall increase as the iteration t continues. The latter (aka. largerρ) forces z to become an increasingly good approximation of x, that we will furtherdiscuss in Chapter 4.As detailed previously, solving the above sub-problems of both ADMM andHQS relies on direct Fourier and Inverse Fourier transformation, and soft-shrinkageoperation. In addition, we have an observation that HQS method typically requiresfewer iterations to converge in practice than ADMM methods.2.4 Natural image deblurringIt is clear that optical aberrations can be designed for easy computational removalin a post-capture step, for instance, by favoring PSFs that preserve high spatialfrequencies. This combination of optics and computation is one of the major de-velopment tendencies in photography industry. Different criteria, e.g., design foot-print, light and computational efficiency, are to be traded off [34, 174]. In thisdissertation we mainly apply image deconvolution in our image restorations.2.4.1 Image deconvolutionBack to earlier time, the classical Richardson-Lucy iterative method for image de-convolution [34, 93, 126] was developed for use in astronomy to compensate forblur in optics. This famous algorithm has been investigated and extended in nu-merous ways, including residual deconvolution [173] and Bilateral Richardson-Lucy [174]. Recent advances in image deconvolution seek to include local statisti-cal prior knowledge, like using the well-known total variation (TV) prior [19, 128],which promotes piece-wise constant images to restore high-frequency informationsuch as edges. Cho et al. [20] proposed to match the gradient distribution for im-age restoration. Alternative non-convex regularization terms are intensively inves-tigated [79, 87], empirically giving improved results at reasonable local optimum.19Non-blind image deconvolution techniques mostly assume the aberration-inducedblur kernel is known [70] or can be easily calibrated. Despite adding generic priorsin the optimization, learning-based methods like filters learning [166] and convo-lutional neural network [133] have been recently reported.We note that a known convolution kernel can be impractical to calibrate, espe-cially where a large field-of-view (FOV) imaging is expected. Accordingly, blinddeconvolution algorithms aim to reconstruct the latent image and estimate the ker-nel simultaneously [89, 136, 175]. Most blind deconvolution algorithms use an ex-pected maximization-type manner to blindly estimating the blur kernel [80, 149].Typically, this is an ill-posed problem to solve. In order to make the non-convexoptimization converge to a good local optimum, we shall discriminatively designgood regularizers [167]. A straight-forward way is to directly apply the image pri-ors from non-blind deconvolution research (e.g. TV). In this work, we exploit thesimilar prior-based deconvolution in both non-blind (Chapter 3, 4, 5) and blindmanner (Chapter 6).2.4.2 Color fidelity preservationColor fidelity has driven many advances in both the optics and graphics com-munity [127, 137]. For conventional image acquisition, it is not a notable issuesince most refractive lenses are carefully designed with spectral invariant focusingpower [145, 169]. Thus, state-of-the-art deconvolution formations can be directlyapplied to computationally recover color images. However, the spectral responseof a DOE affects color fidelity drastically because its highly spectral variant fo-cusing power goes against the blur kernel convolution model in sRGB color space.Consequently, the ill-posed inverse problem may lead to unacceptable color arti-facts, even after the perceptual blur has been mostly removed. Preserving colorfidelity is crucial for any consumer level imaging device [148], which is also thegoal we seek in the diffractive achromat (Chapter 4) design and the encoded lens(Chapter 6) design.Regarding computational imaging with simplified optics, it is necessary to re-move color fringes on sharp edges resulted from different PSFs in channels. Fromthis point of view, many low-level techniques have been applied in complex optical20systems [72]. Schuler et al. [131] proposed an aberration correction algorithm for asingle lens in YUV color space. Chung et al. [22] measured chromatic aberrationsat edges through color differences and compensates them locally. Further, a convexcross-channel prior is developed and efficiently solved [57]. This method modelsthe three color channels separately and exploits the correlation between gradientsin different channels, which provides better localization. Alternatively, the sym-metry of the convolution kernel [132] and the geometric and visual priors [176]are investigated. Very recently, Sun et al. [150] investigated a blind deconvolutionscheme that included cross-channel sharing in the fitting model.State-of-the-art models yield reasonably good results with chromatic aberrationcorrected. Although correcting chromatic aberrations for diffractive imaging is alsoplausible using existing techniques, the residual color artifacts are inevitable dueto the inherent non-uniform distribution of spectral PSFs (Chapter 3), which leadsto a metamerism problem, as discussed in Sec. 2.2.2.5 Matrix factorizationAs a common mathematical discipline, matrix factorization describes a decompo-sition problem in computational design of multiplexing, for instance, our encodedlenses (Chapter 6) and mix-and-match holography (Chapter 7). The forward prob-lem yields a target matrix by multiplying two vectors (matrices), thus the inverseproblem is to find two vectors (matrices) that best approximate the target matrixin a least-squares sense [62, 66]. In graphics, a large number of works have usednon-negative matrix factorization (NMF) to enable a variety of layered displays.We specially analyze this topic in Sec. 2.6. There are many algorithms that havebeen developed aiming for solving NMF problems efficiently, that can be foundin this survey [84]. We here address two kinds of classical update rules includingthe common multiplicative algorithms [62] and projected gradient descent updatealgorithms [23, 91]. Refer to corresponding chapters for the detailed optimization.2.6 Computational image representationAnother active research direction in graphics is to represent the image for visual-ization via some designed optics. From the perspective of optics, both ray optics21theory and wave optics theory can be applied in the design space.2.6.1 Refractive caustic designsThroughout the last few years, a large body of computational caustic designs haveemerged in the graphics community, that mostly focusing on a refractive setting [67].We can distinguish between discrete designs [113, 114, 178] where the surfaceis broken into smaller elements, and continuous designs (see [36, 77] as wellas [134, 179]) that yield a single smooth surface. Discrete designs have the ad-vantage that they allow for discontinuities in the surface, which results in flattergeometries, and more lightweight, compact objects. On the other hand, the dis-continuities also introduce artifacts (including diffraction) that limit the contrastcompared to smooth designs. The design process of all the above work is based onray optics, and thus neglects diffraction.The concept of multiplexing images from two or more portions has been re-ported recently. Fuchs et al. [39] proposed a reflectance field display that com-bines different optical elements and encoded transparencies to allow for projectingdifferent light fields. Refractive optics in combination with multiple attenuationmasks have also been used for generating animations [6]. Other solutions includecombining structured white light with masked prisms [65] or combining singlesource image with a refractive lenslet array [114]. This last work in particular canreveal different hidden images when the lens is placed at prescribed orientations onthe source image or viewed from different angles. These works on steganography,display, and spectral color production with refractive means serve as an inspirationfor our own work with a holographic image formation.Instead of a purely refractive approach, Damberg et al. [27, 28] proposed adesign that uses a spatial phase modulator but also a geometric optics image for-mation model to arrive at a piecewise smooth phase profile with high diffractionefficiency. This work can be interpreted as a hybrid between refractive and diffrac-tive designs. Unfortunately, it remains very challenging to maintain the piecewisesmoothness while multiplexing multiple images into a single phase pattern.222.6.2 Diffractive displaysResearchers in graphics have proposed several designs for exploiting diffractive op-tics in caustic display applications. Recent attempts to fabricate the optical surfaceswith desirable bidirectional reflectance distribution function (BRDF) or BTDF havebeen investigated [49, 90, 170]. Ye et al. [170] proposed to use multiple diffractivelayers to encode a 4D BRDF or BTDF. Finally, Xue et al. [168] multiplexed mul-tiple images into a single diffraction grating at different wavelengths of light basedon an analytic model. These approaches may fail in our framework for two rea-sons: first, they are not intended for broadband imaging DOE designs; second, theyrely on either matured random search algorithms or their extended versions, low-ering the computational efficiency and optimization robustness in our diffractiveachromat or encoded lenses.2.6.3 Computer-generated holographyIn optics community, there has been a considerable amount of work on computergenerated hologram (CGH) for decades. Earlier work includes generating binaryholograms by computer drawing [92] and forming an image by a wavefront recon-struction device called Kinoform [86]. Later development of holography includesthe work on computationally designing complex transmittance or reflectance func-tions for holographic systems (see e.g. [15, 106, 122, 154]). There have been var-ious heuristic algorithms capable of generating desired phase or amplitude distri-butions considering variables in both the spatial and the spectral domain. Iterativemethods based on Gerchberg-Saxton (GS) search, simulated annealing (SA) algo-rithm, and direct binary search (DBS) have been applied to the optimization of bothmonochromatic and broadband DOEs [35, 44, 75, 123, 182].2.6.4 Compressive computational displaysCompressive computational displays are of relevance to our work since they par-tially compensate for the limited data bandwidth of current display hardware (e.g.LCD panels). The fundamental concept behind compressive displays is to encodehigh-dimensional data such as light fields into a small number of either multiplica-tive (e.g. [58, 82, 164]) or additive planes (e.g. [85]). Similar methods have been23applied to achieve super-resolution effects (e.g. [10, 129]). The methodology be-hind these methods is often based on non-negative matrix or tensor factorization,and related optimization problems. This is also the case for our framework, whichis most similar to work by Heide et al. [59], but operates in a diffractive setting.2.6.5 Multiplexed hologramsMulti-layer or volume holograms can provide further increased bandwidth andthus allow for additional degrees-of-freedom in the design process [45]. Usingangular multiplexing, Mok [104] stored many holograms within a single crystalof LiNbO3 to increase the bandwidth. A survey of this inspiring body of workcan be found [121], which also investigates their application potential in data en-cryption and storage. Our mix-and-match design differs from conventional volumeholograms that use reference-recording scheme with either wavelength or angularmultiplexing to encode the data. Alternatively, Brand [14] proposed the specu-lar holography that uses differential geometry to control light steering for severalviewing geometries, which is not derived from diffractive optics.Finally, Shi et al. [138] and Meng et al. [99] investigated similar concepts ofhiding multiple images using two or more layered phase holograms. However,their setting is quite different from ours (Chapter 7) in several ways: first, theirwork results in a static design of a fixed pair of front and rear surface, while weinstead use a combination of phase retrieval with complex matrix factorization toenable mix-and-match designs of many pairings. Second, in their scenario, multi-ple images are represented as essentially a depth-dependent point spread function,i.e. the decoding step involves moving the projection screen. By comparison, ourdesign decodes different holograms and images by changing the relative geometricalignment of the DOEs themselves.We envision that our mix-and-match holography scheme (Chapter 7) can beapplied in a number of scenarios. For instance, one can leverage diffractive prop-agation to obtain high quality imaging or encryption visualization with (virtual)alignment for (digital) security applications. Besides, one can obtain actual dis-play functionality with state-of-the-art hardware capabilities, like using fabricateddiffractive optical element (DOE) or commercial spatial light modulator (SLM).242.7 Fabrication techniques of flat opticsUnlike conventional computer numeric control (CNC) tooling methods [27] thatare generally used in fabricating refractive optics with continuous surface profiles,micro-fabrication techniques usually involve several repeated procedures that areoperated in an alignment accuracy of micro-meter or even sub-micro-meter level.Precisely, either continuous or discrete surface profiles are optional in fabricatingDOEs [31, 183]. One typical realization is to use photolithography technique,which uses light to transfer a geometric pattern from a photomask to a light-sensitive chemical photoresist on a substrate using UV light exposure [141].State-of-the-art photolithography techniques can create the relatively small pat-terns with the feature size going down to 1 µm. The main disadvantages of pho-tolithography currently are the relative complex procedure and the limited patternsize it can create at one time process. More advanced fabrication approaches in-clude grayscale photolithography and nano-imprinting [1]. These two techniquesare able to create micro-structures with more smooth profile, which further lead toimproving diffraction efficiency of some sense. In this thesis work, we repeatedlyapply four iterations the lithography technique followed by a reactive ion etch (RIE)step, that a certain amount of material in the exposed areas on the substrate is re-moved by chemically reactive plasma. The specifications of our on-hand fabri-cation and a diagram of state-of-the-art iterative photolithgraphy fabrication areprovided in Appendix chapter A.25Chapter 3Lightweight Diffractive-refractiveHybrid LensesAs highlighted in Chapter 2, diffractive optical elements (DOEs) show great promisefor imaging optics that are thinner and more lightweight than conventional refrac-tive lenses while preserving their light efficiency. Unfortunately, severe spectraldispersion currently limits their use in consumer-level lens design.In this chapter we jointly design lightweight diffractive-refractive optics andpost-processing algorithms to enable imaging under white light illumination. Us-ing the Fresnel lens as a general platform, we show three phase-plate designs, in-cluding a super-thin stacked plate design, a diffractive-refractive-hybrid lens, and aphase coded-aperture lens. Combined with cross-channel deconvolution algorithm,both spherical and chromatic aberrations are corrected. Experimental results indi-cate that using our computational imaging approach, diffractive-refractive optics isan alternative candidate to build light efficient and thin optics. The primary conceptaddressed in this chapter and the corresponding experimentation together fit as thefirst proof-of-concept solution for the following chapters.3.1 IntroductionModern photography needs an efficient light acquisition procedure, in which highquality lens design is key. In the past decades, the revolution from film-based to26digital photography and steadily increasing processing power in the capturing de-vices have enabled many computational imaging applications. This computationalimaging technology has created a demand for new and flexible lens systems [55].However, designing an excellent digital imaging device requires consideration ofmulti-dimensional factors, such as light acquisition efficiency, material cost, size,design flexibility, and computational complexity. In graphics community, a com-mon approach to handling the imperfections of optical systems is to remove theremaining optical aberrations in a post-capture deconvolution step.Computational imaging with simplified lenses [57] provide the possibility tomotivate inexpensive computer vision devices. From this perspective, flat lenseslike DOEs are promising with three advantages over conventional refractive lenses:first, they can be fabricated on a thin sheet; second, a single DOE can perform manyoptical operations simultaneously, making it promising to act as a light modulationplatform; third, because of the relatively flat structure, diffractive optics reducesdistortions and therefore has much more uniform PSF distributions across the im-age plane. The last point is particularly beneficial for a post-capture deconvolutionstep.Instead of purely relying on refractive optics, as the very first attempt through-out this thesis work, we seek an alternative way, diffractive phase modulation [52],to benefit modern lens design. One the one hand, we could obtain thinner struc-tures, less material cost, better off-axis imaging performance than conventionalrefractive lenses. On the other hand, we could seek better chromatic performancethan a single DOE. We focus on diffractive-refractive optics design jointly withcomputation to form a complete pipeline as shown in Fig. 3.1. Our method al-lows capturing a blurry intermediate image by a super thin phase modulated DOE,and then recovering the sharp image with deconvolution algorithms using cross-channel regularizations. Consequently, compact and lightweight broadband cam-eras could be designed. The results of our method indicate that our algorithmcorrects the residual aberrations and preserves the details. We believe it is of con-siderable significance for follow-up research designs.The main technical contributions include: We propose a white light computa-tional photography method using diffractive-refractive optics to design lenses withthin structures at a controllable cost. We analyze diffractive phase modulation to27Figure 3.1: Compared to the imaging pipeline of conventional imaging withwell corrected complex refractive lenses (top), our diffractive-refractivecomputational imaging allows capturing a blurry intermediate image bya super thin phase modulated DOE, and then recovering the sharp im-age with deconvolution algorithms using cross-channel regularizations.An exemplary blur image (bottom left) and its corresponding recoveredimage (bottom right) are shown. Insets of the images indicate that ouralgorithm corrects the residual aberrations and preserves the details.offer considerable benefits for broadband imaging by allowing flexible compro-mise between diverse efficiencies in different application scenarios, including theuse of a super-thin diffractive optical lens, a diffractive-refractive-hybrid lens, anda phase coded-aperture lens. We employ a fast but robust deconvolution to correctchromatic aberrations introduced by diffraction.3.2 Diffractive-refractive optics imagingThe goal is to provide alternative options for originally complex lens design inmodern photography by implementing phase modulation on a thin substrate.283.2.1 Optical characteristics of Fresnel lensesWe address two optical characteristics of Fresnel zone-plate (FZP), multi-orderdiffraction and wavelength sensitivity, before introducing the imaging model.Multi-order diffraction. The multi-order diffraction is accounted for as the mainreason for low light efficiency of amplitude-type DOEs. Due to the fact that theincident light is distributed into a series of diffraction orders as shown in Fig. 3.2,diffraction efficiency is usually very low. Ideally, if a continuous profile could befabricated in a phase-only DOE, the diffraction efficiency is 100%. In practice,continuous profiles are difficult to fabricate. Instead, they are approximated bymulti-level binary optical elements. We use multi-level phase zone plates (PZPs)in our implementation for high light efficiency. For a N-level DOE, the diffractionefficiency [151] is given in Eq. 3.1:ηmN =sin[pi((n−1)dλ −m)]pi((n−1)dλ −m)2 sin[pi((n−1)dλN)]pi((n−1)dλN)2, (3.1)where λ is the designed principal wavelength,m is the diffraction order, d is totalheight of the micro-structures on the substrate, and n is the refractive index at thedesign wavelength. With as many levels as one can fabricate, a discrete PZP canapproximate a continuous Fresnel lens with a high enough diffraction efficiency. Arealistic design for a DOE has 16 layers, in which case the diffraction efficiency is99%. This number drops to 95% with 8 layers.Wavelength-dependent imaging. PZPs are wavelength-dependent DOEs. Whenapplied in broadband imaging, especially under visible white light illumination,PZPs could be focused for the green wavelength on the image plane, but red wave-length would focus in the front while blue wavelength would focus in back of theimage plane. The difference of the focal lengths for different wavelengths is themain drawback any DOE suffers from. The equivalent Abbe number for a PZPis -3.45 [112]. It differs significantly from refractive lenses (usually a large posi-tive number) both in the chromatic distribution and magnitude, leading to a strong29Figure 3.2: A binary amplitude Fresnel zone plate (top-left) diffracts lightinto multi-orders (0, ±1, ±3, ...), resulting in a low diffraction effi-ciency. With continuous phase-only diffractive lens (top-right), the lightfalls into +1 order could theoretically reach 100%. In practice, multi-level microstructures (bottom) are usually easy to fabricate and approx-imate the continuous profile very well.negative chromatic aberration.3.2.2 Imaging pipelineAs indicated in Fig. 3.1, the imaging pipeline comprises an optical design partand a following optimization procedure. The optical design aims at pre-defininga point spread function (PSF) with expected energy distribution in the spatial andspectral domain. The PSF could be obtained by scalar diffraction theory [52] as theFourier transform of the complex transmittance function of a lens. The resultingPZPs are fabricated by state-of-the-art photolithography techniques. Our optimiza-tion is designed jointly with the optical design. It corrects the residual chromaticaberrations that are introduced by the use of DOEs. We will demonstrate that ourmethod of combining the optical design and algorithm design shows comparableimage quality.3.3 Optical designsWe carry out a study of three designs in the PZP form to make reasonable com-promises among the diverse dimensions of the design space mentioned previously— super-thin diffractive optics, diffractive-refractive-hybrid lens and phase-coded-aperture lens. We fix two focal lengths as 50 mm and 25 mm for each design with30Figure 3.3: Schematic comparisons of a super-thin diffractive optics (top), adiffractive-refractive-hybrid lens (middle), and a simple refractive lens(bottom) for f1 = 50 mm (left) and f2 = 25 mm (right). The R,G,B colorsrepresent the focal powers for different wavelengths. “+” indicates thefocal length for red is longer than that of blue, and “−” for the otherway round.the same aperture diameter 12.7 mm for all designs. The design parameters andcomparisons are illustrated in Fig. 3.3.Super-thin diffractive optics. First, we use a single PZP or stacked multiple PZPsas a lens, of which the maximum thickness could be only 1 mm. We make thefollowing design comparison to illustrate the benefits of the super-thin DOE (seethe top row in Fig. 3.3). By maintaining the same optical power, the phase-platedesign is much thinner and lightweight than the refractive lens. The reductionof thickness and volume will become crucial when designing a lens with a largeraperture but shorter focal length.Although a multi-level PZP gains high light efficiency and high material ef-ficiency, its aberration performance is generally not satisfactory for white light31imaging, especially with respect to (negative) chromatic aberration, but also someresidual spherical aberration. These aberrations increase the computational com-plexity of the reconstruction step, and limit its efficacy.Diffractive-refractive-hybrid lens. Instead of shifting all chromatic aberrations tothe computation side, we design a diffractive-refractive-hybrid lens. Specifically,a simple lens is attached to the above super-thin DOE design, therefore, the hugenegative chromatic gap is partially compensated for by the positive chromatic aber-ration of the refractive lens (see the middle row in Fig. 3.3).Note that the idea is not to design an achromatic doublet, such as the onesin [100, 107]. Instead, we seek a good compromise between obtaining a betterPSF distribution for post-processing and maintaining a thin structure. The attachedsimple lens can be a cheap off-the-shelf lens. Our designed PZP still takes charge ofmajor focusing ability. For maintaining the same optical power as well as keepinga compact structure, one can achieve a compromised optical performance betweena pure phase-plate and a simple lens. Besides, the hybrid design helps weaken thenegative influence of other aberrations, compared with using a simple lens only.Phase-coded-aperture lens. Apart from its super-thin structure benefit, a DOEhas the advantage of design freedom. One can design expected functions on asingle DOE, rather than applying a physical coded-aperture or multiple-capturetechniques on a refractive lens. To modulate the energy distribution of a PSF,further to discrete information of light rays, one can gain additional benefits forcomputational processing [87]. Generally, a designed binary mask is configuredat the aperture plane to selectively block light rays. On the contrary, we yield thecompromise of blocking light rays to encoding the phase distribution of a PZP.Therefore, through interference and diffraction of the sub-structure on the PZP, theexpected PSF shaping is obtained. During this procedure, the total light energy isalmost preserved.According to the above analysis, we designed two 16-level DOEs as shown inFig. 3.4. Referring to Eq. (3.1), a 16-level profile can already provide a satisfactorylight efficiency in common lighting environment. All DOEs are designed at the32Figure 3.4: Left: Microscopic images of our PZP and phase-coded-aperturelens on a Nikon Eclipse L200N 5X microscope. Right: Real color PSFs(normalized) captured are shown with their individual RGB compo-nents. Note that the axises in the plots denote the pixel domain, andthe colorbars indicate the relative intensity.principal wavelength 550 nm, thereby, the red and blue channels suffer from se-vere blurring with considerable noise. Note that here in the figure, the gray scalesindicate the discrete phase distributions in the central area. By making use of theanti-symmetric phase distribution (see Fig. 3.4 on the right), we divide the 2D cir-cular plane by 6 sectors, each of which has a symmetric counterpart with the phaseshifted by pi . The original Gaussian-like PSF can then be encoded into a three-pikestar PSF.3.4 Image reconstructionIn this section, we analyze the post-capture computation to recover the imagescaptured by the presented diffractive-refractive optical designs.3.4.1 Imaging modelWith respect to diffractive phase modulation, the image model (recall Sec. 2, Chap-ter 2) is subject to the variation of wavelength information. The blur kernels Pc forthe RGB color channels c ∈ {1,2,3} are explained as integrals of the PSFs in therespective spectral range λc:Pc =∫λcP(λ )dλ , (3.2)33λ here is the wavelength of input light ray. Thus, in general the PSF is contentdependent. A pre-calibrated PSF can therefore only reduce aberrations to a certainlevel, but never completely.3.4.2 Optimization methodWe next describe the image reconstruction as an optimization solving the inverseproblem. Our approach follows closely the recent work by Heide et al. [60], withseveral improvements to exploit the specific structure of chromatic aberrations inDOE optics, as well as several improvements to computational efficiency.We achieve state-of-the-art deblurring results by using Bayesian optimizationthat exploits patch-based correlation and cross-channel gradient correlation in theunknown latent image, that can be cast as an optimization problem, as follows:ic = argminicµc2‖bc−Kcic‖22+Γ(ic), (3.3)where the first term is a standard least-square data fitting term, µc is the weightof this term for channel c, and the second term Γ(ic) is a collection of priors thatregularize the reconstruction. The regularization terms are taken from [60] and arebriefly summarized in the following.Cross-channel prior. Diffractive optics tend to produce good image quality for asingle design wavelength. In white light applications, chromatic aberration limitsthe overall image quality, but one channel in the RGB image can usually be focusedquite well and is significantly sharper than the other channels. Heide et al. [57]proposed a cross-channel prior that we can modify to take advantage of this fact.The idea is that, in the real world the edges in different color channels appear atthe same place. Furthermore, normalizing with the channel intensity allows forluma gradients, but discourages chroma gradients. For a pair of channels l ,m thecross-channel prior is:∇im · il ≈ ∇il · im⇔ ∇im/im ≈ ∇il/il. (3.4)34Therefore, the regularizer term is rewritten as:Γ(ic) = ∑m 6=l ‖ D · im · il−D · ii · im ‖1 , (3.5)where il,m stand for the images convolved with different blur kernels, again withthe l1 norm. The convolution matrix D{x,y} implements the first derivatives in x andy direction of the image.Patch-based prior. A DOE suffers from more aberrated blur kernels than the onesmeasured for a simple refractive lens. Our reconstructions may be sensitive to noiseand extreme chromatic aberrations. We use a slight modification of the patch-basedBM3D prior [26, 60]. To exploit even more cross-channel information, we per-form multiple BM3D passes with block-matching and 3D thresholding in differentcolor spaces, e.g. the OPP and YUV color space in our work. Using differentcolor spaces in this process has the effect that the distance of different colors areweighted differently according to the color space itself. Therefore using differentcolor spaces affects the similarity based matching and the 3D thresholding of theBM3D method, which gives us significant improvements.3.4.3 Analysis on efficient optimizationInspecting again Eq. (3.3), we have only included the blurry convolution matrixKc in the data term, but no demosaicking mask as proposed in [60]. This hasthe benefit that we can efficiently solve Eq. (3.3) in the Fourier domain. While aspatial domain solver (e.g. conjugate gradient (CG) methods) would add significantcomputational overhead, in general, results can be dramatically improved whenjointly solving demosaicking and deblurring. However, in our specific case, thePSFs may distribute more than a hundred pixels. Sequential demosaicking doesnot change the results much, but can drastically increase the speed. We also foundthat only a small number of iterations of the outer loop are sufficient to convergeto a good image quality. Despite these performance improvements, we note thatour method is effective at recovering good quality images from very challenginginput data: the PSFs in the two blurry channels are of very large size and classicalmethods do fail here. For the data with the largest PSFs – corresponding to a design35with a small f -number – we employ a multi-scale deconvolution scheme [174] tofurther speed up the processing.3.5 ResultsWe show a number of experimental results implementing the above designs. Forthe figures shown, some are cropped regions to highlight the details.3.5.1 PrototypesIn the fabrication of DOEs, either continuous or discrete surface profiles are op-tional if different manufacture methods are adopted [31, 183]. We build the fol-lowing prototypes by repeatedly performing the photolithography and reactive ionetching techniques. The feature size for the fabrication is 1 µm, which is the min-imum size on the edges in all our designs. To simplify the fabrication, we dis-cretized the profile into the 16 levels. The diameter of all the lenses are 8.0 mm.Accordingly, we use off-the-shelf Thorlabs lens tubes and plano-convex lenses forbuilding all the prototypes. Further we consider to use commercial scientific cam-eras (PointGrey GS3-U3-50S5C) to implement all the experiments. The pixel pitchof the sensor is 3.45 µm. It has a full well capacity of 6,402e- and outputs 14-bitimages with maximum SNR 38 dB. In the experiments, we measured the noiselevel under the same illumination setting and exposure setting as in Fig. 3.5. Thesensor characterized noise level is represented by the variance of 8.233×10−5 withno gain added in the sensor setting, and the exposure time set to 100 ms. All imagesin the experiments are based on one single exposure image.Super-thin DOE optics. We weight more on the super-thin and lightweight expec-tation such that we tolerate a certain amount of aberrations in hardware. One single16-level PZP ( f = 100 mm, D = 8 mm) can be exclusively used as an imaging lenswith a thickness of only 0.5 mm. Moreover, two PZPs can be stacked together tomake a doublet lens f = 50 mm with the thickness only 1mm. Figure 3.5 presentsthe results of the single and stacked DOEs (b) and (c). From which we see thatthere is some light efficiency loss due to the fabrication errors in our prototype,reflected by a slight drop of image contrast. However, the deconvolved results still36Figure 3.5: Top of (b-f): The same scene captured by our single PZP ( f =100 mm), stacked PZP ( f = 50 mm), hybrid coded-aperture PZP andrefractive lens ( f = 50 mm), hybrid PZP and refractive lens ( f = 50mm), and a simple lens ( f = 50 mm) respectively. Middle of (b-f):Corresponding deconvolved images for each lens. Bottom of (b-f): R,G and B components of the PSFs for each setup. The scene is projectedonto a board by a projector. Note that we only show the same croppedarea for cross comparison here. The ground truth image captured bya Canon DSLR is presented for reference (a-2), with the PSNR plotprovided (a-1) corresponding to the designs (c-f).make sense in those scenarios where high-quality imaging is not strictly required,but the lens volume and thickness matter significantly.Diffractive-refractive-hybrid lens. For the prototype of the hybrid design, one 16-level PZP ( f = 100 mm, D = 8 mm) is attached to the back side of a plano-convexlens ( f = 100 mm) to make a group lens f = 50 mm with the optics thicknessless than 2.6 mm. Figure 3.5 present the results of the diffractive-refractive-hybridlens, and conventional refractive lens (e) and (f). From this result we see thatthe recovered image is comparable with the one recovered from a simple lens,and even slightly better in off-axis regions. Compared to pure DOE designs, theessential huge negative chromatic gap has been mitigated to a level that the cross-channel deconvolution can handle. Compared to single refractive lens, the off-axisspherical aberration behavior is better with the thickness slightly reduced.37Phase-coded-aperture lens. We fabricate a 6-sector phase plate with anti-symmetricstructures. Same as the above prototypes, we attach a refractive lens for experimen-tal convenience, such that we have the effective f = 50 mm with aperture size D =8 mm. Figure 3.5 presents the result of the phase-coded-aperture lens (d). We en-vision more latent information is preserved due to the shaping of PSFs, especiallyin red and blue channels, which shall benefit the deconvolution procedure. Notethat here we only implement a possible coded-aperture design.In Fig. 3.5, we compare the performance of the three designs for the same housescene, which was projected onto a board. We show a ground truth image capturedby a Canon DSLR in the bottom left for reference. The RGB components of thecorresponding PSFs are shown in the bottom row to demonstrate the individualcharacteristics of our designs. We seek to use PSNR to intuitively present thebenefit and the downside of our designs. Despite some color reproduction artifacts,our designs generate detailed, sharp images that are visually comparable to therefractive simple lens results. The hybrid designs (c) and (d) indeed help to improvethe image quality. These comparisons support our intentions of trading-off theimage quality and lightweight requirement in the design space.3.5.2 Experimental resultsWe then present several pairs of selected deconvolved results of natural scenes tovalidate the cross-channel optimization applicability, see Fig. 3.6, Fig. 3.7, Fig. 3.8,and Fig. 3.9 for the details. In each result, we show the captured blurry image onthe left, and the deconvolved sharp image on the right. Insets of the cropped imagesshow our algorithm recovers the spatial details well.We use non-blind deconvolution, where the spatial variant PSFs are calibratedthrough deconvolving a series of random patterns or using optical simulations. ThePSF estimation procedure is described in detail in [57]. All experiments are per-formed on a core of Intel Core2 i7 CPU with 2.40 GHz and 8GB RAM. The maincomputational cost is the cross-channel based optimization to deblur the image andremove chromatic aberrations. Generally, an image with around 2 megapixel sizeand 100 pixel size of moderate blur kernel needs less than 1 minute to get accept-able recovery with cross-channel prior, and 3 more minutes with one BM3D pass.38Figure 3.6: Top: Images captured and deconvolved for a text scene with asingle PZP ( f = 100 mm). The teapot and white texts portions show theability of our algorithms to remove most of the chromatic aberrations.Bottom: Captured and deconvolved images for a flower and book scenewith stacked PZPs ( f = 50 mm). The yellow and green objects are wellrecovered in the same scene. Note that both scenes are at the samedistance from the lenses, and the field of views differ because of thefocal lengths.As mentioned previously, we implement two priors jointly in the deconvolutionprocedure.Spatial behavior of the PSF. A DOE exhibits lower geometric distortion than asimple lens, as illustrated in Fig. 3.10. Therefore, the monochromatic blur is closerto the spatial invariant assumption, where PSF estimation for large field of view isnot practical. This is a main advantage of introducing diffractive-refractive opticsinto the deconvolution problem. Note that the relatively lower contrast of the lattertwo is partially accounted on the light efficiency issue of our fabricated prototypes.39Figure 3.7: Images captured (a) and deconvolved (b) for a box scene withour diffractive-refractive-hybrid lens ( f = 50 mm). By introducing arefractive lens, chromatic aberrations are smaller than in the single PZPand stacked PZP cases in Fig. 3.6. The inverse problem becomes morewell-conditioned and shaper images could be recovered.3.5.3 DiscussionAcquisition efficiency. We jointly consider three aspects, hardware structure (re-flected by optical thickness and material volume cost), light transport value (re-flected by the NA number), pre-conditioning performance for post-processing (eval-uated by PSF distribution), regarding the efficiency requirements for consumer-level capture devices.The introduced hybrid design helps to maintain a relatively compromised ef-fective f -number, in our case around f -6. With respect to hardware structure, allDOEs we fabricated are with the thickness 0.5 mm. The fabrication technique maybehave high-cost at the very beginning for a research lab, however, once a maskhas been fabricated, mass product manufacturing can operate smoothly. The sub-strate wafer is made from normal glass, thereby diffractive-refractive optics is morematerial efficient compared with modern refractive lenses, which usually require40Figure 3.8: Images captured (a) and deconvolved (b) for a ball scene with ourphase-coded-aperture hybrid lens ( f = 50 mm). By coding the intensitydistribution of the PSF, high frequency features as well as colors can bepreserved better. Note that the net around the ball is recovered whichotherwise would hardly be seen in the blurry image.expensive glasses.Image reconstruction robustness. The cross-channel prior (recall Sec. 3.4.2) addedis itself not guaranteed convergence to a global optimum for all channels. How-ever, in most cases, choosing an appropriate solver and appropriate parameters canstill give reasonable good results.Our design exhibits a robust reconstruction with a certain depth variation, il-lustrated in Fig. 3.11. For the outdoor scene with around 2 to 8 m depth variation,we see that most chromatic aberration has been removed, with the exception ofsmall artifacts in the red component of the background depth. We regard this asa limitation in our experiment due to insufficient PSF estimation. Further, simi-lar to the depth of field (DOF)-extension potential in [24], DOEs focus differentwavelengths at different depths. Therefore, by carefully selecting the weights ofthe cross-channel term for each spatial patch, one can recover an image with ex-41Figure 3.9: Images captured (top) and deconvolved (bottom) for theISO12233 resolution chart with our three designs: (from left to right)hybrid PZP and lens, hybrid phase-coded-aperture PZP and lens, andstacked PZPs. The focal lengths are f = 50 mm for each design.tended DOF. Notably, the deconvolved results show a haze effect that reduces theimage contrast. We analyse this artifact together with other DOE designs in latertext (Sec. 4 in Chapter 4).Limitations and future work. Several essential limitations can not be ignoredwhen applying diffractive phase modulation. Firstly, the number of levels of aDOE in practical manufacturing is dependent on its aperture size and required fo-cal power. Generally, the phase modulation designs with discrete profile behavelower in diffraction efficiency, especially when being operated in harmonic diffrac-tion orders [152], than conventional refractive lenses do. Fortunately, the ease ofbuilding larger apertures without increasing the thickness can compensate for theefficiency loss. The photolithography process used for our prototypes limits thefeature size in the DOEs to 1 µm. This directly affects the quality of the DOE, aswell as the ability to scale down the optical systems to smaller cameras and sensors.Moving to e-beam lithography process would alleviate this limitation.The chromatic performance is not compatible as that of a carefully designedlens currently. Although the cross-channel prior has removed most chromatic aber-rations resulted from the phase modulation, color recovery is not entirely matched42Figure 3.10: Off-axis performance comparison for a simple refractive lens(left), our hybrid PZP lens (center), and stacked PZPs (right). Notethat the first row presents the capture images of a black-white checker-board, while the close-ups illustrate the on-axis and off-axis patches intheir green channel, red channel, and blue channel visualization. Fromthe comparison of patch A and patch B, we observe that our designsshow better spatial uniformity, despite the chromatic aberration. Allthe focal lengths are f = 50 mm and the full field of view is around30◦.to the ground-truth data due to a possible mismatch between color spectrum ofcaptured data and estimated PSFs. From the above deconvolution results, we stillobserve a fair amount of artifacts, including color ringing effects of edges and non-uniformity of less high-frequency containing regions. A color calibration step is tobe added in the post-processing, which is essentially scene-dependent though.The current design is not intended for replacing high-end DSLR lenses, butproviding alternative options for camera designs with special form factor require-ments. Considering the imaging performance of Fresnel lens designs, we regardone promising application is to act as the capture lens of modern mobile device orportable virtual reality device. In this market, the thickness and volume cost weightover the other factors. Referring to the current high-end mobile devices, even 0.5mm thinner can be a breakthrough.43Figure 3.11: Images of captured (left) and deconvolved (right) for an outdoorscene with our hybrid PZP lens (f = 50 mm). The far and near objectsvary from around 8 m (house) to around 2 m (tree) from the lens.3.6 ConclusionWe propose a novel computational camera that bridges diffractive-refractive op-tics and computational algorithms, to build thin, compact and lightweight imagingsystems. Three prototype designs, using a combination of custom nano-fabricatedDOEs and off-the-shelf products, are implemented to validate our concept. The hy-brid designs hold smaller chromatic aberration than pure DOEs, while thinner andlighter, and better off-axis performance than that of a refractive lens. We employ afast image deconvolution with cross-channel BM3D prior to tackle the chromaticaberration and the noise jointly.In combining diffractive-refractive optics and computation algorithms, mod-ern imaging devices can be designed towards larger aperture, thinner structure,more flexible designs while still maintaining visually pleasing imaging perfor-mance. Overall speaking, our approach provides an alternative technical solutionto design and fabricate fashion thin lenses. This chapter acts as the fundamentalbasis in terms of validating the practicability of combining diffractive optics andpost-processing to build compact or even flat lenses. The analysis of diffractionefficiency and fabrication procedure hold for the following chapters.44Chapter 4The Diffractive AchromatWe can observe an insight from the previous chapter that the inherent strong dis-persion of DOEs is a tremendous obstacle that limits their use in full spectrumimaging, causing unacceptable loss of color fidelity in the images. In particular,metamerism introduces a data dependency in the image blur, which has been ne-glected in computational imaging methods so far.In this chapter we introduce both a diffractive achromat based on computa-tional optimization, as well as a corresponding algorithm for correction of residualaberrations. Using this approach, we demonstrate high fidelity color diffractive-only imaging over the full visible spectrum. In the optical design, the height profileof a diffractive lens is optimized to balance the focusing contributions of differentwavelengths for a specific focal length. The spectral point spread functions (PSFs)become nearly identical to each other, creating approximately spectrally invariantblur kernels. This property guarantees good color preservation in the captured im-age and facilitates the correction of residual aberrations in our fast two-step decon-volution without additional color priors. We demonstrate our design of diffractiveachromat on a 0.5mm ultrathin substrate by photolithography techniques. Experi-mental results show that our achromatic diffractive lens produces high color fidelityand better image quality in the full visible spectrum.454.1 IntroductionHigh quality imaging with reduced optical complexity has for a long time been thetarget of investigation in both academic and industrial research. In conventionalimaging systems, ever increasing optical complexity is inevitable because higherand higher sensor resolutions require ever improved correction of aberrations ofall kinds. Recent advances in computational imaging have introduced computationas a virtual component that can shift the burden from optics to algorithms. Thisallows for significantly reduced optical complexity while maintaining high imagefidelity at full sensor resolution and realistic apertures (e.g. [57, 131]). In partic-ular, DOEs have drawn great attention because of their ultrathin and lightweightphysical structure, a large, flexible design space, availability of mature fabrica-tion techniques, as well as better off-axis imaging behavior. Integrating diffractiveimaging elements and computational methods in a single imaging system has re-sulted in several new computational imaging devices with ultra compactness in thepast few years (e.g. [47, 117, 147]).As aforementioned, two main problems still exist for diffractive imaging ap-plied in the full visible spectrum — first, the wavelength-dependency of diffractionleads to strong chromatic aberrations that degrade the image quality with blurs ofa very large diameter; second, the strong chromatic aberrations in addition causea significant wavelength dependency of the point spread functions (PSFs) evenwithin a single color channel. In particular, this wavelength dependency means thatobjects with the same RGB color are blurred differently if the underlying spectraldistributions differ. This metamerism problem is particularly pronounced for thestrong dispersion present in diffractive optics.We aim to overcome the above limitations of diffractive imaging in the full visi-ble spectrum by introducing not only improvements in the deconvolution, but moreimportantly, optimizing the diffractive optics itself. We find that the color fidelityloss in conventional diffractive imaging is caused by the inherent non-uniformityof spectral PSFs, which inspires us to design a diffractive achromat by optimizingthe surface profile to produce nearly identical spectral PSF distribution for a seriesof wavelengths. The benefit of this strategy is twofold. On the one hand, chromaticaberrations are reduced because of the balance among spectral PSFs. On the other46hand, the quasi-uniform spectral PSFs significantly improve the color reproductionof the captured image. Effectively, we sacrifice sharpness in a single channel forspectral uniformity of PSF, which we can exploit to facilitate robust and efficientdeconvolution. In addition to validating the usability of diffractive imaging understate-of-the-art deconvolution schemes, we explore a two-step cross-scale decon-volution scheme to recover sharp and hight color fidelity images.The technical contributions are as follows:• We introduce the diffractive achromat, a diffractive imaging system for full-spectrum visible light imaging that combines optimization in both diffractiveoptical design and post-capture image reconstruction.• We employ an effective optimization method for designing achromatic DOEssubject to diverse achromatic requirements, which rearranges the spatial andspectral distributions of PSFs so that chromatic aberrations and color corrup-tion are largely eliminated in hardware.• We propose a cross-scale prior in the deconvolution to further mitigate theaberrations introduced by the diffractive optics. Benefiting from the opti-mized uniform PSFs, our method is more robust and efficient than state ofthe art.• We build a prototype achromatic diffractive lens on a 0.5 mm ultrathin glassplate to validate the practicality of our ultrathin and lightweight diffractiveimaging in full visible spectrum for different indoor and outdoor scenarios.4.1.1 Diffractive imaging modelImage formation The recorded image in channel c is an integration of spectralimages over the wavelength range Λ, weighted by the spectral response Qc(λ ) ofthe sensors for that channel. Each spectral image reflects the joint modulation ofillumination, surface reflectance and sensor response. This process can be writtenin Eq. 4.1:bc(x,y) =∫ΛQc(λ ) ·A (i(x,y;λ ))dλ , (4.1)47where i(x,y;λ ) is the latent spectral image, and A (·) denotes an operator describ-ing the aberrations of the lens.In Fourier optics, the incoherent imaging process is modeled as a convolutionof the latent image and the system (intensity) PSF, so the aberration operator isdefined in Eq. 4.2:A (i(x,y;λ )) = i(x,y;λ )⊗|g(x,y;λ )|2, (4.2)where ⊗ denotes 2D convolution. g(x,y;λ ) is spectral amplitude PSF, from whichthe spectral intensity PSF can be derived as k(x,y;λ ) = |g(x,y;λ )|2. The amplitudePSF can be further derived from scalar diffraction theory [52] as:g(x,y;λ ) =Aλ zi∫∫P(u,v;λ )exp(−j 2piλ zi(ux+ vy))dudv, (4.3)where A is a constant amplitude, zi is the distance from the lens to the imageplane and (u,v) are coordinates on the lens plane. A generalized pupil functionP(u,v;λ ) accounts for the lens:P(u,v;λ ) = P(u,v)exp(jΦ(u,v)) , (4.4)where the aperture function P(u,v) is usually a circ function. The phase termΦ(u,v) describes the phase retardation of light for each point on the aperture, whichin a general imaging system could be caused by either refractive or diffractiveoptics, or a combination of the two. In our case, Φ(u,v) is the function that willbe optimized to achieve a desired lensing effect and PSF. Using Eq. 4.4, we canrewrite Eq. 4.1 as:bc(x,y) =∫ΛQc(λ )i(x,y;λ )⊗ k(x,y;λ )dλ=∫ ∞−∞∫ΛQc(λ )i(ξ ,η ;λ )k(x−ξ ,y−η ;λ )dλdξdη .(4.5)For a conventional diffractive lens, the spectral PSFs k(x,y;λ ) are highly wave-length dependent. Therefore, the PSF is not separable from the inner integrationover wavelength. This effect is usually neglected in state-of-the-art image forma-48tion models for deblurring, where the blur kernel is assumed to be convolved withthe latent color image. However, this approximation does not hold for large chro-matic aberrations, and therefore current deconvolution algorithms fail to recoverthe latent image with high color fidelity.To guarantee the deconvolution algorithms work in RGB space, we design thespectral PSFs to be nearly wavelength independent, at least over the spectral sup-port of each color channel. Only then can we use the approximation k(x,y;λ ) ≈kc(x,y), resulting in Eq. 4.6:bc(x,y)≈∫ ∞−∞kc(x−ξ ,y−η)∫σQc(λ )i(ξ ,η ;λ )dλ︸ ︷︷ ︸ic(x,y)dξdη= kc(x,y)⊗ ic(x,y),(4.6)where ic(x,y) is the latent color image in channel c.For an RGB image, the vector form of Eq. 4.6 can then be written as:bc = Kcic, c = 1,2,3. (4.7)As long as we can design a diffractive lens with nearly constant spectral PSFs (i.e.a diffractive achromat), the convolutional image formation model holds again.Imaging approach overview Our achromatic diffractive imaging approach con-sists of two main parts — achromatic lens design and post-capture processing.In the lens design, we devise a method that optimizes a DOE to focus a rangeof wavelengths at a certain target focal length while maintaining both a compactPSF as well as spectral uniformity. This optimization is performed for a dense setof discrete wavelengths that are uniformly distributed over the target range.Once we have a diffractive lens that exhibits spectral invariant kernel behavior,the post-capture processing is modeled as an optimization problem to solve theinverse problem of Eq. 4.7. A two-step, cross-scale deconvolution is presented toreconstruct the image.494.2 Optical design optimization4.2.1 Optimization modelWe have seen from Eqs. 4.3 and 4.4 that the spectral PSFs are determined by thephase profile of the diffractive lens. Further, the phase profile is a function of theheight map h(u,v) of a transmissive substrate:Φ(u,v) =2piλ(nλ −1)h(u,v), (4.8)where nλ is the refractive index of the substrate. We can design the height map of adiffractive lens such that the PSFs can fit our target PSFs in the `1 sense. We choosethe `1 norm instead of the least squared error term because penalizing the absolutevalue is more robust to outliers, which in our case are the sparse high frequencycomponents (e.g. glitch intensity) of PSFs. These high spatial frequencies arecaused by optimizing discrete wavelengths individually, which leads to coherentinterference patterns. Note that we are targeting imaging with incoherent light,such that the real PSFs are always smooth in practice. We simulate incoherence bylow pass filtering and outlier removal.The minimization problem is written as:hopt = argminh ∑λi∈Λwi ‖ pi(h)− t ‖1, (4.9)where we have omitted the spatial coordinates for brevity. Here pi(h) are the op-timized PSFs, and t is the wavelength-independent target PSF. The weights wi(wi > 0,∑wi = 1) are assigned to balance relative diffraction efficiencies amongwavelengths (see below). Note that the optimization uses a discrete set of designwavelengths λi which densely sample the target spectral range. This is possible be-cause PSFs vary smoothly with wavelength. The implementation of pi(h) followsdirectly from Eq. 4.3: we first calculate the amplitude PSFs by Fresnel diffractionpropagation and then take the magnitude squared.In this work, we consider only rotationally symmetric patterns. As a result,we can reduce the optimization to a 1D problem. Operating directly on height50Algorithm 1 Optimization on diffractive achromat1: k = 0, h0m = hinit , t = tinit ,w0i = winit ,v0m = vinit2: for iterations do . Repeat till convergence/termination3: for all seeds 1, ..,m do . Implement m seeds in-parallel4: hk+1m = hkm+ vkm . Update height profile5: hk+1m = OPT(hkm,hk+1m ) . Update local optimal6: h¯k+1 = OPT(hk+1m ) . Update global optimal7: vk+1m = vkm+ c1(h¯k+1m −hk+1m )+ c2(h¯k+1−hk+1m ) . Update velocityvector for each seed8: end for9: end for10: function OPT(h) . Store updated height profile11: for i = 1 to N do12: pi = pi(h) · f . Update PSFs13: end for14: hopt = minh ∑i wki ‖ pi− t ‖1 . Evaluate objective15: wk+1i =‖pi−t‖1∑i‖pi−t‖1 . Update weights16: return hopt17: end functionprofile is beneficial as fabrication constraints and other proximity effects can beincorporated into the model. The proposed method is summarized in Alg. 1. Wediscuss the algorithm in detail below.Target PSFs. In optical design, blur kernels are usually represented as Gaussiandistributions with different variances. However, in our design we will end up sacri-ficing resolution in the central wavelengths for improved resolution at both longerand shorter wavelengths, so that the final PSF is achromatic but not as sharp. Weexpect this process to introduce longer tails that are not represented well with asingle Gaussian.To seek an optimal distribution that is feasible with the current physical pro-file, we adaptively tune the target function. Specifically, after a few iterations weaverage the PSF distributions of all wavelengths and fit this average to a mixture51model of three (centered) Gaussians to represent the PSF:t = ∑j=1,2,3a jG (µ,σ j), (4.10)where a j are the weights for each component function and ∑3j=1 a j = 1. The co-efficients and parameters of the fitting model are then tuned (e.g. σ j are shrunkby a few pixels) to generate new target PSFs with a sharper distribution. Via thisupdate strategy, we only need to initialize the target PSFs once based on prelimi-nary simulation. Repeatedly fitting the averaged distribution assures that we haveconsidered the focusing contribution of all wavelengths.Despite its achromatic focusing constraint, this PSF representation helps main-tain a relatively sharp peak in the PSFs so that high-frequency features can bepreserved. This is particularly beneficial for deconvolution under the challenge oflarge kernels.Initialization. The optimization begins by calculating an initial guess for the start-ing height profile. A purely random height profile leaves our input too far awayfrom the solution. We could start from a Fresnel phase plate designed for a sin-gle, central wavelength (e.g. 550nm for the visible spectrum). However, we foundthat a better choice is to use a composite of multiple zone plates. When optimiz-ing for N discrete wavelengths, we divide the aperture into N rings of equal area.Within each ring, we initialize the height field to a Fresnel phase plate for a specificwavelength λi at the target focal length (see Fig. 4.1).Adaptive weights. The final goal of our optimization is to uniformly assign thefocal contribution of each wavelength. However, during the optimization we stilladaptively tune the weights wi in Eq. 4.9 for selected wavelengths according to theirdeviations of fitting errors. Specially, if the current design suffers from a weakeroptical power at one wavelength, which is explicitly reflected by larger fitting error,the weight on this wavelength in the cost function will be adjusted in the nextiteration following a simple rule as wk+1i =‖pi−t‖1∑i‖pi−t‖1 . We see that all weights wiwill be updated toward nearly identical value if the optimization approaches the52Figure 4.1: The initial height profile for the optimization is a mixture of sub-regions screened from Fresnel lenses at the same focal length for dif-ferent wavelengths. The whole area is divided into N subregions for Ndifferent wavelengths. The resulting initial height profile is the super-position of such subregions.optimal.Low-pass filtering. The resulting PSFs after diffraction propagation contain nar-row spectral peaks and valleys (Fig. 4.2). These artifacts are due to two typesof discretization. First, both the DOE plane and the image plane are representedas point samples, which introduces high spatial frequencies that in reality get aver-aged out by integrating over finite pixel areas. The second, and maybe more impor-tant, effect is that our simulation is also discretized along the wavelength direction,which effectively treats the light as coherent and introduces artificial interferencepatterns that are not present in real-world broad-band imaging scenarios. As a re-sult of this analysis, we treat the spectral peaks and valleys as outliers that we filterout by applying a blur along the spectral dimension. Our experience indicates thisfiltering benefits the robustness and convergence speed of the optimization.4.2.2 Stochastic optimization algorithmWe choose the Particle Swarm Optimization (PSO) algorithm to solve our opti-mization [32]. The advantage of the PSO algorithm is its high computing efficiencycompared with other stochastic algorithms, e.g. Genetic Algorithm etc. By imple-53menting a series of seeds (i.e. m seeds in Alg. 1) in parallel at each iteration, theheight profile update is more robust and faster to converge.At each iteration, the seeds update their height profiles and velocities of currentdesign by tracking the optimal solution of their own h and that of the group h¯following a strategy hk+1 = hk+vk, where vk is a velocity vector indicated in line 7of Alg. 1. The two weights c1,c2 are randomly assigned between (0,1). We furtherset a constraint in the implementation that |vm| ≤ 0.25hmax.The idea is to assume individual seeds can evolve according to the informa-tion gathered from their own experience and that of the group, so that the focalpower change for each individual wavelength is not drastic and purely random.This update strategy is beneficial to avoid falling into a local minimum, as well asleveraging parallelism.We borrow the from-coarse-to-fine strategy from multi-scale optimization todivide the N wavelengths to be optimized into several scales. For instance, we startby optimizing N1 = 9 sampled wavelengths in the spectrum. Once the optimizationfor this scale has converged, we increase the number of wavelengths to a secondlevel N2 = 15, and N3 = 29 eventually. We can accept an 10 nm interval in thewavelength sampling to already approximate achromaticity among wavelengths inthe full visible spectrum.Figure 4.2 shows the comparison of simulated PSF cross-sections for threeselected wavelengths λ1 = 650nm, λ2 = 550 nm and λ3 = 450 nm in the full spec-trum. We compare a single wavelength Fresnel phase plate (left) with the multi-ring initialization (center) and the final, optimized result (right). Our algorithmsacrifices the performance of the central wavelength to compromise the spectralfocal contributions. It’s worth noting that the shown PSFs correspond to singlewavelengths, which explains the high spatial frequencies, which are interferencepatterns due to coherence. These patterns average out for incoherent illumination(see Sec. 4.4 and Fig. 4.10).4.3 Image reconstructionIn this section we present the algorithms for solving the inverse problem to recoversharp and high-fidelity color images using a cross-scale prior in a two-step fashion.54Figure 4.2: Cross sections of normalized PSFs for a regular Fresnel lens(left), the initial guess (center), and the optimized diffractive achromat(right) at three selected typical wavelengths λ1 = 650 nm, λ2 = 550nm and λ3 = 450 nm. Note that here we optimize for 29 wavelengthsfrom 410 nm to 690 nm with an interval of 10 nm. We sacrifice theperformance at the central wavelength to equalize PSF distributions.4.3.1 Optimization methodIn order to recover latent images, our approach seeks the solution to Eq. 4.7 bysolving the following minimization problem (as presented in Sec. 2.3 of Chapter 2):ic = argminicµc2‖bc−Kcic‖22+Γ(ic), (4.11)where the first term is a standard least-square data fitting term, µc is the weightof data fitting term for channel c = 1,2,3. The second term Γ(ic) is a regulariza-tion term that enforces natural image priors in the solution, as explained in detailsbelow.We investigate an efficient non-blind deconvolution accounting for the follow-ing properties of our diffractive achromat design. First, despite the large PSF sizeof our lens, they show a preserved central intensity peak as well as quasi-uniformintensity distribution for all the color channels, therefore we do not need additionalcolor priors in our problem. Second, inspired by the image-pyramid strategy, wesolve our problem at two scales, enforcing similar gradient distributions betweenscales. Intuitively, the edges in natural images always exist at the same locations55and are barely affected when the image is downsampled to a lower scale. However,downsampling an image would lead to an improved signal-to-noise ratio, whichhelps improving the conditioning of the problem. The latter point can be particu-larly beneficial in our case.Fast deconvolution and denoising at downsampled scale. We propose to imple-ment the first step deconvolution on a downsampled image, for instance, half thesize of the original image, to deblur large edges and remove strong color corruptednoise. By defining the regularization term Γ(ic) in this scale, the cost function fora single channel in Eq. 4.11 is reformatted as:idc = argminidc µ‖bdc −Kidc‖22+β‖Didc‖22⇔ idopt =F−1(µF (K)∗F(bdc)µF (K)∗F (K)+βF (D)∗F (D)),(4.12)where F (·) represents the Fourier transform, F−1 (·) is its inverse. The super-script ∗ indicates the complex conjugate operation, D is the first-order derivativefilter matrix, and µ , β are the respective weights for each term. The superscriptd denotes that all images are in the downsampled scale. This quadratic problemleads to a closed-form solution in frequency domain, such that we can directly usefast inversion to recover the sharp image idc at downsampled scale. In practice, wesuggest to apply an additional denoising solver here at this scale if the capturedimage suffers from strong noise.Cross-scale deconvolution at full scale. In the second step, we apply a cross-scale prior in our regularization term, which borrows the relatively sharp and de-noised edge information from the upsampled image of the 1st step’s result to benefitthe deconvolution at full scale. Our cross-scale prior is inspired by cross-channelprior [57] and the multi-scale deconvolution scheme [174]. It is reasonable to as-sume large edges and shapes will be located where they are at both the upsampledversion of the downsampled scale image and the original full scale image:∇ic ≈ ∇isc ⇔ Dic ≈ Disc, (4.13)56where ic is the latent image at full scale, and isc is the upsampled version of thedeconvolved image idc in the first step using a bicubic sampling scheme.Note that strong noise in the original blurry image has been smoothed andmostly removed during the processing at downsampled scale, then the second stepcan be run with relatively weak regularization. Our approach differs from themulti-scale deconvolution in [174], which progressively refines details in multi-ple scales such that at each scale iterative residual deconvolution is necessary. Inour approach, we instead directly add the upsampled deconvolved image isc as anadditional prior term to transfer their edge information between two scales. Theremainder of our deconvolution problem still follows the same fast deconvolutionscheme described above. By tuning the weights of two prior terms, we can flexi-bly compromise on sharpness and smoothness of the recovered image in a simplescheme. Then, by rewriting Γ(ic) that includes the gradient prior as well as ourcross-scale prior, the cost function in Eq. 4.11 is reformulated as:ic = argminicµ2‖bc−Kic‖22+β‖Dic‖1+ γ‖Dic−Disc‖1. (4.14)Adding the cross-scale prior results in a non-linear optimization problem that canbe solved by introducing slack variables for the `1 term. Specially, we form theproximal operators [12] for the subproblems, thus turning the `1 terms into shrink-age operators. We define p=Dic as a slack variable, so that we have the followings:proxθ‖·‖1(p) = max(1− θ|p| ,0) pproxθ‖ · −α‖1(p) = max(1− θ|p−α| ,0) p+α,(4.15)where α =Disc. The proximal operators for their convex conjugates can be derivedfrom [12]. We then use a similar half-quadratic penalty scheme as in [79] to solveEq. (4.14).Specifically, Eq. 4.14 can be rewritten by introducing the variable p = Di thathas been denoted. Then, regarding one channel we have:iopt = argminiµ2‖b−Ki‖22+β‖p‖1+ γ‖p−Dis‖1+ρ2‖Di−p‖22, (4.16)57Figure 4.3: Comparisons for different combinations of deblurring step anddenoising step on synthetic dataset (top row) and real capture (bot-tom row), each with input blurry noised image, fast LUT deconvo-lution [79], direct deconvolution + MLP [133], their multi-scale mix-implementations, and our approach with cross-scale prior. For the toprow, the blurry image is synthesized using hyperspectral images with 29wavelengths, blurred by the kernel in Fig. 4.9, and added with σ = 0.005Gaussian white noise. The inset numbers indicate the peak signal-to-noise ratio (PSNR) and runtime on the synthetic image with 1.45megapixels size in Matlab IDE on a commercial PC. Ours is provedto be robust and efficient.where ρ is a weight that shall be gradually increased during the optimization, suchthat once the weight ρ is becoming very large, the modified objective Eq. 4.16 isequivalent to the original objective we have defined in Eq. 4.14. We set the valuesnamely integer powers of√2 between 1 and 256 in all our implementation. Thenwe divide the objective into two subproblems and apply alternating methods tosolve them iteratively. For the two quadratic terms, by fixing p we have a similarclosed-form solution as in the first step deconvolution. Next we fix i, and optimizethe two l1 norm terms that contain p, the proximal operators that have been definedin Eq. A.9 turn them to be pixel-wise shrinkage problems.4.3.2 Efficiency and robustness analysisIn the first step, the direct division in frequency domain is very fast, but may resultin the recovered image being strongly corrupted by noise. One can further apply amulti-layer perceptron (MLP) approach in [133] or similar fast method to denoisethe lowest scale before upsampling. Here we directly use the MLP-based denoiserwith off-line learned data in [133] at the downsampled scale. Although the ker-58nel of our lens differs from the Gaussian kernels with which the system has beentrained the results are still perceptually pleasing for our purpose.In the second step, we can recover the latent image efficiently at the full scaledue to the introduction of the cross-scale prior. We have compared our implemen-tation with recently reported non-blind deconvolution methods in Tab. 4.1. The lastcolumn in Tab. 4.1 indicates the case of a regular Fresnel lens with cross-channeldeconvolution. We find that all the results reconstructed from the diffractive achro-mat show much higher PSNR than those from a regular Fresnel lens. This validatesour design motivation that diffractive achromat preserves higher color fidelity thanconventional diffractive lens. Further, by running our two-step deconvolution, de-noising and cross-scale edge preservation, our results outperform existing methods.Table 4.1: Averaged PSNR comparisons of recovered images from differentdeconvolution schemes. The first 5 columns indicate the results of ourdiffractive achromat using respectively the deconvolution by Krishnan,Schuler, multi-scale Krishnan, and multi-scale Schuler+Krishnan (seeFig. 4.3), while the last column shows the result for a standard Fresnellens with a cross-channel prior.Method 1 2 3 4 Ours FresnelPSNR/dB 25.86 26.05 26.51 26.73 27.10 22.54At either scale, we do not rely on a color regularizer in the cost function, whichsignificantly lowers computing burden. Additionally, the high frequency compo-nents in the compromised PSF is usually of very low intensity. They are actu-ally mixed with additive noise, and are smoothed in the downsampling process.Although the theoretical size of PSFs can be considerably large, in practice, itdoesn’t degrade the image so much. In our experiment, for a 20 megapixel RGBimage and 200-pixel PSFs, the running time of our two-step algorithm in MatlabIDE is around 250 seconds on a PC with Intel Xeon i7@2.70Hz CPU. Furtherperformance improvements could be achieved through GPU optimization and par-allelization.594.4 ImplementationWe show our implementation of the proposed achromatic diffractive imaging method,including the prototype design, fabrication and a number of experimental resultsfor synthetic data as well as real indoor and outdoor scenes at full sensor resolu-tion. For figures some figures we present only cropped regions to highlight detailstructure.Prototype. We designed two types of diffractive lenses, our diffractive achromatand a conventional diffractive lens with the same optical parameters for compari-son. The focal length is designed at f = 100 mm with an aperture diameter 8 mmfor both cases. The conventional diffractive lens is designed at central wavelength550nm. Our diffractive achromat is optimized for wavelengths from 410 nm to 690nm with a 10 nm sampling interval. Both lenses are attached on a DSLR camerabody (Canon 70D) that has 5,740 × 3,648 pixels with the pixel pitch of 4.1 µm.We fabricated the designed diffractive achromat using multi-level photolithographytechniques. Diffractive lenses with 2D micro-structures approximated by 16-levelachieve a theoretical diffraction efficiency of up to 95%, and increasing the num-ber of levels to 32 yields almost no improvement [37]. See Appendix chapter A formore details on the fabrication process.Experimental results. We show in Fig. 4.3 the comparison results using differ-ent deblurring and denoising methods, including single scale fast LUT deconvolu-tion [79], single scale MLP-based deconvolution [133], two-scale fast LUT decon-volution, one scale MLP-based deconvolution followed by one scale fast LUT de-convolution, and our two-step cross-scale deconvolution. Our algorithm producessharp, high color fidelity images in the full visible spectrum due to the achromaticdesign. We run our experiments for a hyperspectral image database with 50 im-ages and the results yield an averaged PSNR 26.2 dB with our proposed algorithm,even the worst result still stays above 25.0 dB (See Tab. 4.1). Figure 4.4 showssome synthetic results of natural scenes from the hyperspectral datasets in [17] and[142].Figures 4.5 and 4.6 show the experimental results captured using our diffractive60Figure 4.4: Illustration of synthetic results with blurred (top row) and de-blurred (bottom) results of selected scenes, with σ = 0.005 Gaussianwhite noise added, the left 3 pairs are synthesized with 29 wavelengths’hyperspectral images, while the right 2 pairs 71 wavelengths. The quan-titative evaluations on full dataset are also provided in the supplementof [118].achromat. We have presented diverse natural scenes, including indoor, outdoor,rich color, high reflection, etc.. The results show that our method is close to spatialand depth invariant achromatic aberration within a long range. Refer to the captionof each figure for scene details.Figure 4.5: Blurred (top) and deblurred (bottom) results of real capturedscenes. All four scenes are captured at different depths with single ex-posure and no gain, under indoor artificial and mixed illumination (lefttwo pairs), and natural sunlight illumination (right two pairs), using asingle 0.5 mm ultrathin diffractive achromat we fabricate. Note that weroughly use a white light source attached with a pinhole to calibrate thePSF only at one depth (2 m), and we use it for all deconvolutions.61Figure 4.6: Blurred and deblurred results of outdoor scene with a large depthvariance and reflection feature (left pair), desktop stuff with rich col-ors (center pair) and natural human face (right pair). The same PSFcalibrated for Fig. 4.5 is used here.Compared to synthetic results, the captured results suffer from an additionalhaze effect, which degrades the image quality. We identify several sources for thesedeviations: First, the discrete height profile and limited spatial resolution of ourphotolithography process reduces the diffraction efficiency of the prototype. Thiscould be alleviated by moving to an electron-beam lithography process. Second,the engineering errors derived from the custom optical mounts, e.g. custom holderand aperture result in some light leaks in the camera. Finally, our prototype isdesigned to be optimal for a spectral range from 410nm to 690nm, while the sensormay have a wider spectral response.4.5 Evaluation and discussionWe have demonstrated our diffractive achromat is able to image natural scenes withcompetitive resolution and color fidelity. In this section, we analyze the imagingperformance from the perspectives of spatial resolution, off-axis behavior and colorperformance. The potential applications and limitations are discussed as well.Resolution measurement. We evaluate the resolution of our achromatic lens bytaking an image of the ISO 12233 resolution chart, shown in Fig. 4.7. The cap-tured image is very blurry due to the large blur kernel from our achromatic lens.62Figure 4.7: Blurred (left) and deblurred (right) results of capturing a standardresolution chart image which is projected by a projector onto a whiteplane. The capture distance is around 1.8m. The same PSF estimationused in Fig. 4.6 is applied here.High-frequency features, such as small edge patterns shown in the close-ups, are in-distinguishable, After the deconvolution step, the recovered image preserves mostof the low-frequency and mid-frequency features. Image contrast is also improved.The resolution comparison with a pure refractive lens is provided in the following.Off-axis behavior. In addition to the benefit of ultrathin and lightweight structure,our diffractive achromat exhibits lower off-axis distortion than that of a refractivelens with the same focal length, as illustrated in Fig. 4.8. Compared to the simplerefractive lens that has highly spatially variant PSFs across the field of view, ourdiffractive achromat exhibits smaller field curvature, thereby the resulting PSFs arealmost uniform across the image plane. This property simplifies PSF calibrationthat otherwise is time-consuming or even impractical for refractive lenses. Specifi-cally, after the deconvolution step, the off-axis image patch of a hybrid lens exhibitssharper edges than that of a pure refractive lens (see the zoomed-in patches). TheMTFs estimated from the slant edge method [130] are presented in the bottom-right of Fig. 4.8. We see that the diffractive achromat results in a good compromisebetween on-axis and off-axis performance. The computational burden for splittingimage patches to account for the spatial variant blur kernels is also eliminated andis further beneficial to the deconvolution step.Owing to the spatial invariant PSFs of our diffractive achromat, we can alsointroduce a hybrid refractive-diffractive lens design as in [117]. Our hybrid de-sign concept differs from conventional hybrid designs in that we do not leverage63Figure 4.8: On-axis and off-axis behavior comparison of an achromatic re-fractive lens (top-left) and a hybrid diffractive-refractive lens (bottom-left). From the blurred and deblurred patch pairs presented in the righthand side, we observe that embedding our DOE design in a lens ex-hibits better spatial uniformity, despite the residual aberration. Herewe assume within each selected patch the PSF is locally invariant. Ac-cordingly, the MTFs estimated from gray-scale slant edges inside eachdeblurred patch are provided (bottom-right). The auxiliary refractiveoptics used are Thorlabs achromatic lenses, with focal length 50 mm,100 mm, and thickness 8.7 mm, 4.7 mm, respectively. The equivalentfocal lengths are all 50 mm and full field of view is around 30◦.the negative chromatic aberration of DOEs to compensate for the positive chro-matic aberration of refractive lenses [100]. One can combine our diffractive achro-mat with any off-the-shelf refractive lens to assemble a hybrid lens that has im-proved spatial uniformity compared to purely refractive designs. See supplemen-tary document for the extensional experimental results of an achromatic diffractive-refractive hybrid lens.Color performance. The advantage of our diffractive achromat is to generatespectrally invariant PSFs, which is usually neglected by conventional diffractiveimaging methods. An ideal diffractive imaging device should preserve high spatialfrequencies across all color channels, similar to a refractive optical system. Our64Figure 4.9: Color performance comparison of conventional diffractive lensand our achromatic diffractive lens. Conventional diffractive lens hasa sharp green channel but severely blurred blue and red channels. Ourachromatic lens balances three channels to show averaged performancefor all. The MTF plots for our design (solid) are closer to each other,compared to those of the conventional diffractive lens (dashed). Theinsets show the respective color PSFs of the two lenses.proposed optimization method achieves this for all the three color channels. Weshow the Modulation Transfer Functions (MTFs) for each channel to illustrate thebenefit of our design in Fig. 4.9. Compared with a conventional diffractive lens,our design shows balanced performance in three channels indicated by the MTFcurves getting close to each other (solid plots), while the MTFs for the conven-tional diffractive lens are separated far away from each other (dashed plots). Thiscan be more clearly seen from the captured color PSFs for both lenses. The con-ventional diffractive lens has a peak for the green channel and the other channelsare of very low intensity, so the PSF looks green. For our diffractive achromat,since we optimize for 29 wavelengths to equalize the intensity distribution, thecolor PSF looks much more natural and closer to that of the refractive lens. Themeasured results of PSFs are presented in Fig. 4.10, from which we see that thequasi-uniform spectral PSF behavior has been established.Limitations. Our prototype suffers from a few shortcomings, which can be tracedback to limitations of the current manufacturing process. Like other DOEs with65Figure 4.10: The measured PSFs on a sRGB color sensor for 5 selected spec-trum bands (left), and the PSF we have calibrated for all above de-convolutions (right-most). Note that for individual measurement, abandpass filter with a 40 nm FWHM is attached in front of our lens.discrete levels of surface relief, our achromatic prototype cannot achieve 100%diffraction efficiency for all the wavelengths. The spread of the lost energy re-sults in slightly foggy appearance in the captured image, which is still difficult toeliminate especially in high contrast scenes.Our prototype is also limited by the resolution of the photolithography processthat we employ for fabrication. This results in a minimum feature size of 1 µm,and places a limit on the aperture sizes and focal lengths that are possible with thisprocess. We note that the resolution could be improved by two orders of magnitudeswitching to an electron-beam lithography process, which, however, was beyondthe scope of this work.Finally, we also note that there is a small residual wavelength dependency leftin the optimized design of the diffractive achromat. For imaging in natural environ-ments with broad-band, incoherent light, we have demonstrated that this is not anissue. However, for partially coherent or fully coherent light, the PSF will containinterference patterns with high spatial frequencies that are sensitive to small shiftsin wavelength. The deconvolution method will not be able to restore a high-qualityimage in this scenario.Potential applications and future work. We have shown in this paper that diffrac-tive achromats are viable for full spectrum imaging in the visible domain. How-66ever, our method is not limited to applications in the visible band. Diffractive lensesare very promising for imaging in the ultraviolet or far infrared spectrum, whererefractive lenses are not able to transmit the desired wavelengths with high effi-ciency [71, 161]. The ultrathin structure of diffractive optics drastically improvestransmission for these wavelengths, while simultaneously reducing weight.Moreover, it is possible to use our approach to design custom optics for spe-cialized purposes that require imaging of multiple discrete wavelengths. By op-timizing only for these wavelengths while neglecting others, better PSFs can beachieved than with a broad band optimization. Due to its flat field property, ourdiffractive achromat design is an alternative option to be combined with off-the-shelf lenses to correct off-axis aberrations flexibly such that the number of lensesin an imaging system can be significantly reduced without decreases in image qual-ity [2]. Additional work on how to apply the known kernel distributions to benefitthe deconvolution process can also be interesting.4.6 ConclusionWe have proposed a novel achromatic diffractive imaging method that bridgesdiffractive optical elements and computational algorithms to build lightweight andthin optics for the full visible spectrum. By introducing optimization to the de-sign of diffractive optics, we develop a diffractive achromat that trades off colorfidelity and overall spatial resolution. The residual aberrations resulted from thiscompromise are tackled by employing a two-step image deconvolution with ourproposed cross-scale prior. The algorithm includes a downsampled scale fast de-convolution and denoising step, and a full scale cross-scale deconvolution. Bothsteps are implemented without additional color priors to rapidly recover high qual-ity color. We envision our method offers an opportunity for the joint design ofultrathin diffractive achromats and computational algorithms, further to boost thepotential application of compact imaging devices in many broadband illuminationscenarios.67Chapter 5Multi-focal Sweep ImagingThe previous chapter addresses design in spectral domain. Starting from this chap-ter we address encodability in spatial domain to enable extra functionalities. A firstapplication example is to extend depth-of-field (DOF) of conventional photogra-phy. Depth-dependent defocus results in a limited depth-of-field in consumer-levelcameras. Computational imaging provides alternative solutions to resolve all-in-focus images with the assistance of designed optics and algorithms.In this chapter, we extend the concept of focal sweep from refractive optics todiffractive optics, where we fuse multiple focal powers onto one single element.In contrast to state-of-the-art sweep models, ours can generate better-conditionedpoint spread function (PSF) distributions along the expected depth range with dras-tically shortened sweep distance. Further, by encoding axially asymmetric PSFssubject to color channels, and then sharing sharp information across channels, wepreserve details as well as color fidelity. We prototype two diffractive imagingsystems that work in the monochromatic and RGB color domain. Experimental re-sults indicate that depth-of-field can be significantly extended with fewer artifactsremaining after the deconvolution.5.1 IntroductionExtending depth-of-field (DOF) is an exciting research direction in computationalimaging [25, 153], particularly for embedded cameras where a large numerical68aperture (aka. a small f -number) is necessary to ensure high light throughput. Re-cent advances seek to design optics in combination with post-processing algorithmsto either preserve more information or enable extra functionality by reducing thecomplexity of lenses.Work on this problem ranges from capturing the entire light field [29, 109] toengineering point spread function (PSF) shapes [87, 181]. Using the prior knowl-edge on the mapping relation between kernel shapes and scene depths, one canrecover all-in-focus images. In addition to engineering PSF shapes on the pupilplane of a lens, another advance is to apply sweep-type solutions such as spatialfocal sweep or focal stack cameras [63, 81, 101].The focal sweeps facilitate nearly depth invariant blur kernel by capturing anintegrated PSF over a time duration, and then apply a deconvolution step to removethe residual blur effect [81, 172]. Sweeping reduces calibration requirements ofdepth-variant PSFs in the capture. This strategy has been applied in not only imag-ing domain but projection display domain [68]. Despite much research, auxiliarymechanics is usually required to sweep either the optics or the sensor for a physicaldistance.One common fact that has not been addressed by state-of-the-art sweep typecameras is that these systems rely on sweeping complex refractive optics. Planaroptics, like diffractive optical elements or metasurface lenses, have recently beenproven effective to shrink the camera lenses in both weight and thickness [43, 117].Although this advantage is prominent for sweep configurations, a regular Fresnellens still requires a considerably large sweep distance as its refractive counterpart.We note that using DOEs as imaging lenses provides flexibility to create mul-tiple foci with only one single element [118]. Despite much research in optics onmulti-focal lenses for ophthalmic and clinical applications [69, 73, 163], existingconsumer-level cameras barely use this kind of design. Theoretically by enablingmultiple focal powers subject to desirable depths, it is possible to reduce the sweepdistance as well as to achieve better conditioned integrating imaging (see Sec. 5.2).In this work, we present the following contributions:• We introduce a multi-focal sweep imaging system for extending depth-of-field from one aggregated image that incorporates optical design and post-69capture image reconstruction.• We propose a diffractive lens design that is fused with multiple focal pow-ers subject to two aspects: expected depth-of-field, and three color channelfidelity. The better-conditioned kernel after sweeping integration leads to anefficient deconvolution to resolve all-in-focus images. Moreover, the colorfidelity is preserved by enforcing cross-channel information sharing.• We present two prototype lenses to validate the concept, with sweeping ultra-thin tri-focal and novem-focal (aka. 9) diffractive lenses. We test our decon-volution on different scenarios with large depth variance. The results exhibitvisually pleasing quality, especially in terms of resolving all-in-focus imageswhile preserving color fidelity and suppressing edging artifacts.5.2 Optical design5.2.1 Multi-focal diffractive lensWe start by investigating the ability of a multi-focal lens on decreasing the sweepdistance of focal sweep imaging. Using geometry optics, the relationship amongsweep distance s, foci number N, focal length difference M f and focused depthrange (on image end) l can be cast as follows:[− l2,l2]⊆N−1⋃n=0[− s2− 2n− (N−1)2(N−1) M f ,s2− 2n− (N−1)2(N−1) M f ]. (5.1)If we consider a lens with the focal length of 50 mm to cover the focused depthrange (on object end) from 1.5 m to 9 m (l = 1.4 mm), and the sweep distance is0.5 mm, then the focal length difference M f should be beyond 0.96 mm and thefoci number N should be at least 3. Subsequently, we can choose M f = 1 mm andN = 3. This means the above requirements can be realized by using a tri-focal lenswith the focal length of 49.5 mm, 50 mm, 50.5 mm, respectively.As shown in Fig. 5.1, it is clear that for a tri-focal lens, due to the approximatelyperiodic distribution of the focal plane along an expected depth range (the green,blue, and red curves), we only need to sweep the image plane from 50.75 mm to70Figure 5.1: Comparison of focused distances subject to object distance andfocal length under the assumption of thin lens model (the math deriva-tion is given). The three color curves visualize the relations of usinglenses with the focal length of 49.5 mm, 50.0 mm, 50.5 mm, respec-tively. s1 and s2 represent the sweep distances needed for a tri-focal lensand a one-focal lens, respectively.51.25 mm (s = 0.5 mm) to cover the focused depth range. However, for a lens withsingle focal length, we need to sweep the image plane from 50.3 mm to 51.7 mm(s = l = 1.4 mm) to cover the same depth range from 1.5m to 9m (the green centercurve). This matches the derivation of Eq. 5.1, indicating that the sweep distancecan be drastically shortened by introducing multi-focal designs.We note that the sweep distance derived above is the minimum sweep distance.In practice, we choose to use a relatively larger sweep distance, e.g. 0.8 mm inthe aforementioned scenario. This is reasonable considering the different defocuseffect of each object plane within the range of the minimum sweep distance. Thefinal kernel that is integrated over a sequence of more uniform depth PSFs leads toa better-conditioned deconvolution.We generate the multi-focal lens by fusing multiple Fresnel lenses onto onesingle element. As mentioned above, we design two lenses. First, we divide theaperture into three rings of equal area. Thus, the monochromatic design is a radialmixture of subregions screened from Fresnel lenses at the wavelength of 550 nm forthree different focal lengths, which we call tri-focal lens. Similarly, the RGB colordesign is an axially asymmetric mixture of three aforementioned monochromatic71Figure 5.2: Comparison of synthetic PSF behavior of sweeping a regularFresnel lens (top) and our tri-focal lens (center and bottom) subject totarget depths. Notice that this design aims for monochromatic imagingthat is integrated over a spectrum of 10nm FWHM. The axes of eachsubfigure represent the size with a pixel pitch of 5µm. The normalizedcross-sections (right-most) indicate that our multi-focal sweep designsexhibit less variance (quantized) regarding PSF distributions. We sacri-fice the peak intensity at the central depth to minimize the variance ofPSF distributions along the full depth range.designs subject to three spectra, namely 640 nm, 550 nm, 460 nm, which we callnovem-focal lens (novem means nine in Latin). The graph fusion schemes and themicroscope images of our prototype lenses are shown top of Fig. 5.3.5.2.2 PSF behaviorsFigure 5.2 visualizes the synthetic PSF behaviors of a regular Fresnel lens and ourtri-focal lens, swept along a distance of 0.8 mm (top and center group) and 0.5mm (bottom group). We see that although none of its PSFs is highly focused, ourtri-focal lens exhibits less variance in terms of the size of peripheral energy distri-bution along the full depth range. This more depth-invariant blur makes it possibleto deconvolve the full image using only one calibrated PSF. We also note that PSFsbecome more depth-invariant when increasing the sweep distance slightly (centerFig. 5.2). This can be justified by the provided quantitative values as well. We willfurther explain our choice in the experiments.Figure 5.3 visualizes the real PSF behaviors of our two multi-focal lensesswept along a distance of 0.5 mm. Concerning our novem-focal lens (shown right72Figure 5.3: Diagrams of graph fusion schemes (top) and cropped microscopeimages of fabricated lenses as well as the experimental PSFs behavior ofsweeping a tri-focal lens (bottom-left) and a novem-focal lens (bottom-right). This experiment aims for RGB color imaging that is integratedon a RGB Bayer sensor.Fig. 5.3), despite the relatively small variance regarding the size of peripheral en-ergy distribution in different channels, the PSFs exhibit axially asymmetric shapewith high-frequency components. As the high-frequency components vary in spa-tial distribution from channel to channel, it is possible to recover that with sharedinformation from different channels.5.3 Image reconstructionWe next derive the image formation model for the focal sweep, as well as presentthe optimization for image reconstruction.5.3.1 Sweeping image formationThe defocus effect can be formulated as a latent image convolved with a blur ker-nel. Similar as the previous chapters, we write the recorded image in channel c invector-matrix format as:bc = Kcic+n, (5.2)where for a channel c, b, K, i, n are the captured image, convolution matrix, sharplatent image and additive noise in the capture procedure, respectively.73Regarding a sweep imaging system with a diffractive lens, Kc can be derivedfrom the PSF Pc integrated over the depth range ∆ and spectrum Λ as follows:Pc(x,y) =∫Λ∫∆Qc(λ ) · (P(x,y,z;λ ))dλdz, (5.3)where P(x,y,z;λ ) is the spatial spectral variant PSF describing the aberrations oflens, which is a function of both spatial position (x,y,z) and spectral component λ .Qc represents the sensor response, which can be reasonably assumed as the con-stant when used under narrowband scenario. As aforementioned, after sweepingintegration the PSF Pc is approximately depth invariant.5.3.2 Optimization methodTo resolve all-in-focus images, we formulate the inverse problem of Eq. 5.2 as anoptimization problem containing a least-squares data fitting term and a collectionof priors that regularize the reconstruction.Deconvolution on individual channels. Regarding the deconvolution of an indi-vidual channel, which is also the application scenario of monochromatic imaging,the prior term Γ(ic) is a total variation prior (i.e. l1-norm on gradients that arederived from multiplying matrix D). The optimization becomes as follows:idc = argminicµc2‖bc−Kcic‖22+‖Dic‖1. (5.4)We can directly use the Split Bregman method [51] to efficiently solve Eq. 5.4. Atrick here is to assign a slightly larger weight of µc so as to yield the deconvolvedresult idc exhibiting sharp edges and features. The intermediate resolved imagesserve as references for further cross-channel processing.Cross-channel regularization. The cross-channel regularization follows closelythe recent work [61]. This is realized by enforcing the gradient information tobe consistent among different color channels. With respect to the color multi-focal sweep scenario, ours differs from state-of-the-art methods in that there is nospecific sharp channel set as a reference. In our case, the images in three channels74shall all serve as references, since the color PSF exhibits different behaviors. Thatsaid, although none of the three channels is sufficiently sharp before processing,each channel shall preserve the details of some sense to the recovery of images inothers. The optimization then becomes as follows:ic = argminicα2‖bc−Kic‖22+β‖Dic‖1+∑m6=lγ‖Dic−Didc‖1, (5.5)where α,β ,γ are the tunable weights for each term.Specially, Eq. 5.5 can be solved by introducing slack variables for the l1 termand then using a similar solver scheme as in [118], which is detailed in Appendixchapter A as well. Although the deblurred images of each individual channel(Eq. 5.4) may suffer from sensor noise, most edges and features can be robustly re-covered from the cross-channel information sharing. These roughly deblurred im-ages idc are used as reference channel images in the cross-channel terms (Eq. 5.5)to iteratively recover the three channel images.5.3.3 Efficiency and robustness analysisWe note that the cross-channel regularizer makes the optimization problem morecomplex and non-linear and the resolved results could be highly dependent on thequality of the reference channel. However, we manage to gain reasonable goodresults with a reasonable amount of tuning effort.Using the color PSFs derived from two real prototype lenses, we have imple-mented simulations on a number of test images (BSDS500 dataset [96]). Extra0.5% Gaussian noise is added. The comparison results are illustrated in Tab. 5.1.Additional visualization results are presented later. With respect to the tri-focallens, we enforce cross-channel sharing only from the green channel image to theother relatively blurred red and blue channel images. With respect to the novem-focal lens, we enforce cross-channel sharing between all three channel images. Forthe latter, the averaged run time for one image with the size of 1,384 by 1,036 pix-els is around 7 seconds on Matlab IDE run on a laptop PC with 2.4 GHz CPU and16 GB ram.We have the two observations. First, enforcing cross-channel information shar-ing contributes to resolving higher quality images in both scenarios. Further, en-75abling graph fusion subject to colors in addition explores cross-channel informa-tion sharing to preserve higher color fidelity.Table 5.1: Comparisons of synthetic image reconstruction with PSNR aver-aged over 100 dataset images. 1 indicates the tri-focal lens and 2 indi-cates the novem-focal lens.w/o. cross. 1 w/. cross. 1 w/o. cross. 2 w/. cross. 220.392 23.299 23.523 24.6255.4 Implementation and discussionIn this section, we present selected experimental results that validate the effective-ness of the proposed approach. Before detailing these experiments, we start byintroducing the parameters of our prototype lenses.5.4.1 Prototype parametersWe have designed two types of multi-focal diffractive lenses, one for monochro-matic imaging while the other for RGB color imaging. The aperture diameter is8mm for both designs. The monochromatic one is designed at the central wave-length of 550nm and fused with 3 Fresnel lens patterns subject to the focal lengthof 49.5 mm, 50.0 mm, 50.5 mm. The color one is designed with fusing afore-designed monochromatic patterns subject to the wavelength of 640 nm, 550 nm,460 nm. Both lenses are attached in front of a PointGrey sensor (GS3-14S5C) thathas 1,384 × 1,036 pixels with the pixel pitch of 6.45 µm. The exposure time is500 ms for lab setting scenes and 650 ms for office setting scenes, during when 0.5mm axial distance is swept. The experimental setup is illustrated in Fig. 5.4.We fabricated our designed lenses by repeatedly applying photolithographyand RIE techniques [118]. Again, refer to Appendix chapter A for fabrication de-tails. Note that in this chapter, we use 4pi phase modulation corresponding to 2.39µm etching depth on the wafer. The higher order diffraction benefits to yielding ashort focal length (aka. small f -number) design with the practical feature size ofstate-of-the-art fabrication techniques.76Figure 5.4: Photograph of our experimental setup. Left: a captured scenewith a large depth range; Right: the prototype lenses are mounted on aholder while the sensor is mounted on a controlled translation stage.5.4.2 Simulation resultsSimulation results of two standard images are presented in Fig. 5.5. From thezoomed-in insets, we observe that our axially asymmetric fusion design preserveshigher color fidelity than a regular symmetric multi-focal design, while its abilityto distinguish fine details is slightly traded off.5.4.3 Real world resultsThe real world results are presented in Fig. 5.6 and in Fig. 5.7. The depth range isset from 1.5 m to 3.5 m for the lab setting scenes (shown left of Fig. 5.4) and setfrom 2 m to 8 m for the office setting scenes. In particular, since the first proto-type aims for monochromatic imaging, the reconstructed results of green channelexhibit decent quality.We first set the green channel as the reference and use a cross-channel priorto deconvolve the images. As shown bottom row of Fig. 5.6, although exhibitingreasonable spatial resolution, its color fidelity is quite low. This is because naiveFresnel lenses suffer from severe chromatic aberration. A regular cross-channelprior is not sufficiently robust to preserve both spatial frequency and color fidelity.In contrast, the second prototype additionally favors axially asymmetric PSFssubject to three color channels. That said, we have a relatively high-intensity peak77Figure 5.5: Selected simulation results: (a) ground truth inputs and kernels;(b) degraded images blurred by corresponding kernels; (c) reconstruc-tion results using TV-based deconvolution on individual channels; (d)reconstruction results using deconvolution with TV and cross-channelregularization. The two color PSFs used to degrade the images are cal-ibrated from the two prototype lenses under white point light source il-lumination. In addition to the background noise in the calibrated PSFs,0.5% white noise is added. Zoom in for fine structures.with high-frequency long tails of PSF in each channel such that in the deconvolu-tion we can preserve color fidelity (shown in Fig. 5.7). However, limited by thedata bandwidth of the DOE, we have to trade off the spatial resolution in somesense. The overall image quality is still visually acceptable. Again, this workaims for extending DOF rather than naively pursuing spatial resolution. From thisperspective, despite the slight image contrast loss due to the fabrication and proto-type issues, our multi-focal lenses outperform off-the-shelf products, as shown inFig. 5.8. To achieve a competitive DOF behavior, one needs to drastically shrinkdown the aperture to at least f /12, which requires much longer exposure in practice.78Figure 5.6: Cropped regions of real world results (tri-focal lens). Top: de-graded inputs; Bottom: reconstruction results using deconvolution withTV and cross-channel regularizations. For experimental convenience,we capture a depth range from 1.5 m to 3.5 m for the left two scenesand from 2 m to 8 m for the right scene with a sweep distance of 0.5mm. We use the single calibrated PSF shown in Fig. 5.3 to deconvolveall images.Figure 5.7: Cropped regions of real world results (novem-focal lens). Top:degraded inputs; Bottom: reconstruction results using deconvolutionwith TV and cross-channel regularizations. We use the same experi-mental setting as in Fig. 5.6.79Figure 5.8: DOF comparison between our tri-focal lens (left) and a standardEF50 (Canon) refractive lens (right), both with a f -number 6.25. Thescene depth range is 1.5 m to 3.5 m, highlighted by different color rect-angles. We here extract the green channel for a fair comparison.5.4.4 DiscussionOn the optics end, the current fusion scheme of multiple foci is heuristic. Ourprototype lenses are just two effective designs. The optimal spatial distributionof PSFs may vary. Designing the fusion scheme in an intelligent way remains anopen but interesting direction. We anticipate learning strategies, like look-up tableor dictionary search, can be introduced to guide the design.Remaining artifacts like low image contrast and residual blur are due to severalengineering factors. Careful readers may observe from the results that a slight hor-izontal shift (2-pixel level) occurs when sweeping the lens. This is mainly becauseour sweeping axis is not strictly perpendicular to the sensor plane. The customizedlens holder and cover may inevitably introduce ambient light that amplifies thenoise. We also note that metamerism issues exists since our scope is not aiming forfull spectrum, thus slight color artifacts may remain when used under white lightillumination.We note that current DOEs with 16-level structure still suffer from a non-trivialdiffraction efficiency loss, especially for high diffraction order designs. This down-side can be observed as low image contrast and additional blur. Due to the inherentlimitation of feature size, it is challenging to create a diffractive lens with a verysmall f -number. We will switch to more advanced techniques like nano-imprintingor grayscale photolithography, to overcome these problems.80On the reconstruction end, the cross-channel regularization can be further ex-ploited. We anticipate there shall exist a better strategy to define reference channelsrather than enforcing the current two-step deconvolution scheme. Additional de-noising solver can be added to obtain better visualization.On application end, the narrowband or single color design is promising insurveillance scenarios where a large FOV and a large DOF are strongly acknowl-edged. In addition, a wide range of depth sensors with active illumination is excel-lent candidates for embedding our multi-focal lenses. Active illumination ensuresthat fusing only a few wavelength designs is reasonable. Thus, one can yield greatdesign freedom to enable the ability to extend depth of field (DOF).5.5 ConclusionWe have proposed a computational imaging approach that jointly considers sweep-ing diffractive design and image reconstruction algorithms, and demonstrated thepracticality of extending depth-of-field with compact lenses. Benefiting from thedesign flexibility of diffractive optics, the proposed design significantly shortensthe required sweep distance meanwhile it exhibits a better conditioned depth-invariantkernel behavior.Moreover, color fidelity is preserved by fusing spectral variant PSF behaviorsin the diffractive lens design and enforcing cross-channel regularization in the de-convolution. We have validated the effectiveness and robustness of our methodwith a variety of captured scenes. Although current experimental results sufferfrom the problems of slight blurry and low contrast that can be resolved via a rea-sonable amount of engineering effort, our approach shall be an effective solutionto extend depth-of-field especially in situations where thin and lightweight opticsis expected.81Chapter 6Encoded Diffractive LensesThe previous chapter addresses an heuristic design in spatial domain to derive thediffractive phase patterns subject to describable imaging scenarios. We note thatthe design scheme can actually be improved to a more specific and intelligent man-ner, for instance, by involving formal optimization.In this chapter we introduce numerically optimized encoded phase masks inwhich different optical parameters such as focus or zoom can be accessed throughchanges in the mechanical alignment of an ultrathin stack of two or more phasemasks (DOEs). Our encoded diffractive designs are combined with a new compu-tational approach for self-calibrating imaging (i.e. blind deconvolution) that canrestore high quality images several orders of magnitude faster than the state-of-the-art without pre-calibration of the optical system. This co-design of optics andcomputation enables tunable, full-spectrum imaging using thin diffractive optics.6.1 IntroductionIn conventional imaging systems great efforts have been made to combat aberra-tions by designing an increasingly large number of refractive or reflective lenses [37,76]. The number of elements in a commercial digital single-lens reflex (DSLR)camera, for example, can easily reach more than a dozen. As another example,groups of doublets or triplets are necessary in high quality microscopes to achievediffraction-limited performance. In addition, more complicated imaging systems,82e.g. zoom lenses, frequently require additional longitudinally moving groups tochange the focal lengths in a specific range [144, 171]. Furthermore, special glasswith extraordinary dispersion properties is incorporated into many designs to min-imize primary and higher order chromatic aberrations. Consequently, conventionalimaging systems are, in most cases, bulky, heavy, inflexible, and costly.A promising avenue for reducing the size and complexity of conventional lensesis to replace complex components with simpler and smaller elements in combina-tion with computational techniques. The aberrations in a simplified optical systemcan be designed for easy computational removal, for example by favoring pointspread functions (PSFs) that preserve high spatial frequencies. Diffractive opticalelements (DOEs) in particular are an interesting replacement for complex imagingsystems. DOEs can be fabricated on ultra-thin plates and are very flexible in mod-ulating light, thus delivering greater freedom in design parameters. One drawback,however, is that they introduce strong chromatic aberrations [112] which limit theirapplication in color imaging. This has been addressed in the previous chapters.Computational imaging techniques have proven to be able to digitally correct opti-cal aberrations [131, 177] including recent work by Heide et al. [57] that proposeda statistical prior for reducing chromatic aberrations from single element opticalsystems.In recent years there have been a large number of computational imagers thatemploy a co-design of refractive or diffractive optics and computation. Examplesinclude cubic phase plates for extended depth of field imaging [95], DOEs produc-ing double-helix PSFs [116] for single-molecule microscopy beyond the diffrac-tion limit, phase microscopy [54], coded amplitude masks instead of lenses [5],and anti-symmetric gratings integrated with a sensor to produce an ultra-miniaturelens-less imager PicoCam [46, 48, 146]. Although these approaches have demon-strated the possibility to develop imaging systems with reduced optical complex-ity by shifting the burden to computational reconstruction, the resulting systemsexhibit little flexibility, such that refocusing or zooming are either not supportedat all, or require tedious re-calibration. Additional shortcomings of existing sys-tems are poor image quality and lack of colour imaging for diffractive designs,as well as high computational cost. In this work, we jointly exploit computation-ally optimized diffractive optics and computational image reconstruction (shown83in Fig. 6.1).Figure 6.1: Overview of the imaging system. We computationally designstacks of DOEs (b) that can encode arbitrary lenses for different geomet-ric configurations, such as relative shift or rotation. Examples includefocus-tunable lenses, tunable cubic phase-plates, and axicons shown in(a) with their corresponding phase transmittance functions. Given a tar-get phase, we design two or more phase plates using complex matrixfactorization (c). If more than two phase plates are used even complexoptics like zoom-systems can be designed. A novel computational ap-proach enables broadband imaging for our encoded lenses (d and e).When a scene is imaged with our diffractive encoded lens, points withdifferent spectral distributions result in significantly different PSFs. Ourreconstruction algorithm (e) jointly self-calibrates the spatially-varying,scene-dependent PSFs and recovers the latent image exploiting cross-channel statistics. Aberrations in the reconstruction are effectively re-moved (actual reconstruction shown here).We computationally design stacks of two or more DOEs in order to encodearbitrary lenses for different geometric configurations of DOEs. The relative posi-tioning of these DOEs (i.e. different rotations or translations) encodes different op-tical designs such as a variety of focal lengths or parameters of a cubic phase plate.Our approach is a generalization of the work [8, 9] which presented an analyticalphase pattern for two rotating DOEs to change focus. We achieve general encod-ings of lenses by formulating the design as a numerical optimization problem. This84optimization problem turns out to be a complex matrix factorization problem. Anefficient optimization method allows us to solve this problem and thus design twoor more DOEs such that, when arranged in different geometric configurations (e.g.rotations or translations), complex transmitances corresponding to different targetlenses can be generated. Other analytic designs have been described and analyzedin the prior work [7, 53, 78]. Note that all of these approaches propose analyticsolutions for selected target designs, while our method numerically optimizes forphase patterns given arbitrary target designs.Aberrations are deliberately tolerated in the images captured using our opti-cal designs and removed by a new self-calibrating computational reconstruction.This blind image reconstruction method jointly estimates the underlying scene-dependent, spatially-varying wavelength dependent point spread functions (PSFs)and the latent sharp image. Our method adaptively self-calibrates the PSFs from theimage itself as well as reconstructs sharp images without any physical calibrationof the imaging system for all possible geometric configurations. Exploiting sta-tistical cross-channel correlation as part of our joint reconstruction, major issuesof current state-of-the-art blind PSF estimation methods can be eliminated. Weshow that our method can restore high quality images exceeding the image qualityof previous fully calibrated (i.e. non-blind) approaches, while being 5 orders ofmagnitude faster than the state of the art.We experimentally show three applications of our flexible encoded diffractivelens — static colour imaging with a single diffractive interface, rotational refocus-ing and rotation-only zooming. In simulation (see Fig. 6.2), we further demonstrateadditional applications including shift-only focus changing, adjustable cubic phaseplates and time-multiplexed diffractive lenses. The proposed technique reveals itsversatility for achieving more compact, flexible, and powerful imaging systems.6.2 Methods6.2.1 Encoded diffractive lens designThe design of encoded diffractive lens is a complex tensor factorization problem.The target transmittance function T in polar coordinates (r,ω) can be described as85the multiplication of two transmittance functions T1 and T2:T(r,ω−θ) = T1(r,ω) ·T2(r,ω−θ), (6.1)where we have encoded optical parameters into the relative rotation angle θ ofthe second element. Our target is to design a varifocal lens whose transmittancefunction is:T(r,ω) = exp(−j2piλ· r22 f (θ)), (6.2)where λ is the design wavelength and f (θ) is the angle-dependent focal length. Itis possible to rewrite the above equation as a complex matrix factorization problemby remapping polar coordinates (r,ω) to linear indices addressing the columns oftwo complex matrices A and B. We solve the following optimization problem:Aopt,Bopt = argminA∈Cm×r|·|=1,B∈Cn×r|·|=112∥∥W◦T−W◦AB†∥∥2F , (6.3)where B† is the complex conjugate of B, ‖.‖F denotes the Frobenius norm, and◦ denotes the Hadamard product. A weighting matrix W is added to select onlythose matrices physically possible over each rotation angle. We refer the reader toHo [62] for a detailed introduction to weighted (non-negative) matrix factorizationmethods. For static DOEs, the A,B matrices are of rank-1, AB† can also be in-terpreted as the inner product of two vectors that correspond to the mask patternof each DOE. For dynamic patterns higher rank factorizations can be realized, asdescribed at the beginning of next section.For our discussion, we assume the static rank-1 setting, but we note that thesame methods can be applied for higher rank factorizations. The bi-convex ma-trix factorization problem can be efficiently solved as a sequence of convex sub-problems using an alternating approach. Note that we aim for converging to astationary point. In each sub-problem we diagonalize the weighting matrix andperform the outer vector product operation followed by vectorization. We can eas-ily derive the gradient and the Hessian results in a diagonal matrix whose inversionusing Newton’s method then becomes a point-wise division. The details will bediscussed later in Chapter 7 where we use a similar complex matrix factorization.86By symmetry, both sub-problems can be solved very efficiently.6.2.2 Reconstruction algorithmRaw images captured by our encoded lens suffer from inherent chromatic aber-rations introduced by the wavelength-dependency property of DOEs. Therefore,computationally reconstructing images without aberrations is a key part in ourimaging system. We devise a fast blind deconvolution process that jointly performsPSF self-calibration and deconvolution.Note that, in general, different spectral distributions of objects and illumina-tion cause spatially varying PSFs for a diffractive lens when measured on a coloursensor. This effect is shown in Fig. 6.1. Given three points with different spec-tra shown in Fig. 6.1 (d), the response of the lens and sensor lead to differentlycoloured PSFs (cyan, yellow, and magenta in Fig. 6.1 (d)) in different sub-regionsof the image. Note that material and illumination properties, on the other hand, arespatially low-frequency which result in a PSF that is spatially invariant in a localneighbourhood.To handle such spatially varying PSFs our reconstruction exploits cross-channelcorrelation as a statistical prior (see Chapter 3). In particular, the empirical distribu-tion of gradient differences between two spectral bands turns out to be heavy-tailed.Intuitively this means gradients of two colour channels are very similar in most ar-eas except for occasional pronounced changes in chroma and luma (e.g. object ormaterial boundaries). A nice property of the encoded diffractive lens is that, whentaking an RGB image, at least one channel can be focused well while the othertwo channels are blurry. This structure is directly exploited by the cross-channelprior. Note that our formulation model differs from that previously proposed [57]which only assumes chroma changes to be sparse. This assumption leads to severeinstabilities of the prior for low intensities.Using this modified cross-channel prior we jointly estimate the spatially vary-ing PSFs and the latent sharp image. Hence, the reconstruction process is self-calibrated and adaptive to the scene. Compared to other common methods whichrequire tedious and complicated physical calibration, our method avoids these stepsfor all possible geometric configurations and is therefore more convenient and ro-87bust. It turns out that adding the prior in the blind deconvolution converges to asignificantly better optimum in drastically less computational time than competingmethods.Self-calibrated blind deconvolution. As discussed, the image reconstruction methodjointly estimates the underlying spatially-varying PSFs and the latent sharp image.Having ensured that our encoded lens focuses well in at least one channel, we cansolve for the PSFs by exploiting cross-channel correlation between the channels.In particular, we solve the following optimization problem:xopt,vopt = argminx,v‖x∗v− j‖22+α3∑a=1‖ha ∗x‖22+β‖x‖22+γ∑a6=r‖ha ∗v−ha ∗ ir‖1+µ3∑a=1‖ha ∗v‖1,(6.4)where j is the captured raw image, v is the unknown latent image, and x is theunknown PSF that needs to be re-normalized at each iteration. The convolution ofthe latent image with the PSF is expressed as x ∗ v using the convolution operator∗. Here ha (a= 1,2,3) are first order spatial gradient filters for RGB channels, andir is a sharp image in reference channel r, and α,β ,γ and µ are weights for theregularization terms.In the first row we have a classical `2 data term, a gradient prior and a low-energy term on the PSF x in the `2 sense. The cross-channel prior on the latentimage v in the second row of Eq. (6.4) penalizes the gradients between blur chan-nels and the reference channel in the `1 sense. A sparsity penalty on the gradientsof images in all channels is further added. Refer to Appendix chapter A for adetailed discussion of all penalty terms from Eq. (6.4).The bi-convex minimization problem in Eq. (6.4) can be solved via coordinatedescent without any additional priors or optimization schedule tricks. We keep oneof the two variables x,v fixed while minimizing the objective with respect to theother in an alternating fashion. This approach leads to the two sub-problems, x-step and v-step, both of which are much simpler to solve than the joint objective.Note that due to a scale ambiguity between x and v, there should be a normalization88after x-step. The full alternating coordinate descent method is given in Alg. 2.Algorithm 2 Blind PSF Estimation Using Coordinate Descent1: x0 = δ , v0 = j2: repeat3: xk+1 = argminx∑a ‖Vkx− j‖22 + . x-stepα∑a ‖Hax‖22+β‖x‖224: xk+1 = Normalization(xk+1) . kernel re-normalization5: vk+1i = argminv∑a ‖Xkv− j‖22+ . v-stepγ∑a ‖Hav−Hair‖1+µ‖Hav‖16: k := k+17: until Optimality achievedMinimizing the quadratic from the x-step is equivalent to solving a system oflinear equations composed of structured matrices, all of which represent convolu-tions. This structure can be exploited by reformulating this linear equation systemin the frequency domain and inverted efficiently by point-wise division. The v-steprequires the solution of a deconvolution problem with known kernel xk. It involvesa quadratic data term, a sparse cross-channel correlation term and a sparse gradientterm. Due to the `1-norm penalty of the last two terms, solving this minimizationproblem does not reduce to a quadratic problem as for the x-step. We solve it us-ing a splitting approach (see Appendix chapter A). Compared with state-of-the-artmethods, ours converges to a significantly better optimum in drastically less com-putational time. Refer to the supplementary of [61] for quantitative evaluations.6.3 ResultsOur encoded diffractive imaging system consists of a numerically optimized op-tical design and self-calibrated reconstruction (aka. blind deconvolution) as ex-plained above. In the following, the computational design of the encoded opticsis discussed first. After that, the reconstruction approach takes the corrupted mea-surements as the input, and we show the experimental results in different imagingscenarios.896.3.1 Encoded diffractive lensConsider a rotational encoding of different phase function designs where each de-sign is implemented by a different rotation angle of two phase gratings. We de-scribe the complex target transmittance function as the product of two transmit-tance functions of the individual DOEs. For special transmittance functions suchas Fresnel lenses, approximations can be found analytically [8, 9], but these meth-ods do not generalize to arbitrary target transmission functions, and the accuracyof the approximation is limited. We instead reformulate the design problem as amatrix factorization problem, whereby a general target transmittance function isrepresented as a matrix that can encode all possible geometric alignments betweenthe DOEs since it contains an entry for every combination of pixels on the first andsecond DOE. A stack of two static DOEs can produce a Rank-1 approximation ofthis complex matrix, which can be found through complex matrix factorization.Note, however, that mutual rotation cannot “perfectly” encode optical elementswith non-rotational symmetric targets in a continuous angular range.The described principle can be generalized in several ways: stacking more thantwo diffractive masks results in a tensor factorization problem rather than a matrixfactorization problem. If one or more of the DOEs are replaced with a dynamicphase modulator that allows for temporally multiplexing different diffractive pat-terns, then the rank is increased to the number of different patterns displayed duringa single exposure. In this setting super-resolved features and higher bit-depth thanfor a single phase modulator might be achieved [59].Any target phase function can be encoded including Fresnel lenses, axicons [98],and cubic phase plates [16]. Figure 6.2 shows the phase profiles of both the tar-get lenses and the approximated results using our factorization methods. Varifocallenses can be achieved by encoding focal lengths into mutual rotation of two layersusing Rank-1 matrix factorization (Fig. 6.2 (a)) for two static elements. Similarly,focal lengths (Fig. 6.2 (b)) or cubic phase functions (Fig. 6.2 (c)) can be approxi-mated by relative translation of two layers using Rank-1 factorization. In the time-multiplexing case, higher rank factorization can be used. A focal length encodingexample of this kind is also shown in Fig. 6.2 (d) using Rank-4 factorization.We built four rotationally encoded lenses for experimental verification. All90Figure 6.2: Encoding examples of the optimization for different target phasefunctions. (a) Focal lengths are encoded into mutual rotation of twolayers using rank-1 complex matrix factorization. (b) Focal lengths areencoded into translation of two layers using Rank-1 complex matrix fac-torization. (c) Cubic phase functions are encoded into translation of twolayers using Rank-1 complex matrix factorization. (d) Focal lengths areencoded into time-multiplexing using Rank-4 complex matrix factoriza-tion. From left to right: four target phase functions, reconstructed phasefunction and their respective optimized phase patterns.the lenses are designed at wavelength λ = 550 nm on 0.5 mm thick Fused Silicasubstrates with an aperture diameter of 8mm. This fabrication setting follows thesame manner as in previous chapters. The focal length ranges are [100mm,200mm],[100mm,300mm], (−∞,−100mm]∪ [100mm,∞), and (−∞,−50mm]∪ [50mm,∞)respectively.6.3.2 Experimental setupOur encoded diffractive lens system is shown in Fig. 6.3 (a). The refractive lenson a Canon EOS Rebel T5 DSLR is replaced by a group of encoded diffractivelenses consisting of two DOEs stacked face-to-face. The encoded lens is placed atthe flange focal distance (44 mm) of the camera which contains a APS-C sensor(22.3 mm × 14.9 mm). The phase patterns are fabricated using photonlithography91Figure 6.3: Experimental set-up. (a) Our encoded imaging system consistsof a encoded lens, mechanical mounts and a camera body. The en-coded lens is mounted on a 3-degree-of-freedom adjustment cage sys-tem. A Canon EOS Rebel T5 camera body is connected to the mount-ing system by a lens tube to prevent stray light. (b) The central ar-eas of the two component DOEs in the design with focal length range(−∞,−100mm] ∪ [100mm,∞) are observed on a microscope NikonEclipse L200N with a 5× objective. The scale bar is 500 µm. (c) Set-upof the static broadband imaging system. Without movement, our systemis able to image the scene as a conventional lens. Our computationalapproach reconstructs high quality output images from the corruptedmeasurements. (d) Set-up of the rotational refocusing system. Withonly mutual rotation, our system is able to refocus on different objectswithin an extremely large depth of field. (e) Set-up of the rotation-onlyzooming system. Employing two encoded lenses, our system is able tozoom in or out with only rotation required to keep the objects sharp.on 4” Fused Silica wafers. Figure 6.3 (b) shows the 5× microscopic images of thecentral areas of both DOEs. When the two elements are rotated relative to eachother, the equivalent focal length varies as a function of the rotation angle.In order to achieve precise alignment (tilt and decenter) between the two DOEs,we first mount the two elements on two cage rotation mounts (CRM1/M, Thorlabs)respectively. The cage systems are then mounted in an XY translator with microm-eter drivers (ST1XY-S/M, Thorlabs). This three degree of freedom mounting sys-tem guarantees accurate alignment as well as smooth 360◦ rotation. The mountedlens is then connected to the camera body with a lens tube.92Three applications of the encoded diffractive lens are illustrated in Fig. 6.3 (c)-(e). First, the static configuration of the system is a traditional diffractive lensat a specific focal length (Fig. 6.3 (c)). Second, when the distance between theencoded lens and the sensor is fixed, the lens can focus on objects at differentdistances (Fig. 6.3 (d)) if one element is rotated with respect to the other. Third,a rotation-only zoom lens can be achieved by two encoded lenses separated by ashort distance. When the two groups are rotated at specific angles respectively, theimage can be zoomed in or out while remaining sharp. No longitudinal movementis involved leaving the overall length of the system unchanged. In all cases, ahigh quality image is reconstructed by the proposed joint blind PSF estimation anddeconvolution approach.6.3.3 Static colour imagingWhen the two DOEs are stationary with respect to each other, the equivalent phaseis a diffractive lens with a specific focal length. The focal length can range froma few millimetres to nearly infinity. In all cases, we have focused our lens for thegreen channel, that is, the central part of the green channel’s spectral response.We show several results for different scenes captured with different focal lengthsettings in Fig. 6.4 to demonstrate the flexibility of encoding focal lengths intorelative rotation.Figure 6.4 shows results of our encoded lens in static colour imaging. The in-door illumination is a fluorescent lamp. Although strong colour fringes occur atedges (Fig. 6.4 (a-b)), especially for white objects, these chromatic aberrations in-troduced by the diffractive lens have been completely removed in the reconstructedimage (see the corresponding insets for highlights). The algorithm is also robustto transparent materials with reflective surfaces (Fig. 6.4 (c)). In addition to imag-ing similarly to a conventional lens, our encoded diffractive lens can also performin macro mode to image very small objects (Fig. 6.4 (d-e)). For natural scenessuch as human faces (Fig. 6.4 (f)), the algorithm is also able to reconstruct a highfidelity image. Outdoor scenes under sunlight illumination are shown in Fig. 6.4 (g-h). Note that, although the aberrations are depth and wavelength dependent, boththe foreground and background objects can be recovered simultaneously to extend93depth of field (DOF).The resolution chart capture in Fig. 6.5 shows that resolution is increased whileeliminating chromatic aberrations. Having focused our lens again for the greenchannel, we can see that severe aberrations occur for the red and blue channel in (b-e). The proposed approach successfully eliminates these aberrations and recoverslost resolution. Especially the high resolution lines are reconstructed accurately forall patches shown in Fig. 6.5. Note that no deconvolution artefacts (e.g. ringing)are introduced. Our method is robust to noise and removes the strong noise in thepatches from Fig. 6.5 which usually causes deconvolution artefacts.6.3.4 Rotational refocusingSince the focal lengths change along with rotation angles, our encoded lens canperform refocusing with only mutual rotation while the distance between the lensand the sensor remains fixed. This is particularly useful when space is limited andlongitudinal movement of the lens is prohibited.Figure 6.6 (a) illustrates a refocusing setup with an extremely large depth offield. A very small screw is located 10mm in front of the lens and a poker cardon a magic cube is located 1500mm further away. Our system can easily focus onboth the far objects (Fig. 6.6 (b)) and near objects (Fig. 6.6 (c)) through rotation.When the poker card is focused, the screw is hardly seen due to strong out of focusblur. Similarly, when the screw is in focus the poker card is severely blurred. Anexample of a more complicated scene is shown in Fig. 6.6 (d). The lens is focusedon the far beam-splitter cage and the near Matryoshka doll is out of focus whilein Fig. 6.6 (f), the near doll is focused instead of the far cage. The sharp imagesare shown in Figs. 6.6 (e) and (g) respectively. Insets show that our algorithm canrestore details very well (e.g. the screw holes on the table). Note that the saturatedarea on the doll and the very dark area on the corner of the black cage indicate thatour algorithm is also robust to different exposure levels.6.3.5 Rotation-only zoomingConventional zoom lens requires two types of motion to change the image mag-nification while maintaining focus. First is the longitudinal moving of an inner94Figure 6.4: Results of static colour imaging with a single encodable lenses.For each individual result, both the blurry observation (top) and recon-struction (bottom) are shown for indoor and outdoor natural scenes. (a-b) Indoor scenes under fluorescent lamp illumination. (c) Transparentobjects with reflective surfaces. (d-e) Small objects imaged in macromode. (f) Human face. (g-h) Outdoor scenes under sunlight illumina-tion. Insets in each result highlight the details before and after compu-tational reconstruction. Chromatic aberrations are completely removedin all scenarios. The focal length range of the encoded lens for experi-mental results shown here is [100mm,300mm].95Figure 6.5: Resolution chart captures with a single encodable lenses. In (a)the blurred measurement (top) and the corresponding reconstruction(bottom) are shown along with a few extracted patches. The red andblue colour channels of these patches are shown in (b - e), the greenchannel was focused. The result (b) shows the red channel with theblurry capture (top) and the reconstruction (bottom). (c) shows the bluechannel for the same patch. (d-e) shows the same channels for a differ-ent patch. We can see that for both blurred channels our reconstructionshave a drastic increase in resolution. The focal length range of the en-coded lens for experimental results shown here is [100mm,200mm].group of lenses within the entire lens to change the focal length, and the secondis longitudinal moving of the whole lens to refocus. Therefore, zoom lenses areusually bulky, heavy and require precise mechanical components to maintain lensalignment while moving [171].Benefiting from the focal length encoding capability of our encoded lens, weare able to realize rotation-only zooming without longitudinal movement. Asshown in Fig.6.3 (e), our zoom lens consists of two encoded lenses separated by afixed distance. This is a two-component system in which the object distance, imagedistance and separation between the two lenses are all fixed. There exists a singlepair of focal lengths that correspond to a specific magnification [145]. Therefore,to achieve a change in magnification (i.e. zooming) the focal lengths f1 and f296Figure 6.6: Results of refocusing with relative rotation between two DOEs.(a) A screw is located 10 mm in front of the lens and a poker card ona magic cube is located 1,500 mm further away. (b) When the lens isfocused on the rear poker card, the front screw is hardly seen due tothe strong out of focus. (c) When the lens is focused on the screw, thepoker card is severely blurred. (d) Captured image when the lens is fo-cused on the rear beam-splitter cage. (e) Reconstructed image of (d).(f) Captured image when the lens is focused on the foreground doll. (g)Reconstructed image of (f). Insets of the results from (d) to (e) indi-cate that our algorithm is robust to both saturated and dark areas in thesame scene. The focal length range of the encoded lens for experimentalresults shown here is (−∞,−50mm]∪ [50mm,∞).must change subject to the following relationship:1f1=ms−md− s′msd1f2=d−ms+ s′ds′.(6.5)where the object distance s, image distance s′ and spacing d between the two groupsof DOEs are all fixed.Our two-group encoded lens system is able to achieve rotation-only zoomingby finding different sets of focal length parameters. All the adjustments involvedare implemented by mutual rotation of the two DOEs in a single group. Figure 6.797shows a set of reconstructed images captured at different magnifications. In eachimage auxiliary arrow marks indicate the magnification change compared to theprevious image. Note that due to the Chromium (Cr) aperture of our DOEs the im-age quality is somewhat degraded compared to the static broadband results becauseof the inter-reflection between the two groups of DOEs.Figure 6.7: Results of rotation-only zooming with two groups of encodedlenses. A sequence of images of the same lamp from (a) to (e) indicatethe magnification change. Arrows in the same colour have the same sizein all the images.6.4 ConclusionThe combination of computationally optimized diffractive optics and computa-tional reconstruction opens up new territory for high quality broadband computa-tional imaging. In this work we have achieved several major advances in imagingtechnology. First, we demonstrate that pairs of diffractive optical elements can becomputationally designed to encode a wide range of optical parameters into sim-ple mechanical rotations and translations. Thin lenses of varying focal lengths, orcubic phase plates and axicons with different parameters can all be achieved withcontinuous rotational or translational changes in alignment. More complex opticalsystems such as zoom lenses can be assembled with two or more pairs of DOEs.These are initial examples of this design approach, which should, however, ap-ply much more widely, thus paving the road for significantly more flexible futureimaging systems.Secondly, we show that, using computational imaging approaches, such en-coded diffractive optical elements can be used for full-spectrum and colour imag-ing. We believe this is the first time high-fidelity colour imaging has been demon-strated using only diffractive optical components.98Finally, we make significant improvements on two salient problems with manycomputational imaging methods: calibration effort and computational expense. Wedemonstrate a fully self-calibrated imaging system, in which image restoration isachieved on individual images without additional calibration information. Thismethod is also significantly faster than existing image restoration methods withoutsacrificing reconstruction quality.Encoded diffractive optics generalizes well to wavelengths outside the visiblespectrum. For imaging in a narrow spectral band at any frequency, the only re-quirement is that diffractive optical elements can be fabricated for that wavelength.For broadband imaging, we additionally require that the image sensor has multi-ple spectral channels, one of which corresponds to the design wavelength of thediffractive optical elements. We believe that this property makes encoded diffrac-tive optics particularly interesting in the UV or THz range, where refractive opticsis either expensive or not available at all.99Chapter 7Mix-and-match HolographyUntil now, the previous chapters present a number of diffractive lens designs onthe capture end. As highlighted at the very beginning of this thesis dissertation,information delivery is also a key component. Computational caustics and lightsteering displays offer a wide range of interesting applications, ranging from artworks and architectural installations to energy efficient high dynamic range (HDR)projection.In this chapter we expand on the concept of computational caustics by encodingseveral target images into pairs of front and rear phase-distorting surfaces. Differ-ent target holograms can be decoded by mixing and matching different front andrear surfaces under specific geometric alignments. Our approach, which we callmix-and-match holography, is made possible by moving from a refractive causticimage formation process to a diffractive, holographic one. This provides the extrabandwidth that is required to multiplex several images into pairing surfaces.Particularly, we derive an image formation model for the setting of holographicprojection displays, as well as a multiplexing method based on a combination ofphase retrieval methods and complex matrix factorization. We demonstrate severalapplication scenarios in both simulation and physical prototypes.1007.1 IntroductionGoal-based caustics and computational light steering have for a long time beena promising research topic in computer graphics, with proposed applications inart and architecture (e.g. [36, 77, 113]) as well as energy efficient high dynamicrange (HDR) projection [28, 64]. The general problem can be stated as designingan optical system that arranges micro-structures, material distributions, or otherphysical parameters of an optical element to produce a target light distributionunder a given illumination condition.Traditionally in graphics, this problem has been tackled using a refractive set-ting, i.e. by optimizing a free-form lens with one flat and one curved or patternedsurface, that generates the desired light distribution under collimated white light.The optics community has also extensively explored the diffractive, or holographicsetting, which offers several advantages. First, refractive designs for computationalcaustics can be quite bulky, while diffractive designs are flat at a macroscopic scale,and can therefore be very thin and lightweight. Second, whereas refractive designsmust aim for relatively smooth surfaces in order to avoid diffraction artifacts, holo-graphic approaches directly incorporate diffraction into the image formation pro-cess, and can therefore make use of smaller features, which drastically increasesthe bandwidth (i.e. the number of degrees of freedom) for the optical design. Fi-nally, we note that, while diffractive optical elements (DOEs) are challenging toprototype, they can be mass-produced at high speed and low cost using methodssuch as nano imprinting [1, 21].In the research of holographic techniques, angular or wavelength multiplexinghave been developed to further increase the bandwidth of holograms so as to en-code a large amount of desired data. In our work we propose mix-and-match holog-raphy, which is a new spatial multiplexing scheme for holography, by which differ-ent target holograms are encoded onto pairs of front and rear DOEs (see Fig. 7.1).These DOEs are designed in such a way that different pairings or geometric align-ments represent different target holograms. By mixing and matching different frontand rear DOEs and alignments, we can decode the different holograms, and thesecan in turn be used to produce different target projected images when illuminatedunder a specific lighting condition. In addition to adding an interesting dynamic101component to existing applications of caustic displays, this mix-and-match designenables new applications in security, similar to the work by Naor and Shamir [108],since the two phase DOEs can be thought of as an encrypted message and a key:both are required to decode the target information, but can be stored separately,unlike, for example, volume holograms.In this work, the technical contributions are as follows:• We re-formulate the goal-driven caustic problem as a holographic recon-struction problem, and derive detailed diffractive image formation modelsfor this scenario.• We exploit the larger design space of diffractive optics by firstly using amodified computer-generated hologram algorithm to retrieve a target phasedistribution, and more importantly, using complex matrix factorization tomultiplex holograms of several target images onto a pair of diffractive sur-faces via mix-and-match holography.• We synthetically and experimentally validate four application scenarios sub-ject to different mix-and-match schemes, including encoding high-contrastholographic image through offset layers, representing multiple target imagesthrough combinatorial pairings, reconstructing animated images through lin-ear translation, and encoding multiview images through a barrier-based dis-play configuration.7.2 Computational designIn this section we discuss in detail the two aspects of the design process — thecomputer-generated holography and the multiplexing step. First, however, weagain briefly review the image formation using Fresnel diffractive propagation.7.2.1 Diffractive image formationRecall the diffractive propagation analysis in Chapter 2. Using the Fresnel approx-imation, the light propagation from the hologram plane to the image planeis then102Figure 7.1: Processing pipeline. Given a sequence of target inputs and theuser defined reconstruction requirements, we first obtain a sequence ofhigh fidelity phase holograms by applying a modified GS-based phaseretrieval algorithm. Following that, we multiplex multiple of thesephase holograms into a latent complex transmittance function, based ona designated alignment scheme and fabrication constraints. Relying oncomplex matrix factorization, two factorized matrices (vectors) can thenbe found, which represent the complex transmittance functions of twodiffractive optical elements. Finally, by pairing two phase DOEs un-der different geometric alignments, we can de-multiplex multiple holo-grams. The corresponding projected images are formed by illuminatingthe setup under a specific lighting condition. (Albert Einstein Photo-graph by Orren Jack Turner)given as:U˜(x, y) =exp(jkz1)jλ z1∫ ∫∑U˜(x1, y1)exp{jk2z1[(x− x1)2+(y− y1)2]}dx1dy1,(7.1)where (x,y) and (x1,y1) refer to positions on the image plane and the hologramplane, respectively. j is the imaginary unit, k the wave number, and z1 stands fordistance from the hologram plane to the image plane. According to Eq. 7.1, onecan calculate the integral over the whole DOE transmittance aperture ∑ to computethe complex amplitude distribution U˜(x,y), including both amplitude and phase in-formation. Please note that Eq. 7.1 represents essentially a Fourier transformation(Fourier transform (FT)).Notice that U˜(x1,y1) is the wave exiting the DOE plane, i.e. after diffractionin the hologram has already occurred. It is related to the incident wave U˜0 at the103hologram via the complex transmittance function t(x1,y1):U˜(x1,y1) = U˜0(x1,y1) · t(x1,y1) (7.2)witht(x1,y1) = A(x1,y1) · exp(jΦ(x1,y1)), (7.3)where A represents the amplitude and Φ represents the phase of each point on thehologram plane.In our case, we design a phase hologram, meaning that the amplitude is con-stant across the plane, and thus the transmittance function is determined exclusivelyby the phase distribution Φ(x1,y1), which is physically represented by the heightprofile h(x1,y1) of a DOE:Φ(x1,y1) = (n−1)h(x1,y1)λ 2pi. (7.4)Putting all equations together, we can now design a height profile h on a sub-strate of refractive index n, that results in a given complex target amplitude U˜(x,y)on the image plane. In the following, both the transmittance function and the heightprofile are pixelated into discrete squares of constant phase retardation and height,respectively. We call these squares cells in order to avoid confusion with pixels inthe target images.7.2.2 Hologram optimizationOur hologram optimization approach (summarized in Fig. 7.2) follows the generalGerchberg-Saxton (GS) framework and its extended versions [35, 44, 139]: westart with a target image in the image plane and convert it to a wave representationU˜ (0)(x1,y1), where the phase at each point is random, and the amplitude repre-sents the image intensity at each point. This complex waveform is then propagatedbackwards from the image plane to the hologram plane using Eq. 7.1 to obtainU˜ (0)(x,y). Since we are designing a phase-only hologram, the amplitude is set to 1across the hologram plane, yielding an updated waveform U˜ (1)(x,y), which is thenforward-propagated to the image plane (again using Eq. 7.1), to give U˜ (1)(x1,y1).In the image plane, an amplitude constraint is applied by enforcing that the ampli-104Figure 7.2: Flowchart of our modified double-sampling GS algorithm. Byrelying on a virtual intermediate sampling plane (see Appendix for de-tails), the sampling rates on the image plane and the hologram plane canbe efficiently decoupled and tuned subject to design intentions. This re-sults in a two-step FFT light propagation process. The input is a gray-scale target image (intensity distribution), while the output is a phasehologram.tude matches the intensity distribution of the target image, before starting the nextiteration.As mentioned before, both the forward and the backward transport are essen-tially Fourier transformations. However, a single-step discrete Fourier transform(DFT) as used in a standard GS algorithm has two constraints. First, the sam-pling rate on the image plane is determined by the wavelength and the propagationdistance. Second, the sampling pattern on the image plane is restricted by thesampling pattern on the hologram plane, and vice versa. We would like to designholograms that expand the projection, such that a small DOE can produce a largeprojected image. This requires different sampling rates and regions of interest onthe hologram plane versus the image plane.Modified GS-based phase retrieval We apply a two-step FFT process to elimi-nate of the sampling limitation imposed by transport with a single discrete Fouriertransform (DFT), similar to the work by Okada et al. [111]. As shown in Fig. 7.3,this procedure uses a combination of two FFTs with a virtual intermediate plane, to105Figure 7.3: Configuration of diffractive wave propagation under a double-sampling scheme, where the red dashed arrows indicate the first (back-ward) propagation, and the blue solid arrows indicate the second (for-ward) propagation. By relying on a virtual intermediate sampling plane,the sampling rate ∆x on image plane and the sampling rate ∆x0 on holo-gram plane can be decoupled and efficiently tuned subject to designintentions.replace the original single forward FFT process [111, 123]. Another modificationis to use a spherical incident wave (i.e. a point light) rather than a regular planarwave, similar to the work by Nagahama et al. [106]. We position the center ofthis spherical wave at the virtual intermediate plane. This strategy is effective toenlarge the projection angle, and offers the side benefit of a relatively uniform lightdistribution that reduces speckle noise.Specifically, during each step, applying single-step DFT gives two samplingrates:∆xv =λ rN∆x0,∆x =λ (r+ z)N∆xv, (7.5)where ∆xv, ∆x0, and ∆x are the sampling rates on virtual intermediate plane, holo-gram plane, and image plane, respectively. r is the distance from sphere wave lightsource to hologram plane, while z is the propagated distance from hologram planeto image plane. N is the original sampling number on hologram plane in DFTprocess.106Figure 7.4: Illustration of the matrix structure for our multiplexing scheme.Left: matrix formation for translational pairing. Center: a sequence oftarget holograms for multiplexing. Right: matrix formation for offsetpairing. Specially, the matrix can be represented as the outer productof two vectors a,b, representing two DOEs, where the physically over-lapped cells are selected via different geometric alignments (as shownin different color cell fragments). fw here indicates an element-wisemapping from one of the target holograms tk to the target matrix T.Thus, we can derive the relationship between ∆x and ∆x0, as follows:∆x =(z+ r)∆x0r. (7.6)Now we successfully derive the sampling configuration without involving λ and N,which gives us flexibility to tune different sampling rates and regions of interest onhologram plane versus image plane. This is particularly beneficial when designingholograms that expand the projection, such that a small DOE can recreate a largeprojected image. The calculation time for the modified double-sampling diffractionis O(N2 log2 N), where a hologram has the size of N by N pixels.7.2.3 Mix-and-match multiplexingAfter the CGH step, we have obtained a target transmittance function tk for eachof the k holograms. We would now like to multiplex the k transmittance functionsonto two DOEs under different alignments. We assume that the two DOEs are indirect contact, so that the alignment has two degrees of freedom for translation, andone for rotation, although in this work we make no use of the rotational degree offreedom. We denote the complex transmittance of cell i in the first DOE as ai, whilesimilarly, b j denotes the complex transmittance of the cell j in the second DOE.107As a result, the synthesized complex transmittance of each paired cell fragment ti, jcan be obtained as:ti, j = wi, j(ai ·b j), (7.7)where wi, j here denotes the physical overlap of cell i and cell j (i∈ [1,n], j ∈ [1,m])in a given geometric configuration.We note that each individual tk can be expressed by a set of ti, j subject to aspecific geometric alignment of the two phase DOEs. As shown in Fig. 7.4, themultiplexing procedure can then be expressed as a mapping via a function fw,defined as:ti, j← fw(tk). (7.8)In vector-matrix form, all the input target holograms can be stacked into a singlematrix T:T = W◦ (AB†), (7.9)where B† is the complex conjugate of B, and ◦ denotes the Hadamard product. Weinclude a sparse weighted matrix W to select only those cells that are physicallyoverlapped.Matrix factorization. With this notation, it is easy to see that the computationaldesign of multiplexing is in fact a complex matrix factorization problem of findingtwo phase distributions that best approximate the target hologram in a least-squaressense. Specifically, we rewrite Eq. 7.9, by addressing the columns of two complexmatrices A and B [59, 62]. Then, the goal becomes to solve the following con-strained least-squares optimization problem:Aopt,Bopt = argminA∈Cm×r|·|=1,B∈Cn×r|·|=112∥∥W◦T−W◦AB†∥∥22 . (7.10)Note that r here is the rank of the approximation, with r = 1 if both DOEs havestatic (time-invariant) transmittance functions.To solve this problem, the following well-known multiplicative update scheme108can be used [62], expressed as:Ak+1← Ak ◦ (W◦T)Bk(W◦ (AkB†k))Bk,Bk+1← Bk ◦ A†k(W◦T)A†k(W◦ (AkB†k)).(7.11)Notice that the division here is the Hadamard (element-wise) matrix division. Forthese update rules to be well-defined, we need to ensure that the denominators arenon-zero by adding a small value ε in the algorithm. The above update rule isstraightforward and effective, however, the convergence of this algorithm is slow,which affects computational efficiency.Alternatively, one can use a gradient descent update scheme [23, 91] to solveEq. 7.10 efficiently. Assuming function g(Ak,Bk) represents the non-negative leastsquares term in Eq. 7.10, the update rules now are defined as follows:Ak+1← Ak−αk∇Ag(Ak,Bk),Bk+1← Bk−αk∇Bg(Ak,Bk),(7.12)where αk is a step size that is adapted in every iteration. Specifically, αk is updatedas αk+1← αk/β or αk+1← αk ·β , as described in detail in [91].During the complex matrix factorization, it is necessary to enforce the physicalconstraints of the optical system. This is analogous to enforcing non-negativityin traditional non-negative matrix factorization (NMF) problems. In our case, thephysical constraint is that we desire a phase-only modulation in both DOEs, so thatthe amplitude of the complex transmittance must equal 1 in each cell. We enforcethis with the following update rules for the individual components of A and B,which we call conditional projection:Ak+1i ←Ak+1i|Ak+1i ||Ak+1i |> 1,Ak+1i else,(7.13)where | · | denotes the amplitude of a complex value.109Note that in theory it is reasonable to force the extracted phase component ofevery pixel to be unified for the amplitude all the time, while in practice, it is notnecessary to normalize it at every step since we only need the final optimum tofit the physical model. Thus, our conditional projection maintains as much of thedata as possible while eliminating ill-conditioned cases, leading to not just fasterbut also better convergence. Once we approach convergence, the complex phaseexpression can then be forced to have unit amplitude.As mentioned, if we assume static DOEs, the matrices A,B are of rank 1.Specifically, A is a n×1 matrix (aka. a vector), B is a m×1 matrix, W is a n×msparse weight matrix, and T is a n×m sparse matrix containing the target complextransmittance of each cell. We can rewrite the model in Eq. 7.10 as:bopt,aopt = argmina∈Cm×1|·|=1,b∈Cn×1|·|=112∥∥T−W◦ab†∥∥2W . (7.14)For the Rank-1 weighted optimization, we use Newton’s updates to enforcea faster convergence. Algorithm 3 shows the specific algorithm we use to solvethis Rank-1 problem. It consists of an unconstrained, bi-convex rank 1 alternatingleast-squared weighted matrix factorization, followed by a constrained projectionoperation after solving each subproblem, as in Eq. 7.13. In each subproblem, wediagonalize the weighting matrix and compute the outer product, followed by vec-torization.Algorithm 3 Rank 1 Alternating Least Squares Weighted Complex Matrix Factor-ization1: k = 0, a0opt = ainit , b0opt = binit2: repeat3: bk+1opt := argminb12∥∥T−ab†∥∥2W . b-step4: bk+1opt := pro j(|bk+1opt |>1)(bk+1opt ) . conditional projection5: ak+1opt := argmina12∥∥T−ab†∥∥2W . a-step6: ak+1opt := pro j(|ak+1opt |>1)(ak+1opt ) . conditional projection7: k := k+18: until Optimality achieved110Newton’s update. Here we derive the b-step in Alg. 3, and the derivations fora-step follow from symmetry. Specifically, we first reformulate the optimizationproblem to derive f as:bopt = argminb12∥∥T−ab†∥∥2W= argminb12‖diag(W) t−diag(W)Oab‖22= argminb12f (b),(7.15)andf (b) =t† diag(W)2 t−2t† diag(W)2 Oab+b†O†a diag(W)Oab,(7.16)where diag(·) puts the matrix from the subscript on the diagonal and O(·) corre-sponds to the outer vector product operation with the vector in the subscript and theright hand side, followed by vectorization. Then we can easily derive the gradientas:∇b f = O†a diag(W)Oab−O†a diag(W)2 t. (7.17)The operator O†a is the same as the outer vector product operation plus a sub-sequent summation over the rows of the resulting matrix. We finally get a diagonalmatrix Hessian as:∂ 2 f∂b2= O†a diag(W)Oa. (7.18)Since Hessian is diagonal, the matrix inversion in Newton’s method becomesa point-wise division, expressed as in Alg. 4. Thus, we have cast the non-convexproblem as a sequence of convex optimization problems. The b-step and a-stepin Alg. 3 can be solved using gradient-based or Hessian-based methods, such asfast Newton’s updates. Our method shares similarities with Encoded DiffractiveOptics [61], although we are now dealing with holograms that exhibit more ran-domness than that of a symmetric imaging lens or a relative smooth phase profile.In addition, we apply a different projection rule (refer to Eq. 7.13) rather thannaively projecting all current optima back to their phase profile. This additional111filtering projection leads to a faster and better convergence than state-of-the-arts.Algorithm 4 Newton Update to solve the b-step1: k = 0,b0opt = ainit ,b0opt = binit2: repeat3: bk+1opt := bkopt− ∇b f∂2 f∂b2. Pointwise division4: bk+1opt := pro j(|bk+1opt |>1)(bk+1opt )5: . Conditional projection6: k := k+17: until Optimality achieved7.3 Applications and resultsThis section describes four specific applications of mix-and-match holography —high-contrast holographic projection via offset pairing, versatile image reconstruc-tions via multiple pairings, animated image reconstruction via translational pair-ings, and barrier-based multi-view projection display.For all our applications we built custom phase DOEs using photolithographyand reactive ion etching (RIE). For our current proof-of-concept stage, four itera-tions of lithography and RIE were repeated to produce 16 discrete height levels ineach DOE [118]. The principal design wavelength is chosen as λ = 550 nm, andthe corresponding total etching depth is 1195 nm. The feature size on the DOE isset to 8 µm. The higher resolution would provide additional degrees of freedomfor multiplexing a larger number of holograms. The right most close-up of Fig. 7.1shows a 20× microscopic image of the central areas of our fabricated DOEs.The pairing of DOEs requires element-wise alignment (i.e. tilt and de-center).Accordingly, alignment marks are added to the individual DOEs. During experi-ments, the DOEs are mounted on a 3 degree-of-freedom alignment system for ac-curate control of position and rotation (see left-most close-up in Fig. 7.6). For theillumination we use a Thorlabs 532 nm laser kit, a 5× beam expander, and a cus-tom lens system to realize different lighting conditions (collimated or spot light).Except for the DOEs, all the auxiliary components are off-the-shelf products.1127.3.1 Application scenariosIn Sec. 7.2, we discussed the general principle of hologram multiplexing in a mix-and match sense. In the following we describe 4 specific scenarios that are enabledby different choices of the layer matrices A and B, as well as the weight matrix W.Offset layered high-contrast projection. The first application is a “super-resolution”example with just one target image. The feature size in the hologram has a directimpact on contrast and image quality of the target image, but depending on thefabrication process it can be difficult or expensive to further reduce the featuresize. This example shows that the factorization from Section 7.2.3 can be used toapproximate a high resolution (4µm feature size) hologram with two lower reso-lution (8µm feature size) DOEs as long as the global alignment can be controlledwith 4µm accuracy.Specifically, in this application, the alignment of the two DOEs is designed foran offset of half the feature size (Fig. 7.4). Accordingly, the two layers A,B areexpressed as column vectors a ∈ Cn×1|·|=1 and b ∈ Cm×1|·|=1, respectively. We can effi-ciently solve this rank 1 optimization problem. Simulation results for this scenarioare shown in Fig. 7.5. From the PSNRs and the close-ups (red and blue rectangles),we see that the offsetting configuration reconstructs an image with higher contrastand less speckle noise than that of a single layer. We also note that the maximumdiffraction angle is dependent on the feature size. Decreasing the feature size al-lows for enlarged diffraction angles, and with our factorization we can achieve thiseffect without the need to use more precise fabrication techniques. Particularly,we note that this design feature of incorporating super-resolution effects into holo-grams can not be achieved by state-of-the-art CGH methods (see Chapter 2).A real result of projecting an image through an offset pairing of DOEs is shownin Fig. 7.6. Illuminated by the designated green incident light, we obtain a high-contrast image at the designated distance (center green close-up), whose quality isvisually matched to the synthetic result. Note that the observed speckle is inherentin the use of highly coherent laser illumination, and can be reduced with a numberof strategies that are orthogonal to our approach [28]. As one would expect, usinganother wavelength of light that differs from the design wavelength by >100 nm113Figure 7.5: Synthetic results of the sub-cell offset design. The images arereconstructed from a single DOE with a cell size of 8 µm (left col-umn), two DOEs with a cell size of 16 µm that are offset by half a cell(middle), and a single DOE with 16 µm cells (right). The PSNR areevaluated. Notice that the images are resampled and integrated over thephysical size of 200 µm per pixel subject to the angular resolution limitof the human eye at the observation distance of 50 cm. The same settingapplies to the followings. (1951 USAF Resolution Chart by MIL-STD-150A Standard, Albert Einstein Photograph by Philippe Halsman)results in a blurred image (right red close-up).Multiple discrete pairings. A first example of creating multiple images via dif-ferent pairings of DOEs, is to have each DOE composed of several discrete blockssuch that each combination of blocks from the two DOEs represents a hologramfor a different target image, as shown in Fig. 7.7. For instance, if we have 9 targetimages, they can be assigned to a matrix T containing 3× 3 blocks. We createthree blocks in the column direction to form a matrix A and three blocks in the rowdirection to form a matrix B. Specially, we put all the elements on the diagonal ofeach individual block. Now we have sparse matrices A ∈ C3m×m|·|=1 , and B ∈ Cm×3m|·|=1 .114Figure 7.6: Synthetic and experimental results of projecting an imagethrough an offset pair of DOEs, with the diagram of set-up indicated(left). We present synthetic results (center left), experimental results il-luminated with green laser light (center right), and experimental resultsilluminated with red laser light (right-most). We present two exampleshere. The size of DOEs are 1.4 cm by 1.4 cm individually, and the pro-jection distance is set at 1m. The final projected image size is around 12cm by 12 cm. (Albert Einstein Photograph by Orren Jack Turner, LenaSo¨derberg Photograph by Dwight Hooker.)Either Eq. 7.11 or Eq. 7.12 in Section 7.2.3 can be used to solve this problem.In this scenario we only need to fabricate 2n DOEs in two groups to representn2 images. This application mode also has applications in privacy or security [108],that differ from prior holographic encryption examples that are discussed in Chap-ter 2. We can interpret one of the DOEs as a set of encoded messages, while theother is a set of keys. Different keys can unlock different messages. For instanceto identify multiple banking accounts under one username, the user holds only onepiece (either digitally or in real plate) while the bank holds multiple paired pieces,or vice versa. Both synthetic and experimental results of this scenario are shownin Figs. 7.7 for images and 7.8 for message encoding.For Fig. 7.8, we show specifically that, if image quality is not of paramountimportance, the feature size of the holograms can be reduced fairly significantly,115Figure 7.7: Synthetic (left) and experimental (right) results for all 9 possiblecombinations of 3 individual patterns on the front DOE and 3 individualpatterns on the back DOE. Each individual image is reconstructed byaligning one pattern from group A and one pattern from group B. Weuse the same incident light source and the same projection distance as inFig. 7.6. The final reconstructed image size is around 3.7 cm by 3.7 cm.(Multiple source images modified from resource in the public domain,specifically, Felix Cat credit at DreamWorks Animation, Hello Kittycredit at SANRIO CO., Speedy Gonzales credit at Warner Brothers.)which makes the mix-and-match patterns more robust to misalignment. The featuresize in this example is 50 µm, which allows the phase DOEs to be mounted into in-dividual frames (whose procedure is a one-time work and can be done practically)which can be snapped together easily with sufficient precision and accuracy to de-code the multiplexed image. A hands-on mix-and-match demo is presented in thesupplementary video. We emphasize that this is a stress test for the image quality,which suffers from the lower diffraction efficiency in DOEs with larger features. Itis certainly possible to building mechanical frames for alignment accuracy muchbetter than 50 µm.116Figure 7.8: Experimental setup (left) and decoded (right) results for all 4 pos-sible combinations. For this demonstration, we have released the featuresize to an extreme value of 50 µm and aligned the DOEs onto mechan-ical frames, allowing for fast hand-on manipulation. The labels F andR indicate the front and rear DOEs, respectively. By snapping togethertwo frames from each group, the messages is decoded. We use the sameincident light source and the same projection distance as in Fig. 7.7. Thefinal reconstructed image size is around 1.6 cm by 1.6 cm.Animated projection under translation. By extending the multiplexing from a dis-crete set of static pairings to a translational alignment, one can generate animations.The concept is to shift a top layer DOE along one direction to overlap the differentregions of a bottom layer DOE (see the left close-up of Fig. 7.9). At each desig-nated alignment, the overlapped region of two layers generates the hologram fora specific frame of the animation. In our example, the bottom layer is only twicethe width of the top layer, but we multiplex 26 frames into this configuration. Inthis scenario, the two layers A,B are expressed as column vectors a ∈ Cn×1|·|=1 andb ∈ C2n×1|·|=1 , respectively. The experimental results of projecting 26 key frames ofanimated images through translation pairing of DOES are shown in Fig. 7.9.Barrier-based multiview projection. Theoretically any subregion of a hologramis able to reconstruct partial frequencies of the full target image. This unique prop-117Figure 7.9: The experimental results of projecting an animation throughtranslation DOE pairings. Left: diagram of experimental set up. Topright: synthetic results of selected frames. Bottom right: experimen-tal results of selected frames. The projection distance is set to 1m as inFig. 7.6, while here we use a planar incident illumination to simplify thesetup. The labels in the bottom-left of each sub-figure indicate the framenumber. Refer to the supplementary document of [119] for the completeframe sequence. (The heavy guy character credit at Vavle [103])erty makes it possible to synthesize a hologram that contains a partial design formultiple perspectives. By incorporating coded illumination via a binary barrier ora lenticular lens sheet, as in conventional auto-stereoscopic displays, one can mul-tiplex holograms of multiview images into a single pair of DOEs. A possible con-figuration is illustrated on the left of Fig. 7.10. In this setting, the illumination fromone direction only reaches certain regions of cells to reconstruct the correspondingperspective image. Note that this application can be combined with sub-cell align-ment to recover some of the bandwidth lost by multiplexing multiple views intothe same DOE pair.The synthetic and experimental results of multiview projection are shown onthe right of Fig. 7.10. One can see the change in occlusion between the fish (front)and the coral (back). Readers may also observe a residual image of the mask in thecenter of the measured image, which does not appear in the synthetic results. Thisis mainly due to a loss of diffraction efficiency caused by fabrication errors as wellas variations in the illumination wavelength.118Figure 7.10: Top: system diagram. Center: synthetic results. Bottom: exper-imental results of projecting multiview images via a combination of aparallax barrier and directional incident illumination. The sizes of theDOEs are 1.4 cm by 3.2 cm individually, where we multiplex 5 per-spective holograms. In our experiments, d = 1.2cm and ∆ = 2.24mm(Fish scene light field images credit at [164]).7.4 AnalysisThe presented prototype architecture and the resulting images that can be achieveddemonstrate the feasibility of the concept. In the following we discuss some of theissues that our initial prototypes suffer from.Image artifacts. The observed artifacts include image blur, noisy energy distribu-tion, and image contrast loss, which can be traced back to the limitations of currentmanufacturing and assembling process. We quantitatively evaluate the perceptualquality of a reconstructed image using the offset pairing design. The variables119added include: random Gaussian white noise to approximate fabrication error, andrelative in-plane shift and rotation to approximate misalignment of two physicallayers. Thus, based on the ground-truth image reconstructed with well-definedconditions (center left subfigure in Fig. 7.6), we present structural similarity in-dex (SSIM) assessment [162] in Tab. 7.1. A full version of comparison visualiza-tion is given in Fig. 7.11.Table 7.1: SSIM comparisons of synthetic images reconstructed from the off-set pairing subject to different experimental conditions. The columns in-dicate the results under the conditions— 1 weak noise, 2 strong noise,3 small rotation, 4 large rotation, 5 one directional shift, 6 dual direc-tional shift, 7 weak noise + small rotation, 8 weak noise + large rotation,9 strong noise + small rotation, 10 strong noise + large rotation, respec-tively.Condition 1 2 3 4 5SSIM 0.994 0.975 0.712 0.525 0.476Condition 6 7 8 9 10SSIM 0.420 0.708 0.523 0.697 0.520We find that relative shift and rotation account for the main factors that degradethe image quality. A direct consequence is the existence of ghost (crosstalk) imagesdue to the overlap of diffraction orders. With careful alignment this artifact canbe drastically reduced. Where precise pixel level alignment is not practical, aneffective solution is to pad extra zeros outside the target image in the phase mapretrieval (see Sec. 7.2). This trick directly enlarges the physical distance of twoneighboring diffraction orders such that the overlap issue can be mitigated withinthe field. In our experiments, misalignment of less than one cell (< 8 µm) will notresult in significant image degradation. Although tedious, this alignment accuracyis quite feasible in modern manufacturing and engineering methods.Hardware limitations. The accuracy of DOEs fabricated by our photolithographyprocess is highly dependent on the alignment accuracy between the used lithogra-phy masks and the stability of etching speed. Empirically, we have designed severalpairs of alignment fiducials (from 20 µm to 1 µm ) to assist the alignment, and haveused a gradual etching method to enforce a stable etching speed. In practice some120Figure 7.11: Comparisons of synthetic images reconstructed from the off-set pairing subject to different implementation conditions— (a) per-fect condition, (b) weak noise, (c) strong noise, (d) small rotation, (e)large rotation, (f) strong noise + small rotation one (g) directional shift,(h) dual directional shift, (i) strong noise + large rotation, respectively.Notice that all images have been resampled and integrated over 200µmper pixel subject to human eye’s angular resolution at a viewing dis-tance of 50cm. (Lena image source from Wikipedia)of DOEs are well-fabricated (Fig. 7.6 and Fig. 7.7) while others (Fig. 7.10) maysuffer from over-etching errors. Like other DOEs with discrete levels of surface re-121lief, our fabrication process cannot guarantee perfect single order diffraction. Thishigher order diffraction artifact is still difficult to eliminate.Our prototype is also limited by the micro-meter level accuracy of the opticalmounts we built. With additional engineering efforts (as demonstrated for the hand-on pairing scenario) or mass-production, most of these alignment issues could beaddressed. We also note that, for mass production, very inexpensive and fast meth-ods exist for reproducing the same high-resolution phase DOE over and over again.For example, nano-imprinting [1] can stamp sub-wavelength features directly intopolymer film at the scale of rolls of material. Therefore, the aforementioned issuesprimarily affect prototyping and research, but not potential commercial applica-tions.Computational performance. Although convergence to a global optimum cannotbe guaranteed for either the modified GS algorithm or the multiplexing method,both problems are in practice well behaved. Empirically, the modified GS algo-rithm converges in around 60 iterations for a zero-padded hologram with a reso-lution of 7,000 × 7,000 pixels. This takes about 15 minutes in our Matlab imple-mentation on a PC with Intel Xeon i7 2.70GHz CPU. The multiplexing step takesanother 2 minutes on the same system.Within our work, we consider two types of efficient factorization algorithms.We analyze their convergence and the run time performance. Some results areshown in Fig. 7.12. We make two observations — first, Newton’s update generallyleads to a faster convergence than the well-known multiplicative update; second,enforcing a conditional projection (recall Sec. 7.2) further outperforms conven-tional naive projection for both methods. Please refer to the supplement for de-tailed quantitative evaluations, including the factorization of a Fresnel lens profilethat is derived from analytic calculation.Notice that we solve a constrained matrix factorization problem over the do-main of complex numbers. The physical constraints (i.e. enforcing projection)are somewhat akin to a traditional non-negative matrix factorization (NMF). Asin NMF, the resulting non-convexity mirrors more or less the situation with NMFalgorithms as well. That said, there may be other formulations, which may yieldeven better results without projection. We leave this open but interesting question122Figure 7.12: Convergence for two factorization methods (rank 1 factorizationon a 1,751 × 1,751 pixel target image) with and without projection(Eq. 7.13). The two images tested are same as in Fig. 7.5.to future work.Spectral performance. While designing at 550 nm, we use easily available 532nm and 650 nm laser light sources in Fig. 7.6. This wavelength deviation may re-sult in diffraction efficiency loss. In Fig. 7.6, we have presented a photographedresult of directly illuminating a green light designed DOE pair with red illumina-tion. Noticeably, the image is degraded because of the wavelength dependency ofdiffraction. We also present another example showing the result illuminated bya white spot light source with 550nm bandpass filter (half width 10 nm), shownin Fig. 7.13 (center). Although the flux is not as strong as coherent laser illumi-nation, we observe a decoded image with less speckle noise (aka. blurred by theintegration over 10 nm wavelength). Readers may notice the “grid” effect, whichis a ghosting effect due to misalignment, as discussed previously. As expected,the full spectrum white example (right) suffers from severe blur because of colordispersion.To obtain a color image, one could encode different holograms for differentcolor primaries into the same DOE pair, similar to the parallax barrier setup [119].Nevertheless, it is challenging to obtain full-color display in high definition whentrading-off light efficiency, image size, spatial resolution, and system complexity.Multiple modulation pairings, similar to multi-modulator projectors, may be analternative. We leave this topic for future work.123Figure 7.13: Results of an offset pairing example illuminated by a coherencegreen light source (left), broadband green light source (center) and afull spectrum white light source (right). (Lena Photograph by DwightHooker)More than two layers. Another extension that we leave for future work is to in-crease the number of layers from two to three or more. This would change thematrix factorization problem in Eq. 7.9 to a tensor factorization problem, similarin spirit to the work by Wetzstein et al. [164] on tensor displays.7.5 ConclusionWe have demonstrated that it is viable to encode the holograms for many targetimages into a pair of diffractive optical elements, where different alignments of theDOEs generate different target images. This is achieved by exploring the largebandwidth of holographic approaches. Computationally, we combine iterativephase retrieval methods and complex matrix factorization for multiplexing holo-grams onto a pair of diffractive optical elements, subject to a designated alignmentscheme. While the method exhibits some artifacts, which we analyze in detail, theoverall image quality is very good for many applications.In the future, our solution can serve as the basis for further generalizations toarbitrary numbers (and thicknesses) of physical or air layers, as well as more gen-eral types of geometric alignment. The unique benefits of multiplexed hologramscan potentially inspire a wide range of application scenarios where the refractivecaustic designs or single layer diffractive designs may fail.124Applications in the security market can be very promising, like public-privateencryption visualization of identity purposes. We note, for this application, the ran-domness of holograms and the effort needed for precise alignment in fact strengththe security and privacy. Our solution can be a trigger for many people to startresearch on a variety of novel encryption modes. The display industry is anotherfield where our solution may pave the way for exploiting diffraction in a varietyof augmented reality (AR) designs. In addition, embedding holographic imagescan make architecture or industrial products more appealing. We believe that ani-mations and other dynamic effects enabled by multiplexed holography enable newand enriched user experiences in many of these fields.125Chapter 8ConclusionAs a classic topic in optics that has been studied for decades or even for cen-turies, imaging systems have evolved drastically from naive pinhole model intomore comprehensive optical solutions to better approximate the mammalian visionsystem. While significant progress has occurred in modern photography technol-ogy, developing more form factor-efficient and more image quality-powerful solu-tions continue to be highly desirable, especially with the rise of mobile imagingsystems such as smartphone cameras. Computational imaging serves as an essen-tial approach to link the two most important subjects — optics and computation toremove blur and recover sharp images during post-processing as well as for inspir-ing new form factors and functionalities of cameras.In this dissertation, we have presented our research on this topic, precisely, thejoint-design of diffractive optics and computational algorithms. With the proposedComputational Imaging with Diffractive Optics, we make steps toward a novelimaging modality to enable a variety of compact and encodable capture and displaysolutions for both scientific and consumer imaging applications.In addition, with respect to optics, we investigate the practicability of intro-ducing diffractive optics into broadband imaging, which provides more flexibilityfor engineering desirable optical functionalities with special form factor hardware.With respect to image reconstruction, throughout this dissertation we exploit differ-ent optimization methods together with different prior-based approaches that aimto remove artifacts and increase accuracy and efficiency. With respect to the intelli-126gent design space, we apply both the heuristic design and the modern optimizationfootprints subject to different design intentions. Specifically, the complex matrixfactorization is exploited to efficiently derive the phase maps of each layer for boththe encoded lenses and the multiplexed holograms.Throughout this work, developing both light and form factor-efficient opticsand computing-efficient optimization solvers for image restoration have been theprinciple tasks for a given computational imaging scenario. In the following, were-iterate the key findings of each chapter. After that, we highlight promising futuredirections that could build on this thesis work.In Chapter 3, we demonstrate the first attempt of introducing diffractive op-tics into consumer level photography solutions. Together with the validation ofphotolithography fabrication technique, this chapter serves as a proof-of-conceptproject under the high-level computational imaging framework.In Chapter 4, we demonstrate a pure diffractive optics camera that incorporatesthe optimization both in optical design and in image reconstruction. This is realizedby engineering desirable point spread function (PSF) behaviors along full visiblespectrum, further leading to obtaining high-quality computational imaging in fullcolor. The optical design step is formed as an optimization problem of updating thephase profile of an ahcromat, meanwhile the image reconstruction step is formedas a deconvolution problem of removing the residual achromatic blur.In Chapter 5 and Chapter 6, we expand the encoding from the spectral domain(Chapter 4) to the spatial domain. To this end, we first demonstrate an extendeddepth of field (DOF) imaging with multi-focal sweep imaging modality. This isrealized by fusing multiple focal powers onto one single element so as to obtainapproximated depth-invariant kernel behavior after the sweeping. The remainingresidual blur is tackled by the deconvolution algorithms. Further, we demonstratean efficient design scheme to encode multiple optical functionalities onto diffrac-tive optics, including refocusing and zooming effect. In particular, we developfocus tunable lenses with a pair of diffractive optics. Despite the optimization (i.e.self-calibrating deconvolution) in the image reconstruction step to remove aberra-tions, a complex matrix factorization scheme is specially introduced in deriving thepattern of each phase layer.Finally in Chapter 7, we extend the partial optimization scheme from Chap-127ter 6, but fits its practicability on display end. Precisely, we demonstrate that thecombination of controlled diffractive propagation and geometrical aligned holo-grams can enable a mix-and-match holography technique so as to deliver multipledifferent target images. The modified phase retrieval methods and the complexmatrix factorization are exploited.Again, this dissertation makes first promising steps toward the vision of thejoint design with diffractive optics in photography and visualization scenarios. Tomake this problem tractable at all, there are many assumptions being made aboutthe capture and the reconstruction in all considered cases. In addition, a part of thedesign scheme in optics relies heavily on the heuristic observations and updates,that can actually be guided by more advanced strategies like introducing machinelearning to exploit a paradigm. This could be a challenging but interesting directionthat we would like to dive into.Note that even with improved inverse methods, resolving imaging quality isfundamentally limited by the light efficiency and the content fidelity that can beeventually preserved by the sensor. We are not intended to completely replace re-fractive lenses in current high-end DSLR cameras, in terms of the cost and theimage quality. Instead, the designs investigated in this thesis dissertation shall pro-vide insights for alternative camera designs with special form factor requirementsor special functionality expectations.Overall, we envision the emergence between flat optics and formal numericaloptimization techniques will stimulate many novel ideas and successes in the fieldof imaging applications, ranging from capture to display end. This joint designtendency is definitely full of significance in both academia and industry, and Iincline to keep investigating on it as future work.128Bibliography[1] S. H. Ahn and L. J. Guo. Large-area roll-to-roll and roll-to-platenanoimprint lithography: a step toward high-throughput application ofcontinuous nanoimprinting. ACS Nano, 3(8):2304–2310, 2009. → pages25, 101, 122[2] F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, andF. Capasso. Aberration-free ultrathin flat lenses and axicons at telecomwavelengths based on plasmonic metasurfaces. Nano Letters, 12(9):4932–4936, 2012. → pages 67[3] F. Aieta, M. A. Kats, P. Genevet, and F. Capasso. Multiwavelengthachromatic metasurfaces by dispersive phase compensation. Science, 347(6228):1342–1345, 2015. → pages 14[4] N. Antipa, G. Kuo, R. Heckel, B. Mildenhall, E. Bostan, R. Ng, andL. Waller. Diffusercam: lensless single-exposure 3d imaging. Optica, 5(1):1–9, 2018. → pages 15[5] M. S. Asif, A. Ayremlou, A. Sankaranarayanan, A. Veeraraghavan, andR. Baraniuk. Flatcam: Thin, bare-sensor cameras using coded aperture andcomputation. arXiv preprint arXiv:1509.00116, 2015. → pages 83[6] I. Baran, P. Keller, D. Bradley, S. Coros, W. Jarosz, D. Nowrouzezahrai,and M. Gross. Manufacturing layered attenuators for multiple prescribedshadow images. In Computer Graphics Forum, volume 31, pages 603–610.Wiley Online Library, 2012. → pages 22[7] I. M. Barton, S. N. Dixit, L. J. Summers, C. A. Thompson, K. Avicola, andJ. Wilhelmsen. Diffractive alvarez lens. Optics Letters, 25(1):1–3, 2000.→ pages 85129[8] S. Bernet and M. Ritsch-Marte. Adjustable refractive power fromdiffractive moire´ elements. Applied Optics, 47(21):3722–3730, 2008. →pages 84, 90[9] S. Bernet, W. Harm, and M. Ritsch-Marte. Demonstration of focus-tunablediffractive moire´-lenses. Optics Express, 21(6):6955–6966, 2013. → pages11, 84, 90[10] F. Berthouzoz and R. Fattal. Resolution enhancement by vibrating displays.ACM Transactions on Graphics (TOG), 31(2):15, 2012. → pages 24[11] P.-A. Blanche, A. Bablumian, R. Voorakaranam, C. Christenson, W. Lin,T. Gu, D. Flores, P. Wang, W.-Y. Hsieh, M. Kathaperumal, et al.Holographic three-dimensional telepresence using large-areaphotorefractive polymer. Nature, 468(7320):80–83, 2010. → pages 3[12] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributedoptimization and statistical learning via the alternating direction method ofmultipliers. Foundations and Trends R© in Machine Learning, 3(1):1–122,2011. → pages 17, 18, 57[13] S. P. Boyd and L. Vandenberghe. Convex optimization. CambridgeUniversity Press, 2004. → pages 16[14] M. Brand. Specular holography. Applied Optics, 50(25):5042–5046, 2011.→ pages 24[15] E. Buckley. Holographic projector using one lens. Optics Letters, 35(20):3399–3401, 2010. → pages 23[16] W. T. Cathey and E. R. Dowski. New paradigm for imaging systems.Applied Optics, 41(29):6080–6092, 2002. → pages 90[17] A. Chakrabarti and T. Zickler. Statistics of real-world hyperspectralimages. In Proc. Computer Vision and Pattern Recognization, pages193–200. IEEE, 2011. → pages 60[18] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convexproblems with applications to imaging. Journal of Mathematical Imagingand Vision, 40(1):120–145, 2011. → pages 17, 153[19] S. H. Chan, R. Khoshabeh, K. B. Gibson, P. E. Gill, and T. Q. Nguyen. Anaugmented lagrangian method for total variation video restoration. IEEETransactions Image Processing (TIP), 20(11):3097–3111, 2011. → pages19130[20] T. S. Cho, C. L. Zitnick, N. Joshi, S. B. Kang, R. Szeliski, and W. T.Freeman. Image restoration by matching gradient distributions. IEEETransactions on Pattern Analysis and Machine Intelligence (TPAMI), 34(4):683–694, 2012. → pages 19[21] S. Y. Chou, P. R. Krauss, and P. J. Renstrom. Nanoimprint lithography.Journal of Vacuum Science & Technology B: Microelectronics andNanometer Structures Processing, Measurement, and Phenomena, 14(6):4129–4133, 1996. → pages 101[22] S.-W. Chung, B.-K. Kim, and W.-J. Song. Detecting and eliminatingchromatic aberration in digital images. In Image Processing (ICIP), 16thIEEE International Conference on, pages 3905–3908. IEEE, 2009. →pages 21[23] A. Cichocki and R. Zdunek. Multilayer nonnegative matrix factorizationusing projected gradient approaches. International Journal of NeuralSystems, 17(06):431–446, 2007. → pages 21, 109[24] O. Cossairt and S. Nayar. Spectral focal sweep: Extended depth of fieldfrom chromatic aberrations. In Computational Photography (ICCP), IEEEInternational Conference on, pages 1–8. IEEE, 2010. → pages 41[25] O. Cossairt, M. Gupta, and S. K. Nayar. When does computational imagingimprove performance? IEEE Transactions Image Processing (TIP), 22(2):447–458, 2013. → pages 68[26] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image denoising bysparse 3-d transform-domain collaborative filtering. IEEE Transactions onImage Processing (TIP), 16(8):2080–2095, 2007. → pages 16, 35[27] G. Damberg and W. Heidrich. Efficient freeform lens optimization forcomputational caustic displays. Optics Express, 23(8):10224–10232, 2015.→ pages 22, 25[28] G. Damberg, J. Gregson, and W. Heidrich. High brightness hdr projectionusing dynamic freeform lensing. ACM Transactions on Graphics (TOG),35(3):24, 2016. → pages 22, 101, 113[29] D. G. Dansereau, O. Pizarro, and S. B. Williams. Linear volumetric focusfor light field cameras. ACM Transactions on Graphics (TOG), 34(2):15–1,2015. → pages 69131[30] E. R. Dowski and W. T. Cathey. Extended depth of field throughwave-front coding. Applied Optics, 34(11):1859–1866, 1995. → pages 2[31] W. Duoshu, C. Luo, Y. Xiong, T. Chen, H. Liu, and J. Wang. Fabricationtechnology of the centrosymmetric continuous relief diffractive opticalelements. Physics Procedia, 18:95–99, 2011. → pages 25, 36[32] R. C. Eberhart and Y. Shi. Particle swarm optimization: developments,applications and resources. In Proc. of Congress on EvolutionaryComputation, pages 81–86. IEEE, 2001. → pages 53[33] D. Elkind, Z. Zalevsky, U. Levy, and D. Mendlovic. Optical transferfunction shaping and depth of focus by using a phase only filter. AppliedOptics, 42(11):1925–1931, 2003. → pages 11[34] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman.Removing camera shake from a single photograph. In ACM Transactionson Graphics (TOG), volume 25, pages 787–794. ACM, 2006. → pages 19[35] J. R. Fienup. Phase retrieval algorithms: a comparison. Applied Optics, 21(15):2758–2769, 1982. → pages 23, 104[36] M. Finckh, H. Dammertz, and H. P. Lensch. Geometry construction fromcaustic images. In Proc. European Conference on Computer Vision(ECCV), pages 464–477, 2010. → pages 22, 101[37] R. E. Fischer, B. Tadic-Galeb, P. R. Yoder, and R. Galeb. Optical systemdesign. McGraw Hill New York, 2000. → pages 60, 82[38] G. R. Fowles. Introduction to modern optics. Courier Dover Publications,2012. → pages 10[39] M. Fuchs, R. Raskar, H.-P. Seidel, and H. Lensch. Towards passive 6dreflectance field displays. ACM Transactions on Graphics (TOG), 27(3):58,2008. → pages 22[40] C. F. Gauss. Dioptrische Untersuchungen von CF Gauss. in derDieterichschen Buchhandlung, 1843. → pages 10[41] J. M. Geary. Introduction to lens design: with practical ZEMAX examples.Willmann-Bell Richmond, 2002. → pages 10[42] D. Geman and C. Yang. Nonlinear image recovery with half-quadraticregularization. IEEE Transactions on Image Processing (TIP), 4(7):932–946, 1995. → pages 17, 18132[43] P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin. Recentadvances in planar optics: from plasmonic to dielectric metasurfaces.Optica, 4(1):139–152, 2017. → pages 69[44] R. W. Gerchberg and W. Saxton. Phase determination from image anddiffraction plane pictures in electron-microscope. Optik, 34(3):275–+,1971. → pages 23, 104[45] T. D. Gerke and R. Piestun. Aperiodic volume optics. Nature Photonics, 4(3):188–193, 2010. → pages 24[46] P. R. Gill. Odd-symmetry phase gratings produce optical nulls uniquelyinsensitive to wavelength and depth. Optics Letters, 38(12):2074–2076,2013. → pages 83[47] P. R. Gill and D. G. Stork. Lensless ultra-miniature imagers usingodd-symmetry spiral phase gratings. In Computational Optical Sensingand Imaging, pages CW4C–3. Optical Society of America, 2013. → pages2, 14, 46[48] P. R. Gill and D. G. Stork. Phase gratings with odd symmetry forhigh-resolution lensed and lensless optical sensing, Aug. 18 2015. USPatent 9,110,240. → pages 83[49] D. Glasner, T. Zickler, and A. Levin. A reflectance display. ACMTransactions on Graphics (TOG), 33(4):61, 2014. → pages 23[50] B. Goldluecke and D. Cremers. Introducing total curvature for imageprocessing. In Computer Vision (ICCV), IEEE International Conferenceon, pages 1267–1274. IEEE, 2011. → pages 16[51] T. Goldstein and S. Osher. The split bregman method for l1-regularizedproblems. SIAM Journal on Imaging Sciences, 2(2):323–343, 2009. →pages 17, 74[52] J. Goodman. Introduction to Fourier optics. McGraw-hill, 2008. → pages12, 27, 30, 48[53] A. Grewe, M. Hillenbrand, and S. Sinzinger. Aberration analysis ofoptimized alvarez–lohmann lenses. Applied Optics, 53(31):7498–7506,2014. → pages 85[54] W. Harm, C. Roider, A. Jesacher, S. Bernet, and M. Ritsch-Marte.Dispersion tuning with a varifocal diffractive-refractive hybrid lens. OpticsExpress, 22(5):5260–5269, 2014. → pages 83133[55] S. W. Hasinoff and K. N. Kutulakos. Light-efficient photography. IEEETransactions on Pattern Analysis and Machine Intelligence, 33(11):2203–2214, 2011. → pages 27[56] F. Heide. Structure-aware computational imaging. PhD thesis, TheUniversity of British Columbia, 2016. → pages 3[57] F. Heide, M. Rouf, M. B. Hullin, B. Labitzke, W. Heidrich, and A. Kolb.High-quality computational imaging through simple lenses. ACMTransactions on Graphics (TOG), 32(5):149, 2013. → pages 2, 21, 27, 34,38, 46, 56, 83, 87, 151, 153[58] F. Heide, J. Gregson, G. Wetzstein, R. Raskar, and W. Heidrich.Compressive multi-mode superresolution display. Optics Express, 22(12):14981–14992, 2014. → pages 23[59] F. Heide, D. Lanman, D. Reddy, J. Kautz, K. Pulli, and D. Luebke.Cascaded displays: spatiotemporal superresolution using offset pixellayers. ACM Transactions on Graphics (TOG), 33(4):60, 2014. → pages24, 90, 108[60] F. Heide, M. Steinberger, Y.-T. Tsai, M. Rouf, D. Pajak, D. Reddy,O. Gallo, J. Liu, W. Heidrich, K. Egiazarian, et al. Flexisp: A flexiblecamera image processing framework. ACM Transactions on Graphics(TOG), 33(6):231, 2014. → pages 34, 35[61] F. Heide, Q. Fu, Y. Peng, and W. Heidrich. Encoded diffractive optics forfull-spectrum computational imaging. Scientific Reports, 6, 2016. → pages74, 89, 111[62] N.-D. Ho. Nonnegative matrix factorization algorithms and applications.PhD thesis, PhD thesis, Universite´ catholique de Louvain, 2008. → pages21, 86, 108, 109[63] S. Honnungar, J. Holloway, A. K. Pediredla, A. Veeraraghavan, andK. Mitra. Focal-sweep for large aperture time-of-flight cameras. In Proc.International Cconference of Image Processing (ICIP), 2016. → pages 69[64] R. Hoskinson, B. Stoeber, W. Heidrich, and S. Fels. Light reallocation forhigh contrast projection using an analog micromirror array. ACMTransactions on Graphics (TOG), 29(6):165, 2010. → pages 101134[65] R. Hostettler, R. Habel, M. Gross, and W. Jarosz. Dispersion-based colorprojection using masked prisms. In Computer Graphics Forum, volume 34,pages 329–338. Wiley Online Library, 2015. → pages 22[66] P. O. Hoyer. Non-negative matrix factorization with sparseness constraints.Journal of Machine Learning Research, 5(Nov):1457–1469, 2004. →pages 21[67] M. B. Hullin, I. Ihrke, W. Heidrich, T. Weyrich, G. Damberg, andM. Fuchs. Computational fabrication and display of material appearance.In Eurographics 2013-State of the Art Reports, pages 137–153, 2012. →pages 22[68] D. Iwai, S. Mihara, and K. Sato. Extended depth-of-field projector by fastfocal sweep projection. IEEE Transactions on Visualization and ComputerGraphics (TVCG), 21(4):462–470, 2015. → pages 69[69] J. C. Javitt and R. F. Steinert. Cataract extraction with multifocalintraocular lens implantation: a multinational clinical trial evaluatingclinical, functional, and quality-of-life outcomes. Ophthalmology, 107(11):2040–2048, 2000. → pages 69[70] N. Joshi, R. Szeliski, and D. Kriegman. Psf estimation using sharp edgeprediction. In Proc. Computer Vision and Pattern Recognization (CVPR),pages 1–8. IEEE, 2008. → pages 20[71] G. Kang, Q. Tan, X. Wang, and G. Jin. Achromatic phase retarder appliedto mwir & lwir dual-band. Optics Express, 18(2):1695–1703, 2010. →pages 67[72] S. B. Kang. Automatic removal of chromatic aberration from a singleimage. In Computer Vision and Pattern Recognition (CVPR), 2007. IEEEConference on, pages 1–8. IEEE, 2007. → pages 21[73] R. H. Keates, J. L. Pearce, and R. T. Schneider. Clinical results of themultifocal lens. Journal of Cataract & Refractive Surgery, 13(5):557–560,1987. → pages 69[74] M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, andF. Capasso. Metalenses at visible wavelengths: Diffraction-limitedfocusing and subwavelength resolution imaging. Science, 352(6290):1190–1194, 2016. → pages 14135[75] G. Kim, J. A. Domı´nguez-Caballero, and R. Menon. Design and analysisof multi-wavelength diffractive optics. Optics Express, 20(3):2814–2823,2012. → pages 23[76] R. Kingslake and R. B. Johnson. Lens design fundamentals. AcademicPress, 2009. → pages 10, 82[77] T. Kiser, M. Eigensatz, M. M. Nguyen, P. Bompas, and M. Pauly.Architectural Caustics—Controlling Light with Geometry. Springer, 2013.→ pages 22, 101[78] A. Kolodziejczyk and Z. Jaroszewicz. Diffractive elements of variableoptical power and high diffraction efficiency. Applied Optics, 32(23):4317–4322, 1993. → pages 85[79] D. Krishnan and R. Fergus. Fast image deconvolution usinghyper-laplacian priors. In Advances in Neural Information ProcessingSystems, pages 1033–1041. NIPS, 2009. → pages 19, 57, 58, 60[80] D. Krishnan, T. Tay, and R. Fergus. Blind deconvolution using anormalized sparsity measure. In Proc. Computer Vision and PatternRecognition, pages 233–240. IEEE, 2011. → pages 20[81] S. Kuthirummal, H. Nagahara, C. Zhou, and S. K. Nayar. Flexible depth offield photography. IEEE Transactions on Pattern Analysis and MachineIntelligence (TPAMI), 33(1):58–71, 2011. → pages 69[82] D. Lanman, G. Wetzstein, M. Hirsch, W. Heidrich, and R. Raskar.Polarization fields: dynamic light field display using multi-layer lcds. ACMTransactions on Graphics (TOG), 30(6):186, 2011. → pages 23[83] M. Lebrun. An analysis and implementation of the bm3d image denoisingmethod. Image Processing On Line, (2012), 2012. → pages 16[84] D. D. Lee and H. S. Seung. Algorithms for non-negative matrixfactorization. In Advances in neural information processing systems, pages556–562, 2001. → pages 21[85] S. Lee, C. Jang, S. Moon, J. Cho, and B. Lee. Additive light field displays:realization of augmented reality with holographic optical elements. ACMTransactions on Graphics (TOG), 35(4):60, 2016. → pages 23[86] L. Lesem, P. Hirsch, and J. Jordan. The kinoform: a new wavefrontreconstruction device. IBM Journal of Research and Development, 13(2):150–155, 1969. → pages 23136[87] A. Levin, R. Fergus, F. Durand, and W. T. Freeman. Image and depth froma conventional camera with a coded aperture. ACM Transactions onGraphics (TOG), 26(3):70, 2007. → pages 1, 19, 32, 69[88] A. Levin, S. W. Hasinoff, P. Green, F. Durand, and W. T. Freeman. 4dfrequency analysis of computational cameras for depth of field extension.In ACM Transactions on Graphics (TOG), volume 28, page 97. ACM,2009. → pages 2[89] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman. Understanding andevaluating blind deconvolution algorithms. In Computer Vision and PatternRecognition (CVPR), pages 1964–1971. IEEE, 2009. → pages 20[90] A. Levin, D. Glasner, Y. Xiong, F. Durand, W. Freeman, W. Matusik, andT. Zickler. Fabricating brdfs at high spatial resolution using wave optics.ACM Transactions on Graphics (TOG), 32(4):144, 2013. → pages 23[91] C.-b. Lin. Projected gradient methods for nonnegative matrix factorization.Neural Computation, 19(10):2756–2779, 2007. → pages 21, 109[92] A. W. Lohmann and D. Paris. Binary fraunhofer holograms, generated bycomputer. Applied Optics, 6(10):1739–1748, 1967. → pages 23[93] L. B. Lucy. An iterative technique for the rectification of observeddistributions. The Astronomical Journal, 79:745, 1974. → pages 19[94] D. Malacara-Herna´ndez and Z. Malacara-Herna´ndez. Handbook of opticaldesign. CRC Press, 2016. → pages 10[95] D. L. Marks, R. A. Stack, D. J. Brady, and J. van der Gracht.Three-dimensional tomography using a cubic-phase plate extendeddepth-of-field system. Optics Letters, 24(4):253–255, 1999. → pages 83[96] D. Martin, C. Fowlkes, D. Tal, and J. Malik. A database of humansegmented natural images and its application to evaluating segmentationalgorithms and measuring ecological statistics. In Proc. of the IEEEInternational Conference on Computer Vision (ICCV), pages 416–423.IEEE, 2001. → pages 75[97] K. Marwah, G. Wetzstein, Y. Bando, and R. Raskar. Compressive light fieldphotography using overcomplete dictionaries and optimized projections.ACM Transactions on Graphics (TOG), 32(4):46, 2013. → pages 2137[98] J. H. McLeod. The axicon: a new type of optical element. Journal of OSA,44(8):592–597, 1954. → pages 90[99] X. Meng, L. Cai, Y. Wang, X. Yang, X. Xu, G. Dong, X. Shen, H. Zhang,and X. Cheng. Hierarchical image encryption based on cascaded iterativephase retrieval algorithm in the fresnel domain. Journal of Optics A: Pureand Applied Optics, 9(11):1070, 2007. → pages 24[100] M. M. Meyers. Hybrid refractive/diffractive achromatic camera lens, Feb. 31998. US Patent 5,715,091. → pages 32, 64[101] D. Miau, O. Cossairt, and S. K. Nayar. Focal sweep videography withdeformable optics. In Proc. International Conference on ComputationalPhotography (ICCP), pages 1–8, 2013. → pages 69[102] P. Milanfar. A tour of modern image filtering: New insights and methods,both practical and theoretical. IEEE Signal Processing Magazine, 30(1):106–128, 2013. → pages 3[103] J. Mitchell, M. Francke, and D. Eng. Illustrative rendering in team fortress2. In Proceedings of the 5th international symposium on Non-photorealisticanimation and rendering, pages 71–76. ACM, 2007. → pages 118[104] F. H. Mok. Angle-multiplexed storage of 5000 holograms in lithiumniobate. Optics Letters, 18(11):915–917, 1993. → pages 24[105] M. Monjur, L. Spinoulas, P. R. Gill, and D. G. Stork. Ultra-miniature,computationally efficient diffractive visual-bar-position sensor. In Proc.SensorComm. IEIFSA, 2015. → pages 14[106] Y. Nagahama, T. Shimobaba, T. Kawashima, T. Kakue, and T. Ito.Holographic multi-projection using the random phase-free method.Applied Optics, 55(5):1118–1123, 2016. → pages 23, 106[107] T. Nakai and H. Ogawa. Research on multi-layer diffractive opticalelements and their application to camera lenses. In Diffractive Optics andMicro-Optics, page DMA2. OSA, 2002. → pages 14, 32[108] M. Naor and A. Shamir. Visual cryptography. In Workshop on the Theoryand Application of of Cryptographic Techniques, pages 1–12. Springer,1994. → pages 102, 115[109] R. Ng, M. Levoy, M. Bre´dif, G. Duval, M. Horowitz, and P. Hanrahan.Light field photography with a hand-held plenoptic camera. CSTR, 2(11):1–11, 2005. → pages 69138[110] A. Nikonorov, R. Skidanov, V. Fursov, M. Petrov, S. Bibikov, andY. Yuzifovich. Fresnel lens imaging with post-capture image processing. InProc. of the IEEE Conference on Computer Vision and Pattern Recognition(CVPR) Workshops, pages 33–41, 2015. → pages 14[111] N. Okada, T. Shimobaba, Y. Ichihashi, R. Oi, K. Yamamoto, M. Oikawa,T. Kakue, N. Masuda, and T. Ito. Band-limited double-step fresneldiffraction and its application to computer-generated holograms. OpticsExpress, 21(7):9192–9197, 2013. → pages 105, 106[112] D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather. Diffractiveoptics: design, fabrication, and test. 2004. → pages 11, 29, 83[113] M. Papas, W. Jarosz, W. Jakob, S. Rusinkiewicz, W. Matusik, andT. Weyrich. Goal-based caustics. In Computer Graphics Forum,volume 30, pages 503–511, 2011. → pages 22, 101[114] M. Papas, T. Houit, D. Nowrouzezahrai, M. Gross, and W. Jarosz. Themagic lens: refractive steganography. ACM Transactions on Graphics(TOG), 31(6):186, 2012. → pages 22[115] N. Parikh, S. Boyd, et al. Proximal algorithms. Foundations and Trends R©in Optimization, 1(3):127–239, 2014. → pages 17[116] S. R. P. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J.Twieg, R. Piestun, and W. Moerner. Three-dimensional, single-moleculefluorescence imaging beyond the diffraction limit by using a double-helixpoint spread function. Proc. of the National Academy of Sciences, 106(9):2995–2999, 2009. → pages 83[117] Y. Peng, Q. Fu, H. Amata, S. Su, F. Heide, and W. Heidrich. Computationalimaging using lightweight diffractive-refractive optics. Optics Express, 23(24):31393–31407, 2015. → pages 2, 11, 14, 46, 63, 69[118] Y. Peng, Q. Fu, F. Heide, and W. Heidrich. The diffractive achromat fullspectrum computational imaging with diffractive optics. ACM Transactionson Graphics (TOG), 35(4):31, 2016. → pages 3, 61, 69, 75, 76, 112[119] Y. Peng, X. Dun, Q. Sun, and W. Heidrich. Mix-and-match holography.ACM Transactions on Graphics (TOG), 2017. → pages 118, 123[120] D. Perrone and P. Favaro. Total variation blind deconvolution: The devil isin the details. In Proc. of the IEEE Conference on Computer Vision andPattern Recognition (CVPR), pages 2909–2916, 2014. → pages 150139[121] D. Psaltis and G. W. Burr. Holographic data storage. Computer, 31(2):52–60, 1998. → pages 24[122] D. Psaltis, K. Curtis, M. Levene, A. Pu, and G. Barbastathis. Holographicstorage using shift multiplexing. Optics Letters, 20(7):782–784, 1995. →pages 23[123] W. Qu, H. Gu, H. Zhang, and Q. Tan. Image magnification in lenslessholographic projection using double-sampling fresnel diffraction. AppliedOptics, 54(34):10018–10021, 2015. → pages 23, 106[124] S. Quirin and R. Piestun. Depth estimation and image recovery usingbroadband, incoherent illumination with engineered point spread functions[invited]. Applied Optics, 52(1):A367–A376, 2013. → pages 14[125] R. Raskar, A. Agrawal, and J. Tumblin. Coded exposure photography:motion deblurring using fluttered shutter. In ACM Transactions onGraphics (TOG), volume 25, pages 795–804. ACM, 2006. → pages 2[126] W. H. Richardson. Bayesian-based iterative method of image restoration.Journal of OSA, 62(1):55–59, 1972. → pages 19[127] G. Rougeron and B. Pe´roche. Color fidelity in computer graphics: Asurvey. In Computer Graphics Forum, volume 17, pages 3–15. WileyOnline Library, 1998. → pages 20[128] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noiseremoval algorithms. Physica D: Nonlinear Phenomena, 60(1-4):259–268,1992. → pages 16, 19[129] B. Sajadi, M. Gopi, and A. Majumder. Edge-guided resolutionenhancement in projectors via optical pixel sharing. ACM Trans. Graph, 31(4):79, 2012. → pages 24[130] E. Samei, M. J. Flynn, and D. A. Reimann. A method for measuring thepresampled mtf of digital radiographic systems using an edge test device.Medical physics, 25(1):102–113, 1998. → pages 63[131] C. J. Schuler, M. Hirsch, S. Harmeling, and B. Scho¨lkopf. Non-stationarycorrection of optical aberrations. In Proc. International Conference onComputer Vision (ICCV), pages 659–666, 2011. → pages 2, 21, 46, 83[132] C. J. Schuler, M. Hirsch, S. Harmeling, and B. Scho¨lkopf. Blind correctionof optical aberrations. In Proc. European Conference on Computer Vision(ECCV), pages 187–200, 2012. → pages 21140[133] C. J. Schuler, H. C. Burger, S. Harmeling, and B. Scholkopf. A machinelearning approach for non-blind image deconvolution. In Proc. ComputerVision and Pattern Recognition (CVPR), pages 1067–1074. IEEE, 2013. →pages 20, 58, 60[134] Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly. High-contrastcomputational caustic design. ACM Transactions on Graphics (TOG), 33(4):74, 2014. → pages 22[135] L. Seidel. Ueber die theorie der fehler. mit welchen die durch optischeInstrumente gesehenen Bilder behaftet sind, und u¨ber die mathematischenBedingungen ihrer Aufhebung. Abhandlungen derNaturwissenschaftlich-Technischen Commission bei der Ko¨nigl.Bayerischen Akademie der Wissenschaften in Mu¨nchen. Cotta, 2:4, 1857.→ pages 10[136] Q. Shan, J. Jia, and A. Agarwala. High-quality motion deblurring from asingle image. ACM Transactions on Graphics (TOG), 27(3):73, 2008. →pages 20[137] G. Sharma. Color fundamentals for digital imaging. Digital color imaginghandbook, 20, 2003. → pages 20[138] Y. Shi, G. Situ, and J. Zhang. Multiple-image hiding in the fresnel domain.Optics Letters, 32(13):1914–1916, 2007. → pages 24[139] T. Shimobaba, M. Makowski, T. Kakue, M. Oikawa, N. Okada, Y. Endo,R. Hirayama, and T. Ito. Lensless zoomable holographic projection usingscaled fresnel diffraction. Optics Express, 21(21):25285–25290, 2013. →pages 104[140] M. Singh, J. Tervo, and J. Turunen. Broadband beam shaping withharmonic diffractive optics. Optics Express, 22(19):22680–22688, 2014.→ pages 14[141] S. Sinzinger, J. Jahns, V. R. Daria, and J. Glu¨ckstad. Planar microopticalsystems for correlation and security applications. In Optical ImagingSensors and Systems for Homeland Security Applications, pages 339–366.2006. → pages 25[142] T. Skauli and J. Farrell. A collection of hyperspectral images for imagingsystems research. In Proc. IS&T/SPIE Electronic Imaging, pages86600C–86600C. SPIE, 2013. → pages 60141[143] G. G. Sliusarev. Aberration and optical design theory. Bristol, England,Adam Hilger, Ltd., 1984, 672 p. Translation., 1984. → pages 1, 10[144] W. J. Smith. Modern optical engineering. Tata McGraw-Hill Education,1966. → pages 83[145] W. J. Smith. Modern lens design. McGraw-Hill, 2005. → pages 20, 96[146] D. G. Stork and P. R. Gill. Lensless ultra-miniature cmos computationalimagers and sensors. Proc. SENSORCOMM, pages 186–190, 2013. →pages 83[147] D. G. Stork and P. R. Gill. Optical, mathematical, and computationalfoundations of lensless ultra-miniature diffractive imagers and sensors.International Journal on Advances in Systems and Measurements, 7(3):4,2014. → pages 2, 46[148] Z. Su, K. Zeng, L. Liu, B. Li, and X. Luo. Corruptive artifacts suppressionfor example-based color transfer. IEEE Transactions on Multimedia, 16(4):988–999, 2014. → pages 20[149] L. Sun, S. Cho, J. Wang, and J. Hays. Edge-based blur kernel estimationusing patch priors. In Proc. International Conference on ComputationalPhotography (ICCP), pages 1–8, 2013. → pages 20, 152[150] T. Sun, Y. Peng, and W. Heidrich. Revisiting cross-channel informationtransfer for chromatic aberration correction. In Proc. InternationalConference on Computer Vision (ICCV), pages 3248–3256, 2017. → pages21[151] G. J. Swanson. Binary optics technology: theoretical limits on thediffraction efficiency of multilevel diffractive optical elements. Technicalreport, Massachusetts Insititute of Tech Lexington Lincoln Lab, 1991. →pages 29[152] D. W. Sweeney and G. E. Sommargren. Harmonic diffractive lenses.Applied Optics, 34(14):2469–2475, 1995. → pages 42[153] M. W. Tao, S. Hadap, J. Malik, and R. Ramamoorthi. Depth fromcombining defocus and correspondence using light-field cameras. In Proc.International Conference on Computer Vision (ICCV), pages 673–680,2013. → pages 68142[154] G. Tricoles. Computer generated holograms: an historical review. AppliedOptics, 26(20):4351–4360, 1987. → pages 23[155] S. C. Tucker, W. T. Cathey, and E. R. Dowski. Extended depth of field andaberration control for inexpensive digital microscope systems. OpticsExpress, 4(11):467–474, 1999. → pages 11[156] A. Veeraraghavan, R. Raskar, A. Agrawal, A. Mohan, and J. Tumblin.Dappled photography: Mask enhanced cameras for heterodyned light fieldsand coded aperture refocusing. ACM Transactions on Graphics (TOG), 26(3):69, 2007. → pages 2[157] B. Wahlberg, S. Boyd, M. Annergren, and Y. Wang. An admm algorithmfor a class of total variation regularized estimation problems. arXivpreprint arXiv:1203.1828, 2012. → pages 18[158] P. Wang and R. Menon. Computational spectrometer based on a broadbanddiffractive optic. Optics Express, 22(12):14575–14587, 2014. → pages 11[159] P. Wang and R. Menon. Ultra-high-sensitivity color imaging via atransparent diffractive-filter array and computational optics. Optica, 2(11):933–939, 2015. → pages 2, 14[160] P. Wang, N. Mohammad, and R. Menon. Chromatic-aberration-correcteddiffractive lenses for ultra-broadband focusing. Scientific Reports, 6, 2016.→ pages 14[161] Y. Wang, W. Yun, and C. Jacobsen. Achromatic fresnel optics for widebandextreme-ultraviolet and x-ray imaging. Nature, 424(6944):50–53, 2003. →pages 67[162] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli. Image qualityassessment: from error visibility to structural similarity. Image Processing,IEEE Transactions on, 13(4):600–612, 2004. → pages 120[163] H. A. Weeber. Diffractive multifocal lens having radially varying lightdistribution, 2011. US Patent 7,871,162. → pages 69[164] G. Wetzstein, D. Lanman, M. Hirsch, and R. Raskar. Tensor displays:compressive light field synthesis using multilayer displays with directionalbacklighting. ACM Transactions on Graphics (TOG), 31(4):80, 2012. →pages 23, 119, 124143[165] G. Wyszecki and W. S. Stiles. Color science, volume 8. Wiley New York,1982. → pages 2[166] L. Xiao, J. Wang, W. Heidrich, and M. Hirsch. Learning high-order filtersfor efficient blind deconvolution of document photographs. In EuropeanConference on Computer Vision, pages 734–749. Springer, 2016. → pages20[167] L. Xu, S. Zheng, and J. Jia. Unnatural l0 sparse representation for naturalimage deblurring. In Proc. of the IEEE conference on Computer Vision andPattern Recognition, pages 1107–1114, 2013. → pages 20[168] G. Xue, J. Liu, X. Li, J. Jia, Z. Zhang, B. Hu, and Y. Wang. Multiplexingencoding method for full-color dynamic 3d holographic display. OpticsExpress, 22(15):18473–18482, 2014. → pages 23[169] M. Yamaguchi, H. Haneishi, and N. Ohyama. Beyond red–green–blue(rgb): Spectrum-based color imaging technology. JIST, 52(1):10201–1,2008. → pages 20[170] G. Ye, S. Jolly, V. M. Bove Jr, Q. Dai, R. Raskar, and G. Wetzstein. Towardbxdf display using multilayer diffraction. ACM Transactions on Graphics(TOG), 33(6):191, 2014. → pages 23[171] P. R. Yoder Jr. Opto-mechanical systems design. CRC press, 2005. →pages 83, 96[172] R. Yokoya and S. K. Nayar. Extended depth of field catadioptric imagingusing focal sweep. In Proc. International Conference on Computer Vision(ICCV), pages 3505–3513, 2015. → pages 69[173] L. Yuan, J. Sun, L. Quan, and H.-Y. Shum. Image deblurring withblurred/noisy image pairs. In ACM Transactions on Graphics (TOG),volume 26, page 1. ACM, 2007. → pages 19[174] L. Yuan, J. Sun, L. Quan, and H.-Y. Shum. Progressive inter-scale andintra-scale non-blind image deconvolution. In Acm Transactions onGraphics (TOG), volume 27, page 74. ACM, 2008. → pages 19, 36, 56, 57[175] T. Yue, S. Cho, J. Wang, and Q. Dai. Hybrid image deblurring by fusingedge and power spectrum information. In European Conference onComputer Vision, pages 79–93. Springer, 2014. → pages 20144[176] T. Yue, J. Suo, J. Wang, X. Cao, and Q. Dai. Blind optical aberrationcorrection by exploring geometric and visual priors. In Proc. ComputerVision and Pattern Recognition, 2015. → pages 21[177] T. Yue, J. Suo, J. Wang, X. Cao, and Q. Dai. Blind optical aberrationcorrection by exploring geometric and visual priors. In Proc. of the IEEEConference on Computer Vision and Pattern Recognition (CVPR), pages1684–1692, 2015. → pages 83[178] Y. Yue, K. Iwasaki, B.-Y. Chen, Y. Dobashi, and T. Nishita. Pixel art withrefracted light by rearrangeable sticks. Computer Graphics Forum, 31(2pt3):575–582, 2012. → pages 22[179] Y. Yue, K. Iwasaki, B.-Y. Chen, Y. Dobashi, and T. Nishita. Poisson-basedcontinuous surface generation for goal-based caustics. ACM Transactionson Graphics (TOG), 2014. → pages 22[180] C. Zhou and S. Nayar. What are good apertures for defocus deblurring? InComputational Photography (ICCP), IEEE International Conference on,pages 1–8. IEEE, 2009. → pages 2[181] C. Zhou, S. Lin, and S. Nayar. Coded aperture pairs for depth fromdefocus. In Proc. International Conference on Computer Vision (ICCV),pages 325–332. IEEE, 2009. → pages 69[182] G. Zhou, Y. Chen, Z. Wang, and H. Song. Genetic local search algorithmfor optimization design of diffractive optical elements. Applied Optics, 38(20):4281–4290, 1999. → pages 23[183] J. Zhou, L. Li, N. Naples, T. Sun, and Y. Y. Allen. Fabrication ofcontinuous diffractive optical elements using a fast tool servo diamondturning process. Journal of Micromechanics and Microengineering, 23(7):075010, 2013. → pages 25, 36145Appendix AAppendicesThis appendix includes the supporting materials that are closely relevant to theindividual chapter in main text.A.1 Supplementary details of fabricationA.1.1 Specifying a Fresnel lens.As has been discussed in Chapter 3, 4, without loss of generality, we take anamplitude-type Fresnel zone plate (FZP) as the instance to specify a Fresnel lens.Using scalar diffraction theory, the complex amplitude of a point at the image planeis the sum-up of all complex amplitudes derived from all sub-waves emitted fromring-zones of a FZP. To obtain constructive interference at the image plane, thezones should switch from the opaque to transparent at the radius given by:Rn =√nλ f +n2λ 24, (A.1)where n is the integer for the nth ring, λ is the wavelength of the illumination thezone plate is meant to focus, and f is the distance from the center of the zone plateto the focal point. In this case, each two neighbouring ring-zones holds an opticalpath difference of half wavelength (corresponding to the phase modulation of 2pi),such that they contribute the positive and negative amplitude to the spatial point,146respectively. Those negative components are blocked.With respect to the phase modulation, we expect the incident light to travelthrough the substrate medium with different thickness, such that create the requiredoptical path difference without blocking the light. That is to say, we can use thesame math described above to derive the radius for each ring zone, and then assignthem with different physical heights, corresponding to the phase delay they modu-late on the incident light rays. These height change can be designed continuous ordiscrete subject to the fabrication method we intend to apply.A.1.2 Fabrication procedure illustration.The DOEs are fabricated by repeatedly applying the lithography technique fol-lowed by a reactive ion etching step. Each round of fabrication results in a binaryprofile on the substrate. Therefore, by repeating the process for N times, a 2N-levelphase plate can be fabricated. The fabrication process is shown in Fig. A.1.We use four inch Fused Silica wafers with 0.5 mm thickness as the substrates.In each fabrication cycle, a 200 nm Cr layer is first deposited on the wafer. A0.6 µm photoresist layer (AZ1505) is then spin-coated on the Cr layer and gainsits shape after a softbake process (120 ◦C for 60 seconds). The designed patternsare transferred from the photomask to the photoresist under ultraviolet (UV) lightexposure. After exposure, the chemical property of the exposed area on the pho-toresist changes and can consequently be removed in the developer (MIL AZ726)for 20 seconds. Subsequently, the open area of the Cr layer is removed in Cr etchantand the patterns are transferred to the wafer. In the RIE step, the material (SiO2) inthe open area is removed by a mixture of Argon and SF6 plasma. Each fabricationcycle doubles the number of microstructure levels on the previous profile. A finalCr layer is deposited and etched to cover the outer area of the lens for preventionof stray light. This receipt holds for DOEs used in all the previous chapters.In our designs, the principle wavelength is λ0 = 550 nm. A total 2pi phasedelay at this wavelength corresponds to a depth of d = λ/(nλ − 1) = 1,195 nm,where nλ = 1.4599 is the refractive index of the wafer at the principle wavelength.Therefore, each level has a depth of 75 nm in our designs. A final depth error of ±10 nm has been achieved. As an instance, the microscope images of the fabricated147Figure A.1: The fabrication of DOEs consists of two major steps — lithogra-phy and etching. In the lithography step, patterns are transferred fromthe mask to photoresist on the wafer through exposure to UV light. Inthe etching step, the transferred patterns are converted to a binary pro-file on the wafer by reactive ion beam bombardment.Figure A.2: Illustration of 5X microscope images of the fabricated Fresnellens (left), and diffractive achromat (right). Note that the bright anddark changes indicate the physical height distribution.148Figure A.3: Illustration of 20X microscope image of a fabricated hologram.DOEs, a regular Fresnel lens (Chapter 3) and our diffractive achromat (Chapter 4)are shown in Fig. A.2, and a hologram (Chapter 7) is shown in Fig. A.3. We seethat a Fresnel lens has more regular periodical ring distribution, while a diffractiveachromat has more non-periodical ring distribution, corresponding to the optimizedsurface profile.A.2 Supplementary details of image reconstructionalgorithmsThis appendix section covers the image reconstruction that jointly estimates theunderlying scene-dependent, spatially-varying PSFs and the latent sharp image,see Chapter 6. Our method can be understood as a self-calibrating reconstructionmethod that adapts to the scene. Note that the cross-scale based non-blind de-convolution in Chapter 4 and the cross-channel based non-blind deconvolution inChapter 5 share a similar solver manner.149A.2.1 Blind deconvolution.Recall the discussion in Chapter 6, the PSFs of our DOEs change based on thereflectance of the objects and illumination spectra in the scene. Having ensuredthat our optics focus at least well in one channel, we can solve for the PSFs ex-ploiting cross-channel correlation between the channels. In particular we solve thefollowing optimization problem:xopt = argminx,v‖Vx− j‖22+α3∑a=1‖Hax‖22+β‖x‖22+γ∑a6=r‖Hav−Hair‖1+µ∑a‖Hav‖1.(A.2)The matrix V expresses here the convolution with the latent image v. Thatmeans the results of the matrix-vector product Vx is the same as the vectorizedconvolution v ∗ x. Convolution is a symmetric operator, hence it can be inter-changeably expressed as either Vx or Xv, where V and X are the correspondingconvolution matrices. The operators Ha (a = 1,2,3) are first order spatial gradientfilters for RGB channels, and ir is a sharp image in reference channel r.In the first row we have a classical `2 dataterm and gradient priors on the PSFx. However, this simple term is not sufficient to get robust blind estimation quality.That is why recent approaches use complex priors on the latent image. They areactually not just solving a simple quadratic problem, but a complex sequence ofindividual heuristic steps. A solid discussion of different heuristics behind blinddeconvolution can be found in [120].Cross-channel prior and algorithm. We add a cross-channel term on the latentimage v. This cross-channel term is derived from prior statistical knowledge aboutthe correlation between different spectral bands, recall the discussion in Chapter 3and Chapter 6. We model gradient differences with a heavy-tailed prior, expressedas follows:p(xopt|γ, ir) ∝ exp(−γ∑a‖Hav−Hair‖1). (A.3)The key observation of our cross-channel prior is that for a latent (sharp) image,changes in chroma and luma are sparsely distributed in natural images. That means,150in most areas the gradients of different channels are close, except for few largechanges in chroma and luma (e.g. commonly at object or material boundaries).Note that this model differs from the one proposed in [57]. The authors as-sume only chroma changes to be sparse. The downside is that this leads to se-vere instabilities of the prior term for low intensities. The cross-channel term inEq. (A.2) follows from Maximum-a-Posteriori (MAP) Estimation with the priorfrom Eq. (A.3). We use the above prior for the designs in Chapter 6 and Chapter 5.Recall Alg. 2 of Chapter 6 both sub-problems contain simple linear operatorsXk and Vk that represent convolutions. This formulation can be used since convo-lutions are commutative:Vx = Xv since x∗v = v∗x. (A.4)The optimization solves the bi-convex minimization problem via coordinatedescent. That means we keep one of the two variables x,v fixed at a time whileminimizing the objective with respect to the other variable in an alternating fashion.This approach leads to the two sub-problems from the steps x-step and v-step ofalgorithm Alg. 2 of Chapter 6, both of which are now much easier to solve. All ofthese sub-steps are now convex optimization problems, while the joint objective isnon-convex, shown as in Eq. A.5:V′ = VE and x′ = E−1x since,V′x′ = VEE−1x = VIx.(A.5)Having transformed our non-convex problem into a sequence of convex sub-problems, we can derive efficient solutions for these individual sub-problems.151Solving the x-step and the v-step. On the one hand, the x-step can be solvedefficiently as follows:xopt = argminx‖Vx− j‖22+α∑a‖Hax‖22+β‖x‖22︸ ︷︷ ︸χ(x)⇔ ∂χ(xopt)∂x= 2(VT Vxopt−VT j)+α∑a2HTa Haxopt+2βxopt!= 0⇔ VT Vxopt+α∑aHTa Haxopt+βxopt = VT j⇔(VT V+α∑aHTa Ha+β I)xopt = VT j⇔ xopt =F−1(F (V)∗F (j)F (V)∗F (V)+α∑aF (Ha)∗F (Ha)+β).(A.6)Note that minimizing the quadratic from the x-step is equivalent to solving thelinear equation system from the second last row in Eq. A.6. The linear systemis composed of the matrices V, H and β I which all are very structured, that isall of these matrices are convolution matrices. This structure can be exploited byreformulating this linear equation system in the frequency domain in the last rowof Eq. A.6. This can then be inverted very efficiently by point-wise division.On the other hand, the v-step of our blind deconvolution requires the solutionof a deconvolution problem with known kernel Xk. It involves a quadratic dataterm, sparse cross channel correlation term and sparse gradient term. Due to the`1-norm penalty of the last two terms solving this minimization problem does notreduce to a quadratic problem as for the x-step in Eq. A.6.Note that (given a kernel estimate) solving the non-blind deconvolution prob-lem is in fact part of our joint blind reconstruction. This is in strong contrast torecent state-of-the-art blind deconvolution methods such as [149].A.2.2 Non-blind deconvolution.The previous subsection has described the overall reconstruction Alg. 2 of Chap-ter 6 that jointly self-calibrates the PSF and estimates the unknown latent image.152This subsection explains how we solve the non-blind deconvolution sub-problem(v-step), that is solving for the latent image given a fixed PSF estimate Xk andprior knowledge of spectral gradient correlation. The cross-channel prior fromSec. A.2.1 allows for formulating an efficient inversion method. To show this, itmakes sense to first revisit the formulation proposed by [57], that we also brieflydiscuss in Chapter 3:vopt = argminv‖Xv− j‖22+µ∑a‖Hav‖1 + γ∑a‖Hav · ir−Hair ·v‖1. (A.7)Note that slack variables are introduced for all the terms Hax and all the Hav ·ir−Hair ·v terms. Two slack variables are necessary since statistics on normalizedgradients are considered (multiplication with ir and v). The resulting optimizationproblem is then solved via Chambolle and Pock’s primal-dual framework [18].Introducing many slack terms results in slow convergence in order to achieve con-sensus between all the consensus constraints and the objective terms.Efficient optimization. Our formulation of the cross-channel prior leads to thefollowing sub-problem (v-step of Alg. 2):vopt = argminv‖Xv− j‖22+µ∑a‖Hav‖1 + γ∑a‖Hav−Hair‖1. (A.8)This formulation only requires a single slack variables for Hav. That is thesame number of slack variables usually needed for `1-deconvolution methods. Tosee this, let us formulate the proximal operators that would be necessary if we wereto apply [18] to Eq. A.8. Two proximal operators, one for the sparse gradient termand the other for the cross-channel term, can be formulated using only p = Hav as153variable:proxθ‖·‖1(p) = max(1− θβ|p| ,0) p Shrinkageproxθ‖ · −α‖1(p) = argminxθ‖x−α‖1+ 12‖x−p‖22= α+ argminzθ‖z‖1+ 12‖z+α−p‖22 with z := x−α= α+ argminzθ‖z‖1+ 12‖z− (p−α)‖22= proxθ‖·‖1(p−α)+α= max(1− θβ|p−α| ,0) p+α Cross-shrinkage(A.9)where here α = Hair. To derive the second proximal operator, we have used asimple substitution trick via substitution of z := x−α . The proximal operator thenreduces to the shrinkage operator. Using this trick, we have defined the proximaloperators for both, gradient and cross channel term, using the same slack variablefor p = Hax. It is easy to see, that for Eq. (A.7) the multiplicative term in thecross-channel term prohibits this approach and thus another slack variable needs tobe introduced.Having defined the objective in a way that it possesses an efficient splitting,we can now choose between different proximal algorithms to minimize it. Sub-sequently, we test two choices, the well known Alternating Direction Method ofMultiplier (ADMM) and the fast approximation using fast Half Quadratic Splitting(HQS) (recall Sec. 2.3 of Chapter 2 for details).154
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Computational imaging with diffractive optics
Text