File isValid.cpp
File List > algorithm > isValid.cpp
Go to the documentation of this file
// Copyright (c) 2012-2013, IGN France.
// Copyright (c) 2012-2022, Oslandia.
// SPDX-License-Identifier: LGPL-2.0-or-later
#include "SFCGAL/algorithm/isValid.h"
#include "SFCGAL/GeometryCollection.h"
#include "SFCGAL/LineString.h"
#include "SFCGAL/MultiLineString.h"
#include "SFCGAL/MultiPoint.h"
#include "SFCGAL/MultiPolygon.h"
#include "SFCGAL/MultiSolid.h"
#include "SFCGAL/Polygon.h"
#include "SFCGAL/PolyhedralSurface.h"
#include "SFCGAL/Triangle.h"
#include "SFCGAL/TriangulatedSurface.h"
#include "SFCGAL/Exception.h"
#include "SFCGAL/Kernel.h"
#include "SFCGAL/algorithm/connection.h"
#include "SFCGAL/algorithm/distance.h"
#include "SFCGAL/algorithm/distance3d.h"
#include "SFCGAL/algorithm/intersection.h"
#include "SFCGAL/algorithm/intersects.h"
#include "SFCGAL/algorithm/length.h"
#include "SFCGAL/algorithm/normal.h"
#include "SFCGAL/algorithm/orientation.h"
#include "SFCGAL/algorithm/plane.h"
#include "SFCGAL/detail/ForceValidityVisitor.h"
#include "SFCGAL/detail/GetPointsVisitor.h"
#include "SFCGAL/detail/algorithm/coversPoints.h"
#include "SFCGAL/detail/tools/Log.h"
#include <boost/format.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/undirected_dfs.hpp>
#include <boost/graph/visitors.hpp>
using namespace SFCGAL::detail::algorithm;
namespace SFCGAL {
void
SFCGAL_ASSERT_GEOMETRY_VALIDITY_(const Geometry &g, const std::string &ctxt)
{
if (!(g).hasValidityFlag()) {
const Validity sfcgalAssertGeometryValidity = algorithm::isValid(g);
if (!sfcgalAssertGeometryValidity) {
throw GeometryInvalidityException(
(boost::format(ctxt + "%s is invalid : %s : %s") %
(g).geometryType() % sfcgalAssertGeometryValidity.reason() %
(g).asText())
.str());
}
}
}
void
SFCGAL_ASSERT_GEOMETRY_VALIDITY(const Geometry &g)
{
SFCGAL_ASSERT_GEOMETRY_VALIDITY_(g, "");
}
void
SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(const Geometry &g)
{
if (!(g).hasValidityFlag()) {
using namespace SFCGAL;
if ((g).is3D()) {
std::unique_ptr<SFCGAL::Geometry> sfcgalAssertGeometryValidityClone(
(g).clone());
algorithm::force2D(*sfcgalAssertGeometryValidityClone);
SFCGAL_ASSERT_GEOMETRY_VALIDITY_((*sfcgalAssertGeometryValidityClone),
"When converting to 2D - ");
} else {
SFCGAL_ASSERT_GEOMETRY_VALIDITY(g);
}
}
}
void
SFCGAL_ASSERT_GEOMETRY_VALIDITY_3D(const Geometry &g)
{
if (!(g).hasValidityFlag()) {
using namespace SFCGAL;
if (!(g).is3D()) {
std::unique_ptr<Geometry> sfcgalAssertGeometryValidityClone((g).clone());
algorithm::force3D(*sfcgalAssertGeometryValidityClone);
SFCGAL_ASSERT_GEOMETRY_VALIDITY_((*sfcgalAssertGeometryValidityClone),
"When converting to 3D - ");
} else {
SFCGAL_ASSERT_GEOMETRY_VALIDITY(g);
}
}
}
void
SFCGAL_ASSERT_GEOMETRY_VALIDITY_ON_PLANE(const Geometry & /*g*/)
{
throw NotImplementedException(
"validation on geometry projected on arbitrary plane is not implemented");
}
namespace algorithm {
// to detect unconnected interior in polygon
struct LoopDetector : public boost::dfs_visitor<> {
LoopDetector(bool &hasLoop) : _hasLoop(hasLoop) {}
template <class Edge, class Graph>
void
back_edge(Edge /*unused*/, const Graph & /*unused*/)
{
_hasLoop = true;
}
private:
bool &_hasLoop;
};
// const Validity isValid( const Coordinate & p)
//{
// BOOST_ASSERT( !p.isEmpty() );
// if ( !CGAL::is_finite(p.x()) || CGAL::is_finite(p.y()) ) return
// Validity::invalid("NaN coordinate");
// //if ( p.x().is_inf() || p.y().is_inf() ) return
// Validity::invalid("infinite coordinate"); return Validity::valid();
//}
//
auto
isValid(const Point &p) -> const Validity
{
if (p.isEmpty()) {
return Validity::valid();
}
(void)p;
// return( isValid(p.coordinate() ) );
return Validity::valid();
}
auto
isValid(const LineString &l, const double &toleranceAbs) -> const Validity
{
if (l.isEmpty()) {
return Validity::valid();
}
// const size_t numPoints = l.numPoints();
// for ( size_t p=0; p!=numPoints; ++p) {
// const Validity v = isValid(l.pointN(p));
// if (!v) return Validity::invalid( ( boost::format("Point %d is
// invalid: %s") % p % v.reason() ).str() );
// }
return length3D(l) > toleranceAbs ? Validity::valid()
: Validity::invalid("no length");
}
auto
isValid(const Polygon &p, const double &toleranceAbs) -> const Validity
{
if (p.isEmpty()) {
return Validity::valid();
}
// Closed simple rings
const size_t numRings = p.numRings();
for (size_t r = 0; r != numRings; ++r) {
if (p.ringN(r).numPoints() < 4) {
return Validity::invalid(
(boost::format("not enough points in ring %d") % r).str());
}
// const Validity v = isValid( p.ringN(r) );
// if (!v) return Validity::invalid((boost::format("ring %d is
// invalid: %s") % r % v.reason()).str() );
const double distanceToClose =
p.is3D() ? distancePointPoint3D(p.ringN(r).startPoint(),
p.ringN(r).endPoint())
: distancePointPoint(p.ringN(r).startPoint(),
p.ringN(r).endPoint());
if (distanceToClose > 0) {
return Validity::invalid(
(boost::format("ring %d is not closed") % r).str());
}
if (p.is3D() ? selfIntersects3D(p.ringN(r)) : selfIntersects(p.ringN(r))) {
return Validity::invalid(
(boost::format("ring %d self intersects") % r).str());
}
}
// When polygon degenerates to a single point, it is not trapped by the rest
// of the code, so we check that here
for (size_t r = 0; r != numRings; ++r) {
const LineString &ring = p.ringN(r);
const Point &start = ring.startPoint();
size_t i = 0;
for (; i < ring.numPoints() && start == ring.pointN(i); i++) {
; // noop
}
if (i == ring.numPoints()) {
return Validity::invalid(
(boost::format("ring %d degenerated to a point") % r).str());
}
}
// Orientation in 2D
if (!p.is3D()) {
// Opposit orientation for interior and exterior rings
const bool extCCWO = isCounterClockWiseOriented(p.exteriorRing());
for (std::size_t r = 0; r < p.numInteriorRings(); ++r) {
if (extCCWO == isCounterClockWiseOriented(p.interiorRingN(r))) {
return Validity::invalid(
(boost::format("exterior ring and interior ring %d have the same "
"orientation") %
r)
.str());
}
}
}
// Orientation in 3D
else {
// Polygone must be planar (all points in the same plane)
if (!isPlane3D<Kernel>(p, toleranceAbs)) {
return Validity::invalid("points don't lie in the same plane");
}
// interior rings must be oriented opposit to exterior;
if (p.hasInteriorRings()) {
const CGAL::Vector_3<Kernel> nExt = normal3D<Kernel>(p.exteriorRing());
for (std::size_t r = 0; r < p.numInteriorRings(); ++r) {
const CGAL::Vector_3<Kernel> nInt =
normal3D<Kernel>(p.interiorRingN(r));
if (nExt * nInt > 0) {
return Validity::invalid(
(boost::format("interior ring %d is oriented in the same "
"direction as exterior ring") %
r)
.str());
}
}
}
}
// Rings must not share more than one point (no intersection)
{
using Edge = std::pair<int, int>;
std::vector<Edge> touchingRings;
for (size_t ri = 0; ri < numRings;
++ri) { // no need for numRings-1, the next loop won't be entered for
// the last ring
for (size_t rj = ri + 1; rj < numRings; ++rj) {
std::unique_ptr<Geometry> inter =
p.is3D() ? intersection3D(p.ringN(ri), p.ringN(rj))
: intersection(p.ringN(ri), p.ringN(rj));
if (!inter->isEmpty() && !inter->is<Point>()) {
return Validity::invalid(
(boost::format("intersection between ring %d and %d") % ri % rj)
.str());
}
if (!inter->isEmpty() && inter->is<Point>()) {
touchingRings.emplace_back(ri, rj);
}
}
}
{
using namespace boost;
using Graph = adjacency_list<vecS, vecS, undirectedS, no_property,
property<edge_color_t, default_color_type>>;
using vertex_t = graph_traits<Graph>::vertex_descriptor;
Graph g(touchingRings.begin(), touchingRings.end(), numRings);
bool hasLoop = false;
LoopDetector const vis(hasLoop);
undirected_dfs(g, root_vertex(vertex_t(0))
.visitor(vis)
.edge_color_map(get(edge_color, g)));
if (hasLoop) {
return Validity::invalid("interior is not connected");
}
}
}
if (p.hasInteriorRings()) {
// Interior rings must be interior to exterior ring
for (size_t r = 0; r < p.numInteriorRings();
++r) { // no need for numRings-1, the next loop won't be entered for
// the last ring
if (p.is3D() ? !coversPoints3D(Polygon(p.exteriorRing()),
Polygon(p.interiorRingN(r)))
: !coversPoints(Polygon(p.exteriorRing()),
Polygon(p.interiorRingN(r)))) {
return Validity::invalid(
(boost::format("exterior ring doesn't cover interior ring %d") % r)
.str());
}
}
// Interior ring must not cover one another
for (size_t ri = 0; ri < p.numInteriorRings();
++ri) { // no need for numRings-1, the next loop won't be entered for
// the last ring
for (size_t rj = ri + 1; rj < p.numInteriorRings(); ++rj) {
if (p.is3D() ? coversPoints3D(Polygon(p.interiorRingN(ri)),
Polygon(p.interiorRingN(rj)))
: coversPoints(Polygon(p.interiorRingN(ri)),
Polygon(p.interiorRingN(rj)))) {
return Validity::invalid(
(boost::format("interior ring %d covers interior ring %d") % ri %
rj)
.str());
}
}
}
}
return Validity::valid();
}
auto
isValid(const Triangle &t, const double &toleranceAbs) -> const Validity
{
return isValid(t.toPolygon(), toleranceAbs);
}
auto
isValid(const MultiLineString &ml, const double &toleranceAbs) -> const Validity
{
if (ml.isEmpty()) {
return Validity::valid();
}
const size_t numLineString = ml.numGeometries();
for (size_t l = 0; l != numLineString; ++l) {
Validity const v = isValid(ml.lineStringN(l), toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("LineString %d is invalid: %s") % l % v.reason())
.str());
}
}
return Validity::valid();
}
auto
isValid(const MultiPolygon &mp, const double &toleranceAbs) -> const Validity
{
if (mp.isEmpty()) {
return Validity::valid();
}
const size_t numPolygons = mp.numGeometries();
for (size_t p = 0; p != numPolygons; ++p) {
Validity const v = isValid(mp.polygonN(p), toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("Polygon %d is invalid: %s") % p % v.reason()).str());
}
}
for (size_t pi = 0; pi != numPolygons; ++pi) {
for (size_t pj = pi + 1; pj < numPolygons; ++pj) {
std::unique_ptr<Geometry> inter =
mp.is3D() ? intersection3D(mp.polygonN(pi), mp.polygonN(pj))
: intersection(mp.polygonN(pi), mp.polygonN(pj));
// intersection can be empty, a point, or a set of points
if (!inter->isEmpty() && inter->dimension() != 0) {
return Validity::invalid(
(boost::format("intersection between Polygon %d and %d") % pi % pj)
.str());
}
}
}
return Validity::valid();
}
auto
isValid(const GeometryCollection &gc, const double &toleranceAbs)
-> const Validity
{
if (gc.isEmpty()) {
return Validity::valid();
}
const size_t numGeom = gc.numGeometries();
for (size_t g = 0; g != numGeom; ++g) {
Validity const v = isValid(gc.geometryN(g), toleranceAbs);
if (!v) {
return Validity::invalid((boost::format("%s %d is invalid: %s") %
gc.geometryN(g).geometryType() % g % v.reason())
.str());
}
}
return Validity::valid();
}
auto
isValid(const TriangulatedSurface &tin, const SurfaceGraph &graph,
const double &toleranceAbs) -> const Validity
{
if (tin.isEmpty()) {
return Validity::valid();
}
size_t const numTriangles = tin.numTriangles();
for (size_t t = 0; t != numTriangles; ++t) {
Validity const v = isValid(tin.triangleN(t), toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("Triangle %d is invalid: %s") % t % v.reason()).str());
}
}
if (!isConnected(graph)) {
return Validity::invalid("not connected");
}
if (tin.is3D() ? selfIntersects3D(tin, graph) : selfIntersects(tin, graph)) {
return Validity::invalid("self intersects");
}
return Validity::valid();
}
auto
isValid(const TriangulatedSurface &tin, const double &toleranceAbs)
-> const Validity
{
if (tin.isEmpty()) {
return Validity::valid();
}
const SurfaceGraph graph(tin);
return graph.isValid() ? isValid(tin, graph, toleranceAbs) : graph.isValid();
}
auto
isValid(const PolyhedralSurface &s, const SurfaceGraph &graph,
const double &toleranceAbs) -> const Validity
{
if (s.isEmpty()) {
return Validity::valid();
}
size_t const numPolygons = s.numPolygons();
for (size_t p = 0; p != numPolygons; ++p) {
Validity const v = isValid(s.polygonN(p), toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("Polygon %d is invalid: %s") % p % v.reason()).str());
}
}
if (!isConnected(graph)) {
return Validity::invalid("not connected");
}
if (s.is3D() ? selfIntersects3D(s, graph) : selfIntersects(s, graph)) {
return Validity::invalid("self intersects");
}
return Validity::valid();
}
auto
isValid(const PolyhedralSurface &s, const double &toleranceAbs)
-> const Validity
{
if (s.isEmpty()) {
return Validity::valid();
}
const SurfaceGraph graph(s);
return graph.isValid() ? isValid(s, graph, toleranceAbs) : graph.isValid();
}
auto
isValid(const Solid &solid, const double &toleranceAbs) -> const Validity
{
if (solid.isEmpty()) {
return Validity::valid();
}
const size_t numShells = solid.numShells();
for (size_t s = 0; s != numShells; ++s) {
const SurfaceGraph graph(solid.shellN(s));
Validity const v = isValid(solid.shellN(s), graph, toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("PolyhedralSurface (shell) %d is invalid: %s") % s %
v.reason())
.str());
}
if (!isClosed(graph)) {
return Validity::invalid(
(boost::format("PolyhedralSurface (shell) %d is not closed") % s)
.str());
}
}
if (solid.numInteriorShells() != 0U) {
BOOST_THROW_EXCEPTION(
Exception("function is not fully implemented (orientation, covering "
"and intersections of interior shells missing"));
}
return Validity::valid();
}
auto
isValid(const MultiSolid &ms, const double &toleranceAbs) -> const Validity
{
if (ms.isEmpty()) {
return Validity::valid();
}
const size_t numMultiSolid = ms.numGeometries();
for (size_t s = 0; s != numMultiSolid; ++s) {
Validity const v = isValid(ms.solidN(s), toleranceAbs);
if (!v) {
return Validity::invalid(
(boost::format("Solid %d is invalid: %s") % s % v.reason()).str());
}
}
return Validity::valid();
}
auto
isValid(const Geometry &g, const double &toleranceAbs) -> const Validity
{
switch (g.geometryTypeId()) {
case TYPE_POINT:
return isValid(g.as<Point>());
case TYPE_LINESTRING:
return isValid(g.as<LineString>(), toleranceAbs);
case TYPE_POLYGON:
return isValid(g.as<Polygon>(), toleranceAbs);
case TYPE_TRIANGLE:
return isValid(g.as<Triangle>(), toleranceAbs);
case TYPE_SOLID:
return isValid(g.as<Solid>(), toleranceAbs);
case TYPE_MULTIPOINT:
return Validity::valid();
case TYPE_MULTILINESTRING:
return isValid(g.as<MultiLineString>(), toleranceAbs);
case TYPE_MULTIPOLYGON:
return isValid(g.as<MultiPolygon>(), toleranceAbs);
case TYPE_MULTISOLID:
return isValid(g.as<MultiSolid>(), toleranceAbs);
case TYPE_GEOMETRYCOLLECTION:
return isValid(g.as<GeometryCollection>(), toleranceAbs);
case TYPE_TRIANGULATEDSURFACE:
return isValid(g.as<TriangulatedSurface>(), toleranceAbs);
case TYPE_POLYHEDRALSURFACE:
return isValid(g.as<PolyhedralSurface>(), toleranceAbs);
}
BOOST_THROW_EXCEPTION(Exception(
(boost::format("isValid( %s ) is not defined") % g.geometryType())
.str()));
return Validity::invalid(
(boost::format("isValid( %s ) is not defined") % g.geometryType())
.str()); // to avoid warning
}
void
propagateValidityFlag(Geometry &g, bool valid)
{
detail::ForceValidityVisitor v(valid);
g.accept(v);
}
} // namespace algorithm
} // namespace SFCGAL