[Liblas-commits] hg-main-tree: reproj and scaling filters now working

liblas-commits at liblas.org liblas-commits at liblas.org
Thu Apr 28 12:57:24 EDT 2011


details:   http://hg.libpc.orghg-main-tree/rev/e0413b5a06e6
changeset: 692:e0413b5a06e6
user:      Michael P. Gerlek <mpg at flaxen.com>
date:      Thu Apr 28 09:40:40 2011 -0700
description:
reproj and scaling filters now working
Subject: hg-main-tree: reproj and scaling filters working

details:   http://hg.libpc.orghg-main-tree/rev/23340cd26242
changeset: 693:23340cd26242
user:      Michael P. Gerlek <mpg at flaxen.com>
date:      Thu Apr 28 09:57:02 2011 -0700
description:
reproj and scaling filters working
Subject: hg-main-tree: merge

details:   http://hg.libpc.orghg-main-tree/rev/89f0502b790f
changeset: 694:89f0502b790f
user:      Michael P. Gerlek <mpg at flaxen.com>
date:      Thu Apr 28 09:57:12 2011 -0700
description:
merge

diffstat:

 apps/pc2pc.cpp                          |    5 +-
 include/libpc/filters/ScalingFilter.hpp |    5 +
 src/filters/ReprojectionFilter.cpp      |    6 +-
 src/filters/ScalingFilter.cpp           |   86 +++++++++++---
 test/unit/ReprojectionFilterTest.cpp    |  192 ++++++++++++++++++-------------
 5 files changed, 190 insertions(+), 104 deletions(-)

diffs (truncated from 413 to 300 lines):

diff -r 2c7b904a79a9 -r 89f0502b790f apps/pc2pc.cpp
--- a/apps/pc2pc.cpp	Thu Apr 28 08:25:32 2011 -0700
+++ b/apps/pc2pc.cpp	Thu Apr 28 09:57:12 2011 -0700
@@ -190,12 +190,13 @@
             if (m_srs.size() > 0)
             {
                 ref.setFromUserInput(m_srs);
-                writer.setSpatialReference(ref);                            
+                writer.setSpatialReference(ref);
             }
         }
         if (hasOption("compress"))
         {
-            writer.setCompressed(true);            
+            if (m_bCompress)
+                writer.setCompressed(true);            
         }
         writer.write(numPoints);
 
diff -r 2c7b904a79a9 -r 89f0502b790f include/libpc/filters/ScalingFilter.hpp
--- a/include/libpc/filters/ScalingFilter.hpp	Thu Apr 28 08:25:32 2011 -0700
+++ b/include/libpc/filters/ScalingFilter.hpp	Thu Apr 28 09:57:12 2011 -0700
@@ -59,7 +59,9 @@
     // notes:
     //   - "forward=true" means doubles --> ints
     //   - "forward=false" means ints --> doubles
+    //   - 1st version uses the scale/offset values already present
     ScalingFilter(const Stage& prevStage, bool forward);
+    ScalingFilter(const Stage& prevStage, double scaleX, double offsetX, double scaleY, double offsetY, double scaleZ, double offsetZ, bool forward);
 
     const std::string& getDescription() const;
     const std::string& getName() const;
@@ -80,6 +82,9 @@
     void checkImpedance();
     void initialize();
     
+    bool m_customScaleOffset;
+    double m_scaleX, m_scaleY, m_scaleZ;
+    double m_offsetX, m_offsetY, m_offsetZ;
     bool m_forward;
 
     ScalingFilter& operator=(const ScalingFilter&); // not implemented
diff -r 2c7b904a79a9 -r 89f0502b790f src/filters/ReprojectionFilter.cpp
--- a/src/filters/ReprojectionFilter.cpp	Thu Apr 28 08:25:32 2011 -0700
+++ b/src/filters/ReprojectionFilter.cpp	Thu Apr 28 09:57:12 2011 -0700
@@ -210,9 +210,9 @@
 
     for (boost::uint32_t pointIndex=0; pointIndex<numPoints; pointIndex++)
     {
-        double x = data.getField<boost::int32_t>(pointIndex, indexX);
-        double y = data.getField<boost::int32_t>(pointIndex, indexY);
-        double z = data.getField<boost::int32_t>(pointIndex, indexZ);
+        double x = data.getField<double>(pointIndex, indexX);
+        double y = data.getField<double>(pointIndex, indexY);
+        double z = data.getField<double>(pointIndex, indexZ);
 
         this->transform(x,y,z);
 
diff -r 2c7b904a79a9 -r 89f0502b790f src/filters/ScalingFilter.cpp
--- a/src/filters/ScalingFilter.cpp	Thu Apr 28 08:25:32 2011 -0700
+++ b/src/filters/ScalingFilter.cpp	Thu Apr 28 09:57:12 2011 -0700
@@ -47,6 +47,32 @@
 
 ScalingFilter::ScalingFilter(const Stage& prevStage, bool forward)
     : Filter(prevStage)
+    , m_customScaleOffset(false)
+    , m_scaleX(0.0)
+    , m_scaleY(0.0)
+    , m_scaleZ(0.0)
+    , m_offsetX(0.0)
+    , m_offsetY(0.0)
+    , m_offsetZ(0.0)
+    , m_forward(forward)
+{
+    checkImpedance();
+
+    initialize();
+
+    return;
+}
+
+
+ScalingFilter::ScalingFilter(const Stage& prevStage, double scaleX, double offsetX, double scaleY, double offsetY, double scaleZ, double offsetZ, bool forward)
+    : Filter(prevStage)
+    , m_customScaleOffset(true)
+    , m_scaleX(scaleX)
+    , m_scaleY(scaleY)
+    , m_scaleZ(scaleZ)
+    , m_offsetX(offsetX)
+    , m_offsetY(offsetY)
+    , m_offsetZ(offsetZ)
     , m_forward(forward)
 {
     checkImpedance();
@@ -90,12 +116,24 @@
         schema.removeDimension(dimYd);
         schema.removeDimension(dimZd);
 
-        dimXi.setNumericScale(dimXd.getNumericScale());
-        dimXi.setNumericOffset(dimXd.getNumericOffset());
-        dimYi.setNumericScale(dimYd.getNumericScale());
-        dimYi.setNumericOffset(dimYd.getNumericOffset());
-        dimZi.setNumericScale(dimZd.getNumericScale());
-        dimZi.setNumericOffset(dimZd.getNumericOffset());
+        if (m_customScaleOffset)
+        {
+            dimXi.setNumericScale(m_scaleX);
+            dimXi.setNumericOffset(m_offsetX);
+            dimYi.setNumericScale(m_scaleY);
+            dimYi.setNumericOffset(m_offsetY);
+            dimZi.setNumericScale(m_scaleZ);
+            dimZi.setNumericOffset(m_offsetZ);
+        }
+        else
+        {
+            dimXi.setNumericScale(dimXd.getNumericScale());
+            dimXi.setNumericOffset(dimXd.getNumericOffset());
+            dimYi.setNumericScale(dimYd.getNumericScale());
+            dimYi.setNumericOffset(dimYd.getNumericOffset());
+            dimZi.setNumericScale(dimZd.getNumericScale());
+            dimZi.setNumericOffset(dimZd.getNumericOffset());
+        }
 
         schema.addDimension(dimXi);
         schema.addDimension(dimYi);
@@ -125,12 +163,24 @@
             throw impedance_invalid("Scaling filter requires X,Y,Z dimensions as int32s not be initially present (reverse direction)");
         }
 
-        dimXd.setNumericScale(dimXi.getNumericScale());
-        dimXd.setNumericOffset(dimXi.getNumericOffset());
-        dimYd.setNumericScale(dimYi.getNumericScale());
-        dimYd.setNumericOffset(dimYi.getNumericOffset());
-        dimZd.setNumericScale(dimZi.getNumericScale());
-        dimZd.setNumericOffset(dimZi.getNumericOffset());
+        if (m_customScaleOffset)
+        {
+            dimXd.setNumericScale(m_scaleX);
+            dimXd.setNumericOffset(m_offsetX);
+            dimYd.setNumericScale(m_scaleY);
+            dimYd.setNumericOffset(m_offsetY);
+            dimZd.setNumericScale(m_scaleZ);
+            dimZd.setNumericOffset(m_offsetY);
+        }
+        else
+        {
+            dimXd.setNumericScale(dimXi.getNumericScale());
+            dimXd.setNumericOffset(dimXi.getNumericOffset());
+            dimYd.setNumericScale(dimYi.getNumericScale());
+            dimYd.setNumericOffset(dimYi.getNumericOffset());
+            dimZd.setNumericScale(dimZi.getNumericScale());
+            dimZd.setNumericOffset(dimZi.getNumericOffset());
+        }
 
         schema.removeDimension(dimXi);
         schema.removeDimension(dimYi);
@@ -205,9 +255,9 @@
             const double yd = srcData.getField<double>(pointIndex, indexYd);
             const double zd = srcData.getField<double>(pointIndex, indexZd);
 
-            const boost::int32_t xi = dimXd.removeScaling<boost::int32_t>(xd);
-            const boost::int32_t yi = dimYd.removeScaling<boost::int32_t>(yd);
-            const boost::int32_t zi = dimZd.removeScaling<boost::int32_t>(zd);
+            const boost::int32_t xi = dimXi.removeScaling<boost::int32_t>(xd);
+            const boost::int32_t yi = dimYi.removeScaling<boost::int32_t>(yd);
+            const boost::int32_t zi = dimZi.removeScaling<boost::int32_t>(zd);
 
             dstData.setField<boost::int32_t>(pointIndex, indexXi, xi);
             dstData.setField<boost::int32_t>(pointIndex, indexYi, yi);
@@ -220,9 +270,9 @@
             const boost::int32_t yi = srcData.getField<boost::int32_t>(pointIndex, indexYi);
             const boost::int32_t zi = srcData.getField<boost::int32_t>(pointIndex, indexZi);
 
-            const double xd = dimXi.applyScaling(xi);
-            const double yd = dimYi.applyScaling(yi);
-            const double zd = dimZi.applyScaling(zi);
+            const double xd = dimXd.applyScaling(xi);
+            const double yd = dimYd.applyScaling(yi);
+            const double zd = dimZd.applyScaling(zi);
 
             dstData.setField<double>(pointIndex, indexXd, xd);
             dstData.setField<double>(pointIndex, indexYd, yd);
diff -r 2c7b904a79a9 -r 89f0502b790f test/unit/ReprojectionFilterTest.cpp
--- a/test/unit/ReprojectionFilterTest.cpp	Thu Apr 28 08:25:32 2011 -0700
+++ b/test/unit/ReprojectionFilterTest.cpp	Thu Apr 28 09:57:12 2011 -0700
@@ -60,88 +60,33 @@
 }
 
 
