[filter] Reverse Morton (#1786)

* Add revert Morton algorithm

* Add test

* Rename option from revert to reverse

* Update documentation

* Some minor changes

doc/stages/filters.mortonorder.rst

@@ -5,8 +5,25 @@ filters.mortonorder
 
 Sorts the XY data using `Morton ordering`_.
 
+It's also possible to compute a reverse Morton code by reading the binary
+representation from the end to the beginning. This way, points are sorted
+with a good dispersement. For example, by successively selecting N
+representative points within tiles:
+
+.. figure:: filters.mortonorder.img1.png
+    :scale: 100 %
+    :alt: Reverse Morton indexing
+
 .. _`Morton ordering`: http://en.wikipedia.org/wiki/Z-order_curve
 
+.. seealso::
+
+    See `LOPoCS`_ and `pgmorton`_ for some use case examples of the
+    Reverse Morton algorithm.
+
+.. _`pgmorton`: https://github.com/Oslandia/pgmorton
+.. _`LOPoCS`: https://github.com/Oslandia/lopocs
+
 .. embed::
 
 Example
@@ -18,7 +35,8 @@ Example
       "pipeline":[
         "uncompressed.las",
         {
-          "type":"filters.mortonorder"
+          "type":"filters.mortonorder",
+          "reverse":"false"
         },
         {
           "type":"writers.las",

filters/MortonOrderFilter.cpp

@@ -53,6 +53,11 @@ CREATE_STATIC_PLUGIN(1, 0, MortonOrderFilter, Filter, s_info)
 
 std::string MortonOrderFilter::getName() const { return s_info.name; }
 
+void MortonOrderFilter::addArgs(ProgramArgs& args)
+{
+    args.add("reverse", "Reverse Morton", m_reverse, false);
+}
+
 //This used to be a lambda, but the VS compiler exploded, I guess.
 typedef std::pair<double, double> Coord;
 namespace
@@ -87,7 +92,91 @@ class CmpZOrder
 };
 }
 
-PointViewSet MortonOrderFilter::run(PointViewPtr inView)
+class ReverseZOrder
+{
+public:
+    static uint32_t encode_morton(uint32_t x, uint32_t y)
+    {
+        return (part1_by1(y) << 1) + part1_by1(x);
+    }
+
+    static uint32_t reverse_morton(uint32_t index)
+    {
+        index = ((index >> 1) & 0x55555555u) | ((index & 0x55555555u) << 1);
+        index = ((index >> 2) & 0x33333333u) | ((index & 0x33333333u) << 2);
+        index = ((index >> 4) & 0x0f0f0f0fu) | ((index & 0x0f0f0f0fu) << 4);
+        index = ((index >> 8) & 0x00ff00ffu) | ((index & 0x00ff00ffu) << 8);
+        index = ((index >> 16) & 0xffffu) | ((index & 0xffffu) << 16);
+        return index;
+    }
+
+private:
+    static uint32_t part1_by1(uint32_t x)
+    {
+        x &= 0x0000ffff;
+        x = (x ^ (x <<  8)) & 0x00ff00ff;
+        x = (x ^ (x <<  4)) & 0x0f0f0f0f;
+        x = (x ^ (x <<  2)) & 0x33333333;
+        x = (x ^ (x <<  1)) & 0x55555555;
+        return x;
+    }
+
+    static uint32_t part1_by2(uint32_t x)
+    {
+        x &= 0x000003ff;
+        x = (x ^ (x << 16)) & 0xff0000ff;
+        x = (x ^ (x <<  8)) & 0x0300f00f;
+        x = (x ^ (x <<  4)) & 0x030c30c3;
+        x = (x ^ (x <<  2)) & 0x09249249;
+        return x;
+    }
+};
+
+PointViewSet MortonOrderFilter::reverseMorton(PointViewPtr inView)
+{
+    const int32_t cell = sqrt(inView->size());
+
+    // compute range
+    BOX2D buffer_bounds;
+    inView->calculateBounds(buffer_bounds);
+    const double xrange = buffer_bounds.maxx - buffer_bounds.minx;
+    const double yrange = buffer_bounds.maxy - buffer_bounds.miny;
+
+    const double cell_width = xrange / cell;
+    const double cell_height = yrange / cell;
+
+    // compute reverse morton code for each point
+    std::multimap<uint32_t, PointId> codes;
+    for (PointId idx = 0; idx < inView->size(); idx++)
+    {
+        const double x = inView->getFieldAs<double>(Dimension::Id::X, idx);
+        const int32_t xpos = floor((x - buffer_bounds.minx) / cell_width);
+
+        const double y = inView->getFieldAs<double>(Dimension::Id::Y, idx);
+        const int32_t ypos = floor((y - buffer_bounds.miny) / cell_height);
+
+        const uint32_t code = ReverseZOrder::encode_morton(xpos, ypos);
+        const uint32_t reverse = ReverseZOrder::reverse_morton( code );
+
+        codes.insert( std::pair<uint32_t, PointId>(reverse, idx) );
+    }
+
+    // a map is yet order by key so its naturally ordered by lod
+    std::multimap<uint32_t, PointId>::iterator it;
+    PointViewPtr outView = inView->makeNew();
+    for (it = codes.begin(); it != codes.end(); ++it)
+    {
+        outView->appendPoint(*inView, it->second);
+    }
+
+    // build output view
+    PointViewSet viewSet;
+    viewSet.insert(outView);
+
+    return viewSet;
+}
+
+PointViewSet MortonOrderFilter::morton(PointViewPtr inView)
 {
     PointViewSet viewSet;
     if (!inView->size())
@@ -121,4 +210,16 @@ PointViewSet MortonOrderFilter::run(PointViewPtr inView)
     return viewSet;
 }
 
+PointViewSet MortonOrderFilter::run(PointViewPtr inView)
+{
+    if ( m_reverse )
+    {
+        return reverseMorton( inView );
+    }
+    else
+    {
+        return morton( inView );
+    }
+}
+
 } // pdal

