CGAL users discussion list

[Claus Volko <>]

Yes, I am doing region growing on a mesh. In the end I would like to check which region intersects with a given line (orthogonal to the xz plane).

I can send you code that is ready to be built. I've attached a file to this email.

The effect of this code should be a ply file in which all regions are black except the one that contains the given point. This last region will be red.

I've already tested the code a bit and it does seem to work, but my tests are not finished yet.

Am Do., 14. Sept. 2023 um 11:32 Uhr schrieb Sebastien Loriot <>:

It is hard to understand what you are trying to do.
If you want to get help, you should post a minimal example showing what
you want to achieve.

Yesterday I thought you were doing region growing on a point set but now
it seems to be on a mesh.

Best,

Sebastien.

On 9/14/23 08:30, Claus Volko ( via cgal-discuss
Mailing List) wrote:
> I've got it now. It is perhaps not the most elegant piece of code but it
> seems to work fine:
>
>    std::list<int> listRegions = {};
>
>    for (std::size_t i = 0; i < regions.size(); i++)
>    {
>        Region_growing::Primitive_and_region region = regions[i];
>        Line_3* l = new Line_3(Point_3(0, -1, 0), Point_3(0, +1, 0));
>        auto intersection =
> boost::get<Point_3>(&*CGAL::intersection(region.first, *l));
>
>        if (intersection && region.first.has_on(*intersection))
>        {
>            auto x11 =
> region.second.front()->halfedge()->vertex()->point().x();
>            auto z11 =
> region.second.front()->halfedge()->vertex()->point().z();
>            auto x12 =
> region.second.front()->halfedge()->next()->vertex()->point().x();
>            auto z12 =
> region.second.front()->halfedge()->next()->vertex()->point().z();
>            auto x13 =
> region.second.front()->halfedge()->next()->next()->vertex()->point().x();
>            auto z13 =
> region.second.front()->halfedge()->next()->next()->vertex()->point().z();
>            auto x1a = x11 < x12 ? x11 < x13 ? x11 : x13 : x12 < x13 ?
> x12 : x13;
>            auto z1a = z11 < z12 ? z11 < z13 ? z11 : z13 : z12 < z13 ?
> z12 : z13;
>            auto x1b = x11 > x12 ? x11 > x13 ? x11 : x13 : x12 > x13 ?
> x12 : x13;
>            auto z1b = z11 > z12 ? z11 > z13 ? z11 : z13 : z12 > z13 ?
> z12 : z13;
>
>            if (intersection->x() > x1a && intersection->x() < x1b &&
> intersection->z() > z1a && intersection->z() < z1b)
>            {
>                std::cout << "Region " << i << " has point on it : " <<
> region.first << " at " << *intersection << std::endl;
>                listRegions.push_back(i);
>            }
>        }
>    }
>
> Am Do., 14. Sept. 2023 um 07:28 Uhr schrieb Claus Volko
> < <mailto:>>:
>
>     I tried coloring the regions with intersections red. The result
>     looked strange. That's why I don't think the code I posted does its job.
>
>     Am Mi., 13. Sept. 2023 um 12:41 Uhr schrieb Claus Volko
>     < <mailto:>>:
>
>         I'm now trying this:
>
>            for (std::size_t i = 0; i < regions.size(); i++)
>            {
>                Region_growing::Primitive_and_region region = regions[i];
>                Line_3* l = new Line_3(Point_3(0, -1, 0), Point_3(0, +1, 0));
>                auto intersection =
>         boost::get<Point_3>(&*CGAL::intersection(region.first, *l));
>
>                if (intersection && region.first.has_on(*intersection))
>         std::cout << "Region has point on it: " << region.first << " at"
>         << *intersection << std::endl;
>            }
>
>         It returns fewer incidences of intersection, but I'm still not
>         sure if it does what it should.
>
>         Am Mi., 13. Sept. 2023 um 11:39 Uhr schrieb Claus Volko
>         < <mailto:>>:
>
>             I've realized the way I did it was not correct because it
>             searched for the intersection of the plane with the line,
>             which always returns a point. But I actually want to check
>             for intersection with the finite region.
>
>             However, I've not understood how region_map() would help me
>             in my endeavour. Could you perhaps point me to an example?
>
>             Thank you!
>
>             Am Mi., 13. Sept. 2023 um 11:03 Uhr schrieb Sebastien Loriot
>             < <mailto:>>:
>
>                 region_map()
>
>                 https://doc.cgal.org/latest/Shape_detection/classCGAL_1_1Shape__detection_1_1Region__growing.html#a60adecff409f4025585c64cd3bd70110 <https://doc.cgal.org/latest/Shape_detection/classCGAL_1_1Shape__detection_1_1Region__growing.html#a60adecff409f4025585c64cd3bd70110>
>
>                 Provides a property map associating indicating the
>                 region index for each
>                 input element. The key of the map depends on your
>                 instantiation. But
>                 everything should be documented in the documentation of
>                 your model.
>
>                 Best,
>
>                 Sebastien.
>
>                 On 9/13/23 09:56, Claus Volko (
>                 <mailto:> via cgal-discuss
>                 Mailing List) wrote:
>                  > OK, it seems I've found it out myself by now:
>                  >
>                  >    for (Region_growing::Primitive_and_region region :
>                 regions)
>                  >    {
>                  >        Line_3* l = new Line_3(Point_3(0, -10, 0),
>                 Point_3(0, +10, 0));
>                  >
>                  >        if (CGAL::intersection(region.first, *l))
>                  >        {
>                  >            std::cout << "Region has point on it" <<
>                 std::endl;
>                  >            break;
>                  >        }
>                  >    }
>                  >
>                  > Am Mi., 13. Sept. 2023 um 09:39 Uhr schrieb Claus Volko
>                  > < <mailto:>
>                 <mailto: <mailto:>>>:
>                  >
>                  >     OK, I've realized by now that the code I posted
>                 was not wrong but
>                  >     that what I really need is to check if a line
>                 intersects with a
>                  >     region, as I don't know the y coordinate of the
>                 point I'm looking
>                  >     for. How can that be done?
>                  >
>                  >     Am Mi., 13. Sept. 2023 um 07:04 Uhr schrieb Claus
>                 Volko
>                  >     < <mailto:>
>                 <mailto: <mailto:>>>:
>                  >
>                  >         I would like to find the region that contains
>                 a given point
>                  >         after using the region growing algorithm for
>                 shape detection. I
>                  >         tried this code but it does not seem to work:
>                  >
>                  >            for (Region_growing::Primitive_and_region
>                 region : regions)
>                  >            {
>                  >                Point_3 *point = new Point_3(0, 0, 0);
>                  >                if (region.first.has_on(*point))
>                  >                    std::cout << "Region has zero
>                 point on it" << std::endl;
>                  >            }
>                  >
>                  >         Also, if I found a region, how can I access
>                 the information
>                  >         about the boundaries of the region?
>                  >
>                  >         --
>                  >         Dipl.-Ing. Dr. Claus D. Volko, BSc
>                  > http://www.cdvolko.net/ <http://www.cdvolko.net/>
>                 <http://www.cdvolko.net/ <http://www.cdvolko.net/>>
>                  >
>                  >
>                  >
>                  >
>                  >
>                  >     --
>                  >     Dipl.-Ing. Dr. Claus D. Volko, BSc
>                  > http://www.cdvolko.net/ <http://www.cdvolko.net/>
>                 <http://www.cdvolko.net/ <http://www.cdvolko.net/>>
>                  >
>                  >
>                  >
>                  >
>                  >
>                  > --
>                  > Dipl.-Ing. Dr. Claus D. Volko, BSc
>                  > http://www.cdvolko.net/ <http://www.cdvolko.net/>
>                 <http://www.cdvolko.net/ <http://www.cdvolko.net/>>
>                  >
>                  >
>                  >
>                  >
>                  > --
>                  > You are currently subscribed to cgal-discuss.
>                  > To unsubscribe or access the archives, go to
>                  > https://sympa.inria.fr/sympa/info/cgal-discuss
>                 <https://sympa.inria.fr/sympa/info/cgal-discuss>
>                  >
>
>                 --
>                 You are currently subscribed to cgal-discuss.
>                 To unsubscribe or access the archives, go to
>                 https://sympa.inria.fr/sympa/info/cgal-discuss
>                 <https://sympa.inria.fr/sympa/info/cgal-discuss>
>
>
>
>
>             --
>             Dipl.-Ing. Dr. Claus D. Volko, BSc
>             http://www.cdvolko.net/ <http://www.cdvolko.net/>
>
>
>
>
>
>         --
>         Dipl.-Ing. Dr. Claus D. Volko, BSc
>         http://www.cdvolko.net/ <http://www.cdvolko.net/>
>
>
>
>
>
>     --
>     Dipl.-Ing. Dr. Claus D. Volko, BSc
>     http://www.cdvolko.net/ <http://www.cdvolko.net/>
>
>
>
>
>
> --
> Dipl.-Ing. Dr. Claus D. Volko, BSc
> http://www.cdvolko.net/ <http://www.cdvolko.net/>
>
>
>
>
> --
> You are currently subscribed to cgal-discuss.
> To unsubscribe or access the archives, go to
> https://sympa.inria.fr/sympa/info/cgal-discuss
>

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss

--

Dipl.-Ing. Dr. Claus D. Volko, BSc
http://www.cdvolko.net/

#define USE_POLYHEDRON

#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Polygon_mesh.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include "include/utils.h"
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <boost/optional/optional_io.hpp>

// Typedefs.
/*
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3                                     Point;
typedef CGAL::Surface_mesh<Point>                           Surface_mesh;
*/

using Kernel  = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT      = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Line_3  = typename Kernel::Line_3;

#ifdef USE_POLYHEDRON
using Polygon_mesh   = CGAL::Polyhedron_3<Kernel>;
#else
using Polygon_mesh   = CGAL::Surface_mesh<Point_3>;
#endif

using Neighbor_query = 
CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query<Polygon_mesh>;
using Region_type    = 
CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<Kernel, 
Polygon_mesh>;
using Sorting        = 
CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting<Kernel, 
Polygon_mesh, Neighbor_query>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, 
Region_type>;

typedef CGAL::Exact_predicates_inexact_constructions_kernel          K;
typedef CGAL::Surface_mesh<K::Point_3>                               Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor                 
vertex_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor               
halfedge_descriptor;

// extract vertices which are at most k (inclusive)
// far from vertex v in the graph of edges
void extract_k_ring(Polygon_mesh::Vertex_handle v,
    int k,
    std::vector<Polygon_mesh::Vertex_handle>& qv)
{
    std::map<Polygon_mesh::Vertex_handle, int>  D;
    qv.push_back(v);
    D[v] = 0;
    std::size_t current_index = 0;
    int dist_v;
    while (current_index < qv.size() && (dist_v = D[qv[current_index]]) < k)
    {
        v = qv[current_index++];
        Polygon_mesh::Halfedge_around_vertex_circulator e(v->vertex_begin()), 
e_end(e);
        do {
            Polygon_mesh::Vertex_handle new_v = e->opposite()->vertex();
            if (D.insert(std::make_pair(new_v, dist_v + 1)).second)
                qv.push_back(new_v);
        } while (++e != e_end);
    }
}

