From 213aee5669aa94dfe8a65e6bc21cf20159797f90 Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Tue, 8 Oct 2013 13:35:35 +0200
Subject: [PATCH] moved eps calculation to AABB
---
GeoLib/AABB.h | 6 +
Gui/DataView/GeoMapper.cpp | 218 ++++++++++++++++++-------------------
2 files changed, 115 insertions(+), 109 deletions(-)
diff --git a/GeoLib/AABB.h b/GeoLib/AABB.h
index c2ffe775f67..0e542b9fa90 100644
--- a/GeoLib/AABB.h
+++ b/GeoLib/AABB.h
@@ -131,6 +131,12 @@ public:
return containsPoint(other_aabb.getMinPoint()) && containsPoint(other_aabb.getMaxPoint());
}
+ double getEpsFromBoundingBox() const
+ {
+ double eps = sqrt(std::numeric_limits<float>::epsilon());
+ return eps * sqrt(MathLib::sqrDist (&(this->getMinPoint()),&(this->getMaxPoint())));
+ }
+
protected:
PNT_TYPE _min_pnt;
PNT_TYPE _max_pnt;
diff --git a/Gui/DataView/GeoMapper.cpp b/Gui/DataView/GeoMapper.cpp
index 9d759d32759..56d098a062a 100644
--- a/Gui/DataView/GeoMapper.cpp
+++ b/Gui/DataView/GeoMapper.cpp
@@ -80,6 +80,113 @@ void GeoMapper::mapOnMesh(const MeshLib::Mesh* mesh)
}
}
+void GeoMapper::mapData()
+{
+ const std::vector<GeoLib::Point*> *points (this->_geo_objects.getPointVec(this->_geo_name));
+ GeoLib::Station* stn_test = dynamic_cast<GeoLib::Station*>((*points)[0]);
+ bool is_borehole(false);
+ if (stn_test != nullptr && static_cast<GeoLib::StationBorehole*>((*points)[0])->type() == GeoLib::Station::StationType::BOREHOLE)
+ is_borehole = true;
+
+ double min_val(0), max_val(0);
+ if (_mesh)
+ {
+ GeoLib::AABB<GeoLib::Point> bounding_box (_mesh->getNodes().begin(), _mesh->getNodes().end());
+ min_val = bounding_box.getMinPoint()[2];
+ max_val = bounding_box.getMaxPoint()[2];
+ }
+ size_t nPoints (points->size());
+
+ if (!is_borehole)
+ {
+ for (unsigned j=0; j<nPoints; ++j)
+ {
+ GeoLib::Point* pnt ((*points)[j]);
+ (*pnt)[2] = (_grid) ? this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val)
+ : this->getDemElevation((*pnt)[0],(*pnt)[1]);
+ }
+ }
+ else
+ {
+ for (unsigned j=0; j<nPoints; ++j)
+ {
+ GeoLib::Point* pnt ((*points)[j]);
+ double offset = (_grid) ? (this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val) - (*pnt)[2])
+ : (this->getDemElevation((*pnt)[0],(*pnt)[1]) - (*pnt)[2]);
+
+ GeoLib::StationBorehole* borehole = static_cast<GeoLib::StationBorehole*>(pnt);
+ const std::vector<GeoLib::Point*> layers = borehole->getProfile();
+ size_t nLayers = layers.size();
+ for (unsigned k=0; k<nLayers; ++k)
+ {
+ GeoLib::Point* layer_pnt = layers[k];
+ (*layer_pnt)[2] = (*layer_pnt)[2] + offset;
+ }
+ }
+ }
+}
+
+float GeoMapper::getDemElevation(double x, double y) const
+{
+ const double origin_x(_raster->getOrigin()[0]);
+ const double origin_y(_raster->getOrigin()[1]);
+ const double cellsize(_raster->getRasterPixelDistance());
+ const std::size_t width(_raster->getNCols());
+ const std::size_t height(_raster->getNRows());
+
+ if ((x<origin_x) || (x>origin_x+(width*cellsize)) || (y<origin_y) || (y>origin_y+(height*cellsize)))
+ return 0;
+
+ const unsigned x_index = static_cast<unsigned>((x-origin_x)/cellsize);
+ const unsigned y_index = static_cast<unsigned>((y-origin_y)/cellsize);
+
+ return static_cast<float>(*(_raster->begin() + (y_index*width+x_index)));
+}
+
+double GeoMapper::getMeshElevation(double x, double y, double min_val, double max_val) const
+{
+ double coords[3] = {x,y,0};
+ const GeoLib::PointWithID* pnt = _grid->getNearestPoint(coords);
+ const std::vector<MeshLib::Element*> elements (_mesh->getNode(pnt->getID())->getElements());
+ GeoLib::Point* intersection (nullptr);
+
+ for (std::size_t i=0; i<elements.size(); ++i)
+ {
+ if (intersection==nullptr)
+ intersection=GeoLib::triangleLineIntersection(*elements[i]->getNode(0), *elements[i]->getNode(1), *elements[i]->getNode(2), GeoLib::Point(x,y,max_val), GeoLib::Point(x,y,min_val));
+ if (intersection==nullptr && elements[i]->getGeomType() == MeshElemType::QUAD)
+ intersection=GeoLib::triangleLineIntersection(*elements[i]->getNode(0), *elements[i]->getNode(2), *elements[i]->getNode(3), GeoLib::Point(x,y,max_val), GeoLib::Point(x,y,min_val));
+ }
+ if (intersection)
+ return (*intersection)[2];
+ // if something goes wrong, simply take the elevation of the nearest mesh node
+ return (*(_mesh->getNode(pnt->getID())))[2];
+}
+
+GeoLib::Grid<GeoLib::PointWithID>* GeoMapper::getFlatGrid(MeshLib::Mesh const*const mesh, std::vector<GeoLib::PointWithID*> sfc_pnts) const
+{
+ if (mesh->getDimension()<3) //much faster
+ {
+ size_t nNodes (mesh->getNNodes());
+ sfc_pnts.resize(nNodes);
+ const std::vector<MeshLib::Node*> nodes (mesh->getNodes());
+ for (unsigned i(0); i<nNodes; ++i)
+ sfc_pnts[i] = new GeoLib::PointWithID(nodes[i]->getCoords(), nodes[i]->getID());
+ }
+ else
+ {
+ double dir[3] = {0,0,1};
+ sfc_pnts = MeshLib::MeshSurfaceExtraction::getSurfaceNodes(*mesh, dir);
+ }
+ size_t nPoints (sfc_pnts.size());
+ for (unsigned i=0; i<nPoints; ++i)
+ {
+ GeoLib::PointWithID* pnt (sfc_pnts[i]);
+ (*pnt)[2] = 0;
+ }
+
+ return new GeoLib::Grid<GeoLib::PointWithID>(sfc_pnts.begin(), sfc_pnts.end());
+}
unsigned getIndexInPntVec(GeoLib::Point const*const pnt, std::vector<GeoLib::Point*> const*const points)
{
@@ -108,9 +215,8 @@ void GeoMapper::advancedMapOnMesh(const MeshLib::Mesh* mesh, const std::string &
const std::vector<GeoLib::Point*> *points (this->_geo_objects.getPointVec(this->_geo_name));
const std::vector<GeoLib::Polyline*> *org_lines (this->_geo_objects.getPolylineVec(this->_geo_name));
- double eps = sqrt(std::numeric_limits<float>::epsilon());
- GeoLib::AABB<GeoLib::Point> aabb(points->begin(), points->end());
- eps *= sqrt(MathLib::sqrDist (&(aabb.getMinPoint()),&(aabb.getMaxPoint())));
+ const GeoLib::AABB<GeoLib::Point> aabb(points->begin(), points->end());
+ const double eps (aabb.getEpsFromBoundingBox());
// copy geometry (and set z=0 for all points)
unsigned nGeoPoints ( points->size() );
@@ -282,113 +388,7 @@ double GeoMapper::getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &line
return max_segment_length;
}
-void GeoMapper::mapData()
-{
- const std::vector<GeoLib::Point*> *points (this->_geo_objects.getPointVec(this->_geo_name));
- GeoLib::Station* stn_test = dynamic_cast<GeoLib::Station*>((*points)[0]);
- bool is_borehole(false);
- if (stn_test != nullptr && static_cast<GeoLib::StationBorehole*>((*points)[0])->type() == GeoLib::Station::StationType::BOREHOLE)
- is_borehole = true;
-
- double min_val(0), max_val(0);
- if (_mesh)
- {
- GeoLib::AABB<GeoLib::Point> bounding_box (_mesh->getNodes().begin(), _mesh->getNodes().end());
- min_val = bounding_box.getMinPoint()[2];
- max_val = bounding_box.getMaxPoint()[2];
- }
- size_t nPoints (points->size());
-
- if (!is_borehole)
- {
- for (unsigned j=0; j<nPoints; ++j)
- {
- GeoLib::Point* pnt ((*points)[j]);
- (*pnt)[2] = (_grid) ? this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val)
- : this->getDemElevation((*pnt)[0],(*pnt)[1]);
- }
- }
- else
- {
- for (unsigned j=0; j<nPoints; ++j)
- {
- GeoLib::Point* pnt ((*points)[j]);
- double offset = (_grid) ? (this->getMeshElevation((*pnt)[0],(*pnt)[1], min_val, max_val) - (*pnt)[2])
- : (this->getDemElevation((*pnt)[0],(*pnt)[1]) - (*pnt)[2]);
-
- GeoLib::StationBorehole* borehole = static_cast<GeoLib::StationBorehole*>(pnt);
- const std::vector<GeoLib::Point*> layers = borehole->getProfile();
- size_t nLayers = layers.size();
- for (unsigned k=0; k<nLayers; ++k)
- {
- GeoLib::Point* layer_pnt = layers[k];
- (*layer_pnt)[2] = (*layer_pnt)[2] + offset;
- }
- }
- }
-}
-
-float GeoMapper::getDemElevation(double x, double y) const
-{
- const double origin_x(_raster->getOrigin()[0]);
- const double origin_y(_raster->getOrigin()[1]);
- const double cellsize(_raster->getRasterPixelDistance());
- const std::size_t width(_raster->getNCols());
- const std::size_t height(_raster->getNRows());
- if ((x<origin_x) || (x>origin_x+(width*cellsize)) || (y<origin_y) || (y>origin_y+(height*cellsize)))
- return 0;
-
- const unsigned x_index = static_cast<unsigned>((x-origin_x)/cellsize);
- const unsigned y_index = static_cast<unsigned>((y-origin_y)/cellsize);
-
- return static_cast<float>(*(_raster->begin() + (y_index*width+x_index)));
-}
-
-double GeoMapper::getMeshElevation(double x, double y, double min_val, double max_val) const
-{
- double coords[3] = {x,y,0};
- const GeoLib::PointWithID* pnt = _grid->getNearestPoint(coords);
- const std::vector<MeshLib::Element*> elements (_mesh->getNode(pnt->getID())->getElements());
- GeoLib::Point* intersection (nullptr);
-
- for (std::size_t i=0; i<elements.size(); ++i)
- {
- if (intersection==nullptr)
- intersection=GeoLib::triangleLineIntersection(*elements[i]->getNode(0), *elements[i]->getNode(1), *elements[i]->getNode(2), GeoLib::Point(x,y,max_val), GeoLib::Point(x,y,min_val));
- if (intersection==nullptr && elements[i]->getGeomType() == MeshElemType::QUAD)
- intersection=GeoLib::triangleLineIntersection(*elements[i]->getNode(0), *elements[i]->getNode(2), *elements[i]->getNode(3), GeoLib::Point(x,y,max_val), GeoLib::Point(x,y,min_val));
- }
- if (intersection)
- return (*intersection)[2];
- // if something goes wrong, simply take the elevation of the nearest mesh node
- return (*(_mesh->getNode(pnt->getID())))[2];
-}
-
-GeoLib::Grid<GeoLib::PointWithID>* GeoMapper::getFlatGrid(MeshLib::Mesh const*const mesh, std::vector<GeoLib::PointWithID*> sfc_pnts) const
-{
- if (mesh->getDimension()<3) //much faster
- {
- size_t nNodes (mesh->getNNodes());
- sfc_pnts.resize(nNodes);
- const std::vector<MeshLib::Node*> nodes (mesh->getNodes());
- for (unsigned i(0); i<nNodes; ++i)
- sfc_pnts[i] = new GeoLib::PointWithID(nodes[i]->getCoords(), nodes[i]->getID());
- }
- else
- {
- double dir[3] = {0,0,1};
- sfc_pnts = MeshLib::MeshSurfaceExtraction::getSurfaceNodes(*mesh, dir);
- }
- size_t nPoints (sfc_pnts.size());
- for (unsigned i=0; i<nPoints; ++i)
- {
- GeoLib::PointWithID* pnt (sfc_pnts[i]);
- (*pnt)[2] = 0;
- }
-
- return new GeoLib::Grid<GeoLib::PointWithID>(sfc_pnts.begin(), sfc_pnts.end());
-}
--
GitLab