Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into issue-2390
Browse files Browse the repository at this point in the history
  • Loading branch information
@abellgithub
abellgithub committed Mar 12, 2019
2 parents d291352 + 3110c00 commit 7854ba1
Show file tree
Hide file tree
Showing 10 changed files with 376 additions and 170 deletions.
11 changes: 8 additions & 3 deletions pdal/compression/ZstdCompression.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,10 @@ class ZstdCompressorImpl
char m_tmpbuf[CHUNKSIZE];
BlockCb m_cb;

ZstdCompressorImpl(BlockCb cb) : m_cb(cb)
ZstdCompressorImpl(BlockCb cb, int compressionLevel) : m_cb(cb)
{
m_strm = ZSTD_createCStream();
ZSTD_initCStream(m_strm, 15);
ZSTD_initCStream(m_strm, compressionLevel);
}

~ZstdCompressorImpl()
Expand Down Expand Up @@ -92,7 +92,12 @@ class ZstdCompressorImpl
};

ZstdCompressor::ZstdCompressor(BlockCb cb) :
m_impl(new ZstdCompressorImpl(cb))
m_impl(new ZstdCompressorImpl(cb, 15))
{}


ZstdCompressor::ZstdCompressor(BlockCb cb, int compressionLevel) :
m_impl(new ZstdCompressorImpl(cb, compressionLevel))
{}


Expand Down
1 change: 1 addition & 0 deletions pdal/compression/ZstdCompression.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ class ZstdCompressor : public Compressor
{
public:
PDAL_DLL ZstdCompressor(BlockCb cb);
PDAL_DLL ZstdCompressor(BlockCb cb, int compressionLevel);
PDAL_DLL ~ZstdCompressor();

PDAL_DLL void compress(const char *buf, size_t bufsize);
Expand Down
165 changes: 71 additions & 94 deletions plugins/i3s/io/EsriReader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,7 @@ void EsriReader::initialize(PointTableRef table)
//create pdal Bounds
m_bounds = createBounds();


//find version
if (jsonBody["store"].isMember("version"))
m_version = Version(jsonBody["store"]["version"].asString());
Expand Down Expand Up @@ -144,21 +145,30 @@ void EsriReader::initialize(PointTableRef table)
"by this driver"));
}

//output arguments for debugging
log()->get(LogLevel::Debug) << "filename: " <<
m_filename << std::endl;
log()->get(LogLevel::Debug) << "threads: " <<
m_args.threads << std::endl;
log()->get(LogLevel::Debug) << "bounds: " <<
m_args.bounds << std::endl;
log()->get(LogLevel::Debug) << "min_density: " <<
m_args.min_density << std::endl;
log()->get(LogLevel::Debug) << "max_density: " <<
m_args.max_density << std::endl;
log()->get(LogLevel::Debug) << "dimensions: " << std::endl;
for (std::string& dim : m_args.dimensions)
//create spatial reference objects
Json::Value spatialJson = m_info["spatialReference"];
std::string spatialStr = "EPSG:" + spatialJson["wkid"].asString();
m_nativeSrs = SpatialReference(spatialStr);
setSpatialReference(m_nativeSrs);

m_ecefSrs = SpatialReference("EPSG:4978");

m_ecefRef = OSRNewSpatialReference(m_ecefSrs.getWKT().c_str());
m_nativeRef = OSRNewSpatialReference(m_nativeSrs.getWKT().c_str());
m_toEcefTransform = OCTNewCoordinateTransformation(m_nativeRef, m_ecefRef);
m_toNativeTransform = OCTNewCoordinateTransformation(m_ecefRef,
m_nativeRef);

//create bounds in ECEF
double minx(m_bounds.minx), maxx(m_bounds.maxx), miny(m_bounds.miny), maxy(m_bounds.maxy), minz(m_bounds.minz), maxz(m_bounds.maxz);

for (std::size_t i(0); i < 8; ++i)
{
log()->get(LogLevel::Debug) << " -" << dim <<std::endl;
double a = (i & 1 ? minx: maxx);
double b = (i & 2 ? miny: maxy);
double c = (i & 4 ? minz: maxz);
OCTTransform(m_toEcefTransform, 1, &a, &b, 0);
m_ecefBounds.grow(a, b, c);
}
}

Expand Down Expand Up @@ -250,18 +260,23 @@ void EsriReader::addDimensions(PointLayoutPtr layout)

