TEXTS IN COMPUTER SCIENCE

Computer Vision

Algorithms and Applications Second Edition

Richard Szeliski

Texts in Computer Science

Series Editors

David Gries, Department of Computer Science, Cornell University, Ithaca, NY, USA

Orit Hazzan $ ^{ID} $, Faculty of Education in Technology and Science, Technion—Israel Institute of Technology, Haifa, Israel

Titles in this series now included in the Thomson Reuters Book Citation Index!

'Texts in Computer Science' (TCS) delivers high-quality instructional content for undergraduates and graduates in all areas of computing and information science, with a strong emphasis on core foundational and theoretical material but inclusive of some prominent applications-related content. TCS books should be reasonably self-contained and aim to provide students with modern and clear accounts of topics ranging across the computing curriculum. As a result, the books are ideal for semester courses or for individual self-study in cases where people need to expand their knowledge. All texts are authored by established experts in their fields, reviewed internally and by the series editors, and provide numerous examples, problems, and other pedagogical tools; many contain fully worked solutions.

The TCS series is comprised of high-quality, self-contained books that have broad and comprehensive coverage and are generally in hardback format and sometimes contain color. For undergraduate textbooks that are likely to be more brief and modular in their approach, require only black and white, and are under 275 pages, Springer offers the flexibly designed Undergraduate Topics in Computer Science series, to which we refer potential authors.

More information about this series at https://link.springer.com/bookseries/3191

Richard Szeliski

Computer Vision

Algorithms and Applications

Second Edition

Richard Szeliski

University of Washington

Seattle, WA, USA

ISSN 1868-0941

ISSN 1868-095X (electronic)

Texts in Computer Science

ISBN 978-3-030-34371-2

ISBN 978-3-030-34372-9 (eBook)

https://doi.org/10.1007/978-3-030-34372-9

1 $ ^{st} $ edition: © Springer-Verlag London Limited 2011

2 $ ^{nd} $ edition: © Springer Nature Switzerland AG 2022

This work is subject to copyright. All rights are reserved by the Publisher, whether the whole or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, reproduction on microfilms or in any other physical way, and transmission or information storage and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or hereafter developed.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Instance segmentation (delineating individual objects in a scene with pixel-accurate masks) using the Mask R-CNN system of He, Gkioxari et al. (2017), © 2017 IEEE (reproduced with permission).

This Springer imprint is published by the registered company Springer Nature Switzerland AG. The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland.

This book is dedicated to my parents,

Zdzistaw and Jadwiga,

and my family,

Lyn, Anne, and Stephen.

1 Introduction

What is computer vision? • A brief history

Book overview • Sample syllabus • Notation

2 Image formation

Geometric primitives and transformations

Photometric image formation • The digital camera

3 Image processing

Point operators • Linear filtering

Non-linear filtering • Fourier transforms

Pyramids and wavelets • Geometric transformations

4 Model fitting and optimization

Scattered data interpolation

Variational methods and regularization

Markov random fields

5 Deep learning

Supervised learning • Unsupervised learning

Deep neural networks • Convolutional networks

More complex models

6 Recognition

Instance recognition • Image classification

Object detection • Semantic segmentation

Video understanding • Vision and language

7 Feature detection and matching

Points and patches • Edges and contours

Contour tracking • Lines and vanishing points

Segmentation

8 Image alignment and stitching 401

Pairwise alignment • Image stitching • Global alignment • Compositing

9 Motion estimation 443

Translational alignment • Parametric motion • Optical flow • Layered motion

10 Computational photography 483

Photometric calibration • High dynamic range imaging • Super-resolution, denoising, and blur removal • Image matting and compositing • Texture analysis and synthesis

11 Structure from motion and SLAM 543

Geometric intrinsic calibration • Pose estimation • Two-frame structure from motion • Multi-frame structure from motion • Simultaneous localization and mapping (SLAM)

12 Depth estimation 595

Epipolar geometry • Sparse correspondence • Dense correspondence • Local methods • Global optimization • Deep neural networks • Multi-view stereo • Monocular depth estimation

13 3D reconstruction 639

Shape from X • 3D scanning • Surface representations • Point-based representations • Volumetric representations • Model-based reconstruction • Recovering texture maps and albedos

14 Image-based rendering 681

View interpolation • Layered depth images • Light fields and Lumigraphs • Environment matting • Video-based rendering • Neural rendering

Preface

The seeds for this book were first planted in 2001 when Steve Seitz at the University of Washington invited me to co-teach a course called “Computer Vision for Computer Graphics”. At that time, computer vision techniques were increasingly being used in computer graphics to create image-based models of real-world objects, to create visual effects, and to merge real-world imagery using computational photography techniques. Our decision to focus on the applications of computer vision to fun problems such as image stitching and photo-based 3D modeling from personal photos seemed to resonate well with our students.

That initial course evolved into a more complete computer vision syllabus and project-oriented course structure that I used to co-teach general computer vision courses both at the University of Washington and at Stanford. (The latter was a course I co-taught with David Fleet in 2003.) Similar curricula were then adopted at a number of other universities and also incorporated into more specialized courses on computational photography. (For ideas on how to use this book in your own course, please see Table 1.1 in Section 1.4.)

This book also reflects my 40 years' experience doing computer vision research in corporate research labs, mostly at Digital Equipment Corporation's Cambridge Research Lab, Microsoft Research, and Facebook. In pursuing my work, I have mostly focused on problems and solution techniques (algorithms) that have practical real-world applications and that work well in practice. Thus, this book has more emphasis on basic techniques that work under real-world conditions and less on more esoteric mathematics that has intrinsic elegance but less practical applicability.

This book is suitable for teaching a senior-level undergraduate course in computer vision to students in both computer science and electrical engineering. I prefer students to have either an image processing or a computer graphics course as a prerequisite, so that they can spend less time learning general background mathematics and more time studying computer vision techniques. The book is also suitable for teaching graduate-level courses in computer vision, e.g., by delving into more specialized topics, and as a general reference to fundamental techniques and the recent research literature. To this end, I have attempted wherever possible to at least cite the newest research in each sub-field, even if the technical details are too complex to cover in the book itself.

In teaching our courses, we have found it useful for the students to attempt a number of small implementation projects, which often build on one another, in order to get them used to working with real-world images and the challenges that these present. The students are then asked to choose an individual topic for each of their small-group, final projects. (Sometimes these projects even turn into conference papers!) The exercises at the end of each chapter contain numerous suggestions for smaller mid-term projects, as well as more open-ended problems whose solutions are still active research topics. Wherever possible, I encourage students to try their algorithms on their own personal photographs, since this better motivates them, often leads to creative variants on the problems, and better acquaints them with the variety and complexity of real-world imagery.

In formulating and solving computer vision problems, I have often found it useful to draw inspi-

ration from four high-level approaches:

• Scientific: build detailed models of the image formation process and develop mathematical techniques to invert these in order to recover the quantities of interest (where necessary, making simplifying assumptions to make the mathematics more tractable).

• Statistical: use probabilistic models to quantify the prior likelihood of your unknowns and the noisy measurement processes that produce the input images, then infer the best possible estimates of your desired quantities and analyze their resulting uncertainties. The inference algorithms used are often closely related to the optimization techniques used to invert the (scientific) image formation processes.

• Engineering: develop techniques that are simple to describe and implement but that are also known to work well in practice. Test these techniques to understand their limitation and failure modes, as well as their expected computational costs (run-time performance).

• Data-driven: collect a representative set of test data (ideally, with labels or ground-truth answers) and use these data to either tune or learn your model parameters, or at least to validate and quantify its performance.

These four approaches build on each other and are used throughout the book.

My personal research and development philosophy (and hence the exercises in the book) have a strong emphasis on testing algorithms. It's too easy in computer vision to develop an algorithm that does something plausible on a few images rather than something correct. The best way to validate your algorithms is to use a three-part strategy.

First, test your algorithm on clean synthetic data, for which the exact results are known. Second, add noise to the data and evaluate how the performance degrades as a function of noise level. Finally, test the algorithm on real-world data, preferably drawn from a wide variety of sources, such as photos found on the web. Only then can you truly know if your algorithm can deal with real-world complexity, i.e., images that do not fit some simplified model or assumptions.

