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 
46 using namespace Eigen;
47 
48 namespace vcg { namespace tri {
49 
50 template<class MeshType>
51 class Fitmaps
52 {
53 public:
54 
55  typedef typename MeshType::FaceType FaceType;
56  typedef typename MeshType::VertexType VertexType;
57  typedef typename MeshType::ScalarType ScalarType;
58  typedef typename MeshType::FaceIterator FaceIterator;
59  typedef typename MeshType::VertexIterator VertexIterator;
60  typedef typename MeshType::CoordType CoordType;
61 
62  typedef vcg::tri::Nring<MeshType> RingWalker;
63 
64 
65  class Bicubic
66  {
67  public:
68 
69  Bicubic() {};
70 
71  Bicubic(vector<double>& input)
72  {
73  data = input;
74  if (input.size() != 16)
75  {
76  assert(0);
77  }
78  }
79 
80 // (u3 u2 u 1) (v3 v2 v 1)
81 //
82 // a u3 v3
83 // b u3 v2
84 // c u3 v1
85 // d u3 1
86 // e u2 v3
87 // f u2 v2
88 // g u2 v1
89 // h u2 1
90 // i u1 v3
91 // l u1 v2
92 // m u1 v1
93 // n u1 1
94 // o 1 v3
95 // p 1 v2
96 // q 1 v1
97 // r 1 1
98 
99  double& a() { return data[0];}
100  double& b() { return data[1];}
101  double& c() { return data[2];}
102  double& d() { return data[3];}
103  double& e() { return data[4];}
104  double& f() { return data[5];}
105  double& g() { return data[6];}
106  double& h() { return data[7];}
107  double& i() { return data[8];}
108  double& l() { return data[9];}
109  double& m() { return data[10];}
110  double& n() { return data[11];}
111  double& o() { return data[12];}
112  double& p() { return data[13];}
113  double& q() { return data[14];}
114  double& r() { return data[15];}
115 
116  vector<double> data;
117 
118  double evaluate(double u, double v)
119  {
120 
121  return
122  a() * u*u*u * v*v*v +
123  b() * u*u*u * v*v +
124  c() * u*u*u * v +
125  d() * u*u*u +
126  e() * u*u * v*v*v +
127  f() * u*u * v*v +
128  g() * u*u * v +
129  h() * u*u +
130  i() * u * v*v*v +
131  l() * u * v*v +
132  m() * u * v +
133  n() * u +
134  o() * v*v*v +
135  p() * v*v +
136  q() * v +
137  r();
138  }
139 
140 
141  double distanceRMS(std::vector<CoordType>& VV)
142  {
143  double error = 0;
144  for(typename std::vector<CoordType>::iterator it = VV.begin(); it != VV.end(); ++it)
145  {
146  double u = it->X();
147  double v = it->Y();
148  double n = it->Z();
149 
150  double temp = evaluate(u,v);
151  error += (n - temp)*(n - temp);
152  }
153  error /= (double) VV.size();
154  return sqrt(error);
155  }
156 
157  static Bicubic fit(std::vector<CoordType> VV)
158  {
159  assert(VV.size() >= 16);
160  Eigen::MatrixXf A(VV.size(),16);
161  Eigen::MatrixXf b(VV.size(),1);
162  Eigen::MatrixXf sol(16,1);
163 
164  for(unsigned int c=0; c < VV.size(); ++c)
165  {
166  double u = VV[c].X();
167  double v = VV[c].Y();
168  double n = VV[c].Z();
169 
170  A(c,0) = u*u*u * v*v*v;
171  A(c,1) = u*u*u * v*v;
172  A(c,2) = u*u*u * v;
173  A(c,3) = u*u*u;
174  A(c,4) = u*u * v*v*v;
175  A(c,5) = u*u * v*v;
176  A(c,6) = u*u * v;
177  A(c,7) = u*u;
178  A(c,8) = u * v*v*v;
179  A(c,9) = u * v*v;
180  A(c,10) = u * v;
181  A(c,11) = u;
182  A(c,12) = v*v*v;
183  A(c,13) = v*v;
184  A(c,14) = v;
185  A(c,15) = 1;
186 
187 
188  b[c] = n;
189  }
190 
191  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A);
192  sol=svd.solve(b);
193 // A.svd().solve(b, &sol);
194 
195  vector<double> r(16);
196 
197  for (int i=0; i < 16; ++i)
198  r.at(i) = sol[i];
199 
200  return Bicubic(r);
201  }
202  };
203 
204  Fitmaps()
205  {}
206 
207  class radSorter
208  {
209  public:
210 
211  radSorter(VertexType* v)
212  {
213  origin = v;
214  }
215 
216  VertexType* origin;
217 
218  bool operator() (VertexType* v1, VertexType* v2)
219  {
220  return (v1->P() - origin->P()).SquaredNorm() < (v2->P() - origin->P()).SquaredNorm();
221  }
222  };
223 
224  float getMeanCurvature(VertexType* vp)
225  {
226  return (vp->K1() + vp->K2())/2.0;
227  }
228 
229  static bool fitBicubicPoints(VertexType* v, std::vector<CoordType>& ref, Bicubic& ret, std::vector<CoordType>& points, std::vector<VertexType*>& ring)
230  {
231  points.clear();
232 
233  if (ring.size() < 16)
234  {
235  return false;
236  }
237 
238  typename std::vector<VertexType*>::iterator b = ring.begin();
239  typename std::vector<VertexType*>::iterator e = ring.end();
240 
241  while(b != e)
242  {
243  CoordType vT = (*b)->P() - v->P();
244 
245  double x = vT * ref[0];
246  double y = vT * ref[1];
247  double z = vT * ref[2];
248 
249  points.push_back(CoordType(x,y,z));
250  ++b;
251  }
252  ret = Bicubic::fit(points);
253  return true;
254  }
255 
256  static double AverageEdgeLenght(MeshType& m)
257  {
258  double doubleA = 0;
259  for (FaceIterator fi = m.face.begin(); fi!=m.face.end(); fi++) if (!fi->IsD()) {
260  doubleA+=vcg::DoubleArea(*fi);
261  }
262  int nquads = m.fn / 2;
263  return sqrt( doubleA/(2*nquads) );
264  }
265 
266  static void computeMFitmap(MeshType& m, float perc, int ringMax = 50)
267  {
270 
273 
275 
276  int countTemp = 0;
277 
278  RingWalker::clearFlags(&m);
279 
280  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
281  {
282  if ((countTemp++ % 100) == 0)
283  cerr << countTemp << "/" << m.vert.size() << endl;
284 
285  RingWalker rw(&*it,&m);
286 
287  CoordType nor = it->N();
288 
289  float okfaces = 0;
290  float flipfaces = 0;
291 
292  int count = 0;
293  do
294  {
295  count++;
296  rw.expand();
297  for(unsigned i=0; i<rw.lastF.size();++i)
298  {
299  CoordType vet1 = nor;
300  CoordType vet2 = rw.lastF[i]->N();
301 
302  vet1.Normalize();
303  vet2.Normalize();
304 
305 
306  double scal = vet1 * vet2;
307  if ((scal) > 0)
308  okfaces += (vcg::DoubleArea(*rw.lastF[i]));
309  else
310  flipfaces += (vcg::DoubleArea(*rw.lastF[i]));
311  }
312  } while ((((flipfaces)/(okfaces + flipfaces))*100.0 < perc) && (count < ringMax));
313 
314  std::sort(rw.lastV.begin(),rw.lastV.end(),radSorter(&*it));
315 
316  it->Q() = ((*rw.lastV.begin())->P() - it->P()).Norm();
317  rw.clear();
318 
319  }
320 
321  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,2);
322  }
323 
324  static vector<VertexType*> gatherNeighsSurface(VertexType* vt, float sigma, MeshType& m)
325  {
326  vector<VertexType*> current;
327 
328  RingWalker rw(vt,&m);
329 
330  bool exit = false;
331 
332  do
333  {
334  rw.expand();
335 
336  exit = true;
337 
338  for(typename vector<VertexType*>::iterator it = rw.lastV.begin(); it != rw.lastV.end(); ++it)
339  {
340  if (((*it)->P() - vt->P()).Norm() < sigma)
341  {
342  current.push_back(*it);
343  exit = false;
344  }
345  }
346 
347  } while (!exit);
348 
349  rw.clear();
350  return current;
351  }
352 
353  static void computeSFitmap(MeshType& m)//, float target = 1000)
354  {
355 
358 
361 
362  // update bounding box
364 
365  int countTemp = 0;
366 
367  double e = AverageEdgeLenght(m);
368 
369  int iteraz = 5; //2.0 * sqrt(m.vert.size()/target);
370 
371  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
372  {
373  if ((countTemp++ % 100) == 0)
374  cerr << countTemp << "/" << m.vert.size() << endl;
375 
376  vector<float> oneX;
377 
378 
379  for (int iteration = 0; iteration<iteraz; ++iteration)
380  {
381  oneX.push_back((iteration+1)*(e));
382  }
383 
384  std::vector<CoordType> ref(3);
385  ref[0] = it->PD1();
386  ref[1] = it->PD2();
387  ref[2] = it->PD1() ^ it->PD2();
388 
389  ref[0].Normalize();
390  ref[1].Normalize();
391  ref[2].Normalize();
392 
393  Bicubic b;
394 
395  RingWalker::clearFlags(&m);
396 
397  std::vector<VertexType*> pointsGlobal = gatherNeighsSurface(&*it,oneX.at(iteraz-1),m);
398 
399  vector<float> onedimensional;
400 
401  for (int iteration = 0; iteration<iteraz; ++iteration)
402  {
403  std::vector<VertexType*> points; // solo quelli nel raggio
404 
405  std::vector<CoordType> projected; // come sopra ma in coord locali
406  for (typename std::vector<VertexType*>::iterator it2 = pointsGlobal.begin(); it2 != pointsGlobal.end(); ++it2)
407  {
408  if (((*it).P() - (*it2)->P()).Norm() < oneX.at(iteration))
409  points.push_back(*it2);
410  }
411 
412  std::vector<VertexType*>& pointsFitting = points;
413 
414 
415  if (!fitBicubicPoints(&*it, ref, b, projected,pointsFitting))
416  {
417  onedimensional.push_back(0);
418  }
419  else
420  {
421  onedimensional.push_back(b.distanceRMS(projected));
422  }
423 
424  }
425 
426 
427  // // vecchio fit ax^4
428  Eigen::MatrixXf Am(onedimensional.size(),1);
429  Eigen::MatrixXf bm(onedimensional.size(),1);
430  Eigen::MatrixXf sol(1,1);
431 
432  for(unsigned int c=0; c < onedimensional.size(); ++c)
433  {
434  double x = oneX.at(c);
435 
436  Am(c,0) = pow(x,4);
437  bm[c] = onedimensional[c];
438  }
439 
440 
441  // Am.svd().solve(bm, &sol);
442  Eigen::JacobiSVD<Eigen::MatrixXd> svd(Am);
443  sol=svd.solve(bm);
444 
445  it->Q() = pow((double)sol[0],0.25);
446 
447  // // nuovo fit ax^4 + b
448  // Eigen::MatrixXf Am(onedimensional.size()+1,2);
449  // Eigen::MatrixXf bm(onedimensional.size()+1,1);
450  // Eigen::MatrixXf sol(2,1);
451  //
452  // Am(0,0) = 0;
453  // Am(0,1) = 0;
454  // bm[0] = 0;
455  //
456  // for(unsigned int c=0; c < onedimensional.size(); ++c)
457  // {
458  // double x = oneX.at(c);
459  //
460  // Am(c,0) = pow(x,4);
461  // Am(c,1) = 1;
462  // bm[c] = onedimensional[c];
463  // }
464  //
465  // //sol = ((Am.transpose()*Am).inverse()*Am.transpose())*bm;
466  // Am.svd().solve(bm, &sol);
467  //
468  // cerr << "------" << sol[0] << " " << sol[1] << endl;
469  // if (sol[0] > 0)
470  // saliency[it] = pow((double)sol[0],0.25);
471  // else
472  // saliency[it] = 0;
473 
474 
475  }
476 
477  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,1);
478  }
479 
480 
481  ~Fitmaps(){};
482 
483 };
484 
485 }} // END NAMESPACES
486 
487 #endif // FITMAPS_H