Skip to content
This repository has been archived by the owner on Jun 3, 2020. It is now read-only.

Commit

Permalink
Merge pull request #4 from Oslandia/polygons_with_constrains_and_holes
Browse files Browse the repository at this point in the history 
Polygons with constrains and holes
  • Loading branch information
@vmora
vmora authored Oct 30, 2019
2 parents 12dd825 + 0c29949 commit 282b115
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 36 deletions.
123 changes: 87 additions & 36 deletions fourmy/_fourmy.cc
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,92 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Polygon_2_algorithms.h>

#include <boost/python.hpp>
#include <boost/python/overloads.hpp>
#include <iostream>
#include <algorithm>
#include <string>

namespace fourmy {

struct FaceInfo2
{
FaceInfo2(){}
int nesting_level;
bool in_domain(){
return nesting_level%2 == 1;
}
};

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_2<K> Vb;
typedef CGAL::Constrained_triangulation_face_base_2<K> Fb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,K> Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<K,Fbb> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
typedef CGAL::Exact_predicates_tag Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag> CDT;
typedef CDT::Point Point;
typedef CDT::Segment Segment;
typedef CGAL::Polygon_2<K> Polygon_2;
typedef std::vector< Polygon_2 > Rings;

void
mark_domains(CDT& ct,
CDT::Face_handle start,
int index,
std::list<CDT::Edge>& border,
const std::vector< Segment >& rings)
{
if(start->info().nesting_level != -1){
return;
}
std::list<CDT::Face_handle> queue;
queue.push_back(start);
while(! queue.empty()){
CDT::Face_handle fh = queue.front();
queue.pop_front();
if(fh->info().nesting_level == -1){
fh->info().nesting_level = index;
for(int i = 0; i < 3; i++){
CDT::Edge e(fh,i);
CDT::Face_handle n = fh->neighbor(i);
if(n->info().nesting_level == -1){
if(ct.is_constrained(e) && std::find(rings.begin(), rings.end(), ct.segment(e)) != rings.end() ){
border.push_back(e);
}
else {
queue.push_back(n);
}
}
}
}
}
}

//explore set of facets connected with non constrained edges,
//and attribute to each such set a nesting level.
//We start from facets incident to the infinite vertex, with a nesting
//level of 0. Then we recursively consider the non-explored facets incident
//to constrained edges bounding the former set and increase the nesting level by 1.
//Facets in the domain are those with an odd nesting level.
void
mark_domains(CDT& cdt, const std::vector< Segment >& rings)
{
for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){
it->info().nesting_level = -1;
}
std::list<CDT::Edge> border;
mark_domains(cdt, cdt.infinite_face(), 0, border, rings);
while(! border.empty()){
CDT::Edge e = border.front();
border.pop_front();
CDT::Face_handle n = e.first->neighbor(e.second);
if(n->info().nesting_level == -1){
mark_domains(cdt, n, e.first->info().nesting_level+1, border, rings);
}
}
}

Polygon_2 polygon_from_ring(const boost::python::object & coords)
{
Expand All @@ -31,25 +98,9 @@ Polygon_2 polygon_from_ring(const boost::python::object & coords)
return std::move(poly);
}

bool is_inside(const Point & point, const Rings & rings)
{
int nesting = 0;
for(const Polygon_2 & ring : rings)
{
if (CGAL::ON_BOUNDED_SIDE == CGAL::bounded_side_2(ring.vertices_begin(), ring.vertices_end(), point, K()))
{
++nesting;
}
}

return nesting == 1;
}


