1 Introduction
Light Detection and Ranging Technology (LiDAR) has been used during the last four decades in Geosciences as a geomatic method to obtain the 3D geometry of the surface of objects [7]. LiDAR uses a laser beam to compute the distance between a device and an object that reflects the beam. When the laser is used to scan the entire field of view at a high speed, the result is a dense cloud of points centered at the device, with each point having additional information such as the intensity of the laser beam reflection. If the device has additional sensors, each point in the cloud will have additional attributes (e.g., if the device includes a camera, each point will have a color value associated).
The decrease in cost of laser scanning devices has helped to drastically increase the application fields for point clouds. For instance, laser scanning has been used to classify and recognize objects in urban environments
[23], in natural environments (e.g., landslides [11], or forests [9]), or even underwater environments [18]. Hence, huge datasets are being produced that require immense computing resources to be processed. To give two examples, the ISPRS benchmark on indoor modelling consists of five point clouds containing points [12] and the mobile laser scanning test data “MLS 1  TUM City Campus” contains more than points collected in 15 minutes.Due to the size of the obtained data, the use of compression is almost mandatory. However, the traditional approach of keeping the data compressed in disk and decompressing the whole dataset before any processing is not effective. Therefore, the classical method for storing LiDAR data (LAZ format) is able to compress the data but, in addition, permits accessing a datum or portions of the data without the need of decompressing the whole dataset. In addition, it is equipped with an index to accelerate the queries. Observe that this setup is very similar to that of many modern compact data structures [17] .
The use of compression for spatial data is not exclusive of LAZ format. In the case of raster data, GeoTiff and NetCDF are able to store the data in compressed form, and in the case of NetCDF, it also permits querying the data directly in that format. Therefore, the application of the knowledge acquired in compact data structures soon led to a new research line.
The wavelet tree [8] was the first compact data structure used in the scope of spatial data. The work in [2] proposes a new point access method based on it. In [16], it is used to represent a set of points in the twodimensional space, each one with an associated value given by an integer function. Another family of compact data structures based on the quadtree also arose. In [5], a compact version of the region quadtree was adapted for storing and indexing compressed rasters. In the same line, Ladra et. al [14] presented an improvement on the previous works.
In this work we continue in this line, now tackling even more complex spatial data. We present a compact data structure to represent LiDAR point clouds, denoted , which compresses and indexes the data. The improvements with respect to the LAZ format are in two aspects. The LAZ format relies on differential encoding plus an entropy encoder, which compresses/decompresses data by blocks, and thus, in order to retrieve a small region, one or more complete blocks must be decompressed. Another drawback of the LAZ format is that it uses a quadtree to index the points, and this only accelerates the queries by the and coordinates. Our new lidar is able to retrieve/decompress a given datum and indexes the three dimensions of the space.
2 Related work: LAS and LAZ format
The American Society for Photogrammetry and Remote Sensing^{1}^{1}1https://www.asprs.org/ defined in 2003 the LASer (LAS) file format, an open data exchange format for LiDAR point data records [21]. The format contains binary data consisting of a header block that describes general information of the point cloud (e.g., number of points, bounding box), variablelength records to describe additional information such as georeferencing information or other metatada, and point records. Each point record consists of a collection of fields describing the point (e.g., the point , and coordinates represented as 4byte integers, the return intensity represented as a 2byte integer, or the laser pulse return number). The standard defines a Point Data Record Format 0, and allows for additional formats to be defined with additional data (e.g., the Point Data Record Format 1 adds the GPS time at which the point was acquired as an 8byte double).
The version 1.4 R14 of the LAS file format specification has been released in 2019 [22]. It defines 11 point data record formats (some of them to provide legacy support), a mechanism to customize the LAS file format to meet applicationspecific needs by adding point classes and attributes, extended variablelength records to carry larger payloads, and additional types of variablelength records (e.g., georeferencing using Well Known Text descriptions, textual description of the LAS file content, or extra bytes for each point record).
Even though the LAS file format avoids using unnecessary space by storing the coordinates as scaled and offset integers, a LAS file requires much storage space (e.g., a 13.2M point cloud requires 254 MB of disk space, see Table 1). The LAZ file format [10], defined by the LASzip lossless compressor for LiDAR, achieves high compression rates supporting streaming, and random access decompression. LASzip encodes the points in the cloud using chunks of 50,000 points. For each chunk, the first point is stored as raw bytes and it is used as the initial value for subsequent prediction schemes. Each additional point is compressed using an entropy coder (an adaptive, contextbased arithmetic coding [20]). These techniques make the LAZ file format very efficient in terms of space. Isenburg [10] showed that the compression ratio is similar or better than generalpurpose compression formats such as ZIP or RAR, while maintaining the possibility of processing the file as a stream of points or directly accessing a particular point without having to decompress the complete file.
The LAZ file format is also very efficient answering range queries on the and dimensions because LAZ files can be indexed using an adaptative quadtree over these coordinates. Each quadtree leaf contains a list of point indexes that can be used to determine the chunk that contains the point. To resolve a range query, the quadtree is first traversed to determine the candidate point indexes, then, the relevant chunks are retrieved, and finally the chunks have to be decompressed and sequentially scanned to determine the points in the result.
Considering additional types of queries, the LAZ file format is highly inefficient on threedimensional queries or queries over attribute data because a sequential scan has to be performed over the points in a chunk. These queries are becoming more common because classification algorithms over point clouds quite often require to locate close points in the threedimensional space, or close points in a twodimensional space with a similar attribute value (e.g., having the same intensity). Both types of queries can be efficiently answered if a threedimensional index is built over the data.
3 Background: trees and trees
The tree [4] is a time and spaceefficient version of a region quadtree [13, 1]. Considering a binary matrix of size , it is divided into quadboxes (submatrices) of size . Each quadbox produces a child of the root node. The label of the node is 1bit if the corresponding quadbox contains at least one 1bit, and 0bit otherwise. The quadboxes having at least one 1bit are divided using the same procedure until reaching a quadbox full of 0bits, or reaching the cells of the original matrix.
The tree, instead of using a classical pointerbased representation of the tree, represents the quadtree using only sequences of bits. More concretely, it uses two bitmaps, denoted as and , where is formed by a breadthfirst traversal of the internal nodes, whereas the is formed by the leaves of the tree. From this basic version, several other improvements yield better space and time performance [4]. Among them, one is that instead of diving each quadbox into quadboxes of size , each division produces quadboxes of size , where is a parameter that can be adapted for each level of the tree. This is usually used to obtain shorter and wider trees that, at the price of a slightly worse space consumption, are faster when querying.
As in the case of the quadtree, where the simple addition of a third dimension produces the octree [15], the tree [5] is simply a 3dimensional tree. Figure 1 shows a 3dimensional binary matrix, its corresponding octree, and the tree represented using bitmaps and .
The tree can be efficiently navigated using rank and select operations^{2}^{2}2Given a bitmap , is the number of occurrences of bit in and is the th occurrence of bit in . over and (see [6, Section 6.2.1]).
4 Our proposal: lidar
In this section, we present the lidar, a data structure to represent LiDAR point clouds in compact space, which allows us to perform efficient queries.
4.1 Conceptual description
Consider a 3D matrix of size that stores a set of LiDAR points. In addition to its coordinates, a point contains several values that correspond to each of its attributes (e.g., intensity, scan angle, color, etc.).
Observe that the tree was designed to store and index a matrix with a bit value for all positions, that is, a 3D binary raster. However, when dealing with LiDAR datasets, we have to change that raster approach to a pointbased data structure. In this scenario, many positions of a LiDAR matrix can be empty, that is, there are no points in those positions.
The lidar recursively divides the matrix into several equalsized submatrices, by following the same strategy used by the tree. Again, this recursive subdivision is represented as a tree, where each node corresponds to a submatrix. As in the original tree, the division of the submatrices stops when the submatrix is empty (full of 0bits, in the case of the tree) and when the process reaches the cells of the original matrix (submatrices of sixe ), placed at level . In addition to these cases, the subdivision of the lidar also stops when the number of points in the processed submatrix is less than or equal to a given threshold .
Regarding the attributes of each point, their values are compactly stored in our structure and can be efficiently retrieved when obtaining a point. Our structure stores the attributes defined in Point Data Record Format 0 (LAS Specification 1.4 [22]) but it can be easily adapted to allow attributes defined in other formats.
4.2 Data structures
We use several data structures to represent the conceptual tree previously described:

