Fixed problems with code of closest.

geomodelr/cpp/basic.hpp

@@ -123,7 +123,7 @@ class add_funcs {
 	}
 
 	bool operator() ( const value_f& v ) const {
-		return this->fa(v) and this->fb(v);
+		return (this->fa(v) and this->fb(v));
 	}
 
 };

geomodelr/cpp/model.cpp

@@ -138,13 +138,18 @@ void Model::set_matches( const vector< std::tuple< std::tuple<wstring, wstring>,
 }
 
 std::pair<int, double> Model::closest_match( bool a, int a_idx, int pol_idx, const point2& pt ) const {
+	// Reference to match.
 	const Match& m = *(this->match[a_idx]);
+	// Reference to section s.
 	const Section& s = (a) ? *(this->sections[a_idx+1]) : *(this->sections[a_idx]);
+	// Get match depending on the direction.
 	const map<int, vector<int>>& mp = (a) ? m.a_to_b : m.b_to_a;
-	auto it = (a) ? mp.find(pol_idx): mp.find(pol_idx);
+	auto it = mp.find(pol_idx);
+	// This polygon does not have matches.
 	if ( it == mp.end() ) {
 		return std::make_pair(-1, std::numeric_limits<double>::infinity());
 	}
+	// In all the matches of this polygon, check the minimum distance.
 	const vector<int>& op = it->second;
 	double mindist = std::numeric_limits<double>::infinity();
 	int minidx = -1;
@@ -156,6 +161,7 @@ std::pair<int, double> Model::closest_match( bool a, int a_idx, int pol_idx, con
 			minidx = pl;
 		}
 	}
+	// Return minimum distance and index of match.
 	return std::make_pair( minidx, mindist );
 }
 

geomodelr/cpp/model.hpp