filters/MortonOrderFilter.hpp

@@ -56,7 +56,13 @@ class PDAL_DLL MortonOrderFilter : public pdal::Filter
     std::string getName() const;
 
 private:
+    virtual void addArgs(ProgramArgs& args);
     virtual PointViewSet run(PointViewPtr view);
+
+    PointViewSet reverseMorton(PointViewPtr view);
+    PointViewSet morton(PointViewPtr view);
+
+    bool m_reverse = false;
 };
 
 } // namespace pdal

test/unit/CMakeLists.txt

@@ -105,6 +105,7 @@ PDAL_ADD_TEST(pdal_filters_ferry_test FILES filters/FerryFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_groupby_test FILES filters/GroupByFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_locate_test FILES filters/LocateFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_merge_test FILES filters/MergeTest.cpp)
+PDAL_ADD_TEST(pdal_morton_order_test FILES filters/MortonOrderTest.cpp)
 PDAL_ADD_TEST(pdal_filters_additional_merge_test FILES
     filters/AdditionalMergeTest.cpp)
 target_include_directories(pdal_filters_additional_merge_test PRIVATE ${PDAL_JSONCPP_INCLUDE_DIR})

test/unit/filters/MortonOrderTest.cpp

@@ -0,0 +1,115 @@
+/******************************************************************************
+* Copyright (c) 2018, Paul Blottiere (blottiere.paul@gmail.com)
+*
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following
+* conditions are met:
+*
+*     * Redistributions of source code must retain the above copyright
+*       notice, this list of conditions and the following disclaimer.
+*     * Redistributions in binary form must reproduce the above copyright
+*       notice, this list of conditions and the following disclaimer in
+*       the documentation and/or other materials provided
+*       with the distribution.
+*     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+*       names of its contributors may be used to endorse or promote
+*       products derived from this software without specific prior
+*       written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+****************************************************************************/
+
+#include <pdal/pdal_test_main.hpp>
+
+#include <io/BufferReader.hpp>
+#include <filters/MortonOrderFilter.hpp>
+
+#include "Support.hpp"
+
+using namespace pdal;
+
+TEST(MortonOrderTest, test_code)
+{
+    PointTable table;
+    table.layout()->registerDim(Dimension::Id::X);
+    table.layout()->registerDim(Dimension::Id::Y);
+
+    PointViewPtr view(new PointView(table));
+
+    int n = 0;
+    for ( int i = 0; i < 4; i++ )
+    {
+        for ( int j = 0; j < 4; j++ )
+        {
+            view->setField(Dimension::Id::X, n, i);
+            view->setField(Dimension::Id::Y, n, j);
+            n += 1;
+        }
+    }
+
+    BufferReader r;
+    r.addView(view);
+
+    MortonOrderFilter filter;
+    Options o;
+    o.add("reverse", "true");
+    filter.setInput(r);
+    filter.setOptions(o);
+
+    filter.prepare(table);
+    PointViewSet s = filter.execute(table);
+    PointViewPtr outView = *s.begin();
+
+    /*
+     ___________________
+    |    |    |    |    |
+    | 0  | 9  | 6  | 2  |
+    |____|____|____|____|
+    |    |    |    |    |
+    | 12 | 15 | 11 | 10 |
+    |____|____|____|____|
+    |    |    |    |    |
+    | 4  | 14 | 8  | 5  |
+    |____|____|____|____|
+    |    |    |    |    |
+    | 1  | 13 | 7  | 3  |
+    |____|____|____|____|
+    */
+
+    // index 0
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 0), 0);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 0), 0);
+
+    // index 1
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 1), 0);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 1), 3);
+
+    // index 2
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 2), 3);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 2), 0);
+
+    // index 3
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 3), 3);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 3), 3);
+
+    // index 4
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 4), 0);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 4), 2);
+
+    // index 5
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::X, 5), 3);
+    EXPECT_EQ(outView->getFieldAs<double>(Dimension::Id::Y, 5), 2);
+}