Not removing polygons, as it has unintended consequences, but polygon…

…s are ok.

geomodelr/cpp/geomodel.cpp

@@ -148,3 +148,37 @@ BOOST_PYTHON_MODULE(cpp)
 						    .add_property("matches", &Model::get_matches, &Model::set_matches);
 }
 
+
+wstring human_failure_type( const geometry::validity_failure_type& fail )
+{
+	switch ( fail ) {
+		case geometry::validity_failure_type::no_failure:
+			return L"no failure";
+		case geometry::validity_failure_type::failure_few_points:
+			return L"failure few points";
+		case geometry::validity_failure_type::failure_wrong_topological_dimension:
+			return L"failure wrong topological dimension";
+		case geometry::validity_failure_type::failure_spikes:
+			return L"failure spikes";
+		case geometry::validity_failure_type::failure_duplicate_points:
+			return L"failure duplicate points";
+		case geometry::validity_failure_type::failure_not_closed:
+			return L"failure not closed";
+		case geometry::validity_failure_type::failure_self_intersections:
+			return L"failure self intersections";
+		case geometry::validity_failure_type::failure_wrong_orientation:
+			return L"failure wrong orientation";
+		case geometry::validity_failure_type::failure_interior_rings_outside:
+			return L"failure interior rings outside";
+		case geometry::validity_failure_type::failure_nested_interior_rings:
+			return L"failure nested interior rings";
+		case geometry::validity_failure_type::failure_disconnected_interior:
+			return L"failure disconnected interior";
+		case geometry::validity_failure_type::failure_intersecting_interiors:
+			return L"failure intersecting interiors";
+		case geometry::validity_failure_type::failure_wrong_corner_order:
+			return L"failure wrong corner order";
+		default:
+			return L"unknown";
+	}
+}

geomodelr/cpp/geomodel.hpp

@@ -35,6 +35,8 @@ using std::string;
 using std::wstring;
 using std::map;
 
+wstring human_failure_type( const geometry::validity_failure_type& fail );
+
 static const double tolerance = 1e-15;
 
 typedef geometry::model::point<double, 2, geometry::cs::cartesian> point2;

geomodelr/cpp/match.cpp

@@ -63,12 +63,15 @@ pylist Match::get() const {
 void Match::match_polygons() {
 	map<wstring, vector<int>> units_a;
 	map<wstring, vector<int>> units_b;
+
 	for ( size_t i = 0; i < this->a->polygons.size(); i++ ) {
 		units_a[this->a->units[i]].push_back(i);
 	}
+
 	for ( size_t i = 0; i < this->b->polygons.size(); i++ ) {
 		units_b[this->b->units[i]].push_back(i);
 	}
+
 	vector<std::pair<int, int>> m;
 	for ( auto it = units_a.begin(); it != units_a.end(); it++ ) {
 		if ( units_b.find(it->first) != units_b.end() ) {
@@ -84,6 +87,7 @@ void Match::match_polygons() {
 			}
 		}
 	}
+
 	this->set( m );
 }
 

geomodelr/cpp/model.cpp

@@ -20,6 +20,7 @@
 #include <algorithm>
 #include <functional>
 #include <cmath>
+#include <iostream>
 
 bool Model::verbose = false;
 
@@ -41,7 +42,7 @@ topography(nullptr)
 	}
 	size_t nsects = python::len(sections);
 	for ( size_t i = 0; i < nsects; i++ ) {
-		const wstring& name    = python::extract<wstring>(sections[i][0]);
+		wstring name    = python::extract<wstring>(sections[i][0]);
 		double cut             = python::extract<double>(sections[i][1]);
 		const pylist& points   = python::extract<pylist>(sections[i][2]);
 		const pylist& polygons = python::extract<pylist>(sections[i][3]);
@@ -50,13 +51,15 @@ topography(nullptr)
 		const pylist& lnames   = python::extract<pylist>(sections[i][6]);
 		this->sections.push_back(new Section(name, cut, points, polygons, units, lines, lnames));
 	}
