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  A.svd().solve(b, &sol);
192 
193  vector<double> r(16);
194 
195  for (int i=0; i < 16; ++i)
196  r.at(i) = sol[i];
197 
198  return Bicubic(r);
199  }
200  };
201 
202  Fitmaps()
203  {}
204 
205  class radSorter
206  {
207  public:
208 
209  radSorter(VertexType* v)
210  {
211  origin = v;
212  }
213 
214  VertexType* origin;
215 
216  bool operator() (VertexType* v1, VertexType* v2)
217  {
218  return (v1->P() - origin->P()).SquaredNorm() < (v2->P() - origin->P()).SquaredNorm();
219  }
220  };
221 
222  float getMeanCurvature(VertexType* vp)
223  {
224  return (vp->K1() + vp->K2())/2.0;
225  }
226 
227  static bool fitBicubicPoints(VertexType* v, std::vector<CoordType>& ref, Bicubic& ret, std::vector<CoordType>& points, std::vector<VertexType*>& ring)
228  {
229  points.clear();
230 
231  if (ring.size() < 16)
232  {
233  return false;
234  }
235 
236  typename std::vector<VertexType*>::iterator b = ring.begin();
237  typename std::vector<VertexType*>::iterator e = ring.end();
238 
239  while(b != e)
240  {
241  CoordType vT = (*b)->P() - v->P();
242 
243  double x = vT * ref[0];
244  double y = vT * ref[1];
245  double z = vT * ref[2];
246 
247  points.push_back(CoordType(x,y,z));
248  ++b;
249  }
250  ret = Bicubic::fit(points);
251  return true;
252  }
253 
254  static double AverageEdgeLenght(MeshType& m)
255  {
256  double doubleA = 0;
257  for (FaceIterator fi = m.face.begin(); fi!=m.face.end(); fi++) if (!fi->IsD()) {
258  doubleA+=vcg::DoubleArea(*fi);
259  }
260  int nquads = m.fn / 2;
261  return sqrt( doubleA/(2*nquads) );
262  }
263 
264  static void computeMFitmap(MeshType& m, float perc, int ringMax = 50)
265  {
268 
271 
273 
274  int countTemp = 0;
275 
276  RingWalker::clearFlags(&m);
277 
278  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
279  {
280  if ((countTemp++ % 100) == 0)
281  cerr << countTemp << "/" << m.vert.size() << endl;
282 
283  RingWalker rw(&*it,&m);
284 
285  CoordType nor = it->N();
286 
287  float okfaces = 0;
288  float flipfaces = 0;
289 
290  int count = 0;
291  do
292  {
293  count++;
294  rw.expand();
295  for(unsigned i=0; i<rw.lastF.size();++i)
296  {
297  CoordType vet1 = nor;
298  CoordType vet2 = rw.lastF[i]->N();
299 
300  vet1.Normalize();
301  vet2.Normalize();
302 
303 
304  double scal = vet1 * vet2;
305  if ((scal) > 0)
306  okfaces += (vcg::DoubleArea(*rw.lastF[i]));
307  else
308  flipfaces += (vcg::DoubleArea(*rw.lastF[i]));
309  }
310  } while ((((flipfaces)/(okfaces + flipfaces))*100.0 < perc) && (count < ringMax));
311 
312  std::sort(rw.lastV.begin(),rw.lastV.end(),radSorter(&*it));
313 
314  it->Q() = ((*rw.lastV.begin())->P() - it->P()).Norm();
315  rw.clear();
316 
317  }
318 
319  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,2);
320  }
321 
322  static vector<VertexType*> gatherNeighsSurface(VertexType* vt, float sigma, MeshType& m)
323  {
324  vector<VertexType*> current;
325 
326  RingWalker rw(vt,&m);
327 
328  bool exit = false;
329 
330  do
331  {
332  rw.expand();
333 
334  exit = true;
335 
336  for(typename vector<VertexType*>::iterator it = rw.lastV.begin(); it != rw.lastV.end(); ++it)
337  {
338  if (((*it)->P() - vt->P()).Norm() < sigma)
339  {
340  current.push_back(*it);
341  exit = false;
342  }
343  }
344 
345  } while (!exit);
346 
347  rw.clear();
348  return current;
349  }
350 
351  static void computeSFitmap(MeshType& m)//, float target = 1000)
352  {
353 
356 
359 
360  // update bounding box
362 
363  int countTemp = 0;
364 
365  double e = AverageEdgeLenght(m);
366 
367  int iteraz = 5; //2.0 * sqrt(m.vert.size()/target);
368 
369  for(VertexIterator it=m.vert.begin(); it!=m.vert.end();++it)
370  {
371  if ((countTemp++ % 100) == 0)
372  cerr << countTemp << "/" << m.vert.size() << endl;
373 
374  vector<float> oneX;
375 
376 
377  for (int iteration = 0; iteration<iteraz; ++iteration)
378  {
379  oneX.push_back((iteration+1)*(e));
380  }
381 
382  std::vector<CoordType> ref(3);
383  ref[0] = it->PD1();
384  ref[1] = it->PD2();
385  ref[2] = it->PD1() ^ it->PD2();
386 
387  ref[0].Normalize();
388  ref[1].Normalize();
389  ref[2].Normalize();
390 
391  Bicubic b;
392 
393  RingWalker::clearFlags(&m);
394 
395  std::vector<VertexType*> pointsGlobal = gatherNeighsSurface(&*it,oneX.at(iteraz-1),m);
396 
397  vector<float> onedimensional;
398 
399  for (int iteration = 0; iteration<iteraz; ++iteration)
400  {
401  std::vector<VertexType*> points; // solo quelli nel raggio
402 
403  std::vector<CoordType> projected; // come sopra ma in coord locali
404  for (typename std::vector<VertexType*>::iterator it2 = pointsGlobal.begin(); it2 != pointsGlobal.end(); ++it2)
405  {
406  if (((*it).P() - (*it2)->P()).Norm() < oneX.at(iteration))
407  points.push_back(*it2);
408  }
409 
410  std::vector<VertexType*>& pointsFitting = points;
411 
412 
413  if (!fitBicubicPoints(&*it, ref, b, projected,pointsFitting))
414  {
415  onedimensional.push_back(0);
416  }
417  else
418  {
419  onedimensional.push_back(b.distanceRMS(projected));
420  }
421 
422  }
423 
424 
425  // // vecchio fit ax^4
426  Eigen::MatrixXf Am(onedimensional.size(),1);
427  Eigen::MatrixXf bm(onedimensional.size(),1);
428  Eigen::MatrixXf sol(1,1);
429 
430  for(unsigned int c=0; c < onedimensional.size(); ++c)
431  {
432  double x = oneX.at(c);
433 
434  Am(c,0) = pow(x,4);
435  bm[c] = onedimensional[c];
436  }
437 
438  Am.svd().solve(bm, &sol);
439 
440  it->Q() = pow((double)sol[0],0.25);
441 
442  // // nuovo fit ax^4 + b
443  // Eigen::MatrixXf Am(onedimensional.size()+1,2);
444  // Eigen::MatrixXf bm(onedimensional.size()+1,1);
445  // Eigen::MatrixXf sol(2,1);
446  //
447  // Am(0,0) = 0;
448  // Am(0,1) = 0;
449  // bm[0] = 0;
450  //
451  // for(unsigned int c=0; c < onedimensional.size(); ++c)
452  // {
453  // double x = oneX.at(c);
454  //
455  // Am(c,0) = pow(x,4);
456  // Am(c,1) = 1;
457  // bm[c] = onedimensional[c];
458  // }
459  //
460  // //sol = ((Am.transpose()*Am).inverse()*Am.transpose())*bm;
461  // Am.svd().solve(bm, &sol);
462  //
463  // cerr << "------" << sol[0] << " " << sol[1] << endl;
464  // if (sol[0] > 0)
465  // saliency[it] = pow((double)sol[0],0.25);
466  // else
467  // saliency[it] = 0;
468 
469 
470  }
471 
472  vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(m,1);
473  }
474 
475 
476  ~Fitmaps(){};
477 
478 };
479 
480 }} // END NAMESPACES
481 
482 #endif // FITMAPS_H