Reading .las version 1.3 using C++, and displaying the results using pcl library

958 Views Asked by At

I'm fairly new to the Point Cloud world. In addition, I'm not so experienced in C++. I need to read .las files and process them using the pcl library. This is a sample file from the dataset that I need to read. I followed this youtube video. However, since the file I'm trying to read is of version 1.3, I followed the corresponding spec file to define the header fields and I used 'Data Record Format 3' which is the data record format mentioned in the file header.

This are my definitions for the header and data record format:

#pragma once
#include <string>
#include <vector>

struct float4
{
    float x, y, z, intensity;
};

class PointCloud
{
    public:

        uint32_t getVertsCount();
        float4* getVertsData();

        template<typename PointT>
        typename pcl::PointCloud<PointT>::Ptr read(const std::string& path);//void read(const std::string &path);

    private:
        std::vector<float4> verts;
        #pragma pack(1)
        struct Header
        {
            char magic[4];
            uint16_t fileSourceID;
            uint16_t globalEncoding;
            uint32_t guidData1;
            uint16_t guidData2;
            uint16_t guidData3;
            uint8_t guidData4[8];
            uint8_t versionMaj, versionMin;
            char systemIdentifier[32];
            char genSoftware[32];
            uint16_t creationDay, creationYear;
            uint16_t headerSize;
            uint32_t pointDataOffset;
            uint32_t numVarLenRecords;
            uint8_t pointDataRecordFormat;
            uint16_t pointDataRecordLen;
            uint32_t numberOfPoints;
            uint32_t numPointsByReturn[5];
            double scaleX, scaleY, scaleZ;
            double offsetX, offsetY, offsetZ;
            double maxX, minX, maxY, minY, maxZ, minZ;
            uint64_t waveform;
        };
        
        //#pragma pack(1)
        struct PointRecord3
        {
            uint32_t x, y, z;
            uint16_t intensity;
            uint8_t flags;
            uint8_t classification;
            uint8_t scanAngleRank;
            uint8_t userData;
            uint16_t pointSourceId;
            double gpsTime;
            uint16_t red;
            uint16_t green;
            uint16_t blue;
        };
};

I used the following code to read the point data, but I failed to get correct points:

template<typename PointT>
typename pcl::PointCloud<PointT>::Ptr PointCloud::read(const string& path)
{
    ifstream inf(path, ios::binary);
    typename pcl::PointCloud<PointT>::Ptr lasCloud(new pcl::PointCloud<PointT>);

    if (inf.is_open())
    {
        Header header;
        inf.read((char*)&header, sizeof(header));
        
        cout << "Signature: " << header.magic << endl;
        cout << "Source ID: " << int(header.fileSourceID) << endl;
        cout << "Global Encoding: " << int(header.globalEncoding) << endl;
        cout << "Guid 1: " << int(header.guidData1) << endl;
        cout << "Guid 2: " << int(header.guidData2) << endl;
        cout << "Guid 3: " << int(header.guidData3) << endl;
        cout << "Guid 4: " << header.guidData4 << endl;
        cout << (int)header.versionMaj << '.' << (int)header.versionMin << endl; 
        cout << "Sys Identifier: " << header.systemIdentifier << endl;
        cout << "Gen Software: " << header.genSoftware << endl;
        cout << "Creation Day: " << header.creationDay << endl;
        cout << "Creation Year: " << header.creationYear << endl;
        cout << header.headerSize << " == " << sizeof(header) << endl;
        cout << "Point Data Offset: " << header.pointDataOffset << endl;
        cout << "Number of Variable Len Records: " << header.numVarLenRecords << endl;
        cout << "point Data Record Format: " << header.pointDataRecordFormat << endl;
        cout << "point Data Record Len: " << header.pointDataRecordLen << endl;
        cout << "Number of Points: " << header.numberOfPoints << endl;
        cout << "Number of Points by Return: " << header.numPointsByReturn << endl;
        cout << "Scales: " << header.scaleX << ", " << header.scaleY << ", " << header.scaleZ << endl;
        cout << "Offsets: " << header.offsetX << ", " << header.offsetY << ", " << header.offsetZ << endl;
        cout << "Xmin = " << header.minX << ", Ymin = " << header.minY << ", Zmin = " << header.minZ << endl;
        cout << "Xmax = " << header.maxX << ", Ymax = " << header.maxY << ", Zmax = " << header.maxZ << endl;
        cout << "Waveform: "<<header.waveform << endl;

        assert(header.versionMaj == 1 && header.versionMin == 3);
        //assert(header.headerSize == sizeof(header));
        assert(header.pointDataRecordFormat == 3);

        //inf.seekg(header.pointDataOffset);
        inf.seekg(sizeof(header));
        //inf.seekg(header.pointDataOffset+sizeof(header.waveform));
        for (uint32_t i = 0; i < header.numberOfPoints; i++)
        {
            //PointRecord1* points = new PointRecord1[header.numberOfPoints];
            PointRecord3 point;
            //inf.read((char*)(points + i), sizeof(PointRecord1));
            //inf.read((char*)&point, sizeof(PointRecord1));
            inf.read((char*)&point, sizeof(PointRecord3));

            PointT cloudPoint;
            cloudPoint.x = (float)(point.x * header.scaleX) + header.offsetX;
            cloudPoint.y = (float)(point.y * header.scaleY) + header.offsetY;
            cloudPoint.z = (float)(point.z * header.scaleZ) + header.offsetZ;
            cloudPoint.intensity = (float)(point.intensity) / 65536.0;
            lasCloud->points.push_back(cloudPoint);
        }
        
        if (!inf.good())
            throw runtime_error("Reading went wrong!");

    }
    else
    {
        throw runtime_error("Can't find any!");
    }

    lasCloud->width = lasCloud->points.size();
    lasCloud->height = 1;
    lasCloud->is_dense = true;
    std::cout << "Cloud size = " << lasCloud->points.size() << endl;
    return lasCloud;
}

int main (int argc, char** argv)
{
    
    std::cout << "starting enviroment" << std::endl;
    pcl::visualization::PCLVisualizer::Ptr viewer (new pcl::visualization::PCLVisualizer ("3D Viewer"));
    CameraAngle setAngle = FPS; //XY, FPS, Side, TopDown
    initCamera(setAngle, viewer);
    
    pcl::PointCloud<pcl::PointXYZI>::Ptr inputCloudI; //
    PointCloud  pcd;
    inputCloudI=pcd.read<pcl::PointXYZI>("C:/Users/hedey/OneDrive/Documents/Research_papers/STDF/10_4231_MFQF-Q141/I-65/LiDAR/RoadSurface/NB/20180524_I65_NB_RoadSurface_1_50.5.las");
    std::cout << "Cloud size = " << inputCloudI->points.size() << endl;

    renderPointCloud(viewer, inputCloudI, "lasCloud");

    while (!viewer->wasStopped())
    {
        viewer->spinOnce();
    }

}

There is a problem that I noticed. The header size is defined to be 227 (this is the value of the header size field). However, there is an 8-byte field named 'Start of Waveform Data Packet Record' at the end of the header, which if included in the header definition will make the header size 235 bytes. Also, the pointDataOffset field which was used to seek the points data in the youtube video is pointing to 227 bytes. When I used it to seek the points data, I got unreasonable point values. My target is to process this point cloud and display it using the pcl cloud, but I'm failing to read the points correctly.

1

There are 1 best solutions below

8
On BEST ANSWER

I have no idea about Lidar and *.las files, but I've noticed that the input file was created with libLAS

> file 20180524_I65_NB_RoadSurface_1_50.5.las 
20180524_I65_NB_RoadSurface_1_50.5.las: LIDAR point data records, version 1.3, SYSID libLAS, Generating Software libLAS 1.6.0

So, why don't you use libLAS to read the data in? https://liblas.org/

libLAS comes with convenient CLI utility programs to handle las files, e.g.:

lasinfo 20180524_I65_NB_RoadSurface_1_50.5.las 
---------------------------------------------------------
  Header Summary
