Patent Publication Number: US-11644682-B2

Title: Systems and methods for diffraction line imaging

Description:
RELATED APPLICATIONS 
     This application claims the benefit of U.S. Provisional Patent Application Ser. No. 63/037,888, filed Jun. 11, 2020, the contents of which are incorporated herein in their entirety. 
    
    
     GOVERNMENT INTEREST 
     This invention was made with government support under Awards Nos. 1730147 and 1900821 issued by the National Science Foundation and Award No. HR0011-16-C-0025 issue by DARPA/Wisconsin (A022264). The government has certain rights in this invention. 
    
    
     BACKGROUND OF THE INVENTION 
     Light passing through a diffraction grating is dispersed as a function of wavelength (measured color). When imaged by a 2D camera, the dispersed light manifests as colorful streaks. Such diffraction gratings are commonly used to create artistic effects like the star-filter commonly used in photography when capturing bright light sources such as vehicle head-lamps and streetlights at night, as shown in  FIG.  1 B . The micro-structure of the grating influences the shapes and number of streaks. The brightness and appearance of the streaks depends on the 2D spatial locations of the sources in the image. Diffraction gratings are used in many imaging and scientific applications, including spectroscopy, multi-spectral sensing and rainbow particle velocimetry. 
     SUMMARY OF THE INVENTION 
     When light sources or bright projections of light, as shown in  FIG.  1 A , are viewed through a transmissive diffraction grating, the incoming light is dispersed as a function of wavelength (measured color). When imaged by a 2D camera, the dispersed light manifests as colorful streaks, creating streaks of light on the camera&#39;s sensor plane, as shown in  FIG.  1 B . If a one-dimensional color sensor, shown as the dotted green line in  FIG.  1 B , or a very few rows in a 2D sensor, intersects with these streaks, the measured colors can be used to efficiently determine the 2D spatial positions for each light source or projected line pixel at high frame rates. For example, the central column in the right-most diffraction image in  FIG.  1 B  can be sensed at 2.6K fps by a camera. 
     Disclosed herein is a novel class of imaging systems that use one or more diffraction gratings in conjunction with line (1D) sensors. Line sensors can yield 2D and 3D information. Line sensors offer several advantages over 2D sensors, namely, high frame rates, high dynamic range, and additional on-chip computation near the sensing area. As a result, this imaging approach results in fast 2D positioning and 3D reconstruction in many applications, including, for example, night imaging, LED-based motion capture, and industrial product scanning. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee. 
         FIGS.  1 A-B  contrast viewing light sources or bright projections of light through standard optics versus through a transmissive diffraction grating. 
         FIGS.  2 A-C  are a graphic representation of various diffraction line imaging applications. 
         FIGS.  3 A-C  are a graphic representation of the diffraction-based 2D positioning method of the present invention. 
         FIGS.  4 A-B  are a schematic representation of the components of the present invention. 
         FIG.  5    shows a method for recovering sparse point sources using a convolutional neural network. 
         FIGS.  6 A-C  show a method utilizing a 2D camera having multiple ROIs; a first alternate embodiment utilizing a horizontal cylindrical lens; and a second alternate embodiment using multi-axis diffraction, respectively. 
         FIG.  7    shows results from the system when applied to high speed operation. 
         FIG.  8    shows results from a double-axis diffraction system having three ROIs. 
         FIG.  9    shows a method for dense scene deconvolution. 
         FIG.  10    shows the process for performing the dense scene deconvolution shown in  FIG.  9   . 
         FIG.  11    is a graphical representation of the convolutional image formation and recovery. 
         FIG.  12    graphically shows the recovery of sources of multiple spectra. 
         FIG.  13    graphically shows the effect of PSF and sensor sampling on position recovery. 
         FIG.  14    graphically shows the process of estimating the uniquely recoverable domain. 
     
    
    
     DEFINITIONS 
     As used herein, the term “light source” refers to an illuminated object in a scene, wherein the object is illuminated either by emitting its own light or by reflecting light from another source. 
     As used herein, the term “X-Y position” or “X-Y coordinates” refers to any orientation of coordinates in an 2d plane and is not meant to be limited to coordinate system having axes oriented horizontally or vertically. 
     DETAILED DESCRIPTION 
     As a starting point in the explanation of the invention, a detailed image formation model for a collimated beam of light incident on a line sensor via a diffraction grating will be presented. This model maps a 3D point source to a 2D X-Y location on a virtual image plane and depends on the relative placement of the line sensor with respect to the diffraction grating and the line sensor&#39;s spectral response. 
     A learning-based algorithm is presented for estimating the 2D locations of a sparse set of light sources or a line illumination projected onto a 3D scene. An analysis is provided of the uncertainty of the 2D positioning and the achieved field-of-view. To improve positioning accuracy significantly, the imaging system is extended to include multiple diffraction gratings and/or an additional cylindrical lens. Finally, the approach is extended to multiple line sensors (or using multiple regions of interest in a 2D sensor) to increase the field-of-view of the imaging system. 
     The approach of the present invention can be thought of as a variant of compressive line sensing with direct decoding of 2D positions, in place of computationally intensive and noise-prone decoding algorithms. 
     The present invention is explained in terms of two applications: (1) the passive imaging of light sources; and (2) Structured Light Line Scanning. 
     Passive Imaging of Light Sources. Night scenes include a variety of light sources such as streetlamps, vehicle headlamps, taillights, turn signals, and bicycle lights, as shown in  FIG.  1 A . Because of their distinctive appearances, light sources can serve as strong features to track, reconstruct and interpret the scene and objects of interest at night. As an example, the flicker of sources may be used to determine the AC phase of electrical circuits and even analyze power grids. In another example, the glows around streetlamps may be used to reveal the weather condition (fog, haze, rain) and visibility. Finally, motion capture systems often attach light sources (LEDs) to estimate the object&#39;s motion in 3D. Capturing sparse light sources with 2D sensors wastes pixel resources and limits the image acquisition rate. The present invention is able to estimate 2D light source positions from line images at high frame rates (up to 2200 fps in the system described herein). In contrast to previous works that use line sensors or spatio-temporal illumination to light a subject, one embodiment of the present invention is based purely on passive observation of light sources without requiring special modulation or synchronization.  FIG.  2 A  shows an example of the passive imaging of light sources, specifically, unmodulated LEDs and car headlights. 
     Structured Light Line Scanning: Structured light 3D scanning often projects and sweeps an illumination line across a scene. A 2D camera observes the intersection of the line with the 3D scene. The 3D scene is then reconstructed by intersecting the pixel ray with the plane of illumination. In this scenario, most of the pixels in the 2D camera are not used, thereby limiting the acquisition rate due to data transfer bandwidth. The scene can be accurately reconstructed using a line light source and diffraction line imaging. The image acquisition of the present invention is limited mainly by the SNR (and hence exposure time) and not by the bandwidth of data transfer from a 2D sensor. This embodiment of the present invention represents a first instance of a structured light system with a 1D source and a 1D camera. This was not possible before and either a 2D projector and/or a 2D camera was required. Further, line illumination and imaging is significantly more energy efficient than 2D illumination and imaging. If the line source is bright enough, it is possible to achieve scan rates of tens of thousands of lines per second, making this approach very useful in industrial/logistics scanning applications.  FIG.  2 B  shows an application of structured light scanning in which objects on a conveyer belt are illuminated by line illumination and scanned. 
     Diffraction Gratings 
     The approach of the present invention to light source positioning starts by exploiting diffraction, which is a wavelength-dependent optical phenomenon that causes light to bend (i.e., diffract) when passing through and around obstacles and slits. A diffraction grating is an optical element having periodic structures forming slits that produces diffraction patterns, useful in scientific applications such as spectroscopy, holography and hyper spectral imaging. 
       FIG.  3 A  shows a grating consisting of a periodic structure repeated every d microns. A monochromatic collimated beam passes through the diffraction grating. The grating diffracts the beam yielding new beam directions according to Eq. (1). In a transmissive grating, incident light diffracts according to:
 
 d [sin(θ m )−sin(θ i )]= mλ   (1)
 
     where: 
     θ m  is the angle of the m th  diffraction order; 
     θ i  is the angle of incident light; 
     d is the distance between the periodic structures in the diffraction grating; 
     m is the diffraction order; and 
     λ is the wavelength, in microns. 
     In Eq. (1), m=0 corresponds to the 0 th -order component, which passes unbent through the grating, while each wavelength disperses into a unique angle for m≠0, producing the characteristic point spread function (PSF) which appears as a rainbow color streak, as shown in  FIG.  1 A . Per diffraction order m, each wavelength is diffracted in a different spatial direction. When imaged by a camera, this results in a vertical rainbow pattern on the image plane. Here, pixel p 0  measures the energy of some wavelength in the green range. 
     For a fixed θ i , a camera imaging the light exiting the diffraction grating maps each wavelength onto unique sensor position. This mapping is the basis for spectroscopy, the purpose of which is the spectral analysis of the incident light. For example, echelle spectroscopes use combinations of high-dispersion and low-dispersion gratings to precisely measure the intensity of thousands of different wavelengths. 
     In contrast, the present invention inverts this process, using color to efficiently and precisely recover the incident direction of light θ i  for unknown light sources.  FIG.  3 C  shows that shifting the incoming beam&#39;s direction results in a spatial shift of the spectral pattern on the image plane. Each incoming beam direction maps a unique wavelength to p 0 . Thus, the color value at p 0  provides information about the incident angle θ i . Consider light exiting the diffraction grating at an angle θ′ belonging to the first diffraction order (i.e., m=1). Suppose that, on the camera image plane, light exiting at θ′ is focused at some corresponding camera pixel coordinate p 0 , as illustrated in  FIG.  3 B . Then, from Eq. (1), the color (wavelength) measured at p 0  is:
 