in order to help students in this process, Appendix C includes pointers to commonly used datasets and software libraries that contain implementations of a wide variety of computer vision algorithms, which can enable you to tackle more ambitious projects (with your instructor's consent).

Notes on the Second Edition

The last decade has seen a truly dramatic explosion in the performance and applicability of computer vision algorithms, much of it engendered by the application of machine learning algorithms to large amounts of visual training data (Su and Crandall 2021).

Deep neural networks now play an essential role in so many vision algorithms that the new edition of this book introduces them early on as a fundamental technique that gets used extensively in subsequent chapters.

The most notable changes in the second edition include:

• Machine learning, deep learning, and deep neural networks are introduced early on in Chapter 5, as they play just as fundamental a role in vision algorithms as more classical techniques, such as image processing, graphical/probabilistic models, and energy minimization, which are introduced in the preceding two chapters.

• The recognition chapter has been moved earlier in the book to Chapter 6, since end-to-end deep learning systems no longer require the development of building blocks such as feature

detection, matching, and segmentation. Many of the students taking vision classes are primarily interested in visual recognition, so presenting this material earlier in the course makes it easier for students to base their final project on these topics. This chapter also includes sections on semantic segmentation, video understanding, and vision and language.

• The application of neural networks and deep learning to myriad computer vision algorithms and applications, including flow and stereo, 3D shape modeling, and newly emerging fields such as neural rendering.

• New technologies such as SLAM (simultaneous localization and mapping) and VIO (visual inertial odometry) that now run reliably and are used in real-time applications such as augmented reality and autonomous navigation.

In addition to these larger changes, the book has been updated to reflect the latest state-of-the-art techniques such as internet-scale image search and phone-based computational photography. The new edition includes over 1500 new citations (papers) and has over 200 new figures.

Acknowledgements

I would like to gratefully acknowledge all of the people whose passion for research and inquiry as well as encouragement have helped me write this book.

Steve Zucker at McGill University first introduced me to computer vision, taught all of his students to question and debate research results and techniques, and encouraged me to pursue a graduate career in this area.

Takeo Kanade and Geoff Hinton, my PhD thesis advisors at Carnegie Mellon University, taught me the fundamentals of good research, writing, and presentation and mentored several generations of outstanding students and researchers. They fired up my interest in visual processing, 3D modeling, and statistical methods, while Larry Matthies introduced me to Kalman filtering and stereo matching. Geoff continues to inspire so many of us with this undiminished passion for trying to figure out “what makes the brain work”. It’s been a delight to see his pursuit of connectionist ideas bear so much fruit in this past decade.

Demetri Terzopoulos was my mentor at my first industrial research job and taught me the ropes of successful publishing. Yvan Leclerc and Pascal Fua, colleagues from my brief interlude at SRI International, gave me new perspectives on alternative approaches to computer vision.

During my six years of research at Digital Equipment Corporation's Cambridge Research Lab, I was fortunate to work with a great set of colleagues, including Ingrid Carlbom, Gudrun Klinker, Keith Waters, William Hsu, Richard Weiss, Stéphane Lavallée, and Sing Bing Kang, as well as to supervise the first of a long string of outstanding summer interns, including David Tonnesen, Sing Bing Kang, James Coughlan, and Harry Shum. This is also where I began my long-term collaboration with Daniel Scharstein.

At Microsoft Research, I had the outstanding fortune to work with some of the world's best researchers in computer vision and computer graphics, including Michael Cohen, Matt Uyttendale, Sing Bing Kang, Harry Shum, Larry Zitnick, Sudipta Sinha, Drew Steedly, Simon Baker, Johannes Kopf, Neel Joshi, Krishnan Ramnath, Anandan, Phil Torr, Antonio Criminisi, Simon Winder, Matthew Brown, Michael Goesele, Richard Hartley, Hugues Hoppe, Stephen Gortler, Steve Shafer, Matthew Turk, Georg Petschnigg, Kentaro Toyama, Ramin Zabih, Shai Avidan, Patrice Simard, Chris Pal, Nebojsa Jojic, Patrick Baudisch, Dani Lischinski, Raanan Fattal, Eric Stollnitz, David Nistér, Blaise Aguera y Arcas, Andrew Fitzgibbon, Jamie Shotton, Wolf Kienzle, Piotr Dollar, and Ross Girshick. I was also lucky to have as interns such great students as Polina Golland, Simon

Baker, Mei Han, Arno Schödl, Ron Dror, Ashley Eden, Jonathan Shade, Jinxiang Chai, Rahul Swaminathan, Yanghai Tsin, Sam Hasinoff, Anat Levin, Matthew Brown, Eric Bennett, Vaibhav Vaish, Jan-Michael Frahm, James Diebel, Ce Liu, Josef Sivic, Grant Schindler, Colin Zheng, Neel Joshi, Sudipta Sinha, Zeev Farbman, Rahul Garg, Tim Cho, Yekeun Jeong, Richard Roberts, Varsha Hedau, Dilip Krishnan, Adarsh Kowdle, Edward Hsiao, Yong Seok Heo, Fabian Langguth, Andrew Owens, and Tianfan Xue. Working with such outstanding students also gave me the opportunity to collaborate with some of their amazing advisors, including Bill Freeman, Irfan Essa, Marc Pollefeys, Michael Black, Marc Levoy, and Andrew Zisserman.

Since moving to Facebook, I've had the pleasure to continue my collaborations with Michael Cohen, Matt Uyttendaele, Johannes Kopf, Wolf Kienzle, and Krishnan Ramnath, and also new colleagues including Kevin Matzen, Bryce Evans, Suhib Alisian, Changil Kim, David Geraghty, Jan Herling, Nils Plath, Jan-Michael Frahm, True Price, Richard Newcombe, Thomas Whelan, Michael Goesele, Steven Lovegrove, Julian Straub, Simon Green, Brian Cabral, Michael Toksvig, Albert Para Pozzo, Laura Sevilla-Lara, Georgia Gkioxari, Justin Johnson, Chris Sweeney, and Vassileios Balntas. I've also had the pleasure to collaborate with some outstanding summer interns, including Tianfan Xue, Scott Wehrwein, Peter Hedman, Joel Janai, Aleksander Holynski, Xuan Luo, Rui Wang, Olivia Wiles, and Yulun Tian. I'd like to thank in particular Michael Cohen, my mentor, colleague, and friend for the last 25 years for his unwavering support of my sprint to complete this second edition.

While working at Microsoft and Facebook, I’ve also had the opportunity to collaborate with wonderful colleagues at the University of Washington, where I hold an Affiliate Professor appointment. I’m indebted to Tony DeRose and David Salesin, who first encouraged me to get involved with the research going on at UW, my long-time collaborators Brian Curless, Steve Seitz, Maneesh Agrawala, Sameer Agarwal, and Yasu Furukawa, as well as the students I have had the privilege to supervise and interact with, including Frédéric Pighin, Yung-Yu Chuang, Doug Zongker, Colin Zheng, Aseem Agarwala, Dan Goldman, Noah Snavely, Ian Simon, Rahul Garg, Ryan Kaminsky, Juliet Fiss, Aleksander Hołyński, and Yifan Wang. As I mentioned at the beginning of this preface, this book owes its inception to the vision course that Steve Seitz invited me to co-teach, as well as to Steve’s encouragement, course notes, and editorial input.

I'm also grateful to the many other computer vision researchers who have given me so many constructive suggestions about the book, including Sing Bing Kang, who was my informal book editor, Vladimir Kolmogorov, Daniel Scharstein, Richard Hartley, Simon Baker, Noah Snavely, Bill Freeman, Svetlana Lazebnik, Matthew Turk, Jitendra Malik, Alyosha Efros, Michael Black, Brian Curless, Sameer Agarwal, Li Zhang, Deva Ramanan, Olga Veksler, Yuri Boykov, Carsten Rother, Phil Torr, Bill Triggs, Bruce Maxwell, Rico Malvar, Jana Košeká, Eero Simoncelli, Aaron Hertzmann, Antonio Torralba, Tomaso Poggio, Theo Pavlidis, Baba Vemuri, Nando de Freitas, Chuck Dyer, Song Yi, Falk Schubert, Roman Pflugfelder, Marshall Tappen, James Coughlan, Sammy Rogmans, Klaus Strobel, Shanmuganathan, Andreas Siebert, Yongjun Wu, Fred Pighin, Juan Cockburn, Ronald Mallet, Tim Soper, Georgios Evangelidis, Dwight Fowler, Itzik Bayaz, Daniel O'Connor, Srikrishna Bhat, and Toru Tamaki, who wrote the Japanese translation and provided many useful errata.

For the second edition, I received significant help and advice from three key contributors. Daniel Scharstein helped me update the chapter on stereo, Matt Deitke contributed descriptions of the newest papers in deep learning, including the sections on transformers, variational autoencoders, and text-to-image synthesis, along with the exercises in Chapters 5 and 6 and some illustrations. Sing Bing Kang reviewed multiple drafts and provided useful suggestions. I'd also like to thank Andrew Glassner, whose book (Glassner 2018) and figures were a tremendous help, Justin Johnson,

Sean Bell, Ishan Misra, David Fouhey, Michael Brown, Abdelrahman Abdelhamed, Frank Dellaert, Xinlei Chen, Ross Girshick, Andreas Geiger, Dmytro Mishkin, Aleksander Hołyński, Joel Janai, Christoph Feichtenhofer, Yuandong Tian, Alyosha Efros, Pascal Fua, Torsten Sattler, Laura Leal-Taixé, Aljosa Osep, Qunjie Zhou, Jiří Matas, Eddy Ilg, Yann LeCun, Larry Jackel, Vasileios Balntas, Daniel DeTone, Zachary Teed, Junhwa Hur, Jun-Yan Zhu, Filip Radenović, Michael Zollhöfer, Matthias Nießner, Andrew Owens, Hervé Jégou, Luowei Zhou, Ricardo Martin Brualla, Pratul Srinivasan, Matteo Poggi, Fabio Tosi, Ahmed Osman, Dave Howell, Holger Heidrich, Howard Yen, Anton Papst, Syamprasad K. Rajagopalan, Abhishek Nagar, Vladimir Kuznetsov, Raphaël Fouque, Marian Ciobanu, Darko Simonovic, and Guilherme Schlinker.

In preparing the second edition, I taught some of the new material in two courses that I helped co-teach in 2020 at Facebook and UW. I'd like to thank my co-instructors Jan-Michael Frahm, Michael Goesele, Georgia Gkioxari, Ross Girshick, Jakob Julian Engel, Daniel Scharstein, Fernando de la Torre, Steve Seitz, and Harpreet Sawhney, from whom I learned a lot about the latest techniques that are included in the new edition. I'd also like to thank the TAs, including David Geraghty, True Price, Kevin Matzen, Akash Bapat, Aleksander Holyński, Keunhong Park, and Svetoslav Kolev, for the wonderful job they did in creating and grading the assignments. I'd like to give a special thanks to Justin Johnson, whose excellent class slides (Johnson 2020), based on earlier slides from Stanford (Li, Johnson, and Yeung 2019), taught me the fundamentals of deep learning and which I used extensively in my own class and in preparing the new chapter on deep learning.

Shena Deuchers and Ian Kingston did a fantastic job copy-editing the first and second editions, respectively and suggesting many useful improvements, and Wayne Wheeler and Simon Rees at Springer were most helpful throughout the whole book publishing process. Keith Price's Annotated Computer Vision Bibliography was invaluable in tracking down references and related work.

If you have any suggestions for improving the book, please send me an e-mail, as I would like to keep the book as accurate, informative, and timely as possible.

The last year of writing this second edition took place during the worldwide COVID-19 pandemic. I would like to thank all of the first responders, medical and front-line workers, and everyone else who helped get us through these difficult and challenging times and to acknowledge the impact that this and other recent tragedies have had on all of us.

Lastly, this book would not have been possible or worthwhile without the incredible support and encouragement of my family. I dedicate this book to my parents, Zdzisław and Jadwiga, whose love, generosity, and accomplishments always inspired me; to my sister Basia for her lifelong friendship; and especially to Lyn, Anne, and Stephen, whose love and support in all matters (including my book projects) makes it all worthwhile.

Lake Wenatchee

May 2021

Contents

Preface ix

1 Introduction 1

1.1 What is computer vision? ..... 3

1.2 A brief history ..... 9

1.3 Book overview ..... 18

1.4 Sample syllabus ..... 24

1.5 A note on notation ..... 25

1.6 Additional reading ..... 26

2 Image formation 27

2.1 Geometric primitives and transformations ..... 29

2.1.1 2D transformations ..... 32

2.1.2 3D transformations ..... 35

2.1.3 3D rotations ..... 36

2.1.4 3D to 2D projections ..... 41

2.1.5 Lens distortions ..... 51

2.2 Photometric image formation ..... 53

2.2.1 Lighting ..... 53

2.2.2 Reflectance and shading ..... 54

2.2.3 Optics ..... 59

2.3 The digital camera ..... 63

2.3.1 Sampling and aliasing ..... 67

2.3.2 Color ..... 69

2.3.3 Compression ..... 78

2.4 Additional reading ..... 80

2.5 Exercises ..... 80

3 Image processing 85

3.1 Point operators ..... 87

3.1.1 Pixel transforms ..... 89

3.1.2 Color transforms ..... 90

3.1.3 Compositing and matting ..... 91

3.1.4 Histogram equalization ..... 92

3.1.5 Application: Tonal adjustment ..... 95

3.2 Linear filtering ..... 95

3.2.1 Separable filtering ..... 99

3.2.2 Examples of linear filtering ..... 100

3.2.3 Band-pass and steerable filters ..... 101

3.3 More neighborhood operators ..... 105

3.3.1 Non-linear filtering ..... 105

3.3.2 Bilateral filtering ..... 107

3.3.3 Binary image processing ..... 110

3.4 Fourier transforms ..... 113

3.4.1 Two-dimensional Fourier transforms ..... 115

3.4.2 Application: Sharpening, blur, and noise removal ..... 118

3.5 Pyramids and wavelets ..... 119

3.5.1 Interpolation ..... 119

3.5.2 Decimation ..... 122

3.5.3 Multi-resolution representations ..... 123

3.5.4 Wavelets ..... 128

3.5.5 Application: Image blending ..... 132

3.6 Geometric transformations ..... 135

3.6.1 Parametric transformations ..... 135

3.6.2 Mesh-based warping ..... 140

3.6.3 Application: Feature-based morphing ..... 142

3.7 Additional reading ..... 143

3.8 Exercises ..... 144

4 Model fitting and optimization ..... 153

4.1 Scattered data interpolation ..... 155

4.1.1 Radial basis functions ..... 157

4.1.2 Overfitting and underfitting ..... 159

4.1.3 Robust data fitting ..... 162

4.2 Variational methods and regularization ..... 163

4.2.1 Discrete energy minimization ..... 166

4.2.2 Total variation ..... 168

4.2.3 Bilateral solver ..... 168

4.2.4 Application: Interactive colorization ..... 169

4.3 Markov random fields ..... 170

4.3.1 Conditional random fields ..... 177

4.3.2 Application: Interactive segmentation ..... 181

4.4 Additional reading ..... 184

4.5 Exercises ..... 185

5 Deep Learning ..... 187

5.1 Supervised learning ..... 191

5.1.1 Nearest neighbors ..... 192

5.1.2 Bayesian classification ..... 194

5.1.3 Logistic regression ..... 198

5.1.4 Support vector machines ..... 199

5.1.5 Decision trees and forests ..... 202

5.2 Unsupervised learning ..... 205

5.2.1 Clustering ..... 205

5.2.2 K-means and Gaussians mixture models ..... 206

Contents

5.2.3 Principal component analysis ..... 209

5.2.4 Manifold learning ..... 211

5.2.5 Semi-supervised learning ..... 212

5.3 Deep neural networks ..... 214

5.3.1 Weights and layers ..... 215

5.3.2 Activation functions ..... 217

5.3.3 Regularization and normalization ..... 219

5.3.4 Loss functions ..... 223

5.3.5 Backpropagation ..... 226

5.3.6 Training and optimization ..... 228

5.4 Convolutional neural networks ..... 231

5.4.1 Pooling and unpooling ..... 234

5.4.2 Application: Digit classification ..... 237

5.4.3 Network architectures ..... 238

5.4.4 Model zoos ..... 242

5.4.5 Visualizing weights and activations ..... 244

5.4.6 Adversarial examples ..... 248

5.4.7 Self-supervised learning ..... 249

5.5 More complex models ..... 252

5.5.1 Three-dimensional CNNs ..... 252

5.5.2 Recurrent neural networks ..... 255

5.5.3 Transformers ..... 257

5.5.4 Generative models ..... 261

5.6 Additional reading ..... 267

5.7 Exercises ..... 268

6 Recognition ..... 273

6.1 Instance recognition ..... 276

6.2 Image classification ..... 278

6.2.1 Feature-based methods ..... 278

6.2.2 Deep networks ..... 285

6.2.3 Application: Visual similarity search ..... 287

6.2.4 Face recognition ..... 289

6.3 Object detection ..... 295

6.3.1 Face detection ..... 295

6.3.2 Pedestrian detection ..... 299

6.3.3 General object detection ..... 301

6.4 Semantic segmentation ..... 307

6.4.1 Application: Medical image segmentation ..... 310

6.4.2 Instance segmentation ..... 311

6.4.3 Panoptic segmentation ..... 312

6.4.4 Application: Intelligent photo editing ..... 314

6.4.5 Pose estimation ..... 315

6.5 Video understanding ..... 316

6.6 Vision and language ..... 319

6.7 Additional reading ..... 326

6.8 Exercises ..... 329

7 Feature detection and matching 333

7.1 Points and patches 335

7.1.1 Feature detectors 337

7.1.2 Feature descriptors 347

7.1.3 Feature matching 352

7.1.4 Large-scale matching and retrieval 358

7.1.5 Feature tracking 361

7.1.6 Application: Performance-driven animation 363

7.2 Edges and contours 364

7.2.1 Edge detection 364

7.2.2 Contour detection 368

7.2.3 Application: Edge editing and enhancement 372

7.3 Contour tracking 373

7.3.1 Snakes and scissors 373

7.3.2 Level Sets 379

7.3.3 Application: Contour tracking and rotoscoping 380

7.4 Lines and vanishing points 381

7.4.1 Successive approximation 381

7.4.2 Hough transforms 381

7.4.3 Vanishing points 384

7.5 Segmentation 386

7.5.1 Graph-based segmentation 388

7.5.2 Mean shift 389

7.5.3 Normalized cuts 391

7.6 Additional reading 393

7.7 Exercises 395

8 Image alignment and stitching 401

8.1 Pairwise alignment 403

8.1.1 2D alignment using least squares 403

8.1.2 Application: Panography 405

8.1.3 Iterative algorithms 406

8.1.4 Robust least squares and RANSAC 408

8.1.5 3D alignment 410

8.2 Image stitching 411

8.2.1 Parametric motion models 412

8.2.2 Application: Whiteboard and document scanning 414

8.2.3 Rotational panoramas 414

8.2.4 Gap closing 416

8.2.5 Application: Video summarization and compression 417

8.2.6 Cylindrical and spherical coordinates 418

8.3 Global alignment 421

8.3.1 Bundle adjustment 421

8.3.2 Parallax removal 424

8.3.3 Recognizing panoramas 425

8.4 Compositing 426

8.4.1 Choosing a compositing surface 426

8.4.2 Pixel selection and weighting (deghosting) 430

Contents

8.4.3 Application: Photomontage ..... 435

8.4.4 Blending ..... 435

8.5 Additional reading ..... 437

8.6 Exercises ..... 438

9 Motion estimation ..... 443

9.1 Translational alignment ..... 445

9.1.1 Hierarchical motion estimation ..... 448

9.1.2 Fourier-based alignment ..... 449

9.1.3 Incremental refinement ..... 451

9.2 Parametric motion ..... 455

9.2.1 Application: Video stabilization ..... 457

9.2.2 Spline-based motion ..... 459

9.2.3 Application: Medical image registration ..... 461

9.3 Optical flow ..... 461

9.3.1 Deep learning approaches ..... 466

9.3.2 Application: Rolling shutter wobble removal ..... 468

9.3.3 Multi-frame motion estimation ..... 468

9.3.4 Application: Video denoising ..... 469

9.4 Layered motion ..... 470

9.4.1 Application: Frame interpolation ..... 473

9.4.2 Transparent layers and reflections ..... 474

9.4.3 Video object segmentation ..... 476

9.4.4 Video object tracking ..... 477

9.5 Additional reading ..... 478

9.6 Exercises ..... 479

10 Computational photography ..... 483

10.1 Photometric calibration ..... 486

10.1.1 Radiometric response function ..... 486

10.1.2 Noise level estimation ..... 488

10.1.3 Vignetting ..... 489

10.1.4 Optical blur (spatial response) estimation ..... 491

10.2 High dynamic range imaging ..... 494

10.2.1 Tone mapping ..... 499

10.2.2 Application: Flash photography ..... 506

10.3 Super-resolution, denoising, and blur removal ..... 508

10.3.1 Color image demosaicing ..... 515

10.3.2 Lens blur (bokeh) ..... 517

10.4 Image matting and compositing ..... 518

10.4.1 Blue screen matting ..... 519

10.4.2 Natural image matting ..... 521

10.4.3 Optimization-based matting ..... 524

10.4.4 Smoke, shadow, and flash matting ..... 527

10.4.5 Video matting ..... 528

10.5 Texture analysis and synthesis ..... 529

10.5.1 Application: Hole filling and inpainting ..... 531

10.5.2 Application: Non-photorealistic rendering ..... 532

10.5.3 Neural style transfer and semantic image synthesis ..... 534

10.6 Additional reading ..... 537

10.7 Exercises ..... 538

11 Structure from motion and SLAM ..... 543

11.1 Geometric intrinsic calibration ..... 545

11.1.1 Vanishing points ..... 547

11.1.2 Application: Single view metrology ..... 548

11.1.3 Rotational motion ..... 549

11.1.4 Radial distortion ..... 550

11.2 Pose estimation ..... 552

11.2.1 Linear algorithms ..... 552

11.2.2 Iterative non-linear algorithms ..... 554

11.2.3 Application: Location recognition ..... 555

11.2.4 Triangulation ..... 558

11.3 Two-frame structure from motion ..... 560

11.3.1 Eight, seven, and five-point algorithms ..... 560

11.3.2 Special motions and structures ..... 564

11.3.3 Projective (uncalibrated) reconstruction ..... 565

11.3.4 Self-calibration ..... 566

11.3.5 Application: View morphing ..... 568

11.4 Multi-frame structure from motion ..... 568

11.4.1 Factorization ..... 568

11.4.2 Bundle adjustment ..... 570

11.4.3 Exploiting sparsity ..... 571

11.4.4 Application: Match move ..... 574

11.4.5 Uncertainty and ambiguities ..... 575

11.4.6 Application: Reconstruction from internet photos ..... 576

11.4.7 Global structure from motion ..... 578

11.4.8 Constrained structure and motion ..... 580

11.5 Simultaneous localization and mapping (SLAM) ..... 583

11.5.1 Application: Autonomous navigation ..... 585

11.5.2 Application: Smartphone augmented reality ..... 587

11.6 Additional reading ..... 588

11.7 Exercises ..... 590

12 Depth estimation ..... 595

12.1 Epipolar geometry ..... 599

12.1.1 Rectification ..... 600

12.1.2 Plane sweep ..... 602

12.2 Sparse correspondence ..... 604

12.2.1 3D curves and profiles ..... 604

12.3 Dense correspondence ..... 606

12.3.1 Similarity measures ..... 607

12.4 Local methods ..... 609

12.4.1 Sub-pixel estimation and uncertainty ..... 610

12.4.2 Application: Stereo-based head tracking ..... 611

12.5 Global optimization ..... 612

Contents

12.5.1 Dynamic programming ..... 614

12.5.2 Segmentation-based techniques ..... 616

12.5.3 Application: Z-keying and background replacement ..... 617

12.6 Deep neural networks ..... 618

12.7 Multi-view stereo ..... 620

12.7.1 Scene flow ..... 624

12.7.2 Volumetric and 3D surface reconstruction ..... 624

12.7.3 Shape from silhouettes ..... 630

12.8 Monocular depth estimation ..... 632

12.9 Additional reading ..... 634

12.10 Exercises ..... 635

13 3D reconstruction ..... 639

13.1 Shape from X ..... 641

13.1.1 Shape from shading and photometric stereo ..... 642

13.1.2 Shape from texture ..... 645

13.1.3 Shape from focus ..... 646

13.2 3D scanning ..... 647

13.2.1 Range data merging ..... 650

13.2.2 Application: Digital heritage ..... 654

13.3 Surface representations ..... 654

13.3.1 Surface interpolation ..... 655

13.3.2 Surface simplification ..... 656

13.3.3 Geometry images ..... 656

13.4 Point-based representations ..... 657

13.5 Volumetric representations ..... 658

13.5.1 Implicit surfaces and level sets ..... 658

13.6 Model-based reconstruction ..... 660

13.6.1 Architecture ..... 660

13.6.2 Facial modeling and tracking ..... 663

13.6.3 Application: Facial animation ..... 665

13.6.4 Human body modeling and tracking ..... 668

13.7 Recovering texture maps and albedos ..... 674

13.7.1 Estimating BRDFs ..... 675

13.7.2 Application: 3D model capture ..... 676

13.8 Additional reading ..... 677

13.9 Exercises ..... 679

14 Image-based rendering ..... 681

14.1 View interpolation ..... 683

14.1.1 View-dependent texture maps ..... 685

14.1.2 Application: Photo Tourism ..... 686

14.2 Layered depth images ..... 688

14.2.1 Impostors, sprites, and layers ..... 688

14.2.2 Application: 3D photography ..... 690

14.3 Light fields and Lumigraphs ..... 693

14.3.1 Unstructured Lumigraph ..... 696

14.3.2 Surface light fields ..... 696

14.3.3 Application: Concentric mosaics ..... 698

14.3.4 Application: Synthetic re-focusing ..... 698

14.4 Environment mattes ..... 699

14.4.1 Higher-dimensional light fields ..... 700

14.4.2 The modeling to rendering continuum ..... 701

14.5 Video-based rendering ..... 701

14.5.1 Video-based animation ..... 702

14.5.2 Video textures ..... 703

14.5.3 Application: Animating pictures ..... 705

14.5.4 3D and free-viewpoint Video ..... 706

14.5.5 Application: Video-based walkthroughs ..... 708

14.6 Neural rendering ..... 711

14.7 Additional reading ..... 718

14.8 Exercises ..... 719

15 Conclusion ..... 723

A Linear algebra and numerical techniques ..... 727

A.1 Matrix decompositions ..... 728

A.1.1 Singular value decomposition ..... 728

A.1.2 Eigenvalue decomposition ..... 729

A.1.3 QR factorization ..... 731

A.1.4 Cholesky factorization ..... 732

A.2 Linear least squares ..... 733

A.2.1 Total least squares ..... 734

A.3 Non-linear least squares ..... 736

A.4 Direct sparse matrix techniques ..... 737

A.4.1 Variable reordering ..... 737

A.5 Iterative techniques ..... 738

A.5.1 Conjugate gradient ..... 739

A.5.2 Preconditioning ..... 740

A.5.3 Multigrid ..... 741

B Bayesian modeling and inference ..... 743

B.1 Estimation theory ..... 745

B.2 Maximum likelihood estimation and least squares ..... 747

B.3 Robust statistics ..... 748

B.4 Prior models and Bayesian inference ..... 750

B.5 Markov random fields ..... 751

B.6 Uncertainty estimation (error analysis) ..... 752

C Supplementary material ..... 755

C.1 Datasets and benchmarks ..... 756

C.2 Software ..... 761

C.3 Slides and lectures ..... 768

References ..... 769

Index ..... 905

Introduction

1.1 What is computer vision? ..... 3

1.2 A brief history ..... 9

1.3 Book overview ..... 18

1.4 Sample syllabus ..... 24

1.5 A note on notation ..... 25

1.6 Additional reading ..... 26

1.1 What is computer vision?

As humans, we perceive the three-dimensional structure of the world around us with apparent ease. Think of how vivid the three-dimensional percept is when you look at a vase of flowers sitting on the table next to you. You can tell the shape and translucency of each petal through the subtle patterns of light and shading that play across its surface and effortlessly segment each flower from the background of the scene (Figure 1.1). Looking at a framed group portrait, you can easily count and name all of the people in the picture and even guess at their emotions from their facial expressions (Figure 1.2a). Perceptual psychologists have spent decades trying to understand how the visual system works and, even though they can devise optical illusions $ ^{1} $ to tease apart some of its principles (Figure 1.3), a complete solution to this puzzle remains elusive (Marr 1982; Wandell 1995; Palmer 1999; Livingstone 2008; Frisby and Stone 2010).

Researchers in computer vision have been developing, in parallel, mathematical techniques for recovering the three-dimensional shape and appearance of objects in imagery. Here, the progress in the last two decades has been rapid. We now have reliable techniques for accurately computing a 3D model of an environment from thousands of partially overlapping photographs (Figure 1.2c). Given a large enough set of views of a particular object or façade, we can create accurate dense 3D surface models using stereo matching (Figure 1.2d). We can even, with moderate success, delineate most of the people and objects in a photograph (Figure 1.2a). However, despite all of these advances, the dream of having a computer explain an image at the same level of detail and causality as a two-year old remains elusive.

Why is vision so difficult? In part, it is because it is an inverse problem, in which we seek to recover some unknowns given insufficient information to fully specify the solution. We must therefore resort to physics-based and probabilistic models, or machine learning from large sets of examples, to disambiguate between potential solutions. However, modeling the visual world in all of its rich complexity is far more difficult than, say, modeling the vocal tract that produces spoken sounds.

The forward models that we use in computer vision are usually developed in physics (radiometry, optics, and sensor design) and in computer graphics. Both of these fields model how objects move and animate, how light reflects off their surfaces, is scattered by the atmosphere, refracted through camera lenses (or human eyes), and finally projected onto a flat (or curved) image plane. While computer graphics are not yet perfect, in many domains, such as rendering a still scene composed of everyday objects or animating extinct creatures such as dinosaurs, the illusion of reality is essentially there.

In computer vision, we are trying to do the inverse, i.e., to describe the world that we see in one or more images and to reconstruct its properties, such as shape, illumination, and color distributions. It is amazing that humans and animals do this so effortlessly, while computer vision algorithms are so error prone. People who have not worked in the field often underestimate the difficulty of the problem. This misperception that vision should be easy dates back to the early days of artificial intelligence (see Section 1.2), when it was initially believed that the cognitive (logic proving and planning) parts of intelligence were intrinsically more difficult than the perceptual components (Boden 2006).

The good news is that computer vision is being used today in a wide variety of real-world applications, which include:

• Optical character recognition (OCR): reading handwritten postal codes on letters (Fig.

1.1 What is computer vision?

ure 1.4a) and automatic number plate recognition (ANPR);

• Machine inspection: rapid parts inspection for quality assurance using stereo vision with specialized illumination to measure tolerances on aircraft wings or auto body parts (Figure 1.4b) or looking for defects in steel castings using X-ray vision;

• Retail: object recognition for automated checkout lanes and fully automated stores (Wingfield 2019);

• Warehouse logistics: autonomous package delivery and pallet-carrying “drives” (Guizzo 2008; O'Brian 2019) and parts picking by robotic manipulators (Figure 1.4c; Ackerman 2020);

• Medical imaging: registering pre-operative and intra-operative imagery (Figure 1.4d) or performing long-term studies of people's brain morphology as they age;

• Self-driving vehicles: capable of driving point-to-point between cities (Figure 1.4e; Monte-merlo, Becker et al. 2008; Urmson, Anhalt et al. 2008; Janai, Güney et al. 2020) as well as autonomous flight (Kaufmann, Gehrig et al. 2019);

• 3D model building (photogrammetry): fully automated construction of 3D models from aerial and drone photographs (Figure 1.4f);

Match move: merging computer-generated imagery (CGI) with live action footage by tracking feature points in the source video to estimate the 3D camera motion and shape of the environment. Such techniques are widely used in Hollywood, e.g., in movies such as Jurassic Park (Roble 1999; Roble and Zafar 2009); they also require the use of precise matting to insert new elements between foreground and background elements (Chuang, Agarwala et al. 2002).

• Motion capture (mocap): using retro-reflective markers viewed from multiple cameras or other vision-based techniques to capture actors for computer animation;

• Surveillance: monitoring for intruders, analyzing highway traffic and monitoring pools for drowning victims (e.g., https://swimeye.com);

• Fingerprint recognition and biometrics: for automatic access authentication as well as forensic applications.

David Lowe's website of industrial vision applications (http://www.cs.ubc.ca/spider/lowe/vision.html) lists many other interesting industrial applications of computer vision. While the above applications are all extremely important, they mostly pertain to fairly specialized kinds of imagery and narrow domains.

In addition to all of these industrial applications, there exist myriad consumer-level applications, such as things you can do with your own personal photographs and videos. These include:

• Stitching: turning overlapping photos into a single seamlessly stitched panorama (Figure 1.5a), as described in Section 8.2;

• Exposure bracketing: merging multiple exposures taken under challenging lighting conditions (strong sunlight and shadows) into a single perfectly exposed image (Figure 1.5b), as described in Section 10.2;

• Morphing: turning a picture of one of your friends into another, using a seamless morph transition (Figure 1.5c);

1.1 What is computer vision?

• 3D modeling: converting one or more snapshots into a 3D model of the object or person you are photographing (Figure 1.5d), as described in Section 13.6;

• Video match move and stabilization: inserting 2D pictures or 3D models into your videos by automatically tracking nearby reference points (see Section 11.4.4) $ ^{2} $ or using motion estimates to remove shake from your videos (see Section 9.2.1);

• Photo-based walkthroughs: navigating a large collection of photographs, such as the interior of your house, by flying between different photos in 3D (see Sections 14.1.2 and 14.5.5);

• Face detection: for improved camera focusing as well as more relevant image searching (see Section 6.3.1);

• Visual authentication: automatically logging family members onto your home computer as they sit down in front of the webcam (see Section 6.2.4).

The great thing about these applications is that they are already familiar to most students; they are, at least, technologies that students can immediately appreciate and use with their own personal media. Since computer vision is a challenging topic, given the wide range of mathematics being covered $ ^{3} $ and the intrinsically difficult nature of the problems being solved, having fun and relevant problems to work on can be highly motivating and inspiring.

The other major reason why this book has a strong focus on applications is that they can be used to formulate and constrain the potentially open-ended problems endemic in vision. Thus, it is better to think back from the problem at hand to suitable techniques, rather than to grab the first technique that you may have heard of. This kind of working back from problems to solutions is typical of an engineering approach to the study of vision and reflects my own background in the field.

First, I come up with a detailed problem definition and decide on the constraints and specifications for the problem. Then, I try to find out which techniques are known to work, implement a few of these, evaluate their performance, and finally make a selection. In order for this process to work, it is important to have realistic test data, both synthetic, which can be used to verify correctness and analyze noise sensitivity, and real-world data typical of the way the system will finally be used. If machine learning is being used, it is even more important to have representative unbiased training data in sufficient quantity to obtain good results on real-world inputs.

However, this book is not just an engineering text (a source of recipes). It also takes a scientific approach to basic vision problems. Here, I try to come up with the best possible models of the physics of the system at hand: how the scene is created, how light interacts with the scene and atmospheric effects, and how the sensors work, including sources of noise and uncertainty. The task is then to try to invert the acquisition process to come up with the best possible description of the scene.

The book often uses a statistical approach to formulating and solving computer vision problems. Where appropriate, probability distributions are used to model the scene and the noisy image acquisition process. The association of prior distributions with unknowns is often called Bayesian modeling (Appendix B). It is possible to associate a risk or loss function with misestimating the answer (Section B.2) and to set up your inference algorithm to minimize the expected risk. (Consider a robot trying to estimate the distance to an obstacle: it is usually safer to underestimate than to overestimate.) With statistical techniques, it often helps to gather lots of training data from which to

Navigation

Fun photos

learn probabilistic models. Finally, statistical approaches enable you to use proven inference techniques to estimate the best answer (or distribution of answers) and to quantify the uncertainty in the resulting estimates.

Because so much of computer vision involves the solution of inverse problems or the estimation of unknown quantities, my book also has a heavy emphasis on algorithms, especially those that are known to work well in practice. For many vision problems, it is all too easy to come up with a mathematical description of the problem that either does not match realistic real-world conditions or does not lend itself to the stable estimation of the unknowns. What we need are algorithms that are both robust to noise and deviation from our models and reasonably efficient in terms of run-time resources and space. In this book, I go into these issues in detail, using Bayesian techniques, where applicable, to ensure robustness, and efficient search, minimization, and linear system solving algorithms to ensure efficiency. $ ^{4} $ Most of the algorithms described in this book are at a high level, being mostly a list of steps that have to be filled in by students or by reading more detailed descriptions elsewhere. In fact, many of the algorithms are sketched out in the exercises.

Now that I’ve described the goals of this book and the frameworks that I use, I devote the rest of this chapter to two additional topics. Section 1.2 is a brief synopsis of the history of computer vision. It can easily be skipped by those who want to get to “the meat” of the new material in this book and do not care as much about who invented what when.

The second is an overview of the book's contents, Section 1.3, which is useful reading for everyone who intends to make a study of this topic (or to jump in partway, since it describes chapter interdependencies). This outline is also useful for instructors looking to structure one or more courses around this topic, as it provides sample curricula based on the book's contents.

1.2 A brief history

In this section, I provide a brief personal synopsis of the main developments in computer vision over the last fifty years (Figure 1.6) with a focus on advances I find personally interesting and that have stood the test of time. Readers not interested in the provenance of various ideas and the evolution of this field should skip ahead to the book overview in Section 1.3.

1970s. When computer vision first started out in the early 1970s, it was viewed as the visual perception component of an ambitious agenda to mimic human intelligence and to endow robots with intelligent behavior. At the time, it was believed by some of the early pioneers of artificial intelligence and robotics (at places such as MIT, Stanford, and CMU) that solving the “visual input” problem would be an easy step along the path to solving more difficult problems such as higher-level reasoning and planning. According to one well-known story, in 1966, Marvin Minsky at MIT asked his undergraduate student Gerald Jay Sussman to “spend the summer linking a camera to a computer and getting the computer to describe what it saw” (Boden 2006, p. 781). $ ^{5} $ We now know that the problem is slightly more difficult than that. $ ^{6} $

What distinguished computer vision from the already existing field of digital image processing (Rosenfeld and Pfaltz 1966; Rosenfeld and Kak 1976) was a desire to recover the three-dimensional

structure of the world from images and to use this as a stepping stone towards full scene understanding. Winston (1975) and Hanson and Riseman (1978) provide two nice collections of classic papers from this early period.

Early attempts at scene understanding involved extracting edges and then inferring the 3D structure of an object or a “blocks world” from the topological structure of the 2D lines (Roberts 1965). Several line labeling algorithms (Figure 1.7a) were developed at that time (Huffman 1971; Clowes 1971; Waltz 1975; Rosenfeld, Hummel, and Zucker 1976; Kanade 1980). Nalwa (1993) gives a nice review of this area. The topic of edge detection was also an active area of research; a nice survey of contemporaneous work can be found in (Davis 1975).

Three-dimensional modeling of non-polyhedral objects was also being studied (Baumgart 1974; Baker 1977). One popular approach used generalized cylinders, i.e., solids of revolution and swept closed curves (Agin and Binford 1976; Nevatia and Binford 1977), often arranged into parts relationships $ ^{7} $ (Hinton 1977; Marr 1982) (Figure 1.7c). Fischler and Elschlager (1973) called such elastic arrangements of parts pictorial structures (Figure 1.7b).

A qualitative approach to understanding intensities and shading variations and explaining them by the effects of image formation phenomena, such as surface orientation and shadows, was championed by Barrow and Tenenbaum (1981) in their paper on intrinsic images (Figure 1.7d), along with the related $ 2\frac{1}{2} $-D sketch ideas of Marr (1982). This approach has seen periodic revivals, e.g., in the work of Tappen, Freeman, and Adelson (2005) and Barron and Malik (2012).

More quantitative approaches to computer vision were also developed at the time, including the first of many feature-based stereo correspondence algorithms (Figure 1.7e) (Dev 1974; Marr and Poggio 1976, 1979; Barnard and Fischler 1982; Ohta and Kanade 1985; Grimson 1985; Pollard, Mayhew, and Frisby 1985) and intensity-based optical flow algorithms (Figure 1.7f) (Horn and Schunck 1981; Huang 1981; Lucas and Kanade 1981; Nagel 1986). The early work in simultaneously recovering 3D structure and camera motion (see Chapter 11) also began around this time (Ullman 1979; Longuet-Higgins 1981).

1.2 A brief history

A lot of the philosophy of how vision was believed to work at the time is summarized in David Marr’s (1982) book.⁸ In particular, Marr introduced his notion of the three levels of description of a (visual) information processing system. These three levels, very loosely paraphrased according to my own interpretation, are:

• Computational theory: What is the goal of the computation (task) and what are the constraints that are known or can be brought to bear on the problem?

• Representations and algorithms: How are the input, output, and intermediate information represented and which algorithms are used to calculate the desired result?

• Hardware implementation: How are the representations and algorithms mapped onto actual hardware, e.g., a biological vision system or a specialized piece of silicon? Conversely, how can hardware constraints be used to guide the choice of representation and algorithm? With the prevalent use of graphics chips (GPUs) and many-core architectures for computer vision, this question is again quite relevant.

As I mentioned earlier in this introduction, it is my conviction that a careful analysis of the problem specification and known constraints from image formation and priors (the scientific and statistical approaches) must be married with efficient and robust algorithms (the engineering approach) to design successful vision algorithms. Thus, it seems that Marr's philosophy is as good a guide to framing and solving problems in our field today as it was 25 years ago.

1980s. In the 1980s, a lot of attention was focused on more sophisticated mathematical techniques for performing quantitative image and scene analysis.

Image pyramids (see Section 3.5) started being widely used to perform tasks such as image blending (Figure 1.8a) and coarse-to-fine correspondence search (Rosenfeld 1980; Burt and Adelson 1983b; Rosenfeld 1984; Quam 1984; Anandan 1989). Continuous versions of pyramids using the concept of scale-space processing were also developed (Witkin 1983; Witkin, Terzopoulos, and Kass 1986; Lindeberg 1990). In the late 1980s, wavelets (see Section 3.5.4) started displacing or augmenting regular image pyramids in some applications (Mallat 1989; Simoncelli and Adelson 1990a; Simoncelli, Freeman et al. 1992).

The use of stereo as a quantitative shape cue was extended by a wide variety of shape-from-X techniques, including shape from shading (Figure 1.8b) (see Section 13.1.1 and Horn 1975; Pentland 1984; Blake, Zisserman, and Knowles 1985; Horn and Brooks 1986, 1989), photometric stereo (see Section 13.1.1 and Woodham 1981), shape from texture (see Section 13.1.2 and Witkin 1981; Pentland 1984; Malik and Rosenholtz 1997), and shape from focus (see Section 13.1.3 and Nayar, Watanabe, and Noguchi 1995). Horn (1986) has a nice discussion of most of these techniques.

Research into better edge and contour detection (Figure 1.8c) (see Section 7.2) was also active during this period (Canny 1986; Nalwa and Binford 1986), including the introduction of dynamically evolving contour trackers (Section 7.3.1) such as snakes (Kass, Witkin, and Terzopoulos 1988), as well as three-dimensional physically based models (Figure 1.8d) (Terzopoulos, Witkin, and Kass 1987; Kass, Witkin, and Terzopoulos 1988; Terzopoulos and Fleischer 1988).

Researchers noticed that a lot of the stereo, flow, shape-from-X, and edge detection algorithms could be unified, or at least described, using the same mathematical framework if they were posed as variational optimization problems and made more robust (well-posed) using regularization (Fig-

1.2 A brief history

ure 1.8e) (see Section 4.2 and Terzopoulos 1983; Poggio, Torre, and Koch 1985; Terzopoulos 1986b; Blake and Zisserman 1987; Bertero, Poggio, and Torre 1988; Terzopoulos 1988). Around the same time, Geman and Geman (1984) pointed out that such problems could equally well be formulated using discrete Markov random field (MRF) models (see Section 4.3), which enabled the use of better (global) search and optimization algorithms, such as simulated annealing.

Online variants of MRF algorithms that modeled and updated uncertainties using the Kalman filter were introduced a little later (Dickmanns and Graefe 1988; Matthies, Kanade, and Szeliski 1989; Szeliski 1989). Attempts were also made to map both regularized and MRF algorithms onto parallel hardware (Poggio and Koch 1985; Poggio, Little et al. 1988; Fischler, Firschein et al. 1989). The book by Fischler and Firschein (1987) contains a nice collection of articles focusing on all of these topics (stereo, flow, regularization, MRFs, and even higher-level vision).

Three-dimensional range data processing (acquisition, merging, modeling, and recognition; see Figure 1.8f) continued being actively explored during this decade (Agin and Binford 1976; Besl and Jain 1985; Faugeras and Hebert 1987; Curless and Levoy 1996). The compilation by Kanade (1987) contains a lot of the interesting papers in this area.

1990s. While a lot of the previously mentioned topics continued to be explored, a few of them became significantly more active.

A burst of activity in using projective invariants for recognition (Mundy and Zisserman 1992) evolved into a concerted effort to solve the structure from motion problem (see Chapter 11). A lot

of the initial activity was directed at projective reconstructions, which did not require knowledge of camera calibration (Faugeras 1992; Hartley, Gupta, and Chang 1992; Hartley 1994a; Faugeras and Luong 2001; Hartley and Zisserman 2004). Simultaneously, factorization techniques (Section 11.4.1) were developed to solve efficiently problems for which orthographic camera approximations were applicable (Figure 1.9a) (Tomasi and Kanade 1992; Poelman and Kanade 1997; Anandan and Irani 2002) and then later extended to the perspective case (Christy and Horaud 1996; Triggs 1996). Eventually, the field started using full global optimization (see Section 11.4.2 and Taylor, Kriegman, and Anandan 1991; Szeliski and Kang 1994; Azarbayejani and Pentland 1995), which was later recognized as being the same as the bundle adjustment techniques traditionally used in photogrammetry (Triggs, McLauchlan et al. 1999). Fully automated 3D modeling systems were built using such techniques (Beardsley, Torr, and Zisserman 1996; Schaffalitzky and Zisserman 2002; Snavely, Seitz, and Szeliski 2006; Agarwal, Furukawa et al. 2011; Frahm, Fite-Georgel et al. 2010).

Work begun in the 1980s on using detailed measurements of color and intensity combined with accurate physical models of radiance transport and color image formation created its own subfield known as physics-based vision. A good survey of the field can be found in the three-volume collection on this topic (Wolff, Shafer, and Healey 1992a; Healey and Shafer 1992; Shafer, Healey, and Wolff 1992).

Optical flow methods (see Chapter 9) continued to be improved (Nagel and Enkelmann 1986; Bolles, Baker, and Marimont 1987; Horn and Weldon Jr. 1988; Anandan 1989; Bergen, Anandan et al. 1992; Black and Anandan 1996; Bruhn, Weickert, and Schnörr 2005; Papenberg, Bruhn et al. 2006), with (Nagel 1986; Barron, Fleet, and Beauchemin 1994; Baker, Scharstein et al. 2011) being good surveys. Similarly, a lot of progress was made on dense stereo correspondence algorithms (see Chapter 12, Okutomi and Kanade (1993, 1994); Boykov, Veksler, and Zabih (1998); Birchfield and Tomasi (1999); Boykov, Veksler, and Zabih (2001), and the survey and comparison in Scharstein and Szeliski (2002)), with the biggest breakthrough being perhaps global optimization using graph cut techniques (Figure 1.9b) (Boykov, Veksler, and Zabih 2001).

Multi-view stereo algorithms (Figure 1.9c) that produce complete 3D surfaces (see Section 12.7) were also an active topic of research (Seitz and Dyer 1999; Kutulakos and Seitz 2000) that continues to be active today (Seitz, Curless et al. 2006; Schöps, Schönberger et al. 2017; Knapitsch, Park et al. 2017). Techniques for producing 3D volumetric descriptions from binary silhouettes (see Section 12.7.3) continued to be developed (Potmesil 1987; Srivasan, Liang, and Hackwood 1990; Szeliski 1993; Laurentini 1994), along with techniques based on tracking and reconstructing smooth occluding contours (see Section 12.2.1 and Cipolla and Blake 1992; Vaillant and Faugeras 1992; Zheng 1994; Boyer and Berger 1997; Szeliski and Weiss 1998; Cipolla and Giblin 2000).

Tracking algorithms also improved a lot, including contour tracking using active contours (see Section 7.3), such as snakes (Kass, Witkin, and Terzopoulos 1988), particle filters (Blake and Isard 1998), and level sets (Malladi, Sethian, and Vemuri 1995), as well as intensity-based (direct) techniques (Lucas and Kanade 1981; Shi and Tomasi 1994; Rehg and Kanade 1994), often applied to tracking faces (Figure 1.9d) (Lanitis, Taylor, and Cootes 1997; Matthews and Baker 2004; Matthews, Xiao, and Baker 2007) and whole bodies (Sidenbladh, Black, and Fleet 2000; Hilton, Fua, and Ronfard 2006; Moeslund, Hilton, and Krüger 2006).

Image segmentation (see Section 7.5) (Figure 1.9e), a topic which has been active since the earliest days of computer vision (Brice and Fennema 1970; Horowitz and Pavlidis 1976; Riseman and Arbib 1977; Rosenfeld and Davis 1979; Haralick and Shapiro 1985; Pavlidis and Liow 1990), was also an active topic of research, producing techniques based on minimum energy (Mumford and Shah 1989) and minimum description length (Leclerc 1989), normalized cuts (Shi and Malik 2000), and mean shift (Comaniciu and Meer 2002).

1.2 A brief history

Statistical learning techniques started appearing, first in the application of principal component eigenface analysis to face recognition (Figure 1.9f) (see Section 5.2.3 and Turk and Pentland 1991) and linear dynamical systems for curve tracking (see Section 7.3.1 and Blake and Isard 1998).

Perhaps the most notable development in computer vision during this decade was the increased interaction with computer graphics (Seitz and Szeliski 1999), especially in the cross-disciplinary area of image-based modeling and rendering (see Chapter 14). The idea of manipulating real-world imagery directly to create new animations first came to prominence with image morphing techniques (Figure 1.5c) (see Section 3.6.3 and Beier and Neely 1992) and was later applied to view interpolation (Chen and Williams 1993; Seitz and Dyer 1996), panoramic image stitching (Figure 1.5a) (see Section 8.2 and Mann and Picard 1994; Chen 1995; Szeliski 1996; Szeliski and Shum 1997; Szeliski 2006a), and full light-field rendering (Figure 1.10a) (see Section 14.3 and Gortler, Grzeszczuk et al. 1996; Levoy and Hanrahan 1996; Shade, Gortler et al. 1998). At the same time, image-based modeling techniques (Figure 1.10b) for automatically creating realistic 3D models from collections of images were also being introduced (Beardsley, Torr, and Zisserman 1996; Debevec, Taylor, and Malik 1996; Taylor, Debevec, and Malik 1996).

2000s. This decade continued to deepen the interplay between the vision and graphics fields, but more importantly embraced data-driven and learning approaches as core components of vision. Many of the topics introduced under the rubric of image-based rendering, such as image stitching (see Section 8.2), light-field capture and rendering (see Section 14.3), and high-dynamic range

(HDR) image capture through exposure bracketing (Figure 1.5b) (see Section 10.2 and Mann and Picard 1995; Debevec and Malik 1997), were re-christened as computational photography (see Chapter 10) to acknowledge the increased use of such techniques in everyday digital photography. For example, the rapid adoption of exposure bracketing to create high dynamic range images necessitated the development of tone mapping algorithms (Figure 1.10c) (see Section 10.2.1) to convert such images back to displayable results (Fattal, Lischinski, and Werman 2002; Durand and Dorsey 2002; Reinhard, Stark et al. 2002; Lischinski, Farbman et al. 2006). In addition to merging multiple exposures, techniques were developed to merge flash images with non-flash counterparts (Eisemann and Durand 2004; Petschnigg, Agrawala et al. 2004) and to interactively or automatically select different regions from overlapping images (Agarwala, Dontcheva et al. 2004).

Texture synthesis (Figure 1.10d) (see Section 10.5), quilting (Efros and Leung 1999; Efros and Freeman 2001; Kwatra, Schödl et al. 2003), and inpainting (Bertalmio, Sapiro et al. 2000; Bertalmio, Vese et al. 2003; Criminisi, Pérez, and Toyama 2004) are additional topics that can be classified as computational photography techniques, since they re-combine input image samples to produce new photographs.

A second notable trend during this decade was the emergence of feature-based techniques (combined with learning) for object recognition (see Section 6.1 and Ponce, Hebert et al. 2006). Some of the notable papers in this area include the constellation model of Fergus, Perona, and Zisserman (2007) (Figure 1.10e) and the pictorial structures of Felzenszwalb and Huttenlocher (2005). Feature-based techniques also dominate other recognition tasks, such as scene recognition (Zhang, Marszalek et al. 2007) and panorama and location recognition (Brown and Lowe 2007; Schindler, Brown, and Szeliski 2007). And while interest in point (patch-based) features tend to dominate current research, some groups are pursuing recognition based on contours (Belongie, Malik, and Puzicha 2002) and region segmentation (Figure 1.10f) (Mori, Ren et al. 2004).

Another significant trend from this decade was the development of more efficient algorithms for complex global optimization problems (see Chapter 4 and Appendix B.5 and Szeliski, Zabih et al. 2008; Blake, Kohli, and Rother 2011). While this trend began with work on graph cuts (Boykov, Veksler, and Zabih 2001; Kohli and Torr 2007), a lot of progress has also been made in message passing algorithms, such as loopy belief propagation (LBP) (Yedidia, Freeman, and Weiss 2001; Kumar and Torr 2006).

The most notable trend from this decade, which has by now completely taken over visual recognition and most other aspects of computer vision, was the application of sophisticated machine learning techniques to computer vision problems (see Chapters 5 and 6). This trend coincided with the increased availability of immense quantities of partially labeled data on the internet, as well as significant increases in computational power, which makes it more feasible to learn object categories without the use of careful human supervision.

2010s. The trend towards using large labeled (and also self-supervised) datasets to develop machine learning algorithms became a tidal wave that totally revolutionized the development of image recognition algorithms as well as other applications, such as denoising and optical flow, which previously used Bayesian and global optimization techniques.

This trend was enabled by the development of high-quality large-scale annotated datasets such as ImageNet (Deng, Dong et al. 2009; Russakovsky, Deng et al. 2015), Microsoft COCO (Common Objects in Context) (Lin, Maire et al. 2014), and LVIS (Gupta, Dollár, and Girshick 2019). These datasets provided not only reliable metrics for tracking the progress of recognition and semantic segmentation algorithms, but more importantly, sufficient labeled data to develop complete solutions based on machine learning.

1.2 A brief history

Another major trend was the dramatic increase in computational power available from the development of general purpose (data-parallel) algorithms on graphical processing units (GPGPU). The breakthrough SuperVision (“AlexNet”) deep neural network (Figure 1.11a; Krizhevsky, Sutskever, and Hinton 2012), which was the first neural network to win the yearly ImageNet large-scale visual recognition challenge, relied on GPU training, as well as a number of technical advances, for its dramatic performance. After the publication of this paper, progress in using deep convolutional architectures accelerated dramatically, to the point where they are now the only architecture considered for recognition and semantic segmentation tasks (Figure 1.11b), as well as the preferred architecture for many other vision tasks (Chapter 5; LeCun, Bengio, and Hinton 2015), including optical flow (Sun, Yang et al. 2018)), denoising, and monocular depth inference (Li, Dekel et al. 2019).

Large datasets and GPU architectures, coupled with the rapid dissemination of ideas through timely publications on arXiv as well as the development of languages for deep learning and the open sourcing of neural network models, all contributed to an explosive growth in this area, both in rapid advances and capabilities, and also in the sheer number of publications and researchers now working on these topics. They also enabled the extension of image recognition approaches to video understanding tasks such as action recognition (Feichtenhofer, Fan et al. 2019), as well as structured regression tasks such as real-time multi-person body pose estimation (Cao, Simon et al. 2017).

Specialized sensors and hardware for computer vision tasks also continued to advance. The

Microsoft Kinect depth camera, released in 2010, quickly became an essential component of many 3D modeling (Figure 1.11d) and person tracking (Shotton, Fitzgibbon et al. 2011) systems. Over the decade, 3D body shape modeling and tracking systems continued to evolve, to the point where it is now possible to infer a person's 3D model with gestures and expression from a single image (Figure 1.11c).

And while depth sensors have not yet become ubiquitous (except for security applications on high-end phones), computational photography algorithms run on all of today's smartphones. Innovations introduced in the computer vision community, such as panoramic image stitching and bracketed high dynamic range image merging, are now standard features, and multi-image low-light denoising algorithms are also becoming commonplace (Liba, Murthy et al. 2019). Lightfield imaging algorithms, which allow the creation of soft depth-of-field effects, are now also becoming more available (Garg, Wadhwa et al. 2019). Finally, mobile augmented reality applications that perform real-time pose estimation and environment augmentation using combinations of feature tracking and inertial measurements are commonplace, and are currently being extended to include pixel-accurate depth occlusion effects (Figure 1.11e).

On higher-end platforms such as autonomous vehicles and drones, powerful real-time SLAM (simultaneous localization and mapping) and VIO (visual inertial odometry) algorithms (Engel, Schöps, and Cremers 2014; Forster, Zhang et al. 2017; Engel, Koltun, and Cremers 2018) can build accurate 3D maps that enable, e.g., autonomous flight through challenging scenes such as forests (Figure 1.11f).

In summary, this past decade has seen incredible advances in the performance and reliability of computer vision algorithms, brought in part by the shift to machine learning and training on very large sets of real-world data. It has also seen the application of vision algorithms in myriad commercial and consumer scenarios as well as new challenges engendered by their widespread use (Su and Crandall 2021).

1.3 Book overview

In the final part of this introduction, I give a brief tour of the material in this book, as well as a few notes on notation and some additional general references. Since computer vision is such a broad field, it is possible to study certain aspects of it, e.g., geometric image formation and 3D structure recovery, without requiring other parts, e.g., the modeling of reflectance and shading. Some of the chapters in this book are only loosely coupled with others, and it is not strictly necessary to read all of the material in sequence.

Figure 1.12 shows a rough layout of the contents of this book. Since computer vision involves going from images to both a semantic understanding as well as a 3D structural description of the scene, I have positioned the chapters horizontally in terms of where in this spectrum they land, in addition to vertically according to their dependence. $ ^{9} $

Interspersed throughout the book are sample applications, which relate the algorithms and mathematical material being presented in various chapters to useful, real-world applications. Many of these applications are also presented in the exercises sections, so that students can write their own.

At the end of each section, I provide a set of exercises that the students can use to implement, test, and refine the algorithms and techniques presented in each section. Some of the exercises are suitable as written homework assignments, others as shorter one-week projects, and still others as

open-ended research problems that make for challenging final projects. Motivated students who implement a reasonable subset of these exercises will, by the end of the book, have a computer vision software library that can be used for a variety of interesting tasks and projects.

If the students or curriculum do not have a strong preference for programming languages, Python, with the NumPy scientific and array arithmetic library plus the OpenCV vision library, are a good environment to develop algorithms and learn about vision. Not only will the students learn how to program using array/tensor notation and linear/matrix algebra (which is a good foundation for later use of PyTorch for deep learning), you can also prepare classroom assignments using Jupyter notebooks, giving you the option to combine descriptive tutorials, sample code, and code to be extended/modified in one convenient location. $ ^{10} $

As this is a reference book, I try wherever possible to discuss which techniques and algorithms work well in practice, as well as provide up-to-date pointers to the latest research results in the areas that I cover. The exercises can be used to build up your own personal library of self-tested and validated vision algorithms, which is more worthwhile in the long term (assuming you have the time) than simply pulling algorithms out of a library whose performance you do not really understand.

The book begins in Chapter 2 with a review of the image formation processes that create the images that we see and capture. Understanding this process is fundamental if you want to take a scientific (model-based) approach to computer vision. Students who are eager to just start implementing algorithms (or courses that have limited time) can skip ahead to the next chapter and dip into this material later. In Chapter 2, we break down image formation into three major components. Geometric image formation (Section 2.1) deals with points, lines, and planes, and how these are mapped onto images using projective geometry and other models (including radial lens distortion). Photometric image formation (Section 2.2) covers radiometry, which describes how light interacts with surfaces in the world, and optics, which projects light onto the sensor plane. Finally, Section 2.3 covers how sensors work, including topics such as sampling and aliasing, color sensing, and in-camera compression.

Chapter 3 covers image processing, which is needed in almost all computer vision applications. This includes topics such as linear and non-linear filtering (Section 3.3), the Fourier transform (Section 3.4), image pyramids and wavelets (Section 3.5), and geometric transformations such as image warping (Section 3.6). Chapter 3 also presents applications such as seamless image blending and image morphing.

Chapter 4 begins with a new section on data fitting and interpolation, which provides a conceptual framework for global optimization techniques such as regularization and Markov random fields (MRFs), as well as machine learning, which we cover in the next chapter. Section 4.2 covers classic regularization techniques, i.e., piecewise-continuous smoothing splines (aka variational techniques) implemented using fast iterated linear system solvers, which are still often the method of choice in time-critical applications such as mobile augmented reality. The next section (4.3) presents the related topic of MRFs, which also serve as an introduction to Bayesian inference techniques, covered at a more abstract level in Appendix B. The chapter also discusses applications to interactive colorization and segmentation.

Chapter 5 is a completely new chapter covering machine learning, deep learning, and deep neural networks. It begins in Section 5.1 with a review of classic supervised machine learning approaches, which are designed to classify images (or regress values) based on intermediate-level features. Section 5.2 looks at unsupervised learning, which is useful for both understanding unlabeled training data and providing models of real-world distributions. Section 5.3 presents the basic elements of

feedforward neural networks, including weights, layers, and activation functions, as well as methods for network training. Section 5.4 goes into more detail on convolutional networks and their applications to both recognition and image processing. The last section in the chapter discusses more complex networks, including 3D, spatio-temporal, recurrent, and generative networks.

Chapter 6 covers the topic of recognition. In the first edition of this book this chapter came last, since it built upon earlier methods such as segmentation and feature matching. With the advent of deep networks, many of these intermediate representations are no longer necessary, since the network can learn them as part of the training process. As so much of computer vision research is now devoted to various recognition topics, I decided to move this chapter up so that students can learn about it earlier in the course.

The chapter begins with the classic problem of instance recognition, i.e., finding instances of known 3D objects in cluttered scenes. Section 6.2 covers both traditional and deep network approaches to whole image classification, i.e., what used to be called category recognition. It also discusses the special case of facial recognition. Section 6.3 presents algorithms for object detection (drawing bounding boxes around recognized objects), with a brief review of older approaches to face and pedestrian detection. Section 6.4 covers various flavors of semantic segmentation (generating per-pixel labels), including instance segmentation (delineating separate objects), pose estimation (labeling pixels with body parts), and panoptic segmentation (labeling both things and stuff). In Section 6.5, we briefly look at some recent papers in video understanding and action recognition, while in Section 6.6 we mention some recent work in image captioning and visual question answering.

In Chapter 7, we cover feature detection and matching. A lot of current 3D reconstruction and recognition techniques are built on extracting and matching feature points (Section 7.1), so this is a fundamental technique required by many subsequent chapters (Chapters 8 and 11) and even in instance recognition (Section 6.1). We also cover edge and straight line detection in Sections 7.2 and 7.4, contour tracking in Section 7.3, and low-level segmentation techniques in Section 7.5.

Feature detection and matching are used in Chapter 8 to perform image alignment (or registration) and image stitching. We introduce the basic techniques of feature-based alignment and show how this problem can be solved using either linear or non-linear least squares, depending on the motion involved. We also introduce additional concepts, such as uncertainty weighting and robust regression, which are essential to making real-world systems work. Feature-based alignment is then used as a building block for both 2D applications such as image stitching (Section 8.2) and computational photography (Chapter 10), as well as 3D geometric alignment tasks such as pose estimation and structure from motion (Chapter 11).

The second part of Chapter 8 is devoted to image stitching, i.e., the construction of large panoramas and composites. While stitching is just one example of computational photography (see Chapter 10), there is enough depth here to warrant a separate section. We start by discussing various possible motion models (Section 8.2.1), including planar motion and pure camera rotation. We then discuss global alignment (Section 8.3), which is a special (simplified) case of general bundle adjustment, and then present panorama recognition, i.e., techniques for automatically discovering which images actually form overlapping panoramas. Finally, we cover the topics of image compositing and blending (Section 8.4), which involve both selecting which pixels from which images to use and blending them together so as to disguise exposure differences.

Image stitching is a wonderful application that ties together most of the material covered in earlier parts of this book. It also makes for a good mid-term course project that can build on previously developed techniques such as image warping and feature detection and matching. Sections 8.2–8.4 also present more specialized variants of stitching such as whiteboard and document scanning, video summarization, panography, full $ 360^{\circ} $ spherical panoramas, and interactive photomontage for

blending repeated action shots together.

In Chapter 9, we generalize the concept of feature-based image alignment to cover dense intensity-based motion estimation, i.e., optical flow. We start with the simplest possible motion models, translational motion (Section 9.1), and cover topics such as hierarchical (coarse-to-fine) motion estimation, Fourier-based techniques, and iterative refinement. We then present parametric motion models, which can be used to compensate for camera rotation and zooming, as well as affine or planar perspective motion (Section 9.2). This is then generalized to spline-based motion models (Section 9.2.2) and finally to general per-pixel optical flow (Section 9.3). We close the chapter in Section 9.4 with a discussion of layered and learned motion models as well as video object segmentation and tracking. Applications of motion estimation techniques include automated morphing, video denoising, and frame interpolation (slow motion).

Chapter 10 presents additional examples of computational photography, which is the process of creating new images from one or more input photographs, often based on the careful modeling and calibration of the image formation process (Section 10.1). Computational photography techniques include merging multiple exposures to create high-dynamic-range images (Section 10.2), increasing image resolution through blur removal and super-resolution (Section 10.3), and image editing and compositing operations (Section 10.4). We also cover the topics of texture analysis, synthesis, and inpainting (hole filling) in Section 10.5, as well as non-photorealistic rendering and style transfer.

Starting in Chapter 11, we delve more deeply into techniques for reconstructing 3D models from images. We begin by introducing methods for intrinsic camera calibration in Section 11.1 and 3D pose estimation, i.e., extrinsic calibration, in Section 11.2. These sections also describe the applications of single-view reconstruction of building models and 3D location recognition. We then cover the topic of triangulation (Section 11.2.4), which is the 3D reconstruction of points from matched features when the camera positions are known.

Chapter 11 then moves on to the topic of structure from motion, which involves the simultaneous recovery of 3D camera motion and 3D scene structure from a collection of tracked 2D features. We begin with two-frame structure from motion (Section 11.3), for which algebraic techniques exist, as well as robust sampling techniques such as RANSAC that can discount erroneous feature matches. We then cover techniques for multi-frame structure from motion, including factorization (Section 11.4.1), bundle adjustment (Section 11.4.2), and constrained motion and structure models (Section 11.4.8). We present applications in visual effects (match move) and sparse 3D model construction for large (e.g., internet) photo collections. The final part of this chapter (Section 11.5) has a new section on simultaneous localization and mapping (SLAM) as well as its applications to autonomous navigation and mobile augmented reality (AR).

In Chapter 12, we turn to the topic of stereo correspondence, which can be thought of as a special case of motion estimation where the camera positions are already known (Section 12.1). This additional knowledge enables stereo algorithms to search over a much smaller space of correspondences to produce dense depth estimates using various combinations of matching criteria, optimization algorithm, and/or deep networks (Sections 12.3–12.6). We also cover multi-view stereo algorithms that build a true 3D surface representation instead of just a single depth map (Section 12.7), as well as monocular depth inference algorithms that hallucinate depth maps from just a single image (Section 12.8). Applications of stereo matching include head and gaze tracking, as well as depth-based background replacement (Z-keying).

Chapter 13 covers additional 3D shape and appearance modeling techniques. These include classic shape-from-X techniques such as shape from shading, shape from texture, and shape from focus (Section 13.1). An alternative to all of these passive computer vision techniques is to use active rangefinding (Section 13.2), i.e., to project patterned light onto scenes and recover the 3D geometry.

through triangulation. Processing all of these 3D representations often involves interpolating or simplifying the geometry (Section 13.3), or using alternative representations such as surface point sets (Section 13.4) or implicit functions (Section 13.5).

The collection of techniques for going from one or more images to partial or full 3D models is often called image-based modeling or 3D photography. Section 13.6 examines three more specialized application areas (architecture, faces, and human bodies), which can use model-based reconstruction to fit parameterized models to the sensed data. Section 13.7 examines the topic of appearance modeling, i.e., techniques for estimating the texture maps, albedos, or even sometimes complete bi-directional reflectance distribution functions (BRDFs) that describe the appearance of 3D surfaces.

In Chapter 14, we discuss the large number of image-based rendering techniques that have been developed in the last three decades, including simpler techniques such as view interpolation (Section 14.1), layered depth images (Section 14.2), and sprites and layers (Section 14.2.1), as well as the more general framework of light fields and Lumigraphs (Section 14.3) and higher-order fields such as environment mattes (Section 14.4). Applications of these techniques include navigating 3D collections of photographs using photo tourism.

Next, we discuss video-based rendering, which is the temporal extension of image-based rendering. The topics we cover include video-based animation (Section 14.5.1), periodic video turned into video textures (Section 14.5.2), and 3D video constructed from multiple video streams (Section 14.5.4). Applications of these techniques include animating still images and creating home tours based on $ 360^{\circ} $ video. We finish the chapter with an overview of the new emerging field of neural rendering.

To support the book's use as a textbook, the appendices and associated website contain more detailed mathematical topics and additional material. Appendix A covers linear algebra and numerical techniques, including matrix algebra, least squares, and iterative techniques. Appendix B covers Bayesian estimation theory, including maximum likelihood estimation, robust statistics, Markov random fields, and uncertainty modeling. Appendix C describes the supplementary material that can be used to complement this book, including images and datasets, pointers to software, and course slides.

1.4 Sample syllabus

Teaching all of the material covered in this book in a single quarter or semester course is a Herculean task and likely one not worth attempting. $ ^{11} $ It is better to simply pick and choose topics related to the lecturer's preferred emphasis and tailored to the set of mini-projects envisioned for the students.

Steve Seitz and I have successfully used a 10-week syllabus similar to the one shown in Table 1.1 as both an undergraduate and a graduate-level course in computer vision. The undergraduate course $ ^{12} $ tends to go lighter on the mathematics and takes more time reviewing basics, while the graduate-level course $ ^{13} $ dives more deeply into techniques and assumes the students already have a decent grounding in either vision or related mathematical techniques. Related courses have also been taught on the topics of 3D photography and computational photography. Appendix C.3 and the book's website list other courses that use this book to teach a similar curriculum.

When Steve and I teach the course, we prefer to give the students several small programming assignments early in the course rather than focusing on written homework or quizzes. With a suitable choice of topics, it is possible for these projects to build on each other. For example, introducing feature matching early on can be used in a second assignment to do image alignment and stitching. Alternatively, direct (optical flow) techniques can be used to do the alignment and more focus can be put on either graph cut seam selection or multi-resolution blending techniques.

In the past, we have also asked the students to propose a final project (we provide a set of suggested topics for those who need ideas) by the middle of the course and reserved the last week of the class for student presentations. Sometimes, a few of these projects have actually turned into conference submissions!

No matter how you decide to structure the course or how you choose to use this book, I encourage you to try at least a few small programming tasks to get a feel for how vision techniques work and how they fail. Better yet, pick topics that are fun and can be used on your own photographs, and try to push your creative boundaries to come up with surprising results.

1.5 A note on notation

For better or worse, the notation found in computer vision and multi-view geometry textbooks tends to vary all over the map (Faugeras 1993; Hartley and Zisserman 2004; Girod, Greiner, and Niemann 2000; Faugeras and Luong 2001; Forsyth and Ponce 2003). In this book, I use the convention I first learned in my high school physics class (and later multi-variate calculus and computer graphics courses), which is that vectors v are lower case bold, matrices M are upper case bold, and scalars $ (T, s) $ are mixed case italic. Unless otherwise noted, vectors operate as column vectors, i.e., they post-multiply matrices, $ \mathbf{M}\mathbf{v} $, although they are sometimes written as comma-separated parent-the-sized lists $ \mathbf{x} = (x, y) $ instead of bracketed column vectors $ \mathbf{x} = [x \ y]^T $. Some commonly used matrices are $ \mathbf{R} $ for rotations, $ \mathbf{K} $ for calibration matrices, and $ \mathbf{I} $ for the identity matrix. Homogeneous coordinates (Section 2.1) are denoted with a tilde over the vector, e.g., $ \tilde{\mathbf{x}} = (\tilde{x}, \tilde{y}, \tilde{w}) = \tilde{w}(x, y, 1) = \tilde{w}\mathbf{x} $ in $ \mathcal{P}^2 $. The cross product operator in matrix form is denoted by $ []\times $.

1.6 Additional reading

This book attempts to be self-contained, so that students can implement the basic assignments and algorithms described here without the need for outside references. However, it does presuppose a general familiarity with basic concepts in linear algebra and numerical techniques, which are reviewed in Appendix A, and image processing, which is reviewed in Chapter 3.

Students who want to delve more deeply into these topics can look in Golub and Van Loan (1996) for matrix algebra and Strang (1988) for linear algebra. In image processing, there are a number of popular textbooks, including Crane (1997), Gomes and Velho (1997), Jähne (1997), Pratt (2007), Russ (2007), Burger and Burge (2008), and Gonzalez and Woods (2017). For computer graphics, popular texts include Hughes, van Dam et al. (2013) and Marschner and Shirley (2015), with Glassner (1995) providing a more in-depth look at image formation and rendering. For statistics and machine learning, Chris Bishop's (2006) book is a wonderful and comprehensive introduction with a wealth of exercises, while Murphy (2012) provides a more recent take on the field and Hastie, Tibshirani, and Friedman (2009) a more classic treatment. A great introductory text to deep learning is Glassner (2018), while Goodfellow, Bengio, and Courville (2016) and Zhang, Lipton et al. (2021) provide more comprehensive treatments. Students may also want to look in other textbooks on computer vision for material that we do not cover here, as well as for additional project ideas (Nalwa 1993; Trucco and Verri 1998; Hartley and Zisserman 2004; Forsyth and Ponce 2011; Prince 2012; Davies 2017).

There is, however, no substitute for reading the latest research literature, both for the latest ideas and techniques and for the most up-to-date references to related literature. $ ^{14} $ In this book, I have attempted to cite the most recent work in each field so that students can read them directly and use them as inspiration for their own work. Browsing the last few years’ conference proceedings from the major vision, graphics, and machine learning conferences, such as CVPR, ECCV, ICCV, SIGGRAPH, and NeurIPS, as well as keeping an eye out for the latest publications on arXiv, will provide a wealth of new ideas. The tutorials offered at these conferences, for which slides or notes are often available online, are also an invaluable resource.

Chapter 2

Image formation

2.1 Geometric primitives and transformations ..... 29

2.1.1 2D transformations ..... 32

2.1.2 3D transformations ..... 35

2.1.3 3D rotations ..... 36

2.1.4 3D to 2D projections ..... 41

2.1.5 Lens distortions ..... 51

2.2 Photometric image formation ..... 53

2.2.1 Lighting ..... 53

2.2.2 Reflectance and shading ..... 54

2.2.3 Optics ..... 59

2.3 The digital camera ..... 63

2.3.1 Sampling and aliasing ..... 67

2.3.2 Color ..... 69

2.3.3 Compression ..... 78

2.4 Additional reading ..... 80

2.5 Exercises ..... 80

2.1 Geometric primitives and transformations

Before we can analyze and manipulate images, we need to establish a vocabulary for describing the geometry of a scene. We also need to understand the image formation process that produced a particular image given a set of lighting conditions, scene geometry, surface properties, and camera optics. In this chapter, we present a simplified model of this image formation process.

Section 2.1 introduces the basic geometric primitives used throughout the book (points, lines, and planes) and the geometric transformations that project these 3D quantities into 2D image features (Figure 2.1a). Section 2.2 describes how lighting, surface properties (Figure 2.1b), and camera optics (Figure 2.1c) interact to produce the color values that fall onto the image sensor. Section 2.3 describes how continuous color images are turned into discrete digital samples inside the image sensor (Figure 2.1d) and how to avoid (or at least characterize) sampling deficiencies, such as aliasing.

The material covered in this chapter is but a brief summary of a very rich and deep set of topics, traditionally covered in a number of separate fields. A more thorough introduction to the geometry of points, lines, planes, and projections can be found in textbooks on multi-view geometry (Hartley and Zisserman 2004; Faugeras and Luong 2001) and computer graphics (Hughes, van Dam et al. 2013). The image formation (synthesis) process is traditionally taught as part of a computer graphics curriculum (Glassner 1995; Watt 1995; Hughes, van Dam et al. 2013; Marschner and Shirley 2015) but it is also studied in physics-based computer vision (Wolff, Shafer, and Healey 1992a). The behavior of camera lens systems is studied in optics (Möller 1988; Ray 2002; Hecht 2015). Some good books on color theory are Healey and Shafer (1992), Wandell (1995), and Wyszecki and Stiles (2000), with Livingstone (2008) providing a more fun and informal introduction to the topic of color perception. Topics relating to sampling and aliasing are covered in textbooks on signal and image processing (Crane 1997; Jähne 1997; Oppenheim and Schafer 1996; Oppenheim, Schafer, and Buck 1999; Pratt 2007; Russ 2007; Burger and Burge 2008; Gonzalez and Woods 2017). The recent book by Ikeuchi, Matsushita et al. (2020) also covers 3D geometry, photometry, and sensor models, with an emphasis on active illumination systems.

A note to students: If you have already studied computer graphics, you may want to skim the material in Section 2.1, although the sections on projective depth and object-centered projection near the end of Section 2.1.4 may be new to you. Similarly, physics students (as well as computer graphics students) will mostly be familiar with Section 2.2. Finally, students with a good background in image processing will already be familiar with sampling issues (Section 2.3) as well as some of the material in Chapter 3.

2.1 Geometric primitives and transformations

In this section, we introduce the basic 2D and 3D primitives used in this textbook, namely points, lines, and planes. We also describe how 3D features are projected into 2D features. More detailed descriptions of these topics (along with a gentler and more intuitive introduction) can be found in textbooks on multiple-view geometry (Hartley and Zisserman 2004; Faugeras and Luong 2001).

Geometric primitives form the basic building blocks used to describe three-dimensional shapes. In this section, we introduce points, lines, and planes. Later sections of the book discuss curves (Sections 7.3 and 12.2), surfaces (Section 13.3), and volumes (Section 13.5).

2D points. 2D points (pixel coordinates in an image) can be denoted using a pair of values, $ \mathbf{x} = (x, y) \in \mathcal{R}^2 $, or alternatively,

$$ \mathbf{x}=\begin{bmatrix}x\\ y\end{bmatrix}. $$

(As stated in the introduction, we use the $ (x_{1}, x_{2}, \ldots) $ notation to denote column vectors.)

2D points can also be represented using homogeneous coordinates, $ \tilde{\mathbf{x}} = (\tilde{x}, \tilde{y}, \tilde{w}) \in \mathcal{P}^2 $, where vectors that differ only by scale are considered to be equivalent. $ \mathcal{P}^2 = \mathcal{R}^3 - (0, 0, 0) $ is called the 2D projective space.

A homogeneous vector $ \tilde{x} $ can be converted back into an inhomogeneous vector x by dividing through by the last element $ \tilde{w} $, i.e.,

$$ \tilde{\mathbf{x}}=(\tilde{x},\tilde{y},\tilde{w})=\tilde{w}(x,y,1)=\tilde{w}\bar{\mathbf{x}}, $$

where $ \bar{\mathbf{x}} = (x, y, 1) $ is the augmented vector. Homogeneous points whose last element is $ \tilde{w} = 0 $ are called ideal points or points at infinity and do not have an equivalent inhomogeneous representation.

2D lines. 2D lines can also be represented using homogeneous coordinates $ \mathbf{\tilde{I}} = (a, b, c) $. The corresponding line equation is

$$ \bar{\mathbf{x}}\cdot\mathbf{\tilde{l}}=a x+b y+c=0. $$

We can normalize the line equation vector so that $\mathbf{l} = (\hat{n}_x, \hat{n}_y, d) = (\hat{\mathbf{n}}, d)$ with $\|\hat{\mathbf{n}}\| = 1$. In this case, $\hat{\mathbf{n}}$ is the normal vector perpendicular to the line and $d$ is its distance to the origin (Figure 2.2). (The one exception to this normalization is the line $at\ infinity\ \mathbf{l} = (0, 0, 1)$, which includes all (ideal) points at infinity.)

We can also express $ \hat{n} $ as a function of rotation angle $ \theta $, $ \hat{n} = (\hat{n}_x, \hat{n}_y) = (\cos\theta, \sin\theta) $ (Figure 2.2a). This representation is commonly used in the Hough transform line-finding algorithm, which is discussed in Section 7.4.2. The combination $ (\theta, d) $ is also known as polar coordinates.

When using homogeneous coordinates, we can compute the intersection of two lines as

$$ \tilde{\mathbf{x}}=\tilde{\mathbf{l}}_{1}\times\tilde{\mathbf{l}}_{2}, $$

where $ \times $ is the cross product operator. Similarly, the line joining two points can be written as

$$ \tilde{1}=\tilde{\mathbf{x}}_{1}\times\tilde{\mathbf{x}}_{2}. $$

When trying to fit an intersection point to multiple lines or, conversely, a line to multiple points, least squares techniques (Section 8.1.1 and Appendix A.2) can be used, as discussed in Exercise 2.1.

2D conics. There are other algebraic curves that can be expressed with simple polynomial homogeneous equations. For example, the conic sections (so called because they arise as the intersection of a plane and a 3D cone) can be written using a quadric equation

$$ \tilde{\mathbf{x}}^{T}\mathbf{Q}\tilde{\mathbf{x}}=0. $$

Quadric equations play useful roles in the study of multi-view geometry and camera calibration (Hartley and Zisserman 2004; Faugeras and Luong 2001) but are not used extensively in this book.

3D points. Point coordinates in three dimensions can be written using inhomogeneous coordinates $ \mathbf{x} = (x, y, z) \in \mathcal{R}^3 $ or homogeneous coordinates $ \tilde{\mathbf{x}} = (\tilde{x}, \tilde{y}, \tilde{z}, \tilde{w}) \in \mathcal{P}^3 $. As before, it is sometimes useful to denote a 3D point using the augmented vector $ \tilde{\mathbf{x}} = (x, y, z, 1) $ with $ \tilde{\mathbf{x}} = \tilde{w}\tilde{\mathbf{x}} $.

3D planes. 3D planes can also be represented as homogeneous coordinates $ \tilde{\mathbf{m}} = (a, b, c, d) $ with a corresponding plane equation

$$ \mathbf{\bar{x}}\cdot\mathbf{\tilde{m}}=a x+b y+c z+d=0. $$

We can also normalize the plane equation as $ \mathbf{m} = (\hat{n}_x, \hat{n}_y, \hat{n}_z, d) = (\hat{\mathbf{n}}, d) $ with $ \|\hat{\mathbf{n}}\| = 1 $. In this case, $ \hat{\mathbf{n}} $ is the normal vector perpendicular to the plane and $ d $ is its distance to the origin (Figure 2.2b). As with the case of 2D lines, the plane at infinity $ \tilde{\mathbf{m}} = (0, 0, 0, 1) $, which contains all the points at infinity, cannot be normalized (i.e., it does not have a unique normal or a finite distance).

We can express $ \hat{n} $ as a function of two angles $ (\theta,\phi) $,

$$ \hat{\mathbf{n}}=(\cos\theta\cos\phi,\sin\theta\cos\phi,\sin\phi), $$

i.e., using spherical coordinates, but these are less commonly used than polar coordinates since they do not uniformly sample the space of possible normal vectors.

3D lines. Lines in 3D are less elegant than either lines in 2D or planes in 3D. One possible representation is to use two points on the line, $ (p, q) $. Any other point on the line can be expressed as a linear combination of these two points

$$ \mathbf{r}=(1-\lambda)\mathbf{p}+\lambda\mathbf{q}, $$

as shown in Figure 2.3. If we restrict $ 0 \leq \lambda \leq 1 $, we get the line segment joining p and q.

If we use homogeneous coordinates, we can write the line as

$$ \tilde{\mathbf{r}}=\mu\tilde{\mathbf{p}}+\lambda\tilde{\mathbf{q}}. $$

A special case of this is when the second point is at infinity, i.e., $ \tilde{\mathbf{q}} = (\hat{d}_x, \hat{d}_y, \hat{d}_z, 0) = (\hat{\mathbf{d}}, 0) $. Here, we see that $ \hat{\mathbf{d}} $ is the direction of the line. We can then re-write the inhomogeneous 3D line equation

$$ \mathbf{r}=\mathbf{p}+\lambda\hat{\mathbf{d}}. $$

A disadvantage of the endpoint representation for 3D lines is that it has too many degrees of freedom, i.e., six (three for each endpoint) instead of the four degrees that a 3D line truly has. However, if we fix the two points on the line to lie in specific planes, we obtain a representation with four degrees of freedom. For example, if we are representing nearly vertical lines, then $z=0$ and $z=1$ form two suitable planes, i.e., the $(x,y)$ coordinates in both planes provide the four coordinates describing the line. This kind of two-plane parameterization is used in the light field and Lumigraph image-based rendering systems described in Chapter 14 to represent the collection of rays seen by a camera as it moves in front of an object. The two-endpoint representation is also useful for representing line segments, even when their exact endpoints cannot be seen (only guessed at).

If we wish to represent all possible lines without bias towards any particular orientation, we can use Plücker coordinates (Hartley and Zisserman 2004, Section 3.2; Faugeras and Luong 2001, Chapter 3). These coordinates are the six independent non-zero entries in the $ 4 \times 4 $ skew symmetric matrix

$$ \mathbf{L}=\tilde{\mathbf{p}}\tilde{\mathbf{q}}^{T}-\tilde{\mathbf{q}}\tilde{\mathbf{p}}^{T}, $$

where $ \tilde{\mathbf{p}} $ and $ \tilde{\mathbf{q}} $ are any two (non-identical) points on the line. This representation has only four degrees of freedom, since $ \mathbf{L} $ is homogeneous and also satisfies $ |\mathbf{L}| = 0 $, which results in a quadratic constraint on the Plücker coordinates.

In practice, the minimal representation is not essential for most applications. An adequate model of 3D lines can be obtained by estimating their direction (which may be known ahead of time, e.g., for architecture) and some point within the visible portion of the line (see Section 11.4.8) or by using the two endpoints, since lines are most often visible as finite line segments. However, if you are interested in more details about the topic of minimal line parameterizations, Förstner (2005) discusses various ways to infer and model 3D lines in projective geometry, as well as how to estimate the uncertainty in such fitted models.

3D quadrics. The 3D analog of a conic section is a quadric surface

$$ \bar{\mathbf{x}}^{T}\mathbf{Q}\bar{\mathbf{x}}=0 $$

(Hartley and Zisserman 2004, Chapter 3). Again, while quadric surfaces are useful in the study of multi-view geometry and can also serve as useful modeling primitives (spheres, ellipsoids, cylinders), we do not study them in great detail in this book.

2.1.1 2D transformations

Having defined our basic primitives, we can now turn our attention to how they can be transformed. The simplest transformations occur in the 2D plane are illustrated in Figure 2.4.

Translation. 2D translations can be written as $ x' = x + t $ or

$$ \mathbf{x}^{\prime}=\begin{bmatrix}\mathbf{I}&\mathbf{t}\end{bmatrix}\bar{\mathbf{x}}, $$

where I is the $ (2 \times 2) $ identity matrix or

$$ \mathbf{\bar{x}}^{\prime}=\left[\begin{matrix}{\mathbf{I}}&{\mathbf{t}}\\ {\mathbf{0}^{T}}&{1}\\ \end{matrix}\right]\mathbf{\bar{x}}, $$

where 0 is the zero vector. Using a $ 2 \times 3 $ matrix results in a more compact notation, whereas using a full-rank $ 3 \times 3 $ matrix (which can be obtained from the $ 2 \times 3 $ matrix by appending a $ [0^T 1] $ row)

makes it possible to chain transformations using matrix multiplication as well as to compute inverse transforms. Note that in any equation where an augmented vector such as $ \bar{x} $ appears on both sides, it can always be replaced with a full homogeneous vector $ \tilde{x} $.

Rotation + translation. This transformation is also known as 2D rigid body motion or the 2D Euclidean transformation (since Euclidean distances are preserved). It can be written as $ \mathbf{x}' = \mathbf{R}\mathbf{x} + \mathbf{t} $ or

$$ \mathbf{x}^{\prime}=\begin{bmatrix}\mathbf{R}&\mathbf{t}\end{bmatrix}\mathbf{\bar{x}}. $$

where

$$ \mathbf{R}=\begin{bmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{bmatrix} $$

is an orthonormal rotation matrix with $ \mathbf{R}\mathbf{R}^T = \mathbf{I} $ and $ |\mathbf{R}| = 1 $.

Scaled rotation. Also known as the similarity transform, this transformation can be expressed as $ \mathbf{x}' = s\mathbf{R}\mathbf{x} + \mathbf{t} $, where $ s $ is an arbitrary scale factor. It can also be written as

$$ \mathbf{x}^{\prime}=\begin{bmatrix}s\mathbf{R}&\mathbf{t}\end{bmatrix}\mathbf{\bar{x}}=\begin{bmatrix}a&-b&t_{x}\\ b&a&t_{y}\end{bmatrix}\mathbf{\bar{x}}, $$

where we no longer require that $ a^{2} + b^{2} = 1 $. The similarity transform preserves angles between lines.

Affine. The affine transformation is written as $ \mathbf{x}' = \mathbf{A}\bar{\mathbf{x}} $, where $ \mathbf{A} $ is an arbitrary $ 2 \times 3 $ matrix, i.e.,

$$ \mathbf{x}^{\prime}=\begin{bmatrix}a_{00}&a_{01}&a_{02}\\ a_{10}&a_{11}&a_{12}\end{bmatrix}\overline{\mathbf{x}}. $$

Parallel lines remain parallel under affine transformations.

Projective. This transformation, also known as a perspective transform or homography, operates on homogeneous coordinates.

$$ \tilde{\mathbf{x}}^{\prime}=\tilde{\mathbf{H}}\tilde{\mathbf{x}}, $$

where $ \tilde{\mathbf{H}} $ is an arbitrary $ 3 \times 3 $ matrix. Note that $ \tilde{\mathbf{H}} $ is homogeneous, i.e., it is only defined up to a scale, and that two $ \tilde{\mathbf{H}} $ matrices that differ only by scale are equivalent. The resulting homogeneous coordinate $ \tilde{\mathbf{x}}' $ must be normalized in order to obtain an inhomogeneous result $ \mathbf{x} $, i.e.,

$$ x^{\prime}=\frac{h_{00}x+h_{01}y+h_{02}}{h_{20}x+h_{21}y+h_{22}}\quad\mathrm{a n d}\quad y^{\prime}=\frac{h_{10}x+h_{11}y+h_{12}}{h_{20}x+h_{21}y+h_{22}}. $$

| Transformation | Matrix | # DoF | Preserves | Icon |

| --- | --- | --- | --- | --- |

| translation | $ [\mathbf{I} \mathbf{t}]_{2\times3} $ | 2 | orientation | ☐ |

| rigid (Euclidean) | $ [\mathbf{R} \mathbf{t}]_{2\times3} $ | 3 | lengths | ◇ |

| similarity | $ [s\mathbf{R} \mathbf{t}]_{2\times3} $ | 4 | angles | ◇ |

| affine | $ [\mathbf{A}]_{2\times3} $ | 6 | parallels | ◇ |

| projective | $ [\tilde{\mathbf{H}}]_{3\times3} $ | 8 | straight lines | ◇ |

Perspective transformations preserve straight lines (i.e., they remain straight after the transformation).

Hierarchy of 2D transformations. The preceding set of transformations are illustrated in Figure 2.4 and summarized in Table 2.1. The easiest way to think of them is as a set of (potentially restricted) $ 3 \times 3 $ matrices operating on 2D homogeneous coordinate vectors. Hartley and Zisserman (2004) contains a more detailed description of the hierarchy of 2D planar transformations.

The above transformations form a nested set of groups, i.e., they are closed under composition and have an inverse that is a member of the same group. (This will be important later when applying these transformations to images in Section 3.6.) Each (simpler) group is a subgroup of the more complex group below it. The mathematics of such Lie groups and their related algebras (tangent spaces at the origin) are discussed in a number of recent robotics tutorials (Dellaert and Kaess 2017; Blanco 2019; Solà, Deray, and Atchuthan 2019), where the 2D rotation and rigid transforms are called SO(2) and SE(2), which stand for the special orthogonal and special Euclidean groups. $ ^{1} $

Co-vectors. While the above transformations can be used to transform points in a 2D plane, can they also be used directly to transform a line equation? Consider the homogeneous equation $ \tilde{1} \cdot \tilde{x} = 0 $. If we transform $ \tilde{x}' = \tilde{H} \tilde{x} $, we obtain

$$ \tilde{\mathbf{l}}^{\prime}\cdot\tilde{\mathbf{x}}^{\prime}=\tilde{\mathbf{l}}^{\prime T}\tilde{\mathbf{H}}\tilde{\mathbf{x}}=(\tilde{\mathbf{H}}^{T}\tilde{\mathbf{l}}^{\prime})^{T}\tilde{\mathbf{x}}=\tilde{\mathbf{l}}\cdot\tilde{\mathbf{x}}=0, $$

i.e., $ \tilde{\mathbf{I}}^{\prime} = \tilde{\mathbf{H}}^{-T}\tilde{\mathbf{I}} $. Thus, the action of a projective transformation on a co-vector such as a 2D line or 3D normal can be represented by the transposed inverse of the matrix, which is equivalent to the adjoint of $ \tilde{\mathbf{H}} $, since projective transformation matrices are homogeneous. Jim Blinn (1998) describes (in Chapters 9 and 10) the ins and outs of notating and manipulating co-vectors.

While the above transformations are the ones we use most extensively, a number of additional transformations are sometimes used.

2.1 Geometric primitives and transformations

Stretch/squash. This transformation changes the aspect ratio of an image.

$$ x^{\prime}=s_{x}x+t_{x} $$

$$ y^{\prime}=s_{y}y+t_{y}, $$

and is a restricted form of an affine transformation. Unfortunately, it does not nest cleanly with the groups listed in Table 2.1.

Planar surface flow. This eight-parameter transformation (Horn 1986; Bergen, Anandan et al. 1992; Girod, Greiner, and Niemann 2000).

$$ \begin{aligned}&x^{\prime}=a_{0}+a_{1}x+a_{2}y+a_{6}x^{2}+a_{7}xy\\&y^{\prime}=a_{3}+a_{4}x+a_{5}y+a_{6}xy+a_{7}y^{2},\\ \end{aligned} $$

arises when a planar surface undergoes a small 3D motion. It can thus be thought of as a small motion approximation to a full homography. Its main attraction is that it is linear in the motion parameters, $ a_{k} $, which are often the quantities being estimated.

Bilinear interpolant. This eight-parameter transform (Wolberg 1990).

$$ x^{\prime}=a_{0}+a_{1}x+a_{2}y+a_{6}xy $$

$$ y^{\prime}=a_{3}+a_{4}x+a_{5}y+a_{7}xy, $$

can be used to interpolate the deformation due to the motion of the four corner points of a square. (In fact, it can interpolate the motion of any four non-collinear points.) While the deformation is linear in the motion parameters, it does not generally preserve straight lines (only lines parallel to the square axes). However, it is often quite useful, e.g., in the interpolation of sparse grids using splines (Section 9.2.2).

2.1.2 3D transformations

The set of three-dimensional coordinate transformations is very similar to that available for 2D transformations and is summarized in Table 2.2. As in 2D, these transformations form a nested set of groups. Hartley and Zisserman (2004, Section 2.4) give a more detailed description of this hierarchy.

Translation. 3D translations can be written as $ x' = x + t $ or

$$ \mathbf{x}^{\prime}=\begin{bmatrix}\mathbf{I}&\mathbf{t}\end{bmatrix}\bar{\mathbf{x}}, $$

where I is the $ (3 \times 3) $ identity matrix.

Rotation + translation. Also known as 3D rigid body motion or the 3D Euclidean transformation or SE(3), it can be written as $ \mathbf{x}' = \mathbf{R}\mathbf{x} + \mathbf{t} $ or

$$ \mathbf{x}^{\prime}=\begin{bmatrix}\mathbf{R}&\mathbf{t}\end{bmatrix}\bar{\mathbf{x}}, $$

where $ \mathbf{R} $ is a $ 3 \times 3 $ orthonormal rotation matrix with $ \mathbf{RR}^T = \mathbf{I} $ and $ |\mathbf{R}| = 1 $. Note that sometimes it is more convenient to describe a rigid motion using

$$ \mathbf{x}^{\prime}=\mathbf{R}(\mathbf{x}-\mathbf{c})=\mathbf{R}\mathbf{x}-\mathbf{R}\mathbf{c}, $$

| Transformation | Matrix | # DoF | Preserves | Icon |

| --- | --- | --- | --- | --- |

| translation | $ [\mathbf{I} \mathbf{t}]_{3\times4} $ | 3 | orientation | ☐ |

| rigid (Euclidean) | $ [\mathbf{R} \mathbf{t}]_{3\times4} $ | 6 | lengths | ◇ |

| similarity | $ [s\mathbf{R} \mathbf{t}]_{3\times4} $ | 7 | angles | ◇ |

| affine | $ [\mathbf{A}]_{3\times4} $ | 12 | parallels | ◇ |

| projective | $ [\tilde{\mathbf{H}}]_{4\times4} $ | 15 | straight lines | ◇ |

where c is the center of rotation (often the camera center).

Compactly parameterizing a 3D rotation is a non-trivial task, which we describe in more detail below.

Scaled rotation. The 3D similarity transform can be expressed as $ \mathbf{x}' = s\mathbf{R}\mathbf{x} + \mathbf{t} $ where $ s $ is an arbitrary scale factor. It can also be written as

$$ \mathbf{x}^{\prime}=\begin{bmatrix}s\mathbf{R}&\mathbf{t}\end{bmatrix}\mathbf{\bar{x}}. $$

This transformation preserves angles between lines and planes.

Affine. The affine transform is written as $ \mathbf{x}' = \mathbf{A}\bar{\mathbf{x}} $, where $ \mathbf{A} $ is an arbitrary $ 3 \times 4 $ matrix, i.e.,

$$ \mathbf{x}^{\prime}=\begin{bmatrix}a_{00}&a_{01}&a_{02}&a_{03}\\ a_{10}&a_{11}&a_{12}&a_{13}\\ a_{20}&a_{21}&a_{22}&a_{23}\end{bmatrix}\mathbf{\bar{x}}. $$

Parallel lines and planes remain parallel under affine transformations.

Projective. This transformation, variously known as a 3D perspective transform, homography, or collineation, operates on homogeneous coordinates.

$$ \tilde{\mathbf{x}}^{\prime}=\tilde{\mathbf{H}}\tilde{\mathbf{x}}, $$

where $ \mathbf{H} $ is an arbitrary $ 4 \times 4 $ homogeneous matrix. As in 2D, the resulting homogeneous coordinate $ \tilde{\mathbf{x}}' $ must be normalized in order to obtain an inhomogeneous result $ \mathbf{x} $. Perspective transformations preserve straight lines (i.e., they remain straight after the transformation).

2.1.3 3D rotations

The biggest difference between 2D and 3D coordinate transformations is that the parameterization of the 3D rotation matrix R is not as straightforward, as several different possibilities exist.

Euler angles

A rotation matrix can be formed as the product of three rotations around three cardinal axes, e.g., $x$, $y$, and $z$, or $x$, $y$, and $x$. This is generally a bad idea, as the result depends on the order in which the transforms are applied.$^{2}$ What is worse, it is not always possible to move smoothly in the parameter space, i.e., sometimes one or more of the Euler angles change dramatically in response to a small change in rotation.$^{3}$ For these reasons, we do not even give the formula for Euler angles in this book—interested readers can look in other textbooks or technical reports (Faugeras 1993; Diebel 2006). Note that, in some applications, if the rotations are known to be a set of uni-axial transforms, they can always be represented using an explicit set of rigid transformations.

Axis/angle (exponential twist)

A rotation can be represented by a rotation axis $ \hat{n} $ and an angle $ \theta $, or equivalently by a 3D vector $ \boldsymbol{\omega} = \theta \hat{\mathbf{n}} $. Figure 2.5 shows how we can compute the equivalent rotation. First, we project the vector $ \mathbf{v} $ onto the axis $ \hat{n} $ to obtain

$$ \mathbf{v}_{\parallel}=\hat{\mathbf{n}}(\hat{\mathbf{n}}\cdot\mathbf{v})=(\hat{\mathbf{n}}\hat{\mathbf{n}}^{T})\mathbf{v}, $$

which is the component of v that is not affected by the rotation. Next, we compute the perpendicular residual of v from $ \hat{n} $,

$$ \mathbf{v}_{\perp}=\mathbf{v}-\mathbf{v}_{\parallel}=(\mathbf{I}-\hat{\mathbf{n}}\hat{\mathbf{n}}^{T})\mathbf{v}. $$

We can rotate this vector by $ 90^{\circ} $ using the cross product.

$$ \mathbf{v}_{\times}=\hat{\mathbf{n}}\times\mathbf{v}_{\perp}=\hat{\mathbf{n}}\times\mathbf{v}=[\hat{\mathbf{n}}]_{\times}\mathbf{v}, $$

where $ [\hat{n}]_\times $ is the matrix form of the cross product operator with the vector $ \hat{n} = (\hat{n}_x, \hat{n}_y, \hat{n}_z) $,

$$ [\mathbf{\hat{n}}]_{\times}=\left[\begin{matrix}{0}&{-\hat{n}_{z}}&{\hat{n}_{y}}\\ {\hat{n}_{z}}&{0}&{-\hat{n}_{x}}\\ {-\hat{n}_{y}}&{\hat{n}_{x}}&{0}\\ \end{matrix}\right]. $$

Note that rotating this vector by another $ 90^{\circ} $ is equivalent to taking the cross product again.

$$ \mathbf{v}_{\times\times}=\hat{\mathbf{n}}\times\mathbf{v}_{\times}=[\hat{\mathbf{n}}]_{\times}^{2}\mathbf{v}=-\mathbf{v}_{\perp}, $$

and hence

$$ \mathbf{v}_{\parallel}=\mathbf{v}-\mathbf{v}_{\perp}=\mathbf{v}+\mathbf{v}_{\times\times}=(\mathbf{I}+[\hat{\mathbf{n}}]_{\times}^{2})\mathbf{v}. $$

We can now compute the in-plane component of the rotated vector u as

$$ \mathbf{u}_{\perp}=\cos\theta\mathbf{v}_{\perp}+\sin\theta\mathbf{v}_{\times}=(\sin\theta[\hat{\mathbf{n}}]_{\times}-\cos\theta[\hat{\mathbf{n}}]_{\times}^{2})\mathbf{v}. $$

Putting all these terms together, we obtain the final rotated vector as

$$ \mathbf{u}=\mathbf{u}_{\perp}+\mathbf{v}_{\parallel}=(\mathbf{I}+\sin\theta[\hat{\mathbf{n}}]_{\times}+(1-\cos\theta)[\hat{\mathbf{n}}]_{\times}^{2})\mathbf{v}. $$

We can therefore write the rotation matrix corresponding to a rotation by $ \theta $ around an axis $ \hat{n} $ as

$$ \mathbf{R}(\hat{\mathbf{n}},\theta)=\mathbf{I}+\sin\theta[\hat{\mathbf{n}}]_{\times}+(1-\cos\theta)[\hat{\mathbf{n}}]_{\times}^{2}, $$

which is known as Rodrigues' formula (Ayache 1989).

The product of the axis $ \hat{\mathbf{n}} $ and angle $ \theta $, $ \boldsymbol{\omega} = \theta \hat{\mathbf{n}} = (\omega_x, \omega_y, \omega_z) $, is a minimal representation for a 3D rotation. Rotations through common angles such as multiples of $ 90^\circ $ can be represented exactly (and converted to exact matrices) if $ \theta $ is stored in degrees. Unfortunately, this representation is not unique, since we can always add a multiple of $ 360^\circ $ ( $ 2\pi $ radians) to $ \theta $ and get the same rotation matrix. As well, $ (\hat{\mathbf{n}}, \theta) $ and $ (-\hat{\mathbf{n}}, -\theta) $ represent the same rotation.

However, for small rotations (e.g., corrections to rotations), this is an excellent choice. In particular, for small (infinitesimal or instantaneous) rotations and $\theta$ expressed in radians, Rodrigues' formula simplifies to

$$ \mathbf{R}(\boldsymbol{\omega})\approx\mathbf{I}+\sin\theta[\hat{\mathbf{n}}]_{\times}\approx\mathbf{I}+[\theta\hat{\mathbf{n}}]_{\times}=\begin{bmatrix}1&-\omega_{z}&\omega_{y}\\ \omega_{z}&1&-\omega_{x}\\ -\omega_{y}&\omega_{x}&1\end{bmatrix}, $$

which gives a nice linearized relationship between the rotation parameters $ \omega $ and $ \mathbf{R} $. We can also write $ \mathbf{R}(\omega)\mathbf{v} \approx \mathbf{v} + \omega \times \mathbf{v} $, which is handy when we want to compute the derivative of $ \mathbf{R}\mathbf{v} $ with respect to $ \omega $,

$$ \frac{\partial\mathbf{R}\mathbf{v}}{\partial\mathbf{\omega}^{T}}=-\left[\mathbf{v}\right]_{\times}=\left[\begin{matrix}{0}&{z}&{-y}\\ {-z}&{0}&{x}\\ {y}&{-x}&{0}\\ \end{matrix}\right]. $$

Another way to derive a rotation through a finite angle is called the exponential twist (Murray, Li, and Sastry 1994). A rotation by an angle $ \theta $ is equivalent to $ k $ rotations through $ \theta/k $. In the limit as $ k \to \infty $, we obtain

$$ \mathbf{R}(\hat{\mathbf{n}},\theta)=\lim_{k\to\infty}(\mathbf{I}+\frac{1}{k}[\theta\hat{\mathbf{n}}]_{\times})^{k}=\exp\left[\boldsymbol{\omega}\right]_{\times}. $$

If we expand the matrix exponential as a Taylor series (using the identity $ [\hat{\mathbf{n}}]_{\times}^{k+2} = -[\hat{\mathbf{n}}]_{\times}^{k}, k > 0 $, and again assuming $ \theta $ is in radians),

$$ \begin{align*}\exp\left[\boldsymbol{\omega}\right]_{\times}&=\mathbf{I}+\theta[\hat{\mathbf{n}}]_{\times}+\frac{\theta^{2}}{2}[\hat{\mathbf{n}}]_{\times}^{2}+\frac{\theta^{3}}{3!}[\hat{\mathbf{n}}]_{\times}^{3}+\cdots\\&=\mathbf{I}+(\theta-\frac{\theta^{3}}{3!}+\cdots)[\hat{\mathbf{n}}]_{\times}+(\frac{\theta^{2}}{2}-\frac{\theta^{4}}{4!}+\cdots)[\hat{\mathbf{n}}]_{\times}^{2}\\&=\mathbf{I}+\sin\theta[\hat{\mathbf{n}}]_{\times}+(1-\cos\theta)[\hat{\mathbf{n}}]_{\times}^{2},\end{align*} $$

which yields the familiar Rodrigues' formula.

In robotics (and group theory), rotations are called SO(3), i.e., the special orthogonal group in 3D. The incremental rotations $\omega$ are associated with a Lie algebra $\mathfrak{se}(3)$ and are the preferred way to formulate rotation derivatives and to model uncertainties in rotation estimates (Blanco 2019; Solà, Deray, and Atchuthan 2019).

Unit quaternions

The unit quaternion representation is closely related to the angle/axis representation. A unit quaternion is a unit length 4-vector whose components can be written as $ \mathbf{q} = (q_x, q_y, q_z, q_w) $ or $ \mathbf{q} = (x, y, z, w) $ for short. Unit quaternions live on the unit sphere $ \|\mathbf{q}\| = 1 $ and antipodal (opposite sign) quaternions, $ \mathbf{q} $ and $ -\mathbf{q} $, represent the same rotation (Figure 2.6). Other than this ambiguity (dual covering), the unit quaternion representation of a rotation is unique. Furthermore, the representation is continuous, i.e., as rotation matrices vary continuously, you can find a continuous quaternion representation, although the path on the quaternion sphere may wrap all the way around before returning to the “origin” $ \mathbf{q}_o = (0, 0, 0, 1) $. For these and other reasons given below, quaternions are a very popular representation for pose and for pose interpolation in computer graphics (Shoemake 1985).

Quaternions can be derived from the axis/angle representation through the formula

$$ \mathbf{q}=(\mathbf{v},w)=(\sin\frac{\theta}{2}\mathbf{\hat{n}},\cos\frac{\theta}{2}), $$

where $ \hat{n} $ and $ \theta $ are the rotation axis and angle. Using the trigonometric identities $ \sin\theta = 2\sin\frac{\theta}{2}\cos\frac{\theta}{2} $ and $ (1 - \cos\theta) = 2\sin^{2}\frac{\theta}{2} $, Rodrigues' formula can be converted to

$$ \begin{aligned}\mathbf{R}(\hat{\mathbf{n}},\theta)&=\mathbf{I}+\sin\theta[\hat{\mathbf{n}}]_{\times}+(1-\cos\theta)[\hat{\mathbf{n}}]_{\times}^{2}\\&=\mathbf{I}+2w[\mathbf{v}]_{\times}+2[\mathbf{v}]_{\times}^{2}.\end{aligned} $$

This suggests a quick way to rotate a vector v by a quaternion using a series of cross products, scalings, and additions. To obtain a formula for $ \mathbf{R}(\mathbf{q}) $ as a function of $ (x, y, z, w) $, recall that

$$ [\mathbf{v}]_{\times}=\left[\begin{matrix}{0}&{-z}&{y}\\ {z}&{0}&{-x}\\ {-y}&{x}&{0}\\ \end{matrix}\right]\quad\mathrm{a n d}\quad[\mathbf{v}]_{\times}^{2}=\left[\begin{matrix}{-y^{2}-z^{2}}&{x y}&{x z}\\ {x y}&{-x^{2}-z^{2}}&{y z}\\ {x z}&{y z}&{-x^{2}-y^{2}}\\ \end{matrix}\right]. $$

We thus obtain

$$ \mathbf{R}(\mathbf{q})=\left[\begin{matrix}{1-2(y^{2}+z^{2})}&{2(x y-z w)}&{2(x z+y w)}\\ {2(x y+z w)}&{1-2(x^{2}+z^{2})}&{2(y z-x w)}\\ {2(x z-y w)}&{2(y z+x w)}&{1-2(x^{2}+y^{2})}\\ \end{matrix}\right]. $$

The diagonal terms can be made more symmetrical by replacing $ 1 - 2(y^2 + z^2) $ with $ (x^2 + w^2 - y^2 - z^2) $, etc.

The nicest aspect of unit quaternions is that there is a simple algebra for composing rotations expressed as unit quaternions. Given two quaternions $ \mathbf{q}_0 = (\mathbf{v}_0, w_0) $ and $ \mathbf{q}_1 = (\mathbf{v}_1, w_1) $, the quaternion multiply operator is defined as

$$ \mathbf{q}_{2}=\mathbf{q}_{0}\mathbf{q}_{1}=\big(\mathbf{v}_{0}\times\mathbf{v}_{1}+w_{0}\mathbf{v}_{1}+w_{1}\mathbf{v}_{0},~w_{0}w_{1}-\mathbf{v}_{0}\cdot\mathbf{v}_{1}\big), $$

with the property that $ \mathbf{R}(\mathbf{q}_2) = \mathbf{R}(\mathbf{q}_0)\mathbf{R}(\mathbf{q}_1) $. Note that quaternion multiplication is not commutative, just as 3D rotations and matrix multiplications are not.

Taking the inverse of a quaternion is easy: Just flip the sign of v or w (but not both!). (You can verify this has the desired effect of transposing the $ \mathbf{R} $ matrix in (2.41).) Thus, we can also define quaternion division as

$$ \mathbf{q}_{2}=\mathbf{q}_{0}/\mathbf{q}_{1}=\mathbf{q}_{0}\mathbf{q}_{1}^{-1}=\big(\mathbf{v}_{0}\times\mathbf{v}_{1}+w_{0}\mathbf{v}_{1}-w_{1}\mathbf{v}_{0},-w_{0}w_{1}-\mathbf{v}_{0}\cdot\mathbf{v}_{1}\big). $$

This is useful when the incremental rotation between two rotations is desired.

In particular, if we want to determine a rotation that is partway between two given rotations, we can compute the incremental rotation, take a fraction of the angle, and compute the new rotation. This procedure is called spherical linear interpolation or slerp for short (Shoemaker 1985) and is given in Algorithm 2.1. Note that Shoemaker presents two formulas other than the one given here. The first exponentiates $ q_r $ by alpha before multiplying the original quaternion,

$$ \mathbf{q}_{2}=\mathbf{q}_{r}^{\alpha}\mathbf{q}_{0}, $$

while the second treats the quaternions as 4-vectors on a sphere and uses

$$ \mathbf{q}_{2}=\frac{\sin(1-\alpha)\theta}{\sin\theta}\mathbf{q}_{0}+\frac{\sin\alpha\theta}{\sin\theta}\mathbf{q}_{1}, $$

where $\theta = \cos^{-1}(\mathbf{q}_0 \cdot \mathbf{q}_1)$ and the dot product is directly between the quaternion 4-vectors. All of these formulas give comparable results, although care should be taken when $\mathbf{q}_0$ and $\mathbf{q}_1$ are close together, which is why I prefer to use an arctangent to establish the rotation angle.

Which rotation representation is better?

The choice of representation for 3D rotations depends partly on the application.

The axis/angle representation is minimal, and hence does not require any additional constraints on the parameters (no need to re-normalize after each update). If the angle is expressed in degrees, it is easier to understand the pose (say, $ 90^\circ $ twist around $ x $-axis), and also easier to express exact rotations. When the angle is in radians, the derivatives of $ \mathbf{R} $ with respect to $ \omega $ can easily be computed (2.36).

Quaternions, on the other hand, are better if you want to keep track of a smoothly moving camera, since there are no discontinuities in the representation. It is also easier to interpolate between rotations and to chain rigid transformations (Murray, Li, and Sastry 1994; Bregler and Malik 1998).

My usual preference is to use quaternions, but to update their estimates using an incremental rotation, as described in Section 11.2.2.

procedure slerp(q₀, q₁, α):

1. $ \mathbf{q}_r = \mathbf{q}_1 / \mathbf{q}_0 = (\mathbf{v}_r, w_r) $

2. if $ w_r < 0 $ then $ \mathbf{q}_r \leftarrow -\mathbf{q}_r $

3. $ \theta_r = 2 \tan^{-1}(\|\mathbf{v}_r\|/w_r) $

4. $ \hat{\mathbf{n}}_r = \mathcal{N}(\mathbf{v}_r) = \mathbf{v}_r / \|\mathbf{v}_r\| $

5. $ \theta_\alpha = \alpha \theta_r $

6. $ \mathbf{q}_\alpha = (\sin \frac{\theta_\alpha}{2} \hat{\mathbf{n}}_r, \cos \frac{\theta_\alpha}{2}) $

7. $ \mathbf{return}\mathbf{q}_2 = \mathbf{q}_\alpha \mathbf{q}_0 $

2.1.4 3D to 2D projections

Now that we know how to represent 2D and 3D geometric primitives and how to transform them spatially, we need to specify how 3D primitives are projected onto the image plane. We can do this using a linear 3D to 2D projection matrix. The simplest model is orthography, which requires no division to get the final (inhomogeneous) result. The more commonly used model is perspective, since this more accurately models the behavior of real cameras.

Orthography and para-perspective

An orthographic projection simply drops the z component of the three-dimensional coordinate p to obtain the 2D point x. (In this section, we use p to denote 3D points and x to denote 2D points.) This can be written as

$$ \mathbf{x}=\left[\mathbf{I}_{2\times2}|\mathbf{0}\right]\mathbf{p}. $$

If we are using homogeneous (projective) coordinates, we can write

$$ \tilde{\mathbf{x}}=\begin{bmatrix}1&0&0&0\\ 0&1&0&0\\ 0&0&0&1\end{bmatrix}\tilde{\mathbf{p}}, $$

i.e., we drop the z component but keep the w component. Orthography is an approximate model for long focal length (telephoto) lenses and objects whose depth is shallow relative to their distance to the camera (Sawhney and Hanson 1991). It is exact only for telecentric lenses (Baker and Nayar 1999, 2001).

In practice, world coordinates (which may measure dimensions in meters) need to be scaled to fit onto an image sensor (physically measured in millimeters, but ultimately measured in pixels). For this reason, scaled orthography is actually more commonly used,

$$ \mathbf{x}=\left[s\mathbf{I}_{2\times2}|\mathbf{0}\right]\mathbf{p}. $$

2.1 Geometric primitives and transformations

This model is equivalent to first projecting the world points onto a local fronto-parallel image plane and then scaling this image using regular perspective projection. The scaling can be the same for all parts of the scene (Figure 2.7b) or it can be different for objects that are being modeled independently (Figure 2.7c). More importantly, the scaling can vary from frame to frame when estimating structure from motion, which can better model the scale change that occurs as an object approaches the camera.

Scaled orthography is a popular model for reconstructing the 3D shape of objects far away from the camera, since it greatly simplifies certain computations. For example, pose (camera orientation) can be estimated using simple least squares (Section 11.2.1). Under orthography, structure and motion can simultaneously be estimated using factorization (singular value decomposition), as discussed in Section 11.4.1 (Tomasi and Kanade 1992).

A closely related projection model is para-perspective (Aloimonos 1990; Poelman and Kanade 1997). In this model, object points are again first projected onto a local reference parallel to the image plane. However, rather than being projected orthogonally to this plane, they are projected parallel to the line of sight to the object center (Figure 2.7d). This is followed by the usual projection onto the final image plane, which again amounts to a scaling. The combination of these two projections is therefore affine and can be written as

$$ \tilde{\mathbf{x}}=\begin{bmatrix}a_{00}&a_{01}&a_{02}&a_{03}\\ a_{10}&a_{11}&a_{12}&a_{13}\\ 0&0&0&1\end{bmatrix}\tilde{\mathbf{p}}. $$

Note how parallel lines in 3D remain parallel after projection in Figure 2.7b–d. Para-perspective provides a more accurate projection model than scaled orthography, without incurring the added complexity of per-pixel perspective division, which invalidates traditional factorization methods (Poelman and Kanade 1997).

Perspective

The most commonly used projection in computer graphics and computer vision is true 3D perspective (Figure 2.7e). Here, points are projected onto the image plane by dividing them by their z component. Using inhomogeneous coordinates, this can be written as

$$ \mathbf{\bar{x}}=\mathcal{P}_{z}(\mathbf{p})=\left[\begin{matrix}{x/z}\\ {y/z}\\ {1}\\ \end{matrix}\right]. $$

In homogeneous coordinates, the projection has a simple linear form.

$$ \mathbf{\widetilde{x}}=\left[\begin{matrix}{1}&{0}&{0}&{0}\\ {0}&{1}&{0}&{0}\\ {0}&{0}&{1}&{0}\\ \end{matrix}\right]\mathbf{\widetilde{p}}, $$

i.e., we drop the w component of p. Thus, after projection, it is not possible to recover the distance of the 3D point from the image, which makes sense for a 2D imaging sensor.

A form often seen in computer graphics systems is a two-step projection that first projects 3D coordinates into normalized device coordinates $ (x, y, z) \in [-1, 1] \times [-1, 1] \times [0, 1] $, and then rescales these coordinates to integer pixel coordinates using a viewport transformation (Watt 1995; OpenGL

ARB 1997). The (initial) perspective projection is then represented using a $ 4 \times 4 $ matrix

$$ \tilde{\mathbf{x}}=\left[\begin{matrix}{1}&{0}&{0}&{0}\\ {0}&{1}&{0}&{0}\\ {0}&{0}&{-z_{\mathrm{f a r}}/z_{\mathrm{r a n g e}}}&{z_{\mathrm{n e a r}}z_{\mathrm{f a r}}/z_{\mathrm{r a n g e}}}\\ {0}&{0}&{1}&{0}\\ \end{matrix}\right]\tilde{\mathbf{p}}, $$

where $ z_{near} $ and $ z_{far} $ are the near and far z clipping planes and $ z_{range} = z_{far} - z_{near} $. Note that the first two rows are actually scaled by the focal length and the aspect ratio so that visible rays are mapped to $ (x, y, z) \in [-1, 1]^2 $. The reason for keeping the third row, rather than dropping it, is that visibility operations, such as z-buffering, require a depth for every graphical element that is being rendered.

If we set $ z_{\text{near}} = 1 $, $ z_{\text{far}} \to \infty $, and switch the sign of the third row, the third element of the normalized screen vector becomes the inverse depth, i.e., the disparity (Okutomi and Kanade 1993). This can be quite convenient in many cases since, for cameras moving around outdoors, the inverse depth to the camera is often a more well-conditioned parameterization than direct 3D distance.

While a regular 2D image sensor has no way of measuring distance to a surface point, range sensors (Section 13.2) and stereo matching algorithms (Chapter 12) can compute such values. It is then convenient to be able to map from a sensor-based depth or disparity value d directly back to a 3D location using the inverse of a $ 4 \times 4 $ matrix (Section 2.1.4). We can do this if we represent perspective projection using a full-rank $ 4 \times 4 $ matrix, as in (2.64).

Camera intrinsics

Once we have projected a 3D point through an ideal pinhole using a projection matrix, we must still transform the resulting coordinates according to the pixel sensor spacing and the relative position of the sensor plane to the origin. Figure 2.8 shows an illustration of the geometry involved. In this section, we first present a mapping from 2D pixel coordinates to 3D rays using a sensor homography $ M_s $, since this is easier to explain in terms of physically measurable quantities. We then relate these quantities to the more commonly used camera intrinsic matrix $ \mathbf{K} $, which is used to map 3D camera-centered points $ \mathbf{p}_c $ to 2D pixel coordinates $ \tilde{\mathbf{x}}_s $.

Image sensors return pixel values indexed by integer pixel coordinates $ (x_s, y_s) $, often with the coordinates starting at the upper-left corner of the image and moving down and to the right. (This convention is not obeyed by all imaging libraries, but the adjustment for other coordinate systems is straightforward.) To map pixel centers to 3D coordinates, we first scale the $ (x_s, y_s) $ values by the pixel spacings $ (s_x, s_y) $ (sometimes expressed in microns for solid-state sensors) and then describe

2.1 Geometric primitives and transformations

the orientation of the sensor array relative to the camera projection center $ \mathbf{O}_c $ with an origin $ \mathbf{c}_s $ and a 3D rotation $ \mathbf{R}_s $ (Figure 2.8).

The combined 2D to 3D projection can then be written as

$$ \mathbf{p}=\begin{bmatrix}\mathbf{R}_{s}&\mathbf{c}_{s}\end{bmatrix}\begin{bmatrix}s_{x}&0&0\\ 0&s_{y}&0\\ 0&0&0\\ 0&0&1\end{bmatrix}\begin{bmatrix}x_{s}\\ y_{s}\\ 1\end{bmatrix}=\mathbf{M}_{s}\mathbf{\bar{x}}_{s}. $$

The first two columns of the $ 3 \times 3 $ matrix $ \mathbf{M}_s $ are the 3D vectors corresponding to unit steps in the image pixel array along the $ x_s $ and $ y_s $ directions, while the third column is the 3D image array origin $ \mathbf{c}_s $.

The matrix $ M_s $ is parameterized by eight unknowns: the three parameters describing the rotation $ R_s $, the three parameters describing the translation $ c_s $, and the two scale factors ( $ s_x, s_y $). Note that we ignore here the possibility of skew between the two axes on the image plane, since solid-state manufacturing techniques render this negligible. In practice, unless we have accurate external knowledge of the sensor spacing or sensor orientation, there are only seven degrees of freedom, since the distance of the sensor from the origin cannot be teased apart from the sensor spacing, based on external image measurement alone.

However, estimating a camera model $ M_{s} $ with the required seven degrees of freedom (i.e., where the first two columns are orthogonal after an appropriate re-scaling) is impractical, so most practitioners assume a general $ 3 \times 3 $ homogeneous matrix form.

The relationship between the 3D pixel center p and the 3D camera-centered point $p_{c}$ is given by an unknown scaling $s$, $\mathbf{p} = s\mathbf{p}_{c}$. We can therefore write the complete projection between $\mathbf{p}_{c}$ and a homogeneous version of the pixel address $\tilde{\mathbf{x}}_{s}$ as

$$ \tilde{\mathbf{x}}_{s}=\alpha\mathbf{M}_{s}^{-1}\mathbf{p}_{c}=\mathbf{K}\mathbf{p}_{c}. $$

The $ 3 \times 3 $ matrix K is called the calibration matrix and describes the camera intrinsics (as opposed to the camera's orientation in space, which are called the extrinsics).

From the above discussion, we see that $ \mathbf{K} $ has seven degrees of freedom in theory and eight degrees of freedom (the full dimensionality of a $ 3 \times 3 $ homogeneous matrix) in practice. Why, then, do most textbooks on 3D computer vision and multi-view geometry (Faugeras 1993; Hartley and Zisserman 2004; Faugeras and Luong 2001) treat $ \mathbf{K} $ as an upper-triangular matrix with five degrees of freedom?

While this is usually not made explicit in these books, it is because we cannot recover the full K matrix based on external measurement alone. When calibrating a camera (Section 11.1) based on external 3D points or other measurements (Tsai 1987), we end up estimating the intrinsic (K) and extrinsic (R, t) camera parameters simultaneously using a series of measurements,

$$ \mathbf{\tilde{x}}_{s}=\mathbf{K}\left[\mathbf{R}\quad\mathbf{t}\right]\mathbf{p}_{w}=\mathbf{P}\mathbf{p}_{w}, $$

where $ p_{w} $ are known 3D world coordinates and

$$ \mathbf{P}=\mathbf{K}[\mathbf{R}|\mathbf{t}] $$

is known as the camera matrix. Inspecting this equation, we see that we can post-multiply K by $ \mathbf{R}_1 $ and pre-multiply $ [\mathbf{R}|t] $ by $ \mathbf{R}_1^T $, and still end up with a valid calibration. Thus, it is impossible based on image measurements alone to know the true orientation of the sensor and the true camera intrinsics.

The choice of an upper-triangular form for $ \mathbf{K} $ seems to be conventional. Given a full $ 3 \times 4 $ camera matrix $ \mathbf{P} = \mathbf{K}[\mathbf{R}|t] $, we can compute an upper-triangular $ \mathbf{K} $ matrix using QR factorization (Golub and Van Loan 1996). (Note the unfortunate clash of terminologies: In matrix algebra textbooks, $ \mathbf{R} $ represents an upper-triangular (right of the diagonal) matrix; in computer vision, $ \mathbf{R} $ is an orthogonal rotation.)

There are several ways to write the upper-triangular form of K. One possibility is

$$ \mathbf{K}=\begin{bmatrix}f_{x}&s&c_{x}\\ 0&f_{y}&c_{y}\\ 0&0&1\end{bmatrix}, $$

which uses independent focal lengths $ f_x $ and $ f_y $ for the sensor x and y dimensions. The entry s encodes any possible skew between the sensor axes due to the sensor not being mounted perpendicular to the optical axis and $ (c_x, c_y) $ denotes the image center expressed in pixel coordinates. The image center is also often called the principal point in the computer vision literature (Hartley and Zisserman 2004), although in optics, the principal points are 3D points usually inside the lens where the principal planes intersect the principal (optical) axis (Hecht 2015). Another possibility is

$$ \mathbf{K}=\begin{bmatrix}f&s&c_{x}\\ 0&a f&c_{y}\\ 0&0&1\end{bmatrix}, $$

where the aspect ratio a has been made explicit and a common focal length f is used.

In practice, for many applications an even simpler form can be obtained by setting a = 1 and s = 0,

$$ \mathbf{K}=\begin{bmatrix}f&0&c_{x}\\ 0&f&c_{y}\\ 0&0&1\end{bmatrix}. $$

Often, setting the origin at roughly the center of the image, e.g., $ (c_x, c_y) = (W/2, H/2) $, where $ W $ and $ H $ are the image width and height, respectively, can result in a perfectly usable camera model with a single unknown, i.e., the focal length $ f $.

Figure 2.9 shows how these quantities can be visualized as part of a simplified imaging model. Note that now we have placed the image plane in front of the nodal point (projection center of the lens). The sense of the y-axis has also been flipped to get a coordinate system compatible with the way that most imaging libraries treat the vertical (row) coordinate.

A note on focal lengths

The issue of how to express focal lengths is one that often causes confusion in implementing computer vision algorithms and discussing their results. This is because the focal length depends on the units used to measure pixels.

If we number pixel coordinates using integer values, say $[0, W) \times [0, H)$, the focal length $f$ and camera center $(c_x, c_y)$ in (2.59) can be expressed as pixel values. How do these quantities relate to the more familiar focal lengths used by photographers?

Figure 2.10 illustrates the relationship between the focal length $f$, the sensor width $W$, and the horizontal field of view $\theta_{H}$, which obey the formula

$$ \tan\frac{\theta_{\mathrm{H}}}{2}=\frac{W}{2f}\qquad\mathrm{o r}\qquad f=\frac{W}{2}\left[\tan\frac{\theta_{\mathrm{H}}}{2}\right]^{-1}. $$

For a traditional 35mm film camera, whose active exposure area is 24mm × 36mm, we have $ W = 36\text{mm} $, and hence $ f $ is also expressed in millimeters. $ ^{4} $ For example, the “stock” lens that often comes with SLR (single lens reflex) cameras is 50mm, which is a good length, whereas 85mm is the standard for portrait photography. Since we work with digital images, however, it is more convenient to express $ W $ in pixels so that the focal length $ f $ can be used directly in the calibration matrix $ \mathbf{K} $ as in (2.59).

Another possibility is to scale the pixel coordinates so that they go from $ [-1, 1) $ along the longer image dimension and $ [-a^{-1}, a^{-1}) $ along the shorter axis, where $ a \geq 1 $ is the image aspect ratio (as opposed to the sensor cell aspect ratio introduced earlier). This can be accomplished using modified normalized device coordinates,

$$ x_{s}^{\prime}=(2x_{s}-W)/S\quad\mathrm{a n d}\quad y_{s}^{\prime}=(2y_{s}-H)/S,\qquad\mathrm{w h e r e}\qquad S=\operatorname*{m a x}(W,H). $$

This has the advantage that the focal length $f$ and image center $(c_x, c_y)$ become independent of the image resolution, which can be useful when using multi-resolution, image-processing algorithms, such as image pyramids (Section 3.5).$^{5}$ The use of $S$ instead of $W$ also makes the focal length the same for landscape (horizontal) and portrait (vertical) pictures, as is the case in 35mm photography. (In some computer graphics textbooks and systems, normalized device coordinates go from $[-1, 1] \times [-1, 1]$, which requires the use of two different focal lengths to describe the camera intrinsics (Watt

1995.) Setting $S = W = 2$ in (2.60), we obtain the simpler (unitless) relationship

$$ f^{-1}=\tan\frac{\theta_H}{2}. $$

The conversion between the various focal length representations is straightforward, e.g., to go from a unitless $f$ to one expressed in pixels, multiply by $W/2$, while to convert from an $f$ expressed in pixels to the equivalent 35mm focal length, multiply by 18mm.

Camera matrix

Now that we have shown how to parameterize the calibration matrix K, we can put the camera intrinsics and extrinsics together to obtain a single $ 3 \times 4 $ camera matrix

$$ \mathbf{P}=\mathbf{K}\left[\mathbf{R}\quad\mathbf{t}\right]. $$

It is sometimes preferable to use an invertible $ 4 \times 4 $ matrix, which can be obtained by not dropping the last row in the P matrix,

$$ \tilde{\mathbf{P}}=\begin{bmatrix}\mathbf{K}&\mathbf{0}\\ \mathbf{0}^{T}&1\end{bmatrix}\begin{bmatrix}\mathbf{R}&\mathbf{t}\\ \mathbf{0}^{T}&1\end{bmatrix}=\tilde{\mathbf{K}}\mathbf{E}, $$

where $\mathbf{E}$ is a 3D rigid-body (Euclidean) transformation and $\mathbf{\tilde{K}}$ is the full-rank calibration matrix. The $4 \times 4$ camera matrix $\tilde{\mathbf{P}}$ can be used to map directly from 3D world coordinates $\tilde{\mathbf{p}}_{w} = (x_{w}, y_{w}, z_{w}, 1)$ to screen coordinates (plus disparity), $\mathbf{x}_{s} = (x_{s}, y_{s}, 1, d)$,

$$ \mathbf{x}_{s}\sim\mathbf{\tilde{P}}\mathbf{\bar{p}}_{w}, $$

where $\sim$ indicates equality up to scale. Note that after multiplication by $\mathbf{P}$, the vector is divided by the third element of the vector to obtain the normalized form $\mathbf{x}_s = (x_s, y_s, 1, d)$.

Plane plus parallax (projective depth)

In general, when using the $ 4 \times 4 $ matrix $ \mathbf{P} $, we have the freedom to remap the last row to whatever suits our purpose (rather than just being the “standard” interpretation of disparity as inverse depth). Let us re-write the last row of $ \tilde{\mathbf{P}} $ as $ \mathbf{p}_3 = s_3[\hat{\mathbf{n}}_0 | c_0] $, where $ \|\hat{\mathbf{n}}_0\| = 1 $. We then have the equation

$$ d=\frac{s_{3}}{z}(\hat{\mathbf{n}}_{0}\cdot\mathbf{p}_{w}+c_{0}), $$

2.1 Geometric primitives and transformations

where $ z = \mathbf{p}_2 \cdot \mathbf{\bar{p}}_w = \mathbf{r}_z \cdot (\mathbf{p}_w - \mathbf{c}) $ is the distance of $ \mathbf{p}_w $ from the camera center $ C $ (2.25) along the optical axis $ Z $ (Figure 2.11). Thus, we can interpret $ d $ as the projective disparity or projective depth of a 3D scene point $ \mathbf{p}_w $ from the reference plane $ \mathbf{\hat{n}}_0 \cdot \mathbf{p}_w + c_0 = 0 $ (Szeliski and Coughlan 1997; Szeliski and Golland 1999; Shade, Gortler et al. 1998; Baker, Szeliski, and Anandan 1998). (The projective depth is also sometimes called $ parallax $ in reconstruction algorithms that use the term plane plus $ parallax $ (Kumar, Anandan, and Hanna 1994; Sawhney 1994).) Setting $ \mathbf{\hat{n}}_0 = \mathbf{0} $ and $ c_0 = 1 $, i.e., putting the reference plane at infinity, results in the more standard $ d = 1/z $ version of disparity (Okutomi and Kanade 1993).

Another way to see this is to invert the P matrix so that we can map pixels plus disparity directly back to 3D points,

$$ \tilde{\mathbf{p}}_{w}=\tilde{\mathbf{P}}^{-1}\mathbf{x}_{s}. $$

In general, we can choose P to have whatever form is convenient, i.e., to sample space using an arbitrary projection. This can come in particularly handy when setting up multi-view stereo reconstruction algorithms, since it allows us to sweep a series of planes (Section 12.1.2) through space with a variable (projective) sampling that best matches the sensed image motions (Collins 1996; Szeliski and Golland 1999; Saito and Kanade 1999).

Mapping from one camera to another

What happens when we take two images of a 3D scene from different camera positions or orientations (Figure 2.12a)? Using the full rank $ 4 \times 4 $ camera matrix $ \tilde{\mathbf{P}} = \tilde{\mathbf{K}}\mathbf{E} $ from (2.64), we can write the projection from world to screen coordinates as

$$ \tilde{\mathbf{x}}_{0}\sim\tilde{\mathbf{K}}_{0}\mathbf{E}_{0}\mathbf{p}=\tilde{\mathbf{P}}_{0}\mathbf{p}. $$

Assuming that we know the z-buffer or disparity value $ d_{0} $ for a pixel in one image, we can compute the 3D point location p using

$$ \mathbf{p}\sim\mathbf{E}_{0}^{-1}\tilde{\mathbf{K}}_{0}^{-1}\tilde{\mathbf{x}}_{0} $$

and then project it into another image yielding

$$ \tilde{\mathbf{x}}_{1}\sim\tilde{\mathbf{K}}_{1}\mathbf{E}_{1}\mathbf{p}=\tilde{\mathbf{K}}_{1}\mathbf{E}_{1}\mathbf{E}_{0}^{-1}\tilde{\mathbf{K}}_{0}^{-1}\tilde{\mathbf{x}}_{0}=\tilde{\mathbf{P}}_{1}\tilde{\mathbf{P}}_{0}^{-1}\tilde{\mathbf{x}}_{0}=\mathbf{M}_{10}\tilde{\mathbf{x}}_{0}. $$

Unfortunately, we do not usually have access to the depth coordinates of pixels in a regular photographic image. However, for a planar scene, as discussed above in (2.66), we can replace the last row of $ \mathbf{P}_0 $ in (2.64) with a general plane equation, $ \hat{\mathbf{n}}_0 \cdot \mathbf{p} + c_0 $, that maps points on the plane to $ d_0 = 0 $ values (Figure 2.12b). Thus, if we set $ d_0 = 0 $, we can ignore the last column of $ \mathbf{M}_{10} $ in (2.70) and also its last row, since we do not care about the final z-buffer depth. The mapping Equation (2.70) thus reduces to

$$ \tilde{\mathbf{x}}_{1}\sim\tilde{\mathbf{H}}_{10}\tilde{\mathbf{x}}_{0}, $$

where $\tilde{\mathbf{H}}_{10}$ is a general $3\times3$ homography matrix and $\tilde{\mathbf{x}}_{1}$ and $\tilde{\mathbf{x}}_{0}$ are now 2D homogeneous coordinates (i.e., 3-vectors) (Szeliski 1996). This justifies the use of the 8-parameter homography as a general alignment model for mosaics of planar scenes (Mann and Picard 1994; Szeliski 1996).

The other special case where we do not need to know depth to perform inter-camera mapping is when the camera is undergoing pure rotation (Section 8.2.3), i.e., when $ t_0 = t_1 $. In this case, we can write

$$ \tilde{\mathbf{x}}_{1}\sim\mathbf{K}_{1}\mathbf{R}_{1}\mathbf{R}_{0}^{-1}\mathbf{K}_{0}^{-1}\tilde{\mathbf{x}}_{0}=\mathbf{K}_{1}\mathbf{R}_{10}\mathbf{K}_{0}^{-1}\tilde{\mathbf{x}}_{0}, $$

which again can be represented with a $ 3 \times 3 $ homography. If we assume that the calibration matrices have known aspect ratios and centers of projection (2.59), this homography can be parameterized by the rotation amount and the two unknown focal lengths. This particular formulation is commonly used in image-stitching applications (Section 8.2.3).

Object-centered projection

When working with long focal length lenses, it often becomes difficult to reliably estimate the focal length from image measurements alone. This is because the focal length and the distance to the object are highly correlated and it becomes difficult to tease these two effects apart. For example, the change in scale of an object viewed through a zoom telephoto lens can either be due to a zoom change or to a motion towards the user. (This effect was put to dramatic use in some scenes of Alfred Hitchcock's film Vertigo, where the simultaneous change of zoom and camera motion produces a disquieting effect.)

This ambiguity becomes clearer if we write out the projection equation corresponding to the simple calibration matrix $ \mathbf{K} $ (2.59),

$$ x_{s}=f\frac{\mathbf{r}_{x}\cdot\mathbf{p}+t_{x}}{\mathbf{r}_{z}\cdot\mathbf{p}+t_{z}}+c_{x} $$

$$ \boldsymbol{y}_{s}=f\frac{\mathbf{r}_{y}\cdot\mathbf{p}+t_{y}}{\mathbf{r}_{z}\cdot\mathbf{p}+t_{z}}+c_{y}, $$

where $ \mathbf{r}_x $, $ \mathbf{r}_y $, and $ \mathbf{r}_z $ are the three rows of $ \mathbf{R} $. If the distance to the object center $ t_z \gg \|\mathbf{p}\| $ (the size of the object), the denominator is approximately $ t_z $ and the overall scale of the projected object depends on the ratio of $ f $ to $ t_z $. It therefore becomes difficult to disentangle these two quantities.

To see this more clearly, let $ \eta_{z}=t_{z}^{-1} $ and $ s=\eta_{z}f $. We can then re-write the above equations as

$$ x_{s}=s\frac{\mathbf{r}_{x}\cdot\mathbf{p}+t_{x}}{1+\eta_{z}\mathbf{r}_{z}\cdot\mathbf{p}}+c_{x} $$

$$ y_{s}=s\frac{\mathbf{r}_{y}\cdot\mathbf{p}+t_{y}}{1+\eta_{z}\mathbf{r}_{z}\cdot\mathbf{p}}+c_{y} $$

(Szeliski and Kang 1994; Pighin, Hecker et al. 1998). The scale of the projection s can be reliably estimated if we are looking at a known object (i.e., the 3D coordinates p are known). The inverse

2.1 Geometric primitives and transformations

distance $ \eta_{z} $ is now mostly decoupled from the estimates of s and can be estimated from the amount of foreshortening as the object rotates. Furthermore, as the lens becomes longer, i.e., the projection model becomes orthographic, there is no need to replace a perspective imaging model with an orthographic one, since the same equation can be used, with $ \eta_{z} \to 0 $ (as opposed to f and $ t_{z} $ both going to infinity). This allows us to form a natural link between orthographic reconstruction techniques such as factorization and their projective/perspective counterparts (Section 11.4.1).

2.1.5 Lens distortions

The above imaging models all assume that cameras obey a linear projection model where straight lines in the world result in straight lines in the image. (This follows as a natural consequence of linear matrix operations being applied to homogeneous coordinates.) Unfortunately, many wide-angle lenses have noticeable radial distortion, which manifests itself as a visible curvature in the projection of straight lines. (See Section 2.2.3 for a more detailed discussion of lens optics, including chromatic aberration.) Unless this distortion is taken into account, it becomes impossible to create highly accurate photorealistic reconstructions. For example, image mosaics constructed without taking radial distortion into account will often exhibit blurring due to the misregistration of corresponding features before pixel blending (Section 8.2).

Fortunately, compensating for radial distortion is not that difficult in practice. For most lenses, a simple quartic model of distortion can produce good results. Let $ (x_c, y_c) $ be the pixel coordinates obtained after perspective division but before scaling by focal length f and shifting by the image center $ (c_x, c_y) $, i.e.,

$$ \begin{aligned}\boldsymbol{x}_{c}&=\frac{\mathbf{r}_{x}\cdot\mathbf{p}+t_{x}}{\mathbf{r}_{z}\cdot\mathbf{p}+t_{z}}\\y_{c}&=\frac{\mathbf{r}_{y}\cdot\mathbf{p}+t_{y}}{\mathbf{r}_{z}\cdot\mathbf{p}+t_{z}}.\end{aligned} $$

The radial distortion model says that coordinates in the observed images are displaced towards (barrel distortion) or away (pincushion distortion) from the image center by an amount proportional to their radial distance (Figure 2.13a–b).⁶ The simplest radial distortion models use low-order polynomials, e.g.,

$$ \begin{align*}\hat{x}_{c}&=x_{c}(1+\kappa_{1}r_{c}^{2}+\kappa_{2}r_{c}^{4})\\\hat{y}_{c}&=y_{c}(1+\kappa_{1}r_{c}^{2}+\kappa_{2}r_{c}^{4}),\end{align*} $$

where $ r_c^2 = x_c^2 + y_c^2 $ and $ \kappa_1 $ and $ \kappa_2 $ are called the radial distortion parameters. $ ^7 $ This model, which also includes a tangential component to account for lens decentering, was first proposed in the photogrammetry literature by Brown (1966), and so is sometimes called the Brown or Brown-Conrady model. However, the tangential components of the distortion are usually ignored because they can lead to less stable estimates (Zhang 2000).

After the radial distortion step, the final pixel coordinates can be computed using

$$ \begin{aligned}\boldsymbol{x}_{s}&=f\hat{\boldsymbol{x}}_{c}+c_{x}\\y_{s}&=f\hat{y}_{c}+c_{y}.\end{aligned} $$

A variety of techniques can be used to estimate the radial distortion parameters for a given lens, as discussed in Section 11.1.4.

Sometimes the above simplified model does not model the true distortions produced by complex lenses accurately enough (especially at very wide angles). A more complete analytic model also includes tangential distortions and decentering distortions (Slama 1980).

Fisheye lenses (Figure 2.13c) require a model that differs from traditional polynomial models of radial distortion. Fisheye lenses behave, to a first approximation, as equi-distance projectors of angles away from the optical axis (Xiong and Turkowski 1997),

$$ r=f\theta, $$

which is the same as the polar projection described by Equations (8.55–8.57). Because of the mostly linear mapping between distance from the center (pixels) and viewing angle, such lenses are sometimes called $f\text{-}\theta\text{-}\lambda\text{lenses}$, which is likely where the popular RICOH THETA $360^{\circ}$ camera got its name. Xiong and Turkowski (1997) describe how this model can be extended with the addition of an extra quadratic correction in $\phi$ and how the unknown parameters (center of projection, scaling factor $s$, etc.) can be estimated from a set of overlapping fisheye images using a direct (intensity-based) non-linear minimization algorithm.

For even larger, less regular distortions, a parametric distortion model using splines may be necessary (Goshtasby 1989). If the lens does not have a single center of projection, it may become necessary to model the 3D line (as opposed to direction) corresponding to each pixel separately (Gremban, Thorpe, and Kanade 1988; Champleboux, Lavallée et al. 1992a; Grossberg and Nayar 2001; Sturm and Ramalingam 2004; Tardif, Sturm et al. 2009). Some of these techniques are described in more detail in Section 11.1.4, which discusses how to calibrate lens distortions.

There is one subtle issue associated with the simple radial distortion model that is often glossed over. We have introduced a non-linearity between the perspective projection and final sensor array projection steps. Therefore, we cannot, in general, post-multiply an arbitrary $ 3 \times 3 $ matrix K with a rotation to put it into upper-triangular form and absorb this into the global rotation. However, this situation is not as bad as it may at first appear. For many applications, keeping the simplified diagonal form of (2.59) is still an adequate model. Furthermore, if we correct radial and other distortions to an accuracy where straight lines are preserved, we have essentially converted the sensor back into a linear imager and the previous decomposition still applies.

2.2 Photometric image formation

In modeling the image formation process, we have described how 3D geometric features in the world are projected into 2D features in an image. However, images are not composed of 2D features. Instead, they are made up of discrete color or intensity values. Where do these values come from? How do they relate to the lighting in the environment, surface properties and geometry, camera optics, and sensor properties (Figure 2.14)? In this section, we develop a set of models to describe these interactions and formulate a generative process of image formation. A more detailed treatment of these topics can be found in textbooks on computer graphics and image synthesis (Cohen and Wallace 1993; Sillion and Puech 1994; Watt 1995; Glassner 1995; Weyrich, Lawrence et al. 2009; Hughes, van Dam et al. 2013; Marschner and Shirley 2015).

2.2.1 Lighting

Images cannot exist without light. To produce an image, the scene must be illuminated with one or more light sources. (Certain modalities such as fluorescence microscopy and X-ray tomography do not fit this model, but we do not deal with them in this book.) Light sources can generally be divided into point and area light sources.

A point light source originates at a single location in space (e.g., a small light bulb), potentially at infinity (e.g., the Sun). (Note that for some applications such as modeling soft shadows (penumbras), the Sun may have to be treated as an area light source.) In addition to its location, a point light source has an intensity and a color spectrum, i.e., a distribution over wavelengths $ L(\lambda) $. The intensity of a light source falls off with the square of the distance between the source and the object being lit, because the same light is being spread over a larger (spherical) area. A light source may also have a directional falloff (dependence), but we ignore this in our simplified model.

Area light sources are more complicated. A simple area light source such as a fluorescent ceiling light fixture with a diffuser can be modeled as a finite rectangular area emitting light equally in all directions (Cohen and Wallace 1993; Sillion and Puech 1994; Glassner 1995). When the distribution is strongly directional, a four-dimensional lightfield can be used instead (Ashdown 1993).

A more complex light distribution that approximates, say, the incident illumination on an object

sitting in an outdoor courtyard, can often be represented using an environment map (Greene 1986) (originally called a reflection map (Blinn and Newell 1976)). This representation maps incident light directions $ \hat{v} $ to color values (or wavelengths, $ \lambda $),

$$ L(\hat{\mathbf{v}};\lambda), $$

and is equivalent to assuming that all light sources are at infinity. Environment maps can be represented as a collection of cubical faces (Greene 1986), as a single longitude–latitude map (Blinn and Newell 1976), or as the image of a reflecting sphere (Watt 1995). A convenient way to get a rough model of a real-world environment map is to take an image of a reflective mirrored sphere (sometimes accompanied by a darker sphere to capture highlights) and to unwrap this image onto the desired environment map (Debevec 1998). Watt (1995) gives a nice discussion of environment mapping, including the formulas needed to map directions to pixels for the three most commonly used representations.

2.2.2 Reflectance and shading

When light hits an object's surface, it is scattered and reflected (Figure 2.15a). Many different models have been developed to describe this interaction. In this section, we first describe the most general form, the bidirectional reflectance distribution function, and then look at some more specialized models, including the diffuse, specular, and Phong shading models. We also discuss how these models can be used to compute the global illumination corresponding to a scene.

The Bidirectional Reflectance Distribution Function (BRDF)

The most general model of light scattering is the bidirectional reflectance distribution function (BRDF). $ ^{8} $ Relative to some local coordinate frame on the surface, the BRDF is a four-dimensional function that describes how much of each wavelength arriving at an incident direction $ \hat{v}_{i} $ is emitted in a reflected direction $ \hat{v}_{r} $ (Figure 2.15b). The function can be written in terms of the angles of the incident and reflected directions relative to the surface frame as

$$ f_{r}(\theta_{i},\phi_{i},\theta_{r},\phi_{r};\lambda). $$

The BRDF is reciprocal, i.e., because of the physics of light transport, you can interchange the roles of $ \hat{v}_i $ and $ \hat{v}_r $ and still get the same answer (this is sometimes called Helmholtz reciprocity).

Most surfaces are isotropic, i.e., there are no preferred directions on the surface as far as light transport is concerned. (The exceptions are anisotropic surfaces such as brushed (scratched) aluminum, where the reflectance depends on the light orientation relative to the direction of the scratches.) For an isotropic material, we can simplify the BRDF to

$$ f_{r}(\theta_{i},\theta_{r},|\phi_{r}-\phi_{i}|;\lambda)\quad\mathrm{o r}\quad f_{r}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{r},\hat{\mathbf{n}};\lambda), $$

as the quantities $ \theta_i, \theta_r $, and $ \phi_r - \phi_i $ can be computed from the directions $ \hat{v}_i, \hat{v}_r $, and $ \hat{n} $.

To calculate the amount of light exiting a surface point p in a direction $ \hat{v}_r $ under a given lighting condition, we integrate the product of the incoming light $ L_i(\hat{v}_i; \lambda) $ with the BRDF (some authors call this step a convolution). Taking into account the foreshortening factor $ \cos^+\theta_i $, we obtain

$$ L_{r}(\hat{\mathbf{v}}_{r};\lambda)=\int L_{i}(\hat{\mathbf{v}}_{i};\lambda)f_{r}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{r},\hat{\mathbf{n}};\lambda)\cos^{+}\theta_{i}d\hat{\mathbf{v}}_{i}, $$

where

$$ \cos^{+}\theta_{i}=\max(0,\cos\theta_{i}). $$

If the light sources are discrete (a finite number of point light sources), we can replace the integral with a summation.

$$ L_{r}(\hat{\mathbf{v}}_{r};\lambda)=\sum_{i}L_{i}(\lambda)f_{r}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{r},\hat{\mathbf{n}};\lambda)\cos^{+}\theta_{i}. $$

BRDFs for a given surface can be obtained through physical modeling (Torrance and Sparrow 1967; Cook and Torrance 1982; Glassner 1995), heuristic modeling (Phong 1975; Lafortune, Foo et al. 1997), or through empirical observation (Ward 1992; Westin, Arvo, and Torrance 1992; Dana, van Ginneken et al. 1999; Marschner, Westin et al. 2000; Matusik, Pfister et al. 2003; Dorsey, Rushmeier, and Sillion 2007; Weyrich, Lawrence et al. 2009; Shi, Mo et al. 2019).⁹ Typical BRDFs can often be split into their diffuse and specular components, as described below.

Diffuse reflection

The diffuse component (also known as Lambertian or matte reflection) scatters light uniformly in all directions and is the phenomenon we most normally associate with shading, e.g., the smooth (non-shiny) variation of intensity with surface normal that is seen when observing a statue (Figure 2.16). Diffuse reflection also often imparts a strong body color to the light, as it is caused by selective absorption and re-emission of light inside the object's material (Shafer 1985; Glassner 1995).

While light is scattered uniformly in all directions, i.e., the BRDF is constant.

$$ f_{d}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{r},\hat{\mathbf{n}};\lambda)=f_{d}(\lambda), $$

the amount of light depends on the angle between the incident light direction and the surface normal $ \theta_{i} $. This is because the surface area exposed to a given amount of light becomes larger at oblique angles, becoming completely self-shadowed as the outgoing surface normal points away from the light (Figure 2.17a). (Think about how you orient yourself towards the Sun or fireplace to get maximum warmth and how a flashlight projected obliquely against a wall is less bright than one pointing directly at it.) The shading equation for diffuse reflection can thus be written as

$$ L_{d}(\hat{\mathbf{v}}_{r};\lambda)=\sum_{i}L_{i}(\lambda)f_{d}(\lambda)\cos^{+}\theta_{i}=\sum_{i}L_{i}(\lambda)f_{d}(\lambda)[\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{n}}]^{+}, $$

where

$$ [\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{n}}]^{+}=\operatorname*{m a x}(0,\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{n}}). $$

