Odm: PDAL pipeline distorts point cloud in certain datasets

[pierotofy]

Input dataset: https://www.dropbox.com/sh/d28k2p7silx5ows/AAD0jglE7C70DuNFT0zu-xEda?dl=0

Problem happens with either OpenSfM or PMVS.

The resulting odm_georeferenced_model.ply point cloud is generated correctly (see below):

The odm_georeferenced_model.ply.las file generated with /code/SuperBuild/build/pdal/bin/pdal pipeline -i /var/www/data/53683562-45b6-495e-86d6-f00a22208e68/odm_georeferencing/pipeline.xml --readers.ply.filename=/var/www/data/53683562-45b6-495e-86d6-f00a22208e68/odm_georeferencing/odm_georeferenced_model.ply --writers.las.filename=/var/www/data/53683562-45b6-495e-86d6-f00a22208e68/odm_georeferencing/odm_georeferenced_model.ply.las however is "striped". See below:

Potree display:

CloudCompare display:

PLY file: https://ufile.io/08728
LAS file: https://ufile.io/d3d511

[pierotofy]

pipeline.xml

<?xml version="1.0" encoding="UTF-8"?>
<Pipeline version="1.0">
   <Writer type="writers.las">
      <Option name="filename">transformed.las</Option>
      <Option name="a_srs">EPSG:32756</Option>
      <Filter type="filters.transformation">
         <Option name="matrix">1  0  0  340695        0  1  0  6259514        0  0  1  0        0  0  0  1</Option>
         <Reader type="readers.ply">
            <Option name="filename">untransformed.ply</Option>
         </Reader>
      </Filter>
   </Writer>
</Pipeline>

[smathermather]

Wacky. Does it do the same with a JSON file instead of XML?

[pierotofy]

Here's another dataset that has the same problem:

https://drive.google.com/open?id=0B4Tf0aIvxRLJTnpKOWlsU3ZaN2c

[smathermather]

Have we tried processing this with a json file instead of the XML? My guess is this is an XML regression that hasn't been caught in PDAL because we are the only ones left using XML in PDAL... . Sorry I haven't had a change to look into this myself.

[dakotabenjamin]

https://github.com/dakotabenjamin/OpenDroneMap/blob/pdal-dem/opendm/types.py I've changed to JSON in this branch. Can someone test?

[pierotofy]

Just tested, I think we should switch to JSON regardless, but it doesn't fix the problem.

I think the problem is due to the limitations of floating point precision.

The matrix in the first dataset for example reads:

"1 0 0 340695 0 1 0 6259514 0 0 1 0 0 0 0 1

Plugging in some coordinates from input files and multiplying them by the matrix yields several points with more than 7 digits of precision.

[dakotabenjamin]

https://github.com/dakotabenjamin/OpenDroneMap/tree/pdal-dem I just pushed an update that forwards the scale. It looks good to me.

[pierotofy]

I confirm that adding "forward": "scale" to the pipeline fixes the problem.

[hblanken]

Works for me! :1st_place_medal:

[dakotabenjamin]

Great. Look forward to a PR in the next week or so.

[hblanken]

Has this now been been integrated into PR #549. Will test both this weekend.

[dakotabenjamin]

Yes it also fixed that by forwarding the scale through the pipeline: https://github.com/OpenDroneMap/OpenDroneMap/pull/549/files#diff-02fe4f960c0e32f739192d0fdfa38756R185

[dakotabenjamin]

@hblanken did you get good results?

[hblanken]

Sorry to only reply now - I pulled your new ODM master and unfortunately have not been able to create a las without stripes. I am now running ODM natively and produce output via plas.io

[pierotofy]

After pulling ODM master, did you remove all previous results from the dataset folder (all the directories expect the images/ directory)?

[hblanken]

@pierotofy Yes, clean dataset folder. I tried today with a dataset with a gcp_list.txt and have similar results. Both in Meshlab and and in plas.io and webODM i get stripes.

[dakotabenjamin]

Have the changes been reflected in WebODM yet @pierotofy ?

[pierotofy]

It has not; the docker image hasn't been upgraded in a while (the build process was stuck, I need to code up some sort of notification system for when this happens in the future).

I'm rebuilding it now, it should be ready in a few hours.

Keep an eye on the "Last pushed" field here https://hub.docker.com/r/opendronemap/opendronemap/ and here https://hub.docker.com/r/opendronemap/node-opendronemap/

Then ./webodm.sh update should do it.

[pierotofy]

Ok, the new images have been built. @hblanken could you confirm that things now work after an update?

[hblanken]

Hi @pierotofy - Thanks for trying to get to the bottom of this. I pulled the latest docker webodm and node and now the potree visual looks better, with quite large bubbles as default (which I know can be adjusted in the potree options), but the stripes are gone. See below. When I pull all the odm assets from docker and view the .las file in plas.io - it still has the stripes. Question: This means, if I run ODM natively (which I tend to do), without WedODM, the stripes will remain in the .las files? (Note for you - the screenshot in your webODM readme also shows the stripes in the sheffield park potree view. Maybe worth updating to the latest commit)

