Code/Sandbox/Plugins/LodGeneratorPlugin/Util/ABF.h

   1 #pragma once
   2
   3 #include "TopologyGraph.h"
   4 #include "Vector.h"
   5 #include "SparseMatrix.h"
   6
   7 namespace LODGenerator
   8 {
   9         class CComponent;
  10     class CABF
  11         {
  12     public:
  13         CABF();
  14
  15         bool do_parameterize_disc(CTopologyGraph* map);
  16                 bool parameterize_disc(CComponent* disc);
  17
  18         // _____________ Tuning parameters _____________________
  19         // Most users will not need to change them.
  20
  21
  22         // Enables Zayer et.al's modified wheel condition
  23         void set_log_mode(bool x) { log_mode_ = x; }
  24         bool get_log_mode() const { return log_mode_; }
  25
  26         // Tolerance for the norm of the gradient
  27         double get_newton_tolf() const { return newton_tolf_; }
  28         void set_newton_tolf(double x) { newton_tolf_ = x; }
  29
  30         // Tolerance for the norm of delta (step vector)
  31         double get_newton_tolx() const { return newton_tolx_; }
  32         void set_newton_tolx(double x) { newton_tolx_ = x; }
  33
  34         // Maximum number of newton iterations (outer loop)
  35         int get_max_newton_iters() const { return max_newton_iter_; }
  36         void set_max_newton_iters(int n) { max_newton_iter_ = n; }
  37
  38         double get_step_length_factor() const { return step_length_factor_; }
  39         void set_step_length_factor(double x) { step_length_factor_ = x; }
  40
  41         public:
  42                 CTopologyGraph* component_to_map(CComponent* comp);
  43
  44     protected:
  45
  46         void allocate_variables();
  47         void deallocate_variables();
  48
  49         // takes into account polygonal facets
  50         static int nb_triangles(CTopologyGraph* map);
  51         static int nb_interior_vertices(CTopologyGraph* map);
  52
  53         // ------------------ Main driver & Utilities
  54
  55         void enumerate_angles();
  56         void compute_beta();
  57         bool solve_angles();
  58
  59         // ------------------ Jacobian of the constraints:
  60
  61         // C1 : alpha_i1 + alpha_i2 + alpha_i3 - PI = 0
  62         //        where Ti = (i1,i2,i3)
  63         // C2 : Sum alpha_ij - 2.PI = 0
  64         //        for alpha_ij incident to vertex i
  65         // C3 : Prod sin(gamma_ij) - Prod sin(gamma'_ij) = 0
  66         //        for gamma_ij and gamma'_ij opposite incident to vertex i
  67         //
  68
  69         void compute_product_sin_angles(
  70            CVertex* v, double& prod_prev_sin, double& prod_next_sin
  71         );
  72
  73         void compute_sum_log_sin_angles(
  74            CVertex* v, double& sum_prev_sin, double& sum_next_sin
  75         );
  76
  77         // JC1 is implicit.
  78         void add_JC2();
  79         void add_JC3();
  80
  81         // ------------------  Right hand side: gradients
  82
  83         // Gradient of the quadratic form
  84         void sub_grad_F();
  85
  86         // For each triangle: sum angles - PI
  87         void sub_grad_C1();
  88
  89         // For each vertex: sum incident angles - 2.PI
  90         void sub_grad_C2();
  91
  92         // For each vertex: prod sin(next angle) - prod sin(prev angle)
  93         void sub_grad_C3();
  94
  95         // ------------------ Solver for one iteration
  96
  97         // computes dalpha, dlambda1 and dlambda2
  98         void solve_current_iteration();
  99
 100         // ------------------ Convergence control
 101
 102         // increases the weights of negative angles.
 103         double compute_step_length_and_update_weights();
 104
 105         // returns the norm of the gradient (quadratic form + cnstrs).
 106         double errf() const;
 107
 108         // returns the norm of the step vector.
 109         double compute_errx_and_update_x(double s);
 110
 111         // ------------------ Utilities
 112
 113         // y += J.D.x
 114         void add_J_D_x(CVectorD& y, const CSparseMatrix& J, CVectorD& D, CVectorD& x);
 115
 116         // M += J.D.J^t
 117         void add_J_D_Jt(CSparseMatrix& M, const CSparseMatrix& J, CVectorD& D);
 118
 119         // M -= J.D.J^t
 120         void sub_J_D_Jt(CSparseMatrix& M, const CSparseMatrix& J, CVectorD& D);
 121
 122     private:
 123
 124         CTopologyGraph* map_;
 125
 126         int nf_;      // Number of facets
 127         int nalpha_;  // Number of angles (= 3.nf)
 128         int nint_;    // Number of interior nodes
 129         int nlambda_; // Number of constraints (= nf+2.nint)
 130         int ntot_;    // Total number of unknowns (= nalpha + nlamda)
 131
 132                 std::map<CHalfedge*,int> angle_index_;
 133                 std::map<CHalfedge*,double> angle_angle_;
 134
 135         // ------ ABF variables & Lagrange multipliers ----------------
 136         CVectorD alpha_;    // Unknown angles. size = nalpha
 137         CVectorD lambda_;   // Lagrange multipliers. size = nlambda
 138         CVectorD beta_;     // Optimum angles.  size = nalpha
 139         CVectorD w_;        // Weights.         size = nalpha
 140
 141         // ------ Step vectors ----------------------------------------
 142         CVectorD dalpha_  ; // size = nalpha; angles
 143         CVectorD dlambda1_; // size = nf    ; C1 part
 144         CVectorD dlambda2_; // size = 2.nint; C2 and C3 part
 145
 146         // ------ Right-hand side ( - gradients ) ---------------------
 147         CVectorD b1_;      // size = nalpha
 148         CVectorD b2_;      // size = nlambda
 149
 150         // ------ Jacobian of the constraints -------------------------
 151         // J1 (Jacobian of constraint 1) is not stored, it is implicit
 152         CSparseMatrix J2_; // size = 2.nint * nalpha
 153
 154
 155         // ------ ABF++ variables -------------------------------------
 156         CVectorD Delta_inv_     ; // size = nalpha
 157         CVectorD Delta_star_inv_; // size = nf;
 158         CSparseMatrix J_star_  ; // size = 2.nint * nf
 159         CVectorD b1_star_       ; // size = nf
 160         CVectorD b2_star_       ; // size = 2.nint
 161
 162         // ------ Final linear system ---------------------------------
 163         CSparseMatrix M_       ; // size = 2.nint * 2.nint
 164         CVectorD       r_       ; // size = 2.nint
 165
 166         CVectorI perm_;   // column permutations (used by SuperLU)
 167
 168         //_______________ constants, tuning parameters ____________________
 169
 170         const double epsilon_;    // threshold for small angles
 171
 172         double newton_tolf_;   // threshold for gradient norm (rhs)
 173         double newton_tolx_;  // threshold for delta norm
 174         int max_newton_iter_;
 175         const double positive_angle_ro_;
 176         double step_length_factor_;
 177
 178         bool low_mem_; // If set, tries to reduce memory footprint
 179
 180         bool log_mode_;
 181     };
 182
 183     //_______________________________________________________________________
 184 }
 185