Tree structures ( and ): We use two data structures to represent the topology of the tree. is a bitmap, similar to that of the tree, but a 0 means that the submatrix is empty or the number of points does not exceed the threshold . The bitmap of tree is not used. This is due to: i) in many cases, the lidar does not reach the last level of subdivision (level ), as the division frequently stops before that level due to empty submatrices or submatrices containing points or less; and ii) in case of reaching the last level of the subdivision, leaf points require a more complex data structure than just a bitmap.
In addition to T, we use another bitmap, , which has one bit for each 0bit in , plus, if the last level is reached, as many bits as cells in that last level. Each bit differentiates empty submatrices with those that contain some points. That is, for each such that , we set iff the submatrix is empty or in other case. In the case of positions corresponding to the last level , a 1bit means that the corresponding cell has points.

Number of points (): We use a bitmap to store the number of points that each leaf contains. When , we store in the number of points, in unary.

Arrays of coordinates (, , and ): Coordinates of a point are stored in arrays X, Y, Z respectively. Points in level are not stored, since their position can be calculated during the descent through the tree. These values are encoded as local coordinates with respect to the node to which they belong. This results in smaller values than the original ones, and thus they are represented with DACs [3]. DACs provide efficient random access to any position and a good compression ratio with small values.

Attributes: Attributes, as intensity, return number, etc., are stored in separated arrays. They follow the same order as the coordinate arrays. The sequence of values are encoded with DACs or with a bitmap when the attribute can be represented with just one bit per point.

Scale factor and Offset: We convert real coordinates into positive integer values. Therefore, we use one scale factor and one different offset for each dimension. The scale factor is a float value that allows us to transform float values into integer values. The offset is an integer value that allows us to translate points to a coordinate system that starts at . In addition to converting negative numbers into positives, the offset also allows us to obtain smaller numbers in the values of the coordinates. For example, if the minimum value in is , we can move all points 1000 positions to the left, that is, the coordinate would be the , the would be the , and so on. These parameters use the same strategy than LAS/LAZ format.
4.3 Construction
Figure 2 shows the recursive division of a cloud of LiDAR points (left), and the lidar representation (righttop). This examples uses and the maximum number of points in a leaf () is 3. Circles represent LiDAR points and they are identified with a unique identifier. In this case, each point only contains an intensity value labeled , where is the identifier of that point. The algorithm starts by dividing the cube in submatrices of equal size. We add a child node of the root for each submatrix. The first child (bottomleft) has 4 points (5, 8, 7 and 6) and, since , we set . This node is then enqueued to be processed later. The second node is empty, thus and . Note that the third submatrix (bottomleft) has 3 points (2, 3 and 1), thus, since , we set , , and add 3 in unary (001) to . Moreover, local coordinates are stored into arrays , , and following a order. Point 2 having global coordinates becomes . Hence, , , and . This point has an intensity value of , so we set . We repeat the same process with nodes 3 and 1, in that order. The algorithm continues with the rest of submatrices until reaching leaf nodes. In case of having more attributes, we would create a sequence for each attribute in an analogous way to the sequence of intensities .
At the bottomright of Figure 2 we include the final structures that represent the lidar of the example.
4.4 Query
In this section, we describe two queries designed for the lidar, which are of interest for LiDAR point clouds.
4.4.1 Obtaining all points of a region (getRegion):
The lidar is able to obtain all points of a given region by performing a topdown traversal of the tree from the root node. We follow the branches corresponding to submatrices that overlap with the region of interest. When a leaf is found, the coordinates of its points are checked and the attributes are only retrieved if the point is within the region. Due to the fact that points follow a order, some mechanisms are included to decrease the number of points checked. For instance, when we reach a point and , the process stops as we can assure that there are no more valid points in that node.
Algorithm 1 shows a pseudocode of the algorithm to solve this query. Let , , , , , , and be global parameters and the size of the dataset. Parameter stores the total number of 1s in bitmap and was calculated previously. The parameter is the list of points returned by the query. Given a region , the first call of the algorithm is getRegion.
Lines 1–3 run through each child node that overlaps the defined region. (Line 4) is the position in of the current child node. The condition in Line 5 determines if the node is at level (which corresponds to submatrices of size ) or not.
Line 6 checks if the node is an internal node or a leaf node. In the first case, a recursive call is invoke (Lines 7–10). Function getLocalCoordinates converts the given region into local coordinates with respect to the current node. The position of its children is calculated as .
Lines 11–22 are executed when the algorithm reaches a leaf node. Line 12 counts the number of 0bits in until position , i.e the number of leaves until that position (). The condition in Line 13 checks if the node is not empty. Line 14 counts the number of leaf nodes containing points () until position .
Lines 15–17 obtain the positions of the first point and the last point of the current child node in arrays , , and . Since the number of points has been inserted in unary code, each node corresponds to a 1bit in the bitmap. With the select operation, we obtain the position of the 1bit corresponding to the previous node and the 1bit of the current node. Intermediate positions are points of the current node. Finally, in Lines 18–22, the algorithm gets the coordinates of each point and checks if it belongs to the region of interest. If affirmative, the corresponding attributes are added and the point is inserted into the final result.
When the algorithm reaches a node in the level , Lines 24–32 are executed. Line 24 calculates the position in , recall that level is not represented in
. Lines 26–29 are equal to lines 14–17. Finally, for each point in the node, the algorithm retrieves its information and adds the point to the list. Observe that it is not necessary to obtain the local coordinates of the vectors
, and .4.4.2 Obtaining all points of a region filtered by attribute value (filterAttRegion):
Given region and a range of values for an attribute , this query obtains all points within the defined region with values for the attribute between and . Again, this query performs a topdown traversal of the tree. Unlike the query getRegion, when the algorithm reaches a leaf node, in addition to the coordinates, it also retrieves the attribute value of the point.
5 Experimental evaluation
We ran some experiments as a proof of concept of the good properties of our proposed data format. More concretely, we compare the space and time results obtained by lidar, to those obtained when using LAS/LAZ formats.
Table 1 shows the description of the LiDAR point clouds used in the experimental evaluation. We use five different datasets coming from two different sources, an airborne LiDAR and a mobile laser scanning:

Three datasets were created from the union of different files of the Plan Nacional de Observación del Territorio^{3}^{3}3http://pnoa.ign.es/productos_lidar (PNOA). Each tile (file) represents an area of Spanish territory of size 22 km with a minimum density of 0.5 points/m. PNOAsmall is composed of 4 tiles and represents an area of 16 km, PNOAmedium is composed of 9 tiles and represents an area of 36 km, and PNOAlarge is composed of 23 tiles and represents an area of 92 km. The number of points can vary from one tile to another. These datasets contain the following attributes: intensity (with values between 0 and 255), return number (with values between 1 and 4), number of returns (with values between 1 and 4), edge of flight line (with values between 0 and 1), scan direction flag (with values between 0 and 1), classification (with values between 1 and 7), scan angle rank (with values between 24 and 28), and point source ID (with values between 175 and 227).

We use datasets TUB1 and FireBrigade from the ISPRS benchmark on indoor modelling^{4}^{4}4http://www2.isprs.org/commissions/comm4/wg5/benchmarkonindoormodelling.html [12]. TUB1 point cloud was captured in one of the buildings of the Technische Universität Braunschweig, Germany, and FireBrigade was captured in the office of fire brigade in Delft, The Netherlands. These datasets do not contain values for any additional attribute, just the point coordinates. The level of clutter, defined as the amount of points belonging to elements that do not constitute the building structures, is low for TUB1 and high for FireBrigade.
In all cases, LiDAR points were converted to Point Data Record Format 0 (LAS Specification 1.4). Then we created indexes for the LAZ files using the lasindex tool of LAStools^{5}^{5}5https://github.com/LAStools/LAStools. These indexes are able to index the and dimensions to improve the query time. LASLib library^{6}^{6}6https://github.com/LAStools/LAStools/tree/master/LASlib was used to execute queries on LAZ files. The lidar was configured with and .
All the experiments were run on an isolated Intel^{®} Core^{TM} i73820 CPU @ 3.60 GHz (4 cores) with 10 MB of cache, and 64 GB of RAM. It ran Debian 9.8 Stretch, using gcc version 6.3.0 with 03 option.
PNOAsmall  PNOAmedium  PNOAlarge  TUB1  FireBridge  
#points  13,265,144  25,108,130  52,627,503  32,597,694  10,406,389 
min x (real)  546,000.00  544,000.00  542,000.00  9.90  44.32 
max x (real)  549,999.99  549,999.99  551,999.99  5.76  58.40 
min y (real)  4,798,000.00  4,798,000.00  4,794,000.00  24.52  23.10 
max y (real)  4,801,999.99  4,803,999.99  4,805,020.45  18.30  77.67 
min z (real)  39.14  162.43  162.43  1.46  1.27 
max z (real)  179.58  1005.05  1005.05  1.10  12.04 
max X (converted)  3,999,990  5,999,990  9,999,990  1,000,000,000  1,000,000,000 
max Y (converted)  3,999,990  5,999,990  11,020,450  1,000,000,000  1,000,000,000 
max Z (converted)  218,720  1,167,480  1,167,480  1,000,000,000  1,000,000,000 
The comparison of space is shown in the first columns of Table 2. LAZ files obtain the best results in the three cases, around 65% less than lidar, which in turn, needs around 53% less space than the uncompressed LAS.
Table 2 also shows the query times. We generated 500 random regions of different sizes. However, the LAZ software failed in many queries (the reported result contains 0 points). Therefore we only considered the times of those techniques that worked properly, that is, LAZ and lidar. Our proposal is around 5 times faster than LAZ in GetRegion queries for PNOA datasets and from 16–23 times faster for ISPRS datasets. FilterAttRegion queries were only executed over PNOA datasets filtering by the intensity attribute, as TUB1 and FireBrigade only contain the coordinates of the points, but do not include any other attribute values. For PNOA datasets, our proposal outperforms LAZ format, as queries are solved 5–10 times faster.
Dataset  # points  Space (MB)  GetRegion (ms)  FilterAttRegion (ms)  

