Skip to content

Commit

Permalink
Use X, Y and Z in geocentric transform.
Browse files Browse the repository at this point in the history 
Remove dead code.
  • Loading branch information
@abellgithub
abellgithub committed May 22, 2019
1 parent 305f56c commit 55a5147
Show file tree
Hide file tree
Showing 5 changed files with 149 additions and 125 deletions.
18 changes: 12 additions & 6 deletions pdal/GDALUtils.hpp
View file
Expand Up @@ -82,7 +82,8 @@ class SpatialRef
SpatialRef(const std::string& srs)
{
newRef(OSRNewSpatialReference(""));
OSRSetFromUserInput(get(), srs.data());
if (OSRSetFromUserInput(get(), srs.data()) != OGRERR_NONE)
m_ref.reset();
}

void setFromLayer(OGRLayerH layer)
Expand All @@ -103,11 +104,16 @@ class SpatialRef
{ return m_ref.get(); }
std::string wkt() const
{
char *pszWKT = NULL;
OSRExportToWkt(m_ref.get(), &pszWKT);
bool valid = (bool)*pszWKT;
std::string output(pszWKT);
CPLFree(pszWKT);
std::string output;

if (m_ref.get())
{
char *pszWKT = NULL;
OSRExportToWkt(m_ref.get(), &pszWKT);
bool valid = (bool)*pszWKT;
output = pszWKT;
CPLFree(pszWKT);
}
return output;
}

Expand Down
162 changes: 92 additions & 70 deletions plugins/i3s/io/EsriReader.cpp
View file
Expand Up @@ -36,6 +36,7 @@

#include <Eigen/Geometry>
#include <ogr_spatialref.h>
#include <pdal/GDALUtils.hpp>

#include "../lepcc/src/include/lepcc_types.h"

Expand All @@ -55,26 +56,21 @@ void EsriReader::addArgs(ProgramArgs& args)
m_args.min_density, -1.0);
args.add("max_density", "Maximum point density",
m_args.max_density, -1.0);

}

void EsriReader::initialize(PointTableRef table)
{
//create proper density if min was set but max wasn't
if(m_args.min_density >= 0 && m_args.max_density < 0)
if (m_args.min_density >= 0 && m_args.max_density < 0)
m_args.max_density = (std::numeric_limits<double>::max)();

//create dimensions map for future lookup
if (!m_args.dimensions.empty())
{
for (std::string& dim : m_args.dimensions)
{
std::transform(
dim.begin(),
dim.end(),
dim.begin(),
[](unsigned char c){ return std::toupper(c); });
if(esriDims.find(dim) != esriDims.end())
dim = Utils::toupper(dim);
if (esriDims.find(dim) != esriDims.end())
m_dimensions[dim] = esriDims.at(dim);
else
m_dimensions[dim];
Expand Down Expand Up @@ -107,12 +103,10 @@ void EsriReader::initialize(PointTableRef table)
//create const for looking into
const NL::json jsonBody = m_info;

//create pdal Bounds
m_bounds = createBounds();

//find version
if (jsonBody["store"].contains("version"))
m_version = Version(jsonBody["store"]["version"].get<std::string>());

if (Version("2.0") < m_version || m_version < Version("1.6"))
log()->get(LogLevel::Warning) << "This version may not work with "
"the current implementation of i3s/slpk reader" << std::endl;
Expand Down Expand Up @@ -142,35 +136,96 @@ void EsriReader::initialize(PointTableRef table)

//create spatial reference objects
NL::json wkid = m_info["spatialReference"]["wkid"];
std::string spatialStr = "EPSG:";
int system(0);
if (wkid.is_string())
spatialStr += wkid.get<std::string>();
else if (wkid.is_number_unsigned())
spatialStr += std::to_string(wkid.get<uint64_t>());
m_nativeSrs = SpatialReference(spatialStr);
{
std::string sval = wkid.get<std::string>();
try
{
system = std::stoi(sval);
}
catch (...)
{}
if (system < 2000)
throwError("Invalid wkid string '" + sval + "' for spatial "
"reference.");
}
else if (wkid.is_number())
{
system = wkid.get<int64_t>();
if (system < 2000)
throwError("Invalid wkid value '" + std::to_string(system) +
"' for spatial reference.");
}
std::string systemString("EPSG:" + std::to_string(system));
m_nativeSrs.set(systemString);
if (m_nativeSrs.empty())
throwError("Unable to create spatial reference for i3s for '" +
systemString + "'.");
setSpatialReference(m_nativeSrs);

m_ecefSrs = SpatialReference("EPSG:4978");
createEcefTransform();
createBounds();
}

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);
void EsriReader::createEcefTransform()
{
// All comparisons of bounding boxes are done in geocentric coords.
gdal::SpatialRef nativeSrs(m_nativeSrs.getWKT());
gdal::SpatialRef ecefSrs("EPSG:4978");
m_toEcefTransform = OCTNewCoordinateTransformation(nativeSrs.get(),
ecefSrs.get());
}

//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)
//Create bounding box that the user specified
void EsriReader::createBounds()
{
const double mn((std::numeric_limits<double>::lowest)());
const double mx((std::numeric_limits<double>::max)());
if (m_args.bounds.is3d())
{
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);
m_bounds = m_args.bounds.to3d();
m_ecefBounds = m_bounds;
OCTTransform(m_toEcefTransform, 1,
&m_ecefBounds.minx, &m_ecefBounds.miny, &m_ecefBounds.minz);
OCTTransform(m_toEcefTransform, 1,
&m_ecefBounds.maxx, &m_ecefBounds.maxy, &m_ecefBounds.maxz);
}
else
{
// Distance to the center of the earth is 6.3 million meters, so
// 10 million should be a reasonable max for ECEF.
BOX2D b = m_args.bounds.to2d();
if (b.empty())
{
m_bounds = BOX3D(mn, mn, mn, mx, mx, mx);
m_ecefBounds = BOX3D(-10e6, -10e6, -10e6, 10e6, 10e6, 10e6);
}
else
{
m_bounds = BOX3D(b); // Will set z values to 0.
m_ecefBounds = m_bounds;
OCTTransform(m_toEcefTransform, 1,
&m_ecefBounds.minx, &m_ecefBounds.miny, &m_ecefBounds.minz);
OCTTransform(m_toEcefTransform, 1,
&m_ecefBounds.maxx, &m_ecefBounds.maxy, &m_ecefBounds.maxz);
m_bounds.minz = mn;
m_bounds.maxz = mx;
m_ecefBounds.minz = -10e6;
m_ecefBounds.maxz = 10e6;
}
}
// Transformation can invert coordinates
if (m_ecefBounds.minx > m_ecefBounds.maxx)
std::swap(m_ecefBounds.minx, m_ecefBounds.maxx);
if (m_ecefBounds.miny > m_ecefBounds.maxy)
std::swap(m_ecefBounds.miny, m_ecefBounds.maxy);
if (m_ecefBounds.minz > m_ecefBounds.maxz)
std::swap(m_ecefBounds.minz, m_ecefBounds.maxz);
}


