diff --git a/Applications/Utils/GeoTools/RemoveUnusedPoints.cpp b/Applications/Utils/GeoTools/RemoveUnusedPoints.cpp
index 4399ed00f1366d68cb83119993a2333a3dd49d90..8de26bf2c5eb313a8dd879441f47586f92cc031f 100644
--- a/Applications/Utils/GeoTools/RemoveUnusedPoints.cpp
+++ b/Applications/Utils/GeoTools/RemoveUnusedPoints.cpp
@@ -85,10 +85,7 @@ int main(int argc, char* argv[])
     {
         for (auto const* polyline : *polylines)
         {
-            for (std::size_t i = 0; i < polyline->getNumberOfPoints(); ++i)
-            {
-                used_points[polyline->getPointID(i)] = true;
-            }
+            GeoLib::markUsedPoints(*polyline, used_points);
         }
     }
     if (surfaces)
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index d79671993339927013e2a9d078fcc66ce74cee20..c41414bc20164cebdf634d89a2665b4409d55a30 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -494,6 +494,21 @@ void resetPointIDs(Polyline& polyline, std::vector<std::size_t> const& mapping)
     }
 }
 
+void markUsedPoints(Polyline const& polyline, std::vector<bool>& used_points)
+{
+    if (polyline.getPointsVec().size() != used_points.size())
+    {
+        OGS_FATAL(
+            "internal error in markUsedPoints(): polyline based on point "
+            "vector of size {}, given used_points has size {}",
+            polyline.getPointsVec().size(), used_points.size());
+    }
+    for (std::size_t i = 0; i < polyline.getNumberOfPoints(); ++i)
+    {
+        used_points[polyline.getPointID(i)] = true;
+    }
+}
+
 bool containsEdge(const Polyline& ply, std::size_t id0, std::size_t id1)
 {
     if (id0 == id1)
diff --git a/GeoLib/Polyline.h b/GeoLib/Polyline.h
index 0f50261e809a712a8531bf67c29b862003f0d524..60790b717f75e4493a7eae7e9836ba9ebcb73907 100644
--- a/GeoLib/Polyline.h
+++ b/GeoLib/Polyline.h
@@ -230,6 +230,9 @@ bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1);
 /// Resets the point IDs of the polyline corresponding to the mapping.
 void resetPointIDs(Polyline& polyline, std::vector<std::size_t> const& mapping);
 
+/// Resets the point IDs of the polyline corresponding to the mapping.
+void markUsedPoints(Polyline const& polyline, std::vector<bool>& used_points);
+
 /**
  * comparison operator
  * \param lhs first polyline