writers.copc: round X/Y/Z values to match LAS

The COPC writer truncated scaled & offset values to int32 during packing. This didn't match the LAS/LAZ writer, which meant that point coordinates could needlessly change during conversion to COPC. Fix this by applying the same rounding approach that is used for LAS.
koordinates · Sep 20, 2023 · c418ae2 · c418ae2
1 parent 810b876
commit c418ae2
Show file tree

Hide file tree

Showing 2 changed files with 74 additions and 4 deletions.
diff --git a/io/private/las/Utils.cpp b/io/private/las/Utils.cpp
@@ -682,11 +682,27 @@ void V14BaseLoader::load(PointRef& point, const char *buf, int bufsize)
 void V14BaseLoader::pack(const PointRef& point, char *buf, int bufsize)
 {
     LeInserter ostream(buf, bufsize);
-    int32_t xi = (int32_t)m_scaling.m_xXform.toScaled(point.getFieldAs<double>(Dimension::Id::X));
-    int32_t yi = (int32_t)m_scaling.m_yXform.toScaled(point.getFieldAs<double>(Dimension::Id::Y));
-    int32_t zi = (int32_t)m_scaling.m_zXform.toScaled(point.getFieldAs<double>(Dimension::Id::Z));
 
-    ostream << xi << yi << zi;
+    auto converter = [](double val, Dimension::Id dim) -> int32_t
+    {
+        int32_t i(0);
+
+        if (!Utils::numericCast(val, i))
+            throw std::runtime_error("Unable to convert scaled value (" +
+                Utils::toString(val) + ") to "
+                "int32 for dimension '" + Dimension::name(dim) );
+        return i;
+    };
+
+    double xOrig = point.getFieldAs<double>(Dimension::Id::X);
+    double yOrig = point.getFieldAs<double>(Dimension::Id::Y);
+    double zOrig = point.getFieldAs<double>(Dimension::Id::Z);
+
+    int32_t x = converter(m_scaling.m_xXform.toScaled(xOrig), Dimension::Id::X);
+    int32_t y = converter(m_scaling.m_yXform.toScaled(yOrig), Dimension::Id::Y);
+    int32_t z = converter(m_scaling.m_zXform.toScaled(zOrig), Dimension::Id::Z);
+
+    ostream << x << y << z;
 
     uint16_t intensity = point.getFieldAs<uint16_t>(Dimension::Id::Intensity);
     int returnNum = point.getFieldAs<int>(Dimension::Id::ReturnNumber);

diff --git a/test/unit/io/CopcWriterTest.cpp b/test/unit/io/CopcWriterTest.cpp
@@ -36,6 +36,8 @@
 
 #include <pdal/pdal_test_main.hpp>
 
+#include <pdal/util/FileUtils.hpp>
+#include <io/BufferReader.hpp>
 #include <io/CopcReader.hpp>
 #include <io/CopcWriter.hpp>
 #include <io/LasReader.hpp>
@@ -195,4 +197,56 @@ TEST(CopcWriterTest, srsUTM)
     EXPECT_TRUE(Utils::startsWith(data, "PROJCS[\"NAD83 / UTM zone 15N\""));
 }
 
+TEST(CopcWriterTest, scaling)
+{
+    using namespace Dimension;
+
+    const std::string FILENAME(Support::temppath("copc_scaling.las"));
+    PointTable table;
+
+    table.layout()->registerDims({Id::X, Id::Y, Id::Z});
+
+    BufferReader bufferReader;
+
+    PointViewPtr view(new PointView(table));
+    view->setField(Id::X, 0, 1406018.497);
+    view->setField(Id::Y, 0, 4917487.174);
+    view->setField(Id::Z, 0, 62.276);
+    bufferReader.addView(view);
+
+    Options writerOps;
+    writerOps.add("filename", FILENAME);
+    writerOps.add("offset_x", "1000000");
+    writerOps.add("scale_x", "0.001");
+    writerOps.add("offset_y", "5000000");
+    writerOps.add("scale_y", "0.001");
+    writerOps.add("offset_z", "0");
+    writerOps.add("scale_z", "0.001");
+
+    CopcWriter writer;
+    writer.setOptions(writerOps);
+    writer.setInput(bufferReader);
+
+    writer.prepare(table);
+    writer.execute(table);
+
+    Options readerOps;
+    readerOps.add("filename", FILENAME);
+
+    PointTable readTable;
+
+    LasReader reader;
+    reader.setOptions(readerOps);
+
+    reader.prepare(readTable);
+    PointViewSet viewSet = reader.execute(readTable);
+    EXPECT_EQ(viewSet.size(), 1u);
+    view = *viewSet.begin();
+    EXPECT_EQ(view->size(), 1u);
+    EXPECT_NEAR(1406018.497, view->getFieldAs<double>(Id::X, 0), .00001);
+    EXPECT_NEAR(4917487.174, view->getFieldAs<double>(Id::Y, 0), .00001);
+    EXPECT_NEAR(62.276, view->getFieldAs<double>(Id::Z, 0), .00001);
+    FileUtils::deleteFile(FILENAME);
+}
+
 } // namespace pdal