From 0eb30408453121c0627c533c9c75aa561e3c9652 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Mon, 1 Nov 2021 15:47:00 -0400 Subject: [PATCH 1/4] Modify offset calculation to take scale factor into account. --- epf/Epf.cpp | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/epf/Epf.cpp b/epf/Epf.cpp index 654f8a3..6ba4342 100644 --- a/epf/Epf.cpp +++ b/epf/Epf.cpp @@ -203,10 +203,6 @@ void Epf::fillMetadata(const pdal::PointLayoutPtr layout) m_b.pointSize += pdal::Dimension::size(di.type); m_b.dimInfo.push_back(di); } - m_b.offset[0] = m_b.bounds.maxx / 2 + m_b.bounds.minx / 2; - m_b.offset[1] = m_b.bounds.maxy / 2 + m_b.bounds.miny / 2; - m_b.offset[2] = m_b.bounds.maxz / 2 + m_b.bounds.minz / 2; - auto calcScale = [](double scale, double low, double high) { if (scale > 0) @@ -221,10 +217,23 @@ void Epf::fillMetadata(const pdal::PointLayoutPtr layout) return std::pow(10, (std::max)(power, -4.0)); }; - m_b.scale[0] = calcScale(m_b.scale[0], m_b.bounds.minx, m_b.bounds.maxx); - m_b.scale[1] = calcScale(m_b.scale[1], m_b.bounds.miny, m_b.bounds.maxy); - m_b.scale[2] = calcScale(m_b.scale[2], m_b.bounds.minz, m_b.bounds.maxz); + m_b.scale[0] = calcScale(m_b.scale[0], m_b.trueBounds.minx, m_b.trueBounds.maxx); + m_b.scale[1] = calcScale(m_b.scale[1], m_b.trueBounds.miny, m_b.trueBounds.maxy); + m_b.scale[2] = calcScale(m_b.scale[2], m_b.trueBounds.minz, m_b.trueBounds.maxz); + + // Find an offset such that (offset - min) / scale is close to an integer. + auto calcOffset = [](double minval, double maxval, double scale) + { + double interval = maxval - minval; + double spacings = interval / scale; // Number of quantized values in our range. + double halfspacings = spacings / 2; // Half of that number. + double offset = (int32_t)halfspacings * scale; // Round to an even value and scale down. + return minval + offset; // Add the base (min) value. + }; + m_b.offset[0] = calcOffset(m_b.trueBounds.minx, m_b.trueBounds.maxx, m_b.scale[0]); + m_b.offset[1] = calcOffset(m_b.trueBounds.miny, m_b.trueBounds.maxy, m_b.scale[1]); + m_b.offset[2] = calcOffset(m_b.trueBounds.minz, m_b.trueBounds.maxz, m_b.scale[2]); } PointCount Epf::createFileInfo(const StringList& input, StringList dimNames, From 2e6eb181a761360f48e7731a54efbdb5f5beed19 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Mon, 1 Nov 2021 15:53:09 -0400 Subject: [PATCH 2/4] Update comment. --- epf/Epf.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/epf/Epf.cpp b/epf/Epf.cpp index 6ba4342..4f7f34e 100644 --- a/epf/Epf.cpp +++ b/epf/Epf.cpp @@ -221,7 +221,12 @@ void Epf::fillMetadata(const pdal::PointLayoutPtr layout) m_b.scale[1] = calcScale(m_b.scale[1], m_b.trueBounds.miny, m_b.trueBounds.maxy); m_b.scale[2] = calcScale(m_b.scale[2], m_b.trueBounds.minz, m_b.trueBounds.maxz); - // Find an offset such that (offset - min) / scale is close to an integer. + // Find an offset such that (offset - min) / scale is close to an integer. This helps + // to eliminate warning messages in lasinfo that complain because of being unable + // to write nominal double values precisely using a 32-bit integer. + // The hope is also that raw input values are written as the same raw values + // on output. This may not be possible if the input files have different scaling or + // incompatible offsets. auto calcOffset = [](double minval, double maxval, double scale) { double interval = maxval - minval; From 0c8b849becfd3ed14d613d523d98e4ca19a06a75 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Tue, 2 Nov 2021 15:16:26 -0400 Subject: [PATCH 3/4] Fix comment. --- epf/Epf.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/epf/Epf.cpp b/epf/Epf.cpp index 4f7f34e..47fc811 100644 --- a/epf/Epf.cpp +++ b/epf/Epf.cpp @@ -232,7 +232,7 @@ void Epf::fillMetadata(const pdal::PointLayoutPtr layout) double interval = maxval - minval; double spacings = interval / scale; // Number of quantized values in our range. double halfspacings = spacings / 2; // Half of that number. - double offset = (int32_t)halfspacings * scale; // Round to an even value and scale down. + double offset = (int32_t)halfspacings * scale; // Round to an int value and scale down. return minval + offset; // Add the base (min) value. }; From dbe9bd0ba2dee4517300155577289aabe933c3b2 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Mon, 24 Jan 2022 13:28:18 -0500 Subject: [PATCH 4/4] Use pre-existing offset if possible. --- epf/Epf.cpp | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/epf/Epf.cpp b/epf/Epf.cpp index c7319dd..5265f5c 100644 --- a/epf/Epf.cpp +++ b/epf/Epf.cpp @@ -285,6 +285,10 @@ PointCount Epf::createFileInfo(const StringList& input, StringList dimNames, filenames.push_back(filename); } + std::vector xOffsets; + std::vector yOffsets; + std::vector zOffsets; + // Determine a driver for each file and get a preview of the file. If we couldn't // Create a FileInfo object containing the file bounds, dimensions, filename and // associated driver. Expand our grid by the bounds and file point count. @@ -315,6 +319,15 @@ PointCount Epf::createFileInfo(const StringList& input, StringList dimNames, m = root.findChild("scale_z"); if (m.valid()) m_b.scale[2] = (std::max)(m_b.scale[2], m.value()); + m = root.findChild("offset_x"); + if (m.valid()) + xOffsets.push_back(m.value()); + m = root.findChild("offset_y"); + if (m.valid()) + yOffsets.push_back(m.value()); + m = root.findChild("offset_z"); + if (m.valid()) + zOffsets.push_back(m.value()); FileInfo fi; fi.bounds = qi.m_bounds; @@ -340,9 +353,26 @@ PointCount Epf::createFileInfo(const StringList& input, StringList dimNames, totalPoints += fi.numPoints; } + // If we had an offset from the input, choose one in the middle of the list of offsets. + if (xOffsets.size()) + { + std::sort(xOffsets.begin(), xOffsets.end()); + m_b.offset[0] = xOffsets[xOffsets.size() / 2]; + } + if (yOffsets.size()) + { + std::sort(yOffsets.begin(), yOffsets.end()); + m_b.offset[1] = yOffsets[yOffsets.size() / 2]; + } + if (zOffsets.size()) + { + std::sort(zOffsets.begin(), zOffsets.end()); + m_b.offset[2] = zOffsets[zOffsets.size() / 2]; + } + // If we have LAS start capability, break apart file infos into chunks of size 5 million. #ifdef PDAL_LAS_START - PointCount ChunkSize = 5000000; + PointCount ChunkSize = 5'000'000; for (const FileInfo& fi : tempFileInfos) { if (fi.driver != "readers.las" || fi.numPoints < ChunkSize)