-static void getPoint(const libpc::PointBuffer& data, boost::uint32_t pointIndex, double& x, double& y, double& z)
+static void getPoint(const libpc::PointBuffer& data, double& x, double& y, double& z, double scaleX, double scaleY, double scaleZ)
 {
     using namespace libpc;
 
     const Schema& schema = data.getSchema();
 
-#if 1
     const int indexX = schema.getDimensionIndex(Dimension::Field_X, Dimension::Int32);
     const int indexY = schema.getDimensionIndex(Dimension::Field_Y, Dimension::Int32);
     const int indexZ = schema.getDimensionIndex(Dimension::Field_Z, Dimension::Int32);
-    const Dimension& xDim = schema.getDimension(indexX);
-    const Dimension& yDim = schema.getDimension(indexY);
-    const Dimension& zDim = schema.getDimension(indexZ);
 
-    const boost::int32_t xraw = data.getField<boost::int32_t>(pointIndex, indexX);
-    const boost::int32_t yraw = data.getField<boost::int32_t>(pointIndex, indexY);
-    const boost::int32_t zraw = data.getField<boost::int32_t>(pointIndex, indexZ);
+    const boost::int32_t xraw = data.getField<boost::int32_t>(0, indexX);
+    const boost::int32_t yraw = data.getField<boost::int32_t>(0, indexY);
+    const boost::int32_t zraw = data.getField<boost::int32_t>(0, indexZ);
 
-    x = xDim.applyScaling(xraw);
-    y = yDim.applyScaling(yraw);
-    z = zDim.applyScaling(zraw);
-#else
-    const int indexX = schema.getDimensionIndex(Dimension::Field_X, Dimension::Double);
-    const int indexY = schema.getDimensionIndex(Dimension::Field_Y, Dimension::Double);
-    const int indexZ = schema.getDimensionIndex(Dimension::Field_Z, Dimension::Double);
-    const Dimension& xDim = schema.getDimension(indexX);
-    const Dimension& yDim = schema.getDimension(indexY);
-    const Dimension& zDim = schema.getDimension(indexZ);
-
-    x = data.getField<double>(pointIndex, indexX);
-    y = data.getField<double>(pointIndex, indexY);
-    z = data.getField<double>(pointIndex, indexZ);
-#endif
+    x = (double)xraw * scaleX;
+    y = (double)yraw * scaleY;
+    z = (double)zraw * scaleZ;
 
     return;
 }
 
 