[pierotofy]

Thanks for pointing out the screenshot in the readme. If you run ODM without WebODM, the stripes should be gone as well. This is a patch we've added in ODM core.

[pierotofy]

Actually, I'm going to backtrack a moment and I think that this is still a problem. The reason the potree point cloud looks fine is because there have been some filename changes (odm_georeferenced_model.ply.las is no longer created, but odm_georeferenced_model.las is), which triggered node-odm's failsafe to generate a las file from the original merged.ply point cloud (which is not georeferenced).

Long story short, this is still a problem, and adding "forward": "scale" does not solve the issue. After reading more, the forward directive only applies when reading from a las/laz source, not ply. So in essence adding "forward" : "scale" doesn't do anything.

The problem is still the matrix transformation:

{
  "pipeline": [
    "untransformed.ply",
    {
      "matrix": "1 0 0 576705 0 1 0 5188170 0 0 1 0 0 0 0 1",
      "type": "filters.transformation"
    },
    {
      "a_srs": "EPSG:32615",
      "filename": "transformed.las"
    }
  ]
}

We lose precision after the transformation is applied. We need to calculate appropriate scale and offset parameters (see https://www.pdal.io/stages/writers.las.html scale_x, scale_y, offset_x, offset_y).

I found this discussion interesting for the purpose: https://gis.stackexchange.com/questions/132632/reprojecting-lidar-data-with-liblas/132635

[pierotofy]

[pierotofy]

@hblanken could you confirm that there's the following line in the last output log for the dataset you processed?

odm_georeferenced_model.ply.las is missing, will attempt to generate it from

@dakotabenjamin perhaps instead of a matrix transformation, we just need to specify offset_x and offset_y?

https://github.com/OpenDroneMap/OpenDroneMap/blob/master/opendm/types.py#L170

[hblanken]

@pierotofy - I believe you mean the geo_referencing logs? I attach both files in the directory:odm_georeferencing_log.txt
odm_georeferencing_utm_log.txt

[pierotofy]

Sorry, no I mean this output log from the WebODM interface:

[hblanken]

Sure - this is what I get, see below: I don't see "odm_georeferenced_model.ply.las is missing"

spacing calculated from diagonal: 0.594157
READING: data/4bc923fd-cbd3-4d72-8331-76411fa99970/odm_georeferencing/odm_georeferenced_model.ply.las
INDEXING: 1,000,000 points processed; 1,000,000 points written; 1.107 seconds passed
closing writer

conversion finished
1,425,761 points were processed and 1,425,761 points ( 100% ) were written to the output.
duration: 2.383s
Compressing all.zip

[dakotabenjamin]

@pierotofy is this an issue with core ODM or node-ODM?

[pierotofy]

ODM. I just processed this with the latest master.

Las file (noticed the stripes):

Ply file (correct point distrubition):

What is the reason for using the matrix transform in the PDAL pipeline?

[dakotabenjamin]

@smathermather would be able to answer that.

[pierotofy]

So I digged through the code, the idea of the matrix transform is to add a offset to the coordinates (x + utm east offset = geo referenced x coordinate, same for y, z is already in meters above the geoid).

Setting an offset_x and offset_y property in the las writer options should have worked, but the problem did not go away.

The reason this is not working is because of the PLY reader!

If you look at the header of our input PLY files:

property float x
property float y
property float z

X,Y,Z coordinates are defined as floats. So PDAL (correctly) creates Point objects that have a float data type: https://github.com/PDAL/PDAL/blob/dd97a033c5d487000e200304a7619e9978de6ae2/filters/TransformationFilter.cpp#L100

    point.setField(Dimension::Id::X,
        x * m_matrix[0] + y * m_matrix[1] + z * m_matrix[2] + m_matrix[3]);

setField will apply the proper offset, but precision is lost in the process because points are treated as floats, not doubles.

Changing the PLY header to:

property double x
property double y
property double z

Yields proper results.

[hblanken]

This error was bugging me for the last weeks with processing several striped Las files. Big breakthrough finding @pierotofy @dakotabenjamin ! So if I understand you correctly the .las is created based on the PLY and if the proper header is passed on for las processing the Las will show up correctly as well? Is this sth I can implement in the code somehow?

[pierotofy]

If the PLY header specifies its dimensions data type as double instead of float, the las file is generated properly. I'm still trying to figure out the best way to fix this, as I'm not sure I'd want to modify the PLY headers for the sake of executing the conversion.

[smathermather]

If PLY is written by OpenSfM, maybe a pull request against that project to correctly specify type in the header is required.

[pierotofy]

This might get fixed in PDAL: https://github.com/PDAL/PDAL/issues/1584

[pierotofy]

It should be fixed in the latest PDAL master: https://github.com/PDAL/PDAL/commit/aea5bb0cacc64b91d626eff491fbdbb5668c06d7

I'll make a PR that uses the latest PDAL version.

[smathermather]

W00t