pythonc++binary-datalas

Opening a LAS file in C++ outputting seemingly incorrect values


I will preface by saying that I am very new to C++, please forgive me if I have made any egregious errors.

I am trying to write a c++ program to input a las file and output each point's coordinate vector. I followed this youtube video (https://www.youtube.com/watch?v=vdaWvEH2xNU) for assistance and I have gotten so close (I think) to getting the correct values.

When I make a python in python to do this, it gives me the correct coordinates in ENU:

import laspy
import numpy as np
las = laspy.read('ppk_cloud_1.las')
point_data = np.stack([las.x, las.y, las.z], axis=0).transpose((1, 0))
print(point_data)

--outputs--

[[2.73720833e+05 4.33655312e+06 1.55323500e+02]
 [2.73720833e+05 4.33655312e+06 1.55323500e+02]
 [2.73720742e+05 4.33655314e+06 1.55200500e+02]
 ...

However my program in c++ outputs these values:

[0]:{x=273720.844 y=8631520.00 z=4295122.50 }
[1]:{x=273720.844 y=8631520.00 z=4295122.50 }
[2]:{x=273720.750 y=8631520.00 z=4295122.50 }
[3]:{x=273720.750 y=8631520.00 z=4295122.50 }
...
[27]:{x=273728.188 y=8631521.00 z=4295116.50 }
[28]:{x=273727.719 y=4336556.00 z=4295117.50 }
[29]:{x=273727.719 y=4336556.00 z=4295117.50 }
[30]:{x=273727.625 y=4336556.00 z=4295117.50 }

The x value is pretty much there (not identical but close), however the other two are not. The y value starts very off and then jumps to close to the correct values (it jumps between these values a few times...really confused by this). The z value is just wrong. (Values should be in meters)

My thoughts as to a possible issue:

(1) Endianess - maybe cpp is reading the binary file in little endian (I am on an Intel computer) and the las file was written in big endian. (ruled this out by realizing that the header file offsets (header.x_off,y_off,z_off ~= 273664.46,4336554.51, 161.24) are pretty much what I want.

(2) It is reading the actual points shifted (not the header) - I double checked the header.point_data_offset value and I am pretty sure that is correct.

(3) In the video he does not use (float)((point.x * header.x_scale) + header.x_off) and instead uses (float)((point.x * header.x_scale) + header.x_off), however in the documentation (https://www.asprs.org/a/society/committees/lidar/Downloads/ASPRS_LAS%201_2.pdf), page 5, it says this "So to scale a given X from the point record, take the point record X multiplied by the X scale factor, and then add the X offset". Thus I think that my formula should be correct (however I believe this to be the source of my issue).

(4) Maybe the values are correct and just in a different coordinate system? (Seems unlikely that the header and points would be in different coordinate systems).

Any help would be GREATLY APPRECIATED. THANK YOU! :)

For reference here is my cpp file:

#include <iostream>
#include <fstream>
#include <stdexcept>
#include <assert.h>
#include "point_cloud.h"


PointCloud::PointCloud(const std::string &path)
{
    read(path);
}

void PointCloud::read(const std::string &path)
{
    std::ifstream inf(path, std::ios::binary);

    if (inf.is_open())
    {
        Header header;
        inf.read((char *)&header, sizeof(header));

        assert(header.vers_major == 1 && header.vers_minor == 2);
        assert(header.header_size == sizeof(header));
        assert(header.point_data_record_format == 1);

        inf.seekg(header.point_data_offset);

        for (uint32_t i = 0; i < header.n_points; i++)
        {
            PointFormat1 point;
            inf.read((char*)&point, sizeof(PointFormat1));

            assert(sizeof(point) == sizeof(PointFormat1));

            xyzfloats r = 
            {
                (float)((point.x * header.x_scale) + header.x_off),
                (float)((point.y * header.y_scale) + header.y_off),
                (float)((point.z * header.z_scale) + header.z_off)
            };

            //std::cout << r.x << "," << r.y << "," << r.z << std::endl;
            
            coordinates.push_back(r);
        }
        std::cout << "done." << std::endl;
    }
    else
    {
        throw std::runtime_error("Can't find file.");
    }
}

uint32_t PointCloud::get_n_points()
{
    return (uint32_t)coordinates.size();
}

xyzfloats *PointCloud::getcoordinates()
{
    return coordinates.data();
}

and here is my header file:

#pragma once

#include <vector>
#include <string>

//https://www.asprs.org/wp-content/uploads/2010/12/LAS_1_4_r13.pdf

struct xyzfloats
{
    float x, y, z;
};


class PointCloud
{
public:

    PointCloud(const std::string &path);

    uint32_t get_n_points();
    xyzfloats *getcoordinates();
    

private:

    std::vector<xyzfloats> coordinates;

#pragma pack(push, 1)
    struct Header
    {
        char     file_sig[4];
        uint16_t file_source_id;
        uint16_t global_encoding;
        uint32_t guid_data_1;
        uint16_t guid_data_2;
        uint16_t guid_data_3;
        uint8_t  guid_data_4[8];
        uint8_t  vers_major, vers_minor;
        char     sys_identifier[32];
        char     generating_software[32];
        uint16_t creation_day, creation_year;
        uint16_t header_size;
        uint32_t point_data_offset;
        uint32_t n_var_length_records;
        uint8_t  point_data_record_format;
        uint16_t point_data_record_length;
        uint32_t n_points;
        uint32_t n_points_by_return[5];
        double   x_scale, y_scale, z_scale;
        double   x_off, y_off, z_off;
        double   x_min, y_min, z_min;
        double   x_max, y_max, z_max;
    };
#pragma pack(pop)

#pragma pack(push, 1)
    struct PointFormat1
    {
        uint32_t x, y, z;
        uint16_t intensity;
        uint8_t  return_info;
        uint8_t  classification;
        uint8_t  scan_angle_rank;
        uint8_t  user_data;
        uint16_t point_src_id;
        double   gps_time;
    };
#pragma pack(pop)

    void read(const std::string &path);

};

PS: I am aware of libLAS, and I am trying to get it working, however I cannot get it installed at the moment.


Solution

  • The LAS format specification uses a signed 32-bit integer as the x/y/z values in the LAS Point Data Format Record 1: http://www.asprs.org/wp-content/uploads/2019/07/LAS_1_4_r15.pdf

    In your PointFormat1 struct, switch to int32_t and it should then work.