void EsriReader::ready(PointTableRef)
{
Json::Value spatialJson = m_info["spatialReference"];
std::string spatialStr = "EPSG:" + spatialJson["wkid"].asString();
m_nativeSrs = SpatialReference(spatialStr);
setSpatialReference(m_nativeSrs);

m_ecefSrs = SpatialReference("EPSG:4978");

m_ecefRef = OSRNewSpatialReference(m_ecefSrs.getWKT().c_str());
m_nativeRef = OSRNewSpatialReference(m_nativeSrs.getWKT().c_str());
m_toEcefTransform = OCTNewCoordinateTransformation(m_nativeRef, m_ecefRef);
m_toNativeTransform = OCTNewCoordinateTransformation(m_ecefRef,
m_nativeRef);
//output arguments for debugging
log()->get(LogLevel::Debug) << "filename: " <<
m_filename << std::endl;
log()->get(LogLevel::Debug) << "threads: " <<
m_args.threads << std::endl;
log()->get(LogLevel::Debug) << "bounds: " <<
m_args.bounds << std::endl;
log()->get(LogLevel::Debug) << "min_density: " <<
m_args.min_density << std::endl;
log()->get(LogLevel::Debug) << "max_density: " <<
m_args.max_density << std::endl;
log()->get(LogLevel::Debug) << "dimensions: " << std::endl;
for (std::string& dim : m_args.dimensions)
{
log()->get(LogLevel::Debug) << " -" << dim <<std::endl;
}
}


Expand All @@ -279,6 +294,7 @@ point_count_t EsriReader::read(PointViewPtr view, point_count_t count)
// Build the node list that will tell us which nodes overlap with bounds
std::vector<int> nodes;
const Json::Value initJson = fetchJson(m_filename + "/nodepages/0");
log()->get(LogLevel::Debug) << "Traversing metadata" << std::endl;
traverseTree(initJson, 0, nodes, 0, 0);

// Create view with overlapping nodes at desired depth
Expand All @@ -292,10 +308,11 @@ point_count_t EsriReader::read(PointViewPtr view, point_count_t count)
std::string localUrl;

localUrl = m_filename + "/nodes/" + std::to_string(nodes[i]);
int nodeIndex(nodes[i]);

p.add([localUrl, this, &view]()
p.add([localUrl, nodeIndex, this, &view]()
{
createView(localUrl, *view);
createView(localUrl, nodeIndex, *view);
});
}
p.await();
Expand Down Expand Up @@ -333,14 +350,16 @@ void EsriReader::traverseTree(Json::Value page, int index,
m_maxNode = firstChild + cCount - 1;
}

// create box from OBB information and check for overlaps
BOX3D nodeBox = parseBox(page["nodes"][index]);
bool overlap = m_bounds.overlaps(nodeBox);

//if it doesn't overlap, then none of the nodes in this subtree will
BOX3D nodeBox = createCube(page["nodes"][index]);
bool overlap = m_ecefBounds.overlaps(nodeBox);

// if it doesn't overlap, then none of the nodes in this subtree will
if (!overlap)
return;
// if it's a child node and the density hasn't been set, add leaf nodes


// if it's a child node and we're fetching full density, add leaf nodes
if (m_args.max_density == -1 && m_args.min_density == -1 && cCount == 0)
{
nodes.push_back(name);
Expand Down Expand Up @@ -388,80 +407,36 @@ void EsriReader::traverseTree(Json::Value page, int index,
}
}


// Create the BOX3D for the node. This will be used for testing overlap.
BOX3D EsriReader::parseBox(Json::Value base)
//Finds a sphere from the given center and halfsize vector of the OBB
//and makes a cube around it. This should help with collision detection
//and pruning of nodes before fetching binaries.
BOX3D EsriReader::createCube(Json::Value base)
{
// Pull XYZ in lat/lon.
Json::Value center = base["obb"]["center"];
double x = center[0].asDouble();
double y = center[1].asDouble();
double z = center[2].asDouble();

// Convert to ECEF so the OBB offsets in meters will match up.
OCTTransform(m_toEcefTransform, 1, &x, &y, &z);

Json::Value hsize = base["obb"]["halfSize"];
const double hx = hsize[0].asDouble();
const double hy = hsize[1].asDouble();
const double hz = hsize[2].asDouble();

Json::Value quat = base["obb"]["quaternion"];
const double w = quat[0].asDouble();
const double i = quat[1].asDouble();
const double j = quat[2].asDouble();
const double k = quat[3].asDouble();

// Create quat object and normalize it for use
Eigen::Quaterniond q(w, i, j, k);
q.normalize();

// Here we'll convert our halfsizes to x, y, and z vectors in respective
// directions, which will give us new bounding planes
Eigen::Vector3d vecx(hx, 0, 0);
Eigen::Vector3d vecy(0, hy, 0);
Eigen::Vector3d vecz(0, 0, hz);

// Create quaternion-like vectors
Eigen::Quaterniond quatVecX, pxmin, quatVecY, pymin, quatVecZ, pzmin;
quatVecX.w() = 0;
quatVecY.w() = 0;
quatVecZ.w() = 0;
quatVecX.vec() = vecx;
quatVecY.vec() = vecy;
quatVecZ.vec() = vecz;

// Rotate all the individual vectors
// gives us offset for the of the new x/y/zmax/min planes
Eigen::Quaterniond rotx = q * quatVecX * q.inverse();
Eigen::Quaterniond roty = q * quatVecY * q.inverse();
Eigen::Quaterniond rotz = q * quatVecZ * q.inverse();

BOX3D bounds;
for (std::size_t i(0); i < 8; ++i)
{
double a(x), b(y), c(z);

a += rotx.vec()[0] * (i & 1 ? 1.0 : -1.0);
b += rotx.vec()[1] * (i & 1 ? 1.0 : -1.0);
c += rotx.vec()[2] * (i & 1 ? 1.0 : -1.0);

a += roty.vec()[0] * (i & 2 ? 1.0 : -1.0);
b += roty.vec()[1] * (i & 2 ? 1.0 : -1.0);
c += roty.vec()[2] * (i & 2 ? 1.0 : -1.0);

a += rotz.vec()[0] * (i & 4 ? 1.0 : -1.0);
b += rotz.vec()[1] * (i & 4 ? 1.0 : -1.0);
c += rotz.vec()[2] * (i & 4 ? 1.0 : -1.0);
double hx = hsize[0].asDouble();
double hy = hsize[1].asDouble();
double hz = hsize[2].asDouble();

OCTTransform(m_toNativeTransform, 1, &a, &b, &c);
bounds.grow(a, b, c);
}
return bounds;
//transform (x,y,z) to ECEF to match the half sizes in meters.
OCTTransform(m_toEcefTransform, 1, &x, &y, &z);
//take half size vector and find magnitude of it multiplied by sqrt(2)
double r = std::sqrt(2) * std::sqrt(std::pow(hx, 2) + std::pow(hy, 2) + std::pow(hz, 2));
//create cube around this radius
double maxx(x+r), maxy(y+r), maxz(z+r), minx(x-r), miny(y-r), minz(z-r);
BOX3D nodeBox(minx,miny,minz,maxx,maxy,maxz);

//returning in ecef
return nodeBox;
}