void EsriReader::addDimensions(PointLayoutPtr layout)
{
if (!m_info.contains("attributeStorageInfo"))
Expand Down Expand Up @@ -259,7 +314,6 @@ void EsriReader::addDimensions(PointLayoutPtr layout)

void EsriReader::ready(PointTableRef)
{

//output arguments for debugging
log()->get(LogLevel::Debug) << "filename: " <<
m_filename << std::endl;
Expand All @@ -272,10 +326,9 @@ void EsriReader::ready(PointTableRef)
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 Down Expand Up @@ -323,7 +376,7 @@ point_count_t EsriReader::read(PointViewPtr view, point_count_t count)
//the tree and test if it overlaps with the bounds created by user.
//If it's a leaf node(the highest resolution) and it overlaps, add
//it to the list of nodes to be pulled later.
void EsriReader::traverseTree(NL::json& page, int index,
void EsriReader::traverseTree(NL::json page, int index,
std::vector<int>& nodes, int depth, int pageIndex)
{
// find node information
Expand Down Expand Up @@ -407,6 +460,7 @@ void EsriReader::traverseTree(NL::json& page, int index,
}
}


//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.
Expand All @@ -428,16 +482,15 @@ BOX3D EsriReader::createCube(const NL::json& base)
//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;
return BOX3D(minx, miny, minz, maxx, maxy, maxz);
}


void EsriReader::createView(std::string localUrl, int nodeIndex, PointView& view)
void EsriReader::createView(std::string localUrl, int nodeIndex,
PointView& view)
{
// pull the geometries to start
const std::string geomUrl = localUrl + "/geometries/";
Expand Down Expand Up @@ -574,7 +627,6 @@ void EsriReader::createView(std::string localUrl, int nodeIndex, PointView& view
view.setField(Dimension::Id::NumberOfReturns,
startId + i, numReturns);
}

}
else
{
Expand All @@ -598,41 +650,11 @@ void EsriReader::createView(std::string localUrl, int nodeIndex, PointView& view
}


//Create bounding box that the user specified
BOX3D EsriReader::createBounds()
{
if (m_args.bounds.is3d())
return m_args.bounds.to3d();

// If empty, return maximal extents to select everything.
const BOX2D b(m_args.bounds.to2d());
if (b.empty())
{
double mn((std::numeric_limits<double>::lowest)());
double mx((std::numeric_limits<double>::max)());
return BOX3D(mn, mn, mn, mx, mx, mx);
}

// Finally if 2D and non-empty, convert to 3D with all-encapsulating
// Z-values.
return BOX3D(
b.minx, b.miny, (std::numeric_limits<double>::lowest)(),
b.maxx, b.maxy, (std::numeric_limits<double>::max)());
}


void EsriReader::done(PointTableRef)
{
m_stream.reset();

if (m_nativeRef)
OSRDestroySpatialReference(m_nativeRef);
if (m_ecefRef)
OSRDestroySpatialReference(m_ecefRef);
if (m_toEcefTransform)
OCTDestroyCoordinateTransformation(m_toEcefTransform);
if (m_toNativeTransform)
OCTDestroyCoordinateTransformation(m_toNativeTransform);
}

} //namespace pdal
Expand Down
14 changes: 5 additions & 9 deletions plugins/i3s/io/EsriReader.hpp
View file
Expand Up @@ -77,9 +77,6 @@ std::map<std::string, pdal::Dimension::Type> const dimTypes

class PDAL_DLL EsriReader : public Reader
{
public:
BOX3D createBounds();

protected:
virtual void initInfo() = 0;
virtual std::vector<char> fetchBinary(std::string url, std::string attNum,
Expand Down Expand Up @@ -160,12 +157,7 @@ class PDAL_DLL EsriReader : public Reader

typedef void* ReferencePtr;
typedef void* TransformPtr;
ReferencePtr m_nativeRef;
ReferencePtr m_ecefRef;

TransformPtr m_toEcefTransform;
TransformPtr m_toNativeTransform;


struct dimData
{
Expand All @@ -186,8 +178,12 @@ class PDAL_DLL EsriReader : public Reader
void createView(std::string localUrl, int nodeIndex, PointView& view);
BOX3D createCube(const NL::json& base);
BOX3D parseBox(const NL::json& base);
void traverseTree(NL::json& page, int index, std::vector<int>& nodes,
void traverseTree(NL::json page, int index, std::vector<int>& nodes,
int depth, int pageIndex);

private:
void createBounds();
void createEcefTransform();
};

} // namespace pdal
Expand Down

0 comments on commit 55a5147

Please sign in to comment.