Skip to content

Commit fea94f6

Browse files
runetterunette
authored
Add Meshes to the pdal API (#82)
* first commit * pymesh typos * pymesh typos * pymesh bugs and added python defs * pypipeline bugs * Update PyPipeline.hpp * testing * testing * fixing the zero data bug * Update PyMesh.cpp * Update PyMesh.cpp * Update PyMesh.cpp * Update PyMesh.cpp * Update PyMesh.cpp * Update PyMesh.cpp * Added Meshio integration * meshio bug fixes * Delete .DS_Store * Delete .DS_Store * Update .gitignore * Delete .DS_Store * tests * test fix * Updates Docs * typos * typos * Update README.rst
1 parent 763d5c3 commit fea94f6

File tree

12 files changed

+535
-7
lines changed

12 files changed

+535
-7
lines changed

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@ dist/*
99
*.o
1010
*.so
1111
*.dylib
12+
.DS_Store

‎PKG-INFO‎

Lines changed: 37 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -149,16 +149,48 @@ Description: ===================================================================
149149
count = p.execute()
150150
print (count)
151151

152+
Accessing Mesh Data
153+
................................................................................
154+
155+
Some PDAL stages (for instance ``filters.delaunay``) create TIN type mesh data.
156+
157+
This data can be accessed in Python using the ``Pipeline.meshes`` property, which returns ``numpy.ndarray``
158+
of shape (1,n) where n is the number of Triangles in the mesh. If the PointView contains no mesh data, then n = 0.
159+
160+
Each Triangle is a tuple ``(A,B,C)`` where A, B and C are indicees into the PointView for the point that is the vertex for the Triangle.
161+
162+
Meshio Integration
163+
................................................................................
164+
165+
The meshes property provides the face data but is not easy to use as a mesh. Therefore, we have provided optional Integration
166+
into the `Meshio <https://github.com/nschloe/meshio>`__ library.
167+
168+
The ``pdal.Pipeline`` class provides the ``get_meshio(idx: int) -> meshio.Mesh`` method. This
169+
method creates a `Mesh` object from the `PointView` array and mesh properties.
170+
171+
.. note:: The meshio integration requires that meshio is installed (e.g. ``pip install meshio``). If it
172+
is not, then the method fails with an informative RuntimeError.
173+
174+
Simple use of the functionality could be as follows:
175+
176+
-- code-block:: python
177+
import pdal
178+
179+
...
180+
pl = pdal.Pipeline(pipeline)
181+
pl.execute()
152182

183+
mesh = pl.get_meshio(0)
184+
mesh.write('test.obj')
153185

154186

155-
.. _`Numpy`: http://www.numpy.org/
156-
.. _`schema`: http://www.pdal.io/dimensions.html
157-
.. _`metadata`: http://www.pdal.io/development/metadata.html
187+
.. _`Numpy`: http://www.numpy.org/
188+
.. _`schema`: http://www.pdal.io/dimensions.html
189+
.. _`metadata`: http://www.pdal.io/development/metadata.html
158190

159191

160-
.. image:: https://ci.appveyor.com/api/projects/status/of4kecyahpo8892d
161-
:target: https://ci.appveyor.com/project/hobu/python/
192+
.. image:: https://ci.appveyor.com/api/projects/status/of4kecyahpo8892d
193+
:target: https://ci.appveyor.com/project/hobu/python/
162194

163195
Requirements
164196
================================================================================

‎README.rst‎

Lines changed: 79 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ sorts it by the ``X`` dimension:
6262
metadata = pipeline.metadata
6363
log = pipeline.log
6464
65-
Reading using Numpy Arrays
65+
Reading using Numpy Arrays
6666
................................................................................
6767

6868
The following more complex scenario demonstrates the full cycling between
@@ -148,6 +148,84 @@ PDAL and Python:
148148
print (count)
149149
150150
151+
Accessing Mesh Data
152+
................................................................................
153+
154+
Some PDAL stages (for instance ``filters.delaunay``) create TIN type mesh data.
155+
156+
This data can be accessed in Python using the ``Pipeline.meshes`` property, which returns a ``numpy.ndarray``
157+
of shape (1,n) where n is the number of Triangles in the mesh.
158+
159+
If the PointView contains no mesh data, then n = 0.
160+
161+
Each Triangle is a tuple ``(A,B,C)`` where A, B and C are indices into the PointView identifying the point that is the vertex for the Triangle.
162+
163+
Meshio Integration
164+
................................................................................
165+
166+
The meshes property provides the face data but is not easy to use as a mesh. Therefore, we have provided optional Integration
167+
into the `Meshio <https://github.com/nschloe/meshio>`__ library.
168+
169+
The ``pdal.Pipeline`` class provides the ``get_meshio(idx: int) -> meshio.Mesh`` method. This
170+
method creates a `Mesh` object from the `PointView` array and mesh properties.
171+
172+
.. note:: The meshio integration requires that meshio is installed (e.g. ``pip install meshio``). If it is not, then the method fails with an informative RuntimeError.
173+
174+
Simple use of the functionality could be as follows:
175+
176+
.. code-block:: python
177+
178+
import pdal
179+
180+
...
181+
pl = pdal.Pipeline(pipeline)
182+
pl.execute()
183+
184+
mesh = pl.get_meshio(0)
185+
mesh.write('test.obj')
186+
187+
Advanced Mesh Use Case
188+
................................................................................
189+
190+
USE-CASE : Take a LiDAR map, create a mesh from the ground points, split into tiles and store the tiles in PostGIS.
191+
192+
.. note:: Like ``Pipeline.arrays``, ``Pipeline.meshes`` returns a list of ``numpy.ndarray`` to provide for the case where the output from a Pipeline is multiple PointViews
193+
194+
(example using 1.2-with-color.las and not doing the ground classification for clarity)
195+
196+
.. code-block:: python
197+
198+
import pdal
199+
import json
200+
import psycopg2
201+
import io
202+
203+
pipe = [
204+
'.../python/test/data/1.2-with-color.las',
205+
{"type": "filters.splitter", "length": 1000},
206+
{"type": "filters.delaunay"}
207+
]
208+
209+
pl = pdal.Pipeline(json.dumps(pipe))
210+
pl.execute()
211+
212+
conn = psycopg(%CONNNECTION_STRING%)
213+
buffer = io.StringIO
214+
215+
for idx in range(len(pl.meshes)):
216+
m = pl.get_meshio(idx)
217+
if m:
218+
m.write(buffer, file_format = "wkt")
219+
with conn.cursor() as curr:
220+
curr.execute(
221+
"INSERT INTO %table-name% (mesh) VALUES (ST_GeomFromEWKT(%(ewkt)s)",
222+
{ "ewkt": buffer.getvalue()}
223+
)
224+
225+
conn.commit()
226+
conn.close()
227+
buffer.close()
228+
151229
152230
153231
.. _`Numpy`: http://www.numpy.org/

‎pdal/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
set(EXTENSION_SRC
22
PyArray.cpp
33
PyArray.hpp
4+
PyMesh.cpp
5+
PyMesh.hpp
46
PyDimension.hpp
57
PyPipeline.cpp
68
PyPipeline.hpp)

‎pdal/PyMesh.cpp‎

Lines changed: 216 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,216 @@
1+
/******************************************************************************
2+
* Copyright (c) 2021, Runette Software Ltd (www.runette.co.uk)
3+
*
4+
* All rights reserved.
5+
*
6+
* Redistribution and use in source and binary forms, with or without
7+
* modification, are permitted provided that the following
8+
* conditions are met:
9+
*
10+
* * Redistributions of source code must retain the above copyright
11+
* notice, this list of conditions and the following disclaimer.
12+
* * Redistributions in binary form must reproduce the above copyright
13+
* notice, this list of conditions and the following disclaimer in
14+
* the documentation and/or other materials provided
15+
* with the distribution.
16+
* * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
17+
* names of its contributors may be used to endorse or promote
18+
* products derived from this software without specific prior
19+
* written permission.
20+
*
21+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22+
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23+
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24+
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25+
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26+
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27+
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
28+
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
29+
* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30+
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31+
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32+
* OF SUCH DAMAGE.
33+
****************************************************************************/
34+
35+
#include "PyMesh.hpp"
36+
37+
#include <numpy/arrayobject.h>
38+
39+
40+
namespace pdal
41+
{
42+
namespace python
43+
{
44+
45+
// Create new empty mesh
46+
//
47+
Mesh::Mesh() : m_mesh(nullptr)
48+
{
49+
hasMesh = false;
50+
if (_import_array() < 0)
51+
throw pdal_error("Could not import numpy.core.multiarray.");
52+
}
53+
54+
Mesh::~Mesh()
55+
{
56+
if (m_mesh)
57+
Py_XDECREF((PyObject *)m_mesh);
58+
}
59+
60+
// Update from a PointView
61+
//
62+
void Mesh::update(PointViewPtr view)
63+
{
64+
npy_intp size;
65+
66+
if (m_mesh)
67+
Py_XDECREF((PyObject *)m_mesh);
68+
m_mesh = nullptr; // Just in case of an exception.
69+
70+
TriangularMesh* mesh = view->mesh();
71+
if (mesh)
72+
{
73+
hasMesh = true;
74+
size = mesh->size();
75+
} else {
76+
hasMesh = false;
77+
size = 0;
78+
}
79+
80+
PyObject *dtype_dict = (PyObject*)buildNumpyDescription(view);
81+
if (!dtype_dict)
82+
throw pdal_error("Unable to build numpy dtype "
83+
"description dictionary");
84+
85+
PyArray_Descr *dtype = nullptr;
86+
if (PyArray_DescrConverter(dtype_dict, &dtype) == NPY_FAIL)
87+
throw pdal_error("Unable to build numpy dtype");
88+
Py_XDECREF(dtype_dict);
89+
90+
// This is a 1 x size array.
91+
m_mesh = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype,
92+
1, &size, 0, nullptr, NPY_ARRAY_CARRAY, nullptr);
93+
94+
for (PointId idx = 0; idx < size; idx++)
95+
{
96+
char* p = (char *)PyArray_GETPTR1(m_mesh, idx);
97+
const Triangle& t = (*mesh)[idx];
98+
uint32_t a = (uint32_t)t.m_a;
99+
std::memcpy(p, &a, 4);
100+
uint32_t b = (uint32_t)t.m_b;
101+
std::memcpy(p + 4, &b, 4);
102+
uint32_t c = (uint32_t)t.m_c;
103+
std::memcpy(p + 8, &c, 4);
104+
}
105+
}
106+
107+
108+
109+
PyArrayObject *Mesh::getPythonArray() const
110+
{
111+
return m_mesh;
112+
}
113+
114+
PyObject* Mesh::buildNumpyDescription(PointViewPtr view) const
115+
{
116+
// Build up a numpy dtype dictionary
117+
//
118+
// {'formats': ['f8', 'f8', 'f8', 'u2', 'u1', 'u1', 'u1', 'u1', 'u1',
119+
// 'f4', 'u1', 'u2', 'f8', 'u2', 'u2', 'u2'],
120+
// 'names': ['X', 'Y', 'Z', 'Intensity', 'ReturnNumber',
121+
// 'NumberOfReturns', 'ScanDirectionFlag', 'EdgeOfFlightLine',
122+
// 'Classification', 'ScanAngleRank', 'UserData',
123+
// 'PointSourceId', 'GpsTime', 'Red', 'Green', 'Blue']}
124+
//
125+
126+
Dimension::IdList dims = view->dims();
127+
128+
PyObject* dict = PyDict_New();
129+
PyObject* formats = PyList_New(3);
130+
PyObject* titles = PyList_New(3);
131+
132+
PyList_SetItem(titles, 0, PyUnicode_FromString("A"));
133+
PyList_SetItem(formats, 0, PyUnicode_FromString("u4"));
134+
PyList_SetItem(titles, 1, PyUnicode_FromString("B"));
135+
PyList_SetItem(formats, 1, PyUnicode_FromString("u4"));
136+
PyList_SetItem(titles, 2, PyUnicode_FromString("C"));
137+
PyList_SetItem(formats, 2, PyUnicode_FromString("u4"));
138+
139+
140+
PyDict_SetItemString(dict, "names", titles);
141+
PyDict_SetItemString(dict, "formats", formats);
142+
143+
return dict;
144+
}
145+
146+
bool Mesh::rowMajor() const
147+
{
148+
return m_rowMajor;
149+
}
150+
151+
Mesh::Shape Mesh::shape() const
152+
{
153+
return m_shape;
154+
}
155+
156+
157+
MeshIter& Mesh::iterator()
158+
{
159+
MeshIter *it = new MeshIter(*this);
160+
m_iterators.push_back(std::unique_ptr<MeshIter>(it));
161+
return *it;
162+
}
163+
164+
MeshIter::MeshIter(Mesh& mesh)
165+
{
166+
m_iter = NpyIter_New(mesh.getPythonArray(),
167+
NPY_ITER_EXTERNAL_LOOP | NPY_ITER_READONLY | NPY_ITER_REFS_OK,
168+
NPY_KEEPORDER, NPY_NO_CASTING, NULL);
169+
if (!m_iter)
170+
throw pdal_error("Unable to create numpy iterator.");
171+
172+
char *itererr;
173+
m_iterNext = NpyIter_GetIterNext(m_iter, &itererr);
174+
if (!m_iterNext)
175+
{
176+
NpyIter_Deallocate(m_iter);
177+
throw pdal_error(std::string("Unable to create numpy iterator: ") +
178+
itererr);
179+
}
180+
m_data = NpyIter_GetDataPtrArray(m_iter);
181+
m_stride = NpyIter_GetInnerStrideArray(m_iter);
182+
m_size = NpyIter_GetInnerLoopSizePtr(m_iter);
183+
m_done = false;
184+
}
185+
186+
MeshIter::~MeshIter()
187+
{
188+
NpyIter_Deallocate(m_iter);
189+
}
190+
191+
MeshIter& MeshIter::operator++()
192+
{
193+
if (m_done)
194+
return *this;
195+
196+
if (--(*m_size))
197+
*m_data += *m_stride;
198+
else if (!m_iterNext(m_iter))
199+
m_done = true;
200+
return *this;
201+
}
202+
203+
MeshIter::operator bool () const
204+
{
205+
return !m_done;
206+
}
207+
208+
char * MeshIter::operator * () const
209+
{
210+
return *m_data;
211+
}
212+
213+
214+
} // namespace python
215+
} // namespace pdal
216+

0 commit comments

Comments
 (0)