Merge remote-tracking branch 'origin/master' into ept

CMakeLists.txt

@@ -50,7 +50,6 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
 # Name of C++ library
-
 #
 # To facilitate one-library linking, we do special things for various platforms.
 #

doc/stages/filters.optimalneighborhood.rst

@@ -44,4 +44,4 @@ min_k
 
 max_k
   The maximum number of k nearest neighbors to consider for optimal
-  neighborhood selection. [Default: 100]
+  neighborhood selection. [Default: 14]

doc/stages/filters.sample.rst

@@ -3,8 +3,8 @@
 filters.sample
 ===============================================================================
 
-The **Sample Filter** performs Poisson sampling of the input ``PointView``. The
-The practice of performing Poisson sampling via "Dart Throwing" was introduced
+The **Sample Filter** performs Poisson sampling of the input ``PointView``. The 
+practice of performing Poisson sampling via "Dart Throwing" was introduced
 in the mid-1980's by [Cook1986]_ and [Dippe1985]_, and has been applied to
 point clouds in other software [Mesh2009]_.
 

doc/stages/filters.transformation.rst

@@ -3,12 +3,9 @@
 filters.transformation
 ======================
 
-The transformation filter applies an arbitrary rotation+translation
+The transformation filter applies an arbitrary homography
 transformation, represented as a 4x4 matrix_, to each xyz triplet.
 
-The filter does *no* checking to ensure the matrix is a valid affine
-transformation.
-
 .. note::
 
     The transformation filter does not apply or consider any spatial
@@ -59,10 +56,11 @@ A full tutorial about transformation matrices is beyond the scope of this
 documentation. Instead, we will provide a few pointers to introduce core
 concepts, especially as pertains to PDAL's handling of the ``matrix`` argument.
 
-Transformations in a 3-dimensional coordinate system can be represented as an
-affine transformation using homogeneous coordinates. This 4x4 matrix can
-represent transformations describing operations like translation, rotation, and
-scaling of coordinates.
+Transformations in a 3-dimensional coordinate system can be represented
+as a homography transformation using homogeneous coordinates. This 4x4
+matrix can represent affine transformations describing operations like
+translation, rotation, and scaling of coordinates.  In addition it can
+represent perspective transformations modeling a pinhole camera.
 
 The transformation filter's ``matrix`` argument is a space delimited, 16
 element string. This string is simply a row-major representation of the 4x4

doc/stages/readers.rst

@@ -32,7 +32,6 @@ like :ref:`readers.pgpointcloud`, or a network service like :ref:`readers.ept`.
    readers.mrsid
    readers.nitf
    readers.numpy
-   readers.oci
    readers.optech
    readers.pcd
    readers.pgpointcloud
@@ -106,9 +105,6 @@ like :ref:`readers.pgpointcloud`, or a network service like :ref:`readers.ept`.
 :ref:`readers.numpy`
     Read point cloud data from Numpy ``.npy`` files.
 
-:ref:`readers.oci`
-    Read data from Oracle point cloud databases. [deprecated]
-
 :ref:`readers.optech`
     Read Optech Corrected Sensor Data (.csd) files.
 

doc/stages/writers.rst

@@ -25,7 +25,6 @@ dimension type, while others only understand fixed dimension names.
    writers.matlab
    writers.nitf
    writers.null
-   writers.oci
    writers.ogr
    writers.pcd
    writers.pgpointcloud
@@ -66,9 +65,6 @@ dimension type, while others only understand fixed dimension names.
     Provides a sink for points in a pipeline. It's the same as sending pipeline
     output to /dev/null.
 
-:ref:`writers.oci`
-    Write data to Oracle point cloud databases. [deprecated]
-
 :ref:`writers.ogr`
     Write a point cloud as a set of OGR points/multipoints
 

filters/IterativeClosestPoint.cpp

