diff --git a/CMakeLists.txt b/CMakeLists.txt index c5c6610419..8a7d1fa595 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -287,6 +287,13 @@ if(WITH_GEOTIFF) endif() find_package(Flann 1.7) +if (FLANN_FOUND) + mark_as_advanced(CLEAR FLANN_INCLUDE_DIR) + mark_as_advanced(CLEAR FLANN_LIBRARY) + include_directories(${FLANN_INCLUDE_DIR}) + set(PDAL_HAVE_FLANN 1) +endif() + if(WITH_GDAL) if (NOT GEOTIFF_FOUND) diff --git a/include/pdal/filters/Index.hpp b/include/pdal/filters/Index.hpp index 4df8ac6dce..ed431981d0 100644 --- a/include/pdal/filters/Index.hpp +++ b/include/pdal/filters/Index.hpp @@ -38,8 +38,11 @@ #include #include -#include -#include +#include + +#ifdef PDAL_HAVE_FLANN +#include +#endif namespace pdal @@ -88,15 +91,37 @@ class PDAL_DLL Index : public pdal::FilterSequentialIterator { public: Index(const pdal::filters::Index& filter, PointBuffer& buffer); - + + std::vector query(double const& x, double const& y, double const& z, double distance, boost::uint32_t count=1); + protected: virtual void readBufferBeginImpl(PointBuffer&); + virtual void readEndImpl(); private: boost::uint64_t skipImpl(boost::uint64_t); boost::uint32_t readBufferImpl(PointBuffer&); bool atEndImpl() const; + + double getScaledValue(PointBuffer& data, + Dimension const& d, + std::size_t pointIndex) const; const pdal::filters::Index& m_stage; + + boost::scoped_array m_query_data; + boost::scoped_array m_distance_data; + boost::scoped_array m_indice_data; + + flann::Index >* m_index; + flann::Matrix* m_dataset; + flann::Matrix* m_indices; + flann::Matrix* m_query; + flann::Matrix* m_distances; + + pdal::Dimension const* m_xDim; + pdal::Dimension const* m_yDim; + pdal::Dimension const* m_zDim; + }; diff --git a/pdal_defines.h.in b/pdal_defines.h.in index a9f4e70281..7df67bad53 100644 --- a/pdal_defines.h.in +++ b/pdal_defines.h.in @@ -24,17 +24,17 @@ /* * availability of 3rd-party libraries */ -#cmakedefine PDAL_HAVE_LIBLAS -#cmakedefine PDAL_HAVE_LASZIP -#cmakedefine PDAL_HAVE_ORACLE -#cmakedefine PDAL_HAVE_GDAL #cmakedefine PDAL_HAVE_GEOS +#cmakedefine PDAL_HAVE_GDAL +#cmakedefine PDAL_HAVE_FLANN +#cmakedefine PDAL_HAVE_LASZIP +#cmakedefine PDAL_HAVE_LIBLAS #cmakedefine PDAL_HAVE_LIBXML2 #cmakedefine PDAL_HAVE_LIBGEOTIFF #cmakedefine PDAL_HAVE_MRSID +#cmakedefine PDAL_HAVE_ORACLE #cmakedefine PDAL_HAVE_P2G #cmakedefine PDAL_HAVE_PYTHON -#cmakedefine PDAL_EMBED_BOOST /* @@ -56,6 +56,12 @@ * platform compiler */ #cmakedefine PDAL_COMPILER_MSVC + +/* + * Use the embedded boost tree instead of system one + */ +#cmakedefine PDAL_EMBED_BOOST + #cmakedefine PDAL_COMPILER_VC10 #cmakedefine PDAL_COMPILER_VC9 #cmakedefine PDAL_COMPILER_VC8 diff --git a/src/StageFactory.cpp b/src/StageFactory.cpp index 445437e1dd..43248d5d8b 100644 --- a/src/StageFactory.cpp +++ b/src/StageFactory.cpp @@ -73,6 +73,7 @@ #include #include #include +#include #include #include @@ -154,6 +155,7 @@ MAKE_FILTER_CREATOR(Color, pdal::filters::Color) MAKE_FILTER_CREATOR(Colorization, pdal::filters::Colorization) MAKE_FILTER_CREATOR(Crop, pdal::filters::Crop) MAKE_FILTER_CREATOR(Decimation, pdal::filters::Decimation) +MAKE_FILTER_CREATOR(Index, pdal::filters::Index) MAKE_FILTER_CREATOR(InPlaceReprojection, pdal::filters::InPlaceReprojection) #ifdef PDAL_HAVE_PYTHON @@ -347,6 +349,7 @@ void StageFactory::registerKnownFilters() REGISTER_FILTER(Crop, pdal::filters::Crop); REGISTER_FILTER(Decimation, pdal::filters::Decimation); REGISTER_FILTER(Reprojection, pdal::filters::Reprojection); + REGISTER_FILTER(Index, pdal::filters::Index); REGISTER_FILTER(InPlaceReprojection, pdal::filters::InPlaceReprojection); #ifdef PDAL_HAVE_PYTHON diff --git a/src/filters/Index.cpp b/src/filters/Index.cpp index 21cf47e88d..f60da3239f 100644 --- a/src/filters/Index.cpp +++ b/src/filters/Index.cpp @@ -40,11 +40,6 @@ #include -#ifdef PDAL_HAVE_GDAL -#include -#include -#include -#endif namespace pdal { @@ -73,40 +68,15 @@ const Options Index::getDefaultOptions() const Option x("x_dim", std::string("X"), "Dimension name to use for 'X' data"); Option y("y_dim", std::string("Y"), "Dimension name to use for 'Y' data"); + Option z("z_dim", std::string("Z"), "Dimension name to use for 'Z' data"); + + Option filename("filename", "", "Filename to store the index in"); - pdal::Option red("dimension", "Red", ""); - pdal::Option b0("band",1, ""); - pdal::Option s0("scale", 1.0f, "scale factor for this dimension"); - pdal::Options redO; - redO.add(b0); - redO.add(s0); - red.setOptions(redO); - - pdal::Option green("dimension", "Green", ""); - pdal::Option b1("band",2, ""); - pdal::Option s1("scale", 1.0f, "scale factor for this dimension"); - pdal::Options greenO; - greenO.add(b1); - greenO.add(s1); - green.setOptions(greenO); - - pdal::Option blue("dimension", "Blue", ""); - pdal::Option b2("band",3, ""); - pdal::Option s2("scale", 1.0f, "scale factor for this dimension"); - pdal::Options blueO; - blueO.add(b2); - blueO.add(s2); - blue.setOptions(blueO); - - pdal::Option reproject("reproject", false, "Reproject the input data into the same coordinate system as the raster?"); options.add(x); options.add(y); - options.add(red); - options.add(green); - options.add(blue); - options.add(reproject); - + options.add(z); + options.add(filename); return options; } @@ -134,7 +104,16 @@ namespace sequential Index::Index(const pdal::filters::Index& filter, PointBuffer& buffer) : pdal::FilterSequentialIterator(filter, buffer) , m_stage(filter) + , m_index(0) + , m_dataset(0) + , m_indices(0) + , m_query(0) + , m_distances(0) + , m_xDim(0) + , m_yDim(0) + , m_zDim(0) { + return; } @@ -143,19 +122,156 @@ void Index::readBufferBeginImpl(PointBuffer& buffer) // Cache dimension positions pdal::Schema const& schema = buffer.getSchema(); - - + + std::string x_name = m_stage.getOptions().getValueOrDefault("x_dim", "X"); + std::string y_name = m_stage.getOptions().getValueOrDefault("x_dim", "Y"); + std::string z_name = m_stage.getOptions().getValueOrDefault("x_dim", "Z"); + + m_xDim = &schema.getDimension(x_name); + m_yDim = &schema.getDimension(y_name); + m_zDim = &schema.getDimension(z_name); + + + if (!m_stage.getNumPoints()) + throw pdal_error("Unable to create index from pipeline that has an indeterminate number of points!"); + + + boost::scoped_array data(new float[ m_stage.getNumPoints() *3 ]); + m_query_data.swap(data); + m_dataset = new flann::Matrix(m_query_data.get(), m_stage.getNumPoints(), 3); + } +std::vector Index::query(double const& x, double const& y, double const& z, double distance, boost::uint32_t count) +{ + + boost::scoped_array distances(new float[ count *3 ]); + m_distance_data.swap(distances); + m_distances = new flann::Matrix(m_distance_data.get(), count, 3); + + boost::scoped_array indices(new int[ m_stage.getNumPoints() *3 ]); + m_indice_data.swap(indices); + m_indices = new flann::Matrix(m_indice_data.get(), m_stage.getNumPoints(), 3); + + std::vector output; + + return output; +} boost::uint32_t Index::readBufferImpl(PointBuffer& data) { const boost::uint32_t numRead = getPrevIterator().read(data); + for (boost::uint32_t pointIndex=0; pointIndex(getScaledValue(data, *m_xDim, pointIndex)); + float y = static_cast(getScaledValue(data, *m_yDim, pointIndex)); + float z = static_cast(getScaledValue(data, *m_zDim, pointIndex)); + + m_query_data[pointIndex] = x; + m_query_data[pointIndex+1] = y; + m_query_data[pointIndex+2] = z; + + } + return numRead; } +double Index::getScaledValue(PointBuffer& data, + Dimension const& d, + std::size_t pointIndex) const +{ + double output(0.0); + + float flt(0.0); + boost::int8_t i8(0); + boost::uint8_t u8(0); + boost::int16_t i16(0); + boost::uint16_t u16(0); + boost::int32_t i32(0); + boost::uint32_t u32(0); + boost::int64_t i64(0); + boost::uint64_t u64(0); + + boost::uint32_t size = d.getByteSize(); + switch (d.getInterpretation()) + { + case dimension::Float: + if (size == 4) + { + flt = data.getField(d, pointIndex); + output = static_cast(flt); + } + if (size == 8) + { + output = data.getField(d, pointIndex); + } + break; + + case dimension::SignedInteger: + case dimension::SignedByte: + if (size == 1) + { + i8 = data.getField(d, pointIndex); + output = d.applyScaling(i8); + } + if (size == 2) + { + i16 = data.getField(d, pointIndex); + output = d.applyScaling(i16); + } + if (size == 4) + { + i32 = data.getField(d, pointIndex); + output = d.applyScaling(i32); + } + if (size == 8) + { + i64 = data.getField(d, pointIndex); + output = d.applyScaling(i64); + } + break; + + case dimension::UnsignedInteger: + case dimension::UnsignedByte: + if (size == 1) + { + u8 = data.getField(d, pointIndex); + output = d.applyScaling(u8); + } + if (size == 2) + { + u16 = data.getField(d, pointIndex); + output = d.applyScaling(u16); + } + if (size == 4) + { + u32 = data.getField(d, pointIndex); + output = d.applyScaling(u32); + } + if (size == 8) + { + u64 = data.getField(d, pointIndex); + output = d.applyScaling(u64); + } + break; + + case dimension::Pointer: // stored as 64 bits, even on a 32-bit box + case dimension::Undefined: + throw pdal_error("Dimension data type unable to be reprojected"); + } + + return output; +} +void Index::readEndImpl() +{ + + // Build the index + m_index = new flann::Index >(*m_dataset, flann::KDTreeIndexParams(4)); + m_index->buildIndex(); +} + boost::uint64_t Index::skipImpl(boost::uint64_t count) { diff --git a/test/data/pipeline/pipeline_index.xml b/test/data/pipeline/pipeline_index.xml new file mode 100644 index 0000000000..45273b8c50 --- /dev/null +++ b/test/data/pipeline/pipeline_index.xml @@ -0,0 +1,10 @@ + + + + + + + + diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 241f35faf1..7109dbe92a 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -22,6 +22,7 @@ SET(PDAL_UNITTEST_TEST_SRC filters/ChipperTest.cpp filters/ColorFilterTest.cpp filters/ColorizationFilterTest.cpp + filters/IndexFilterTest.cpp ConfigTest.cpp filters/CropFilterTest.cpp filters/DecimationFilterTest.cpp diff --git a/test/unit/filters/IndexFilterTest.cpp b/test/unit/filters/IndexFilterTest.cpp new file mode 100644 index 0000000000..33bd5b999a --- /dev/null +++ b/test/unit/filters/IndexFilterTest.cpp @@ -0,0 +1,75 @@ +/****************************************************************************** +* Copyright (c) 2011, Michael P. Gerlek (mpg@flaxen.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 + +#include +#include +#include +#include +#include +#include + +#include "Support.hpp" + +using namespace pdal; + +BOOST_AUTO_TEST_SUITE(IndexFilterTest) + +BOOST_AUTO_TEST_CASE(test_1) +{ + pdal::Option option("filename", Support::datapath("pipeline/pipeline_index.xml")); + pdal::Options options(option); + + pdal::drivers::pipeline::Reader reader(options); + reader.initialize(); + + pdal::filters::Index const* filter = static_cast(reader.getManager().getStage()); + pdal::Options opt = filter->getOptions(); + // std::cout << "filter ops: " << opt << std::endl; + + const pdal::Schema& schema = filter->getSchema(); + pdal::PointBuffer data(schema, 1); + + pdal::StageSequentialIterator* iter = filter->createSequentialIterator(data); + + boost::uint32_t numRead = iter->read(data); + BOOST_CHECK(numRead == 1); + +} + + + + +BOOST_AUTO_TEST_SUITE_END()