template<typename Polygon_mesh, typename Region>
void save_polygon_mesh_regions(
    Polygon_mesh& polygon_mesh,
    const Region& regions,
    const std::string fullpath,
    std::list<int> listRegions) {

    using Color = CGAL::IO::Color;
    srand(static_cast<unsigned int>(time(NULL)));

    using Face_property_color = CGAL::dynamic_face_property_t<Color>;
    using Face_color_map = typename boost::property_map<Polygon_mesh, 
Face_property_color>::type;
    Face_color_map face_color = get(Face_property_color(), polygon_mesh);

    // Iterate through all regions.
    std::ofstream out(fullpath);
    for (std::size_t i = 0; i < regions.size(); i++)
    {
        Region_growing::Primitive_and_region region = regions[i];

/*
        // Generate a random color.
        Color color(
            static_cast<unsigned char>(rand() % 256),
            static_cast<unsigned char>(rand() % 256),
            static_cast<unsigned char>(rand() % 256));
*/

        Color color(0, 0, 0, 0);

        if (std::find(listRegions.begin(), listRegions.end(), i) != 
listRegions.end())
        {
            color = Color(255, 0, 0, 255);
//            std::cout << "Region " << i << " is red." << std::endl;
        }

        // Iterate through all region items.
        for (const auto& item : region.second) {
            put(face_color, item, color);
        }
    }
    CGAL::IO::write_PLY(out, polygon_mesh, 
CGAL::parameters::face_color_map(face_color));
}