@@ -165,7 +165,11 @@ PointViewPtr IterativeClosestPoint::icp(PointViewPtr fixed,
         fixed_idx.reserve(tempMovingTransformed->size());
         moving_idx.reserve(tempMovingTransformed->size());
         double mse(0.0);
-        double sqr_maxdist = m_maxdist * m_maxdist;
+        double sqr_maxdist;
+        if (m_maxdistArg->set())
+            sqr_maxdist = m_maxdist * m_maxdist;
+        else
+            sqr_maxdist = (std::numeric_limits<double>::max)();
 
         // For every point in the centered, moving PointView, find the nearest
         // neighbor in the centered fixed PointView. Record the indices of each
@@ -178,18 +182,14 @@ PointViewPtr IterativeClosestPoint::icp(PointViewPtr fixed,
             std::vector<double> sqr_dists(1);
             kd_fixed.knnSearch(p, 1, &indices, &sqr_dists);
 
-            // In the PCL code, there would've been a check that the square
-            // distance did not exceed a threshold value.
-            if (m_maxdistArg->set())
+            // Check that the square distance does not exceed threshold value.
+            if (sqr_dists[0] < sqr_maxdist)
             {
-                if (sqr_dists[0] > sqr_maxdist)
-                    continue;
+                // Store the indices of the correspondence and update the MSE.
+                moving_idx.push_back(p.pointId());
+                fixed_idx.push_back(indices[0]);
+                mse += std::sqrt(sqr_dists[0]);
             }
-
-            // Store the indices of the correspondence and update the MSE.
-            moving_idx.push_back(p.pointId());
-            fixed_idx.push_back(indices[0]);
-            mse += std::sqrt(sqr_dists[0]);
         }
 
         // Finalize and log the MSE.
@@ -281,15 +281,27 @@ PointViewPtr IterativeClosestPoint::icp(PointViewPtr fixed,
     // Compute the MSE one last time, using the unaltered, fixed PointView and
     // the transformed, moving PointView.
     double mse(0.0);
+    size_t mse_n(0);
+    double sqr_maxdist;
+    if (m_maxdistArg->set())
+        sqr_maxdist = m_maxdist * m_maxdist;
+    else
+        sqr_maxdist = (std::numeric_limits<double>::max)();
     KD3Index& kd_fixed_orig = fixed->build3dIndex();
     for (PointRef p : *moving)
     {
         PointIdList indices(1);
         std::vector<double> sqr_dists(1);
         kd_fixed_orig.knnSearch(p, 1, &indices, &sqr_dists);
-        mse += std::sqrt(sqr_dists[0]);
+
+        // Check that the square distance does not exceed threshold value.
+        if (sqr_dists[0] < sqr_maxdist)
+        {
+            mse_n++;
+            mse += std::sqrt(sqr_dists[0]);
+        }
     }
-    mse /= moving->size();
+    mse /= mse_n;
     log()->get(LogLevel::Debug2) << "MSE: " << mse << std::endl;
 
     // Transformation to demean coords

filters/TransformationFilter.cpp

@@ -190,15 +190,16 @@ bool TransformationFilter::processOne(PointRef& point)
     double x = point.getFieldAs<double>(Dimension::Id::X);
     double y = point.getFieldAs<double>(Dimension::Id::Y);
     double z = point.getFieldAs<double>(Dimension::Id::Z);
+    double s = x * matrix[12] + y * matrix[13] + z * matrix[14] + matrix[15];
 
     point.setField(Dimension::Id::X,
-        x * matrix[0] + y * matrix[1] + z * matrix[2] + matrix[3]);
+       (x * matrix[0] + y * matrix[1] + z * matrix[2] + matrix[3]) / s);
 
     point.setField(Dimension::Id::Y,
-        x * matrix[4] + y * matrix[5] + z * matrix[6] + matrix[7]);
+       (x * matrix[4] + y * matrix[5] + z * matrix[6] + matrix[7]) / s);
 
     point.setField(Dimension::Id::Z,
-        x * matrix[8] + y * matrix[9] + z * matrix[10] + matrix[11]);
+       (x * matrix[8] + y * matrix[9] + z * matrix[10] + matrix[11]) / s);
     return true;
 }