Skip to content

Commit

Permalink
Merge branch 'extended-local-minimum' into derivative-refactor
Browse files Browse the repository at this point in the history
  • Loading branch information
@chambbj
chambbj committed Nov 17, 2016
2 parents 626b0a4 + 5bdc780 commit 1565af6
Show file tree
Hide file tree
Showing 6 changed files with 226 additions and 241 deletions.
112 changes: 8 additions & 104 deletions filters/mongus/MongusFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,6 @@

#include <Eigen/Dense>

#include "gdal_priv.h" // For File I/O
#include "gdal_version.h" // For version info
#include "ogr_spatialref.h" //For Geographic Information/Transformations

namespace pdal
{

Expand Down Expand Up @@ -189,98 +185,6 @@ Eigen::MatrixXd MongusFilter::computeSpline(Eigen::MatrixXd x_prev,
return S;
}

void MongusFilter::writeMatrix(Eigen::MatrixXd data, std::string filename, double cell_size, PointViewPtr view)
{
int cols = data.cols();
int rows = data.rows();

GDALAllRegister();

GDALDataset *mpDstDS(0);

char **papszMetadata;

// parse the format driver, hardcoded for the time being
std::string tFormat("GTIFF");
const char *pszFormat = tFormat.c_str();
GDALDriver* tpDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

// try to create a file of the requested format
if (tpDriver != NULL)
{
papszMetadata = tpDriver->GetMetadata();
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
{
char **papszOptions = NULL;

mpDstDS = tpDriver->Create(filename.c_str(), cols, rows, 1,
GDT_Float32, papszOptions);

// set the geo transformation
double adfGeoTransform[6];
adfGeoTransform[0] = m_bounds.minx; // - 0.5*m_GRID_DIST_X;
adfGeoTransform[1] = cell_size;
adfGeoTransform[2] = 0.0;
adfGeoTransform[3] = m_bounds.maxy; // + 0.5*m_GRID_DIST_Y;
adfGeoTransform[4] = 0.0;
adfGeoTransform[5] = -1 * cell_size;
mpDstDS->SetGeoTransform(adfGeoTransform);

// set the projection
mpDstDS->SetProjection(view->spatialReference().getWKT().c_str());
}
}

// if we have a valid file
if (mpDstDS)
{
// loop over the raster and determine max slope at each location
int cs = 0, ce = cols;
int rs = 0, re = rows;
float *poRasterData = new float[cols*rows];
for (auto i=0; i<cols*rows; i++)
{
poRasterData[i] = std::numeric_limits<float>::min();
}

#pragma omp parallel for
for (auto c = cs; c < ce; ++c)
{
for (auto r = rs; r < re; ++r)
{
if (data(r, c) == 0.0 || std::isnan(data(r, c)) || data(r, c) == std::numeric_limits<double>::max())
continue;
poRasterData[(r * cols) + c] =
data(r, c);
}
}

// write the data
if (poRasterData)
{
GDALRasterBand *tBand = mpDstDS->GetRasterBand(1);

tBand->SetNoDataValue(std::numeric_limits<float>::min());

if (cols > 0 && rows > 0)
#if GDAL_VERSION_MAJOR <= 1
tBand->RasterIO(GF_Write, 0, 0, cols, rows,
poRasterData, cols, rows,
GDT_Float32, 0, 0);
#else

int ret = tBand->RasterIO(GF_Write, 0, 0, cols, rows,
poRasterData, cols, rows,
GDT_Float32, 0, 0, 0);
#endif
}

GDALClose((GDALDatasetH) mpDstDS);

delete [] poRasterData;
}
}

void MongusFilter::writeControl(Eigen::MatrixXd cx, Eigen::MatrixXd cy, Eigen::MatrixXd cz, std::string filename)
{
using namespace Dimension;
Expand Down Expand Up @@ -561,7 +465,7 @@ std::vector<PointId> MongusFilter::processGround(PointViewPtr view)
char buffer[256];
sprintf(buffer, "interp_surface_%d.laz", l);
std::string name(buffer);
// writeMatrix(surface, name, cur_cell_size, view);
// eigen::writeMatrix(surface, name, cur_cell_size, view, m_bounds);
writeControl(x_samp, y_samp, surface, name);

char bufm[256];
Expand All @@ -578,13 +482,13 @@ std::vector<PointId> MongusFilter::processGround(PointViewPtr view)
char rbuf[256];
sprintf(rbuf, "residual_%d.laz", l);
std::string rbufn(rbuf);
// writeMatrix(R, rbufn, cur_cell_size, view);
// eigen::writeMatrix(R, rbufn, cur_cell_size, view, m_bounds);
writeControl(x_samp, y_samp, R, rbufn);

char mbuf[256];
sprintf(mbuf, "median_%d.laz", l);
std::string mbufn(mbuf);
// writeMatrix(M, mbufn, cur_cell_size, view);
// eigen::writeMatrix(M, mbufn, cur_cell_size, view, m_bounds);
writeControl(x_samp, y_samp, M, mbufn);

char buf2[256];
Expand All @@ -607,27 +511,27 @@ std::vector<PointId> MongusFilter::processGround(PointViewPtr view)
char buffer[256];
sprintf(buffer, "final_surface.tif");
std::string name(buffer);
writeMatrix(surface, name, m_cellSize, view);
eigen::writeMatrix(surface, name, m_cellSize, view, m_bounds);
//
// char rbuf[256];
// sprintf(rbuf, "final_residual.tif");
// std::string rbufn(rbuf);
// writeMatrix(R, rbufn, cur_cell_size, view);
// eigen::writeMatrix(R, rbufn, cur_cell_size, view, m_bounds);
//
// char obuf[256];
// sprintf(obuf, "final_opened.tif");
// std::string obufn(obuf);
// writeMatrix(maxZ, obufn, cur_cell_size, view);
// eigen::writeMatrix(maxZ, obufn, cur_cell_size, view, m_bounds);
//
// char Tbuf[256];
// sprintf(Tbuf, "final_tophat.tif");
// std::string Tbufn(Tbuf);
// writeMatrix(T, Tbufn, cur_cell_size, view);
// eigen::writeMatrix(T, Tbufn, cur_cell_size, view, m_bounds);
//
// char tbuf[256];
// sprintf(tbuf, "final_thresh.tif");
// std::string tbufn(tbuf);
// writeMatrix(t, tbufn, cur_cell_size, view);
// eigen::writeMatrix(t, tbufn, cur_cell_size, view, m_bounds);
}

// apply final filtering (top hat) using raw points against TPS
Expand Down
2 changes: 0 additions & 2 deletions filters/mongus/MongusFilter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,6 @@ class PDAL_DLL MongusFilter : public Filter
virtual void addArgs(ProgramArgs& args);
int getColIndex(double x, double cell_size);
int getRowIndex(double y, double cell_size);
void writeMatrix(Eigen::MatrixXd data, std::string filename,
double cell_size, PointViewPtr view);
Eigen::MatrixXd computeSpline(Eigen::MatrixXd x_prev,
Eigen::MatrixXd y_prev,
Eigen::MatrixXd z_prev,
Expand Down
125 changes: 14 additions & 111 deletions filters/smrf/SMRFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,10 +44,6 @@
#include <Eigen/Dense>
#include <Eigen/Sparse>

#include "gdal_priv.h" // For File I/O
#include "gdal_version.h" // For version info
#include "ogr_spatialref.h" //For Geographic Information/Transformations

namespace pdal
{
using namespace Eigen;
Expand Down Expand Up @@ -264,12 +260,12 @@ std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
// neighbor (KNN) approach.
MatrixXd ZImin = eigen::createDSM(*view.get(), m_numRows, m_numCols,
m_cellSize, m_bounds);
writeMatrix(ZImin, "zimin.tif", m_cellSize, view);

eigen::writeMatrix(ZImin, "zimin.tif", m_cellSize, view, m_bounds);
// MatrixXd ZImin_painted = inpaintKnn(cx, cy, ZImin);
// MatrixXd ZImin_painted = TPS(cx, cy, ZImin);
MatrixXd ZImin_painted = expandingTPS(cx, cy, ZImin);
writeMatrix(ZImin_painted, "zimin_painted.tif", m_cellSize, view);
eigen::writeMatrix(ZImin_painted, "zimin_painted.tif", m_cellSize, view, m_bounds);

ZImin = ZImin_painted;

Expand All @@ -290,7 +286,7 @@ std::vector<PointId> SMRFilter::processGround(PointViewPtr view)

// paper has low point happening later, i guess it doesn't matter too much, this is where he does it in matlab code
MatrixXi Low = progressiveFilter(-ZImin, m_cellSize, 5.0, 1.0);
writeMatrix(Low.cast<double>(), "zilow.tif", m_cellSize, view);
eigen::writeMatrix(Low.cast<double>(), "zilow.tif", m_cellSize, view, m_bounds);

// matlab code has net cutting occurring here
MatrixXd ZInet = ZImin;
Expand Down Expand Up @@ -320,12 +316,12 @@ std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
ZInet(r, c) = bigOpen(r, c);
}
}
writeMatrix(ZInet, "zinet.tif", m_cellSize, view);
eigen::writeMatrix(ZInet, "zinet.tif", m_cellSize, view, m_bounds);
}

// and finally object detection
MatrixXi Obj = progressiveFilter(ZInet, m_cellSize, m_percentSlope, m_maxWindow);
writeMatrix(Obj.cast<double>(), "ziobj.tif", m_cellSize, view);
eigen::writeMatrix(Obj.cast<double>(), "ziobj.tif", m_cellSize, view, m_bounds);

// STEP 3:

Expand All @@ -343,12 +339,12 @@ std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
if (Obj(i) == 1 || Low(i) == 1 || isNetCell(i) == 1)
ZIpro(i) = std::numeric_limits<double>::quiet_NaN();
}
writeMatrix(ZIpro, "zipro.tif", m_cellSize, view);
eigen::writeMatrix(ZIpro, "zipro.tif", m_cellSize, view, m_bounds);

// MatrixXd ZIpro_painted = inpaintKnn(cx, cy, ZIpro);
// MatrixXd ZIpro_painted = TPS(cx, cy, ZIpro);
MatrixXd ZIpro_painted = expandingTPS(cx, cy, ZIpro);
writeMatrix(ZIpro_painted, "zipro_painted.tif", m_cellSize, view);
eigen::writeMatrix(ZIpro_painted, "zipro_painted.tif", m_cellSize, view, m_bounds);

ZIpro = ZIpro_painted;

Expand Down Expand Up @@ -425,21 +421,21 @@ std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
};

MatrixXd gx = diffX(ZIpro / m_cellSize);
writeMatrix(gx, "gx.tif", m_cellSize, view);
eigen::writeMatrix(gx, "gx.tif", m_cellSize, view, m_bounds);
MatrixXd gy = diffY(ZIpro / m_cellSize);
writeMatrix(gy, "gy.tif", m_cellSize, view);
eigen::writeMatrix(gy, "gy.tif", m_cellSize, view, m_bounds);
MatrixXd gsurfs = (gx.cwiseProduct(gx) + gy.cwiseProduct(gy)).cwiseSqrt();
writeMatrix(gsurfs, "gsurfs.tif", m_cellSize, view);
eigen::writeMatrix(gsurfs, "gsurfs.tif", m_cellSize, view, m_bounds);

// MatrixXd gsurfs_painted = inpaintKnn(cx, cy, gsurfs);
// MatrixXd gsurfs_painted = TPS(cx, cy, gsurfs);
MatrixXd gsurfs_painted = expandingTPS(cx, cy, gsurfs);
writeMatrix(gsurfs_painted, "gsurfs_painted.tif", m_cellSize, view);
eigen::writeMatrix(gsurfs_painted, "gsurfs_painted.tif", m_cellSize, view, m_bounds);

gsurfs = gsurfs_painted;

MatrixXd thresh = (m_threshold + 1.2 * gsurfs.array()).matrix();
writeMatrix(thresh, "thresh.tif", m_cellSize, view);
eigen::writeMatrix(thresh, "thresh.tif", m_cellSize, view, m_bounds);

for (PointId i = 0; i < view->size(); ++i)
{
Expand Down Expand Up @@ -527,7 +523,7 @@ MatrixXi SMRFilter::progressiveFilter(MatrixXd const& ZImin, double cell_size,
if (diff(i) > threshold)
Obj(i) = 1;
}
// writeMatrix(Obj, "obj.tif", m_cellSize, view);
// eigen::writeMatrix(Obj, "obj.tif", m_cellSize, view, m_bounds);

// The algorithm then proceeds to the next window radius (up to the
// maximum), and proceeds as above with the last opened surface acting
Expand Down Expand Up @@ -923,97 +919,4 @@ MatrixXd SMRFilter::expandingTPS(MatrixXd cx, MatrixXd cy, MatrixXd cz)
return S;
}

void SMRFilter::writeMatrix(MatrixXd data, std::string filename,
double cell_size, PointViewPtr view)
{
int cols = data.cols();
int rows = data.rows();

GDALAllRegister();

GDALDataset *mpDstDS = NULL;

char **papszMetadata;

// parse the format driver, hardcoded for the time being
std::string tFormat("GTIFF");
const char *pszFormat = tFormat.c_str();
GDALDriver* tpDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

// try to create a file of the requested format
if (tpDriver != NULL)
{
papszMetadata = tpDriver->GetMetadata();
if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
{
char **papszOptions = NULL;

mpDstDS = tpDriver->Create(filename.c_str(), cols, rows, 1,
GDT_Float32, papszOptions);

// set the geo transformation
double adfGeoTransform[6];
adfGeoTransform[0] = m_bounds.minx; // - 0.5*m_GRID_DIST_X;
adfGeoTransform[1] = cell_size;
adfGeoTransform[2] = 0.0;
adfGeoTransform[3] = m_bounds.maxy; // + 0.5*m_GRID_DIST_Y;
adfGeoTransform[4] = 0.0;
adfGeoTransform[5] = -1 * cell_size;
mpDstDS->SetGeoTransform(adfGeoTransform);

// set the projection
mpDstDS->SetProjection(view->spatialReference().getWKT().c_str());
}
}

// if we have a valid file
if (mpDstDS)
{
// loop over the raster and determine max slope at each location
int cs = 1, ce = cols - 1;
int rs = 1, re = rows - 1;
float *poRasterData = new float[cols*rows];
for (auto i=0; i<cols*rows; i++)
{
poRasterData[i] = std::numeric_limits<float>::min();
}

// #pragma omp parallel for
for (auto c = cs; c < ce; ++c)
{
for (auto r = rs; r < re; ++r)
{
if (data(r, c) == 0.0 || std::isnan(data(r, c)) || data(r, c) == std::numeric_limits<double>::max())
continue;
poRasterData[(r * cols) + c] =
data(r, c);
}
}

// write the data
if (poRasterData)
{
GDALRasterBand *tBand = mpDstDS->GetRasterBand(1);

tBand->SetNoDataValue(std::numeric_limits<float>::min());

if (cols > 0 && rows > 0)
#if GDAL_VERSION_MAJOR <= 1
tBand->RasterIO(GF_Write, 0, 0, cols, rows,
poRasterData, cols, rows,
GDT_Float32, 0, 0);
#else

int ret = tBand->RasterIO(GF_Write, 0, 0, cols, rows,
poRasterData, cols, rows,
GDT_Float32, 0, 0, 0);
#endif
}

GDALClose((GDALDatasetH) mpDstDS);

delete [] poRasterData;
}
}

} // namespace pdal
4 changes: 0 additions & 4 deletions filters/smrf/SMRFilter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,10 +101,6 @@ class PDAL_DLL SMRFilter : public Filter
MatrixXd TPS(MatrixXd cx, MatrixXd cy, MatrixXd cz);
MatrixXd expandingTPS(MatrixXd cx, MatrixXd cy, MatrixXd cz);

// writeMatrix writes out Eigen matrices to GeoTIFFs for debugging.
void writeMatrix(MatrixXd data, std::string filename,
double cell_size, PointViewPtr view);

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

0 comments on commit 1565af6

Please sign in to comment.