Specular reflection

The second major component of a typical BRDF is specular (gloss or highlight) reflection, which depends strongly on the direction of the outgoing light. Consider light reflecting off a mirrored surface (Figure 2.17b). Incident light rays are reflected in a direction that is rotated by $ 180^{\circ} $ around the surface normal $ \hat{n} $. Using the same notation as in Equations (2.29–2.30), we can compute the specular reflection direction $ \hat{s}_{i} $ as

$$ \mathbf{\hat{s}}_{i}=\mathbf{v}_{\parallel}-\mathbf{v}_{\perp}=(2\mathbf{\hat{n}}\mathbf{\hat{n}}^{T}-\mathbf{I})\mathbf{v}_{i}. $$

2.2 Photometric image formation

The amount of light reflected in a given direction $ \hat{\mathbf{v}}_{r} $ thus depends on the angle $ \theta_{s} = \cos^{-1}(\hat{\mathbf{v}}_{r} \cdot \hat{\mathbf{s}}_{i}) $ between the view direction $ \hat{\mathbf{v}}_{r} $ and the specular direction $ \hat{\mathbf{s}}_{i} $. For example, the Phong (1975) model uses a power of the cosine of the angle,

$$ f_{s}(\theta_{s};\lambda)=k_{s}(\lambda)\cos^{k_{e}}\theta_{s}, $$

