LASLib C++; Creating new custom las files from mulitple LAS files; how to get accurate & correct LASpoint x,y and z?

[Joep Lijnen]

Hello,

I am working on some C++ code to generate custom tiles containing points of interest, the short version of my question would be:

"How can I correctly and preferably with minimal loss of accuracy represent the x, y and z coordinates an original LAS point to new  LAS file and header?"

I run into a situation I do not know how to solve; I want existing points from one LAS to be correctly represented for another LAS file.
When I just 'copy' points from one LAS to another with the same header, this is not an issue. (as per the classical LASlib example)
When I create new points, I can use the .init() function for the header I create; however; this erases all prior information in the point.


The principle is as follows:

I have a set of LAS input tiles, and a set of predicates that describe a polygon which are to be populated with points from the LAS input tiles.

However, I notice (screenshot attached) that upon writing my LAS file, the points that came from different source-LAS files, are shifted.
(The screenshot depicts a QGIS screenshot, where my 4 artificially constructed original LAS tiles contain 'source' points, the green polygon named 905 represents the target tile to generate, however, the green, outputted points, are not correctly positioned after writing.)

This is the header with which I try to write my points:

LASheader lasheader;
lasheader.x_scale_factor = 0.00001;
lasheader.y_scale_factor = 0.00001;
lasheader.z_scale_factor = 0.00001;
lasheader.x_offset = predicate->boundingBox.averageX();
lasheader.y_offset = predicate->boundingBox.averageY();
lasheader.z_offset = 0.0; 
lasheader.point_data_format = same_data_format;
lasheader.point_data_record_length = same_data_length;

I choose a nice origin (x,y,z_offset), for the best numerical stability of the points.

The //Todo in the code below is likely where my problem and solution are to be found.

while (lasreader->read_point()) //I read the point from one the source las file
{
  const auto &point = &lasreader->point;
  
  if (predicate->operator()(static_cast<point3>(*point))) //If the point is inside interest polygon.
  {
    //Todo: how do I get the point to contain the correct representation for its X/Y in the new header, which lasWriter uses? - the point comes from a file with a different x/y/z_offset.
    lasWriter->write_point(point);   //Here I write the point.
    lasWriter->update_inventory(point);

    pointNumber++;
  }
}

How can I correctly and preferably with minimal loss of accuracy represent the x, y and z coordinates of the original LAS point to another header?

Possibly I can use the raw integer point information, as George Vosselman has done/attempted in the past?, then change the Point's quantizer to the new header and put it back in? 

Ideally however, I would like to see a simple method to accurately 'reproject' a LASpoint to another 'LASheader', perhaps I have overlooked it.

Thank you for any insight/help.

Kind Regards,

Joep Lijnen

[Evon Silvia]

Basically, I think you need to convert the integer-form XYZ values to double-floats by applying the scale/offset from the reader. Then you'd apply the scale-offset from the writer to convert them to new integers. In my own code I do this myself by following the order of operations in the LAS specification and assigning the new XYZ values to point directly, but you could also use the LASlib functions.

If you're using the LASzip DLL, I'd encourage you to look at EXAMPLE_SIX as it does most of what you're trying to do. Take particular note of the use of laszip_auto_offset on line 1338, and laszip_set_coordinates on line 1405. There was a laszip function to get the transformed XYZ values directly, but now I can't find it. Hopefully someone else can. 

Evon

--

Evon Silvia PLS

QSI Solutions Developer

ASPRS LAS Working Group Chair

Quantum Spatial

517 SW 2nd Street, Suite 400, Corvallis, OR 97333

P: (541) 452-8502

E: esi...@quantumspatial.com

--
Download LAStools at
http://lastools.org
http://rapidlasso.com
Be social with LAStools at
http://facebook.com/LAStools
http://twitter.com/LAStools
http://linkedin.com/groups/LAStools-4408378
Manage your settings at
http://groups.google.com/group/lastools/subscribe

[Martin Isenburg]

Hello,

Here is some sample code (not tested but some variation of this should work) using LASlib (LGPL 2.1). This assumed the points have the same type and size. Also see the las2las.cpp code how the reprojection is implemented. Or the lasreadermerged.cpp code. If the quantizers are drastically different this will not work. Note that your precision chosen below is in increments of 10th of a micrometer (0.00001). An individual human hair has a diameter of 40 to 120 micrometer  (maybe depending on which shampoo you use - at least that's what those commercials on TV say ... (-;). So unless you are scanning with some form of super high-accuracy ranging device I suggest to go up to centimeter (scale factors of 0.01) or millimeters (scale factors of 0.001). For regular airborne scans integer overflows are likely.with micrometer scale factors ...

======================

LASreadOpener lasreadopener;

LASwriteOpener laswriteopener;

LASreader* lasreader1 =  lasreadopener.open("file1.laz);

LASreader* lasreader2 =  lasreadopener.open("file2.laz);

laswriteopener set_file_name("merged.laz"

LASwriter* laswriter = laswriteopener.open(&lasreader1->header);

// just write points of first file

while (lasreader1->read_point())

{

   laswriter->write_point(&lasreader1->point);

   laswriter->update_invenory( &lasreader1->point);

}

lasreader1->close();

// requantize points of second file before writing

while (lasreader2->read_point())

{

   // store the current coordinates computed with the quantizer of lasreader2 as F64 internally in the point

   lasreader2->point.compute_coordinates();

   // recompute the I32 coordinates using the quantizer of lasreader1 from the F64 coordinates computed in the previous step

   lasreader2->point.compute_XYZ(&lasreader1->header);

   laswriter->write_point(&lasreader2->point);

   laswriter->update_inventory(&lasreader2->point);

}

lasreader2->close();

laswriter->update_inventory();

laswriter->close();

delete laswriter;

delete lasreader1;

delete lasreader2;

======================

[Jérôme Chatillon]

Hi Joep,

I'll have to check my code, but i think I did something like that:

// create and initialize the point you want to write
LASpoint writterLASpoint;
writterLASpoint.init( &lasWritterHeader, lasWritterHeader.point_data_format, lasWritterHeader.point_data_record_length);
// copy information from the lasreader point
writterLASpoint = lasreader->point;
// reassign position since your laswritter has different scale and offset [set x/y scale to 0.001 for airborne]
writterLASpoint.set_x( lasreader->point.get_x() );
writterLASpoint.set_y( lasreader->point.get_y() );
writterLASpoint.set_z( lasreader->point.get_z() );

Jerome

[Joep Lijnen]

Ah, thank you everyone for your suggestions, apologies for not responding;
I believed that Google had failed to post my message, I only just found out that it had indeed successfully posted.

Thank you for bringing to my attention the LAS Specification, Evon.
To conclude my story, after the 'failed' post, went on to search a solution, which I did find after analyzing the LASPoint anatomy; I am happy to see it matches the ideas of the ones suggested by Jérome and Martin;

(some code omitted regarding initialization/cleanup):

 while (lasreader->read_point()) //process all points from source tile las/laz.

{
    const auto &point = &lasreader->point;

    
    //The 2 lines below get the RAW integer coordinates, saves these in a temporary variable, then converts
    //back to raw using the SPECIFIED las header. This BREAKS the point afterward, its data is now incompatible
    //with the original header; it is now valid for the updated header.
    point->compute_coordinates();
    point->compute_XYZ(&lasheader);
    
    lasWriter->write_point(point);   //Write the point, valid for the updated header.
    lasWriter->update_inventory(point);
}

@Martin thanks for your detailed reply. Las in question is from a mobile mapping platform. I figured I had an 'overkill' in accuracy. After getting my code to work, I verified the accuracy-loss due to changing the header offsets was very small; which it was; 4 micrometers.
Is there any associated downside to picking a scale factor which is too small? - Will this e.g. make it impossible to represent a point correctly that is too far away from the center? - Currently tiles I generate are less than 100 meters in diameter..
In any case, I suppose I could safely remove one or two zeroes from the scale.

Thanks again for responding so quickly, I'll have to be more careful with Google groups, before giving up on it working :)

Kind Regards,

Joep Lijnen

ITC Developer 

Department of Earth Observation Science (ITC-EOS)

On Tuesday, 10 April 2018 14:07:10 UTC+2, Joep Lijnen wrote: