Abstract:
An apparatus and method for image processing of XYZ point clouds obtained from a GmAPD LADAR using low-pass filtering followed by high-pass filtering and deconvolution. Preferably, the low-pass filter parameters are developed numerically utilizing Weiner-Levinson Deconvolution (WLD).

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
       [0001]    This application is a non-provisional of U.S. patent application Ser. No. 61/511,004, filed Jul. 22, 2011, the disclosure of which is incorporated by reference herein in its entirety. 
     
    
     BACKGROUND 
       [0002]    This disclosure relates generally to the field of imaging and more particularly to enhancing images obtained from Geiger mode Avalanche Photo Diode detectors using three-dimensional statistical differencing. 
         [0003]    Imaging sensors such as laser radar sensors (LADARs) acquire point clouds of a scene. The point clouds of the scene are then image processed to generate three dimensional (3D) models of the actual environment of the scene. The image processing of the 3D models enhances the visualization and interpretation of the scene. Typical applications include surface measurements in airborne and ground-based industrial, commercial and military scanning applications such as site surveillance, terrain mapping, reconnaissance, bathymetry, autonomous control navigation and collision avoidance and the detection, ranging and recognition of remote military targets. 
         [0004]    Presently there exist many types of LADARs for acquiring point clouds of a scene. A point cloud acquired by a LADAR typically comprise x, y &amp; z data points from which range to target, two spatial angular measurements and strength (i.e., intensity) may be computed. However, the origins of many of the individual data points in the point cloud are indistinguishable from one another. As a result, most computations employed to generate the 3D models treat all of the points in the point cloud the same, thereby resulting in indistinguishable “humps/bumps” on the 3D surface model of the scene. 
         [0005]    Various imaging processing techniques have been employed to reconstruct the blurred image of the scene. The blurring or convolution of the image is a result of the low resolution (i.e., the number of pixels/unit area) of the intensity images at longer distances and of distortion of the intensity image by the LADAR optics and by data processing. Accordingly, the image must be de-blurred (deconvolved). 
         [0006]    Relevant herein, LADARs may comprise arrays of avalanche photodiode (APD) detectors operating in Geiger-mode (hereinafter “GmAPD”) that are capable of detecting a single photons incident onto one of the detectors.  FIG. 1  diagrammatically depicts a typical GmAPD LADAR  10  including focal plane arrays  12  of avalanche photodiode (APD) detectors  14  operating in Geiger-mode. Integrated timing and readout circuitry (not shown) is provided for each detector  14 . In typical operation, a laser pulse emitted from a microchip laser  16  passes through a band-pass filter  18 , variable divergence optics  20 , a half-wave plate  22 , a polarizing beam splitter  24 , and is then directed via mirrors  26  and  28  through a beam expander  30  and a quarter wave plate  32 . Scanning mirrors  34  then steer the laser pulses to scan the scene  36  of interest. It is noted that the scanning mirrors  34  may allow the imaging of large areas from a single angle of incidence or small areas imaged from a variety of angles on a single pass. Return reflections of the pulse from objects in the scene  36  (e.g., tree and tank) pass in the opposite direction through the polarizing beam splitter  24 , a narrow band filter  38 , and then through a zoom lens  40  onto the detector array  12 . The outputs of the detector array  12  forming a point cloud  42  of XYZ data are then provided to an image processor  44  for viewing on a display  46 . 
         [0007]    More particularly, the operation of a GmAPD LADAR occurs as follows. After the transmit laser pulse leaves the GmAPD LADAR, the detectors  14  are over-biased into Geiger-mode for a short time, corresponding to the expected time of arrival of the return pulse. The window in time when the GmAPD is armed to receive the return pulse is known as the range gate. During the range gate, the GmAPD and its integrated readout circuitry is sensitive to single photons. The high quantum efficiency in the GmAPD results in a high probability of generating a photoelectron. The few volts of overbias ensure that each free electron has a high probability of creating the growing avalanche which produces the volt-level pulse that is detected by the CMOS readout circuitry. This operation is more particularly described in U.S. Pat. No. 7,301,608, the disclosure of which is hereby incorporated by reference herein. 
         [0008]    Unfortunately, during photon detection, the GmAPD does not distinguish among free electrons generated from laser pulses, background light, and thermal excitations within the absorber region (dark counts). High background and dark count rates are directly detrimental because they introduce noise (see, e.g., FIG. 7 of U.S. Pat. No. 7,301,608) and are indirectly detrimental because they reduce the effective sensitivity to signal photons that arrive later in the range gate. See generally, M. Albota, “Three-dimensional imaging laser radar with a photon-counting avalanche photodiode array and microstrip laser”, Applied Optics, Vol. 41, No. 36, Dec. 20, 2002, the disclosure of which is hereby incorporated by reference herein. Nevertheless, single photon counting GmAPDs are favored due to efficient use of the power-aperture. 
         [0009]    There presently exist several techniques for extracting the desired signal from the noise in a point cloud acquired by a GmAPD LADAR. Representative techniques include Z-Coincidence Processing (ZCP) that counts the number of points in fixed-size voxels to determine if a single return point is noise or a true return, Neighborhood Coincidence Processing (NCP) that considers points in neighboring voxels, and various hybrids thereof (NCP/ZCP). See P. Ramaswami, “Coincidence Processing of Geiger-Mode 3D Laser Radar Data”, Optical, Society of America, 2006, the disclosure of which is hereby incorporated by reference herein. 
         [0010]    In addition to removal of noise from a point cloud through the use of NCP or ZCP techniques, it is often desirable to enhance the resulting image. Prior art image enhancement techniques include un-sharp masking techniques using a high-pass filter, techniques for emphasizing medium-contrast details more than large-contrast details using adaptive filters and statistical differential techniques that provide high enhancement in edges while presenting a low effect on homogenous areas. 
       SUMMARY 
       [0011]    According to one embodiment, a method for processing XYZ point cloud of a scene acquired by a GmAPD LADAR is disclosed. The method of this embodiment includes: applying low-pass filtering utilizing Deconvolution to the XYZ point cloud to produce a D point cloud; and displaying an image of the D point cloud. 
         [0012]    According to another embodiment, a method for processing a XYZ point cloud of a scene acquired by a GmAPD LADAR that includes: Z-clipping the XYZ point cloud adaptive histogramming to produce a Z-clipped point cloud; applying low-pass filtering utilizing Weiner-Levinson Deconvolution to the XYZ point cloud, utilizing a Deconvolution Matrix having at least one parameter that is operator-selectable to produce a WLD point cloud; thresholding the WSD point cloud to produce a first thresholded point cloud; sharpening the WLD point cloud in the X-Y plane by highpass filtering to produce a sharpened point cloud; thresholding the sharpened point cloud to produce a second thresholded point cloud; mitigating timing uncertainty in the second thresholded point cloud by deconvolving the second thresholded point cloud in the vertical direction to produce a deconvolved point cloud; thresholding and cleansing the deconvolved point cloud in the vertical direction to produce a thresholded/cleansed point cloud; and displaying an image of the thresholded/cleansed point cloud by counting photons at points in the thresholded/cleansed point cloud is disclosed. 
         [0013]    According to another embodiment, a system for processing a XYZ point cloud of a scene acquired by a GmAPD LADAR is disclosed. The system includes an image processor that performs low-pass filtering utilizing Deconvolution to the XYZ point cloud to produce a D point cloud and a display for displaying an image of the D point cloud. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0014]    For a fuller understanding of the present disclosure and its advantages, reference is now made to the following description, taken in conjunction with the accompanying drawings, in which: 
           [0015]      FIG. 1  is a diagrammatic view of a typical GmAPD LADAR that may be employed by the present invention to acquire an XYZ point cloud representing the image of the scene of interest; 
           [0016]      FIG. 2  is a process flow diagram of the method of the invention implemented on an image processor for display or further processing; 
           [0017]      FIG. 3  is a diagrammatic view of adaptive histogramming; 
           [0018]      FIG. 4  is a diagrammatic view of the Sharpen (high-pass) Matrix employed in the method of the invention; and 
           [0019]      FIG. 5  is a diagrammatic view of the Deconvolution Matrix from the Weiner-Levinson Deconvolution (WLD) employed in the method of the invention. 
       
    
    
       [0020]    Similar reference characters refer to similar parts throughout the several views of the drawings. 
       DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
       [0021]    The following description is of the best mode presently contemplated for carrying out the invention. This description is not to be taken in a limiting sense, but is made merely for the purpose of describing one or more preferred embodiments of the invention. The scope of the invention should be determined with reference to the claims. 
         [0022]    The apparatus and method of the invention comprises a typical GmAPD LADAR  10  described above in connection with  FIG. 1  to acquire a point cloud  42 A of XYZ data of a scene of interest  36  that is provided to an image processor  44 . It shall be understood without departing from the spirit and scope of the invention, that neither the apparatus nor method of the invention is limited to any particular type or brand of GmAPD LADARs  10 . 
         [0023]    The image processor  44  may be embodied in a general purpose computer with a conventional operating system or may constitute a specialized computer without a conventional operating system so long as it is capable of processing the XYZ point cloud  42 A in accordance with the process flow diagram of  FIG. 2 . Further, it shall be understood without departing from the spirit and scope of the invention, that neither the apparatus nor the method of the invention is limited to any particular type or brand of image processor  44 . 
         [0024]    As shown in  FIG. 2 , a method according to one embodiment includes storing the XYZ point cloud  42 A of data into the memory of the image processor  44  at block  202 . The memory may comprise any type or form of memory. The image processor  44  may comprise a computational device such as application specific integrated circuits (ASIC), or a central processing unit (CPU), digital signal processor (DSP) or field-programmable gate arrays (FGPA) containing firmware or software, that sequentially performs the following computations on the XYZ point cloud  42 A. 
         [0025]    After being stored, the XYZ point cloud  42 A is Z-clipped based on adaptive histogramming at block  202  to form a Z-clipped point cloud  42 B. The Z-clipping performed at block  202  can include, for example, applying histogram equalization in a window sliding over the image pixel-by-pixel to transform the grey level of the central window pixel. However, to reduce the noise enhancement and distortion of the field edge, as shown in  FIG. 3 , a contrast-limited adaptive histogram equalization is preferably performed in the Z-direction to clip histograms from the contextual regions before equalization, thereby diminishing the influence of dominate grey levels. 
         [0026]    At block  204 , the reference “waveform” generated by histogramming photon return times, is used in the Weiner-Levinson Deconvolution (WLD) to “flatten” the response into an impulse and form WLD point cloud  42 c. The Deconvolution Matrix is derived as follows:
       (1) A portion of the waveform that is known to contain the impulse of interest is auto-correlated. The auto correlation Rxx(t) of a function x(t) is defined as       
 
         [0000]    
       
         
           
             
               
                 R 
                 xx 
               
                
               
                 ( 
                 t 
                 ) 
               
             
             = 
             
               
                 
                   x 
                    
                   
                     ( 
                     t 
                     ) 
                   
                 
                 ⊗ 
                 
                   x 
                    
                   
                     ( 
                     t 
                     ) 
                   
                 
               
               = 
               
                 
                   ∫ 
                   
                     - 
                     ∞ 
                   
                   ∞ 
                 
                  
                 
                   
                     x 
                      
                     
                       ( 
                       r 
                       ) 
                     
                   
                    
                   
                     x 
                      
                     
                       ( 
                       
                         l 
                         + 
                         r 
                       
                       ) 
                     
                   
                    
                   
                      
                     t 
                   
                 
               
             
           
         
       
     
         [0000]    where the symbol           denotes correlation.
 