while the Torrance and Sparrow (1967) micro-facet model uses a Gaussian.

$$ f_{s}(\theta_{s};\lambda)=k_{s}(\lambda)\exp(-c_{s}^{2}\theta_{s}^{2}). $$

Larger exponents $ k_{e} $ (or inverse Gaussian widths $ c_{s} $) correspond to more specular surfaces with distinct highlights, while smaller exponents better model materials with softer gloss.

Phong shading

Phong (1975) combined the diffuse and specular components of reflection with another term, which he called the ambient illumination. This term accounts for the fact that objects are generally illuminated not only by point light sources but also by a general diffuse illumination corresponding to inter-reflection (e.g., the walls in a room) or distant sources, such as the blue sky. In the Phong model, the ambient term does not depend on surface orientation, but depends on the color of both the ambient illumination $ L_a(\lambda) $ and the object $ k_a(\lambda) $.

$$ f_{a}(\lambda)=k_{a}(\lambda)L_{a}(\lambda). $$

Putting all of these terms together, we arrive at the Phong shading model.

$$ L_{r}(\hat{\mathbf{v}}_{r};\lambda)=k_{a}(\lambda)L_{a}(\lambda)+k_{d}(\lambda)\sum_{i}L_{i}(\lambda)[\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{n}}]^{+}+k_{s}(\lambda)\sum_{i}L_{i}(\lambda)(\hat{\mathbf{v}}_{r}\cdot\hat{\mathbf{s}}_{i})^{k_{e}}. $$

Figure 2.18 shows a typical set of Phong shading model components as a function of the angle away from the surface normal (in a plane containing both the lighting direction and the viewer).

Typically, the ambient and diffuse reflection color distributions $ k_{a}(\lambda) $ and $ k_{d}(\lambda) $ are the same, since they are both due to sub-surface scattering (body reflection) inside the surface material (Shafer 1985). The specular reflection distribution $ k_{s}(\lambda) $ is often uniform (white), since it is caused by interface reflections that do not change the light color. (The exception to this is emphetallic materials, such as copper, as opposed to the more common dielectric materials, such as plastics.)

The ambient illumination $ L_a(\lambda) $ often has a different color cast from the direct light sources $ L_i(\lambda) $, e.g., it may be blue for a sunny outdoor scene or yellow for an interior lit with candles or incandescent lights. (The presence of ambient sky illumination in shadowed areas is what often causes shadows to appear bluer than the corresponding lit portions of a scene). Note also that the diffuse component of the Phong model (or of any shading model) depends on the angle of the incoming light source $ \hat{v}_i $, while the specular component depends on the relative angle between the viewer $ v_r $ and the specular reflection direction $ \hat{s}_i $ (which itself depends on the incoming light direction $ \hat{v}_i $ and the surface normal $ \hat{n} $).

