Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into grid-1.4.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@hobu
hobu committed Jun 25, 2015
2 parents 05b861e + 6c1031f commit 9842033
Show file tree
Hide file tree
Showing 11 changed files with 580 additions and 214 deletions.
300 changes: 256 additions & 44 deletions include/pdal/KDIndex.hpp
View file
Expand Up @@ -34,14 +34,15 @@

#pragma once

#include "nanoflann.hpp"

#include <memory>
#include <pdal/PointView.hpp>

namespace nanoflann
{
template<typename Distance, class DatasetAdaptor, int DIM,
typename IndexType>
class KDTreeSingleIndexAdaptor;
typename IndexType> class KDTreeSingleIndexAdaptor;

template<class T, class DataSource, typename _DistanceType>
struct L2_Adaptor;
Expand All @@ -50,71 +51,282 @@ namespace nanoflann
namespace pdal
{

template<int DIM>
class PDAL_DLL KDIndex
{
public:
KDIndex(const PointView& buf);
~KDIndex();
protected:
KDIndex(const PointView& buf) : m_buf(buf)
{}

~KDIndex()
{}

std::size_t kdtree_get_point_count() const;
public:
std::size_t kdtree_get_point_count() const
{ return m_buf.size(); }

double kdtree_get_pt(const PointId idx, int dim) const;
double kdtree_distance(const double *p1, const PointId p2_idx,
size_t /*numDims*/) const;
template <class BBOX> bool kdtree_get_bbox(BBOX& bb) const;
void build()
{
m_index.reset(new my_kd_tree_t(DIM, *this,
nanoflann::KDTreeSingleIndexAdaptorParams(10, DIM)));
m_index->buildIndex();
}

double kdtree_distance(
const PointId idx_p2,
point_count_t size) const;
protected:
const PointView& m_buf;

template <class BBOX> bool kdtree_get_bbox(BBOX &bb) const
typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<
double, KDIndex, double>, KDIndex, -1, std::size_t> my_kd_tree_t;

std::unique_ptr<my_kd_tree_t> m_index;

private:
KDIndex(const KDIndex&);
KDIndex& operator=(KDIndex&);
};

class PDAL_DLL KD2Index : public KDIndex<2>
{
public:
KD2Index(const PointView& buf) : KDIndex<2>(buf)
{
BOX3D bounds = m_buf.calculateBounds();
if (bounds.empty())
return false;
if (!buf.hasDim(Dimension::Id::X))
throw pdal_error("KD2Index: point view missing 'X' dimension.");
if (!buf.hasDim(Dimension::Id::Y))
throw pdal_error("KD2Index: point view missing 'Y' dimension.");
}

bb[0].low = bounds.minx;
bb[0].high = bounds.maxx;
bb[1].low = bounds.miny;
bb[1].high = bounds.maxy;
PointId neighbor(double x, double y)
{
std::vector<PointId> ids = neighbors(x, y, 1);
return (ids.size() ? ids[0] : 0);
}

if (m_3d)
{
bb[2].low = bounds.minx;
bb[2].high = bounds.maxx;
}
return true;
std::vector<PointId> neighbors(double x, double y, point_count_t k)
{
k = std::min(m_buf.size(), k);
std::vector<PointId> output(k);
std::vector<double> out_dist_sqr(k);
nanoflann::KNNResultSet<double, PointId, point_count_t> resultSet(k);

resultSet.init(&output[0], &out_dist_sqr[0]);

std::vector<double> pt;
pt.push_back(x);
pt.push_back(y);
m_index->findNeighbors(resultSet, &pt[0], nanoflann::SearchParams(10));
return output;
}

std::vector<PointId> radius(double const& x, double const& y,
double const& r) const
{
std::vector<PointId> output;
std::vector<std::pair<std::size_t, double>> ret_matches;
nanoflann::SearchParams params;
params.sorted = true;

std::vector<double> pt;
pt.push_back(x);
pt.push_back(y);

// Our distance metric is square distance, so we use the square of
// the radius.
const std::size_t count =
m_index->radiusSearch(&pt[0], r * r, ret_matches, params);

for (std::size_t i = 0; i < count; ++i)
output.push_back(ret_matches[i].first);
return output;
}
};

std::vector<std::size_t> radius(
double const& x,
double const& y,
double const& z,
double const& r) const;
class PDAL_DLL KD3Index : public KDIndex<3>
{
public:
KD3Index(const PointView& buf) : KDIndex<3>(buf)
{
if (!buf.hasDim(Dimension::Id::X))
throw pdal_error("KD3Index: point view missing 'X' dimension.");
if (!buf.hasDim(Dimension::Id::Y))
throw pdal_error("KD3Index: point view missing 'Y' dimension.");
if (!buf.hasDim(Dimension::Id::Z))
throw pdal_error("KD3Index: point view missing 'Z' dimension.");
}

PointId neighbor(double const& x, double const& y, double const& z) const
PointId neighbor(double x, double y, double z)
{
std::vector<PointId> ids = neighbors(x, y, z, 1);
return ids[0];
return (ids.size() ? ids[0] : 0);
}

std::vector<PointId> neighbors(
double const& x,
double const& y,
double const& z,
point_count_t count) const;
std::vector<PointId> neighbors(double x, double y, double z,
point_count_t k)
{
k = std::min(m_buf.size(), k);
std::vector<PointId> output(k);
std::vector<double> out_dist_sqr(k);
nanoflann::KNNResultSet<double, PointId, point_count_t> resultSet(k);

void build(bool b3d = true);
resultSet.init(&output[0], &out_dist_sqr[0]);

private:
const PointView& m_buf;
bool m_3d;
std::vector<double> pt;
pt.push_back(x);
pt.push_back(y);
pt.push_back(z);
m_index->findNeighbors(resultSet, &pt[0], nanoflann::SearchParams(10));
return output;
}

typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Adaptor<
double, KDIndex, double>, KDIndex, -1, std::size_t> my_kd_tree_t;
std::vector<PointId> radius(double x, double y, double z, double r) const
{
std::vector<PointId> output;
std::vector<std::pair<std::size_t, double>> ret_matches;
nanoflann::SearchParams params;
params.sorted = true;

std::unique_ptr<my_kd_tree_t> m_index;
std::vector<double> pt;
pt.push_back(x);
pt.push_back(y);
pt.push_back(z);

KDIndex(const KDIndex&);
KDIndex& operator=(KDIndex&);
// Our distance metric is square distance, so we use the square of
// the radius.
const std::size_t count =
m_index->radiusSearch(&pt[0], r * r, ret_matches, params);

for (std::size_t i = 0; i < count; ++i)
output.push_back(ret_matches[i].first);
return output;
}
};

template<>
inline
double KDIndex<2>::kdtree_get_pt(const PointId idx, int dim) const
{
if (idx >= m_buf.size())
return 0.0;

Dimension::Id::Enum id = Dimension::Id::Unknown;
switch (dim)
{
case 0:
id = Dimension::Id::X;
break;
case 1:
id = Dimension::Id::Y;
break;
default:
throw pdal_error("kdtree_get_pt: Request for invalid dimension "
"from nanoflann");
}
return m_buf.getFieldAs<double>(id, idx);
}

template<>
inline
double KDIndex<3>::kdtree_get_pt(const PointId idx, int dim) const
{
if (idx >= m_buf.size())
return 0.0;

Dimension::Id::Enum id = Dimension::Id::Unknown;
switch (dim)
{
case 0:
id = Dimension::Id::X;
break;
case 1:
id = Dimension::Id::Y;
break;
case 2:
id = Dimension::Id::Z;
break;
default:
throw pdal_error("kdtree_get_pt: Request for invalid dimension "
"from nanoflann");
}
return m_buf.getFieldAs<double>(id, idx);
}

// nanoflann hands us a vector that represents the position of p1. We fetch
// the position of p2 and and compute the square distance.
template<>
inline double KDIndex<2>::kdtree_distance(const double *p1, const PointId idx,
size_t /*numDims*/) const
{
double d0 = p1[0] - m_buf.getFieldAs<double>(Dimension::Id::X, idx);
double d1 = p1[1] - m_buf.getFieldAs<double>(Dimension::Id::Y, idx);

return (d0 * d0 + d1 * d1);
}

template<>
inline double KDIndex<3>::kdtree_distance(const double *p1, const PointId idx,
size_t /*numDims*/) const
{
double d0 = p1[0] - m_buf.getFieldAs<double>(Dimension::Id::X, idx);
double d1 = p1[1] - m_buf.getFieldAs<double>(Dimension::Id::Y, idx);
double d2 = p1[2] - m_buf.getFieldAs<double>(Dimension::Id::Z, idx);

return (d0 * d0 + d1 * d1 + d2 * d2);
}


template<>
template <class BBOX>
bool KDIndex<2>::kdtree_get_bbox(BBOX& bb) const
{
if (m_buf.empty())
{
bb[0].low = 0.0;
bb[0].high = 0.0;
bb[1].low = 0.0;
bb[1].high = 0.0;
}
else
{
BOX3D bounds = m_buf.calculateBounds();

bb[0].low = bounds.minx;
bb[0].high = bounds.maxx;
bb[1].low = bounds.miny;
bb[1].high = bounds.maxy;
}
return true;
}

template<>
template <class BBOX>
bool KDIndex<3>::kdtree_get_bbox(BBOX& bb) const
{
if (m_buf.empty())
{
bb[0].low = 0.0;
bb[0].high = 0.0;
bb[1].low = 0.0;
bb[1].high = 0.0;
bb[2].low = 0.0;
bb[2].high = 0.0;
}
else
{
BOX3D bounds = m_buf.calculateBounds();

bb[0].low = bounds.minx;
bb[0].high = bounds.maxx;
bb[1].low = bounds.miny;
bb[1].high = bounds.maxy;
bb[2].low = bounds.minz;
bb[2].high = bounds.maxz;
}
return true;
}

} // namespace pdal

8 changes: 7 additions & 1 deletion io/ply/PlyReader.cpp
View file
Expand Up @@ -90,6 +90,7 @@ int readPlyCallback(p_ply_argument argument)
long numToRead;
p_ply_property property;
const char * propertyName;

if (!ply_get_argument_element(argument, &element, &index))
{
std::stringstream ss;
Expand Down Expand Up @@ -229,9 +230,14 @@ point_count_t PlyReader::read(PointViewPtr view, point_count_t num)
context.view = view;
context.dimensionMap = m_vertexDimensions;

// It's possible that point_count_t holds a value that's larger than the
// long that is the maximum rply (don't know about ply) point count.
long cnt;
cnt = Utils::inRange<long>(num) ? num : (std::numeric_limits<long>::max)();
for (auto it : m_vertexDimensions)
{
ply_set_read_cb(m_ply, "vertex", it.first.c_str(), readPlyCallback, &context, num);
ply_set_read_cb(m_ply, "vertex", it.first.c_str(), readPlyCallback,
&context, cnt);
}
if (!ply_read(m_ply))
{
Expand Down
8 changes: 4 additions & 4 deletions kernels/delta/DeltaKernel.cpp
View file
Expand Up @@ -152,8 +152,8 @@ int DeltaKernel::execute()
}

// Index the candidate data.
KDIndex index(*candView);
index.build(m_3d);
KD3Index index(*candView);
index.build();

MetadataNode root;

Expand All @@ -168,7 +168,7 @@ int DeltaKernel::execute()


MetadataNode DeltaKernel::dump(PointViewPtr& srcView, PointViewPtr& candView,
KDIndex& index, DimIndexMap& dims)
KD3Index& index, DimIndexMap& dims)
{
MetadataNode root;

Expand Down Expand Up @@ -217,7 +217,7 @@ void DeltaKernel::accumulate(DimIndex& d, double v)


MetadataNode DeltaKernel::dumpDetail(PointViewPtr& srcView,
PointViewPtr& candView, KDIndex& index, DimIndexMap& dims)
PointViewPtr& candView, KD3Index& index, DimIndexMap& dims)
{
MetadataNode root;

Expand Down

0 comments on commit 9842033

Please sign in to comment.