Patent Publication Number: US-10776111-B2

Title: Point cloud data method and apparatus

Description:
CROSS REFERENCE TO RELATED APPLICATION 
     This application claims the benefit of U.S. Provisional Application No. 62/531,495, filed Jul. 12, 2017, the disclosure of which is herein incorporated by reference in its entirety. 
    
    
     TECHNICAL FIELD 
     The present invention relates to point cloud processing, and, more particularly, to point cloud rendering for real-time point cloud data collected from a variety of sensor types. 
     BACKGROUND 
     The use of three-dimensional (3D) sensing systems and 3D imaging data is commonly used for generating 3D images of a location for use in various applications. For example, such 3D images are commonly used for generating topographical maps or for surveillance of a location, and such sensing systems typically operate by capturing elevation data associated with the location. 
     An example of one well-known 3D imaging system that generates 3D point cloud data is a Light Detection and Ranging (LIDAR) system. Typically, such LIDAR systems use a high-energy laser, optical detector, and timing circuitry to determine the distance to a target. In such LIDAR systems, data is generated by recording multiple range echoes from a single pulse of laser light to generate a frame (also referred to as an image frame). Each frame of LIDAR data is comprised of a collection of points in three dimensions (i.e., a 3D point cloud) which correspond to multiple range echoes within a sensor aperture. These points may be organized into so-called “voxels” which represent values on a regular grid in a 3D space. Voxels used in 3D imaging are akin to pixels used in a two-dimensional (2D) imaging device context. These frames may be processed to reconstruct a 3D image of a location where each point in the 3D point cloud has an individual (x, y, z) value representing the actual surface within the 3D scene under investigation. 
     As will be appreciated, LIDAR sensors collect vast amounts of data with scan rates approaching one million measurements per second. As such, these 3D sensing systems make efficient storage, data processing, and data visualization challenging given the data set sizes collected. 
     Therefore, a need exists for a point cloud rendering and point cloud filtering technique for real-time point cloud data collection from a variety of sensor types that delivers enhanced performance including reducing processing requirements, limiting local memory consumption and optimizing overall data visualization. 
     BRIEF SUMMARY OF THE EMBODIMENTS 
     In accordance with various embodiments, a point cloud rendering method and apparatus for real-time point cloud data collection from a variety of sensor types is provided that delivers enhanced performance including reducing processing requirements, limiting local memory consumption and optimizing overall data visualization. 
     More particularly, in accordance with an embodiment, the processing of large 3D point cloud data sets collected from a variety of sensors (e.g., mobile mapping systems, terrestrial laser scanners, and unmanned aircraft systems, to name just a few) is facilitated using a particular tree traversing library. In the core of the library, in accordance with the embodiment is a succinct representation of full binary tree(s) and algorithm(s) for operating with the succinct representation for tree traversing and executing so-called range queries on the binary tree. This representation allows for an efficient and generic way for assigning attributes of any type to inner and leaf nodes of full binary trees. 
     In conjunction with this full binary succinct representation, in accordance with the embodiment, a provably succinct representation and traversing algorithm is employed for arbitrary binary trees and generic storage for attributes assigned to inner and leaf nodes of binary trees. This provably succinct binary tree representation is utilized for building binary space partition trees and to apply spatial indexing, executing spatial queries, and applying nearest neighbor searches in point clouds. 
     The tree traversing library, in accordance with an embodiment, provides a way to store and to efficiently traverse extremely large (e.g., billions of nodes) binary trees, retrieve and process assigned node attribute data, and run queries; all completely out of core, without the requirement to load and keep all the data in memory. In this way, the tree traversing library, associated data structures, algorithms and/or file formats can be used generically for efficient storing and running queries on a collection of arbitrary structured data types of fixed and variable length. The tree traversing library can further be considered as a lossless controlled lossy data compression technique providing efficient arbitrary access to compressed data which does not require full decompression. This enables efficient point cloud storage with improved data compression, high performance queries, and flexibility in supporting attributes to each point(s). 
     These and other advantages of the embodiments will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows a flowchart of illustrative operations for out-of-core rendering of large point cloud data sets in accordance with an embodiment; 
         FIG. 2  shows an illustrative BSP tree constructed in accordance with an embodiment; 
         FIG. 3  shows an illustrative compacted prefix tree in accordance with the embodiment; 
         FIG. 4  shows a flowchart of illustrative operations for frustum culling using the compacted prefix tree of  FIG. 3  in accordance with an embodiment; 
         FIG. 5  shows illustrative results obtained in performing the operations set forth in  FIG. 4 ; 
         FIG. 6  is a high-level block diagram of an exemplary computer in accordance with an embodiment; 
         FIG. 7  shows a flowchart of illustrative operations for point cloud filtering in accordance with an embodiment; 
         FIG. 8  shows an illustrative point cloud filtered in accordance with the operation of  FIG. 7  in accordance with an embodiment; 
         FIG. 9  shows an illustrative example of the principal component grids having overlapping boxes for a point cloud in accordance with an embodiment; and 
         FIG. 10  shows an illustrative Gauss grid for a point cloud in accordance with an embodiment. 
     
    
    
     DETAILED DESCRIPTION 
     In accordance with various embodiments, a point cloud rendering method and apparatus for real-time point cloud data collection from a variety of sensor types is provided that delivers enhanced performance including reducing processing requirements, limiting local memory consumption and optimizing overall data visualization. 
     More particularly, in accordance with an embodiment, the processing of large 3D point cloud data collected from a variety of sensors (e.g., mobile mapping systems, terrestrial laser scanners, and unmanned aircraft systems, to name just a few) is facilitated using a particular tree traversing library. In the core of the library, in accordance with an embodiment is a succinct representation of full binary tree(s) and algorithm(s) for operating with the succinct representation for tree traversing and executing range queries on the binary tree. This representation allows for an efficient and generic way for assigning attributes of any type to inner and leaf nodes of full binary trees. 
     In conjunction with this full binary succinct representation, in accordance with an embodiment, a provably succinct representation and traversing algorithm is employed for arbitrary binary trees and generic storage for attributes assigned to inner and leaf nodes of binary trees. This provably succinct binary tree representation is utilized for building binary space partition trees and to apply spatial indexing, executing spatial queries, and applying nearest neighbor searches in point clouds. 
     As will be appreciated, succinct representations of binary trees are well-known. Existing representations allow for efficient pre-order or level-order tree traversal and provide a way to efficiently assign data attributes of the same type to tree nodes. However, known succinct implementations do not guarantee good locality of data and do not allow adding node attributes of arbitrary-type separately to inner nodes and leaves. By good locality, it is meant herein that accessing node data and traversing to nodes that are close in the tree requires access to bytes that reside close(er) in memory. On the other hand, the embodiments herein are both succinct and allow fast random traversing and separate attributes for inner-nodes and leaves of the tree. 
     The embodiments herein exploit certain succinct binary tree representations described for compact data storage and out-of-core (i.e., without loading the whole data set into memory) data rendering and processing, and facilitates a way to store and access point cloud coordinates and attribute data in a generic way allowing for an enhanced and unified application programming interface (API) for various versions of data storage formats. This further allows for efficient arbitrary tree traversals (including pre-order, in-order, and/or post-order traversals), and for the efficient assigning of different types of data attributes to inner tree nodes and to tree leaves. 
       FIG. 1  shows a flowchart of illustrative operations  100  for out-of-core rendering of large point cloud data sets (e.g., LIDAR data sets) in accordance with an embodiment. At step  105 , a set of n floating-point coordinates is received defined as:
 
{[ x   0   float   ,y   0   float   ,z   0   float ] . . . [ x   n   float   ,y   n   float   ,z   n   float ]}
 
At step  110 , the received floating-point coordinates are received and translated to an origin given by:
 
[   j   float ,   j   float ,   j   float ]=[ x   j   float   ,y   j   float   ,z   j   float ]−[ x   0   float   ,y   0   float   ,z   0   float ].
 
These coordinates are encoded, at step  115 , as integers with a specified precision p (e.g., p=1000 for 1 millimeter precision or p=100 for 1 centimeter precision) as follows:
 
[ x   j   int   ,y   j   int   ,z   j   int ]=[   j   float ,   j   float ,   j   float ]× p.  
 
     Next, at step  120 , there is an application of shift coordinates as unsigned integers given by:
 
[ x   j   uint   ,y   j   uint   ,z   j   uint ]=[ x   j   int   ,y   j   int   ,z   j   int ]+INT_MAX/2,
 
and m bits of (x, y, z) unsigned integer coordinates of each point are scrambled, at step  125 , as follows:
 
     
       
         
           
             
               { 
               
                 
                   
                     
                       [ 
                       
                         
                           x 
                           0 
                         
                         ⁢ 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         ⁢ 
                         
                           x 
                           m 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       [ 
                       
                         
                           y 
                           0 
                         
                         ⁢ 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         ⁢ 
                         
                           y 
                           m 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       [ 
                       
                         
                           z 
                           0 
                         
                         ⁢ 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         . 
                         
                             
                         
                         ⁢ 
                         
                           z 
                           m 
                         
                       
                       ] 
                     
                   
                 
               
               } 
             
             ⇒ 
             
               
                 { 
                 
                   
                     [ 
                     
                       
                         x 
                         0 
                       
                       , 
                       
                         y 
                         0 
                       
                       , 
                       
                         z 
                         0 
                       
                     
                     ] 
                   
                   ⁢ 
                   
                       
                   
                   . 
                   
                       
                   
                   . 
                   
                       
                   
                   . 
                   
                       
                   
                   ⁢ 
                   
                     [ 
                     
                       
                         x 
                         m 
                       
                       , 
                       
                         y 
                         m 
                       
                       , 
                       
                         z 
                         m 
                       
                     
                     ] 
                   
                 
                 } 
               
               . 
             
           
         
       
     
     At step  130 , the scrambled coordinates are lexicographically sorted as unsigned integers employing any number of conventional sorting algorithms to construct and output the constructed prefix tree, at steps  135  and  140 , respectively. The constructed prefix tree is a binary space partitioning (BSP) tree constructed using well-known Morton encoding in which tree nodes represent hyper-rectangular regions of space and tree branches represent hyper-planar partitions of these hyper-rectangles into halves. In this way, the constructed BSP prefix tree is one in which tree nodes represent hyper-rectangular regions of space and tree branches represent hyper-planar partitions of these hyper-rectangles into halves. 
       FIG. 2  shows an illustrative BSP tree  200  constructed in accordance with an embodiment and the aforementioned operations. As shown, BSP tree  200  is comprised of root node  205  and tree nodes  210 - 1  through  210 - 16 . The portioning is shown in graph  220 . For example, given a two-dimensional point cloud with σ points (shown as circles) with a 4 by 4 bounding box, and where X and Y coordinates are given in binary code the process of building the tree (shown on the right hand side of  FIG. 2 ) by dividing the bounding box in half along the X-axis (where x 0  represents the most significant bit of the X coordinate). Each half is then divided along the Y-axis where y 0  represents the most significant bit of the Y coordinate. Note that branch x 0 =0, y 0 =1 is not created because the bounding box for this branch is empty. The operations are then repeated for the second-significant bit in X and Y coordinate 
       FIG. 3  shows an illustrative compacted prefix tree  300  in accordance with the embodiment. As shown, compacted prefix tree  300  comprises root node  305  and child nodes  310 - 1  through  310 - 11 . Compacted prefix tree  300  has many useful properties, for example, these tree types have fewer nodes and are a fully binary tree which are required by the operations herein. 
     In accordance with an embodiment, compacted prefix tree  300  is employed to store large amounts of data (e.g., large numbers of nodes) and is traversed in order to retrieve and process assigned node attribute data, and execute any number of different queries, all completely out-of-core Illustratively, well-known techniques such as spatial indexing, spatial queries, viewing frustum, and nearest neighbor searching in point clouds are applied to compacted prefix tree  300 . For example, view frustum culling is the process of removing objects that lie completely outside the viewing frustum from the rendering process. Rending such objects from the collected data set would waste precious processing cycles given that these objects are not directly visible. Typically, to speed up the culling, this is usually done using bounding volumes surrounding the objects rather than the objects themselves, as will now be discussed. 
       FIG. 4  shows a flowchart of illustrative operations  400  for a spatial query (i.e., enumerate all points within given query region of space) using the prefix tree of  FIG. 3  in accordance with an embodiment. In particular, every query starts, at step  410 , as a query for the root node. At step  420 , a prefix tree query is executed (and the axis-aligned node bounding box is calculated) in order to determine, at step  425 , whether there is an intersection of the query region with the axis-aligned bounding box (AABB) of the root node (e.g., root node  305 ). If there is no intersection, the process is aborted. If there is an intersection then, at step  430 , a determination is made if the query region completely contains the AABB. If so, the process short-circuits to the leaf descendants (i.e., directly enumerate all leaves under this node) at step  435 . Otherwise, if the query region partially contains the AABB, as determined in step  430 , there is a recursive processing of the node&#39;s children, if any, at steps  440  and  445  for the left and right children. Recursive processing of the node&#39;s children means that the process node query steps (i.e., steps  420 - 445 ) are executed for both the left child and right child of the current node. Of course, processing child nodes can invoke processing of grand-children nodes and so on. Eventually the process reaches a point where leaf nodes do not have children, and the bounding box is a point for them, such that partial intersection is impossible for them, so such node is either in or out. 
     For example,  FIG. 5  shows illustrative results  500  obtained in performing the operations set forth in  FIG. 4 . As shown, query region  505  results in root node  510  having children nodes  515 - 1  through  515 - 7 . Query region  505  partially intersects the bounding box of the full data set (i.e., AABB of root node), so the children are processed. Further, query region  505  does not intersect the AABB of the left child (i.e., x 0 =0), so this node can be skipped, and query region  505  partially intersects the AABB of the right child (x 0 =1), so these children nodes are processes. As shown, left child (i.e., x 0 =1, y 0 =0, x 1 =0, y 1 =0) is a leaf node which is within the query region, so this node is enumerated as a result of the query. Query region  505  does not intersect the AABB of the right child (x 0 =1, y 0 =1), so this node is skipped. Therefore, this illustrative query returns just one node (i.e., x=10, y=00). 
     As detailed above, the various embodiments herein can be embodied in the form of methods and apparatuses for practicing those methods. The disclosed methods may be performed by a combination of hardware, software, firmware, middleware, and computer-readable medium (collectively “computer”) installed in and/or communicatively connected to a user device.  FIG. 6  is a high-level block diagram of an exemplary computer  600  that may be used for implementing point cloud rendering for real-time point cloud data collection from a variety of sensor types in accordance with the various embodiments herein. Computer  600  comprises a processor  610  operatively coupled to a data storage device  620  and a memory  630 . Processor  610  controls the overall operation of computer  600  by executing computer program instructions that define such operations. Communications bus  660  facilitates the coupling and communication between the various components of computer  600 . Of course, computer  600  may be any type of computing device such a computer, tablet, server, mobile device, smart phone, to name just a few. The computer program instructions may be stored in data storage device  620 , or a non-transitory computer readable medium, and loaded into memory  630  when execution of the computer program instructions is desired. 
     Thus, the steps of the disclosed method (see, e.g.,  FIGS. 1   4 , and  7 ) and the associated discussion herein above) can be defined by the computer program instructions stored in memory  630  and/or data storage device  620  and controlled by processor  610  executing the computer program instructions. For example, the computer program instructions can be implemented as computer executable code programmed by one skilled in the art to perform the illustrative operations defined by the disclosed method. Accordingly, by executing the computer program instructions, processor  610  executes an algorithm defined by the disclosed method. Computer  600  also includes one or more communication interfaces  650  for communicating with other devices via a network (e.g., a wireless communications network) or communications protocol (e.g., Bluetooth®). For example, such communication interfaces may be a receiver, transceiver or modem for exchanging wired or wireless communications in any number of well-known fashions. Computer  600  also includes one or more input/output devices  640  that enable user interaction with the user device (e.g., camera, display, keyboard, mouse, speakers, microphone, buttons, etc.). 
     Processor  610  may include both general and special purpose microprocessors, and may be the sole processor or one of multiple processors of computer  600 . Processor  610  may comprise one or more central processing units (CPUs), for example. Processor  610 , data storage device  620 , and/or memory  630  may include, be supplemented by, or incorporated in, one or more application-specific integrated circuits (ASICs) and/or one or more field programmable gate arrays (FPGAs). 
     Data storage device  620  and memory  630  each comprise a tangible non-transitory computer readable storage medium. Data storage device  620 , and memory  630 , may each include high-speed random access memory, such as dynamic random access memory (DRAM), static random access memory (SRAM), double data rate synchronous dynamic random access memory (DDR RAM), or other random access solid state memory devices, and may include non-volatile memory, such as one or more magnetic disk storage devices such as internal hard disks and removable disks, magneto-optical disk storage devices, optical disk storage devices, flash memory devices, semiconductor memory devices, such as erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), compact disc read-only memory (CD-ROM), digital versatile disc read-only memory (DVD-ROM) disks, or other non-volatile solid state storage devices. 
     Input/output devices  640  may include peripherals, such as a camera, printer, scanner, display screen, etc. For example, input/output devices  640  may include a display device such as a cathode ray tube (CRT), plasma or liquid crystal display (LCD) monitor for displaying information to the user, a keyboard, and a pointing device such as a mouse or a trackball by which the user can provide input to computer  600 . 
     In accordance with various embodiments, a point cloud filter method and apparatus for use in point cloud rendering from real-time point cloud data collection (e.g., as detailed herein above) from a variety of sensor types is provided that delivers enhanced performance including reducing processing requirements, limiting local memory consumption, and optimizing overall data visualization. 
     More particularly, in accordance with an embodiment, a point cloud filter is employed to smooth and/or resample collected point cloud data (e.g., as detailed herein above) in order to smooth the associated point cloud thereby increasing precision while preserving the details and features from the original point cloud data set. In accordance with the embodiment, flat surfaces defined in the point cloud data set are detected and adjustments are made in individual point position (e.g., towards the associated surface) in order to reduce the overall noise level in the collected point cloud data. 
     In the event of a particular dense collected point cloud data set, in accordance with an embodiment, an estimation is made with respect to the overall feature amount in the point cloud data set and resampling is applied to those particular areas requiring less feature definition (e.g., planar surface) in order to preserve the original density for areas with increased feature definitions. In the embodiment, at least three parameters are utilized in the removal of redundant data points: (i) Sigma (σ): defined as the size of the Gaussian filter aperture. All feature details smaller than σ will be treated as noise and removed. The larger the σ value the larger the smoothing and, illustratively, optimal σ ranges for the subject point cloud filter are in the range of 0.05-0.2 meters; (ii) feature preservation: for very large objects (i.e., larger than the specific σ value), sharp edges are detected and preserved from smoothing. The overall feature preservation applied is governed by a specified feature parameter in the range of 1 to 1000 (where 1 means the entire data set is smoothed, and 1000 means no smoothing is applied whatsoever and everything in the point cloud data set is a feature). Illustratively, in accordance with an embodiment, optimal feature preservation values are in the range of 10 to 16; and (iii) caching type: this parameter dictates whether filtering caching will be applied or not which translates to the overall processing speed (i.e., either precise or fast) applied for feature detection in the point cloud filter. In accordance with the embodiments herein, two computational speeds are facilitated and defined as “precise” and “fast”, respectively, wherein the precise speed applies certain filtering caching and the fast speed does not apply any such caching. 
       FIG. 7  shows a flowchart of illustrative operations  700  for point cloud filtering in accordance with an embodiment. The operations encompass both of the aforementioned precise and fast embodiments which are differentiated, in part, by the application of the filter cache, as will be detailed below. At step  705 , the point cloud data set is received and, at step  710 , redundant points are removed. As detailed above, at least three parameters (i.e., 6, feature preservation, and caching) are defined in order to facilitate the removal of redundant points, and these parameters may be, illustratively, user-defined or automatically assigned by the filtering system/processor. At step  715 , a determination is made as to whether caching will be applied. In accordance with an embodiment, caching may be applied to accelerate the overall processing speed of the filtering operations (i.e., the fast speed embodiment) or not applied for more precision (i.e., the precision speed embodiment). If caching is to be applied, the data points are stored in a filter cache, at step  720 , and are processed until the full point cloud data set is processed. 
     At step  725 , the filtering operations continue such that for every point P i  in the point cloud, all points P j  that belong to some area (A) around point P i  are processed. In particular, area A is defined as a 4σ*4σ*4σ cube in well-known earth-centered, earth-fixed (ECEF) coordinates. Next, at step  730 , well-known principal component analysis (PCA) is applied to identify a covariance matrix S defined as: 
     
       
         
           
             
               S 
               = 
               
                 
                   
                     [ 
                     
                       
                         
                           
                             S 
                             xx 
                           
                         
                         
                           
                             S 
                             xy 
                           
                         
                         
                           
                             S 
                             xz 
                           
                         
                       
                       
                         
                           
                             S 
                             xy 
                           
                         
                         
                           
                             S 
                             yy 
                           
                         
                         
                           
                             S 
                             yz 
                           
                         
                       
                       
                         
                           
                             S 
                             xz 
                           
                         
                         
                           
                             S 
                             xz 
                           
                         
                         
                           
                             S 
                             xz 
                           
                         
                       
                     
                     ] 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     where 
                     : 
                     
                       
 
                     
                     ⁢ 
                     
                       S 
                       x 
                     
                   
                 
                 = 
                 
                   
                     1 
                     N 
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       P 
                       
                         j 
                         , 
                         x 
                       
                     
                   
                 
               
             
             , 
             
               
                 S 
                 y 
               
               = 
               
                 
                   1 
                   N 
                 
                 ⁢ 
                 
                   
                     ∑ 
                     
                       
                         P 
                         j 
                       
                       ∈ 
                       A 
                     
                   
                   ⁢ 
                   
                     P 
                     
                       j 
                       , 
                       y 
                     
                   
                 
               
             
             , 
             
               
                 S 
                 z 
               
               = 
               
                 
                   1 
                   N 
                 
                 ⁢ 
                 
                   
                     ∑ 
                     
                       
                         P 
                         j 
                       
                       ∈ 
                       A 
                     
                   
                   ⁢ 
                   
                     P 
                     
                       j 
                       , 
                       z 
                     
                   
                 
               
             
           
         
       
       
         
           
             
               
                 S 
                 xx 
               
               = 
               
                 
                   - 
                   
                     S 
                     x 
                     2 
                   
                 
                 + 
                 
                   
                     1 
                     N 
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       
                         P 
                         
                           j 
                           , 
                           x 
                         
                       
                       2 
                     
                   
                 
               
             
             , 
             
               
                 S 
                 yy 
               
               = 
               
                 
                   - 
                   
                     S 
                     y 
                     2 
                   
                 
                 + 
                 
                   
                     1 
                     N 
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       
                         P 
                         
                           j 
                           , 
                           y 
                         
                       
                       2 
                     
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 S 
                 zz 
               
               = 
               
                 
                   - 
                   
                     S 
                     y 
                     2 
                   
                 
                 + 
                 
                   
                     1 
                     N 
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       
                         P 
                         
                           j 
                           , 
                           z 
                         
                       
                       2 
                     
                   
                 
               
             
           
         
       
       
         
           
             
               
                 
                   S 
                   xy 
                 
                 = 
                 
                   
                     
                       - 
                       
                         S 
                         x 
                       
                     
                     ⁢ 
                     
                       S 
                       y 
                     
                   
                   + 
                   
                     
                       1 
                       N 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           
                             P 
                             j 
                           
                           ∈ 
                           A 
                         
                       
                       ⁢ 
                       
                         
                           P 
                           
                             j 
                             , 
                             x 
                           
                         
                         ⁢ 
                         
                           P 
                           
                             j 
                             , 
                             y 
                           
                         
                       
                     
                   
                 
               
               , 
               
                 
                   S 
                   xz 
                 
                 = 
                 
                   
                     
                       - 
                       
                         S 
                         y 
                       
                     
                     ⁢ 
                     
                       S 
                       z 
                     
                   
                   + 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       
                         P 
                         
                           j 
                           , 
                           x 
                         
                       
                       ⁢ 
                       
                         P 
                         
                           j 
                           , 
                           z 
                         
                       
                     
                   
                 
               
               , 
               
                 
 
               
               ⁢ 
               
                 
                   S 
                   yz 
                 
                 = 
                 
                   
                     
                       - 
                       
                         S 
                         y 
                       
                     
                     ⁢ 
                     
                       S 
                       z 
                     
                   
                   + 
                   
                     
                       ∑ 
                       
                         
                           P 
                           j 
                         
                         ∈ 
                         A 
                       
                     
                     ⁢ 
                     
                       
                         P 
                         
                           j 
                           , 
                           y 
                         
                       
                       ⁢ 
                       
                         P 
                         
                           j 
                           , 
                           z 
                         
                       
                     
                   
                 
               
             
             ⁢ 
             
                 
             
           
         
       
       
         
           
             
               and 
               ⁢ 
               
                   
               
               ⁢ 
               
                 P 
                 
                   j 
                   , 
                   x 
                 
               
             
             , 
             
               P 
               
                 j 
                 , 
                 x 
               
             
             , 
             
               P 
               
                 j 
                 , 
                 z 
               
             
             , 
             
               are 
               ⁢ 
               
                   
               
               ⁢ 
               x 
             
             , 
             y 
             , 
             
               z 
               ⁢ 
               
                   
               
               ⁢ 
               coordinates 
               ⁢ 
               
                   
               
               ⁢ 
               of 
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   P 
                   j 
                 
                 . 
               
             
           
         
       
     
     At step  735 , using the covariance matrix S, eigenvalues λ 1 , λ 2 , λ 3  (λ 1 &gt;λ 2 &gt;λ 3 ) and corresponding eigenvectors v 1 , v 2 , v 3  are determined. For example,  FIG. 8  shows an illustrative point cloud  800  filtered in accordance with an embodiment. As shown, point cloud  800  is filtered, as detailed above, such that for every point P i  in the point cloud, all points P j    810 - 1  that belong to area (A)  805  around point P i    810 - 2  are processed with eigenvalues  815 - 1 ,  815 - 2  and  815 - 3 , and eigenvectors  820 - 1 ,  820 - 2  and  820 - 3 , respectively. In this way,  FIG. 8  illustrates a two-dimensional projection whereby v 2  is perpendicular to the other respective eigenvectors (i.e., v 1 , and v 3 ) At step  740 , Gaussian filtration (a well-known technique) for point P i    810 - 1  is performed in order to find position correction vector F i  as follows: 
     
       
         
           
             
               
                 F 
                 i 
               
               = 
               
                 
                   
                     ∑ 
                     
                       
                         P 
                         j 
                       
                       ∈ 
                       A 
                     
                   
                   ⁢ 
                   
                     
                       P 
                       i 
                     
                     ⁢ 
                     
                       P 
                       j 
                     
                     ⁢ 
                     
                       w 
                       ij 
                     
                   
                 
                 
                   
                     ∑ 
                     j 
                   
                   ⁢ 
                   
                     w 
                     ij 
                   
                 
               
             
             , 
             
               
                 where 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   w 
                   ij 
                 
               
               = 
               
                 exp 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         P 
                         i 
                       
                       ⁢ 
                       
                         
                           P 
                           j 
                         
                         · 
                         
                           P 
                           i 
                         
                       
                       ⁢ 
                       
                         P 
                         j 
                       
                     
                     
                       2 
                       ⁢ 
                       
                         σ 
                         2 
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
     At step  745 , filtering is performed only in a direction perpendicular to the surface so that in-plane point position are left intact. So a new correction G i  is found as follows:
 
 G   i =( F   i   ·v   3 ) v   3  
 
At step  750 , a determination is made as to whether the points in point cloud  800  in area A  805  represent a surface, and if so, at step  755 , a surface detection function f(λ 1 , λ 2 , λ 3 ) is applied such that: 0≤f(λ 1 , λ 2 , λ 3 )≤1, f(λ 1 , λ 2 , λ 3 )=1 when points in area A  805  represent a surface, and f(λ 1 , λ 2 , λ 3 )=0 when points in area A  805  do not represent surface. Illustratively, the surface detection function is given by:
 
 f (λ 1 ,λ 2 ,λ 3 )=1−(αλ 3 /λ 2 ) 2  capped to [0,1]
 
     At step  760 , a correction factor to P i  is applied to the extent allowed by f(λ 1 , λ 2 , λ 3 ) to determine a corrected point P′ i  as follows:
 
 P′   i   =P   i   +f (λ 1 ,λ 2 ,λ 3 )* G   i  
 
At step  765 , the filtered point cloud data set results are output.
 
     In accordance with the embodiment set forth in  FIG. 7 , for every point P i  in the point cloud this point&#39;s neighbors P j  within 4σ*4σ*4σ cube are found. Points are processed through points P j  (illustratively, the operations in  FIG. 7  are conducted on a single pass through the collected data points of the point cloud) and a calculation of S and F i  is made. The overall complexity of these operations is given by the formula: O(N*M avg ), where N is the total number of points in the point cloud, and M avg  is the average number of points per query. Assuming the point cloud comprises mainly surfaces (e.g., terrain, buildings, etc.) with density given by: Density [points/m 2 ], then the overall complexity is O(N*Density*σ 2 ). As such, processing may become slower with the square of the sigma parameter setting (i.e., for deep filtration). It will also be noted that this is just the complexity of filtration itself, and the complexity of the point cloud database random access is different (O(N*log N)). As will be appreciated, the complexity of any filtration operation that uses a binary search tree (BST) cannot be less than the full enumeration of the BST which is (O(N*log N)), where N is the total number of points in the point cloud, and log N is the average complexity of BST access. In this way, (O(N*log N) represents the theoretical maximum for filtration speed. 
     As such, the complexity associated with the “precise” speed embodiment herein is equal to the complexity where calculating S and F i  is equal to O(N*M avg ), and M avg  is the average number of points per query (at steps  825  and  835 , as detailed above). For example, if the point cloud mainly comprises surfaces (e.g., terrain, buildings, etc.) with a “Density” given by points/m 2  and a total area A in m 2 , then the overall complexity is O(N*Density*σ 2 )=O(N 2 *σ 2 /A). As such the complexity may become O(N 2 ) if the density is increasing at a constant area and this is significantly slower that the theoretical maximum O(N)*log N). 
     As noted above, processing speed may be varied, in accordance with an embodiment, by employing filter caching (i.e., fast speed) such that the filtering operations build several grids on top of the filter cache every time the cache is reloads. In this way, all the point calculations detailed above will use only these several grids such the overall complexity of the filtering operations does not depend upon the total number of points in the filter cache and the complexity approaches O(N). To reduce the filtration operational complexity and increase filtration speed for dense data sets the “fast” speed embodiment is used which processes points from the data set in batches. Each batch takes points that fall into the box and builds three (3) different grids on top of box (i.e., initial box plus 2a from all sides). Then all per point calculations use only these grids so that overall complexity does not depend on the number of points in filter cache and thus is close to O(N). 
     In the fast speed embodiment the 3 grids are defined as: big, small, and Gauss. Big and small grids are used to store principal component analysis, and the Gauss grid for Gauss filter calculations. Big and small grids have 9*9*9 cells, each cell represents a box of 4σ*4σ*4σ (i.e., big grid) or 2σ*2σ*2σ (i.e., small grid) centered around (x 0 +i*σ,y 0 +j*σ,y 0 +k*σ), where (x 0 ,y 0 ,y 0 ) is a center of the filter cache, and i, j, k are integers ranging from −4 to +4. The boxes for grid cells are overlapping and when filter cache is loaded principal component analysis (i.e., calculate the S matrix, vectors and lambdas) is performed and with the storage of eigenvalues and eigenvectors separately for each cell(s) in these grids, with an illustrative example shown in  FIG. 9 . As shown, in  FIG. 9  for point cloud  900  the three principal component grids have overlapping boxes (i.e., box  910 ,  920 , and  930 ) with eigenvalues (i.e., eigenvalue  940  and  950 ) and eigenvectors (i.e., eigenvectors  960  and  970 ) stored from the principal component analysis. 
     The Gauss grid has 24*24*24 cells, ½σ*½σ*½σ each, no overlaps, and covers the whole 12σ*12σ*12σ of the filter cache.  FIG. 10  shows an illustrative Gauss grid  1040  for a point cloud in accordance with an embodiment. Each cell (i.e., cell  1005 ,  1010 ,  1015 ,  1020 ) stores a number m of points of the point cloud (i.e., point  1030 - 1 ,  1030 - 2 ,  1030 - 3 ,  1030 -N) that fall in this cell and their center mass (i.e., center mass  1050 - 1 ,  1050 - 2 ,  1050 - 3 , and  1050 - 4 ), as shown in  FIG. 10 . Calculating the correction vector F i , precedes as detailed above (see steps  735 - 740  in  FIG. 7 ) but in this embodiment the centers of mass  1015 - 1 ,  1052 - 2 ,  1050 - 3 , and  1050 - 4  from Gauss grid  10140  are utilized instead of actual points. As such, the Gaussian filter calculation always takes a number of steps defined by 7*7*7 regardless of the density and Sigma parameters. For example, the correction vector F i  calculation is given by: 
     
       
         
           
             
               
                 F 
                 i 
               
               = 
               
                 
                   
                     ∑ 
                     
                       
                         P 
                         j 
                       
                       ∈ 
                       A 
                     
                   
                   ⁢ 
                   
                     
                       P 
                       i 
                     
                     ⁢ 
                     
                       P 
                       j 
                     
                     ⁢ 
                     
                       w 
                       ij 
                     
                   
                 
                 
                   
                     ∑ 
                     j 
                   
                   ⁢ 
                   
                     w 
                     ij 
                   
                 
               
             
             , 
             
               
                 where 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   w 
                   ij 
                 
               
               = 
               
                 m 
                 · 
                 
                   exp 
                   ⁡ 
                   
                     ( 
                     
                       
                         
                           P 
                           i 
                         
                         ⁢ 
                         
                           
                             P 
                             j 
                           
                           · 
                           
                             P 
                             i 
                           
                         
                         ⁢ 
                         
                           P 
                           j 
                         
                       
                       
                         2 
                         ⁢ 
                         
                           σ 
                           2 
                         
                       
                     
                     ) 
                   
                 
               
             
           
         
       
     
     In instances where determining a surface detection function f(λ 1 , λ 2 , λ 3 ) for filtration is desired, the smallest of the aforementioned stored values for the small and big grid are used. In special cases when too few points are available in the grid cell for reliable analysis the following is applied: if the small grid cell has less than twelve (12) points then f(λ 1 , λ 2 , λ 3 ) from the big grid is used; or if the big grid cell has less than twenty (20) points then f(λ 1 , λ 2 , λ 3 )=0, and no filtration occurs. As will be appreciated, the small grids help to detect smaller objects (e.g., objects that could be lost in the big grid if the rest of the surface is perfectly flat), and the big grids help with a big object with many points that can disturb a flat surface even from further away such that the small grid does not see the object at all. In accordance with the fast speed embodiment, the points in the input filter cache are enumerated one time to build the respective grids. A such, all grids are built for the full point cloud by enumerating around (12*12*12/(8*8*8))*N=3.375*N. Further, filtration also takes a fixed number of steps per point, and as mentioned above, O(N) defines the complexity of the filtration itself and the complexity of point cloud database random access is at least O(N*log N). 
     It should be noted that for clarity of explanation, the illustrative embodiments described herein may be presented as comprising individual functional blocks or combinations of functional blocks. The functions these blocks represent may be provided through the use of either dedicated or shared hardware, including, but not limited to, hardware capable of executing software. Illustrative embodiments may comprise digital signal processor (“DSP”) hardware and/or software performing the operation described herein. Thus, for example, it will be appreciated by those skilled in the art that the block diagrams herein represent conceptual views of illustrative functions, operations and/or circuitry of the principles described in the various embodiments herein. Similarly, it will be appreciated that any flowcharts, flow diagrams, state transition diagrams, pseudo code, program code and the like represent various processes which may be substantially represented in computer readable medium and so executed by a computer, machine or processor, whether or not such computer, machine or processor is explicitly shown. One skilled in the art will recognize that an implementation of an actual computer or computer system may have other structures and may contain other components as well, and that a high level representation of some of the components of such a computer is for illustrative purposes. 
     The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.