The Phong shading model has been superseded in terms of physical accuracy by newer models in computer graphics, including the model developed by Cook and Torrance (1982) based on the original micro-facet model of Torrance and Sparrow (1967). While, initially, computer graphics hardware implemented the Phong model, the advent of programmable pixel shaders has made the use of more complex models feasible.

Di-chromatic reflection model

The Torrance and Sparrow (1967) model of reflection also forms the basis of Shafer's (1985) dichromatic reflection model, which states that the apparent color of a uniform material lit from a single source depends on the sum of two terms,

$$ L_{r}(\hat{\mathbf{v}}_{r};\lambda)=L_{i}(\hat{\mathbf{v}}_{r},\hat{\mathbf{v}}_{i},\hat{\mathbf{n}};\lambda)+L_{b}(\hat{\mathbf{v}}_{r},\hat{\mathbf{v}}_{i},\hat{\mathbf{n}};\lambda) $$

$$ =c_{i}(\lambda)m_{i}(\hat{\mathbf{v}}_{r},\hat{\mathbf{v}}_{i},\hat{\mathbf{n}})+c_{b}(\lambda)m_{b}(\hat{\mathbf{v}}_{r},\hat{\mathbf{v}}_{i},\hat{\mathbf{n}}), $$

i.e., the radiance of the light reflected at the interface, $ L_i $, and the radiance reflected at the surface body, $ L_b $. Each of these, in turn, is a simple product between a relative power spectrum $ c(\lambda) $, which depends only on wavelength, and a magnitude $ m(\hat{\mathbf{v}}_r, \hat{\mathbf{v}}_i, \hat{\mathbf{n}}) $, which depends only on geometry. (This model can easily be derived from a generalized version of Phong’s model by assuming a single light source and no ambient illumination, and rearranging terms.) The di-chromatic model has been successfully used in computer vision to segment specular colored objects with large variations in shading (Klinker 1993) and has inspired local two-color models for applications such as Bayer pattern demosaicing (Bennett, Uyttendaele et al. 2006).

Global illumination (ray tracing and radiosity)

The simple shading model presented thus far assumes that light rays leave the light sources, bounce off surfaces visible to the camera, thereby changing in intensity or color, and arrive at the camera. In reality, light sources can be shadowed by occluders and rays can bounce multiple times around a scene while making their trip from a light source to the camera.