int main(int argc, char *argv[]) {

  // Load data either from a local folder or a user-provided file.
  const bool is_default_input = argc > 1 ? false : true;
  const std::string filename = is_default_input ? 
CGAL::data_file_path("meshes/building.off") : argv[1];
  std::ifstream in(filename);
  CGAL::IO::set_ascii_mode(in);

  Polygon_mesh polygon_mesh;
  if (!CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, 
polygon_mesh)) {
    std::cerr << "ERROR: cannot read the input file!" << std::endl;
    return EXIT_FAILURE;
  }

  for (auto f : faces(polygon_mesh))
  {
      auto w = CGAL::Polygon_mesh_processing::compute_face_normal(f, 
polygon_mesh);
      if (abs(w.y()) < .5)
      {
          CGAL::Euler::remove_face(f->halfedge(), polygon_mesh);
      }
  }

  for (auto f1 : faces(polygon_mesh))
  {
      for (auto f2 : faces(polygon_mesh))
      {
          if (f1 == f2) continue;

          if (f1->plane().point().y() == f2->plane().point().y())
          {
              auto x11 = f1->halfedge()->vertex()->point().x();
              auto z11 = f1->halfedge()->vertex()->point().z();
              auto x12 = f1->halfedge()->next()->vertex()->point().x();
              auto z12 = f1->halfedge()->next()->vertex()->point().z();
              auto x13 = 
f1->halfedge()->next()->next()->vertex()->point().x();
              auto z13 = 
f1->halfedge()->next()->next()->vertex()->point().z();
              auto x1a = x11 < x12 ? x11 < x13 ? x11 : x13 : x12 < x13 ? x12 
: x13;
              auto z1a = z11 < z12 ? z11 < z13 ? z11 : z13 : z12 < z13 ? z12 
: z13;
              auto x1b = x11 > x12 ? x11 > x13 ? x11 : x13 : x12 > x13 ? x12 
: x13;
              auto z1b = z11 > z12 ? z11 > z13 ? z11 : z13 : z12 > z13 ? z12 
: z13;

              auto x21 = f2->halfedge()->vertex()->point().x();
              auto z21 = f2->halfedge()->vertex()->point().z();
              auto x22 = f2->halfedge()->next()->vertex()->point().x();
              auto z22 = f2->halfedge()->next()->vertex()->point().z();
              auto x23 = 
f2->halfedge()->next()->next()->vertex()->point().x();
              auto z23 = 
f2->halfedge()->next()->next()->vertex()->point().z();
              auto x2a = x21 < x22 ? x21 < x23 ? x21 : x23 : x22 < x23 ? x22 
: x23;
              auto z2a = z21 < z22 ? z21 < z23 ? z21 : z23 : z22 < z23 ? z22 
: z23;
              auto x2b = x21 > x22 ? x21 > x23 ? x21 : x23 : x22 > x23 ? x22 
: x23;
              auto z2b = z21 > z22 ? z21 > z23 ? z21 : z23 : z22 > z23 ? z22 
: z23;

              if (x2a > x1a && x2b < x1b && z2a > z1a && z2b < z1b)
                  CGAL::Euler::remove_face(f2->halfedge(), polygon_mesh);
          }
      }
  }

  for (Polygon_mesh::Vertex_iterator v1 = polygon_mesh.vertices_begin(); v1 
!= polygon_mesh.vertices_end(); v1++)
  {
      for (Polygon_mesh::Vertex_iterator v2 = polygon_mesh.vertices_begin(); 
v2 != polygon_mesh.vertices_end(); v2++)
      {
          if (v1 == v2) continue;
          if ((abs(v1->point().x() - v2->point().x()) < .01)
              && (abs(v1->point().y() - v2->point().y()) < .01)
              && (abs(v1->point().z() - v2->point().z()) < .01))
          {
              //std::cout << "merging vertices " << v1->point() << " and " << 
v2->point() << std::endl;
              v2->point() = v1->point();
          }
      }
  }

  std::cout << "Before stitching : " << std::endl;
  std::cout << "\t Number of vertices  :\t" << 
