Skip to content

File covers.cpp

File List > algorithm > covers.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/covers.h"
#include "SFCGAL/Geometry.h"
#include "SFCGAL/Kernel.h"
#include "SFCGAL/algorithm/intersection.h"
#include "SFCGAL/algorithm/intersects.h"
#include "SFCGAL/detail/GeometrySet.h"
#include "SFCGAL/detail/TypeForDimension.h"

#include <CGAL/box_intersection_d.h>

using namespace SFCGAL::detail;

namespace SFCGAL::algorithm {

auto
covers(const PrimitiveHandle<3> & /*unused*/,
       const PrimitiveHandle<3> & /*unused*/) -> bool
{
  return false;
}

template <int Dim>
auto
segmentsLength(const GeometrySet<Dim> &gs) -> double
{
  double result = 0.0;

  for (auto it = gs.segments().begin(); it != gs.segments().end(); ++it) {
    result = result + sqrt(CGAL::to_double(it->primitive().squared_length()));
  }

  return result;
}

auto
solidsVolume(const GeometrySet<3> &gs, bool planarSurface = false) -> double
{
  double result = 0.0;

  for (const auto &it : gs.volumes()) {
    // TODO : we use areas of surfaces here instead of volumes
    const MarkedPolyhedron &poly = it.primitive();

    if (poly.is_closed() && planarSurface) {
      continue;
    }

    if (!poly.is_closed() && !planarSurface) {
      continue;
    }

    BOOST_ASSERT(poly.is_pure_triangle());

    CGAL::Point_3<Kernel> p1;
    CGAL::Point_3<Kernel> p2;
    CGAL::Point_3<Kernel> p3;

    for (MarkedPolyhedron::Facet_const_iterator fit = poly.facets_begin();
         fit != poly.facets_end(); ++fit) {
      MarkedPolyhedron::Halfedge_around_facet_const_circulator cit =
          fit->facet_begin();
      p1 = cit->vertex()->point();
      cit++;
      p2 = cit->vertex()->point();
      cit++;
      p3 = cit->vertex()->point();
      CGAL::Triangle_3<Kernel> const tri(p1, p2, p3);
      result = result + sqrt(CGAL::to_double(tri.squared_area()));
    }
  }

  return result;
}

auto
surfacesArea(const GeometrySet<2> &gs) -> double
{
  Kernel::FT result = 0.0;

  for (const auto &it : gs.surfaces()) {
    const CGAL::Polygon_with_holes_2<Kernel> &polygon = it.primitive();
    result = result + CGAL::abs(polygon.outer_boundary().area());

    for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) {
      result = result - CGAL::abs(hit->area());
    }
  }

  return CGAL::to_double(result);
}

auto
surfacesArea(const GeometrySet<3> &gs) -> double
{
  double result = 0.0;

  if (gs.surfaces().empty() && !gs.volumes().empty()) {
    result = solidsVolume(gs, /* planarSurface = */ true);
  }

  for (const auto &it : gs.surfaces()) {
    result = result + sqrt(CGAL::to_double(it.primitive().squared_area()));
  }

  return result;
}

auto
solidsVolume(const GeometrySet<2> & /*unused*/) -> double
{
  return 0.0;
}

template <int Dim>
auto
equalLength(const GeometrySet<Dim> &a, const GeometrySet<Dim> &b, int dim)
    -> bool
{
  // compare 'length' of primitives in A with 'length' of primitives in B
  // 'length' is :
  // - number of elements for points
  // - length for segments
  // - area for surfaces
  // - should be volume for volumes. We use area here

  double const tol = 1e-9;

  switch (dim) {
  case 0: {
    if (a.points().size() != b.points().size()) {
      return false;
    }
  } break;

  case 1: {

    //
    // Compare lengths
    double const lengthA = segmentsLength(a);
    double const lengthB = segmentsLength(b);
    double const cmp     = (lengthA - lengthB) * (lengthA - lengthB);

    if (cmp > tol) {
      return false;
    }
  } break;

  case 2: {
    //
    // Compare areas
    double const areaA = surfacesArea(a);
    double const areaB = surfacesArea(b);
    double const cmp   = (areaA - areaB) * (areaA - areaB);

    if (cmp > tol) {
      return false;
    }
  } break;

  case 3: {
    // Compare volumes
    double const volA = solidsVolume(a);
    double const volB = solidsVolume(b);
    double const cmp  = (volA - volB) * (volA - volB);

    if (cmp > tol) {
      return false;
    }
  } break;
  }

  return true;
}

template <int Dim>
auto
covers(const GeometrySet<Dim> &a, const GeometrySet<Dim> &b) -> bool
{
  int const dimA = a.dimension();
  int const dimB = b.dimension();

  if (dimA == -1 || dimB == -1) {
    return false;
  }

  if (dimB > dimA) {
    return false;
  }

  //
  // This is a very naive (not efficient) implementation of covers() !
  //
  // covers(A,B) <=> A inter B == B
  // '==' is here implemented with comparison of length, area and volumes
  // TODO use only predicates if possible
  GeometrySet<Dim> inter;
  algorithm::intersection(a, b, inter);

  if (b.hasPoints() && !equalLength(b, inter, 0)) {
    return false;
  }

  if (b.hasSegments() && !equalLength(b, inter, 1)) {
    return false;
  }

  if (b.hasSurfaces() && !equalLength(b, inter, 2)) {
    return false;
  }

  if (b.hasVolumes() && !equalLength(b, inter, 3)) {
    return false;
  }

  return true;
}

template bool
covers<2>(const GeometrySet<2> &a, const GeometrySet<2> &b);
template bool
covers<3>(const GeometrySet<3> &a, const GeometrySet<3> &b);

auto
covers(const Geometry &ga, const Geometry &gb) -> bool
{
  if (ga.isEmpty() || gb.isEmpty()) {
    return false;
  }

  GeometrySet<2> const gsa(ga);
  GeometrySet<2> const gsb(gb);

  return covers(gsa, gsb);
}

auto
covers3D(const Geometry &ga, const Geometry &gb) -> bool
{
  if (ga.isEmpty() || gb.isEmpty()) {
    return false;
  }

  GeometrySet<3> const gsa(ga);
  GeometrySet<3> const gsb(gb);

  return covers(gsa, gsb);
}
} // namespace SFCGAL::algorithm