+
 	std::sort(this->sections.begin(), this->sections.end(), [](const Section* a, const Section* b){ return a->cut < b->cut; });
 	for ( size_t i = 0; i < this->sections.size(); i++ ) {
 		this->cuts.push_back(this->sections[i]->cut);
 	}
 }
 
 Model::~Model(){
+	/* clears matches and sections from memory */
 	if ( this->topography != nullptr ) {
 		delete this->topography;
 	}
@@ -66,6 +69,7 @@ Model::~Model(){
 	}
 }
 void Model::clear_matches() {
+	/* clears all the matches from memory */
 	for ( Match * c: this->match ) {
 		delete c;
 	}
@@ -74,6 +78,7 @@ void Model::clear_matches() {
 
 pydict Model::make_matches() {
 	this->clear_matches();
+
 	map<wstring, vector<triangle_pt>> faults;
 	auto add_to_faults = [&](const map<wstring, vector<triangle_pt>>& m) {
 		for ( auto it = m.begin(); it != m.end(); it++ ) {
@@ -82,19 +87,23 @@ pydict Model::make_matches() {
 			f.insert(f.end(), it->second.begin(), it->second.end());
 		}
 	};
-	// Get all the matching faults.
+
 	for ( size_t i = 1; i < this->sections.size(); i++ ) {
+
 		this->match.push_back(new Match(this->sections[i-1], this->sections[i]));
+		// Match the polygons.
 		this->match.back()->match_polygons();
+		// Get the matching faults.
 		map<wstring, vector<triangle_pt>> m = this->match.back()->match_lines();
 		add_to_faults(m);
 	}
-
+	
 	// Get the extended faults from the begining.
 	if ( this->sections.size() and (this->sections[0]->cut - this->cuts_range.first) > tolerance ) {
 		map<wstring, vector<triangle_pt>> m = this->sections[0]->last_lines(true, this->cuts_range.first);
 		add_to_faults(m);
 	}
+
 	// Get the extended faults from the end.
 	if ( this->sections.size() and (this->cuts_range.second - this->sections.back()->cut) > tolerance ) {
 		// Get the extended faults from the end.
@@ -107,12 +116,12 @@ pydict Model::make_matches() {
 		point3 glp = this->to_inverse_point(p, gz(pt));
 		return python::make_tuple(gx(glp), gy(glp), gz(glp));
 	};
-
+	
 	auto triangle_coords = [point_coords]( const triangle_pt& tr ) {
 		return python::make_tuple(point_coords(g0(tr)), point_coords(g1(tr)), point_coords(g2(tr)));
 	};
-
-	// Now convert to python and return.
+	
+	// Now convert faults to python and return.
 	pydict ret;
 	for ( auto it = faults.begin(); it != faults.end(); it++ ) {
 		pylist tris;
@@ -267,6 +276,7 @@ vector<Model::Possible> Model::get_candidates(size_t a_idx, const point2& pt ) c
 			possible.insert(a_possible(c));
 		}
 	}
+
 	return vector<Model::Possible>(possible.begin(), possible.end());
 }
 
@@ -359,6 +369,7 @@ vector<Model::Possible> Model::all_closest( size_t a_idx, const point2& pt ) con
 }
 
 std::tuple<int, int, double> Model::closest_to( size_t a_idx, const point2& pt, double cut ) const {
+	/* Returns the closest formation to a model. */
 	double s_dist = this->cuts[a_idx+1]-this->cuts[a_idx];
 	double cut_rem = cut - this->cuts[a_idx];
 	double mult_b = cut_rem/s_dist;
@@ -457,6 +468,7 @@ pylist Model::possible_closest(const pyobject& pypt) const {
 	auto it = std::upper_bound(this->cuts.begin(), this->cuts.end(), mp.second);
 	size_t a_idx = it - this->cuts.begin();
 	// If it's behind the last or above the first, return the closest in the section.
+
 	if ( a_idx <= 0 or a_idx >= this->sections.size() ) {
 		const Section& s = ( a_idx <=0 ) ? *(this->sections.front()) : *(this->sections.back());
 		std::pair<int, double> cls = s.closest_to(mp.first, geometry::index::satisfies(always_true));
@@ -465,15 +477,20 @@ pylist Model::possible_closest(const pyobject& pypt) const {
 		ret.append(python::make_tuple(unit, cls.second, cls.second));
 		return ret;
 	}
+
 	a_idx--;
+
 	vector<Model::Possible> ap = this->all_closest(a_idx, mp.first);
+
 	pylist ret;
+
 	for ( size_t i = 0; i < ap.size(); i++ ) {
 		int a_match = ap[i].a_match;
 		int b_match = ap[i].b_match;
 		wstring unit = ( a_match != -1 ) ? this->sections[a_idx]->units[a_match] : this->sections[a_idx+1]->units[b_match];
 		ret.append(python::make_tuple(unit, ap[i].a_dist, ap[i].b_dist));
 	}
+
 	return ret;
 }
 
@@ -498,6 +515,7 @@ pytuple Model::closest_topo( const pyobject& pypt ) const {
 }
 
 pytuple Model::closest(const pyobject& pypt) const {
+	/* python exposed functions to search for the closest formation */
 	const double& p0 = python::extract<double>(pypt[0]);
 	const double& p1 = python::extract<double>(pypt[1]);
 	const double& p2 = python::extract<double>(pypt[2]);
@@ -512,9 +530,10 @@ pytuple Model::closest(const pyobject& pypt) const {
 
 	std::pair<point2, double> mp = this->to_model_point(pt);
 	auto it = std::upper_bound(this->cuts.begin(), this->cuts.end(), mp.second);
+
 	size_t a_idx = it - this->cuts.begin();
 	// If it's behind the last or above the first, return the closest in the section.
-
+	
 	auto closest_single =[&](const Section& s) {
 		std::pair<int, double> cls = s.closest_to(mp.first, geometry::index::satisfies(always_true));
 		if ( cls.first == -1 ) {
@@ -523,11 +542,13 @@ pytuple Model::closest(const pyobject& pypt) const {
 		wstring unit = s.units[cls.first];
 		return python::make_tuple(unit, cls.second);
 	};
+
 	// For a cut below the lowest or above the highest.
 	if ( a_idx <= 0 or a_idx >= this->sections.size() ) {
 		const Section& s = ( a_idx <=0 ) ? *(this->sections.front()) : *(this->sections.back());
 		return closest_single(s);
 	}
+
 	// Check if it crosses a fault and then only evaluate the cross in the actual site.
 	a_idx--;
 	int crosses = this->match[a_idx]->crosses_triangles(mp.first, mp.second);

geomodelr/cpp/section.cpp

@@ -58,17 +58,22 @@ Section::Section(const wstring& name, double cut, const pylist& points,
 				}
 			}
 		}
-
-		string message;
-		if ( not geometry::is_valid(pol, message) or not geometry::is_simple(pol) ) {
-			if ( Model::verbose ) {
-				string su(unit.begin(), unit.end());
-				std::cerr << "removed polygon with unit " << su << " valid: " << (geometry::is_valid(pol)?"true":"false")
-					  <<" simple: " << (geometry::is_simple(pol)?"true":"false") << "\n";
+
+		geometry::validity_failure_type failure;
+		if ( not geometry::is_valid(pol, failure) ) {
+			geometry::correct(pol);
+			string reason;
+			if ( not geometry::is_valid(pol, reason) ) {
+				if ( Model::verbose ) {
+					std::wcerr << L"non valid polygon in section " << name << L" from with unit " << unit << L" valid: " << (geometry::is_valid(pol)?L"true":L"false")
+						   << L" simple: " << (geometry::is_simple(pol)?L"true":L"false") << "\n";
+				}
+				// continue; not avoiding non valid polygons, as they have been validated by shapely. Somehow these polygons get wronget.
 			}
+		}
+		if ( not geometry::is_simple(pol) ) {
 			continue;
 		}
-
 		// Calculate the envelope and add it to build the rtree layer.
 		box env;
 		geometry::envelope(pol, env);

geomodelr/model.py

@@ -131,14 +131,15 @@ def __init__( self, geolojson ):
 
         super(GeologicalModel, self).__init__(cuts, list(base_point), list(direction), geomap, topography, sections)
         self.make_matches()
-
+        
     def make_matches(self):
         """ Prepares the model to query by matching polygons and lines.
             It finds which polygons, when projected to the next cross section,
             intersect. After that, it tries to match faults with the same name
             by triangulating them and trying to find a continuous set of triangles
             between the two lines that go from the ends to the other side.
         """
+
         self.joined_faults = super(GeologicalModel, self).make_matches()
 
     def print_information( self, verbose=False ):

geomodelr/test.py

@@ -33,17 +33,19 @@
 
 class TestGeoModelR(unittest.TestCase):
     def setUp(self):
+        pass
         # Profile GeoModelR
-        self.pr = cProfile.Profile()
-        self.pr.enable()
+        # self.pr = cProfile.Profile()
+        # self.pr.enable()
     def tearDown(self):
+        pass
         # Print profiling of GeoModelR.
-        self.pr.disable()
-        s = StringIO.StringIO()
-        sortby = 'cumulative'
-        ps = pstats.Stats(self.pr, stream=s).sort_stats(sortby)
-        ps.print_stats()
-        print >> sys.stderr, s.getvalue()
+        # self.pr.disable()
+        # s = StringIO.StringIO()
+        # sortby = 'cumulative'
+        # ps = pstats.Stats(self.pr, stream=s).sort_stats(sortby)
+        # ps.print_stats()
+        # print >> sys.stderr, s.getvalue()
 
     # Tests function fault plane for lines.
     def test_faultplane_for_lines(self):
@@ -215,6 +217,27 @@ def test_model_matching(self):
             model.matches = [((u'B-B', u'C-C'), [(0, 0), (2, 2), (3, 4)])]
         with self.assertRaises(Exception):
             model.matches = [((u'C-C', u'D-D'), [(0, 0), (2, 2), (3, 4)])]
+
+    def test_closest_single(self):
+        # Evaluate a model that has a hole in the middle.
+        points_1 = [[0, 0], [3, 0], [3, 3], [0, 3], [1, 1], [2, 1], [2, 2], [1, 2]]
+        points_2 = [[0, 0], [3, 0], [3, 3], [0, 3]]
+
+        polygons_1 = [[[0, 1, 2, 3], [7, 6, 5, 4]], [[4, 5, 6, 7]]]
+        polygons_2 = [[[0, 1, 2, 3]]] 
+
+        units_1 = ["unit1", "unit2"]
+        units_2 = ["unit1"]
+
+        model = cpp.Model([0,0,3,3],[0, 0], [1, 0], {}, {}, [("A-A", 1, points_1, polygons_1, units_1, [], []), ("B-B", 2, points_2, polygons_2, units_2, [], [])])
+        model.make_matches()
+        self.assertEqual(model.model_point([1.5, 0.1, 1.5]), (1.5, 1.5, 0.1))
+        pos_cls = model.possible_closest((1.5, 1.1, 1.5))
+        self.assertEqual(pos_cls[0][0], 'unit1')
+        cls = model.closest((1.5, 1.1, 1.5))
+        self.assertEqual(cls[0], 'unit2')
+
+
     # Test the possible closest, all the units that can be a match given a line.
     def test_possible_closest(self):
         # Evaluate simple models.
@@ -452,7 +475,7 @@ def query_func(p):
         self.assertGreaterEqual(ref_grid['grid'].shape[2], ref_grid['units'].shape[2])
 
         grid = utils.generate_octtree_grid(query_func, bbox, 3, 3, 3)
-        vols, elems = utils.octtree_volume_calculation(query_func, bbox, 10, 3)
+        vols, elems = utils.octtree_volume_calculation(query_func, bbox, 10, 2)
 
 def main(args=None):
     unittest.main()