For the discrete implementation, let Y represent a sequence whose indexing can be negative, let n be the number of elements in the input sequence x, and assume that the indexed elements of x that lie outside its range are equal to zero:
 
         [0000]        X   j =0, j&lt; 0 or j≧n
 
         [0000]    then obtain the elements of Y using: 
         [0000]    
       
         
           
             
               y 
               j 
             
             = 
             
               
                 ∑ 
                 
                   k 
                   - 
                   0 
                 
                 
                   n 
                   - 
                   1 
                 
               
                
               
                 
                   X 
                   k 
                 
                  
                 
                   X 
                   
                     j 
                     + 
                     k 
                   
                 
               
             
           
         
       
     
         [0000]    for j=−(n−1), (n−2), . . . , 2, −1, 0, 1, 2, . . . , n−1. The elements of the output sequence Rxx are related to the elements in the sequence Y by: 
         [0000]      Rxx i =y i−(n−1 ) 
         [0000]    for i=0, 1, 2, . . . , 2n−2. The number of elements in the output sequence Rxx is 2n−1. Thus, Rxx n−1 =Rxx n−1 *1.01 (to add 1% white noise to the peak). 
         [0028]    (2) An m×m matrix A is constructed from the sequence above as follows: 
         [0000]    
       
         
           
             A 
             = 
             
               [ 
               
                 
                   
                     
                       Rxx 
                       
                         n 
                         - 
                         1 
                       
                     
                   
                   
                     
                       Rxx 
                       n 
                     
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     
                       Rxx 
                       
                         n 
                         - 
                         1 
                         + 
                         m 
                       
                     
                   
                 
                 
                   
                     
                       Rxx 
                       
                         n 
                         - 
                         2 
                       
                     
                   
                   
                     
                       Rxx 
                       
                         n 
                         - 
                         1 
                       
                     
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     
                       Rxx 
                       
                         n 
                         - 
                         2 
                         + 
                         m 
                       
                     
                   
                 
                 
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                 
                 
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     … 
                   
                 
                 
                   
                     
                       Rxx 
                       
                         B 
                         - 
                         1 
                         + 
                         m 
                       
                     
                   
                   
                     
                       Rxx 
                       
                         B 
                         - 
                         m 
                       
                     
                   
                   
                     … 
                   
                   
                     … 
                   
                   
                     
                       Rxx 
                       
                         B 
                         - 
                         1 
                       
                     
                   
                 
               
               ] 
             
           
         
       
     
         [0000]    A delayed impulse vector V of length m is defined as: 
         [0000]    
       
         
           
             V 
             = 
             
               [ 
               
                 
                   
                     0 
                   
                 
                 
                   
                     1 
                   
                 
                 
                   
                     ⋮ 
                   
                 
                 
                   
                     ⋮ 
                   
                 
                 
                   
                     ⋮ 
                   
                 
                 
                   
                     0 
                   
                 
                 
                   
                     0 
                   
                 
               
               ] 
             
           
         
       
     
         [0029]    (3) The solution vector C of the Weiner Levinson coefficients is then computer using: 
         [0000]      C=A −1 V 
         [0000]    and then C is normalized to C 0 .
       (4) The solution vector C (which is a z-directional vector) is then expanded in the x and y directions as far out as necessary using a distance relationship. This results in a 3-dimensional matrix, D, the deconvolution matrix.              (5) The deconvolution matrix is then scaled, rounded and converted to integer format.   (6) The coefficient vector D is used as the coefficients of a FIR filter on the image F. That is:       
 
         [0000]    
       
      
       G=F 
               
       D.  
      
       
         
           
             where the symbol           denotes 3-dimensional cross-correlation.
 
 FIG. 5  illustrates the Deconvolution Matrix from the WLD. Notably, the voxelizing and deconvolving in three dimensions eliminates (or substantially reduces) noise and distributes energy to accommodate dispersive targets. The resulting WLD point cloud  42 C is saved in memory for further processing according to the method of the invention.
 
           
         
       
     
         [0034]    Referring again to  FIG. 2 , the resulting WLD point cloud  42 C is thresholded at block  206  to reduce processing time. The resulting thresholded point cloud  42 D is saved in memory for further processing according to the method of the invention. To reduce processing time, the thresholded point cloud  42 D is sharpened in the X-Y plane by a refocus (high-pass) matrix as illustrated in  FIG. 4  at block  208 . The resulting sharpened point cloud  42 E can then be thresholded again at block  210  (producing thresholded point cloud  42 F) to reduce additional noise around the edges of the scene thereby sharpening the image. 
         [0035]    The resulting thresholded point cloud  42 F can then be deconvolved at block  212  in the vertical Z direction { . . . , −d2, −d1, −d0, +d0, +d1, +d2, . . . } using a spiking function to mitigate timing uncertainty. The resulting deconvolved point cloud  42 G can then by thresholded and cleansed downwardly in the Z direction at block  214  to minimize processing. The result is thresholded/cleansed point cloud  42 H that represents the photons returned from the scene. 
         [0036]    At block  216 , thresholded/cleansed point cloud  42 H representing the photons returned from the scene, are counted at each point in the scene  46  and the resulting image is displayed via display  46  at block  218 . It shall be understood that in various embodiments any of the previously described point clouds could have their photons counted and be displayed. 
         [0037]    The present disclosure includes that contained in the appended claims, as well as that of the foregoing description. Although this invention has been described in its preferred form with a certain degree of particularity, it is understood that the present disclosure of the preferred form has been made only by way of example and that numerous changes in the details of construction and the combination and arrangement of parts may be resorted to without departing from the spirit and scope of the invention. 
         [0038]    Now that the invention has been described,