boost::python::object tessellate(const boost::python::object & polygon, const boost::python::object & lines = boost::python::object(), const boost::python::object & points = boost::python::object())
{
using namespace boost::python;
Rings rings;

const char * type = extract<const char *>(polygon.attr("type"));
if (type != std::string("Polygon"))
Expand All @@ -60,26 +111,21 @@ boost::python::object tessellate(const boost::python::object & polygon, const bo
{
Polygon_2 poly(polygon_from_ring(polygon.attr("exterior").attr("coords")));
cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
rings.push_back(poly);
}
const size_t nrings = extract<size_t>(polygon.attr("interiors").attr("__len__")());
for (size_t r=0; r<nrings; r++)
{
Polygon_2 poly(polygon_from_ring(polygon.attr("interiors")[r].attr("coords")));
cdt.insert_constraint(poly.vertices_begin(), poly.vertices_end(), true);
rings.push_back(poly);
}


// insert points
if (points != boost::python::object())
{
const size_t sz = extract<size_t>(points.attr("__len__")());
for (size_t i=0; i<sz; i++)
cdt.push_back(Point(extract<double>(points[i].attr("coords")[0][0]), extract<double>(points[i].attr("coords")[0][1])));
std::vector< Segment > border;
for (CDT::Constrained_edges_iterator e=cdt.constrained_edges_begin(); e != cdt.constrained_edges_end(); e++){
const Segment s(cdt.segment(*e));
border.push_back(s);
border.push_back(s.opposite());
}


// insert line constrains
if (lines != boost::python::object())
{
Expand All @@ -96,17 +142,22 @@ boost::python::object tessellate(const boost::python::object & polygon, const bo
}
}

//mark_domains(cdt);
// insert points
if (points != boost::python::object())
{
const size_t sz = extract<size_t>(points.attr("__len__")());
for (size_t i=0; i<sz; i++)
cdt.push_back(Point(extract<double>(points[i].attr("coords")[0][0]), extract<double>(points[i].attr("coords")[0][1])));
}

mark_domains(cdt, border);

object shapely_geom = import("shapely.geometry");
list triangles;
for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin();
fit!=cdt.finite_faces_end();++fit)
{
const Point c((fit->vertex(0)->point().x()+fit->vertex(1)->point().x()+fit->vertex(2)->point().x())/3,
(fit->vertex(0)->point().y()+fit->vertex(1)->point().y()+fit->vertex(2)->point().y())/3);

if (is_inside(c, rings))
if (fit->info().in_domain())
{
list ring;
for (size_t i=0; i<4; i++)
Expand All @@ -116,7 +167,8 @@ boost::python::object tessellate(const boost::python::object & polygon, const bo
}
return shapely_geom.attr("MultiPolygon")(triangles);
}
}

} // namespace fourmy

BOOST_PYTHON_FUNCTION_OVERLOADS(tessellate_overloads, fourmy::tessellate, 1, 3)

Expand All @@ -125,4 +177,3 @@ BOOST_PYTHON_MODULE(_fourmy)
using namespace boost::python;
def("tessellate", &fourmy::tessellate, tessellate_overloads(args("polygon", "lines", "points"), "This is tesselate's docstring"));
}

11 changes: 11 additions & 0 deletions test/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,3 +15,14 @@

t = tessellate(Polygon([(0,0), (1,0), (1,1), (0,1)]), lines=None, points=[Point(.9, .9)])
print(t.wkt)

t = tessellate(Polygon(((-4.165589, -29.100525), (8.623957000000001, -28.461553), (21.413503, -27.822581), (10.706928, -13.90117), (0.000353, 0.020242), (-2.082618, -14.540141), (-4.165589, -29.100525))))
print(t.wkt)
assert(len(t) == 4)
print(len(t))
for e in t:
print(e.area, e.wkt)

t = tessellate(Polygon([(0,0), (1,0), (1,1), (0,1)], [[(.2,.2), (.2,.8), (.8,.8), (.8, .2)]]), lines=[LineString([(.1, .1), (.9,.1)]), LineString([(.9, .1), (.9,.9)]), LineString([(.9, .9), (.1,.9)]), LineString([(.1, .9), (.1,.1)])])
print(len(t))
print(t.wkt)

0 comments on commit 282b115

Please sign in to comment.