polygon_mesh.size_of_vertices() << std::endl;
  std::cout << "\t Number of halfedges :\t" << 
polygon_mesh.size_of_halfedges() << std::endl;
  std::cout << "\t Number of facets    :\t" << polygon_mesh.size_of_facets() 
<< std::endl;

  CGAL::Polygon_mesh_processing::stitch_borders(polygon_mesh);

  std::cout << "Stitching done : " << std::endl;
  std::cout << "\t Number of vertices  :\t" << 
polygon_mesh.size_of_vertices() << std::endl;
  std::cout << "\t Number of halfedges :\t" << 
polygon_mesh.size_of_halfedges() << std::endl;
  std::cout << "\t Number of facets    :\t" << polygon_mesh.size_of_facets() 
<< std::endl;


  /*
  std::vector<Polygon_mesh::Facet_handle>  new_facets;
  std::vector<Polygon_mesh::Vertex_handle> new_vertices;
  CGAL::Polygon_mesh_processing::refine(polygon_mesh, faces(polygon_mesh),
      std::back_inserter(new_facets),
      std::back_inserter(new_vertices),
      CGAL::parameters::density_control_factor(2.));

  int vertexCount = 0;
  for (Polygon_mesh::Vertex_iterator v = polygon_mesh.vertices_begin(); v != 
polygon_mesh.vertices_end(); v++)
  {
      vertexCount++;
  }
  std::vector<Polygon_mesh::Vertex_handle> region1;
  std::vector<Polygon_mesh::Vertex_handle> region2;
  int i = 0;
  for (Polygon_mesh::Vertex_iterator v = polygon_mesh.vertices_begin(); v != 
polygon_mesh.vertices_end(); v++, i++)
  {
      if (i < vertexCount / 2)
          region1.push_back(v);
      else
          region2.push_back(v);
  }
  CGAL::Polygon_mesh_processing::fair(polygon_mesh, region1);
  CGAL::Polygon_mesh_processing::fair(polygon_mesh, region2);


  std::cout << "Before stitching : " << std::endl;
  std::cout << "\t Number of vertices  :\t" << 
polygon_mesh.size_of_vertices() << std::endl;
  std::cout << "\t Number of halfedges :\t" << 
polygon_mesh.size_of_halfedges() << std::endl;
  std::cout << "\t Number of facets    :\t" << polygon_mesh.size_of_facets() 
<< std::endl;
  
  CGAL::Polygon_mesh_processing::stitch_borders(polygon_mesh);
  CGAL::Polygon_mesh_processing::stitch_borders(polygon_mesh);
  CGAL::Polygon_mesh_processing::stitch_borders(polygon_mesh);
  CGAL::Polygon_mesh_processing::stitch_borders(polygon_mesh);

  std::cout << "Stitching done : " << std::endl;
  std::cout << "\t Number of vertices  :\t" << 
polygon_mesh.size_of_vertices() << std::endl;
  std::cout << "\t Number of halfedges :\t" << 
polygon_mesh.size_of_halfedges() << std::endl;
  std::cout << "\t Number of facets    :\t" << polygon_mesh.size_of_facets() 
<< std::endl;
  */

  const auto& face_range = faces(polygon_mesh);
  std::cout << "* number of input faces: " << face_range.size() << std::endl;
  assert(!is_default_input || face_range.size() == 32245);

  // Testing the usage of clip and copy_face_graph (Claus Volko)