+// Test reprojecting UTM 15 to DD with a filter
 BOOST_AUTO_TEST_CASE(test_1)
 {
-    // Test reprojecting UTM 15 to DD with a filter
-    libpc::drivers::las::LasReader reader1(Support::datapath("utm15.las"));
-    libpc::drivers::las::LasReader reader2(Support::datapath("utm15.las"));
-        
-    const libpc::SpatialReference& in_ref = reader1.getSpatialReference();
-        
     const char* epsg4326_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]";
-    libpc::SpatialReference out_ref(epsg4326_wkt);
-
-    {
-        const char* utm15_wkt = "PROJCS[\"NAD83 / UTM zone 15N\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4269\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-93],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"26915\"]]";
-        BOOST_CHECK(in_ref.getWKT() == utm15_wkt);
-        BOOST_CHECK(out_ref.getWKT() == epsg4326_wkt);
-    }
-
-    libpc::filters::ScalingFilter scalingFilter(reader2, false);
-
-    libpc::filters::ReprojectionFilter reprojectionFilter(scalingFilter, in_ref, out_ref);
-    
-    libpc::filters::ScalingFilter descalingFilter(reprojectionFilter, true);
-
-    const libpc::Schema& schema1 = reader1.getSchema();
-    const libpc::SchemaLayout layout1(schema1);
-    libpc::PointBuffer data1(layout1, 1);
-
-    const libpc::Schema& schema2 = descalingFilter.getSchema();
-    const libpc::SchemaLayout layout2(schema2);
-    libpc::PointBuffer data2(layout2, 1);
-
-    {
-        libpc::SequentialIterator* iter1 = reader1.createSequentialIterator();
-        boost::uint32_t numRead = iter1->read(data1);
-        BOOST_CHECK(numRead == 1);
-        delete iter1;
-    }
-
-    {
-        libpc::SequentialIterator* iter2 = descalingFilter.createSequentialIterator();
-        boost::uint32_t numRead = iter2->read(data2);
-        BOOST_CHECK(numRead == 1);
-        delete iter2;
-    }
+    const char* utm15_wkt = "PROJCS[\"NAD83 / UTM zone 15N\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.2572221010002,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4269\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-93],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"26915\"]]";
 
     const double preX = 470692.447538;
     const double preY = 4602888.904642;
@@ -150,27 +95,112 @@
     const double postY = 41.577148;
     const double postZ = 16.000000;
 
-    // note this file has only 1 points, so yes, the extent's mins and maxes are the same
-    const libpc::Bounds<double> oldBounds_ref(preX, preY, preZ, preX, preY, preZ);
-    const libpc::Bounds<double> newBounds_ref(postX, postY, postZ, postX, postY, postZ);
-    const libpc::Bounds<double>& oldBounds = reader1.getBounds();
-    const libpc::Bounds<double>& newBounds = descalingFilter.getBounds();
-    compareBounds(oldBounds_ref, oldBounds);
-    compareBounds(newBounds_ref, newBounds);
+    // we compute three answers:
+    //   (1) w/out reprojection
+    //   (2) with scaling and reproj
+    //   (3) with custom scaling and reproj
 
-    double x1, x2, y1, y2, z1, z2;


More information about the Liblas-commits mailing list