Skip to content

Commit

Permalink
Use neighbors() instead of knnSearch() to get only valid neighbors. (#…
Browse files Browse the repository at this point in the history 
…3072)

* Use neighbors() instead of knnSearch() to get only valid neighbors.
Compute filled vector in place.

* Merge fix.
  • Loading branch information
@abellgithub
abellgithub committed May 26, 2020
1 parent 9e00fde commit ea30f6b
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 25 deletions.
62 changes: 41 additions & 21 deletions filters/SMRFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,11 @@ void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
gsurfs = (gx.cwiseProduct(gx) + gy.cwiseProduct(gy)).cwiseSqrt();
std::vector<double> gsurfsV(gsurfs.data(),
gsurfs.data() + gsurfs.size());
std::vector<double> gsurfs_fillV = knnfill(view, gsurfsV);

//ABELL - We can eliminate this copy if we're OK with not writing
// both the filled and non-filled array to output.
std::vector<double> gsurfs_fillV = gsurfsV;
knnfill(view, gsurfs_fillV);
gsurfs = Map<MatrixXd>(gsurfs_fillV.data(), m_rows, m_cols);
thresh =
(m_args->m_threshold + m_args->m_scalar * gsurfs.array()).matrix();
Expand Down Expand Up @@ -349,6 +353,8 @@ void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
int c = static_cast<int>(floor((x - m_bounds.minx) / m_args->m_cell));
int r = static_cast<int>(floor((y - m_bounds.miny) / m_args->m_cell));

size_t cell = c * m_rows + r;

// TODO(chambbj): We don't quite do this by the book and yet it seems to
// work reasonably well:
// "The calculation requires that both elevation and slope are
Expand All @@ -358,7 +364,7 @@ void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
// DEM nearly corresponds to the resolution of the LIDAR data. Based on
// these results, we find that a splined cubic interpolation provides
// the best results."
if (std::isnan(ZIpro[c * m_rows + r]))
if (std::isnan(ZIpro[cell]))
continue;

if (std::isnan(gsurfs(r, c)))
Expand All @@ -368,7 +374,7 @@ void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
// ground/object LIDAR points. This is accomplished by measuring the
// vertical distance between each LIDAR point and the provisional
// DEM, and applying a threshold calculation."
if (std::fabs(ZIpro[c * m_rows + r] - z) > thresh(r, c))
if (std::fabs(ZIpro[cell] - z) > thresh(r, c))
view->setField(Id::Classification, i, ClassLabel::Unclassified);
else
view->setField(Id::Classification, i, ClassLabel::Ground);
Expand Down Expand Up @@ -470,15 +476,20 @@ std::vector<double> SMRFilter::createZImin(PointViewPtr view)
int c = static_cast<int>(floor((x - m_bounds.minx) / m_args->m_cell));
int r = static_cast<int>(floor((y - m_bounds.miny) / m_args->m_cell));

if (z < ZIminV[c * m_rows + r] || std::isnan(ZIminV[c * m_rows + r]))
ZIminV[c * m_rows + r] = z;
size_t cell = c * m_rows + r;
if (z < ZIminV[cell] || std::isnan(ZIminV[cell]))
ZIminV[cell] = z;
}

// "...some grid points of ZImin will go unfilled. To fill these values, we
// rely on computationally inexpensive image inpainting techniques. Image
// inpainting involves the replacement of the empty cells in an image (or
// matrix) with values calculated from other nearby values."
std::vector<double> ZImin_fillV = knnfill(view, ZIminV);

//ABELL - We can eliminate this copy if we're OK with not writing
// both the filled and non-filled array to output.
std::vector<double> ZImin_fillV = ZIminV;
knnfill(view, ZImin_fillV);

if (!m_args->m_dir.empty())
{
Expand Down Expand Up @@ -556,7 +567,10 @@ std::vector<double> SMRFilter::createZIpro(PointViewPtr view,

// "These cells are then inpainted according to the same process described
// previously, producing a provisional DEM (ZIpro)."
std::vector<double> ZIpro_fillV = knnfill(view, ZIproV);
//ABELL - We can eliminate this copy if we're OK with not writing
// both the filled and non-filled array to output.
std::vector<double> ZIpro_fillV = ZIproV;
knnfill(view, ZIpro_fillV);

if (!m_args->m_dir.empty())
{
Expand All @@ -575,9 +589,14 @@ std::vector<double> SMRFilter::createZIpro(PointViewPtr view,
}

// Fill voids with the average of eight nearest neighbors.
std::vector<double> SMRFilter::knnfill(PointViewPtr view,
std::vector<double> const& cz)
void SMRFilter::knnfill(PointViewPtr view, std::vector<double>& cz)
{
//ABELL - This potentially means moving a lot of data from the raster
// to the temporary view. This can be improved by either
// 1) using some method other than a KD tree to find neighbors
// 2) build a KDtree from the raster data directly, rather than moving it
// to a view.

// Create a temporary PointView that encodes our raster values so that we
// can construct a 2D KDIndex and perform nearest neighbor searches.
PointViewPtr temp = view->makeNew();
Expand All @@ -586,37 +605,41 @@ std::vector<double> SMRFilter::knnfill(PointViewPtr view,
{
for (int r = 0; r < m_rows; ++r)
{
if (std::isnan(cz[c * m_rows + r]))
size_t cell = c * m_rows + r;
double val = cz[cell];
if (std::isnan(val))
continue;

temp->setField(Id::X, i,
m_bounds.minx + (c + 0.5) * m_args->m_cell);
temp->setField(Id::Y, i,
m_bounds.miny + (r + 0.5) * m_args->m_cell);
temp->setField(Id::Z, i, cz[c * m_rows + r]);
temp->setField(Id::Z, i, val);
i++;
}
}

// https://github.com/PDAL/PDAL/issues/2794#issuecomment-625297062
if (!temp->size())
return;

KD2Index& kdi = temp->build2dIndex();

// Where the raster has voids (i.e., NaN), we search for that cell's eight
// nearest neighbors, and fill the void with the average value of the
// neighbors.
std::vector<double> out = cz;
for (int c = 0; c < m_cols; ++c)
{
for (int r = 0; r < m_rows; ++r)
{
if (!std::isnan(out[c * m_rows + r]))
size_t cell = c * m_rows + r;
if (!std::isnan(cz[cell]))
continue;

double x = m_bounds.minx + (c + 0.5) * m_args->m_cell;
double y = m_bounds.miny + (r + 0.5) * m_args->m_cell;
int k = 8;
PointIdList neighbors(k);
std::vector<double> sqr_dists(k);
kdi.knnSearch(x, y, k, &neighbors, &sqr_dists);
const int k = 8;
PointIdList neighbors = kdi.neighbors(x, y, k);

double M1(0.0);
size_t j(0);
Expand All @@ -626,12 +649,9 @@ std::vector<double> SMRFilter::knnfill(PointViewPtr view,
double delta = temp->getFieldAs<double>(Id::Z, n) - M1;
M1 += (delta / j);
}

out[c * m_rows + r] = M1;
cz[cell] = M1;
}
}

return out;
}

// Iteratively open the estimated surface. progressiveFilter can be used to
Expand Down
7 changes: 3 additions & 4 deletions filters/SMRFilter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ class PDAL_DLL SMRFilter : public Filter
public:
SMRFilter();
~SMRFilter();
SMRFilter& operator=(const SMRFilter&) = delete;
SMRFilter(const SMRFilter&) = delete;

std::string getName() const;

Expand Down Expand Up @@ -76,12 +78,9 @@ class PDAL_DLL SMRFilter : public Filter
std::vector<int> const&,
std::vector<int> const&,
std::vector<int> const&);
std::vector<double> knnfill(PointViewPtr, std::vector<double> const&);
void knnfill(PointViewPtr, std::vector<double>&);
std::vector<int> progressiveFilter(std::vector<double> const&, double,
double);

SMRFilter& operator=(const SMRFilter&); // not implemented
SMRFilter(const SMRFilter&); // not implemented
};

} // namespace pdal

0 comments on commit ea30f6b

Please sign in to comment.