diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index 33412a2f8b66b3f2516dd29986da9d5fc9cc4899..c9c9aa03ffd38e904324894d8deb235c27787723 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -217,146 +217,6 @@ double GeoMapper::getMeshElevation(
return (*(_surface_mesh->getNode(pnt->getID())))[2];
}
-std::unique_ptr<std::vector<GeoLib::Polyline*>> copyPolylinesVector(
- std::vector<GeoLib::Polyline*> const& polylines,
- std::vector<GeoLib::Point*> const& points)
-{
- std::size_t nLines = polylines.size();
- auto new_lines = std::unique_ptr<std::vector<GeoLib::Polyline*>>(
- new std::vector<GeoLib::Polyline*>(nLines));
-
- for (std::size_t i=0; i<nLines; ++i)
- {
- (*new_lines)[i] = new GeoLib::Polyline(points);
- std::size_t nLinePnts (polylines[i]->getNumberOfPoints());
- for (std::size_t j=0; j<nLinePnts; ++j)
- (*new_lines)[i]->addPoint(polylines[i]->getPointID(j));
- }
- return new_lines;
-}
-
-
-void GeoMapper::advancedMapOnMesh(
- MeshLib::Mesh const* mesh, std::string const& new_geo_name)
-{
- const std::vector<GeoLib::Point*> *points(_geo_objects.getPointVec(_geo_name));
- const std::vector<GeoLib::Polyline*> *org_lines(_geo_objects.getPolylineVec(_geo_name));
-
- const GeoLib::AABB aabb(points->begin(), points->end());
- const double eps = sqrt(std::numeric_limits<float>::epsilon()) *
- sqrt(MathLib::sqrDist(aabb.getMinPoint(),aabb.getMaxPoint())) ;
-
- // copy geometry (and set z=0 for all points)
- auto new_points = std::unique_ptr<std::vector<GeoLib::Point*>>(
- new std::vector<GeoLib::Point*>);
- new_points->reserve(points->size());
- std::transform(points->cbegin(), points->cend(), std::back_inserter(*new_points),
- [](GeoLib::Point* p) { return new GeoLib::Point((*p)[0],(*p)[1],0.0); });
-
- auto new_lines = copyPolylinesVector(*_geo_objects.getPolylineVec(_geo_name),
- *new_points);
-
- GeoLib::Grid<GeoLib::Point> grid(new_points->begin(), new_points->end());
- double max_segment_length(getMaxSegmentLength(*new_lines));
- // squared so it can be compared to the squared distances calculated later
- max_segment_length *= max_segment_length;
-
- const unsigned nMeshNodes ( mesh->getNumberOfNodes() );
- // index of closest geo point for each mesh node in (x,y)-plane
- std::vector<int> closest_geo_point(nMeshNodes, -1);
- // distance between geo points and mesh nodes in (x,y)-plane
- std::vector<double> dist(nMeshNodes);
- auto zero_coords = GeoLib::Point{}; // All coordinates zero.
- for (std::size_t i=0; i<nMeshNodes; ++i) {
- zero_coords[0] = (*mesh->getNode(i))[0];
- zero_coords[1] = (*mesh->getNode(i))[1];
- GeoLib::Point* pnt = grid.getNearestPoint(zero_coords);
- dist[i] = MathLib::sqrDist(*pnt, zero_coords);
- closest_geo_point[i] = (dist[i]<=max_segment_length) ? pnt->getID() : -1;
- }
-
- // store for each point the line segment to which it was added.
- const std::size_t nLines(new_lines->size());
- std::vector< std::vector<unsigned> > line_segment_map(nLines);
- for (std::size_t i=0; i<nLines; ++i)
- {
- line_segment_map[i] = std::vector<unsigned>((*new_lines)[i]->getNumberOfPoints(),0);
- std::iota(line_segment_map[i].begin(), line_segment_map[i].end(), 0);
- }
-
- for (std::size_t i=0; i<nMeshNodes; ++i)
- {
- // if mesh node too far away or exactly at point position
- if (closest_geo_point[i] == -1 || dist[i] < eps) continue;
-
- const MeshLib::Node* node (mesh->getNode(i));
- for (std::size_t l=0; l<nLines; ++l)
- {
- // find relevant polylines
- if (!(*org_lines)[l]->isPointIDInPolyline(closest_geo_point[i])) continue;
-
- // find point position of closest geo point in original polyline
- GeoLib::Polyline* ply ((*org_lines)[l]);
- std::size_t nLinePnts ( ply->getNumberOfPoints() );
- std::size_t node_index_in_ply (0);
- for (node_index_in_ply=0; node_index_in_ply<nLinePnts; ++node_index_in_ply)
- if (ply->getPoint(node_index_in_ply) == (*points)[closest_geo_point[i]])
- break;
- const GeoLib::Point* geo_point (ply->getPoint(node_index_in_ply));
-
- // check if line segments connected to closest geo point intersect connected elements of current node
- const std::vector<MeshLib::Element*>& elements (node->getElements());
- const std::size_t nElems = elements.size();
- for (std::size_t e=0; e<nElems; ++e)
- {
- const unsigned nEdges (elements[e]->getNumberOfEdges());
- unsigned intersection_count (0);
-
- for (unsigned n=0; n<nEdges; ++n)
- {
- if (intersection_count>1) break; //already two intersections
-
- const MeshLib::Element* line = elements[e]->getEdge(n);
- unsigned index_offset(0); // default: add to first line segment
- GeoLib::Point* intersection (nullptr);
- if (node_index_in_ply>0) // test line segment before closest point
- intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply-1));
- if (intersection == nullptr && node_index_in_ply<(nLinePnts-1)) // test line segment after closest point
- {
- intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply+1));
- index_offset = 1; // add to second segment
- }
- if (intersection)
- {
- intersection_count++;
- unsigned start_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset-1;})));
- unsigned end_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset;})));
- std::size_t pos = getPointPosInLine((*new_lines)[l], start_point_idx, end_point_idx, intersection, eps);
-
- if (pos)
- {
- const std::size_t pnt_pos (new_points->size());
- new_points->push_back(intersection);
- (*new_lines)[l]->insertPoint(pos, pnt_pos);
- line_segment_map[l].insert(line_segment_map[l].begin()+pos, node_index_in_ply+index_offset-1);
- }
- }
- }
- }
- }
- }
-
- _geo_objects.addPointVec(std::move(new_points), const_cast<std::string&>(new_geo_name));
- std::vector<std::size_t> pnt_id_map = this->_geo_objects.getPointVecObj(new_geo_name)->getIDMap();
- for (auto & new_line : *new_lines)
- new_line->updatePointIDs(pnt_id_map);
- _geo_objects.addPolylineVec(std::move(new_lines), new_geo_name);
-
- // map new geometry incl. additional point using the normal mapping method
- this->_geo_name = new_geo_name;
- this->mapOnMesh(mesh);
-}
-
/// Find the 2d-element within the \c elements that contains the given point \c p.
///
/// The algorithm projects every element of the elements vector and the point
@@ -705,66 +565,4 @@ void GeoMapper::advancedMapOnMesh(MeshLib::Mesh const& mesh)
}
}
-GeoLib::Point* GeoMapper::calcIntersection(MathLib::Point3d const*const p1, MathLib::Point3d const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const
-{
- const double x1 = (*p1)[0], x2 = (*p2)[0], x3 = (*q1)[0], x4 = (*q2)[0];
- const double y1 = (*p1)[1], y2 = (*p2)[1], y3 = (*q1)[1], y4 = (*q2)[1];
-
- const double det = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
- if (std::abs(det) < std::numeric_limits<double>::epsilon()) return nullptr;
-
- const double pre = (x1*y2 - y1*x2);
- const double post = (x3*y4 - y3*x4);
- const double x = ( pre * (x3 - x4) - (x1 - x2) * post ) / det;
- const double y = ( pre * (y3 - y4) - (y1 - y2) * post ) / det;
-
- // Check if the x and y coordinates are within both line segments
- if (isPntInBoundingBox(x1,y1,x2,y2,x,y) && isPntInBoundingBox(x3,y3,x4,y4,x,y))
- return new GeoLib::Point(x, y, 0);
- return nullptr;
-}
-
-unsigned GeoMapper::getPointPosInLine(GeoLib::Polyline const*const line, unsigned start, unsigned end, GeoLib::Point const*const point, double eps) const
-{
- const double* first = line->getPoint(start)->getCoords();
- const double* pnt = point->getCoords();
- const double max_dist = MathLib::sqrDist(first, pnt);
-
- // if point is at start or end of line segment
- if (max_dist<eps && MathLib::sqrDist(pnt, line->getPoint(end)->getCoords())) return 0;
-
- for (std::size_t i=start+1; i<end; ++i)
- {
- const double* current = (*line->getPoint(i)).getCoords();
- if (MathLib::sqrDist(pnt, current) < eps) return 0;
- if (MathLib::sqrDist(first, current) > max_dist) return i;
- }
- return end; // last point of segment
-}
-
-bool GeoMapper::isPntInBoundingBox(double ax, double ay, double bx, double by, double px, double py) const
-{
- if ( px < (std::min(ax, bx)-std::numeric_limits<double>::epsilon()) || px > (std::max(ax, bx)+std::numeric_limits<double>::epsilon()) ||
- py < (std::min(ay, by)-std::numeric_limits<double>::epsilon()) || py > (std::max(ay, by)+std::numeric_limits<double>::epsilon()) )
- return false;
- return true;
-}
-
-double GeoMapper::getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &lines) const
-{
- double max_segment_length (0);
- const std::size_t nPlys ( lines.size() );
- for (std::size_t i=0; i<nPlys; ++i)
- {
- const GeoLib::Polyline* line = lines[i];
- const std::size_t nPlyPoints = line->getNumberOfPoints();
- for (std::size_t j=1; j<nPlyPoints; ++j)
- {
- const double dist (line->getLength(j)-line->getLength(j-1));
- if (dist>max_segment_length)
- max_segment_length=dist;
- }
- }
- return max_segment_length;
-}
} // end namespace MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/GeoMapper.h b/MeshGeoToolsLib/GeoMapper.h
index 222b0e91f0aac2035fa63ea3ff49829a6e6b603d..09e4b0542a7ccfeb6c32bf9160a03a2061ba4476 100644
--- a/MeshGeoToolsLib/GeoMapper.h
+++ b/MeshGeoToolsLib/GeoMapper.h
@@ -56,16 +56,14 @@ public:
void mapToConstantValue(double value);
/**
- * Maps the geometry based on the given mesh file. A new geometry is created where all
- * geometric points are assigned an elevation value on the mesh surface. Additional points are
- * inserted whenever a polyline from the original geometry intersects a mesh node or the edge
- * of a mesh element.
+ * Maps the geometry based on the given mesh file. I.e., all geometric
+ * points are assigned an elevation value on the mesh surface. Additional
+ * points are inserted whenever a polyline from the original geometry
+ * intersects a mesh node or the edge of a mesh element.
* \param mesh Mesh the geometry is mapped on
- * \param new_geo_name Name of the newly created geometry
- * \result A new geometry with the given name is inserted into _geo_objects
+ * \result A new geometry with the given name is inserted into
+ * _geo_objects
*/
- void advancedMapOnMesh(const MeshLib::Mesh* mesh, const std::string &new_geo_name);
-
void advancedMapOnMesh(MeshLib::Mesh const& mesh);
private:
@@ -85,18 +83,6 @@ private:
/// Returns the elevation at Point (x,y) based on a raster
float getDemElevation(GeoLib::Point const& pnt) const;
- /// Calculates the intersection of two lines embedded in the xy-plane
- GeoLib::Point* calcIntersection(MathLib::Point3d const*const p1, MathLib::Point3d const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const;
-
- /// Returns the position of a point within a line-segment
- unsigned getPointPosInLine(GeoLib::Polyline const*const line, unsigned start, unsigned end, GeoLib::Point const*const point, double eps) const;
-
- /// Returns the maximum segment length in a polyline vector
- double getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &lines) const;
-
- /// Returns if a point p is within a bounding box defined by a and b
- bool isPntInBoundingBox(double ax, double ay, double bx, double by, double px, double py) const;
-
GeoLib::GEOObjects& _geo_objects;
std::string& _geo_name;