Two methods have traditionally been used to model such effects. If the scene is mostly specular (the classic example being scenes made of glass objects and mirrored or highly polished balls), the preferred approach is ray tracing or path tracing (Glassner 1995; Akenine-Möller and Haines 2002; Marschner and Shirley 2015), which follows individual rays from the camera across multiple bounces towards the light sources (or vice versa). If the scene is composed mostly of uniform albedo simple geometry illuminators and surfaces, radiosity (global illumination) techniques are preferred (Cohen and Wallace 1993; Sillion and Puech 1994; Glassner 1995). Combinations of the two techniques have also been developed (Wallace, Cohen, and Greenberg 1987), as well as more general light transport techniques for simulating effects such as the caustics cast by rippling water.

The basic ray tracing algorithm associates a light ray with each pixel in the camera image and finds its intersection with the nearest surface. A primary contribution can then be computed using

the simple shading equations presented previously (e.g., Equation (2.94)) for all light sources that are visible for that surface element. (An alternative technique for computing which surfaces are illuminated by a light source is to compute a shadow map, or shadow buffer, i.e., a rendering of the scene from the light source's perspective, and then compare the depth of pixels being rendered with the map (Williams 1983; Akenine-Möller and Haines 2002).) Additional secondary rays can then be cast along the specular direction towards other objects in the scene, keeping track of any attenuation or color change that the specular reflection induces.

Radiosity works by associating lightness values with rectangular surface areas in the scene (including area light sources). The amount of light interchanged between any two (mutually visible) areas in the scene can be captured as a form factor, which depends on their relative orientation and surface reflectance properties, as well as the $ 1/r^{2} $ fall-off as light is distributed over a larger effective sphere the further away it is (Cohen and Wallace 1993; Sillion and Puech 1994; Glassner 1995). A large linear system can then be set up to solve for the final lightness of each area patch, using the light sources as the forcing function (right-hand side). Once the system has been solved, the scene can be rendered from any desired point of view. Under certain circumstances, it is possible to recover the global illumination in a scene from photographs using computer vision techniques (Yu, Debevec et al. 1999).

The basic radiosity algorithm does not take into account certain near field effects, such as the darkening inside corners and scratches, or the limited ambient illumination caused by partial shadowing from other surfaces. Such effects have been exploited in a number of computer vision algorithms (Nayar, Ikeuchi, and Kanade 1991; Langer and Zucker 1994).

While all of these global illumination effects can have a strong effect on the appearance of a scene, and hence its 3D interpretation, they are not covered in more detail in this book. (But see Section 13.7.1 for a discussion of recovering BRDFs from real scenes and objects.)

2.2.3 Optics

Once the light from a scene reaches the camera, it must still pass through the lens before reaching the analog or digital sensor. For many applications, it suffices to treat the lens as an ideal pinhole that simply projects all rays through a common center of projection (Figures 2.8 and 2.9).

However, if we want to deal with issues such as focus, exposure, vignetting, and aberration, we need to develop a more sophisticated model, which is where the study of optics comes in (Möller 1988; Ray 2002; Hecht 2015).

Figure 2.19 shows a diagram of the most basic lens model, i.e., the thin lens composed of a single piece of glass with very low, equal curvature on both sides. According to the lens law (which can be derived using simple geometric arguments on light ray refraction), the relationship between the distance to an object $ z_{o} $ and the distance behind the lens at which a focused image is formed $ z_{i} $ can be expressed as

$$ \frac{1}{z_{o}}+\frac{1}{z_{i}}=\frac{1}{f}\ , $$

where $f$ is called the focal length of the lens. If we let $z_o \to \infty$, i.e., we adjust the lens (move the image plane) so that objects at infinity are in focus, we get $z_i = f$, which is why we can think of a lens of focal length $f$ as being equivalent (to a first approximation) to a pinhole at a distance $f$ from the focal plane (Figure 2.10), whose field of view is given by (2.60).

If the focal plane is moved away from its proper in-focus setting of $ z_i $ (e.g., by twisting the focus ring on the lens), objects at $ z_o $ are no longer in focus, as shown by the gray plane in Figure 2.19. The amount of misfocus is measured by the circle of confusion $ c $ (shown as short thick blue line segments on the gray plane). $ ^{10} $ The equation for the circle of confusion can be derived using similar triangles; it depends on the distance of travel in the focal plane $ \Delta z_i $ relative to the original focus distance $ z_i $ and the diameter of the aperture $ d $ (see Exercise 2.4).

The allowable depth variation in the scene that limits the circle of confusion to an acceptable number is commonly called the depth of field and is a function of both the focus distance and the aperture, as shown diagrammatically by many lens markings (Figure 2.20). Since this depth of field depends on the aperture diameter $d$, we also have to know how this varies with the commonly displayed $f$-number, which is usually denoted as $f/\#$ or $N$ and is defined as

$$ f/\#\equiv N=\frac{f}{d}\mathrm{~,~} $$

where the focal length $f$ and the aperture diameter $d$ are measured in the same unit (say, millimeters).

The usual way to write the f-number is to replace the # in $f/\#$ with the actual number, i.e., $f/1.4$, $f/2$, $f/2.8$, ..., $f/22$. (Alternatively, we can say $N = 1.4$, etc.) An easy way to interpret

these numbers is to notice that dividing the focal length by the f-number gives us the diameter d, so these are just formulas for the aperture diameter. $ ^{11} $

Notice that the usual progression for f-numbers is in full stops, which are multiples of $ \sqrt{2} $, since this corresponds to doubling the area of the entrance pupil each time a smaller f-number is selected. (This doubling is also called changing the exposure by one exposure value or EV. It has the same effect on the amount of light reaching the sensor as doubling the exposure duration, e.g., from $ ^{1}/_{250} $ to $ ^{1}/_{125} $; see Exercise 2.5.)

Now that you know how to convert between f-numbers and aperture diameters, you can construct your own plots for the depth of field as a function of focal length $f$, circle of confusion $c$, and focus distance $z_{o}$, as explained in Exercise 2.4, and see how well these match what you observe on actual lenses, such as those shown in Figure 2.20.

Of course, real lenses are not infinitely thin and therefore suffer from geometric aberrations, unless compound elements are used to correct for them. The classic five Seidel aberrations, which arise when using third-order optics, include spherical aberration, coma, astigmatism, curvature of field, and distortion (Möller 1988; Ray 2002; Hecht 2015).

Chromatic aberration

Because the index of refraction for glass varies slightly as a function of wavelength, simple lenses suffer from chromatic aberration, which is the tendency for light of different colors to focus at slightly different distances (and hence also with slightly different magnification factors), as shown in Figure 2.21. The wavelength-dependent magnification factor, i.e., the transverse chromatic aberration, can be modeled as a per-color radial distortion (Section 2.1.5) and, hence, calibrated using the techniques described in Section 11.1.4. The wavelength-dependent blur caused by longitudinal chromatic aberration can be calibrated using techniques described in Section 10.1.4. Unfortunately, the blur induced by longitudinal aberration can be harder to undo, as higher frequencies can get strongly attenuated and hence hard to recover.

To reduce chromatic and other kinds of aberrations, most photographic lenses today are compound lenses made of different glass elements (with different coatings). Such lenses can no longer be modeled as having a single nodal point $P$ through which all of the rays must pass (when approximating the lens with a pinhole model). Instead, these lenses have both a front nodal point, through which the rays enter the lens, and a rear nodal point, through which they leave on their way to the

sensor. In practice, only the location of the front nodal point is of interest when performing careful camera calibration, e.g., when determining the point around which to rotate to capture a parallax-free panorama (see Section 8.2.3 and Littlefield (2006) and Houghton (2013)).

Not all lenses, however, can be modeled as having a single nodal point. In particular, very wide-angle lenses such as fisheye lenses (Section 2.1.5) and certain catadioptric imaging systems consisting of lenses and curved mirrors (Baker and Nayar 1999) do not have a single point through which all of the acquired light rays pass. In such cases, it is preferable to explicitly construct a mapping function (look-up table) between pixel coordinates and 3D rays in space (Gremban, Thorpe, and Kanade 1988; Champleboux, Lavallée et al. 1992a; Grossberg and Nayar 2001; Sturm and Ramalingam 2004; Tardif, Sturm et al. 2009), as mentioned in Section 2.1.5.

Vignetting

Another property of real-world lenses is vignetting, which is the tendency for the brightness of the image to fall off towards the edge of the image.

Two kinds of phenomena usually contribute to this effect (Ray 2002). The first is called natural vignetting and is due to the foreshortening in the object surface, projected pixel, and lens aperture, as shown in Figure 2.22. Consider the light leaving the object surface patch of size $ \delta o $ located at an off-axis angle $ \alpha $. Because this patch is foreshortened with respect to the camera lens, the amount of light reaching the lens is reduced by a factor $ \cos \alpha $. The amount of light reaching the lens is also subject to the usual $ 1/r^2 $ fall-off; in this case, the distance $ r_o = z_o / \cos \alpha $. The actual area of the aperture through which the light passes is foreshortened by an additional factor $ \cos \alpha $, i.e., the aperture as seen from point $ O $ is an ellipse of dimensions $ d \times d \cos \alpha $. Putting all of these factors together, we see that the amount of light leaving $ O $ and passing through the aperture on its way to the image pixel located at $ I $ is proportional to

$$ \frac{\delta o\cos\alpha}{r_{o}^{2}}\pi\left(\frac{d}{2}\right)^{2}\cos\alpha=\delta o\frac{\pi}{4}\frac{d^{2}}{z_{o}^{2}}\cos^{4}\alpha. $$

Since triangles $ \Delta OPQ $ and $ \Delta IPJ $ are similar, the projected areas of the object surface $ \delta o $ and image pixel $ \delta i $ are in the same (squared) ratio as $ z_o : z_i $,

$$ \frac{\delta o}{\delta i}=\frac{z_{o}^{2}}{z_{i}^{2}}. $$

2.3 The digital camera

Putting these together, we obtain the final relationship between the amount of light reaching pixel i and the aperture diameter d, the focusing distance $ z_i \approx f $, and the off-axis angle $ \alpha $,

$$ \delta o\frac{\pi}{4}\frac{d^{2}}{z_{o}^{2}}\cos^{4}\alpha=\delta i\frac{\pi}{4}\frac{d^{2}}{z_{i}^{2}}\cos^{4}\alpha\approx\delta i\frac{\pi}{4}\left(\frac{d}{f}\right)^{2}\cos^{4}\alpha, $$

which is called the fundamental radiometric relation between the scene radiance L and the light (irradiance) E reaching the pixel sensor,

$$ E=L\frac{\pi}{4}\left(\frac{d}{f}\right)^{2}\cos^{4}\alpha, $$

(Horn 1986; Nalwa 1993; Ray 2002; Hecht 2015). Notice in this equation how the amount of light depends on the pixel surface area (which is why the smaller sensors in point-and-shoot cameras are so much noisier than digital single lens reflex (SLR) cameras), the inverse square of the f-stop $ N = f/d $ (2.98), and the fourth power of the $ \cos^4 \alpha $ off-axis fall-off, which is the natural vignetting term.

The other major kind of vignetting, called mechanical vignetting, is caused by the internal occlusion of rays near the periphery of lens elements in a compound lens, and cannot easily be described mathematically without performing a full ray-tracing of the actual lens design. $ ^{12} $ However, unlike natural vignetting, mechanical vignetting can be decreased by reducing the camera aperture (increasing the f-number). It can also be calibrated (along with natural vignetting) using special devices such as integrating spheres, uniformly illuminated targets, or camera rotation, as discussed in Section 10.1.3.

2.3 The digital camera

After starting from one or more light sources, reflecting off one or more surfaces in the world, and passing through the camera's optics (lenses), light finally reaches the imaging sensor. How are the photons arriving at this sensor converted into the digital (R, G, B) values that we observe when we look at a digital image? In this section, we develop a simple model that accounts for the most important effects, such as exposure (gain and shutter speed), non-linear mappings, sampling and aliasing, and noise. Figure 2.23, which is based on camera models developed by Healey and Kondepudy (1994), Tsin, Ramesh, and Kanade (2001), and Liu, Szeliski et al. (2008), shows a simple version of the processing stages that occur in modern digital cameras. Chakrabarti, Scharstein, and Zickler (2009) developed a sophisticated 24-parameter model that is an even better match to the processing performed in digital cameras, while Kim, Lin et al. (2012), Hasinoff, Sharlet et al. (2016), and Karaimer and Brown (2016) provide more recent models of modern in-camera processing pipelines. Most recently, Brooks, Mildenhall et al. (2019) have developed detailed models of in-camera image processing pipelines to invert (unprocess) noisy JPEG images into their RAW originals, so that they can be better denoised, while Tseng, Yu et al. (2019) develop a tunable model of camera processing pipelines that can be used for image quality optimization.

Light falling on an imaging sensor is usually picked up by an active sensing area, integrated for the duration of the exposure (usually expressed as the shutter speed in a fraction of a second, e.g., $ \frac{1}{125} $, $ \frac{1}{60} $, $ \frac{1}{30} $), and then passed to a set of sense amplifiers. The two main kinds of sensor used in digital still and video cameras today are charge-coupled device (CCD) and complementary metal oxide on silicon (CMOS).

In a CCD, photons are accumulated in each active well during the exposure time. Then, in a transfer phase, the charges are transferred from well to well in a kind of “bucket brigade” until they are deposited at the sense amplifiers, which amplify the signal and pass it to an analog-to-digital converter (ADC). $ ^{13} $ Older CCD sensors were prone to blooming, when charges from one over-exposed pixel spilled into adjacent ones, but most newer CCDs have anti-blooming technology (“troughs” into which the excess charge can spill).

In CMOS, the photons hitting the sensor directly affect the conductivity (or gain) of a photodetector, which can be selectively gated to control exposure duration, and locally amplified before being read out using a multiplexing scheme. Traditionally, CCD sensors outperformed CMOS in quality-sensitive applications, such as digital SLRs, while CMOS was better for low-power applications, but today CMOS is used in most digital cameras.

The main factors affecting the performance of a digital image sensor are the shutter speed, sampling pitch, fill factor, chip size, analog gain, sensor noise, and the resolution (and quality) of the analog-to-digital converter. Many of the actual values for these parameters can be read from the EXIF tags embedded with digital images, while others can be obtained from the camera manufacturers' specification sheets or from camera review or calibration websites. $ ^{14} $

Shutter speed. The shutter speed (exposure time) directly controls the amount of light reaching the sensor and hence determines if images are under- or over-exposed. (For bright scenes, where a large aperture or slow shutter speed is desired to get a shallow depth of field or motion blur, neutral density filters are sometimes used by photographers.) For dynamic scenes, the shutter speed

also determines the amount of motion blur in the resulting picture. Usually, a higher shutter speed (less motion blur) makes subsequent analysis easier (see Section 10.3 for techniques to remove such blur). However, when video is being captured for display, some motion blur may be desirable to avoid stroboscopic effects.

Sampling pitch. The sampling pitch is the physical spacing between adjacent sensor cells on the imaging chip (Figure 2.24). A sensor with a smaller sampling pitch has a higher sampling density and hence provides a higher resolution (in terms of pixels) for a given active chip area. However, a smaller pitch also means that each sensor has a smaller area and cannot accumulate as many photons; this makes it not as light sensitive and more prone to noise.

Fill factor. The fill factor is the active sensing area size as a fraction of the theoretically available sensing area (the product of the horizontal and vertical sampling pitches). Higher fill factors are usually preferable, as they result in more light capture and less aliasing (see Section 2.3.1). While the fill factor was originally limited by the need to place additional electronics between the active sensing areas, modern backside illumination (or back-illuminated) sensors, coupled with efficient microlens designs, have largely removed this limitation (Fontaine 2015).$^{15}$ The fill factor of a camera can be determined empirically using a photometric camera calibration process (see Section 10.1.4).

Chip size. Video and point-and-shoot cameras have traditionally used small chip areas ( $ \frac{1}{4} $-inch to $ \frac{1}{2} $-inch sensors $ ^{16} $), while digital SLR cameras try to come closer to the traditional size of a 35mm film frame. $ ^{17} $ When overall device size is not important, having a larger chip size is preferable, since each

sensor cell can be more photo-sensitive. (For compact cameras, a smaller chip means that all of the optics can be shrunk down proportionately.) However, larger chips are more expensive to produce, not only because fewer chips can be packed into each wafer, but also because the probability of a chip defect goes up exponentially with the chip area.

Analog gain. Before analog-to-digital conversion, the sensed signal is usually boosted by a sense amplifier. In video cameras, the gain on these amplifiers was traditionally controlled by automatic gain control (AGC) logic, which would adjust these values to obtain a good overall exposure. In newer digital still cameras, the user now has some additional control over this gain through the ISO setting, which is typically expressed in ISO standard units such as 100, 200, or 400. Since the automated exposure control in most cameras also adjusts the aperture and shutter speed, setting the ISO manually removes one degree of freedom from the camera's control, just as manually specifying aperture and shutter speed does. In theory, a higher gain allows the camera to perform better under low light conditions (less motion blur due to long exposure times when the aperture is already maxed out). In practice, however, higher ISO settings usually amplify the sensor noise.

Sensor noise. Throughout the whole sensing process, noise is added from various sources, which may include fixed pattern noise, dark current noise, shot noise, amplifier noise, and quantization noise (Healey and Kondepudy 1994; Tsin, Ramesh, and Kanade 2001). The final amount of noise present in a sampled image depends on all of these quantities, as well as the incoming light (controlled by the scene radiance and aperture), the exposure time, and the sensor gain. Also, for low light conditions where the noise is due to low photon counts, a Poisson model of noise may be more appropriate than a Gaussian model (Alter, Matsushita, and Tang 2006; Matsushita and Lin 2007a; Wilburn, Xu, and Matsushita 2008; Takamatsu, Matsushita, and Ikeuchi 2008).

As discussed in more detail in Section 10.1.1, Liu, Szeliski et al. (2008) use this model, along with an empirical database of camera response functions (CRFs) obtained by Grossberg and Nayar (2004), to estimate the noise level function (NLF) for a given image, which predicts the overall noise variance at a given pixel as a function of its brightness (a separate NLF is estimated for each color channel). An alternative approach, when you have access to the camera before taking pictures, is to pre-calibrate the NLF by taking repeated shots of a scene containing a variety of colors and luminances, such as the Macbeth Color Chart shown in Figure 10.3b (McCamy, Marcus, and Davidson 1976). (When estimating the variance, be sure to throw away or downweight pixels with large gradients, as small shifts between exposures will affect the sensed values at such pixels.) Unfortunately, the pre-calibration process may have to be repeated for different exposure times and gain settings because of the complex interactions occurring within the sensing system.

In practice, most computer vision algorithms, such as image denoising, edge detection, and stereo matching, all benefit from at least a rudimentary estimate of the noise level. Barring the ability to pre-calibrate the camera or to take repeated shots of the same scene, the simplest approach is to look for regions of near-constant value and to estimate the noise variance in such regions (Liu, Szeliski et al. 2008).

ADC resolution. The final step in the analog processing chain occurring within an imaging sensor is the analog to digital conversion (ADC). While a variety of techniques can be used to implement this process, the two quantities of interest are the resolution of this process (how many bits it yields) and its noise level (how many of these bits are useful in practice). For most cameras, the number of bits quoted (eight bits for compressed JPEG images and a nominal 16 bits for the RAW formats provided by some DSLRs) exceeds the actual number of usable bits. The best way to

2.3 The digital camera

tell is to simply calibrate the noise of a given sensor, e.g., by taking repeated shots of the same scene and plotting the estimated noise as a function of brightness (Exercise 2.6).

Digital post-processing. Once the irradiance values arriving at the sensor have been converted to digital bits, most cameras perform a variety of digital signal processing (DSP) operations to enhance the image before compressing and storing the pixel values. These include color filter array (CFA) demosaicing, white point setting, and mapping of the luminance values through a gamma function to increase the perceived dynamic range of the signal. We cover these topics in Section 2.3.2 but, before we do, we return to the topic of aliasing, which was mentioned in connection with sensor array fill factors.

Newer imaging sensors. The capabilities of imaging sensor and related technologies such as depth sensors continue to evolve rapidly. Conferences that track these developments include the IS&T Symposium on Electronic Imaging Science and Technology sponsored by the Society for Imaging Science and Technology and the Image Sensors World blog.

2.3.1 Sampling and aliasing

What happens when a field of light impinging on the image sensor falls onto the active sense areas in the imaging chip? The photons arriving at each active cell are integrated and then digitized, as shown in Figure 2.24. However, if the fill factor on the chip is small and the signal is not otherwise band-limited, visually unpleasing aliasing can occur.

To explore the phenomenon of aliasing, let us first look at a one-dimensional signal (Figure 2.25), in which we have two sine waves, one at a frequency of $ f = \frac{3}{4} $ and the other at $ f = \frac{5}{4} $. If we sample these two signals at a frequency of $ f = 2 $, we see that they produce the same samples (shown in black), and so we say that they are $ aliased.^{18} $ Why is this a bad effect? In essence, we can no longer reconstruct the original signal, since we do not know which of the two original frequencies was present.

In fact, Shannon's Sampling Theorem shows that the minimum sampling (Oppenheim and Schafer 1996; Oppenheim, Schafer, and Buck 1999) rate required to reconstruct a signal from its

instantaneous samples must be at least twice the highest frequency, $ ^{19} $

$$ f_{s}\geq2f_{max}. $$

The maximum frequency in a signal is known as the Nyquist frequency and the inverse of the minimum sampling frequency $ r_{s}=1/f_{s} $ is known as the Nyquist rate.

However, you may ask, as an imaging chip actually averages the light field over a finite area, are the results on point sampling still applicable? Averaging over the sensor area does tend to attenuate some of the higher frequencies. However, even if the fill factor is 100%, as in the right image of Figure 2.25, frequencies above the Nyquist limit (half the sampling frequency) still produce an aliased signal, although with a smaller magnitude than the corresponding band-limited signals.

A more convincing argument as to why aliasing is bad can be seen by downsampling a signal using a poor quality filter such as a box (square) filter. Figure 2.26 shows a high-frequency chirp image (so called because the frequencies increase over time), along with the results of sampling it with a 25% fill-factor area sensor, a 100% fill-factor sensor, and a high-quality 9-tap filter. Additional examples of downsampling (decimation) filters can be found in Section 3.5.2 and Figure 3.29.

The best way to predict the amount of aliasing that an imaging system (or even an image processing algorithm) will produce is to estimate the point spread function (PSF), which represents the response of a particular pixel sensor to an ideal point light source. The PSF is a combination (convolution) of the blur induced by the optical system (lens) and the finite integration area of a chip sensor. $ ^{20} $

If we know the blur function of the lens and the fill factor (sensor area shape and spacing) for the imaging chip (plus, optionally, the response of the anti-aliasing filter), we can convolve these (as described in Section 3.2) to obtain the PSF. Figure 2.27a shows the one-dimensional cross-section of a PSF for a lens whose blur function is assumed to be a disc with a radius equal to the pixel spacing s plus a sensing chip whose horizontal fill factor is 80%. Taking the Fourier transform of this PSF (Section 3.4), we obtain the modulation transfer function (MTF), from which we can estimate the amount of aliasing as the area of the Fourier magnitude outside the $ f \leq f_s $ Nyquist frequency. $ ^{21} $

2.3 The digital camera

If we defocus the lens so that the blur function has a radius of 2s (Figure 2.27c), we see that the amount of aliasing decreases significantly, but so does the amount of image detail (frequencies closer to $ f = f_{s} $).

Under laboratory conditions, the PSF can be estimated (to pixel precision) by looking at a point light source such as a pinhole in a black piece of cardboard lit from behind. However, this PSF (the actual image of the pinhole) is only accurate to a pixel resolution and, while it can model larger blur (such as blur caused by defocus), it cannot model the sub-pixel shape of the PSF and predict the amount of aliasing. An alternative technique, described in Section 10.1.4, is to look at a calibration pattern (e.g., one consisting of slanted step edges (Reichenbach, Park, and Narayanswamy 1991; Williams and Burns 2001; Joshi, Szeliski, and Kriegman 2008)) whose ideal appearance can be re-synthesized to sub-pixel precision.

In addition to occurring during image acquisition, aliasing can also be introduced in various image processing operations, such as resampling, upsampling, and downsampling. Sections 3.4 and 3.5.2 discuss these issues and show how careful selection of filters can reduce the amount of aliasing.

2.3.2 Color

In Section 2.2, we saw how lighting and surface reflections are functions of wavelength. When the incoming light hits the imaging sensor, light from different parts of the spectrum is somehow

integrated into the discrete red, green, and blue (RGB) color values that we see in a digital image. How does this process work and how can we analyze and manipulate color values?

You probably recall from your childhood days the magical process of mixing paint colors to obtain new ones. You may recall that blue+yellow makes green, red+blue makes purple, and red+green makes brown. If you revisited this topic at a later age, you may have learned that the proper subtractive primaries are actually cyan (a light blue-green), magenta (pink), and yellow (Figure 2.28b), although black is also often used in four-color printing (CMYK). $ ^{22} $ If you ever subsequently took any painting classes, you learned that colors can have even more fanciful names, such as alizarin crimson, cerulean blue, and chartreuse. The subtractive colors are called subtractive because pigments in the paint absorb certain wavelengths in the color spectrum.

Later on, you may have learned about the additive primary colors (red, green, and blue) and how they can be added (with a slide projector or on a computer monitor) to produce cyan, magenta, yellow, white, and all the other colors we typically see on our TV sets and monitors (Figure 2.28a).

Through what process is it possible for two different colors, such as red and green, to interact to produce a third color like yellow? Are the wavelengths somehow mixed up to produce a new wavelength?

You probably know that the correct answer has nothing to do with physically mixing wavelengths. Instead, the existence of three primaries is a result of the tri-stimulus (or tri-chromatic) nature of the human visual system, since we have three different kinds of cells called cones, each of which responds selectively to a different portion of the color spectrum (Glassner 1995; Wandell 1995; Wyszecki and Stiles 2000; Livingstone 2008; Frisby and Stone 2010; Reinhard, Heidrich et al. 2010; Fairchild 2013). Note that for machine vision applications, such as remote sensing and terrain classification, it is preferable to use many more wavelengths. Similarly, surveillance applications can often benefit from sensing in the near-infrared (NIR) range.

CIE RGB and XYZ

To test and quantify the tri-chromatic theory of perception, we can attempt to reproduce all monochromatic (single wavelength) colors as a mixture of three suitably chosen primaries. (Pure wavelength light can be obtained using either a prism or specially manufactured color filters.) In the

2.3 The digital camera

1930s, the Commission Internationale d'Eclairage (CIE) standardized the RGB representation by performing such color matching experiments using the primary colors of red (700.0nm wavelength), green (546.1nm), and blue (435.8nm).

Figure 2.29 shows the results of performing these experiments with a standard observer, i.e., averaging perceptual results over a large number of subjects. $ ^{24} $ You will notice that for certain pure spectra in the blue–green range, a negative amount of red light has to be added, i.e., a certain amount of red has to be added to the color being matched to get a color match. These results also provided a simple explanation for the existence of metamers, which are colors with different spectra that are perceptually indistinguishable. Note that two fabrics or paint colors that are metamers under one light may no longer be so under different lighting.

Because of the problem associated with mixing negative light, the CIE also developed a new color space called XYZ, which contains all of the pure spectral colors within its positive octant. (It also maps the Y axis to the luminance, i.e., perceived relative brightness, and maps pure white to a diagonal (equal-valued) vector.) The transformation from RGB to XYZ is given by

$$ \begin{aligned}\begin{bmatrix}X\\Y\\Z\end{bmatrix}&=\frac{1}{0.17697}\begin{bmatrix}0.49&0.31&0.20\\0.17697&0.81240&0.01063\\0.00&0.01&0.99\end{bmatrix}\begin{bmatrix}R\\G\\B\end{bmatrix}.\end{aligned} $$

While the official definition of the CIE XYZ standard has the matrix normalized so that the Y value corresponding to pure red is 1, a more commonly used form is to omit the leading fraction, so that the second row adds up to one, i.e., the RGB triplet $ (1,1,1) $ maps to a Y value of 1. Linearly blending the $ (\bar{r}(\lambda),\bar{g}(\lambda),\bar{b}(\lambda)) $ curves in Figure 2.29a according to (2.104), we obtain the resulting $ (\bar{x}(\lambda),\bar{y}(\lambda),\bar{z}(\lambda)) $ curves shown in Figure 2.29b. Notice how all three spectra (color matching functions) now have only positive values and how the $ \bar{y}(\lambda) $ curve matches that of the luminance perceived by humans.

If we divide the XYZ values by the sum of $ X+Y+Z $, we obtain the chromaticity coordinates

$$ x=\frac{X}{X+Y+Z},\quad y=\frac{Y}{X+Y+Z},\quad z=\frac{Z}{X+Y+Z}, $$

which sum to 1. The chromaticity coordinates discard the absolute intensity of a given color sample and just represent its pure color. If we sweep the monochromatic color $\lambda$ parameter in Figure 2.29b from $\lambda = 380\mathrm{nm}$ to $\lambda = 800\mathrm{nm}$, we obtain the familiar chromaticity diagram shown in Figure 2.30a. This figure shows the $(x,y)$ value for every color value perceivable by most humans. (Of course, the CMYK reproduction process in this book does not actually span the whole gamut of perceivable colors.) The outer curved rim represents where all of the pure monochromatic color values map in $(x,y)$ space, while the lower straight line, which connects the two endpoints, is known as the purple line. The inset triangle spans the red, green, and blue single-wavelength primaries used in the original color matching experiments, while E denotes the white point.

A convenient representation for color values, when we want to tease apart luminance and chromaticity, is therefore Yxy (luminance plus the two most distinctive chrominance components).

L*a*b* color space

While the XYZ color space has many convenient properties, including the ability to separate luminance from chrominance, it does not actually predict how well humans perceive differences in color or luminance.

Because the response of the human visual system is roughly logarithmic (we can perceive relative luminance differences of about 1%), the CIE defined a non-linear re-mapping of the XYZ space called L*a*b* (also sometimes called CIELAB), where differences in luminance or chrominance are more perceptually uniform, as shown in Figure 2.30b.²⁵

2.3 The digital camera

The L* component of lightness is defined as

$$ L^{*}=116f\left(\frac{Y}{Y_{n}}\right), $$

where $ Y_{n} $ is the luminance value for nominal white (Fairchild 2013) and

$$ f(t)=\left\{\begin{array}{ll}t^{1/3}&t>\delta^{3}\\t/(3\delta^{2})+2\delta/3&else,\end{array}\right. $$

is a finite-slope approximation to the cube root with $ \delta = 6/29 $. The resulting 0...100 scale roughly measures equal amounts of lightness perceptibility.

In a similar fashion, the a* and b* components are defined as

$$ a^{*}=500\left[f\left(\frac{X}{X_{n}}\right)-f\left(\frac{Y}{Y_{n}}\right)\right]\quad\mathrm{a n d}\quad b^{*}=200\left[f\left(\frac{Y}{Y_{n}}\right)-f\left(\frac{Z}{Z_{n}}\right)\right], $$

where again, $ (X_n, Y_n, Z_n) $ is the measured white point. Figure 2.33i–k show the $ L^*a^*b^* $ representation for a sample color image.

Color cameras

While the preceding discussion tells us how we can uniquely describe the perceived tri-stimulus description of any color (spectral distribution), it does not tell us how RGB still and video cameras actually work. Do they just measure the amount of light at the nominal wavelengths of red (700.0nm), green (546.1nm), and blue (435.8nm)? Do color monitors just emit exactly these wavelengths and, if so, how can they emit negative red light to reproduce colors in the cyan range?

In fact, the design of RGB video cameras has historically been based around the availability of colored phosphors that go into television sets. When standard-definition color television was invented (NTSC), a mapping was defined between the RGB values that would drive the three color guns in the cathode ray tube (CRT) and the XYZ values that unambiguously define perceived color (this standard was called ITU-R BT.601). With the advent of HDTV and newer monitors, a new standard called ITU-R BT.709 was created, which specifies the XYZ values of each of the color primaries,

$$ \begin{aligned}\begin{bmatrix}X\\Y\\Z\end{bmatrix}&=\begin{bmatrix}0.412453&0.357580&0.180423\\0.212671&0.715160&0.072169\\0.019334&0.119193&0.950227\end{bmatrix}\begin{bmatrix}R_{709}\\G_{709}\\B_{709}\end{bmatrix}.\end{aligned} $$

In practice, each color camera integrates light according to the spectral response function of its red, green, and blue sensors.

$$ \begin{aligned}R&=\int L(\lambda)S_{R}(\lambda)d\lambda,\\G&=\int L(\lambda)S_{G}(\lambda)d\lambda,\\B&=\int L(\lambda)S_{B}(\lambda)d\lambda,\end{aligned} $$

where $L(\lambda)$ is the incoming spectrum of light at a given pixel and $\{S_{R}(\lambda), S_{G}(\lambda), S_{B}(\lambda)\}$ are the red, green, and blue spectral sensitivities of the corresponding sensors.

Can we tell what spectral sensitivities the cameras actually have? Unless the camera manufacturer provides us with these data or we observe the response of the camera to a whole spectrum of

monochromatic lights, these sensitivities are not specified by a standard such as BT.709. Instead, all that matters is that the tri-stimulus values for a given color produce the specified RGB values. The manufacturer is free to use sensors with sensitivities that do not match the standard XYZ definitions, so long as they can later be converted (through a linear transform) to the standard colors.

Similarly, while TV and computer monitors are supposed to produce RGB values as specified by Equation (2.109), there is no reason that they cannot use digital logic to transform the incoming RGB values into different signals to drive each of the color channels. $ ^{26} $ Properly calibrated monitors make this information available to software applications that perform color management, so that colors in real life, on the screen, and on the printer all match as closely as possible.

Color filter arrays

While early color TV cameras used three vidicons (tubes) to perform their sensing and later cameras used three separate RGB sensing chips, most of today's digital still and video cameras use a color filter array (CFA), where alternating sensors are covered by different colored filters (Figure 2.24). $ ^{27} $

The most commonly used pattern in color cameras today is the Bayer pattern (Bayer 1976), which places green filters over half of the sensors (in a checkerboard pattern), and red and blue filters over the remaining ones (Figure 2.31). The reason that there are twice as many green filters as red and blue is because the luminance signal is mostly determined by green values and the visual system is much more sensitive to high-frequency detail in luminance than in chrominance (a fact that is exploited in color image compression—see Section 2.3.3). The process of interpolating the missing color values so that we have valid RGB values for all the pixels is known as demosaicing and is covered in detail in Section 10.3.1.

Similarly, color LCD monitors typically use alternating stripes of red, green, and blue filters placed in front of each liquid crystal active area to simulate the experience of a full color display. As before, because the visual system has higher resolution (acuity) in luminance than chrominance, it is possible to digitally prefilter RGB (and monochrome) images to enhance the perception of crispness (Betrisey, Blinn et al. 2000; Platt 2000b).

2.3 The digital camera

Color balance

Before encoding the sensed RGB values, most cameras perform some kind of color balancing operation in an attempt to move the white point of a given image closer to pure white (equal RGB values). If the color system and the illumination are the same (the BT.709 system uses the daylight illuminant $ D_{65} $ as its reference white), the change may be minimal. However, if the illuminant is strongly colored, such as incandescent indoor lighting (which generally results in a yellow or orange hue), the compensation can be quite significant.

A simple way to perform color correction is to multiply each of the RGB values by a different factor (i.e., to apply a diagonal matrix transform to the RGB color space). More complicated transforms, which are sometimes the result of mapping to XYZ space and back, actually perform a color twist, i.e., they use a general $ 3 \times 3 $ color transform matrix. $ ^{28} $ Exercise 2.8 has you explore some of these issues.

Gamma

In the early days of black and white television, the phosphors in the CRT used to display the TV signal responded non-linearly to their input voltage. The relationship between the voltage and the resulting brightness was characterized by a number called $ \gamma $ (γ), since the formula was roughly

$$ B=V^{\gamma}, $$

with a $ \gamma $ of about 2.2. To compensate for this effect, the electronics in the TV camera would pre-map the sensed luminance Y through an inverse gamma,

$$ Y^{\prime}=Y^{\frac{1}{\gamma}}, $$

with a typical value of $ \frac{1}{\gamma}=0.45 $.

The mapping of the signal through this non-linearity before transmission had a beneficial side effect: noise added during transmission (remember, these were analog days!) would be reduced (after applying the gamma at the receiver) in the darker regions of the signal where it was more

visible (Figure 2.32). $ ^{29} $ (Remember that our visual system is roughly sensitive to relative differences in luminance.)

When color television was invented, it was decided to separately pass the red, green, and blue signals through the same gamma non-linearity before combining them for encoding. Today, even though we no longer have analog noise in our transmission systems, signals are still quantized during compression (see Section 2.3.3), so applying inverse gamma to sensed values remains useful.

Unfortunately, for both computer vision and computer graphics, the presence of gamma in images is often problematic. For example, the proper simulation of radiometric phenomena such as shading (see Section 2.2 and Equation (2.88)) occurs in a linear radiance space. Once all of the computations have been performed, the appropriate gamma should be applied before display. Unfortunately, many computer graphics systems (such as shading models) operate directly on RGB values and display these values directly. (Fortunately, newer color imaging standards such as the 16-bit scRGB use a linear space, which makes this less of a problem (Glassner 1995).)

In computer vision, the situation can be even more daunting. The accurate determination of surface normals, using a technique such as photometric stereo (Section 13.1.1) or even a simpler operation such as accurate image deblurring, require that the measurements be in a linear space of intensities. Therefore, it is imperative when performing detailed quantitative computations such as these to first undo the gamma and the per-image color re-balancing in the sensed color values. Chakrabarti, Scharstein, and Zickler (2009) develop a sophisticated 24-parameter model that is a good match to the processing performed by today's digital cameras; they also provide a database of color images you can use for your own testing.

For other vision applications, however, such as feature detection or the matching of signals in stereo and motion estimation, this linearization step is often not necessary. In fact, determining whether it is necessary to undo gamma can take some careful thinking, e.g., in the case of compensating for exposure variations in image stitching (see Exercise 2.7).

If all of these processing steps sound confusing to model, they are. Exercise 2.9 has you try to tease apart some of these phenomena using empirical investigation, i.e., taking pictures of color charts and comparing the RAW and JPEG compressed color values.

Other color spaces

While RGB and XYZ are the primary color spaces used to describe the spectral content (and hence tri-stimulus response) of color signals, a variety of other representations have been developed both in video and still image coding and in computer graphics.

The earliest color representation developed for video transmission was the YIQ standard developed for NTSC video in North America and the closely related YUV standard developed for PAL in Europe. In both of these cases, it was desired to have a luma channel Y (so called since it only roughly mimics true luminance) that would be comparable to the regular black-and-white TV signal, along with two lower frequency chroma channels.

In both systems, the Y signal (or more appropriately, the Y' luma signal since it is gamma compressed) is obtained from

$$ Y_{601}^{\prime}=0.299R^{\prime}+0.587G^{\prime}+0.114B^{\prime}, $$

where $ R^{\prime}G^{\prime}B^{\prime} $ is the triplet of gamma-compressed color components. When using the newer color definitions for HDTV in BT.709, the formula is

$$ Y^{\prime}_{709}=0.2125R^{\prime}+0.7154G^{\prime}+0.0721B^{\prime}. $$

2.3 The digital camera

The UV components are derived from scaled versions of $ (B^{\prime}-Y^{\prime}) $ and $ (R^{\prime}-Y^{\prime}) $, namely,

$$ U=0.492111(B^{\prime}-Y^{\prime})\quad\mathrm{a n d}\quad V=0.877283(R^{\prime}-Y^{\prime}), $$

whereas the IQ components are the UV components rotated through an angle of $ 33^{\circ} $. In composite (NTSC and PAL) video, the chroma signals were then low-pass filtered horizontally before being modulated and superimposed on top of the $ Y' $ luma signal. Backward compatibility was achieved by having older black-and-white TV sets effectively ignore the high-frequency chroma signal (because of slow electronics) or, at worst, superimposing it as a high-frequency pattern on top of the main signal.

While these conversions were important in the early days of computer vision, when frame grabers would directly digitize the composite TV signal, today all digital video and still image compression standards are based on the newer YCbCr conversion. YCbCr is closely related to YUV (the $ C_b $ and $ C_r $ signals carry the blue and red color difference signals and have more useful mnemonics than UV) but uses different scale factors to fit within the eight-bit range available with digital signals.

For video, the $ Y' $ signal is re-scaled to fit within the [16...235] range of values, while the Cb and Cr signals are scaled to fit within [16...240] (Gomes and Velho 1997; Fairchild 2013). For still images, the JPEG standard uses the full eight-bit range with no reserved values,

$$ \begin{aligned}\begin{bmatrix}Y^{\prime}\\C_{b}\\C_{r}\end{bmatrix}&=\begin{bmatrix}0.299&0.587&0.114\\-0.168736&-0.331264&0.5\\0.5&-0.418688&-0.081312\end{bmatrix}\begin{bmatrix}R^{\prime}\\G^{\prime}\\B^{\prime}\end{bmatrix}+\begin{bmatrix}0\\128\\128\end{bmatrix},\end{aligned} $$

where the R'G'B' values are the eight-bit gamma-compressed color components (i.e., the actual RGB values we obtain when we open up or display a JPEG image). For most applications, this formula is not that important, since your image reading software will directly provide you with the eight-bit gamma-compressed R'G'B' values. However, if you are trying to do careful image deblocking (Exercise 4.3), this information may be useful.

Another color space you may come across is hue, saturation, value (HSV), which is a projection of the RGB color cube onto a non-linear chroma angle, a radial saturation percentage, and a luminance-inspired value. In more detail, value is defined as either the mean or maximum color value, saturation is defined as scaled distance from the diagonal, and hue is defined as the direction around a color wheel (the exact formulas are described by Hall (1989), Hughes, van Dam et al. (2013), and Brown (2019)). Such a decomposition is quite natural in graphics applications such as color picking (it approximates the Munsell chart for color description). Figure 2.331–n shows an HSV representation of a sample color image, where saturation is encoded using a gray scale (saturated = darker) and hue is depicted as a color.

If you want your computer vision algorithm to only affect the value (luminance) of an image and not its saturation or hue, a simpler solution is to use either the $ Y_{xy} $ (luminance + chromaticity) coordinates defined in (2.105) or the even simpler color ratios,

$$ r=\frac{R}{R+G+B},\quad g=\frac{G}{R+G+B},\quad b=\frac{B}{R+G+B} $$

(Figure 2.33e–h). After manipulating the luma (2.113), e.g., through the process of histogram equalization (Section 3.1.4), you can multiply each color ratio by the ratio of the new to old luma to obtain an adjusted RGB triplet.

While all of these color systems may sound confusing, in the end, it often may not matter that much which one you use. Poynton, in his Color FAQ, https://www.poynton.com/ColorFAQ.

html, notes that the perceptually motivated L*a*b* system is qualitatively similar to the gamma-compressed R'/G'/B' system we mostly deal with, since both have a fractional power scaling (which approximates a logarithmic response) between the actual intensity values and the numbers being manipulated. As in all cases, think carefully about what you are trying to accomplish before deciding on a technique to use.

2.3.3 Compression

The last stage in a camera's processing pipeline is usually some form of image compression (unless you are using a lossless compression scheme such as camera RAW or PNG).

All color video and image compression algorithms start by converting the signal into YCbCr (or some closely related variant), so that they can compress the luminance signal with higher fidelity than the chrominance signal. (Recall that the human visual system has poorer frequency response to color than to luminance changes.) In video, it is common to subsample Cb and Cr by a factor

2.3 The digital camera

of two horizontally; with still images (JPEG), the subsampling (averaging) occurs both horizontally and vertically.

Once the luminance and chrominance images have been appropriately subsampled and separated into individual images, they are then passed to a block transform stage. The most common technique used here is the discrete cosine transform (DCT), which is a real-valued variant of the discrete Fourier transform (DFT) (see Section 3.4.1). The DCT is a reasonable approximation to the Karhunen–Loève or eigenvalue decomposition of natural image patches, i.e., the decomposition that simultaneously packs the most energy into the first coefficients and diagonalizes the joint covariance matrix among the pixels (makes transform coefficients statistically independent). Both MPEG and JPEG use $ 8 \times 8 $ DCT transforms (Wallace 1991; Le Gall 1991), although newer variants, including the new AV1 open standard, $ ^{30} $ use smaller $ 4 \times 4 $ or even $ 2 \times 2 $ blocks. Alternative transformations, such as wavelets (Taubman and Marcellin 2002) and lapped transforms (Malvar 1990, 1998, 2000) are used in compression standards such as JPEG 2000 and JPEG XR.

After transform coding, the coefficient values are quantized into a set of small integer values that can be coded using a variable bit length scheme such as a Huffman code or an arithmetic code (Wallace 1991; Marpe, Schwarz, and Wiegand 2003). (The DC (lowest frequency) coefficients are also adaptively predicted from the previous block's DC values. The term "DC" comes from "direct current", i.e., the non-sinusoidal or non-alternating part of a signal.) The step size in the quantization is the main variable controlled by the quality setting on the JPEG file (Figure 2.34).

With video, it is also usual to perform block-based motion compensation, i.e., to encode the difference between each block and a predicted set of pixel values obtained from a shifted block in the previous frame. (The exception is the motion-JPEG scheme used in older DV camcorders, which is nothing more than a series of individually JPEG compressed image frames.) While basic MPEG uses $ 16 \times 16 $ motion compensation blocks with integer motion values (Le Gall 1991), newer standards use adaptively sized blocks, sub-pixel motions, and the ability to reference blocks from older frames (Sullivan, Ohm et al. 2012). In order to recover more gracefully from failures and to allow for random access to the video stream, predicted P frames are interleaved among independently coded I frames. (Bi-directional B frames are also sometimes used.)

The quality of a compression algorithm is usually reported using its peak signal-to-noise ratio (PSNR), which is derived from the average mean square error,

$$ M S E=\frac{1}{n}\sum_{\mathbf{x}}\left[I(\mathbf{x})-\hat{I}(\mathbf{x})\right]^{2}, $$

where $ I(\mathbf{x}) $ is the original uncompressed image and $ \hat{I}(\mathbf{x}) $ is its compressed counterpart, or equivalently, the root mean square error (RMS error), which is defined as

$$ RMS=\sqrt{MSE}. $$

The PSNR is defined as

$$ PSNR=10\log_{10}\frac{I^{2}_{\max}}{MSE}=20\log_{10}\frac{I_{\max}}{RMS}, $$

where $ I_{max} $ is the maximum signal extent, e.g., 255 for eight-bit images.

While this is just a high-level sketch of how image compression works, it is useful to understand so that the artifacts introduced by such techniques can be compensated for in various computer vision applications. Note also that researchers are currently developing novel image and video compression algorithms based on deep neural networks, e.g., (Rippel and Bourdev 2017; Mentzer, Agustsson et al. 2019; Rippel, Nair et al. 2019) and https://www.compression.cc. It will be interesting to see what kinds of different artifacts these techniques produce.

2.4 Additional reading

As we mentioned at the beginning of this chapter, this book provides but a brief summary of a very rich and deep set of topics, traditionally covered in a number of separate fields.

A more thorough introduction to the geometry of points, lines, planes, and projections can be found in textbooks on multi-view geometry (Faugeras and Luong 2001; Hartley and Zisserman 2004) and computer graphics (Watt 1995; OpenGL-ARB 1997; Hughes, van Dam et al. 2013; Marschner and Shirley 2015). Topics covered in more depth include higher-order primitives such as quadrics, conics, and cubics, as well as three-view and multi-view geometry.

The image formation (synthesis) process is traditionally taught as part of a computer graphics curriculum (Glassner 1995; Watt 1995; Hughes, van Dam et al. 2013; Marschner and Shirley 2015) but it is also studied in physics-based computer vision (Wolff, Shafer, and Healey 1992a). The behavior of camera lens systems is studied in optics (Möller 1988; Ray 2002; Hecht 2015).

Some good books on color theory have been written by Healey and Shafer (1992), Wandell (1995), Wyszecki and Stiles (2000), and Fairchild (2013), with Livingstone (2008) providing a more fun and informal introduction to the topic of color perception. Mark Fairchild's page of color books and links $ ^{31} $ lists many other sources.

Topics relating to sampling and aliasing are covered in textbooks on signal and image processing (Crane 1997; Jähne 1997; Oppenheim and Schafer 1996; Oppenheim, Schafer, and Buck 1999; Pratt 2007; Russ 2007; Burger and Burge 2008; Gonzalez and Woods 2017).

Two courses that cover many of the above topics (image formation, lenses, color and sampling theory) in wonderful detail are Marc Levoy's Digital Photography course at Stanford (Levoy 2010) and Michael Brown's tutorial on the image processing pipeline at ICCV 2019 (Brown 2019). The recent book by Ikeuchi, Matsushita et al. (2020) also covers 3D geometry, photometry, and sensor models, but with an emphasis on active illumination systems.

2.5 Exercises

A note to students: This chapter is relatively light on exercises since it contains mostly background material and not that many usable techniques. If you really want to understand multi-view geometry

2.5 Exercises

in a thorough way, I encourage you to read and do the exercises provided by Hartley and Zisserman (2004). Similarly, if you want some exercises related to the image formation process, Glassner's (1995) book is full of challenging problems.

Ex 2.1: Least squares intersection point and line fitting—advanced. Equation (2.4) shows how the intersection of two 2D lines can be expressed as their cross product, assuming the lines are expressed as homogeneous coordinates.

1. If you are given more than two lines and want to find a point $ \tilde{x} $ that minimizes the sum of squared distances to each line,

$$ D=\sum_{i}(\tilde{\mathbf{x}}\cdot\mathbf{\tilde{l}}_{i})^{2}, $$

how can you compute this quantity? (Hint: Write the dot product as $ \tilde{\mathbf{x}}^{T}\tilde{\mathbf{l}}_{i} $ and turn the squared quantity into a quadratic form, $ \tilde{\mathbf{x}}^{T}\mathbf{A}\tilde{\mathbf{x}} $.)

2. To fit a line to a bunch of points, you can compute the centroid (mean) of the points as well as the covariance matrix of the points around this mean. Show that the line passing through the centroid along the major axis of the covariance ellipsoid (largest eigenvector) minimizes the sum of squared distances to the points.

3. These two approaches are fundamentally different, even though projective duality tells us that points and lines are interchangeable. Why are these two algorithms so apparently different? Are they actually minimizing different objectives?

Ex 2.2: 2D transform editor. Write a program that lets you interactively create a set of rectangles and then modify their "pose" (2D transform). You should implement the following steps:

1. Open an empty window ("canvas").

2. Shift drag (rubber-band) to create a new rectangle.

3. Select the deformation mode (motion model): translation, rigid, similarity, affine, or perspective.

4. Drag any corner of the outline to change its transformation.

This exercise should be built on a set of pixel coordinate and transformation classes, either implemented by yourself or from a software library. Persistence of the created representation (save and load) should also be supported (for each rectangle, save its transformation).

Ex 2.3: 3D viewer. Write a simple viewer for 3D points, lines, and polygons. Import a set of point and line commands (primitive) as well as a viewing transform. Interactively modify the object or camera transform. This viewer can be an extension of the one you created in Exercise 2.2. Simply replace the viewing transformations with their 3D equivalents.

(Optional) Add a z-buffer to do hidden surface removal for polygons.

(Optional) Use a 3D drawing package and just write the viewer control.

Ex 2.4: Focus distance and depth of field. Figure out how the focus distance and depth of field indicators on a lens are determined.

1. Compute and plot the focus distance $ z_o $ as a function of the distance traveled from the focal length $ \Delta z_i = f - z_i $ for a lens of focal length $ f $ (say, 100mm). Does this explain the hyperbolic progression of focus distances you see on a typical lens (Figure 2.20)?

2. Compute the depth of field (minimum and maximum focus distances) for a given focus setting $ z_{o} $ as a function of the circle of confusion diameter c (make it a fraction of the sensor width), the focal length f, and the f-stop number N (which relates to the aperture diameter d). Does this explain the usual depth of field markings on a lens that bracket the in-focus marker, as in Figure 2.20a?

3. Now consider a zoom lens with a varying focal length $f$. Assume that as you zoom, the lens stays in focus, i.e., the distance from the rear nodal point to the sensor plane $z_{i}$ adjusts itself automatically for a fixed focus distance $z_{o}$. How do the depth of field indicators vary as a function of focal length? Can you reproduce a two-dimensional plot that mimics the curved depth of field lines seen on the lens in Figure 2.20b?

Ex 2.5: F-numbers and shutter speeds. List the common f-numbers and shutter speeds that your camera provides. On older model SLRs, they are visible on the lens and shutter speed dials. On newer cameras, you have to look at the electronic viewfinder (or LCD screen/indicator) as you manually adjust exposures.

1. Do these form geometric progressions; if so, what are the ratios? How do these relate to exposure values (EVs)?

2. If your camera has shutter speeds of $ \frac{1}{60} $ and $ \frac{1}{125} $, do you think that these two speeds are exactly a factor of two apart or a factor of 125/60 = 2.083 apart?

3. How accurate do you think these numbers are? Can you devise some way to measure exactly how the aperture affects how much light reaches the sensor and what the exact exposure times actually are?

Ex 2.6: Noise level calibration. Estimate the amount of noise in your camera by taking repeated shots of a scene with the camera mounted on a tripod. (Purchasing a remote shutter release is a good investment if you own a DSLR.) Alternatively, take a scene with constant color regions (such as a color checker chart) and estimate the variance by fitting a smooth function to each color region and then taking differences from the predicted function.

1. Plot your estimated variance as a function of level for each of your color channels separately.

2. Change the ISO setting on your camera; if you cannot do that, reduce the overall light in your scene (turn off lights, draw the curtains, wait until dusk). Does the amount of noise vary a lot with ISO/gain?

3. Compare your camera to another one at a different price point or year of make. Is there evidence to suggest that “you get what you pay for”? Does the quality of digital cameras seem to be improving over time?

Ex 2.7: Gamma correction in image stitching. Here's a relatively simple puzzle. Assume you are given two images that are part of a panorama that you want to stitch (see Section 8.2). The two images were taken with different exposures, so you want to adjust the RGB values so that they match along the seam line. Is it necessary to undo the gamma in the color values in order to achieve this?

Ex 2.8: White point balancing—tricky. A common (in-camera or post-processing) technique for performing white point adjustment is to take a picture of a white piece of paper and to adjust the RGB values of an image to make this a neutral color.

2.5 Exercises

1. Describe how you would adjust the RGB values in an image given a sample “white color” of $ (R_{w}, G_{w}, B_{w}) $ to make this color neutral (without changing the exposure too much).

2. Does your transformation involve a simple (per-channel) scaling of the RGB values or do you need a full $ 3 \times 3 $ color twist matrix (or something else)?

3. Convert your RGB values to XYZ. Does the appropriate correction now only depend on the XY (or xy) values? If so, when you convert back to RGB space, do you need a full $ 3 \times 3 $ color twist matrix to achieve the same effect?

4. If you used pure diagonal scaling in the direct RGB mode but end up with a twist if you work in XYZ space, how do you explain this apparent dichotomy? Which approach is correct? (Or is it possible that neither approach is actually correct?)

If you want to find out what your camera actually does, continue on to the next exercise.

Ex 2.9: In-camera color processing—challenging. If your camera supports a RAW pixel mode, take a pair of RAW and JPEG images, and see if you can infer what the camera is doing when it converts the RAW pixel values to the final color-corrected and gamma-compressed eight-bit JPEG pixel values.

1. Deduce the pattern in your color filter array from the correspondence between co-located RAW and color-mapped pixel values. Use a color checker chart at this stage if it makes your life easier. You may find it helpful to split the RAW image into four separate images (subsampling even and odd columns and rows) and to treat each of these new images as a “virtual” sensor.

2. Evaluate the quality of the demosaicing algorithm by taking pictures of challenging scenes which contain strong color edges (such as those shown in Section 10.3.1).

3. If you can take the same exact picture after changing the color balance values in your camera, compare how these settings affect this processing.

4. Compare your results against those presented in (Chakrabarti, Scharstein, and Zickler 2009), Kim, Lin et al. (2012), Hasinoff, Sharlet et al. (2016), Karaimer and Brown (2016), and Brooks, Mildenhall et al. (2019) or use the data available in their database of color images.