Working vertical point table.

filters/SMRFilter.cpp

@@ -83,14 +83,15 @@ struct SMRArgs
 {
     double m_cell;
     double m_slope;
-    double m_window;
     double m_scalar;
     double m_threshold;
     double m_cut;
     std::string m_dir;
     std::vector<DimRange> m_ignored;
     StringList m_returns;
     Segmentation::PointClasses m_classbits;
+    double m_window;
+    Arg *m_windowArg;
 };
 
 SMRFilter::SMRFilter() : m_args(new SMRArgs) {}
@@ -106,7 +107,6 @@ void SMRFilter::addArgs(ProgramArgs& args)
 {
     args.add("cell", "Cell size?", m_args->m_cell, 1.0);
     args.add("slope", "Percent slope?", m_args->m_slope, 0.15);
-    args.add("window", "Max window size?", m_args->m_window, 18.0);
     args.add("scalar", "Elevation scalar?", m_args->m_scalar, 1.25);
     args.add("threshold", "Elevation threshold?", m_args->m_threshold, 0.5);
     args.add("cut", "Cut net size?", m_args->m_cut, 0.0);
@@ -116,6 +116,8 @@ void SMRFilter::addArgs(ProgramArgs& args)
              {"last", "only"});
     args.add("classbits", "Ignore synthetic|keypoint|withheld "
         "classification bits?", m_args->m_classbits);
+    m_args->m_windowArg = &args.add("window", "Max window size?",
+        m_args->m_window);
 }
 
 void SMRFilter::addDimensions(PointLayoutPtr layout)
@@ -156,6 +158,9 @@ void SMRFilter::prepared(PointTableRef table)
             m_args->m_returns.clear();
         }
     }
+
+    if (!m_args->m_windowArg->set())
+        m_args->m_window = 18 * m_args->m_cell;
 }
 
 void SMRFilter::ready(PointTableRef table)
@@ -245,6 +250,10 @@ PointViewSet SMRFilter::run(PointViewPtr view)
         ((m_bounds.maxx - m_bounds.minx) / m_args->m_cell) + 1);
     m_rows = static_cast<int>(
         ((m_bounds.maxy - m_bounds.miny) / m_args->m_cell) + 1);
+    if (m_cols * m_rows < 10000)
+        log()->get(LogLevel::Warning) << "SMRF running with a small number "
+            "of cells (" << (m_cols * m_rows) << ").  Consider changing "
+            "cell size.\n";
 
     // Create raster of minimum Z values per element.
     std::vector<double> ZImin = createZImin(inlierView);
@@ -389,7 +398,7 @@ std::vector<int> SMRFilter::createLowMask(std::vector<double> const& ZImin)
     std::vector<double> negZImin;
     std::transform(ZImin.begin(), ZImin.end(), std::back_inserter(negZImin),
                    [](double v) { return -v; });
-    std::vector<int> LowV = progressiveFilter(negZImin, 5.0, 1.0);
+    std::vector<int> LowV = progressiveFilter(negZImin, 5.0, m_args->m_cell);
 
     if (!m_args->m_dir.empty())
     {
@@ -598,14 +607,19 @@ void SMRFilter::knnfill(PointViewPtr view, std::vector<double>& cz)
     // can construct a 2D KDIndex and perform nearest neighbor searches.
     PointViewPtr temp = view->makeNew();
     PointId i(0);
+    size_t nancells = 0;
+    size_t cells = m_cols * m_rows;
     for (int c = 0; c < m_cols; ++c)
     {
         for (int r = 0; r < m_rows; ++r)
         {
             size_t cell = c * m_rows + r;
             double val = cz[cell];
             if (std::isnan(val))
+            {
+                nancells++;
                 continue;
+            }
 
             temp->setField(Id::X, i,
                            m_bounds.minx + (c + 0.5) * m_args->m_cell);
@@ -615,6 +629,9 @@ void SMRFilter::knnfill(PointViewPtr view, std::vector<double>& cz)
             i++;
         }
     }
+    std::cerr << "Size = " << cells << "!\n";
+    std::cerr << "Nancells = " << nancells << "!\n";
+    std::cerr << "Filled cells = " << (cells - nancells) << "!\n";
 
     // https://github.com/PDAL/PDAL/issues/2794#issuecomment-625297062
     if (!temp->size())

pdal/DimDetail.hpp

@@ -44,7 +44,7 @@ namespace Dimension
 class Detail
 {
 public:
-    Detail() : m_id(Id::Unknown), m_offset(-1), m_type(Type::None)
+    Detail() : m_id(Id::Unknown), m_offsetOrOrder(-1), m_type(Type::None)
     {}
     //NOTE - This is strange, but for some reason things run faster with
     // this NOOP virtual dtor.  Perhaps it has something to do with
@@ -55,15 +55,19 @@ class Detail
     {}
 
     void setOffset(int offset)
-        { m_offset = offset; }
+        { m_offsetOrOrder = offset; }
+    void setOrder(int order)
+        { m_offsetOrOrder = order; }
     void setType(Type type)
         { m_type = type; }
     void setId(Id id)
         { m_id = id; }
     Id id() const
         { return m_id; }
     int offset() const
-        { return m_offset; }
+        { return m_offsetOrOrder; }
+    int order() const
+        { return m_offsetOrOrder; }
     Type type() const
         { return m_type; }
     size_t size() const
@@ -73,7 +77,7 @@ class Detail
 
 private:
     Id m_id;
-    int m_offset;
+    int m_offsetOrOrder;
     Type m_type;
 };
 typedef std::vector<Detail> DetailList;

pdal/EigenUtils.cpp

@@ -39,6 +39,7 @@
 #include <pdal/util/Bounds.hpp>
 #include <pdal/util/Utils.hpp>
 
+#include <array>
 #include <cfloat>
 #include <numeric>
 #include <vector>
@@ -254,7 +255,7 @@ Eigen::MatrixXd extendedLocalMinimum(const PointView& view, int rows, int cols,
 std::vector<double> dilateDiamond(std::vector<double> data, size_t rows, size_t cols, int iterations)
 {
     std::vector<double> out(data.size(), std::numeric_limits<double>::lowest());
-    std::vector<size_t> idx(5);
+    std::array<size_t, 5> idx;
 
     for (int iter = 0; iter < iterations; ++iter)
     {
@@ -294,7 +295,7 @@ std::vector<double> dilateDiamond(std::vector<double> data, size_t rows, size_t
 std::vector<double> erodeDiamond(std::vector<double> data, size_t rows, size_t cols, int iterations)
 {
     std::vector<double> out(data.size(), (std::numeric_limits<double>::max)());
-    std::vector<size_t> idx(5);
+    std::array<size_t, 5> idx;
 
     for (int iter = 0; iter < iterations; ++iter)
     {

pdal/PipelineManager.cpp

@@ -46,7 +46,7 @@ namespace pdal
 
 PipelineManager::PipelineManager(point_count_t streamLimit) :
     m_factory(new StageFactory),
-    m_tablePtr(new PointTable()), m_table(*m_tablePtr),
+    m_tablePtr(new ColumnPointTable()), m_table(*m_tablePtr),
     m_streamTablePtr(new FixedPointTable(streamLimit)),
     m_streamTable(*m_streamTablePtr),
     m_progressFd(-1), m_input(nullptr)

pdal/PipelineManager.hpp

@@ -156,7 +156,7 @@ class PDAL_DLL PipelineManager
     Options stageOptions(Stage& stage);
 
     std::unique_ptr<StageFactory> m_factory;
-    std::unique_ptr<PointTable> m_tablePtr;
+    std::unique_ptr<SimplePointTable> m_tablePtr;
     PointTableRef m_table;
     std::unique_ptr<FixedPointTable> m_streamTablePtr;
     StreamPointTable& m_streamTable;

pdal/PointLayout.hpp

@@ -229,6 +229,16 @@ class  PointLayout
         return &(m_detail[Utils::toNative(id)]);
     }
 
+    /**
+      Change from offset to order.  Only call during finalize.
+    */
+    void orderDimensions()
+    {
+        int order = 0;
+        for (Dimension::Id id : m_used)
+            m_detail[Utils::toNative(id)].setOrder(order++);
+    }
+
     /**
         Convert the point layout to a metadata format.
 

pdal/PointTable.cpp

@@ -125,23 +125,6 @@ char *PointTable::getPoint(PointId idx)
 }
 
 
-ContiguousPointTable::~ContiguousPointTable()
-{}
-
-
-PointId ContiguousPointTable::addPoint()
-{
-    m_buf.resize(pointsToBytes(m_numPts + 1));
-    return m_numPts++;
-}
-
-
-char *ContiguousPointTable::getPoint(PointId idx)
-{
-    return m_buf.data() + pointsToBytes(idx);
-}
-
-
 MetadataNode BasePointTable::toMetadata() const
 {
     return layout()->toMetadata();

pdal/PointTable.hpp

@@ -98,6 +98,7 @@ class PDAL_DLL BasePointTable : public PointContainer
 private:
     // Point data operations.
     virtual PointId addPoint() = 0;
+    virtual char *getDimension(const Dimension::Detail *d, PointId idx) = 0;
 
 protected:
     virtual char *getPoint(PointId idx) = 0;
@@ -166,25 +167,42 @@ class PDAL_DLL PointTable : public SimplePointTable
     PointLayout m_layout;
 };
 
-class PDAL_DLL ContiguousPointTable : public SimplePointTable
+// This provides a context for processing a set of points and allows the library
+// to be used to process multiple point sets simultaneously.
+class PDAL_DLL ColumnPointTable : public SimplePointTable
 {
 private:
-    std::vector<char> m_buf;
+    // Point storage.
+    using DimBlockList = std::vector<char *>;
+    using MemBlocks = std::vector<DimBlockList>;
+
+    // List of dimension memory block lists.
+    MemBlocks m_blocks;
     point_count_t m_numPts;
+    static const point_count_t m_blockPtCnt = 16384;
 
 public:
-    ContiguousPointTable() : SimplePointTable(m_layout), m_numPts(0)
+    ColumnPointTable() : SimplePointTable(m_layout), m_numPts(0)
         {}
-    virtual ~ContiguousPointTable();
+    virtual ~ColumnPointTable();
     virtual bool supportsView() const
         { return true; }
-
-protected:
-    virtual char *getPoint(PointId idx);
+    virtual void finalize();
+    virtual char *getPoint(PointId idx)
+        { return nullptr; }
 
 private:
+    virtual void setFieldInternal(Dimension::Id id, PointId idx,
+        const void *value);
+    virtual void getFieldInternal(Dimension::Id id, PointId idx,
+        void *value) const;
+
     virtual PointId addPoint();
 
+    // Hide base class calls for now.
+    const char *getDimension(const Dimension::Detail *d, PointId idx) const;
+    char *getDimension(const Dimension::Detail *d, PointId idx);
+
     PointLayout m_layout;
 };
 

pdal/PointView.cpp

@@ -45,13 +45,14 @@ namespace pdal
 int PointView::m_lastId = 0;
 
 PointView::PointView(PointTableRef pointTable) : m_pointTable(pointTable),
-m_size(0), m_id(0)
+    m_layout(pointTable.layout()), m_size(0), m_id(0)
 {
 	m_id = ++m_lastId;
 }
 
 PointView::PointView(PointTableRef pointTable, const SpatialReference& srs) :
-	m_pointTable(pointTable), m_size(0), m_id(0), m_spatialReference(srs)
+	m_pointTable(pointTable), m_layout(pointTable.layout()), m_size(0),
+    m_id(0), m_spatialReference(srs)
 {
 	m_id = ++m_lastId;
 }
@@ -73,6 +74,22 @@ PointViewIter PointView::end()
 }
 
 
+PointId PointView::tableId(PointId idx)
+{
+    if (idx > size())
+        throw pdal_error("Point index must increment.");
+    if (idx == size())
+    {
+        PointId rawId = m_pointTable.addPoint();
+        m_index.push_back(rawId);
+        m_size++;
+        assert(m_temps.empty());
+        return rawId;
+    }
+    return m_index[idx];
+}
+
+
 void PointView::setFieldInternal(Dimension::Id dim, PointId idx,
     const void *buf)
 {