LAS  LAZ  lidar  LAZ  lidar  LAZ  lidar  
PNOAsmall  13,265,144  254  43  119  1,524  249  1,517  145 
PNOAmedium  25,108,130  479  80  225  2,521  424  2,655  374 
PNOAlarge  52,627,503  1004  173  471  6,859  1,189  6,283  1,264 
TUB1  32,597,694  622  196  304  6,145  383  –  – 
FireBrigade  10,406,389  199  77  100  1,717  74  –  – 
6 Conclusions
In this work, we address the main drawback of the LAZ format for LiDAR data, which is its high executing times when answering to queries that retrieve a subset of points using constraints over the third dimension. LAZ is penalized by the fact that decompression is performed by blocks and the index only covers the X and Y coordinates.
We propose a new representation for LiDAR point clouds, denoted lidar, which is able to decompress random points of the cloud and, as it is based on a compact version of an octree, the tree, it can index the three dimensions. This implies significant improvements in the querying times, ranging from 5 times to more than one order of magnitude faster.
The future work will cover two main lines. First, regarding the space requirements, lidar compresses 65% less than LAZ. Therefore, the use of previous ideas from the original tree and tree, such as using different values in different levels, or further compressing frequent submatrices, will be studied. The other line is the indexation of more dimensions, to include the attribute values stored at the points. The raster is a compact data structure designed for raster data that not only indexes the data spatially, but it also indexes the values stored at the cells of the raster. Our aim is to apply similar ideas to address the indexation of LiDAR data.
References
 [1] Cited by: §3.
 [2] (2009) A new point access method based on wavelet trees. In Proc. 3rd International Workshop on Semantic and Conceptual Issues in GIS (SeCoGIS), LNCS 5833, pp. 297–306. Cited by: §1.
 [3] (2013) DACs: bringing direct access to variablelength codes. Information Processing and Management, pp. 392–404. Cited by: 3rd item.
 [4] (2014) Compact representation of web graphs with extended functionality. Information Systems, pp. 152–174. Cited by: §3, §3.
 [5] (2013) Compact querieable representations of raster data. In Proc. 20th SPIRE, pp. 96–108. Cited by: §1, §3.
 [6] (2014) New data structures and algorithms for the efficient management of large spatial datasets. Ph.D. Thesis, Universidade da Coruña. Cited by: §3.
 [7] (2017) LiDAR remote sensing and applications. CRC Press. Cited by: §1.
 [8] (20031214 January) Highorder entropycompressed text indexes. In Proceedings of SODA 2003, Baltimore, Maryland, pp. 841–850. External Links: ISBN 0898715385 Cited by: §1.
 [9] (2018) Forest inventory using laser scanning. In Topographic Laser Ranging and Scanning, pp. 379–412. Cited by: §1.
 [10] (2013) LASzip: lossless compression of lidar data. Photogrammetric Engineering and Remote Sensing 79 (2), pp. 209–217. Cited by: §2.
 [11] (2012) Use of lidar in landslide investigations: a review. Natural hazards 61 (1), pp. 5–28. Cited by: §1.
 [12] (2017) THE isprs benchmark on indoor modelling.. International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences 42. Cited by: §1, 2nd item.
 [13] (1971) Pattern and search statistics. Academic Press. External Links: Document, ISBN 9780126045505 Cited by: §3.
 [14] (2017) Scalable and queryable compressed storage structure for raster data. Information Systems 72, pp. 179–204. Cited by: §1.
 [15] (1982) Geometric modeling using octree encoding. Computer Graphics and Image Processing 19 (2), pp. 129 – 147. External Links: ISSN 0146664X, Document, Link Cited by: §3.
 [16] (2013) Spaceefficient dataanalysis queries on grids. Theoretical Computer Science 482, pp. 60–72. Cited by: §1.
 [17] (2016) Compact data structures: a practical approach. Cambridge University Press.. Cited by: §1.
 [18] (2018) 3D laser scanner for underwater manipulation. Sensors 18 (4), pp. 1086. Cited by: §1.
 [19] (201104) Eyedome lighting: a nonphotorealistic shading technique. Technical report Cited by: Figure 4.
 [20] (2002) Arithmetic coding. In Lossless compression handbook, Cited by: §2.
 [21] (2003) ASPRS LIDAR Data Exchange Format Standard. Version 1.0. Note: Format Specification Cited by: §2.
 [22] (2019) LAS Specification 1.4  R14. Note: Format Specification Cited by: §2, §4.1.
 [23] (2018) LiDAR point clouds to 3d urban models a review. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 11 (2), pp. 606–627. Cited by: §1.
Appendix
To better understand the nature of the datasets, we show a visualization of PNOAlarge in Figure 3, and visualizations of the point clouds TUB1 and FireBrigade in Figure 4.
Comments
There are no comments yet.