#ifdef USE_POLYHEDRON
  /*
  const Kernel::Plane_3& plane = Kernel::Plane_3(CGAL::Point_3<Kernel>(0, 0, 
0),
      CGAL::Point_3<Kernel>(0, 1, 0),
      CGAL::Point_3<Kernel>(1, 1, 0));
  CGAL::Polygon_mesh_processing::clip(polygon_mesh, plane, 
CGAL::parameters::default_values());

  CGAL::Surface_mesh<Point_3> surface_mesh;
  copy_face_graph(polygon_mesh, surface_mesh);
*/
#endif

  // Default parameter values for the data file building.off.
  /*
  const FT          max_distance    = FT(1);
  const FT          max_angle       = FT(45);
  const std::size_t min_region_size = 5;
  */
  const FT          max_distance = FT(1);
  const FT          max_angle = FT(45);
  const std::size_t min_region_size = 1;

  // Create instances of the classes Neighbor_query and Region_type.
  Neighbor_query neighbor_query(polygon_mesh);

  Region_type region_type(
    polygon_mesh,
    CGAL::parameters::
    maximum_distance(max_distance).
    maximum_angle(max_angle).
    minimum_region_size(min_region_size));

  // Sort face indices.
  Sorting sorting(
    polygon_mesh, neighbor_query);
  sorting.sort();

  // Create an instance of the region growing class.
  Region_growing region_growing(
    face_range, sorting.ordered(), neighbor_query, region_type);

  // Run the algorithm.
  std::vector<typename Region_growing::Primitive_and_region> regions;
  region_growing.detect(std::back_inserter(regions));
  std::cout << "* number of found planes: " << regions.size() << std::endl;
  assert(!is_default_input || regions.size() == 365);

  
  std::list<int> listRegions = {};

  for (std::size_t i = 0; i < regions.size(); i++)
  {
      Region_growing::Primitive_and_region region = regions[i];
      Line_3* l = new Line_3(Point_3(0, -1, 0), Point_3(0, +1, 0));
      auto intersection = 
boost::get<Point_3>(&*CGAL::intersection(region.first, *l));

      if (intersection && region.first.has_on(*intersection))
      {
          auto x11 = region.second.front()->halfedge()->vertex()->point().x();
          auto z11 = region.second.front()->halfedge()->vertex()->point().z();
          auto x12 = 
region.second.front()->halfedge()->next()->vertex()->point().x();
          auto z12 = 
region.second.front()->halfedge()->next()->vertex()->point().z();
          auto x13 = 
region.second.front()->halfedge()->next()->next()->vertex()->point().x();
          auto z13 = 
region.second.front()->halfedge()->next()->next()->vertex()->point().z();
          auto x1a = x11 < x12 ? x11 < x13 ? x11 : x13 : x12 < x13 ? x12 : 
x13;
          auto z1a = z11 < z12 ? z11 < z13 ? z11 : z13 : z12 < z13 ? z12 : 
z13;
          auto x1b = x11 > x12 ? x11 > x13 ? x11 : x13 : x12 > x13 ? x12 : 
x13;
          auto z1b = z11 > z12 ? z11 > z13 ? z11 : z13 : z12 > z13 ? z12 : 
z13;

          if (intersection->x() > x1a && intersection->x() < x1b && 
intersection->z() > z1a && intersection->z() < z1b)
          {
              std::cout << "Region " << i << " has point on it : " << 
region.first << " at " << *intersection << std::endl;
              listRegions.push_back(i);
          }
      }
  }
  

  const Region_growing::Region_map& map = region_growing.region_map();

  for (std::size_t i = 0; i < regions.size(); i++)
    for (auto& item : regions[i].second) {
      if (i != get(map, item)) {
        std::cout << "Region map incorrect" << std::endl;
      }
    }

  std::vector<typename Region_growing::Item> unassigned;
  region_growing.unassigned_items(face_range, std::back_inserter(unassigned));

  for (auto& item : unassigned) {
    if (std::size_t(-1) != get(map, item)) {
      std::cout << "Region map for unassigned incorrect" << std::endl;
    }
  }

  // Save regions to a file.
  const std::string fullpath = (argc > 2 ? argv[2] : 
"planes_polygon_mesh.ply");
  save_polygon_mesh_regions(polygon_mesh, regions, fullpath, listRegions);
  //utils::save_polygon_mesh_regions(polygon_mesh, regions, fullpath);

  return EXIT_SUCCESS;
}