@@ -57,11 +57,17 @@ class Model {
 
 	struct Possible {
 		Possible(int a_match, int b_match, double a_dist, double b_dist);
+		// The match in the section a.
 		int a_match;
+		// The match in the section b.
 		int b_match;
+		// The polygon distance in section a.
 		double a_dist;
+		// The polygon distance in section b.
 		double b_dist;
+		// Order to be able to put in a map.
 		bool operator<( const Possible& other ) const;
+		// The given distance, at the given c.
 		double distance( double c ) const;
 	};
 
@@ -76,7 +82,7 @@ class Model {
 
 		// Given a match in section a, find a match in section b
 		double inf = std::numeric_limits<double>::infinity();
-		auto a_possible = [&](const std::pair<int, double>& cls) { 
+		auto a_possible = [&](const std::pair<int, double>& cls) {
 			// If it didn't find a match, just return no match.
 			if ( cls.first == -1 ) {
 				return Possible(-1, -1, inf, inf);
@@ -90,7 +96,7 @@ class Model {
 			}
 		};
 
-		// Given a match in section b, find a match in section b.
+		// Given a match in section b, find a match in section a.
 		auto b_possible = [&](const std::pair<int, double>& cls) { 
 			if ( cls.first == -1 ) {
 				return Possible(-1, -1, inf, inf);
@@ -107,13 +113,14 @@ class Model {
 		std::set<Possible> possible;
 		// Find the closest polygon to a.
 		auto closest_in_a = s_a.closest(pt_a, predicates);
+
 		// Find the closest polygon in s_a => clst_a.
 		Possible clst_a = a_possible(closest_in_a);
 
 		// Find all polygons closer than the one we just found.
 		if ( clst_a.a_match != -1 ) {
 			possible.insert(clst_a);
-			auto not_first = [&](const value_f& v) { return g2(v) != clst_a.b_match; };
+			auto not_first = [&](const value_f& v) { return g2(v) != clst_a.b_match; }; 
 			// Find the pols in s_b that are closer mindist => all_cls_a.
 			vector<std::pair<int, double>> clsr_b = s_b.closer_than( pt_b, clst_a.b_dist, add_funcs<decltype(not_first), decltype(predicates)>(not_first, predicates) );
 			if ( not clsr_b.size() ) 
@@ -131,9 +138,9 @@ class Model {
 
 		if ( clst_b.b_match != -1 ) {
 			possible.insert(clst_b);
-			auto not_first = [&](const value_f& v) { return g2(v) != clst_b.b_match; };
+			auto not_first = [&](const value_f& v) { return g2(v) != clst_b.a_match; }; // CHECK it's b_match.
 			// Find the pols in s_b that are closer mindist => all_cls_a.
-			vector<std::pair<int, double>> clsr_a = s_a.closer_than( pt_b, clst_b.a_dist, add_funcs<decltype(not_first), decltype(predicates)>(not_first, predicates) );
+			vector<std::pair<int, double>> clsr_a = s_a.closer_than( pt_a, clst_b.a_dist, add_funcs<decltype(not_first), decltype(predicates)>(not_first, predicates) ); // CHECK it's pt_b
 			for ( const auto& c: clsr_a ){
 				possible.insert(a_possible(c));
 			}
@@ -143,8 +150,6 @@ class Model {
 	}
 
 
-
-
 	vector<Possible> all_closest( size_t a_idx, const point2& pt_a, const point2& pt_b ) const;
 
 	template<typename Predicates>
@@ -153,6 +158,7 @@ class 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;
+		// assert 0.0 <= mult_b <= 1.0
 		double mult_a = 1.0-mult_b;
 
 		vector<Possible> possible = this->get_candidates(a_idx, pt_a, pt_b, predicates);

geomodelr/test.py

@@ -437,7 +437,7 @@ def test_aburra_valley(self):
         m = model.GeologicalModel(json.loads(f.read()))
         m.height([813487.938015, 1164500.0])
         m.height([80000, 1100000])
-        bbox = m.geojson['bbox']
+        bbox = m.bbox
         t = {}
         def query_func(p):
             q = m.closest(p)[0]

geomodelr/utils.py

@@ -16,6 +16,7 @@
 
 # The utils file contains scripts that can be used by users to 
 # make calculations of their models.
+
 from model import GeologicalModel
 from shared import ModelException, TaskException
 import random
@@ -32,16 +33,13 @@
     import matplotlib.pyplot as plt
     from mpl_toolkits.mplot3d import Axes3D
 except ImportError:
+    print "Plotting solids is not supported"
     pass
 
 import itertools
 import gc
 
-
-def calculate_isosurface(model, unit, grid_divisions):
-
-    bbox = list(model.bbox) # Copy so original is not modified.
-
+def calculate_isovalues(model, unit, grid_divisions, bbox):
     dx = (bbox[3]-bbox[0])/grid_divisions
     dy = (bbox[4]-bbox[1])/grid_divisions
     dz = (bbox[5]-bbox[2])/grid_divisions
@@ -58,6 +56,47 @@ def calculate_isosurface(model, unit, grid_divisions):
     signed_distance = lambda x,y,z: model.signed_distance_bounded(unit, (x,y,z))
     vsigned_distance = np.vectorize(signed_distance, otypes=[np.float])
     sd = vsigned_distance( X, Y, Z )
+    return (X, Y, Z, sd)
+
+def check_bbox_surface( sd ):
+    """
+    Checks the bbox of the object.
+    """
+    mn = [ float("inf"),  float("inf"),  float("inf")]
+    mx = [-float("inf"), -float("inf"), -float("inf")]
+
+    for i in xrange(sd.shape[0]):
+        for j in xrange(sd.shape[1]):
+            for k in xrange(sd.shape[2]):
+                if sd[i,j,k] < 0:
+                    mn[0] = min( i, mn[0] )
+                    mn[1] = min( j, mn[1] )
+                    mn[2] = min( k, mn[2] )
+
+                    mx[0] = max( i, mx[0] )
+                    mx[1] = max( j, mx[1] )
+                    mx[2] = max( k, mx[2] )
+
+    if mn[0] > mx[0]:
+        raise TaskException("This model does not contain the unit or the sample is too coarse")
+    # The bbox is the first possitive.
+    mn = [ max(mn[i]-1, 0) for i in xrange(3) ]
+    mx = [ min(mx[i]+1, sd.shape[i]-1) for i in xrange(3) ]
+    return mn + mx
+
+def calculate_isosurface(model, unit, grid_divisions):
+
+    bbox = list(model.bbox) # Copy so original is not modified.
+    X, Y, Z, sd = calculate_isovalues( model, unit, grid_divisions, bbox )
+
+    # Check if the surface covers a very small volume, then reduce the bbox to calculate surface.
+    bb = check_bbox_surface( sd )
+    total_cells = grid_divisions ** 3
+    obj_cells = (bb[3]-bb[0])*(bb[4]-bb[1])*(bb[5]-bb[2])
+    # If the object is at least 8 times smaller than the full bbox, it will benefit lots from a thinner sample.
+    if ( float(total_cells) / float(obj_cells) ) > 8.0:
+        bbox = [X[bb[0], 0, 0], Y[0, bb[1], 0], Z[0, 0, bb[2]], X[bb[3], 0, 0], Y[0, bb[4], 0], Z[0, 0, bb[5]]]
+        X, Y, Z, sd = calculate_isovalues( model, unit, grid_divisions, bbox )
     try:
         vertices, simplices, normals, values = measure.marching_cubes(sd, 0)
     except ValueError: