1 // @licstart revoropt 2 // This file is part of Revoropt, a library for the computation and 3 // optimization of restricted Voronoi diagrams. 4 // 5 // Copyright (C) 2013 Vincent Nivoliers <vincent.nivoliers@univ-lyon1.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 // @licend revoropt 11 #ifndef _REVOROPT_CVT_COMPUTER_H_ 12 #define _REVOROPT_CVT_COMPUTER_H_ 13 14 #include <stddef.h> 15 #include <Eigen/Dense> 16 17 #include <Revoropt/RVD/rvd.hpp> 18 #include <Revoropt/RVD/polygon.hpp> 19 #include <Revoropt/Mesh/all_def.hpp> 20 #include <Revoropt/Mesh/measure_def.hpp> 21 22 namespace Revoropt { 23 namespace CVT { 24 25 template<typename _Triangulation> 26 class BaseComputer {// : public RVD::Action<_Triangulation> { 27 28 public : 29 30 /* Typedefs and template arguments */ 31 typedef _Triangulation Triangulation ; 32 enum { Dim = Triangulation::VertexDim } ; 33 typedef Eigen::Matrix<double,Dim,1> Vector ; 34 BaseComputer()35 BaseComputer() : g_(NULL), fx_(0), 36 mesh_(NULL), sites_(NULL), 37 triangle_weights_(NULL), 38 mesh_area_(0), mesh_boundary_length_(0), 39 factor_(1), anisotropy_(1), 40 border_anisotropy_(1), border_weight_(0), 41 tolerance_(0) { 42 } 43 44 /* LBFGS variables to compute */ set_g(double * g)45 void set_g( double* g ) { g_ = g ; } set_fx(double val)46 void set_fx( double val ) { fx_ = val ; } get_fx()47 double get_fx() { return fx_ ; } 48 49 /* Mesh and sites */ set_mesh(const Triangulation * mesh)50 void set_mesh( const Triangulation* mesh ) { 51 mesh_ = mesh ; 52 mesh_area_ = mesh_area(mesh) ; 53 //FIXME 54 //mesh_boundary_length_ = mesh->boundary_length() ; 55 } 56 get_mesh()57 const Triangulation* get_mesh() { 58 return mesh_ ; 59 } set_sites(const double * sites)60 void set_sites( const double* sites ) { 61 sites_ = sites ; 62 } set_triangle_weights(const double * triangle_weights)63 void set_triangle_weights( const double* triangle_weights ) { 64 triangle_weights_ = triangle_weights ; 65 } 66 67 /* Parameter setting */ set_factor(double factor)68 void set_factor( double factor ) { 69 factor_ = factor ; 70 } 71 set_anisotropy(double anisotropy)72 void set_anisotropy( double anisotropy ) { 73 anisotropy_ = anisotropy ; 74 } 75 set_border_anisotropy(double border_anisotropy)76 void set_border_anisotropy( double border_anisotropy ) { 77 border_anisotropy_ = border_anisotropy ; 78 } 79 set_border_weight(double border_weight)80 void set_border_weight( double border_weight ) { 81 border_weight_ = border_weight ; 82 } 83 set_tolerance(double tolerance)84 void set_tolerance( double tolerance ) { 85 tolerance_ = tolerance ; 86 } 87 88 protected : 89 90 /* LBFGS variables to compute */ 91 double* g_ ; 92 double fx_ ; 93 94 /* Mesh and sites */ 95 const Triangulation* mesh_ ; 96 const double* sites_ ; 97 const double* triangle_weights_ ; 98 99 double mesh_area_ ; 100 double mesh_boundary_length_ ; 101 102 /* Parameter variables */ 103 104 double factor_ ; 105 double anisotropy_ ; 106 double border_anisotropy_ ; 107 double border_weight_ ; 108 double tolerance_ ; 109 } ; 110 111 /** Direct computer **/ 112 113 /* Compute the value and gradient of the direct CVT objective function when 114 * passed to the corresponding RVD computer. The direct CVT objective function 115 * optimizes a set of samples to sample uniformly a mesh. */ 116 117 template<typename _Triangulation> 118 class DirectComputer : public BaseComputer<_Triangulation> { 119 120 public : 121 122 /* Typedefs and template arguments */ 123 typedef _Triangulation Triangulation ; 124 enum { Dim = Triangulation::VertexDim } ; 125 typedef Eigen::Matrix<double,Dim,1> Vector ; 126 typedef BaseComputer<Triangulation> Base ; 127 128 typedef RVDEdge<Triangulation> Edge ; 129 DirectComputer()130 DirectComputer() : BaseComputer<Triangulation>() {} 131 operator ()(unsigned int site,unsigned int face,const RVDPolygon<Triangulation> & polygon)132 void operator()( unsigned int site, 133 unsigned int face, 134 const RVDPolygon<Triangulation>& polygon 135 ) { 136 //number of vertices of the polygon 137 const unsigned int size = polygon.size() ; 138 139 if(size < 3) return ; 140 141 //send triangles to the triangle backend 142 for(unsigned int i=1; i<size-1; ++i) { 143 triangle_compute( site, 144 face, 145 polygon[0], 146 polygon[i], 147 polygon[i+1] 148 ) ; 149 } 150 151 //std::cout << "after triangle " << fx_ << std::endl ; 152 153 //compute a point inside the polygon for boundary normal computations 154 Vector bar ; 155 bar.setZero() ; 156 for(unsigned int i=0; i<size; ++i) { 157 bar += polygon[i].vertex ; 158 } 159 bar /= size ; 160 161 //send bisector edges to the edge backend 162 for(unsigned int i=0; i<size; ++i) { 163 if(polygon[i].config == Edge::BISECTOR_EDGE) { 164 edge_compute( site, 165 face, 166 polygon[i], 167 polygon[(i+1)%size], 168 bar 169 ) ; 170 } 171 172 /* FIXME 173 if(border_weight_ != 0 && mesh_boundary_length_ != 0) { 174 //one dimensional CVT of the boundary of the mesh 175 if(polygon[i].config == RVD::FACE_EDGE) { 176 //get edge index 177 const unsigned int edge_index = polygon[i].combinatorics ; 178 //get the neighbours of the triangle 179 const unsigned int* triangle_neighbours 180 = mesh_->face_neighbours(triangle) ; 181 //check whether the edge is a boundary edge 182 if( triangle_neighbours[edge_index] 183 == Triangulation::NO_NEIGHBOUR 184 ) { 185 border_compute(site, polygon[i], polygon[(i+1)%size]) ; 186 } 187 } 188 } 189 */ 190 } 191 192 //std::cout << "after edges " << fx_ << std::endl ; 193 } 194 195 private : 196 197 using Base::factor_ ; 198 using Base::anisotropy_ ; 199 using Base::tolerance_ ; 200 using Base::mesh_area_ ; 201 using Base::fx_ ; 202 using Base::g_ ; 203 using Base::sites_ ; 204 using Base::mesh_ ; 205 using Base::triangle_weights_ ; 206 triangle_compute(const unsigned int site,const unsigned int triangle,const RVDEdge<Triangulation> & e1,const RVDEdge<Triangulation> & e2,const RVDEdge<Triangulation> & e3)207 void triangle_compute( const unsigned int site, 208 const unsigned int triangle, 209 const RVDEdge<Triangulation>& e1, 210 const RVDEdge<Triangulation>& e2, 211 const RVDEdge<Triangulation>& e3 212 ) { 213 //vertices of the triangle 214 const Vector& c1 = e1.vertex ; 215 const Vector& c2 = e2.vertex ; 216 const Vector& c3 = e3.vertex ; 217 218 //Triangle basis and area 219 const Vector base = (c2-c1) ; 220 const double base_len = base.norm() ; 221 if(base_len <= tolerance_) return ; 222 const Vector tb1 = base/base_len ; 223 224 Vector height = (c3-c1) ; 225 height -= height.dot(tb1)*tb1 ; 226 const double height_len = height.norm() ; 227 if(height_len <= tolerance_) return ; 228 const Vector tb2 = height/height_len ; 229 230 //triangle weighting 231 const double tweight = triangle_weights_ == NULL ? 232 1 : 233 triangle_weights_[triangle] ; 234 const double area = 0.5*base_len*height_len*tweight ; 235 236 if(area <= tolerance_) return ; 237 238 //site position 239 Eigen::Map<const Vector> v(sites_ + Dim*site) ; 240 241 //function value 242 243 //site vectors 244 Vector sv[3] ; 245 sv[0] = v-c1 ; 246 sv[1] = v-c2 ; 247 sv[2] = v-c3 ; 248 249 //compute value 250 double local_fx = 0; 251 for (int i = 0; i < 3; ++i) { 252 for (int j = 0; j <= i; ++j) { 253 local_fx += sv[i].dot(sv[j]) ; 254 } 255 } 256 local_fx /= 6 ; 257 258 //compute gradient 259 260 //centroid vector 261 Vector g_vect = (c1+c2+c3)/3. ; 262 g_vect = (v-g_vect) ; 263 264 //anisotropy 265 if(anisotropy_ > 1) { 266 //normal component of the face-site vectors 267 Vector normal_component = v-c1 ; 268 normal_component -= normal_component.dot(tb1)*tb1 ; 269 normal_component -= normal_component.dot(tb2)*tb2 ; 270 ////FIXME keep component in normal space ? 271 //for(int i = 3; i < Dim; ++i) { 272 // normal_component[i] = 0 ; 273 //} 274 275 //value modification 276 local_fx += (anisotropy_*anisotropy_-1) 277 * normal_component.dot(normal_component) ; 278 279 //gradient modification 280 g_vect += (anisotropy_*anisotropy_-1) 281 * normal_component ; 282 } 283 284 //store result, normalize by mesh area 285 fx_ += (factor_/mesh_area_)*area*local_fx ; 286 Eigen::Map<Vector> grad(g_ + Dim*site) ; 287 grad += factor_*area*2*g_vect/mesh_area_ ; 288 } 289 edge_compute(const unsigned int site,const unsigned int face,const RVDEdge<Triangulation> & e1,const RVDEdge<Triangulation> & e2,const Vector & inner_point)290 void edge_compute( const unsigned int site, 291 const unsigned int face, 292 const RVDEdge<Triangulation>& e1, 293 const RVDEdge<Triangulation>& e2, 294 const Vector& inner_point 295 ) { 296 if(anisotropy_<=1) return ; 297 298 //bisector containing the edge 299 const unsigned int neighbour = e1.combinatorics ; 300 301 //handling edges only once 302 if (neighbour<site) return ; 303 304 //sites 305 Eigen::Map<const Vector> v1(sites_ + Dim*site) ; 306 Eigen::Map<const Vector> v2(sites_ + Dim*neighbour) ; 307 308 //vertices of the edge 309 const Vector& c1 = e1.vertex ; 310 const Vector& c2 = e2.vertex ; 311 312 //edge vector 313 Vector edge = c2-c1 ; 314 double edge_len = edge.norm() ; 315 if(edge_len <= tolerance_) return ; 316 edge /= edge_len ; 317 318 //triangle weighting 319 const double tweight = triangle_weights_ == NULL ? 320 1 : 321 triangle_weights_[face] ; 322 edge_len *= tweight ; 323 324 //edge normal 325 Vector edge_normal = c1-inner_point ; 326 edge_normal -= edge_normal.dot(edge)*edge ; 327 const double edge_normal_len = edge_normal.norm() ; 328 if(edge_normal_len <= tolerance_) return ; 329 edge_normal /= edge_normal_len ; 330 331 const double boundary_speed_denom = edge_normal.dot(v1-v2) ; 332 const double abs_bsd = boundary_speed_denom < 0 ? 333 - boundary_speed_denom : 334 boundary_speed_denom ; 335 if(abs_bsd <= tolerance_) return ; 336 337 //site vector normal component 338 Vector nc1 = v1-c1 ; 339 nc1 -= nc1.dot(edge)*edge ; 340 nc1 -= nc1.dot(edge_normal)*edge_normal ; 341 Vector nc2 = v2-c1 ; 342 nc2 -= nc2.dot(edge)*edge ; 343 nc2 -= nc2.dot(edge_normal)*edge_normal ; 344 ////FIXME keep component in normal space ? 345 //for(int i = 3; i < Dim; ++i) { 346 // nc1[i] = 0 ; 347 // nc2[i] = 0 ; 348 //} 349 350 //gradient common part 351 const double tmp = (anisotropy_*anisotropy_-1) 352 * (nc1.dot(nc1)-nc2.dot(nc2)) 353 / boundary_speed_denom 354 * edge_len 355 //normalize by the mesh area FIXME for triangle weights 356 * factor_/mesh_area_ ; 357 358 //edge center 359 const Vector edge_center = 0.5*(c1+c2) ; 360 361 //gradient for v1 362 Eigen::Map<Vector> grad1(g_ + Dim*site) ; 363 grad1 += tmp*(v1-edge_center) ; 364 365 //gradient for v2 366 Eigen::Map<Vector> grad2(g_ + Dim*neighbour) ; 367 grad2 -= tmp*(v2-edge_center) ; 368 } 369 370 /* FIXME 371 void border_compute( const unsigned int site, 372 const RVD::RVDEdge& e1, 373 const RVD::RVDEdge& e2 374 ) ; 375 376 void border_vertex_compute( const unsigned int site, 377 const RVD::RVDVertex& c, 378 const Vector3d& edge 379 ) ; 380 */ 381 } ; 382 383 /** Inverse computer **/ 384 385 /* Compute the value and gradient of the inverse CVT objective function when 386 * passed to the corresponding RVD computer. The inverse CVT objective function 387 * optimizes a mesh such that a given sampling regularly samples it.*/ 388 389 template<typename _Triangulation> 390 class InverseComputer : public BaseComputer<_Triangulation> { 391 392 public : 393 394 /* Typedefs and template arguments */ 395 typedef _Triangulation Triangulation ; 396 enum { Dim = Triangulation::VertexDim } ; 397 typedef Eigen::Matrix<double,Dim,1> Vector ; 398 typedef Eigen::Matrix<double,Dim,Dim> Matrix ; 399 400 typedef RVDVertex<Triangulation> Vertex ; 401 402 typedef BaseComputer<Triangulation> Base ; 403 InverseComputer()404 InverseComputer() : BaseComputer<Triangulation>() {} 405 operator ()(unsigned int site,unsigned int face,const RVDPolygon<Triangulation> & polygon)406 void operator()( unsigned int site, 407 unsigned int face, 408 const RVDPolygon<Triangulation>& polygon 409 ) { 410 //number of vertices of the polygon 411 const unsigned int size = polygon.size() ; 412 413 if(size < 3) return ; 414 415 //send triangles to the triangle backend 416 for(unsigned int i=1; i<size-1; ++i) { 417 triangle_compute(site, face, polygon[0], polygon[i], polygon[i+1]) ; 418 } 419 } 420 add_mesh_area_gradient()421 void add_mesh_area_gradient() { 422 //add the contribution of the mesh area gradient 423 mesh_area_gradient(mesh_, g_, -fx_ / mesh_area_) ; 424 } 425 finalize()426 void finalize() { 427 add_mesh_area_gradient() ; 428 } 429 430 private : 431 432 using Base::factor_ ; 433 using Base::anisotropy_ ; 434 using Base::tolerance_ ; 435 using Base::mesh_area_ ; 436 using Base::fx_ ; 437 using Base::g_ ; 438 using Base::sites_ ; 439 using Base::mesh_ ; 440 441 std::vector<double> area_gradient_ ; 442 triangle_compute(const unsigned int site,const unsigned int face,const RVDEdge<Triangulation> & e1,const RVDEdge<Triangulation> & e2,const RVDEdge<Triangulation> & e3)443 void triangle_compute( const unsigned int site, 444 const unsigned int face, 445 const RVDEdge<Triangulation>& e1, 446 const RVDEdge<Triangulation>& e2, 447 const RVDEdge<Triangulation>& e3 448 ) { 449 //vertices of the triangle 450 const RVDVertex<Triangulation>* c[3] ; 451 c[0] = &e1.vertex ; 452 c[1] = &e2.vertex ; 453 c[2] = &e3.vertex ; 454 455 //normalized edges of the triangle 456 Vector ed[3] ; 457 double ed_len[3] ; 458 for(int i=0; i<3; ++i) { 459 ed[i] = *(c[(i+2)%3]) - *(c[(i+1)%3]) ; 460 ed_len[i] = ed[i].norm() ; 461 if(ed_len[i] <= tolerance_) return ; 462 ed[i] /= ed_len[i] ; 463 } 464 465 //normalized heights of the triangle 466 Vector h[3] ; 467 double h_len[3] ; 468 for(int i=0; i<3; ++i) { 469 h[i] = (ed[(i+1)%3] - (ed[(i+1)%3].dot(ed[i])*ed[i]))*ed_len[(i+1)%3] ; 470 h_len[i] = h[i].norm() ; 471 if(h_len[i] <= tolerance_) return ; 472 h[i] /= h_len[i] ; 473 } 474 475 //area of the triangle 476 double area = 0.5*h_len[0]*ed_len[0] ; 477 if(area <= tolerance_) return ; 478 479 //site 480 Eigen::Map<const Vector> v(sites_ + Dim*site) ; 481 482 //site vectors 483 Vector sv[3] ; 484 for(int i = 0; i < 3; ++i) { 485 sv[i] = v-*(c[i]) ; 486 } 487 488 //unweighted objective function value 489 double local_fx = 0 ; 490 for(int i=0; i<3; ++i) { 491 for(int j=0; j<=i; ++j) { 492 local_fx += sv[i].dot(sv[j]) ; 493 } 494 } 495 local_fx /= 6. ; 496 497 //normal component of the site vector 498 const Vector nv = *(c[0]) - v - v.dot(ed[0])*ed[0] - v.dot(h[0])*h[0] ; 499 500 if(anisotropy_ > 1) { 501 local_fx += (anisotropy_*anisotropy_ - 1)*nv.dot(nv) ; 502 } 503 504 //weighted should add local_fx *= area which is done later to use this 505 //intermediate result for gradient computations. 506 507 //gradient of the objective function 508 Vector vertex_gradient ; 509 const Vector g = *(c[0]) + *(c[1]) + *(c[2]) - 4*v ; 510 for(int i=0; i<3; ++i) { 511 //gradient wrt. the vertex c[i] of the triangle 512 513 //variations of the integrand wrt. c[i] 514 vertex_gradient = *(c[i]) + g ; 515 516 if(anisotropy_ > 1) { 517 vertex_gradient -= 518 12*(anisotropy_*anisotropy_-1)/h_len[i]*(*(c[(i+1)%3])-v).dot(h[i])*nv ; 519 } 520 521 vertex_gradient *= area/6. ; 522 523 //variation of the triangle area wrt. c[i] 524 vertex_gradient += 0.5*ed_len[i]*h[i]*local_fx ; 525 526 //composing to get the gradient wrt. the mesh vertices 527 //number of mesh vertices involved in this RVD vertex 528 if(c[i]->config() == Vertex::ORIGINAL) { 529 //original vertex of the mesh 530 //grab the portion of the gradient wrt. the mesh vertex 531 Eigen::Map<Vector> vert_g(g_ + Dim*c[i]->combinatorics()[0]) ; 532 vert_g += factor_ * vertex_gradient / mesh_area_ ; 533 } else { 534 //get the number of mesh vertices involved 535 int comb_size = (c[i]->config() == Vertex::FACE_VERTEX) ? 3 : 2 ; 536 537 //for each vertex in the combinatorics, get the derivative of the RVD 538 //Vertex wrt. the mesh vertex, and compose the derivatives. 539 Matrix diff_compose ; 540 for(int j=0; j<comb_size; ++j) { 541 //get the derivative 542 c[i]->derivative(j, mesh_, sites_, diff_compose.data()) ; 543 //grab the portion of the gradient wrt. the mesh vertex 544 Eigen::Map<Vector> vert_g(g_ + Dim*c[i]->combinatorics()[j]) ; 545 //compose the gradient 546 vert_g += factor_*diff_compose.transpose()*vertex_gradient / mesh_area_ ; 547 } 548 } 549 } 550 551 //finishing the objective function value computation 552 fx_ += factor_ * area * local_fx / mesh_area_ ; 553 } 554 } ; 555 556 } //end of namespace CVT 557 } //end of namespace Revoropt 558 559 #endif 560