Merge remote-tracking branch 'origin/master' into matlab-updates

cmake/options.cmake

@@ -8,9 +8,9 @@ add_feature_info("Bash completion" WITH_COMPLETION
     "completion for PDAL command line")
 
 option(BUILD_PLUGIN_CPD
-    "Choose if Coherent Point Drift kernel is built" FALSE)
+    "Choose if the cpd filter should be built" FALSE)
 add_feature_info("CPD plugin" BUILD_PLUGIN_CPD
-    "run Coherent Point Drift on two datasets")
+    "Coherent Point Drift (CPD) computes rigid or nonrigid transformations between point sets")
 
 option(BUILD_PLUGIN_GEOWAVE
     "Choose if GeoWave support should be built" FALSE)

doc/stages/filters.cpd.rst

@@ -0,0 +1,73 @@
+.. _filters.cpd:
+
+filters.cpd
+==============
+
+The CPD filter uses the Coherent Point Drift :cite:`Myronenko` algorithm to compute a rigid, nonrigid, or affine transformation between datasets.
+The rigid and affine are what you'd expect; the nonrigid transformation uses Motion Coherence Theory :cite:`Yuille1998` to "bend" the points to find a best alignment.
+
+.. note::
+
+    CPD is computationally intensive and can be slow when working with many points (i.e. > 10 000).
+    Nonrigid is significatly slower than rigid and affine.
+
+The first input to the change filter are considered the "fixed" points, and all subsequent inputs are "moving" points.
+The output from the change filter are the "moving" points after the calculated transformation has been applied, one point view per input.
+Any additional information about the cpd registration, e.g. the rigid transformation matrix, will be placed in the stage's metadata.
+
+Examples
+--------
+
+.. code-block:: json
+
+    {
+        "pipeline": [
+            "fixed.las",
+            "moving.las",
+            {
+                "type": "filters.rigid",
+                "method": "rigid"
+            },
+            "output.las"
+        ]
+    }
+
+If "method" is not provided, the cpd filter will default to using the rigid registration method.
+To get the transform matrix, you'll need to use the ``--metadata`` option:
+
+::
+
+    $ pdal pipeline cpd-pipeline.json --metadata cpd-metadata.json
+
+The metadata output might start something like:
+
+.. code-block:: json
+
+    {
+        "stages":
+        {
+            "filters.cpd":
+            {
+                "iterations": 10,
+                "method": "rigid",
+                "runtime": 0.003839,
+                "sigma2": 5.684342128e-16,
+                "transform": "           1 -6.21722e-17  1.30104e-18  5.29303e-11-8.99346e-17            1  2.60209e-18 -3.49247e-10 -2.1684e-19  1.73472e-18            1 -1.53477e-12           0            0            0            1"
+            },
+        },
+
+.. seealso::
+
+    :ref:`filters.transformation` to apply a transform to other points.
+
+Options
+--------
+
+method
+    Change detection method to use.
+    Valid values are "rigid", "affine", and "nonrigid".
+    [Default: **rigid**]
+
+.. _Coherent Point Drift (CPD): https://github.com/gadomski/cpd
+
+.. bibliography:: references.bib

doc/stages/filters.icp.rst

@@ -0,0 +1,56 @@
+.. _filters.icp:
+
+filters.icp
+==============
+
+The ICP filter uses the `PCL's Iterative Closest Point (ICP)`_ algorithm to calculate a rigid (rotation and translation) transformation that best aligns two datasets.
+The first input to the ICP filter is considered the "fixed" points, and all subsequent points are "moving" points.
+The output from the change filter are the "moving" points after the calculated transformation has been applied, one point view per input.
+The transformation matrix is inserted into the stage's metadata.
+
+Examples
+--------
+
+.. code-block:: json
+
+    {
+        "pipeline": [
+            "fixed.las",
+            "moving.las",
+            {
+                "type": "filters.icp"
+            },
+            "output.las"
+        ]
+    }
+
+To get the transform matrix, you'll need to use the ``--metadata`` option:
+
+::
+
+    $ pdal pipeline icp-pipeline.json --metadata icp-metadata.json
+
+The metadata output might start something like:
+
+.. code-block:: json
+
+    {
+        "stages":
+        {
+            "filters.icp":
+            {
+                "converged": true,
+                "fitness": 0.01953125097,
+                "transform": "           1  2.60209e-18 -1.97906e-09       -0.375  8.9407e-08            1  5.58794e-09      -0.5625 6.98492e -10 -5.58794e-09            1   0.00411987           0            0            0            1"
+            }
+
+.. seealso::
+
+    :ref:`filters.transformation` to apply a transform to other points.
+
+Options
+--------
+
+None.
+
+.. _PCL's Iterative Closest Point (ICP): http://docs.pointclouds.org/trunk/classpcl_1_1_iterative_closest_point.html

doc/stages/references.bib

@@ -0,0 +1,21 @@
+@article{Myronenko,
+abstract = {Point set registration is a key component in many computer vision tasks. The goal of point set registration is to assign correspondences between two sets of points and to recover the transformation that maps one point set to the other. Multiple factors, including an unknown nonrigid spatial transformation, large dimensionality of point set, noise, and outliers, make the point set registration a challenging problem. We introduce a probabilistic method, called the Coherent Point Drift (CPD) algorithm, for both rigid and nonrigid point set registration. We consider the alignment of two point sets as a probability density estimation problem. We fit the Gaussian mixture model (GMM) centroids (representing the first point set) to the data (the second point set) by maximizing the likelihood. We force the GMM centroids to move coherently as a group to preserve the topological structure of the point sets. In the rigid case, we impose the coherence constraint by reparameterization of GMM centroid locations with rigid parameters and derive a closed form solution of the maximization step of the EM algorithm in arbitrary dimensions. In the nonrigid case, we impose the coherence constraint by regularizing the displacement field and using the variational calculus to derive the optimal transformation. We also introduce a fast algorithm that reduces the method computation complexity to linear. We test the CPD algorithm for both rigid and nonrigid transformations in the presence of noise, outliers, and missing points, where CPD shows accurate results and outperforms current state-of-the-art methods.},
+author = {Myronenko, Andriy and Song, Xubo},
+issn = {1939-3539},
+journal = {IEEE transactions on pattern analysis and machine intelligence},
+month = {dec},
+number = {12},
+pages = {2262--75},
+pmid = {20975122},
+publisher = {Institute of Electrical {\&} Electronics Engineers (IEEE)},
+title = {{Point set registration: coherent point drift.}},
+volume = {32},
+year = {2010}
+}
+
+@article{Yuille1998,
+author = {Yuille, Alan L. and Grzywacz, Norberto M.},
+journal = {Second International Conference on Computer Vision},
+title = {{The Motion Coherence Theory}},
+year = {1988}
+}

pdal/EigenUtils.hpp

@@ -36,6 +36,7 @@
 
 #include <pdal/pdal_internal.hpp>
 #include <pdal/util/Bounds.hpp>
+#include <pdal/Metadata.hpp>
 
 #include <Eigen/Dense>
 
@@ -949,4 +950,49 @@ PDAL_DLL Eigen::MatrixXd computeSpline(Eigen::MatrixXd x, Eigen::MatrixXd y,
 
 } // namespace eigen
 
+namespace Utils
+{
+
+template <>
+inline bool fromString<Eigen::MatrixXd>(const std::string& s, Eigen::MatrixXd& matrix) {
+    std::stringstream ss(s);
+    std::string line;
+    std::vector<std::vector<double>> rows;
+    while (std::getline(ss, line)) {
+        std::vector<double> row;
+        std::stringstream ss(line);
+        double n;
+        while (ss >> n) {
+            row.push_back(n);
+            if (ss.peek() == ',' || ss.peek() == ' ') {
+                ss.ignore();
+            }
+        }
+        if (!rows.empty() && rows.back().size() != row.size()) {
+            return false;
+        }
+        rows.push_back(row);
+    }
+    if (rows.empty()) {
+        return true;
+    }
+    size_t nrows = rows.size();
+    size_t ncols = rows[0].size();
+    matrix.resize(nrows, ncols);
+    for (size_t i = 0; i < nrows; ++i) {
+        for (size_t j = 0; j < ncols; ++j) {
+            matrix(i, j) = rows[i][j];
+        }
+    }
+    return true;
+}
+}
+
+template <>
+inline void MetadataNodeImpl::setValue(const Eigen::MatrixXd& matrix)
+{
+    m_type = "matrix";
+    m_value = Utils::toString(matrix);
+}
+
 } // namespace pdal

pdal/Metadata.hpp

@@ -488,7 +488,7 @@ class PDAL_DLL MetadataNode
 
         std::string v(Utils::escapeJSON(value()));
         if (m_impl->m_type == "string" || m_impl->m_type == "base64Binary" ||
-            m_impl->m_type == "uuid")
+            m_impl->m_type == "uuid" || m_impl->m_type == "matrix")
         {
             std::string val("\"");
             val += escapeQuotes(v) + "\"";

plugins/cpd/CMakeLists.txt

@@ -1,14 +1,21 @@
-find_package(Cpd 0.5 REQUIRED)
+set(Cpd_VERSION 0.5)
 
-set_package_properties(Cpd PROPERTIES
-    DESCRIPTION "Coherent Point Drift"
-    URL "https://github.com/gadomski/cpd"
-    TYPE OPTIONAL
-    PURPOSE "Register two point sets using the Coherent Point Drift algorithm"
-    )
+find_package(Cpd ${Cpd_VERSION} REQUIRED)
+option(BUILD_PLUGIN_CPD "Build Coherent Point Drift support" ${Cpd_FOUND})
+
+set(files filters/CpdFilter.cpp)
+set(include_dirs "${CMAKE_CURRENT_LIST_DIR}" "${PDAL_VENDOR_DIR}/eigen")
 
-pdal_add_plugin(cpd_kernel_lib_name kernel cpd
-    FILES kernel/CpdKernel.cpp
+PDAL_ADD_PLUGIN(filter_libname filter cpd
+    FILES filters/CpdFilter.cpp
     LINK_WITH Cpd::Library-C++
     )
-target_include_directories(${cpd_kernel_lib_name} PRIVATE ${CMAKE_CURRENT_LIST_DIR})
+target_include_directories(${filter_libname} PRIVATE "${include_dirs}")
+
+if(${WITH_TESTS})
+    PDAL_ADD_TEST(pdal_filters_cpd_test
+        FILES test/CpdFilterTest.cpp
+        LINK_WITH ${filter_libname}
+        )
+    target_include_directories(pdal_filters_cpd_test PRIVATE "${include_dirs}")
+endif()