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 1_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.
CodePudding user response：
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.