$darkmode
VCG Library
fitmaps.h
1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 
24 #ifndef FITMAPS_H
25 #define FITMAPS_H
26 
27 #include <vcg/math/histogram.h>
28 
29 #include <vcg/simplex/face/jumping_pos.h>
30 #include <vcg/complex/algorithms/update/normal.h>
31 #include <vcg/complex/algorithms/update/curvature.h>
32 #include <vcg/complex/algorithms/update/topology.h>
33 #include <vcg/complex/algorithms/update/bounding.h>
34 #include "vcg/complex/algorithms/update/curvature_fitting.h"
35 
36 #include <Eigen/Core>
37 #include <Eigen/QR>
38 #include <Eigen/LU>
39 #include <Eigen/SVD>
40 
41 #include <vcg/complex/algorithms/nring.h>
42 
43 #include <vcg/complex/algorithms/smooth.h>
44 
45 namespace vcg { namespace tri {
46 
47 template<class MeshType>
48 class Fitmaps
49 {
50 public:
51 
52  typedef typename MeshType::FaceType FaceType;
53  typedef typename MeshType::VertexType VertexType;
54  typedef typename MeshType::ScalarType ScalarType;
55  typedef typename MeshType::FaceIterator FaceIterator;
56  typedef typename MeshType::VertexIterator VertexIterator;
57  typedef typename MeshType::CoordType CoordType;
58 
59  typedef vcg::tri::Nring<MeshType> RingWalker;
60 
61 
62  class Bicubic
63  {
64  public:
65 
66  Bicubic() {};
67 
68  Bicubic(vector<double>& input)
69  {
70  data = input;
71  if (input.size() != 16)
72  {
73  assert(0);
74  }
75  }
76 
77 // (u3 u2 u 1) (v3 v2 v 1)
78 //
79 // a u3 v3
80 // b u3 v2
81 // c u3 v1
82 // d u3 1
83 // e u2 v3
84 // f u2 v2
85 // g u2 v1
86 // h u2 1
87 // i u1 v3
88 // l u1 v2
89 // m u1 v1
90 // n u1 1
91 // o 1 v3
92 // p 1 v2
93 // q 1 v1
94 // r 1 1
95 
96  double& a() { return data[0];}
97  double& b() { return data[1];}
98  double& c() { return data[2];}
99  double& d() { return data[3];}
100  double& e() { return data[4];}
101  double& f() { return data[5];}
102  double& g() { return data[6];}
103  double& h() { return data[7];}
104  double& i() { return data[8];}
105  double& l() { return data[9];}
106  double& m() { return data[10];}
107  double& n() { return data[11];}
108  double& o() { return data[12];}
109  double& p() { return data[13];}
110  double& q() { return data[14];}
111  double& r() { return data[15];}
112 
113  vector<double> data;
114 
115  double evaluate(double u, double v)
116  {
117 
118  return
119  a() * u*u*u * v*v*v +
120  b() * u*u*u * v*v +
121  c() * u*u*u * v +
122  d() * u*u*u +
123  e() * u*u * v*v*v +
124  f() * u*u * v*v +
125  g() * u*u * v +
126  h() * u*u +
127  i() * u * v*v*v +
128  l() * u * v*v +
129  m() * u * v +
130  n() * u +
131  o() * v*v*v +
132  p() * v*v +
133  q() * v +
134  r();
135  }
136 
137 
138  double distanceRMS(std::vector<CoordType>& VV)
139  {
140  double error = 0;
141  for(typename std::vector<CoordType>::iterator it = VV.begin(); it != VV.end(); ++it)
142  {
143  double u = it->X();
144  double v = it->Y();
145  double n = it->Z();
146 
147  double temp = evaluate(u,v);
148  error += (n - temp)*(n - temp);
149  }
150  error /= (double) VV.size();
151  return sqrt(error);
152  }
153 
154  static Bicubic fit(std::vector<CoordType> VV)
155  {
156  assert(VV.size() >= 16);
157  Eigen::MatrixXf A(VV.size(),16);
158  Eigen::MatrixXf b(VV.size(),1);
159  Eigen::MatrixXf sol(16,1);
160 
161  for(unsigned int c=0; c < VV.size(); ++c)
162  {
163  double u = VV[c].X();
164  double v = VV[c].Y();
165  double n = VV[c].Z();
166 
167  A(c,0) = u*u*u * v*v*v;
168  A(c,1) = u*u*u * v*v;
169  A(c,2) = u*u*u * v;
170  A(c,3) = u*u*u;
171  A(c,4) = u*u * v*v*v;
172  A(c,5) = u*u * v*v;
173  A(c,6) = u*u * v;
174  A(c,7) = u*u;
175  A(c,8) = u * v*v*v;
176  A(c,9) = u * v*v;
177  A(c,10) = u * v;
178  A(c,11) = u;
179  A(c,12) = v*v*v;
180  A(c,13) = v*v;
181  A(c,14) = v;
182  A(c,15) = 1;
183 
184 
185  b[c] = n;
186  }
187 
188  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A);
189  sol=svd.solve(b);
190 // A.svd().solve(b, &sol);
191 
192  vector<double> r(16);
193 
194  for (int i=0; i < 16; ++i)
195  r.at(i) = sol[i];
196 
197  return Bicubic(r);
198  }
199  };
200 
201  Fitmaps()
202  {}
203 
204  class radSorter
205  {
206  public:
207 
208  radSorter(VertexType* v)
209  {
210  origin = v;
211  }
212 
213  VertexType* origin;
214 
215  bool operator() (VertexType* v1, VertexType* v2)
216  {
217  return (v1->P() - origin->P()).SquaredNorm() < (v2->P() - origin->P()).SquaredNorm();
218  }
219  };
220 
221  float getMeanCurvature(VertexType* vp)
222  {
223  return (vp->K1() + vp->K2())/2.0;
224  }
225 
226  static bool fitBicubicPoints(VertexType* v, std::vector<CoordType>& ref, Bicubic& ret, std::vector<CoordType>& points, std::vector<VertexType*>& ring)
227  {
228  points.clear();
229 
230  if (ring.size() < 16)
231  {
232  return false;
233  }
234 
235  typename std::vector<VertexType*>::iterator b = ring.begin();
236  typename std::vector<VertexType*>::iterator e = ring.end();
237 
238  while(b != e)
239  {
240  CoordType vT = (*b)->P() - v->P();
241 
242  double x = vT * ref[0];
243  double y = vT * ref[1];
244  double z = vT * ref[2];
245 
246  points.push_back(CoordType(x,y,z));
247  ++b;
248  }
249  ret = Bicubic::fit(points);
250  return true;
251  }
252 
253  static double AverageEdgeLenght(MeshType& m)
254  {
255  double doubleA = 0;
256  for (FaceIterator fi = m.face.begin(); fi!=m.face.end(); fi++) if (!fi->IsD()) {
257  doubleA+=vcg::DoubleArea(*fi);
258  }
259  int nquads = m.fn / 2;
260  return sqrt( doubleA/(2*nquads) );
261  }
262 
263  static void computeMFitmap(MeshType& m, float perc, int ringMax = 50)
264  {
267 
270 
272 
273  int countTemp = 0;
274 
275  RingWalker::clearFlags(&m);
276 
277  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
278  {
279  if ((countTemp++ % 100) == 0)
280  cerr << countTemp << "/" << m.vert.size() << endl;
281 
282  RingWalker rw(&*it,&m);
283 
284  CoordType nor = it->N();
285 
286  float okfaces = 0;
287  float flipfaces = 0;
288 
289  int count = 0;
290  do
291  {
292  count++;
293  rw.expand();
294  for(unsigned i=0; i<rw.lastF.size();++i)
295  {
296  CoordType vet1 = nor;
297  CoordType vet2 = rw.lastF[i]->N();
298 
299  vet1.Normalize();
300  vet2.Normalize();
301 
302 
303  double scal = vet1 * vet2;
304  if ((scal) > 0)
305  okfaces += (vcg::DoubleArea(*rw.lastF[i]));
306  else
307  flipfaces += (vcg::DoubleArea(*rw.lastF[i]));
308  }
309  } while ((((flipfaces)/(okfaces + flipfaces))*100.0 < perc) && (count < ringMax));
310 
311  std::sort(rw.lastV.begin(),rw.lastV.end(),radSorter(&*it));
312 
313  it->Q() = ((*rw.lastV.begin())->P() - it->P()).Norm();
314  rw.clear();
315 
316  }
317 
318  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,2);
319  }
320 
321  static vector<VertexType*> gatherNeighsSurface(VertexType* vt, float sigma, MeshType& m)
322  {
323  vector<VertexType*> current;
324 
325  RingWalker rw(vt,&m);
326 
327  bool exit = false;
328 
329  do
330  {
331  rw.expand();
332 
333  exit = true;
334 
335  for(typename vector<VertexType*>::iterator it = rw.lastV.begin(); it != rw.lastV.end(); ++it)
336  {
337  if (((*it)->P() - vt->P()).Norm() < sigma)
338  {
339  current.push_back(*it);
340  exit = false;
341  }
342  }
343 
344  } while (!exit);
345 
346  rw.clear();
347  return current;
348  }
349 
350  static void computeSFitmap(MeshType& m)//, float target = 1000)
351  {
352 
355 
358 
359  // update bounding box
361 
362  int countTemp = 0;
363 
364  double e = AverageEdgeLenght(m);
365 
366  int iteraz = 5; //2.0 * sqrt(m.vert.size()/target);
367 
368  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
369  {
370  if ((countTemp++ % 100) == 0)
371  cerr << countTemp << "/" << m.vert.size() << endl;
372 
373  vector<float> oneX;
374 
375 
376  for (int iteration = 0; iteration<iteraz; ++iteration)
377  {
378  oneX.push_back((iteration+1)*(e));
379  }
380 
381  std::vector<CoordType> ref(3);
382  ref[0] = it->PD1();
383  ref[1] = it->PD2();
384  ref[2] = it->PD1() ^ it->PD2();
385 
386  ref[0].Normalize();
387  ref[1].Normalize();
388  ref[2].Normalize();
389 
390  Bicubic b;
391 
392  RingWalker::clearFlags(&m);
393 
394  std::vector<VertexType*> pointsGlobal = gatherNeighsSurface(&*it,oneX.at(iteraz-1),m);
395 
396  vector<float> onedimensional;
397 
398  for (int iteration = 0; iteration<iteraz; ++iteration)
399  {
400  std::vector<VertexType*> points; // solo quelli nel raggio
401 
402  std::vector<CoordType> projected; // come sopra ma in coord locali
403  for (typename std::vector<VertexType*>::iterator it2 = pointsGlobal.begin(); it2 != pointsGlobal.end(); ++it2)
404  {
405  if (((*it).P() - (*it2)->P()).Norm() < oneX.at(iteration))
406  points.push_back(*it2);
407  }
408 
409  std::vector<VertexType*>& pointsFitting = points;
410 
411 
412  if (!fitBicubicPoints(&*it, ref, b, projected,pointsFitting))
413  {
414  onedimensional.push_back(0);
415  }
416  else
417  {
418  onedimensional.push_back(b.distanceRMS(projected));
419  }
420 
421  }
422 
423 
424  // // vecchio fit ax^4
425  Eigen::MatrixXf Am(onedimensional.size(),1);
426  Eigen::MatrixXf bm(onedimensional.size(),1);
427  Eigen::MatrixXf sol(1,1);
428 
429  for(unsigned int c=0; c < onedimensional.size(); ++c)
430  {
431  double x = oneX.at(c);
432 
433  Am(c,0) = pow(x,4);
434  bm[c] = onedimensional[c];
435  }
436 
437 
438  // Am.svd().solve(bm, &sol);
439  Eigen::JacobiSVD<Eigen::MatrixXd> svd(Am);
440  sol=svd.solve(bm);
441 
442  it->Q() = pow((double)sol[0],0.25);
443 
444  // // nuovo fit ax^4 + b
445  // Eigen::MatrixXf Am(onedimensional.size()+1,2);
446  // Eigen::MatrixXf bm(onedimensional.size()+1,1);
447  // Eigen::MatrixXf sol(2,1);
448  //
449  // Am(0,0) = 0;
450  // Am(0,1) = 0;
451  // bm[0] = 0;
452  //
453  // for(unsigned int c=0; c < onedimensional.size(); ++c)
454  // {
455  // double x = oneX.at(c);
456  //
457  // Am(c,0) = pow(x,4);
458  // Am(c,1) = 1;
459  // bm[c] = onedimensional[c];
460  // }
461  //
462  // //sol = ((Am.transpose()*Am).inverse()*Am.transpose())*bm;
463  // Am.svd().solve(bm, &sol);
464  //
465  // cerr << "------" << sol[0] << " " << sol[1] << endl;
466  // if (sol[0] > 0)
467  // saliency[it] = pow((double)sol[0],0.25);
468  // else
469  // saliency[it] = 0;
470 
471 
472  }
473 
474  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,1);
475  }
476 
477 
478  ~Fitmaps(){};
479 
480 };
481 
482 }} // END NAMESPACES
483 
484 #endif // FITMAPS_H
Definition: fitmaps.h:63
Definition: fitmaps.h:205
Definition: fitmaps.h:49
static void Box(ComputeMeshType &m)
Calculates the bounding box of the given mesh m.
Definition: bounding.h:45
Computation of per-vertex directions and values of curvature.
Definition: curvature_fitting.h:62
static void PerVertexAngleWeighted(ComputeMeshType &m)
Calculates the vertex normal as an angle weighted average. It does not need or exploit current face n...
Definition: normal.h:124
static void VertexFace(MeshType &m)
Update the Vertex-Face topological relation.
Definition: topology.h:467
static void FaceFace(MeshType &m)
Update the Face-Face topological relation by allowing to retrieve for each face what other faces shar...
Definition: topology.h:395
Definition: namespaces.dox:6