[geogram] Update to 1.6.1

CMakeLists.txt

@@ -39,7 +39,7 @@ include(cmake/geogram.cmake)
 
 set(VORPALINE_VERSION_MAJOR 1)
 set(VORPALINE_VERSION_MINOR 6)
-set(VORPALINE_VERSION_PATCH 0)
+set(VORPALINE_VERSION_PATCH 1)
 set(VORPALINE_VERSION ${VORPALINE_VERSION_MAJOR}.${VORPALINE_VERSION_MINOR}.${VORPALINE_VERSION_PATCH})
 
 set(VORPALINE_INCLUDE_SUBPATH geogram${VORPALINE_VERSION_MAJOR})
@@ -96,8 +96,6 @@ configure_file(
 add_custom_target(uninstall
     COMMAND ${CMAKE_COMMAND} -P ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake)
 
-
-
 ##############################################################################
 # Geogram/Vorpaline sources
 

CMakeOptions.txt

@@ -1 +1,3 @@
 set(CPACK_GENERATOR RPM)
+set(CMAKE_INSTALL_PREFIX /Users/sgaspari/dev/code/popart/libs/sandbox/geogram/build/Darwin-clang-dynamic-Release/install)
+

cmake/platforms/Darwin-clang.cmake

@@ -34,7 +34,7 @@ add_flags(CMAKE_CXX_FLAGS -msse3)
 add_flags(CMAKE_C_FLAGS -msse3)
 
 # C++11 standard
-add_flags(CMAKE_CXX_FLAGS -Qunused-arguments -std=c++11 -Wno-c++98-compat)
+add_flags(CMAKE_CXX_FLAGS -Qunused-arguments -std=c++11 -stdlib=libc++ -Wno-c++98-compat)
 
 # Enable glibc parallel mode
 #add_flags(CMAKE_CXX_FLAGS -D_GLIBCXX_PARALLEL)

cmake/platforms/Darwin-clang/config.cmake

@@ -1,5 +1,5 @@
 set(VORPALINE_ARCH_64 true)
-include(${CMAKE_SOURCE_DIR}/cmake/platforms/Darwin-clang.cmake)
+include(${GEOGRAM_SOURCE_DIR}/cmake/platforms/Darwin-clang.cmake)
 add_flags(CMAKE_CXX_FLAGS -m64)
 add_flags(CMAKE_C_FLAGS -m64)
 

src/lib/exploragram/hexdom/FF.cpp

@@ -55,14 +55,14 @@ namespace GEO {
         m = p_m;
         compute_tet_edge_graph(m,v2e, true); // here  need a bidirectionl edge graph to speed up the LBFGS part
         Attribute<vec3> lockB(m->vertices.attributes(), "lockB");
-        num_l_v = 0;
-        num_ln_v = 0;
-        FOR(v, m->vertices.nb()) {
-            if (lockB[v][2] == 1 && lockB[v][0] == 0 && v > 0 && lockB[v-1][0] != 0) num_l_v = v - 1;
-            if (lockB[v][2] == 0 && v > 0 && num_ln_v == 0) num_ln_v = v - 1;
+        num_l_v = m->vertices.nb();
+        num_ln_v = m->vertices.nb();
+        FOR(inv_v, m->vertices.nb()) {
+	    index_t v = m->vertices.nb()-1 - inv_v;
+	    if (lockB[v][0] <.5) num_l_v = v;
+	    if (lockB[v][2] <.5) num_ln_v = v;
         }
-		if (num_ln_v == 0) num_ln_v = m->vertices.nb();
-
+	if (num_ln_v == 0) num_ln_v = m->vertices.nb();
     }
 
     FFopt::~FFopt() {
@@ -151,17 +151,17 @@ namespace GEO {
                 if (generate_sh) sh[v] = fv;
                 if (v >= num_l_v) {
                     vec3 oldz = col(B[v], 2);
-                    if (v > start) {
+					if (v > start && v > num_l_v) {
                         vec3  prev = mat3_to_euler(normalize_columns(B[v - 1]));
-                        B[v] = fv.project_mat3(1e-3, 1e-5, &prev);
+						B[v] = fv.project_mat3(1e-3, 1e-5, &prev);
                     } else 
                         B[v] = fv.project_mat3(1e-3, 1e-5, NULL);
-                    if (v < num_ln_v) {
+                    if (v <= num_ln_v) {
                         AxisPermutation ap;
                         ap.make_col2_equal_to_z(B[v], normalize(oldz));
                         B[v] = Frame(B[v]).apply_permutation(ap);
-                        FOR(d, 3) B[v](d, 2) = oldz[d];// restore size as well
-                    }
+						FOR(d, 3) B[v](d, 2) = oldz[d];// restore size as well
+					}
                 }
             }
         }
@@ -218,6 +218,12 @@ namespace {
         index_t nverts = FF_LBFGS::ffopt_ptr->m->vertices.nb();
         geo_assert(N == 3 * (nverts - FF_LBFGS::Num_ln_v) + FF_LBFGS::Num_ln_v);
 
+
+		Attribute<bool> border_vertex(FF_LBFGS::ffopt_ptr->m->vertices.attributes(), "border_vertex");
+		FOR(v, FF_LBFGS::ffopt_ptr->m->vertices.nb()) border_vertex[v] = false;
+		FOR(c, FF_LBFGS::ffopt_ptr->m->cells.nb())  FOR(cf, 4) if (FF_LBFGS::ffopt_ptr->m->cells.adjacent(c, cf) == NOT_AN_ID)
+			FOR(cfv, 3) border_vertex[FF_LBFGS::ffopt_ptr->m->cells.facet_vertex(c, cf, cfv)] = true;
+
 #ifdef GEO_OPENMP
         int max_threads = omp_get_max_threads();
 #else
@@ -270,21 +276,25 @@ namespace {
                     mSinv = mSinv.transpose();
                     mPst = mSinv* mR; // Pst = S^{-1} * R
 
-
-                    if (v1 > v2) FOR(i, 3) 
-                        f_chunks[thread_id] += 10. / 3.*(pow(mPst(0,i) * mPst(1,i), 2) + pow(mPst(0,i) * mPst(2,i), 2) + pow(mPst(1,i )* mPst(2,i), 2));
+					double scale = 1.;
+					if (HexdomParam::FF.rigid_border) {
+						if (border_vertex[v1])scale += 100.;
+						if (border_vertex[v2])scale += 100.;
+					}
+					if (v1 > v2) FOR(i, 3)
+                        f_chunks[thread_id] += scale *(10. / 3.*(pow(mPst(0,i) * mPst(1,i), 2) + pow(mPst(0,i) * mPst(2,i), 2) + pow(mPst(1,i )* mPst(2,i), 2)));
 
                     if (v1 >= FF_LBFGS::Num_ln_v) {
                         index_t idx = FF_LBFGS::Num_ln_v + (v1 - FF_LBFGS::Num_ln_v) * 3;
                         FOR(d,3) {                         
                             mJPst[d] = mSinv* mJR[d]; // JPst[d] = S^{-1} * JR[d]
                             FOR(i,3)FOR(j,3)
-                                g[idx + d] += 20. / 3.*mPst(i,j) * (pow(mPst(i, (j + 1)%3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[d](i , j);     
+                                g[idx + d] += scale *(20. / 3.*mPst(i,j) * (pow(mPst(i, (j + 1)%3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[d](i , j));
                         }
                     } else if (v1 >= FF_LBFGS::Num_l_v) {
                         mJPst[0] = mSinv* mJR[0]; // JPst[0] = S^{-1} * JR[0] ; JPst[1] and JPst[2] are not initialized
                         FOR(i, 3)FOR(j, 3)
-                            g[v1] += 20. / 3.*mPst(i ,j) * (pow(mPst(i ,(j + 1) % 3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[0](i , j);
+                            g[v1] += scale *(20. / 3.*mPst(i ,j) * (pow(mPst(i ,(j + 1) % 3), 2) + pow(mPst(i , (j + 2) % 3), 2))*mJPst[0](i , j));
                     }
                 } // v2
             } // v1
@@ -395,6 +405,7 @@ namespace GEO {
                         B[v] = B[v]* M.get_mat();
                         FOR(d, 3) if (std::abs(lockU[v][d]) < .1) lockU[v][d] = 0;
                     }
+					Q.push_back(v);
                 }
             }
         }

src/lib/exploragram/hexdom/PGP.cpp

@@ -45,6 +45,9 @@
 #include <exploragram/hexdom/time_log.h>
 
 
+
+
+#include <exploragram/hexdom/quadmesher.h>
 #include <geogram/NL/nl.h>
 
 #include <algorithm>
@@ -77,7 +80,11 @@ namespace GEO {
 			AxisPermutation r = Rij(m, B, v[i], v[(i + 1) % 3]);
 			bool inv;
 			index_t e = edge_from_vertices(v[i], v[(i + 1) % 3], inv);
-			t += R*(inv ? -(r.inverse()*tij[e]) : tij[e]);
+			int i1 = int(tij[e][0]);
+			int i2 = int(tij[e][1]);
+			int i3 = int(tij[e][2]);
+			vec3i ltij(i1,i2,i3);
+			t += R*(inv ? -(r.inverse()*ltij) : ltij);
 			R = R* r.inverse().get_mat();
 		}
 
@@ -245,8 +252,7 @@ namespace GEO {
 			}
 
 
-			if (org != NOT_AN_ID)
-				FOR(i, 4) {
+			if (org != NOT_AN_ID) FOR(i, 4) {
 				index_t corner = m->cells.corner(c, i);
 				UC[corner] = U[m->cells.vertex(c, i)];
 				if (i != org) {
@@ -255,8 +261,9 @@ namespace GEO {
 					bool inv;
 					index_t e = edge_from_vertices(m->cells.vertex(c, org), m->cells.vertex(c, i), inv);
 					geo_assert(e != NOT_AN_ID);
-					vec3i t2 = !inv ? tij[e] : -(change.inverse()*tij[e]);
+					vec3 t2 = !inv ? tij[e] : -(change.inverse()*tij[e]);
 					UC[corner] += vec3(t2[0], t2[1], t2[2]);
+
 				}
 			}
 			else geo_assert(!has_param[m->cells.facet(c, 0)] && !has_param[m->cells.facet(c, 1)] && !has_param[m->cells.facet(c, 2)] && !has_param[m->cells.facet(c, 3)]);
@@ -322,7 +329,6 @@ namespace GEO {
 
 		nlDeleteContext(nlGetCurrent());
 
-//		snap_U_to_round();
 
 		FOR(e, m->edges.nb()) {
 			index_t i = m->edges.vertex(e, 0);

src/lib/exploragram/hexdom/PGP.h

@@ -99,6 +99,7 @@ namespace GEO {
 
 /*
         void snap_U_to_round(double eps = 0.05) {
+			return;
             FOR(v, m->vertices.nb()) FOR(d, 3) {
                 if (std::abs(U[v][d] - round(U[v][d])) < eps) {
                     U[v][d] = round(U[v][d]);
@@ -116,7 +117,8 @@ namespace GEO {
         // A PGP solution is not only a mesh + attrib: it also requires these datas
 
         Attribute<vec3> corr;
-        Attribute<vec3i> tij;
+		//Attribute<vec3i> tij;
+		Attribute<vec3> tij;
 
         vector<index_t> v2e;
         vector<vector<index_t> > v2eopp;

src/lib/exploragram/hexdom/basic.cpp

@@ -41,6 +41,10 @@
 #include <geogram/basic/file_system.h>
 #include <geogram/basic/string.h>
 
+GEO::FF_param GEO::HexdomParam::FF;
+GEO::FF_param::FF_param() { rigid_border = true; }
+
+
 
 std::string plop_file(const char* file_in, int line) {
     std::string file(file_in);
@@ -49,4 +53,3 @@ std::string plop_file(const char* file_in, int line) {
 }
 
 
-

src/lib/exploragram/hexdom/basic.h

@@ -54,6 +54,19 @@
 
 namespace GEO {
 
+	struct FF_param {
+		FF_param();
+		bool rigid_border;
+	};
+	struct HexdomParam {
+		static FF_param  FF;
+	};
+
+	template<class T> void min_equal(T& A, T B) { if (A > B) A = B; }
+	template<class T> void max_equal(T& A, T B) { if (A < B) A = B; }
+
+	inline double nint(double x) { return floor(x + .5); }
+
     inline index_t next_mod(index_t i, index_t imax) {
 	return (i + 1) % imax;
     }
@@ -153,8 +166,29 @@ namespace GEO {
 	while (id >= data.size()) id -= data.size();
 	return data[id];
     }
+
+
+
+
+struct EXPLORAGRAM_API BBox1 {
+	BBox1() { min = 1e20; max = -1e20; }
+	double length() { return max - min; }
+	bool intersect(const BBox1& b) const { return contains(b.min) || contains(b.max) || b.contains(min) || b.contains(max); }
+	bool contains(const double& v) const { return v > min && v < max; }
+	bool is_null() const;
+	void add(const BBox1& b);
+	void add(const double& P) { min_equal(min, P); max_equal(max, P);}
+	void dilate(double eps) { min -= eps; max += eps; }
+	double bary() const { return (max + min) / 2.; }
+
+	double min;
+	double max;
+};
 }
 
+
+
+
 /**************** debugging and logging ***************************/
 
 #define reach(x){GEO::Logger::out("HexDom") <<"\n=========    mark    ===>    "<<#x<<"    =============\n"<< std::endl; }

src/lib/exploragram/hexdom/cavity.cpp

@@ -29,16 +29,9 @@ namespace GEO {
 		if (hex->cells.nb() == 0) return;
 		double eps = 1e-3*get_cell_average_edge_size(hex);
 
-		Attribute<index_t> hexvid(quad->vertices.attributes(), "hexvid");
-		FOR(v, quad->vertices.nb()) hexvid[v] = NOT_AN_ID;
-
-
-
-
 		// add hex surface to quadt
 		index_t off_v = quad->vertices.create_vertices(hex->vertices.nb());
 		FOR(v, hex->vertices.nb()) X(quad)[off_v + v] = X(hex)[v];
-		FOR(v, hex->vertices.nb()) hexvid[off_v + v] = v;
 
 		index_t off_f = quad->facets.create_facets(hex->facets.nb(), 4);
 		FOR(f, hex->facets.nb()) FOR(fc, 4)
@@ -59,7 +52,6 @@ namespace GEO {
 				index_t nearest = NN->get_nearest_neighbor(X(quad)[v].data());
 				if ((X(quad)[v] - X(quad)[nearest]).length2() == 0) {
 					old2new[v] = nearest;
-					hexvid[nearest] = hexvid[v];
 				}
 			}
 

src/lib/exploragram/hexdom/extra_connectivity.cpp

@@ -105,6 +105,90 @@ namespace GEO {
 
 
 
+
+
+
+
+
+
+		FacetsExtraConnectivityWithInvalidFacets::FacetsExtraConnectivityWithInvalidFacets(Mesh * p_m) {
+			m = p_m;
+			facet_is_valid.bind(m->facets.attributes(), "is_valid");
+			FOR(f, m->facets.nb()) facet_is_valid[f] = true;
+			reset();
+		}
+
+		void FacetsExtraConnectivityWithInvalidFacets::reset() {
+			//plop(m->facets.nb());
+			index_t nbc = m->facet_corners.nb();
+			c2f.resize(nbc,NOT_AN_ID);
+			c2c.resize(nbc, NOT_AN_ID);
+			v2c.resize(m->vertices.nb(), NOT_AN_ID);
+			FOR(f, m->facets.nb()) FOR(fc, m->facets.nb_corners(f)) {
+				index_t c = m->facets.corner(f, fc);
+				if (facet_is_valid[f]) { // everything is NOT_AN_ID for invalid facet and associated corner 
+					c2f[c] = f;
+					c2c[c] = c;
+					v2c[m->facets.vertex(f, fc)] = c;
+				}
+			}
+			FOR(f, m->facets.nb()) {
+				if (!facet_is_valid[f]) continue;
+				FOR(fc, m->facets.nb_corners(f)) {
+					index_t c = m->facets.corner(f, fc);
+					c2c[c] = v2c[m->facets.vertex(f, fc)];
+					v2c[m->facets.vertex(f, fc)] = c;
+				}
+			}
+		}
+
+		index_t FacetsExtraConnectivityWithInvalidFacets::org(index_t corner_id) { return m->facet_corners.vertex(corner_id); }
+		index_t FacetsExtraConnectivityWithInvalidFacets::dest(index_t corner_id) { return m->facet_corners.vertex(next(corner_id)); }
+
+		index_t FacetsExtraConnectivityWithInvalidFacets::opposite(index_t corner_id) {
+			index_t cir = corner_id;
+			index_t result = NOT_AN_ID; // not found
+			do {
+				index_t candidate = prev(cir);
+				if ((org(candidate) == dest(corner_id)) && (dest(candidate) == org(corner_id))) {
+					if (result == NOT_AN_ID) result = candidate;
+					else return NOT_AN_ID; // found more than one
+				}
+				if (cir != corner_id && dest(corner_id) == dest(cir))
+					return NOT_AN_ID; // the edge is non manifold
+				cir = c2c[cir];
+			} while (cir != corner_id);
+			return result;
+		}
+		//index_t FacetsExtraConnectivityWithInvalidFacets::next_around_vertex(index_t  cir) { return opposite(prev(cir)); }
+
+		index_t FacetsExtraConnectivityWithInvalidFacets::facet(index_t corner_id) { return c2f[corner_id]; }
+		index_t FacetsExtraConnectivityWithInvalidFacets::local_id(index_t corner_id) { return corner_id - m->facets.corners_begin(c2f[corner_id]); }
+
+		index_t FacetsExtraConnectivityWithInvalidFacets::next(index_t corner_id) {
+			index_t fc = local_id(corner_id);
+			index_t offset = corner_id - fc;
+			return offset + next_mod(fc, m->facets.nb_corners(c2f[corner_id]));
+		}
+		index_t FacetsExtraConnectivityWithInvalidFacets::prev(index_t corner_id) {
+			index_t fc = local_id(corner_id);
+			index_t offset = corner_id - fc;
+			return offset + prev_mod(fc, m->facets.nb_corners(c2f[corner_id]));
+		}
+
+		vec3 FacetsExtraConnectivityWithInvalidFacets::geom(index_t corner_id) {
+			return X(m)[dest(corner_id)] - X(m)[org(corner_id)];
+		}
+
+
+
+
+
+
+
+
+
+
         void halfedge_manip_example(Mesh* m){
                 FacetsExtraConnectivity fec(m);
                 FOR(c, m->facet_corners.nb()){// we can directly loop over each halfedges