Python filter with MaskedArray.

pdal/PipelineManager.cpp

@@ -45,7 +45,8 @@ namespace pdal
 {
 
 PipelineManager::PipelineManager() : m_factory(new StageFactory),
-    m_tablePtr(new PointTable()), m_table(*m_tablePtr),
+//    m_tablePtr(new PointTable()), m_table(*m_tablePtr),
+    m_tablePtr(new ContiguousPointTable()), m_table(*m_tablePtr),
     m_progressFd(-1), m_input(nullptr)
 {}
 

pdal/PipelineManager.hpp

@@ -145,7 +145,8 @@ class PDAL_DLL PipelineManager
     Options stageOptions(Stage& stage);
 
     std::unique_ptr<StageFactory> m_factory;
-    std::unique_ptr<PointTable> m_tablePtr;
+//    std::unique_ptr<PointTable> m_tablePtr;
+    std::unique_ptr<SimplePointTable> m_tablePtr;
     PointTableRef m_table;
     Options m_commonOptions;
     OptionsMap m_stageOptions;

pdal/PointTable.hpp

@@ -178,6 +178,10 @@ class PDAL_DLL ContiguousPointTable : public SimplePointTable
     virtual ~ContiguousPointTable();
     virtual bool supportsView() const
         { return true; }
+    char *data()
+        { return getPoint(0); }
+    point_count_t size() const
+        { return m_numPts; }
 
 protected:
     virtual char *getPoint(PointId idx);

plugins/python/PythonPointView.hpp

@@ -0,0 +1,69 @@
+/******************************************************************************
+* Copyright (c) 2018, Hobu Inc.
+*
+* 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.
+****************************************************************************/
+
+#pragma once
+
+#include <memory>
+
+#include <pdal/PointView.hpp>
+
+namespace pdal
+{
+
+class PythonPointView;
+
+typedef std::shared_ptr<PythonPointView> PythonPointViewPtr;
+
+// The intent of this is NOT for it to be constructed, but to provide a
+// way, through a pointer or reference, to gain access to otherwise
+// inaccessible PointView data.
+class PDAL_DLL PythonPointView : public PointView
+{
+public:
+    PythonPointView() = delete;
+
+    // For testing only.
+    PointId index(PointId id) const
+        { return m_index[id]; }
+
+    static PythonPointViewPtr convert(PointViewPtr& v)
+    {
+        // Not until C++17
+        //return std::reinterpret_pointer_cast<PythonPointView>(v);
+        PointView *t = v.get();
+        return PythonPointViewPtr(v, reinterpret_cast<PythonPointView *>(t));
+    }
+};
+
+} // namespace pdal

plugins/python/filters/PythonFilter.cpp

@@ -76,20 +76,30 @@ void PythonFilter::ready(PointTableRef table)
 {
     if (m_source.empty())
         m_source = FileUtils::readFileIntoString(m_scriptFile);
-    static_cast<plang::Environment*>(plang::Environment::get())->set_stdout(log()->getLogStream());
+    plang::Environment::get()->set_stdout(log()->getLogStream());
     m_script = new plang::Script(m_source, m_module, m_function);
     m_pythonMethod = new plang::Invocation(*m_script);
     m_pythonMethod->compile();
     m_totalMetadata = table.metadata();
+
+    ContiguousPointTable *cpt = dynamic_cast<ContiguousPointTable *>(&table);
+    m_contiguous = (bool)cpt;
+    if (cpt)
+        m_pythonMethod->createArray(cpt->data(), cpt->size(), cpt->layout());
 }
 
 
 PointViewSet PythonFilter::run(PointViewPtr view)
 {
     log()->get(LogLevel::Debug5) << "filters.python " << *m_script <<
         " processing " << view->size() << " points." << std::endl;
-    m_pythonMethod->resetArguments();
-    m_pythonMethod->begin(*view, m_totalMetadata);
+    if (m_contiguous)
+        m_pythonMethod->setMask(PythonPointView::convert(view));
+    else
+    {
+        m_pythonMethod->resetArguments();
+        m_pythonMethod->begin(*view, m_totalMetadata);
+    }
 
     if (!m_pdalargs.empty())
     {
@@ -106,8 +116,8 @@ PointViewSet PythonFilter::run(PointViewPtr view)
         PointViewPtr outview = view->makeNew();
 
         size_t arrSize(0);
-        void *pydata =
-            m_pythonMethod->extractResult("Mask", Dimension::Type::Unsigned8, arrSize);
+        void *pydata = m_pythonMethod->extractResult("Mask",
+            Dimension::Type::Unsigned8, arrSize);
         char *ok = (char *)pydata;
         for (PointId idx = 0; idx < arrSize; ++idx)
             if (*ok++)

plugins/python/filters/PythonFilter.hpp

@@ -60,6 +60,7 @@ class PDAL_DLL PythonFilter : public Filter
     std::string m_module;
     std::string m_function;
     StringList m_addDimensions;
+    bool m_contiguous;
 
     virtual void addArgs(ProgramArgs& args);
     virtual void addDimensions(PointLayoutPtr layout);

plugins/python/plang/Environment.cpp

@@ -309,7 +309,7 @@ void addMetadata(PyObject *dict, MetadataNode m)
     }
 }
 
-int Environment::getPythonDataType(Dimension::Type t)
+int Environment::pythonType(Dimension::Type t)
 {
     using namespace Dimension;
 

plugins/python/plang/Environment.hpp

@@ -73,7 +73,7 @@ class PDAL_DLL Environment
 
     static EnvironmentPtr get();
 
-    static int getPythonDataType(Dimension::Type t);
+    static int pythonType(Dimension::Type t);
     static pdal::Dimension::Type getPDALDataType(int t);
 
 private:

plugins/python/plang/Invocation.cpp

@@ -69,6 +69,7 @@ Invocation::Invocation(const Script& script)
     , m_module(NULL)
     , m_dictionary(NULL)
     , m_function(NULL)
+    , m_directArray(NULL)
     , m_varsIn(NULL)
     , m_varsOut(NULL)
     , m_scriptArgs(NULL)
@@ -153,7 +154,7 @@ void Invocation::insertArgument(std::string const& name, uint8_t* data,
     int flags = NPY_CARRAY;
 #endif
 
-    const int pyDataType = plang::Environment::getPythonDataType(t);
+    const int pyDataType = Environment::pythonType(t);
 
     PyObject* pyArray = PyArray_New(&PyArray_Type, nd, dims, pyDataType,
         strides, data, 0, flags, NULL);
@@ -167,15 +168,15 @@ void *Invocation::extractResult(std::string const& name,
 {
     PyObject* xarr = PyDict_GetItemString(m_varsOut, name.c_str());
     if (!xarr)
-        throw pdal::pdal_error("plang output variable '" + name + "' not found.");
+        throw pdal_error("plang output variable '" + name + "' not found.");
     if (!PyArray_Check(xarr))
         throw pdal::pdal_error("Plang output variable  '" + name +
             "' is not a numpy array");
 
     PyArrayObject* arr = (PyArrayObject*)xarr;
 
     npy_intp one = 0;
-    const int pyDataType = pdal::plang::Environment::getPythonDataType(t);
+    const int pyDataType = Environment::pythonType(t);
     PyArray_Descr *dtype = PyArray_DESCR(arr);
 
     npy_intp nDims = PyArray_NDIM(arr);
@@ -263,45 +264,52 @@ bool Invocation::execute()
     m_scriptArgs = PyTuple_New(numArgs);
 
     if (numArgs > 2)
-        throw pdal::pdal_error("Only two arguments -- ins and outs numpy arrays -- can be passed!");
+        throw pdal_error("Only two arguments -- ins and outs numpy "
+            "arrays -- can be passed!");
 
-    PyTuple_SetItem(m_scriptArgs, 0, m_varsIn);
-    if (numArgs > 1)
+    if (m_directArray)
     {
-        Py_INCREF(m_varsOut);
-        PyTuple_SetItem(m_scriptArgs, 1, m_varsOut);
+        PyTuple_SetItem(m_scriptArgs, 0, m_directArray);
+        if (numArgs > 1)
+        {
+            Py_INCREF(m_varsOut);
+            PyTuple_SetItem(m_scriptArgs, 1, m_varsOut);
+        }
+    }
+    else
+    {
+        PyTuple_SetItem(m_scriptArgs, 0, m_varsIn);
+        if (numArgs > 1)
+        {
+            Py_INCREF(m_varsOut);
+            PyTuple_SetItem(m_scriptArgs, 1, m_varsOut);
+        }
     }
-
-    int success(0);
 
     if (m_metadata_PyObject)
     {
-        success = PyModule_AddObject(m_module, "metadata", m_metadata_PyObject);
-        if (success)
+        if (PyModule_AddObject(m_module, "metadata", m_metadata_PyObject))
             throw pdal::pdal_error("unable to set metadata global");
         Py_INCREF(m_metadata_PyObject);
     }
 
     if (m_schema_PyObject)
     {
-        success = PyModule_AddObject(m_module, "schema", m_schema_PyObject);
-        if (success)
+        if (PyModule_AddObject(m_module, "schema", m_schema_PyObject))
             throw pdal::pdal_error("unable to set schema global");
         Py_INCREF(m_srs_PyObject);
     }
 
     if (m_srs_PyObject)
     {
-        success = PyModule_AddObject(m_module, "spatialreference", m_srs_PyObject);
-        if (success)
+        if (PyModule_AddObject(m_module, "spatialreference", m_srs_PyObject))
             throw pdal::pdal_error("unable to set spatialreference global");
         Py_INCREF(m_schema_PyObject);
     }
 
     if (m_pdalargs_PyObject)
     {
-        success = PyModule_AddObject(m_module, "pdalargs", m_pdalargs_PyObject);
-        if (success)
+        if (PyModule_AddObject(m_module, "pdalargs", m_pdalargs_PyObject))
             throw pdal::pdal_error("unable to set pdalargs global");
         Py_INCREF(m_pdalargs_PyObject);
     }
@@ -321,6 +329,7 @@ bool Invocation::execute()
     return (m_scriptResult == Py_True);
 }
 
+
 PyObject* getPyJSON(std::string const& str)
 {
 
@@ -358,6 +367,77 @@ void Invocation::setKWargs(std::string const& s)
     m_pdalargs_PyObject = getPyJSON(s);
 }
 
+
+void Invocation::createArray(char *data, point_count_t count,
+    PointLayoutPtr layout)
+{
+    DimTypeList dims = layout->dimTypes();
+
+    PyObject *fields = PyDict_New();
+    PyObject *nameslist = PyTuple_New(dims.size());
+    int i = 0;
+    for (auto& dim : dims)
+    {
+        PyObject *tup = PyTuple_New(2);
+        PyObject *name = PyUnicode_FromString(layout->dimName(dim.m_id).data());
+        PyArray_Descr *sub =
+            PyArray_DescrFromType(Environment::pythonType(dim.m_type));
+
+        PyTuple_SET_ITEM(tup, 0, (PyObject *)sub);
+        PyTuple_SET_ITEM(tup, 1, PyLong_FromLong(layout->pointSize()));
+
+        PyDict_SetItem(fields, name, tup);
+        PyTuple_SET_ITEM(nameslist, i++, name);
+    }
+
+    PyArray_Descr *descr = PyArray_DescrFromType(NPY_VOID);
+    npy_intp pointSize = layout->pointSize();
+    descr->fields = fields;
+    descr->names = nameslist;
+    descr->elsize = pointSize;
+    descr->flags = NPY_FROM_FIELDS;
+
+    npy_intp longCount(count);
+    PyObject *arr = PyArray_NewFromDescr(&PyArray_Type, descr, 1,
+        &longCount, &pointSize, data, 0, nullptr);
+
+    // Now wrap the created array in a MaskedArray.
+    PyObject *module = PyImport_ImportModule("numpy.ma");
+    PyObject *dict = PyModule_GetDict(module);
+    PyObject *keys = PyDict_Keys(dict);
+    PyObject *o = PyDict_GetItemString(dict, "MaskedArray");
+
+    PyObject *args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, arr);
+
+    m_directArray = PyObject_CallObject(o, args);
+    if (!m_directArray)
+        std::cerr << "Couldn't call object!\n";
+}
+
+
+// IMPORTANT - This only works with data that's contiguous, which is what's
+//  provided with ContiguousPointTable.
+void Invocation::setMask(PythonPointViewPtr view)
+{
+    PyArrayObject *mask =
+        (PyArrayObject *)PyObject_GetAttrString(m_directArray, "mask");
+    size_t bytes = PyArray_NBYTES(mask);
+    size_t stride = PyArray_STRIDE(mask, 0);
+    char *data = (char *)PyArray_DATA(mask);
+
+    // First fill everything with 1's (invalid)
+    // Then fill the fields for each dim in points in the view to 0 to
+    // indicate valid.
+    std::fill(data, data + bytes, 1);
+    for (PointId id = 0; id < view->size(); ++id)
+    {
+        char *pt = data + (view->index(id) * stride);
+        std::fill(pt, pt + stride, 0);
+    }
+}
+
+
 void Invocation::begin(PointView& view, MetadataNode m)
 {
     PointLayoutPtr layout(view.m_pointTable.layout());

plugins/python/plang/Invocation.hpp

@@ -40,7 +40,7 @@
 #include "Environment.hpp"
 
 #include <pdal/Dimension.hpp>
-#include <pdal/PointView.hpp>
+#include "../PythonPointView.hpp"
 
 namespace pdal
 {
@@ -59,6 +59,8 @@ class PDAL_DLL Invocation
     void resetArguments();
 
 
+    void createArray(char *data, point_count_t count, PointLayoutPtr layout);
+    void setMask(PythonPointViewPtr view);
     // creates a Python variable pointing to a (one dimensional) C array
     // adds the new variable to the arguments dictionary
     void insertArgument(std::string const& name,
@@ -97,6 +99,7 @@ class PDAL_DLL Invocation
     PyObject* m_module;
     PyObject* m_dictionary;
     PyObject* m_function;
+    PyObject *m_directArray;
 
     PyObject* m_varsIn;
     PyObject* m_varsOut;