λ(θ i )= d [sin(θ′)−sin(θ i )]  (2)
 
     and depends on the incident light direction θ i . Thus, for a fixed θ′, the measured color at p 0  indicates the direction of incident light. 
     Diffraction-Based 2D Positioning 
     The diffraction equation (Eq. (2)) has two key properties that enable computing the direction of scene light sources. First, the function λ(θ i ) is an injective (i.e., one-to-one) function across the domain θ i ∈[−p i /2,θ′]. The inverse of λ(θ i ) is therefore given by: 
     
       
         
           
             
               
                 
                   
                     
                       θ 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       λ 
                       ) 
                     
                   
                   = 
                   
                     
                       sin 
                       
                         - 
                         1 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           sin 
                           ⁡ 
                           
                             ( 
                             
                               θ 
                               ′ 
                             
                             ) 
                           
                         
                         - 
                         
                           λ 
                           d 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Second, Eq. (2) and its inverse do not depend on the intensity of light. This property makes the imaging system of the present invention effectively invariant to the incident light intensity and, more specifically, its spectral characteristics. Note that this assumes no under- and over-saturated pixel measurements. 
     A basic schematic of the optical system of an embodiment of the present invention is shown in  FIG.  4 A  (top view) and  FIG.  4 B  (front view). The system consists of a diffraction grating and a color line scan camera. The line scan camera is positioned vertically, containing pixels along the Y-direction (in the camera&#39;s reference frame) denoted by the single coordinate y. Note that the terms “line scan camera” and “1D camera” are used interchangeably herein. For ease of explanation, the line scan camera&#39;s image plane is sometimes referred to as the full 2D image, had the camera been equipped with a standard 2D sensor. 
     With reference to  FIG.  3 A , the first order diffraction is measured by an RGB line scan camera. The measured colors are used to estimate the 2D light source positions. Light incident on the diffraction grating is dispersed as a rainbow streak onto the sensor.  FIG.  3 (B)  shows the situation when the scene contains sparse light-sources, as would be the case, for example, for motion capture. In this case, the line scan camera contains several spikes. Spike location along the vertical line scan camera encodes the vertical coordinate in the virtual image plane, while the spike color encodes the horizontal coordinate.  FIG.  3 (C)  shows 3D scanning using line projection. The projected line yields a piece-wise continuous signal in the line scan camera. 
     The system can then track both 2D spatial coordinates. For each point source, the vertical source coordinate is trivially measured by the rainbow&#39;s streak position along the 1D vertical sensor (i.e., y 1  in  FIG.  4 A ). Computing the horizontal coordinate of the source amounts to computing the θ i  angle, which involves three simple steps: (1) measure the response in RGB space; (2) compute the dominant wavelength A that produces the RGB signal; and (3) evaluate the inverse function θ i (λ). Note that this procedure requires the light sources to be sufficiently broadband, such that the RGB measurement is non-zero for the given incident angle θ i  and the corresponding wavelength λ(θ i ). 
     As previously mentioned, two types of imaging regimes are addressed: (1) sparse light sources. That is, scenes having a sparse set of point sources distributed across the image plane, as shown in  FIG.  4 B ; and (2) structured light line scanning. That is, scenes illuminated by a vertically projected line, which yields a vertical curve on the image plane, as shown in  FIG.  4 C . 
     3D scanning using a projected line can be implemented in various ways. The projected line (i.e., a plane in space) can either be swept across a static object by a projector or can have a fixed direction while the object moves through the plane, for example, by placing the object on a conveyor belt or turn table. 
     Image Formation Model 
     The model connecting the projection of light from 3D positions in space onto a 2D virtual camera image plane is derived here. The line scan camera, modeled as a pinhole, is set to image the first-order diffracted light from the scene. Let O cam  denote the position and orientation of the line scan camera with respect to the diffraction grating. As shown in  FIG.  4 A , a virtual image plane U is defined with pixel coordinates p v =(x v , y v ), that belong to a virtual 2D camera positioned at O vir . 
     The model for a single point source indexed by n=1, having homogeneous world coordinates w 4 ≡[X, Y, Z, 1] T  is first described. On the image plane of the line scan camera, the PSF of the point source creates a holographic image in the shape of a horizontal rainbow line, as shown in  FIG.  1 B  and  FIG.  5   . The line sensor of the camera is orthogonal to the rainbow line, intersecting it over a few pixels centered at y 1 . The geometric relationship between y 1  and the corresponding virtual image coordinate y 1   v  is:
 
 y   1   v   ≡c ( y   1 )  (4)
 
     Neglecting distortions that arise from light entering at highly oblique angles, the y 1  coordinate of the rainbow line is given by standard projective geometry: 
     
       
         
           
             
               
                 
                   
                     γ 
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             
                               y 
                               1 
                             
                           
                         
                         
                           
                             1 
                           
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ( 
                       
                         
                           
                             0 
                           
                           
                             1 
                           
                           
                             0 
                           
                         
                         
                           
                             0 
                           
                           
                             0 
                           
                           
                             1 
                           
                         
                       
                       ) 
                     
                     ⁢ 
                     
                       Pw 
                       1 
                     
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     where: 
     P is the projection matrix of the camera; and 
     γ is an arbitrary scale factor. 
     In turn, since the Y-axes of both the virtual and line sensors are identical, G can be approximated by an affine transformation. 
     As seen in of  FIG.  4 A , the angle of the incident light with respect to the grating can be expressed as: 
     
       
         
           
             
               
                 
                   
                     
                       θ 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       
                         x 
                         i 
                         v 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       tan 
                       
                         - 
                         1 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           [ 
                           
                             
                               x 
                               t 
                               v 
                             
                             - 
                             
                               W 
                               2 
                             
                           
                           ] 
                         
                         f 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     where: 
     ƒ is the focal length of the virtual camera in pixel units; and 
     W is the virtual image plane width. 
     Combining this with Eq. (3) yields: 
     
       
         
           
             
               
                 
                   
                     λ 
                     ⁡ 
                     
                       ( 
                       
                         x 
                         i 
                         v 
                       
                       ) 
                     
                   
                   = 
                   
                     d 
                     ⁡ 
                     
                       ( 
                       
                         
                           sin 
                           ⁡ 
                           
                             ( 
                             
                               θ 
                               ′ 
                             
                             ) 
                           
                         
                         - 
                         
                           sin 
                           ⁡ 
                           
                             [ 
                             
                               
                                 tan 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     [ 
                                     
                                       
                                         x 
                                         t 
                                         v 
                                       
                                       - 
                                       
                                         W 
                                         2 
                                       
                                     
                                     ] 
                                   
                                   f 
                                 
                                 ) 
                               
                             
                             ] 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     The RGB intensity measured by the camera is modeled as:
 
 I   σ ( y   1 )= Tα   σ [λ( x   i   v )] s [λ( x   i   v )]  (8)
 
     where: 
     T is the exposure time; 
     α σ (λ) is the camera&#39;s spectral response function (grayscale/Joule) for every color channel; 
     σ∈{R, G, B}; and 
     s(λ) represents the spectral radiant flux of the source (Joule/sec) falling on y 1 . 
     In Eq. (8) it is assumed that each pixel y integrates a very narrow band of wavelengths due to the very narrow field-of-view. 
     Normalizing the 3×1 vector I σ (y 1 ) by its L 2  norm removes the dependence on exposure time and source spectral flux: 
     
       
         
           
             
               
                 
                   
                     
                       
                         I 
                         _ 
                       
                       σ 
                     
                     ⁡ 
                     
                       ( 
                       
                         y 
                         1 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         
                           I 
                           σ 
                         
                         ⁡ 
                         
                           ( 
                           
                             y 
                             1 
                           
                           ) 
                         
                       
                       
                         
                            
                           
                             
                               I 
                               σ 
                             
                             ⁡ 
                             
                               ( 
                               
                                 y 
                                 1 
                               
                               ) 
                             
                           
                            
                         
                         2 
                       
                     
                     = 
                     
                       
                         
                           
                             α 
                             σ 
                           
                           ⁡ 
                           
                             [ 
                             
                               λ 
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     θ 
                                     ′ 
                                   
                                   ⁢ 
                                   
                                     x 
                                     t 
                                     v 
                                   
                                 
                                 ) 
                               
                             
                             ] 
                           
                         
                         
                           
                              
                             
                               
                                 α 
                                 σ 
                               
                               ⁡ 
                               
                                 [ 
                                 
                                   λ 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         θ 
                                         ′ 
                                       
                                       ⁢ 
                                       
                                         x 
                                         t 
                                         v 
                                       
                                     
                                     ) 
                                   
                                 
                                 ] 
                               
                             
                              
                           
                           2 
                         
                       
                       ≡ 
                       
                         
                           H 
                           σ 
                         
                         ⁡ 
                         
                           ( 
                           
                             x 
                             i 
                             v 
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     Finally, x i   v  is given by inverting Eq. (9):
 
 x   i   v   =H   −1 [ Ī   σ ( y   1 )]  (10)
 
     The model in Eqs. (6-10) assumes that only a single source is predominately projected onto every line camera pixel. Recovering x i   v  is described next. 
     Recovering Horizontal Coordinates 
     Computing H −1  in practice requires accounting for various factors—point sources not truly at infinity, 1D camera-diffraction geometry, the line scan camera&#39;s radiometric response, sensor saturation, and more. Due to camera lens distortions, sensor-diffraction grating misalignment, and deviations to Eq. (1) due to incident light in oblique angles, H −1  generally depends on y 1 . Note that Eq. (1) strictly holds when the grating groves are perpendicular to the incidence plane. The incidence plane is the plane containing the beam direction and the normal to diffraction grating plane. These factors make computing H −1  directly a laborious task. Instead, a learning-based approach may be used. 
     A neural network can be trained to approximate H −1 , denoted by Ĥ −1 . The network receives sensor color measurements in RGB space along with y, termed RGBy, and outputs the x v  coordinates. In practice, different network architectures may be used for sparse points and vertical line projection. 
     Recovering Sparse Point Sources 
     For sparse sources, a convolutional neural network (CNN) can be trained. The network, denoted by maps H point   −1  from a small patch into a predefined discrete set of possible x v  coordinates. 
     The process is shown in  FIG.  5   . In  FIG.  5 ( a ) , a plurality of rainbow streaks are incident on the line scan camera image plane.  FIG.  5 ( b )  shows an average profile of a subset of columns shown in  FIG.  5 ( a ) , representing the signal received by the line scan camera.  FIG.  5 ( c )  shows small 8×9 patches capturing the peak response within a narrow region of interest (ROI) of the sensor. The position of this patch is given by coordinate y n .  FIG.  5 ( d )  shows each detected patch and the coordinate y n  is processed by the trained CNN to obtain a horizontal coordinate x n   v . The 2D color glove image from a helper camera is shown for reference only. 
     In a preferred embodiment of the invention, the neural network consists of three convolutional layers (channel widths 20, 50, 50 respectively) followed by two fully connected layers, as shown in  FIG.  5 ( c ) . Note that the described CNN is exemplary only and that other architectures for the CNN are possible without deviating from the scope of the invention. 
     Consider a scene with N point sources, indexed by n=1, 2, . . . , N. Let Q denote the number of vertical pixels in the line scan camera. A line scan camera can be a few pixels wide in the horizontal direction as well (e.g., due to RGB Bayer color filter). Supposed the line scan camera has a width of 8 pixels. Thus, given an 8×Q color image from the line scan camera, the first step is to identify the rainbow line positions from peaks in the sensor measurement I(y). Then, for every detected rainbow with coordinate y n , provide the CNN with a normalized small 8×9 image patch Ω(y n ) vertically centered at y n , and concatenate the coordinate y n  with the input to the first fully connected layer. Note that normalizing an 8×9 patch consists of dividing all RGB pixel values by a scalar. For example, for L ∞ , Ī[Ω(y m )]=I[Ωy m )]/max RGB  where max RGB  is the maximum value across all patch pixels and color channels. 
     The network outputs a W×1 vector s n  with scores per possible horizontal virtual image location. Then, are recovered as:
 
 x   n   v =arg max( S [ s   n ])  (11)
 
     where: 
     S is a softmax function. 
     The training dataset consists of ground truth 2D coordinates captured with a 2D helper camera, line scan camera coordinates y m , and the corresponding image patches Ī[Ω(y m )]: 
     
       
         
           
             
               
                 
                   
                     
                       { 
                       
                         x 
                         m 
                         
                           v 
                           , 
                           GT 
                         
                       
                       } 
                     
                     
                       m 
                       = 
                       1 
                     
                     M 
                   
                   ↔ 
                   
                     
                       { 
                       
                         
                           
                             I 
                             ¯ 
                           
                           ⁡ 
                           
                             [ 
                             
                               Ω 
                               ⁡ 
                               
                                 ( 
                                 
                                   y 
                                   m 
                                 
                                 ) 
                               
                             
                             ] 
                           
                         
                         , 
                         
                           y 
                           m 
                         
                       
                       } 
                     
                     
                       m 
                       = 
                       1 
                     
                     M 
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     where: 
     M is the number of training examples. 
     The training loss is given by: 
     
       
         
           
             
               
                 
                   L 
                   = 
                   
                     
                       1 
                       M 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           m 
                           = 
                           1 
                         
                         M 
                       
                       ⁢ 
                       
                         BCE 
                         ⁡ 
                         
                           ( 
                           
                             
                               S 
                               ⁡ 
                               
                                 [ 
                                 
                                   s 
                                   m 
                                 
                                 ] 
                               
                             
                             , 
                             
                               D 
                               ⁡ 
                               
                                 [ 
                                 
                                   
                                     x 
                                     m 
                                     
                                       v 
                                       , 
                                       GT 
                                     
                                   
                                   , 
                                   σ 
                                 
                                 ] 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     where: 
     BCE(⋅,⋅) is the binary cross entropy function; and 
     D(⋅) generates a Gaussian probability distribution with mean x m   v,GT  and standard deviation σ. 
     Intuitively, the binary cross entropy function drives the output distribution S[s m ] to match the desired distribution, which is a narrow Gaussian centered at x m   v,GT . Using a Gaussian instead of a Delta function for D(⋅) provides a tolerance for small deviation due to image noise.  FIG.  4    shows an example result. 
     Recovering a Vertically Projected Line 
     3D reconstruction relies on an aggregate of many continuous measurements through interpolation. Thus, here the network implements regression. Let the 3×1 RGB vector u y  denote the mean color of patch Ω(y). u y  is computed by averaging Ω(y) over the spatial domain. 
     The network H line   −1  is fed with a normalized u y  along with y and outputs a continuous scalar approximation for x v : 
     
       
         
           
             
               
                 
                   
                     x 
                     v 
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               
                                 H 
                                 line 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     
                                       u 
                                       y 
                                     
                                     
                                       
                                          
                                         
                                           u 
                                           y 
                                         
                                          
                                       
                                       2 
                                     
                                   
                                   , 
                                   y 
                                 
                                 ) 
                               
                             
                             , 
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 u 
                                 y 
                               
                             
                             ∈ 
                             A 
                           
                         
                       
                       
                         
                           
                             none 
                             , 
                           
                         
                         
                           otherwise 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     where 
        is a subspace of “valid” RGB input values such that:
 
 ={ u   y   :∥u   y ∥ inf   &gt;t,∥u   y ∥ 0 &gt;1}
 
     where: 
     t is a predefined intensity threshold. 
     Subspace   excludes low intensity (and thus noisy) measurements as well as wavelengths that map to only a single-color channel and thus can&#39;t be distinguished after normalization. The latter occurs at the edges of the rainbow streak (i.e., deep reds and blues). 
     In preferred embodiments, the network consists of four fully connected layers of size 300 and is trained using the Huber loss. The training dataset is: 
     
       
         
           
             
               
                 
                   
                     
                       { 
                       
                         x 
                         m 
                         
                           v 
                           , 
                           GT 
                         
                       
                       } 
                     
                     
                       m 
                       = 
                       1 
                     
                     M 
                   
                   ↔ 
                   
                     
                       { 
                       
                         
                           u 
                           m 
                         
                         , 
                         
                           y 
                           m 
                         
                       
                       } 
                     
                     
                       m 
                       = 
                       1 
                     
                     M 
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     and was obtained by scanning a vertical line on a white wall. 
     Expanding the FOV using Multiple ROIs 
     Using multiple 1D regions-of-interests (ROIs) in a 2D camera has the advantage of increasing the field-of-view (FOV) without significant sacrifice of acquisition speed. In one embodiment, the camera readout may be set to a few columns (rows in a 90° rotated camera) at the center of the image plane of the camera. However, reading multiple columns positioned at different spatial locations, as shown in  FIG.  6 A , increases the horizontal FOV by catching rainbow streaks whose visible spectrum do not intersect with the center column (i.e., the domain of the central ROI). Multiple ROIs also reduces recovery uncertainty since they may increase the signal per point. This occurs when a single rainbow streak has sufficient signal over several ROIs. 
     In  FIG.  6 B , the additional left camera images the zeroth-order diffracted light through a cylindrical lens which focuses light in the horizontal axis (producing white streaks). The intersection of the white streaks with the left line sensor yields additional horizontal point coordinate measurements and improves source positioning. In  FIG.  6 C , a multi-axis diffraction grating additionally yields horizontal rainbow streaks. Capturing these streaks using an additional line sensor mitigates the vertically-overlapping ambiguity of the sources. 
     Let I σ   r  denote the line scan camera image from ROIs indexed r=1, 2, . . . , R. When R&gt;1, the network is fed with a concatenation of measurements from all ROIs. Specifically, for sparse points, the R patches I σ   r [Ω(y n )] from y n  are concatenated to form an extended 8×9 (R) color patch and feed it to the network, along with y n . For vertical line projection, the RGB measurements u y   r  are similarly concatenated. The resulting (3R+1)×1 vector is fed to the network. As in Eq. (14), only “valid” y&#39;s are considered, where now a measurement at y is valid if any one of the R terms u y   r  is valid. To preserve the network&#39;s invariance to object albedo, each individual u y   r  is augmented during training (before concatenation). The augmentation consists of a random scale factor, followed by adding simulated sensor noise. 
     Reducing Sparse Point Uncertainty by Incorporating Additional Sensors 
     The diffraction-based point positioning system described so far has two limitations: (1) uncertainty in the horizontal coordinates is typically higher than the vertical one; and (2) ambiguities in position estimates for certain source alignments (vertically overlapping sources). Two variations of the hardware system that mitigate these limitations, for the sparse point case will now be described. 
     Two hardware augmentations are presented: (1) using two line sensors with a diffraction grating and a horizontal cylindrical lens, respectively; and (2) using two line sensors with a double-axis diffraction gratings. 
     Hardware Augmentation 1: Horizontal Cylindrical Lens 
     To improve the accuracy of x v , an additional horizontal line scan camera can be positioned such as to image the zero-order component through a cylindrical lens, as shown in  FIG.  6 B . The cylindrical lens focuses scene points onto vertical lines which intersect the horizontal line sensor. These white streaks yield precise measurement for coordinates x v . However, unlike the rainbow streaks, the white streaks provide no information about y v . Nevertheless, this data can be used to improve 2D positioning. 
     To accomplish this, let x k   v , denote the recovered cylindrical lens coordinates, indexed by k=1, 2, . . . , K. For every recovered point (x n   v , y n   v ) the distance of its x n   v  coordinates to all x k   v  are computed. If for any k, x k   v  is within a predefined threshold distance to x n   v  (e.g., below four pixels in one embodiment), it is assumed that measurement k originated from point n, and thus replace x n   v ↔x k   v . On the other hand, if no close match is found for x n   v , (x n   v ,y n   v ) is discarded as an inaccurate measurement. 
       FIG.  7    shows an example of this result. A computer fan mounted with two LEDs at different radii is rotated at 780 RPM. A diffraction-cylindrical system images the fan at 2220 FPS.  FIG.  7 ( a )  is a spectral space-time diagram for the diffraction camera where columns show the 1D signal at different times (rows).  FIG.  7 ( b )  is a space-time diagram for the additional cylindrical lens camera.  FIG.  7 ( c )  shows recovered source positions and trajectories at three times within the first revolution of the fan. The top row of  FIG.  7 ( c )  shows results using only the diffraction signal, while the bottom row of  FIG.  7 ( c )  shows the improvement by incorporating cylindrical measurements. 
     Hardware Augmentation #2: Double-Axis Diffraction 
     The transmissive grating previously discussed has spatial variation in only one axis, thus it diffracts light along that axis. A double-axis grating has spatial variation in both axes (e.g., a star filter) and thus creates diffraction orders in both horizontal and vertical directions. Using a double-axis grating allows for a straightforward generalization of diffraction-based 2D positioning. 
     The single-axis grating, shown in  FIG.  6 A  can be replaced with a double-axis grating and an additional line sensor camera may be added to image the vertical first-order diffraction component as well, as shown in  FIG.  6 C . The vertical diffraction component creates a vertical rainbow streak on the second line sensor camera (which now requires no tilting). Each line sensor now provides a separate measurement for (x v ,y v ) using the diffraction-based 2D positioning process described above. The process of merging the pair of measurements is straightforward. If two measurements fall within a pre-defined distance E on the virtual image plane, they are assumed to originate from the same scene point and the two measurements are merged by taking the x v  coordinate from the sensor sensing the vertical rainbow streak, and taking the y v  coordinate from the horizontal sensor. 
       FIG.  8    shows an example result.  FIG.  8 ( a )  shows a double-axis diffraction system with 3 ROIs.  FIG.  8 ( b ) , shows the signal from the vertical diffraction (cam  1 ).  FIG.  8 ( c )  shows recovered points from the vertical diffraction superimposed on a 2D helper camera (ground truth, cam  3 ). Here, the three top points share the same y and thus are not fully recovered.  FIG.  8 ( d )  shows the signal from the horizontal diffraction (cam  2 ).  FIG.  8 ( e )  shows recovered points from the horizontal diffraction. Here the right vertical points share the same x and thus are not fully recovered.  FIG.  8 ( f )  shows the merging of the double-axis measurements, yielding correct positions for all points. 
     Hardware and Calibration 
     A prototype of the invention is shown in  FIG.  8 ( a ) . The setup consists of three-color cameras. Camera  1  and  2  capture the vertical and horizontal diffraction streaks, respectively. Camera  3  is the helper camera, used for gathering ground-truth training data and visualizations. Camera  3  is also used in conjunction with the horizontal cylindrical lens embodiment. Preferably, the cameras support up to 64 ROIs, thus effectively acting as line scan cameras. At least two rows are required per ROI to yield color measurements (due to Bayer filter). Camera  1  is rotated by 90° to obtain vertical 1D measurements. Cameras  1  and  2  were stripped of their built-in IR filters and mounted with 1.5 MP 9 mm lenses. Camera  3  was mounted with an 8 mm M12 lens using an adapter. The transmission gratings, in one embodiment, are 50 mm, 1200 grooves/mm. In the double-axis diffraction embodiment, all cameras viewed the scene through two stacked gratings, where one is rotated 90° with respect to the other. Otherwise, a single grating was used. 
     The motion-capture suit and glove prototypes utilize 5 mm 12V white diffuse LEDs. A second glove was fitted with 20 mm 3.3V white LEDs. For 3D scanning, a short throw DLP projector was used. 
     Calibration is performed using camera  3 , whose image plane serves as the virtual image plane shown in  FIG.  4   . For sparse points, approximately 40K samples of a single moving source using all cameras were simultaneously recorded 
     For sources, three types of LEDs (two clear and one diffuse) were used. This procedure yields the training set of Eq. (12) and is used to compute H −1  and G. For 3D scanning, a vertical line on a white wall was swept and imaged. For each projected line, its line equation y v =cx v +d is computed and used to determine the ground truth disparity: 
     
       
         
           
             
               
                 
                   
                     
                       x 
                       
                         v 
                         , 
                         GT 
                       
                     
                     ⁡ 
                     
                       ( 
                       y 
                       ) 
                     
                   
                   = 
                   
                     
                       [ 
                       
                         
                           G 
                           ⁡ 
                           
                             ( 
                             y 
                             ) 
                           
                         
                         - 
                         d 
                       
                       ] 
                     
                     c 
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     Together, all lines yield almost 700K RGBy/position samples. Intrinsic and extrinsic calibration of the helper camera and projector is performed using a checkerboard. 
     Deconvolving Diffraction for Fast Imaging of Sparse Scenes 
     The first embodiment previously discussed uses a diffraction grating to spread incident light and encode positions using the spectrum with a line sensor. Diffraction gratings are easily available and have been used for artistic effect, spectroscopy, multi-spectral sensing, and rainbow particle velocimetry. The first embodiment does not require special modulation of sources but instead relies on purely passive observation of sparse scenes and works with traditional CMOS sensors. However, the first embodiment has three important limitations: (a) it focuses only on a line diffraction pattern; (b) its reconstruction works for points that are very sparse and fails in many situations (e.g., points on a uniform grid) where multiple scene points contribute different spectral intensities to the same pixel; and (c) it requires multiple cameras to resolve the position of points reliably. 
     The first embodiment can be generalized to arbitrary diffraction patterns and arbitrary sensor sampling of the sparse scene point distributions, while using only a single camera. The image formation can be modeled as a spatially-invariant convolution of the diffraction pattern and the sparse scene. The point-spread function (PSF) of the diffraction pattern can be measured for any type of scene point or source. The image is acquired using an arbitrary configuration of ROIs and is deconvolved by enforcing sparsity. 
     Disclosed herein is a second embodiment comprising a model for recoverability of the sparse scene and an approach to minimize the ROI configuration (and maximize frame rate) based on the diffraction PSF. This embodiment is then extended to scenes with multiple spectra and the sparse scene is estimated using a known dictionary of PSFs. In contrast to the first embodiment, the second embodiment does not require any training data, does not require multiple cameras, applies to more dense scenes, and is not limited to line diffraction gratings/sensors or particular sparse scene configurations. 
     The second embodiment is applicable in several domains, such as motion capture of fast-moving markers and particle image velocimetry (PIV). The system is composed of an off-the-shelf diffraction grating mounted in front of a single, ordinary 120 Hz camera and is able to accurately capture high speed motions (kHz) such as skipping ropes, bullets from guns, complex fluid vortices, and even sparklers. 
     While this disclosure primarily demonstrates 2D position estimates of the sparse scene points at high speeds, the second embodiment can be used with multiple cameras and standard 3D estimation pipelines. 
     As an example, consider an application to directly recover points whose signal is linearly mixed on the line sensor (multiple overlaying points, e.g., light sources that lie on the same row, as shown in  FIG.  9   ). Different recovery methods which rely on a known system point spread function (PSF), such as non-blind deconvolution or convolutional sparse coding, can be used. Handling denser scenes can open the door for additional applications such as fast particle image velocimetry (See  FIG.  2 C ). 
       FIG.  9    shows dense scene deconvolution. A scene with 400 points is imaged as shown in  FIG.  9 ( a ) . The ground truth point location is marked by white squares.  FIG.  9 ( b )  shows images on the 2D image plane showing the rainbow streaks generated by the points. Here, the density causes many points to lie on the same row, yielding color sensor mixtures.  FIG.  9 ( c )  shows signals captured using five ROIs, each having a width of 10 pixels. The recovered point positions from  FIG.  9 ( c )  are superimposed in  FIG.  9 ( a )  as red squares.  FIG.  9 ( d )  shows a 2D image sensor plane reconstructed from the recovered points. 
       FIG.  10    shows the process for performing the dense scene deconvolution shown in  FIG.  9   . The 2D image plane in  FIG.  9 ( c )  is obtained by performing a convolution between the dense point field x of  FIG.  9 ( a )  and the point spread function k in  FIG.  9 ( b ) . The 2D image plane is then sub-sampled for multiple regions of interest to obtain measured signal/shown in  FIG.  9 ( c ) . 
     Convolutional Image Formation Model 
     A diffraction grating is a thin optical element that disperses light as a function of wavelength. When illuminated with a wide spectrum (e.g., white light), the periodic micro-structure of a diffraction grating produces rainbow patterns similar to the effect that a prism has on light. Simple periodic structures result in simple patterns that consist of lines, such as the single rainbow streak or a rainbow-colored star, as shown in  FIG.  11 (B) . Presented below is a novel convolutional image formation model for a sparse scene when viewed through such diffractive optics. The model can also support diffractive optical elements (DOEs), which can produce almost any desired light distribution by carefully designing the micro-structure of the optical element. 
       FIG.  11    shows the convolutional image formation and recovery. In  FIG.  11 ( a ) , an image of a 4×4 LED matrix is captured using a standard 2D camera.  FIG.  11 ( b )  shows that, when the image is imaged through a double-axis diffraction grating (inset), a single LED&#39;s point spread function (PSF) has the shape of a rainbow-colored star spanning a large part of the image domain.  FIG.  11 ( c )  shows the LED matrix of FIG.  11 ( a ) imaged through the diffraction star grating. The resulting image can be modeled as a convolution of the standard 2D frame in  FIG.  11 ( a )  with the individual LED PSF in  FIG.  11 ( b ) .  FIG.  11 ( d )  shows the 2D image sub-sampled by taking eight narrow slices of size 2056×2 regions of interest (marked by yellow lines) to yield this sub-sampled image.  FIG.  11 ( e )  shows the sub-sampled image of  FIG.  11 ( c )  vertically stretched for visualization. In  FIG.  11 ( f ) , the sub-sampled image is used with the LED PSF to yield a reconstruction of the full 2D frame. 
     Let δ(s) denote the image of a single bright scene point mapped to a single pixel s. Note that δ(s) is an intensity-only image of the point, where the intensity value at s can take arbitrary values. If the scene consists of multiple points that emit the same normalized spectrum of light, the resulting image, shown in  FIG.  11 ( a ) , is given by: 
     
       
         
           
             
               
                 
                   I 
                   = 
                   
                     
                       ∑ 
                       s 
                     
                     ⁢ 
                     
                       δ 
                       ⁡ 
                       
                         ( 
                         s 
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     When placed in front of a camera, the diffraction grating blurs the image with a rainbow-colored point spread function (PSF). If s is the center of the frame, the result of imaging δ(s) is the PSF itself, as example of which is shown in  FIG.  11 ( b ) . When viewing the scene through the diffraction grating with a plurality of points, the image formation process can be modeled as a convolution between I and the PSF h σ   diff  of the diffraction grating:
 
 I   σ   diff   =I*h   σ   diff    (18)
 
     where σ represents the different color channels of a camera. As shown in  FIG.  11 ( c ) , the resulting “diffraction” image consists of a superposition of multiple copies of the diffraction PSF, shifted to the scene point positions and scaled by the intensity of each point. As a comparison, the first embodiment assumed a scene with no such superposition, limiting the possible scenarios to which the first embodiment may be applied. 
     Capturing a full-frame video stream with a conventional sensor at high frame rates is limited by bandwidth. Instead, the approach of the second embodiment is to recover a full frame image/from a sparse set of measurements on the sensor, modeled as:
 
 I   σ   sparse   =W   σ   ⊙I   σ   diff   =W   σ ⊙( I*h   σ   diff )   (19)
 
     where: 
     the operator ⊙ denotes an element-wise Hadamard product; and the matrix is a binary sampling matrix with ones for sampled pixels and zeros otherwise. 
     Note that W σ  can vary per color channel σ (e.g., to encode the Bayer color filter array associated with most cameras). Although both I σ   sparse  and I are represented here as having the same size, the image I σ   sparse  requires far less bandwidth because W σ  is sparse, and only a small subset of camera pixels are needed to populate I σ   sparse . 
     For example,  FIGS.  11   ( c - d ) sub-samples 16 of 1542 rows, reducing the bandwidth of the video stream by a factor of nearly 100. 
     Recovering the Sparse Scene 
     I can now be recovered from the sparsely sampled sensor measurements I σ   sparse . Here, it is assumed that the PSF h σ   diff  is known and may be acquired through a calibration procedure. Image I encodes two pieces of information about each scene point: intensity and position on the image plane. Some of the applications, for example, motion capture with markers and PIV, are concerned with only the point positions. In this section, a general method for recovering I is first provided, followed by a discussion about the conditions for which the approach can yield precise point positioning. 
     In the general case, a plurality of scene points yield a mixture of point PSFs on the sensor. Coupled with image noise, this makes a direct deconvolution of Eq. (19) an under-determined problem. Therefore, the reconstruction can be formulated as an optimization: 
     
       
         
           
             
               
                 
                   
                     
                       
                         arg 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         min 
                       
                       I 
                     
                     ⁢ 
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         σ 
                       
                       ⁢ 
                       
                         
                            
                           
                             
                               
                                 W 
                                 σ 
                               
                               ⊙ 
                               
                                 ( 
                                 
                                   I 
                                   * 
                                   
                                     h 
                                     σ 
                                     diff 
                                   
                                 
                                 ) 
                               
                             
                             - 
                             
                               I 
                               σ 
                               sparse 
                             
                           
                            
                         
                         2 
                         2 
                       
                     
                   
                   + 
                   
                     γ 
                     ⁢ 
                     
                       
                          
                         I 
                          
                       
                       1 
                     
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
           
         
       
     
     where, the first term is the data fidelity term, the second term is a sparsity enforcing regularization term, and the scalar controls the contribution of the regularizer. This is a convex optimization problem that can be solved with readily available software packages. 
     Extension to Multiple Spectra 
     It was previously assumed that all scene points share the same reflectance spectrum. This allows the image formation model to be the convolution between an intensity-only image and a diffraction PSF, as shown in Eq. (18). However, general scenes can contain a variety of sources, each emitting or reflecting a unique spectrum. The image formation can therefore be extended model to handle multiple spectra and generalize the corresponding reconstruction procedure. 
     Suppose the reflectance spectrum of each scene point can be described as a linear combination of K reflectance spectra. For example, the scene shown in  FIG.  12 ( a )  contains LEDs of four different colors, each having a unique spectrum. Let h σ,k   diff  denote the PSF of a scene point having reflectance spectra k=1, 2, . . . , K. Then the resulting sub-sampled image can be modeled as: 
     
       
         
           
             
               
                 
                   
                     I 
                     σ 
                     sparse 
                   
                   = 
                   
                     
                       W 
                       σ 
                     
                     ⊙ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                       ⁢ 
                       
                         ( 
                         
                           
                             I 
                             k 
                           
                           * 
                           
                             h 
                             
                               σ 
                               , 
                               k 
                             
                             diff 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     The image I k  denotes the contribution that spectrum k has to the signal at every pixel. In other words, each camera pixel can be dominated by a single spectrum coefficient k or be a mixture of several coefficients. The signal can then be reconstructed by optimizing the following: 
     
       
         
           
             
               
                 
                   
                     
                       
                         arg 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         min 
                       
                       
                         
                           I 
                           1 
                         
                         , 
                         … 
                         , 
                         
                           I 
                           K 
                         
                       
                     
                     ⁢ 
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         σ 
                       
                       ⁢ 
                       
                         
                            
                           
                             
                               
                                 W 
                                 σ 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     ∑ 
                                     
                                       k 
                                       = 
                                       1 
                                     
                                     K 
                                   
                                   ⁢ 
                                   
                                     
                                       I 
                                       k 
                                     
                                     * 
                                     
                                       h 
                                       
                                         σ 
                                         , 
                                         k 
                                       
                                       diff 
                                     
                                   
                                 
                                 ) 
                               
                             
                             - 
                             
                               I 
                               σ 
                               sparse 
                             
                           
                            
                         
                         2 
                         2 
                       
                     
                   
                   + 
                   
                     γ 
                     ⁢ 
                     
                       
                         ∑ 
                         
                           k 
                           = 
                           1 
                         
                         K 
                       
                       ⁢ 
                       
                         
                            
                           
                             I 
                             k 
                           
                            
                         
                         1 
                       
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
     An experiment showing a scene having K=4 sources with measured PSFs is shown in  FIG.  12   . 
       FIG.  12    shows the recovery of sources of multiple spectra.  FIG.  12 ( a )  shows wherein a scene with 16 LEDs of four different colors is imaged using a double-axis diffraction grating. In  FIG.  12 ( b ) , the 2D sensor plane is sampled using 14 two-pixel wide regions of interest (ROIs) illustrated by the yellow horizontal lines.  FIG.  12 ( c )  shows that each LED yields a star PSF in the image of  FIG.  11 ( b )  whose spectrum depends on the LED color. A close-up of a single diffraction mode of the four scene spectra is shown on the right.  FIG.  11 ( d )  shows scene recovery using Eq. (20) and a single PSF.  FIG.  11 ( e )  shows scene recovery using a dictionary of the four scene PSFs with Eq. (22). 
     Note that the single-dictionary recovery yields degraded results, while using multiple dictionaries improves the recovery. A by-product of Eq. (22) are the individual coefficient weights of each PSF, visualized using the appropriate LED color on the right. 
     Domain of Position Recoverability 
     The above reconstruction algorithm handles an arbitrary sampling matrix, as well as an arbitrary PSF. However, not all PSF and sampling matrix combinations yield the same accuracy when recovering point positions from I. How the properties of the PSF constrain the sensor sampling choice is discussed below. 
     Consider a scene containing a single point mapping to pixel s. The point&#39;s image plane position can be uniquely recovered for an arbitrary PSF and sensor sampling matrix over an image domain, if there exists an injective function ƒ that maps the sampled and shifted PSF to its positional value s in this domain:
 
ƒ( W   σ ⊙[δ( s )* h   σ   diff ])   s,∇   s ∈Ω   (23)
 
       FIG.  13    shows the effect of PSF and sensor sampling on position recovery. The image domain Ω for which point positions can be recovered uniquely depends on the PSF and the sensor sampling matrix. In  FIG.  13 ( a ) , a vertical one-pixel wide rainbow PSF sampled using a single camera pixel defines a vertical recoverable domain with the same spatial support as the non-zero PSF pixels.  FIG.  13 ( b )  shows that by increasing the PSF width to three pixels, a single pixel sampling can no longer determine s precisely anywhere in the image domain. This is because the vertical (y-axis) position can still be determined with certainty, but now a two-pixel ambiguity exists in the horizontal position.  FIG.  13 ( c )  shows that sampling the entire row resolves the horizontal position ambiguity, showing that different W σ  yield different Ω for the same PSF.  FIG.  13 ( d )  shows that the PSF need not be a single rainbow line. In this example, the PSF shift can be computed using the ratios between the different RGB channels as in  FIGS.  13   ( a - c ). However, here these measurements are spatially separated across different columns of the sampling row. 
       FIG.  13    also illustrates various cases of the dependence of Ω on the PSF and W σ .  FIGS.  13 ( b, c )  illustrate that the same PSF can yield an empty uniquely-recoverable domain or a uniquely-recoverable domain that spans the entire image plane. Specifically, the single pixel sampling in  FIG.  13 ( b )  can yield point positioning with an ambiguity of two horizontal pixels, while the line sampling shown in  FIG.  13 ( c )  yields no horizontal or vertical position ambiguity. A PSF and a sampling matrix that maximize the support of Ω. Eq. (23) states an implicit condition for defining Ω, but does not provide a method for computing it given a PSF and Next, a method for estimating Ω for a subset of PSFs that fulfill a stronger condition is provided, the condition being that the PSF image itself contains individual pixels that define an injective function from their intensity-scaled color channel (e.g. RGB) values to their position. The method is illustrated in  FIG.  14   , which shows estimating the uniquely recoverable domain. Given a PSF and a sensor sampling matrix W σ , the uniquely-recoverable domain can be approximated by the following procedure: First, estimate the PSF&#39;s injective pixel domain h inj  as shown in  FIG.  14 ( a ) . Next, correlate h inj  with W σ  to yield {tilde over (Ω)}.  FIG.  14 ( b )  illustrates this process for choosing the ROIs&#39; pitch P for a uniform horizontal sampling of the image plane. 
     Let h inj  denote a binary image with a domain Y over which there exists an injective function g:
 
 g (α h   σ   diff [ x ])   x,∇   x ∈Γ   (24)
 
     Here, x is an image pixel, and a is an arbitrary scaling factor representing the intensity of the system&#39;s PSF. For example, the injective PSF domain of a thin diffraction rainbow line, shown in  FIG.  13 ( a ) , is approximately an image line where the camera response function yields unique color channel values in response to the incident light spectrum. Thus, given h inj , the recoverable domain can be estimated by computing its image domain support using W σ :
 
{tilde over (Ω)}{ x |( W   σ   *h   inj )[ x ]&gt;0}   (25)
 
     where “*” denotes the correlation operator. Note that Eq. (25) constitutes a sufficient but not necessary condition for Ω and thus {tilde over (Ω)}⊆Ω. As shown in  FIG.  13 ( c ) , a PSF may have Γ=∅ (because no single pixel in the PSF directly yields a precise shift position) yet still yield unique recoverability over the entire image domain. 
     Choosing W σ  Using Camera ROIs 
     Many machine vision cameras allow defining multiple regions of interests (ROIs) for which to perform sensor read-out. The output image is then a concatenation of the ROIs of choice, as shown in  FIGS.  11   ( c - d ). A ROI is defined as a rectangle in the image plane. As most cameras perform readout sequentially one row at a time, the capture speed usually depends on the total number of rows in all the defined ROIs. 
     Given these hardware constraints, the sampling matrix W σ  consists of multiple ROIs of size 2056×M, where M is the vertical width of each ROI. The multiple ROIs are repeated vertically every P pixels. Thus, the image sub-sampling factor is P/M. Higher sub-sampling factors result in faster frame rates. 
     Maximizing the sub-sampling factor requires maximizing P while minimizing M. This can be achieved by setting M to its minimal allowable value (e.g. 2 pixels for some cameras) while finding the highest P=P max  such that Ω covers the entire image domain (field-of-view). To find P max , the PSF&#39;s injective domain h inj  is estimated and Eq. (25) is used to search for the maximal pitch P max . 
     The diffraction line imaging of the first embodiment has been generalized to handle an arbitrary sensor sampling combined with an arbitrary diffractive element. A system implementing the second embodiment is easily assembled, requiring just a single diffractive optical element mounted in front of one camera. Calibration is fast and straightforward, requiring a single capture of the scene&#39;s PSF. The reconstruction algorithm requires no prior data-set for training. There is a trade-off between the different choices of PSF, created by using different diffractive optics. A single-axis diffraction grating is energy efficient because most of the light&#39;s power is concerted in a single diffraction mode (i.e., the rainbow streak or line). However, the single-axis PSF of the first embodiment yields higher positioning uncertainty in the vertical axis, and has a limited injective domain, which reduces the sub-sampling factor. Conversely, the “star” PSF shown in  FIG.  11 ( b )  and  FIG.  12 ( c )  has a large spatial support and contains rainbow streaks in both axes. It therefore provides accurate localization in both spatial coordinates and enables high sensor sub-sampling factors, typically upwards of 100. However, its spatial spread comes at the cost of light efficiency, making it mostly useful for high-speed imaging of relatively bright sources such as LEDs and car headlights. 
     Diffraction-based imaging assumes sparsity in the scene points or sources. As the scene becomes dense a lower sub-sampling factor is required for accurate point positioning. A higher sensor sampling increases the signal, making the optimization more robust to noise and signal saturation. In some embodiments, 8 ROIs were used. The sub-sampling factor was decreased by increasing the width of each ROI while still maintaining the desired capture speed. 
     The performance of the second embodiment depends on the system&#39;s hardware components: camera, diffractive optics, and illumination source. Recent experimental cameras that allow an arbitrary sampling of the sensor plane may yield higher performance by designing more efficient sampling matrices. Conversely, improved performance can be obtained by designing custom diffractive optical elements that yield a PSF with superior performance. For instance, the star PSF used herein can be improved by concentrating most of the energy in just two diffraction modes (rainbow lines): one horizontal line and one vertical line, with minimal energy along the rest of the directions. Finally, incorporating a learning-based component to position recovery may further improve point positioning by learning and denoising the PSF-specific artifacts that result in the reconstructed frames. In conclusion, the second embodiment has the potential to allow high-speed capture of sparse scenes in various applications. 
     EXTENSIONS AND IMPROVEMENTS 
     The underlying principle for diffraction line sensing can be further improved by using different imaging systems. The principle can be applied to light sources at non-visible wavelengths coupled with the appropriate diffractive optics and cameras. Using prisms instead of diffraction gratings can improve light efficiency for the single-diffraction system. Moreover, designing custom diffractive elements whose PSF is more sensitive to horizontal motion may yield better accuracy for recovering the virtual x-coordinates. 
     Diffraction line imaging is a novel imaging principle relying on light diffraction and 1D sensing. Using the principle, proof of concepts for applications such as motion capture of uncontrolled and unmodulated LED lights, tracking car headlights and structured light 3D line scanning have been demonstrated. A key advantage of using line sensors is speed. Speed is crucial for motion capture because it greatly eases tracking. For 3D scanning, speed is essential when scanning moving objects (e.g., an industrial conveyor belt). The invention uses multiple rows from a conventional 2D sensor by configuring the ROI to mimic a line sensor, resulting in fast readout rates of up to 2220 FPS. Higher-speed line sensors could increase the frame rate up to 45,000 FPS (×20 Faster than the current embodiments) and improve light efficiency with large pixels up to 4×32 microns in size (×10 larger). 
     Conversely, using 2D sensors with multiple ROIs allows for smooth control over the speed vs. quality trade-off, namely, additional ROIs reduce speed but increase accuracy. As with any sensor, performance is dependent on the available signal-to-noise ratio (SNR). The SNR of systems depends on the sensor&#39;s spectral response. 
     Thus, hyper-spectral sensors whose pixel&#39;s spectral response is less limited than that of the standard RGB sensor can boost performance. For 3D scanning, bright broadband sources, such as LEDs or supercontinuum/swept-frequency lasers can additionally increase SNR many-fold. Learning-based or dictionary-based approaches may be employed to improve reconstruction quality by extracting multiple vertically overlapping points, which yield a linear color mixture on the 1D line scan sensor. 
     The calculations required to implement the above-described invention can be accomplished using a processor coupled with one or more light sensors and executing software comprising instructions to read the necessary data from the one or more light sensors and calculate the results in accordance with the equations contained in the description of the invention. Specific implementations of the software are not provided herein as any architecture or organization for the software producing the desired results can be used. 
     In another embodiment of the invention, the methods disclosed herein may be used to detect vibrational modes of one or more objects-of-interest in the scene. The row delay in the rolling shutter of a camera may be used as an alternative way for generating a burst of very high-speed samples. In a rolling shutter sensor, the sensor rows are exposed one after the other with a very short delay. Thus, when capturing a scene having sparse sources through a diffraction grating with a rolling shutter sensor, the motion of the sources will be encoded in the individual sensor rows of a single rolling shutter frame. This setup can be used for vibration detection by imaging one or more retro-reflective markers, placed on objects-of-interest in the scene. The vibration of each marker can be detected by observing the deformation of its respective point spread function. Specifically, the imaged point spread function will exhibit a periodic deformation which encodes the frequency and amplitude of vibrations of the object-of-interest to which the retro-reflective marker that generated the point spread function is attached. 
     The deformation of the point spread function can be optically magnified with cylindrical optics, wherein each image axis is optically magnified by a different factor. For example, for a vertical diffraction streak point spread function, the horizontal image dimension can be magnified to increase the PSF deformation along the horizontal axis while keeping the vertical axis magnification as is. To acquire the deformation of multiple scene markers using the non-uniform magnification described above, additional optics can be used to overlay the magnified point-spread functions from different markers onto the single 2d image plane, yielding a mixture of the magnified deformed PSFs in each frame. 
     The invention has been explained in terms of specific embodiments, including specific methods, components and configurations of components. As may be realized, the scope of the invention is not meant to be limited by the exemplary embodiments described herein but is set forth in the following claims. Further, although the invention has been explained largely in terms of using diffraction gratings, it should be realized that prisms could be used in lieu of diffraction gratings to generate the component spectrums of the light sources.