---------------------------------------------------------

  Version:                     1.3
  Source ID:                   0
  Reserved:                    0
  Project ID/GUID:             '00000000-0000-0000-0000-000000000000'
  System ID:                   'libLAS'
  Generating Software:         'libLAS 1.6.0'
  File Creation Day/Year:      144/2018
  Header Byte Size             227
  Data Offset:                 227
  Header Padding:              0
  Number Var. Length Records:  None
  Point Data Format:           3
  Number of Point Records:     22017565
  Compressed:                  False
  Number of Points by Return:  0 0 0 0 0 
  Scale Factor X Y Z:          0.00100000000000 0.00100000000000 0.00100000000000
  Offset X Y Z:                590284.000 4339456.000 157.000
  Min X Y Z:                   589879.772 4338728.975 149.667
  Max X Y Z:                   590334.248 4339568.021 178.397
  Spatial Reference:           None

---------------------------------------------------------
  Schema Summary
---------------------------------------------------------
  Point Format ID:             3
  Number of dimensions:        16
  Custom schema?:              false
  Size in bytes:               34

  Dimensions
---------------------------------------------------------
  'X'                            --  size: 32 offset: 0
  'Y'                            --  size: 32 offset: 4
  'Z'                            --  size: 32 offset: 8
  'Intensity'                    --  size: 16 offset: 12
  'Return Number'                --  size: 3 offset: 14
  'Number of Returns'            --  size: 3 offset: 14
  'Scan Direction'               --  size: 1 offset: 14
  'Flightline Edge'              --  size: 1 offset: 14
  'Classification'               --  size: 8 offset: 15
  'Scan Angle Rank'              --  size: 8 offset: 16
  'User Data'                    --  size: 8 offset: 17
  'Point Source ID'              --  size: 16 offset: 18
  'Time'                         --  size: 64 offset: 20
  'Red'                          --  size: 16 offset: 28
  'Green'                        --  size: 16 offset: 30
  'Blue'                         --  size: 16 offset: 32
  
---------------------------------------------------------
  Point Inspection Summary
---------------------------------------------------------
  Header Point Count: 22017565
  Actual Point Count: 22017565

  Minimum and Maximum Attributes (min,max)
---------------------------------------------------------
  Min X, Y, Z:      600191.027, 4313816.564, 148.621
  Max X, Y, Z:      600212.594, 4314678.007, 156.632
  Bounding Box:     600191.027, 4313816.564, 600212.594, 4314678.007
  Time:         55449082.421688, 55488872.904376
  Return Number:    0, 0
  Return Count:     0, 0
  Flightline Edge:  0, 0
  Intensity:        0, 255
  Scan Direction Flag:  0, 0
  Scan Angle Rank:  0, 0
  Classification:   1, 3
  Point Source Id:  0, 31
  User Data:        0, 0
  Minimum Color (RGB):  0 0 0 
  Maximum Color (RGB):  0 0 0 

  Number of Points by Return
---------------------------------------------------------
    (1) 22017565

  Number of Returns by Pulse
---------------------------------------------------------
    (0) 22017565

  Point Classifications
---------------------------------------------------------
    7187055 Unclassified (1) 
    8128678 Ground (2) 
    6701832 Low Vegetation (3) 
  -------------------------------------------------------
    0 withheld
    0 keypoint
    0 synthetic
  -------------------------------------------------------

and

las2txt 20180524_I65_NB_RoadSurface_1_50.5.las  qq.txt && head qq.txt
600209.243,4313837.086,155.155
600209.342,4313839.620,155.191
600209.232,4313836.806,155.154
600209.338,4313839.516,155.197
600209.221,4313836.523,155.165
600209.333,4313839.398,155.194
600209.206,4313836.177,155.158
600209.328,4313839.285,155.200
600209.189,4313835.778,155.145
600209.322,4313839.152,155.193

This means that the file is OK, the library works, and you'll save yourself lots of time learning the basic library usage rather then trying to reimplement it yourself (think of various data formats etc., error handling etc., testing etc., getting feedback in case of troubles, etc., and your time)