[geogram] Update to 1.6.3

CMakeLists.txt

@@ -39,7 +39,7 @@ include(cmake/geogram.cmake)
 
 set(VORPALINE_VERSION_MAJOR 1)
 set(VORPALINE_VERSION_MINOR 6)
-set(VORPALINE_VERSION_PATCH 2)
+set(VORPALINE_VERSION_PATCH 3)
 set(VORPALINE_VERSION ${VORPALINE_VERSION_MAJOR}.${VORPALINE_VERSION_MINOR}.${VORPALINE_VERSION_PATCH})
 
 set(VORPALINE_INCLUDE_SUBPATH geogram${VORPALINE_VERSION_MAJOR})

cmake/geogram.cmake

@@ -2,7 +2,7 @@
 if(NOT DEFINED GEOGRAM_SOURCE_DIR)
    set(GEOGRAM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 endif()
-message(STATUS "GEOGRAM_SOURCE_DIR: ${GEOGRAM_SOURCE_DIR}")
+
 if(EXISTS ${GEOGRAM_SOURCE_DIR}/CMakeOptions.txt)
    message(STATUS "Using local options file: ${GEOGRAM_SOURCE_DIR}/CMakeOptions.txt")
    include(${GEOGRAM_SOURCE_DIR}/CMakeOptions.txt)

src/examples/graphics/GLUP_basic_example/main.cpp

@@ -852,7 +852,7 @@ namespace {
 
 static void overlay() {
     // Displays ImGui demo window.
-    ImGui::SetNextWindowPos(ImVec2(650, 20), ImGuiSetCond_FirstUseEver);
+    ImGui::SetNextWindowPos(ImVec2(650, 20), ImGuiCond_FirstUseEver);
     ImGui::ShowDemoWindow();
 }
 

src/examples/graphics/demo_Delaunay2d/main.cpp

@@ -447,8 +447,8 @@ namespace {
     void overlay() {
         static int border_shape = 0;
 
-        ImGui::SetNextWindowPos(ImVec2(20, 20), ImGuiSetCond_Once);
-        ImGui::SetNextWindowSize(ImVec2(180, 320), ImGuiSetCond_Once);
+        ImGui::SetNextWindowPos(ImVec2(20, 20), ImGuiCond_Once);
+        ImGui::SetNextWindowSize(ImVec2(180, 320), ImGuiCond_Once);
 
         ImGui::Begin("Delaunay");
 

src/lib/exploragram/optimal_transport/optimal_transport.cpp

@@ -133,7 +133,7 @@ namespace GEO {
 	linesearch_maxiter_ = 100;
 	linesearch_init_iter_ = 0;
 
-	use_direct_solver_ = false;
+	linear_solver_ = OT_PRECG;
 
 	nb_air_particles_ = 0;
 	air_particles_ = nil;
@@ -823,12 +823,18 @@ namespace GEO {
     void OptimalTransportMap::new_linear_system(index_t n, double* x) {
         nlNewContext();
 
-        if(use_direct_solver_) {
-            use_direct_solver_ = (
-		use_direct_solver_ && (nlInitExtension("SUPERLU") == NL_TRUE)
-	    );
-	    if(!use_direct_solver_) {
-		Logger::warn("OTM") << "Could not initialize SUPERLU OpenNL extension"
+        if(linear_solver_ != OT_PRECG) {
+	    if(linear_solver_ == OT_SUPERLU) {
+		if(nlInitExtension("SUPERLU") != NL_TRUE) {
+		    linear_solver_ = OT_PRECG;
+		}
+	    } else if(linear_solver_ == OT_CHOLMOD) {
+		if(nlInitExtension("CHOLMOD") != NL_TRUE) {
+		    linear_solver_ = OT_PRECG;
+		}
+	    }
+	    if(linear_solver_ == OT_PRECG) {
+		Logger::warn("OTM") << "Could not initialize OpenNL extension"
 				    << std::endl;
 		Logger::warn("OTM") << "Falling back to conjugate gradient"
 				    << std::endl;
@@ -840,15 +846,21 @@ namespace GEO {
 	}
 
         nlSolverParameteri(NL_NB_VARIABLES, NLint(n));
-        if(use_direct_solver_) {
-            nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT);
-        } else {
-            nlSolverParameteri(NL_SOLVER, NL_CG);
-            nlSolverParameteri(NL_PRECONDITIONER, NL_PRECOND_JACOBI);
-            nlSolverParameteri(NL_SYMMETRIC, NL_TRUE);
-            nlSolverParameterd(NL_THRESHOLD, linsolve_epsilon_);
-            nlSolverParameteri(NL_MAX_ITERATIONS, NLint(linsolve_maxiter_));                
-        }
+	switch(linear_solver_) {
+	    case OT_PRECG:
+		nlSolverParameteri(NL_SOLVER, NL_CG);
+		nlSolverParameteri(NL_PRECONDITIONER, NL_PRECOND_JACOBI);
+		nlSolverParameteri(NL_SYMMETRIC, NL_TRUE);
+		nlSolverParameterd(NL_THRESHOLD, linsolve_epsilon_);
+		nlSolverParameteri(NL_MAX_ITERATIONS, NLint(linsolve_maxiter_));                
+		break;
+	    case OT_SUPERLU:
+		nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT);		
+		break;
+	    case OT_CHOLMOD:
+		nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);		
+		break;
+	}
 
 	nlEnable(NL_VARIABLES_BUFFER);
         nlBegin(NL_SYSTEM);

src/lib/exploragram/optimal_transport/optimal_transport.h

@@ -41,6 +41,19 @@
 #define H_EXPLORAGRAM_OPTIMAL_TRANSPORT_OPTIMAL_TRANSPORT_H
 
 #include <exploragram/basic/common.h>
+
+namespace GEO {
+    /**
+     * \brief Specifies the linear solver to be used
+     *  with OptimalTransport.
+     */
+    enum OTLinearSolver {
+	OT_PRECG, OT_SUPERLU, OT_CHOLMOD
+    };
+}
+
+#ifndef GOMGEN
+
 #include <geogram/mesh/mesh.h>
 #include <geogram/voronoi/RVD.h>
 #include <geogram/delaunay/delaunay.h>
@@ -55,6 +68,7 @@
 namespace GEO {
     class CentroidalVoronoiTesselation;
 
+
     /**
      * \brief Computes semi-discrete optimal transport maps.
      * \details Computes an optimal transport map between two
@@ -255,16 +269,13 @@ namespace GEO {
 
 	/**
 	 * \brief Specifies whether a direct solver should be used.
-	 * \param[in] x true if a direct solver should be used, false
-	 *  otherwise. Default is false.
-	 * \details The used direct solver is SUPERLU (if OpenNL extension
-	 *  is available) and the iterative solver the Jacobi-preconditioned
-	 *  conjugate gradient. The direct solver is recommended only for
+	 * \param[in] solver one of OT_PRECG (default), OT_SUPERLU, OT_CHOLMOD.
+	 * \details The direct solvers (OT_SUPERLU, OT_CHOLMOD) are recommended only for
 	 *  surfacic data, since the sparse factors become not so sparse when
 	 *  volumetric meshes are considered.
 	 */
-	void set_use_direct_solver(bool x) {
-	    use_direct_solver_ = x;
+	void set_linear_solver(OTLinearSolver solver) {
+	    linear_solver_ = solver;
 	}
 
         /**
@@ -773,8 +784,8 @@ namespace GEO {
 	/** \brief starting number of steplength divisions */
 	index_t linesearch_init_iter_;
 
-	/** \brief if true, uses a direct solver (SUPERLU) */
-	bool use_direct_solver_;
+	/** \brief one of OT_PRECG, OT_SUPERLU, OT_CHOLMOD. */
+	OTLinearSolver linear_solver_;
 
 	/** \brief if set, pointer to the air particles. */
 	const double* air_particles_;
@@ -803,3 +814,5 @@ namespace GEO {
 }
 
 #endif
+
+#endif

src/lib/exploragram/optimal_transport/optimal_transport_2d.cpp

@@ -682,8 +682,17 @@ namespace GEO {
 	    false
 	);
 
-        //OTM.set_regularization(1e-3);
-	//OTM.set_use_direct_solver(true);
+	static bool initialized = false;
+	static bool has_cholmod = false;
+	if(!initialized) {
+	    initialized = true;
+	    has_cholmod = (nlInitExtension("CHOLMOD") == NL_TRUE);
+	}
+	if(has_cholmod) {
+	    OTM.set_regularization(1e-3);
+	    OTM.set_linear_solver(OT_CHOLMOD);
+	}
+
         OTM.set_Newton(true);
 	OTM.set_air_particles(
 	    nb_air_particles, air_particles, air_particles_stride, air_fraction

src/lib/exploragram/optimal_transport/optimal_transport_on_surface.cpp

@@ -44,6 +44,8 @@
 #include <geogram/basic/geometry_nd.h>
 #include <geogram/basic/stopwatch.h>
 
+#include <geogram/bibliography/bibliography.h>
+
 namespace {
     using namespace GEO;
 
@@ -87,7 +89,6 @@ namespace {
 		return;
 	    }
 
-	    geo_argused(t);
 	    double m, mgx, mgy, mgz;
 	    compute_m_and_mg(P, m, mgx, mgy, mgz);
 
@@ -116,7 +117,7 @@ namespace {
 
 	    if(Newton_step_) {
 		// Spinlocks are managed internally by update_Hessian().
-		update_Hessian(P, v);
+		update_Hessian(P, v, t);
 	    }
 
 
@@ -175,15 +176,39 @@ namespace {
 	    }
 	}	
 
+	/**
+	 * \brief Gets the unit normal vector to a facet.
+	 */
+	vec3 facet_normal(index_t t) const {
+	    Mesh& M = OTM_->mesh();
+	    index_t i = M.facets.vertex(t,0);
+	    index_t j = M.facets.vertex(t,1);
+	    index_t k = M.facets.vertex(t,2);	    
+	    const double* p0 = M.vertices.point_ptr(i);
+	    const double* p1 = M.vertices.point_ptr(j);
+	    const double* p2 = M.vertices.point_ptr(k);
+	    double x1 = p1[0] - p0[0];
+	    double y1 = p1[1] - p0[1];
+	    double z1 = p1[2] - p0[2];
+	    double x2 = p2[0] - p0[0];
+	    double y2 = p2[1] - p0[1];
+	    double z2 = p2[2] - p0[2];
+	    double x = y1*z2 - z1*y2;
+	    double y = z1*x2 - x1*z2;
+	    double z = x1*y2 - y1*x2;
+	    return vec3(x,y,z);
+	}
+
 	/**
 	 * \brief Updates the Hessian according to the current intersection
 	 *  polygon.
 	 * \param[in] P a const reference to the current intersection polygon.
 	 * \param[in] i the current seed
+	 * \param[in] t the current mesh facet
 	 */
-	void update_Hessian(
-	    const GEOGen::Polygon& P, index_t i
-	) const {
+	void update_Hessian(const GEOGen::Polygon& P, index_t i, index_t t) const {
+
+	    vec3 N = normalize(facet_normal(t));
 
 	    // The coefficient of the Hessian associated to a pair of
 	    // adjacent cells Lag(i),Lag(j) is :
@@ -202,9 +227,13 @@ namespace {
 		index_t j = index_t(P.vertex(k2).adjacent_seed());
 		if(j != index_t(-1)) {
 		    const double* pj = OTM_->point_ptr(j);
-		    double hij =
-			edge_mass(P.vertex(k1), P.vertex(k2)) /
-			(2.0 * GEO::Geom::distance(pi,pj,3)) ;
+
+		    vec3 pij(pj[0] - pi[0],pj[1] - pi[1],pj[2] - pi[2]);
+		    // Remove from pij the part that is normal to the current polygon.
+		    pij -= dot(pij,N)*N;
+		    double lij = length(pij);
+
+		    double hij = edge_mass(P.vertex(k1), P.vertex(k2)) / (2.0 * lij);
 
 		    // -hij because we maximize F <=> minimize -F
 		    if(hij != 0.0) {
@@ -369,9 +398,10 @@ namespace GEO {
 	    mesh,
 	    (delaunay == "default") ? "BPOW" : delaunay,
 	    BRIO
-        ) {
+    ) {
 	callback_ = new SurfaceOTMPolygonCallback(this);
 	total_mass_ = total_mesh_mass();
+        geo_cite("DBLP:journals/corr/MerigotMT17");
     }
 
     OptimalTransportMapOnSurface::~OptimalTransportMapOnSurface() {
@@ -491,8 +521,17 @@ namespace GEO {
 	    false
 	);
 
-	//OTM.set_regularization(1e-3);
-	//OTM.set_use_direct_solver(true);
+	static bool initialized = false;
+	static bool has_cholmod = false;
+	if(!initialized) {
+	    initialized = true;
+	    has_cholmod = (nlInitExtension("CHOLMOD") == NL_TRUE);
+	}
+	if(has_cholmod) {
+	    OTM.set_regularization(1e-3);
+	    OTM.set_linear_solver(OT_CHOLMOD);
+	}
+
         OTM.set_Newton(true);
         OTM.set_points(nb_points, points);
         OTM.set_epsilon(0.01);

src/lib/geogram/NL/nl.h

@@ -733,7 +733,7 @@ typedef void* NLContext;
  *  alternative solvers, such as SuperLU (sparse direct solver) and CNC 
  *  (iterative solver on the GPU). This function tests whether an extension 
  *  is supported, and initializes what needs to be initialized in the extention.
- * \retval NL_TRUE if the extension is supported and could be sucessfully 
+ * \retval NL_TRUE if the extension is supported and could be successfully 
  *   initialized
  * \retval NL_FALSE otherwise
  */

src/lib/geogram/NL/nl_arpack.h

@@ -56,7 +56,7 @@
  * \brief Initializes the ARPACK extension
  * \details This dynamically loads the SuperLU 
  *  library available in the system (if available) and
- *  retreives the symbols in there. It supports SuperLU 5.x.
+ *  retrieves the symbols in there. It supports SuperLU 5.x.
  * \retval NL_TRUE if ARPACK could be successfully
  *   dynamically loaded and all functions could be
  *   found in it.

src/lib/geogram/NL/nl_cholmod.h

@@ -73,7 +73,7 @@ NLAPI NLMatrix NLAPIENTRY nlMatrixFactorize_CHOLMOD(
  * \brief Initializes the CHOLMOD extension
  * \details This dynamically loads the Cholmod 
  *  library available in the system (if available) and
- *  retreives the symbols in there. 
+ *  retrieves the symbols in there. 
  * \retval NL_TRUE if CHOLMOD could be successfully
  *   dynamically loaded and all functions could be
  *   found in it.

src/lib/geogram/NL/nl_cuda.h

@@ -59,7 +59,7 @@
  * \brief Initializes the CUDA extension
  * \details This dynamically loads the CUDA
  *  libraries available in the system (if available) and
- *  retreives the symbols in there.
+ *  retrieves the symbols in there.
  * \retval NL_TRUE if CUDA could be successfully
  *   dynamically loaded and all functions could be
  *   found in it.

src/lib/geogram/NL/nl_matrix.c

@@ -1035,7 +1035,7 @@ typedef struct {
     NLMatrix M;
 
     /**
-     * \brief NL_TRUE if memory ownership was transfered,
+     * \brief NL_TRUE if memory ownership was transferred,
      *  NL_FALSE otherwise.
      */
     NLboolean owns_M;
@@ -1046,7 +1046,7 @@ typedef struct {
     NLMatrix N;
 
     /**
-     * \brief NL_TRUE if memory ownership was transfered,
+     * \brief NL_TRUE if memory ownership was transferred,
      *  NL_FALSE otherwise.
      */
     NLboolean owns_N;

src/lib/geogram/NL/nl_matrix.h

@@ -123,7 +123,7 @@ NLAPI void NLAPIENTRY nlDeleteMatrix(NLMatrix M);
 /**
  * \brief Computes a matrix x vector product
  * \param[in] M the matrix
- * \param[in] x the vector to be multipled by the matrix
+ * \param[in] x the vector to be multiplied by the matrix
  * \param[out] y the result
  */
 NLAPI void NLAPIENTRY nlMultMatrixVector(
@@ -686,7 +686,7 @@ NLAPI NLuint NLAPIENTRY nlMatrixNNZ(NLMatrix M);
 
 /**
  * \brief Factorizes a matrix.
- * \details The corresponding extension needs to be sucessfully initialized
+ * \details The corresponding extension needs to be successfully initialized
  *  before calling this function.
  * \param[in] M the input matrix
  * \param[in] solver a direct solver, i.e., one of:

src/lib/geogram/NL/nl_superlu.h

@@ -75,7 +75,7 @@ NLAPI NLMatrix NLAPIENTRY nlMatrixFactorize_SUPERLU(
  * \brief Initializes the SUPERLU extension
  * \details This dynamically loads the SuperLU 
  *  library available in the system (if available) and
- *  retreives the symbols in there. It supports SuperLU 5.x.
+ *  retrieves the symbols in there. It supports SuperLU 5.x.
  * \retval NL_TRUE if SUPERLU could be successfully
  *   dynamically loaded and all functions could be
  *   found in it.