void EsriReader::createView(std::string localUrl, PointView& view)
void EsriReader::createView(std::string localUrl, int nodeIndex, PointView& view)
{
// pull the geometries to start
const std::string geomUrl = localUrl + "/geometries/";
Expand All @@ -483,6 +458,7 @@ void EsriReader::createView(std::string localUrl, PointView& view)
std::lock_guard<std::mutex> lock(m_mutex);
startId = view.size();


for (uint64_t j = 0; j < xyz.size(); ++j)
{
double x = xyz[j].x;
Expand Down Expand Up @@ -572,6 +548,7 @@ void EsriReader::createView(std::string localUrl, PointView& view)
}
else if (dimId == Dimension::Id::NumberOfReturns)
{

const std::vector<char> data = fetchBinary(
attrUrl, std::to_string(key), ".bin.gz");

Expand Down
7 changes: 6 additions & 1 deletion plugins/i3s/io/EsriReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@

#include <vector>
#include <map>
#include <math.h>

#include <json/json.h>
#include <arbiter/arbiter.hpp>
Expand Down Expand Up @@ -144,10 +145,12 @@ class PDAL_DLL EsriReader : public Reader
std::unique_ptr<arbiter::Arbiter> m_arbiter;
arbiter::gzip::Decompressor m_decomp;


EsriArgs m_args;
Json::Value m_info;
std::mutex m_mutex;
BOX3D m_bounds;
BOX3D m_ecefBounds;
std::map<std::string, Dimension::Id> m_dimensions;
int m_nodeCap;
int m_maxNode = 0;
Expand All @@ -165,6 +168,7 @@ class PDAL_DLL EsriReader : public Reader
TransformPtr m_toEcefTransform;
TransformPtr m_toNativeTransform;


struct dimData
{
int key;
Expand All @@ -181,7 +185,8 @@ class PDAL_DLL EsriReader : public Reader
virtual void ready(PointTableRef table) override;
virtual point_count_t read(PointViewPtr view, point_count_t count) override;
virtual void done(PointTableRef table) override;
void createView(std::string localUrl, PointView& view);
void createView(std::string localUrl, int nodeIndex, PointView& view);
BOX3D createCube(Json::Value base);
BOX3D parseBox(Json::Value base);
void traverseTree(Json::Value page, int index, std::vector<int>& nodes,
int depth, int pageIndex);
Expand Down
3 changes: 2 additions & 1 deletion plugins/python/filters/PythonFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,8 @@ void PythonFilter::ready(PointTableRef table)
m_source = FileUtils::readFileIntoString(m_scriptFile);
plang::Environment::get()->set_stdout(log()->getLogStream());
m_script = new plang::Script(m_source, m_module, m_function);
m_pythonMethod = new plang::Invocation(*m_script);

m_pythonMethod = new plang::Invocation(*m_script);
m_pythonMethod->compile();
m_totalMetadata = table.metadata();
}
Expand Down

0 comments on commit 7854ba1

Please sign in to comment.