VCG Library
quality.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 #ifndef __VCG_TRI_UPDATE_QUALITY
24 #define __VCG_TRI_UPDATE_QUALITY
25 #include <vcg/complex/algorithms/stat.h>
26 
27 namespace vcg {
28 namespace tri {
30 
32 
34 
43 template <class UpdateMeshType>
45 {
46 public:
47  typedef UpdateMeshType MeshType;
48  typedef typename MeshType::ScalarType ScalarType;
49  typedef typename MeshType::CoordType CoordType;
50  typedef typename MeshType::VertexType VertexType;
51  typedef typename MeshType::VertexPointer VertexPointer;
52  typedef typename MeshType::VertexIterator VertexIterator;
53  typedef typename MeshType::FaceType FaceType;
54  typedef typename MeshType::FacePointer FacePointer;
55  typedef typename MeshType::FaceIterator FaceIterator;
56  typedef typename MeshType::VertexType::QualityType VertexQualityType;
57  typedef typename MeshType::FaceType::QualityType FaceQualityType;
58  typedef typename MeshType::TetraType TetraType;
59  typedef typename MeshType::TetraPointer TetraPointer;
60  typedef typename MeshType::TetraIterator TetraIterator;
61  typedef typename MeshType::TetraType::QualityType TetraQualityType;
62 
63 
64 
67 static void VertexConstant(MeshType &m, VertexQualityType q)
68 {
69  tri::RequirePerVertexQuality(m);
70  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
71  (*vi).Q()=q;
72 }
73 
76 static void VertexValence(UpdateMeshType &m)
77 {
78  tri::RequirePerVertexQuality(m);
79  VertexConstant(m,0);
80  for (size_t i=0;i<m.face.size();i++)
81  {
82  if (m.face[i].IsD())continue;
83 
84  for (int j=0;j<m.face[i].VN();j++)
85  {
86  VertexType *v=m.face[i].V(j);
87  v->Q()+=1;
88  }
89  }
90 }
91 
94 static void VertexClamp(MeshType &m,
95  VertexQualityType qmin,
96  VertexQualityType qmax)
97 {
98  tri::RequirePerVertexQuality(m);
99  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
100  (*vi).Q()=std::min(qmax, std::max(qmin,(*vi).Q()));
101 }
102 
105 static void VertexNormalize(MeshType &m, VertexQualityType qmin=0.0, VertexQualityType qmax=1.0)
106 {
107  tri::RequirePerVertexQuality(m);
108  ScalarType deltaRange = qmax-qmin;
109  std::pair<ScalarType,ScalarType> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax(m);
110  VertexIterator vi;
111  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
112  (*vi).Q() = qmin+deltaRange*((*vi).Q() - minmax.first)/(minmax.second - minmax.first);
113 }
114 
117 static void FaceNormalize(MeshType &m, FaceQualityType qmin=0.0, FaceQualityType qmax=1.0)
118 {
119  tri::RequirePerFaceQuality(m);
120  FaceQualityType deltaRange = qmax-qmin;
121  std::pair<FaceQualityType,FaceQualityType> minmax = tri::Stat<MeshType>::ComputePerFaceQualityMinMax(m);
122  for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
123  (*fi).Q() = qmin+deltaRange*((*fi).Q() - minmax.first)/(minmax.second - minmax.first);
124 }
125 
128 static void FaceConstant(MeshType &m, FaceQualityType q)
129 {
130  tri::RequirePerFaceQuality(m);
131  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
132  (*fi).Q()=q;
133 }
134 
137 static void FaceArea(MeshType &m)
138 {
139  tri::RequirePerFaceQuality(m);
140  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
141  (*fi).Q()=FaceQualityType(vcg::DoubleArea(*fi)/ScalarType(2.0));
142 }
143 
144 static void TetraConstant(MeshType & m, const TetraQualityType q)
145 {
146  tri::RequirePerTetraQuality(m);
147  ForEachTetra(m, [&q] (TetraType & t) {
148  t.Q() = q;
149  });
150 }
151 static void TetraFromVolume(MeshType & m)
152 {
153  tri::RequirePerTetraQuality(m);
154  ForEachTetra(m, [] (TetraType & t) {
155  t.Q() = TetraQualityType(vcg::Tetra::ComputeVolume(t));
156  });
157 }
158 
159 static void TetraFromAspectRatio(MeshType & m)
160 {
161  tri::RequirePerTetraQuality(m);
162  ForEachTetra(m, [] (TetraType & t) {
163  t.Q() = TetraQualityType(vcg::Tetra::AspectRatio(t));
164  });
165 }
166 
167 static void VertexFromFace( MeshType &m, bool areaWeighted=true)
168 {
169  tri::RequirePerFaceQuality(m);
170  tri::RequirePerVertexQuality(m);
171  SimpleTempData<typename MeshType::VertContainer, ScalarType> TQ(m.vert,0);
172  SimpleTempData<typename MeshType::VertContainer, ScalarType> TCnt(m.vert,0);
173 
174  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
175  if(!(*fi).IsD())
176  {
177  VertexQualityType weight=1.0;
178  if(areaWeighted) weight = vcg::DoubleArea(*fi);
179  for(int j=0;j<3;++j)
180  {
181  TQ[(*fi).V(j)]+=(*fi).Q()*weight;
182  TCnt[(*fi).V(j)]+=weight;
183  }
184  }
185 
186  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
187  if(!(*vi).IsD() && TCnt[*vi]>0 )
188  {
189  (*vi).Q() = TQ[*vi] / TCnt[*vi];
190  }
191 }
192 
193 static void VertexFromTetra(MeshType & m, bool volumeWeighted = true)
194 {
195  tri::RequirePerTetraQuality(m);
196  tri::RequirePerVertexQuality(m);
197 
198  SimpleTempData<typename MeshType::VertContainer, ScalarType> TQ(m.vert, 0);
199  SimpleTempData<typename MeshType::VertContainer, ScalarType> TCnt(m.vert, 0);
200 
201  ForEachTetra(m, [&] (TetraType & t) {
202  TetraQualityType w = 1.;
203  if (volumeWeighted)
204  w = vcg::Tetra::ComputeVolume(t);
205 
206  for (int i = 0; i < 4; ++i)
207  {
208  TQ[t.V(i)] += t.Q() * w;
209  TCnt[t.V(i)] += w;
210  }
211  });
212 
213  ForEachVertex(m, [&] (VertexType & v) {
214  v.Q() = TQ[v] / TCnt[v];
215  });
216 }
217 
218 template <class HandleScalar>
219 static void VertexFromAttributeHandle(MeshType &m, typename MeshType::template PerVertexAttributeHandle<HandleScalar> &h)
220 {
221  tri::RequirePerVertexQuality(m);
222  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
223  if(!(*vi).IsD())
224  (*vi).Q()=VertexQualityType(h[vi]);
225 }
226 
227 template <class HandleScalar>
228 static void FaceFromAttributeHandle(MeshType &m, typename MeshType::template PerFaceAttributeHandle<HandleScalar> &h)
229 {
230  tri::RequirePerFaceQuality(m);
231  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
232  (*fi).Q() =FaceQualityType(h[fi]);
233 }
234 
235 static void FaceFromVertex( MeshType &m)
236 {
237  tri::RequirePerFaceQuality(m);
238  tri::RequirePerVertexQuality(m);
239  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
240  {
241  (*fi).Q() =0;
242  for (int i=0;i<(*fi).VN();i++)
243  (*fi).Q() += (*fi).V(i)->Q();
244  (*fi).Q()/=(FaceQualityType)(*fi).VN();
245  }
246 }
247 
248 static void VertexFromPlane(MeshType &m, const Plane3<ScalarType> &pl)
249 {
250  tri::RequirePerVertexQuality(m);
251  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
252  (*vi).Q() =SignedDistancePlanePoint(pl,(*vi).cP());
253 }
254 
255 static void VertexFromGaussianCurvatureHG(MeshType &m)
256 {
257  tri::RequirePerVertexQuality(m);
258  tri::RequirePerVertexCurvature(m);
259  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
260  (*vi).Q() = (*vi).Kg();
261 }
262 
263 static void VertexFromMeanCurvatureHG(MeshType &m)
264 {
265  tri::RequirePerVertexQuality(m);
266  tri::RequirePerVertexCurvature(m);
267  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
268  (*vi).Q() = (*vi).Kh();
269 }
270 
271 static void VertexFromGaussianCurvatureDir(MeshType &m)
272 {
273  tri::RequirePerVertexQuality(m);
274  tri::RequirePerVertexCurvatureDir(m);
275  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
276  (*vi).Q() = (*vi).K1()*(*vi).K2();
277 }
278 
279 static void VertexFromMeanCurvatureDir(MeshType &m)
280 {
281  tri::RequirePerVertexQuality(m);
282  tri::RequirePerVertexCurvatureDir(m);
283  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
284  (*vi).Q() = ((*vi).K1()+(*vi).K2())/2.0f;
285 }
286 static void VertexFromMinCurvatureDir(MeshType &m)
287 {
288  tri::RequirePerVertexQuality(m);
289  tri::RequirePerVertexCurvatureDir(m);
290  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
291  (*vi).Q() = (*vi).K1();
292 }
293 static void VertexFromMaxCurvatureDir(MeshType &m)
294 {
295  tri::RequirePerVertexQuality(m);
296  tri::RequirePerVertexCurvatureDir(m);
297  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
298  (*vi).Q() = (*vi).K2();
299 }
300 
312 static void VertexFromShapeIndexCurvatureDir(MeshType &m)
313 {
314  tri::RequirePerVertexQuality(m);
315  tri::RequirePerVertexCurvatureDir(m);
316  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
317  {
318  ScalarType k1=(*vi).K1();
319  ScalarType k2=(*vi).K2();
320  if(k1<k2) std::swap(k1,k2);
321  (*vi).Q() = (2.0/M_PI)*atan2(k1+k2,k1-k2);
322  }
323 }
334 static void VertexFromCurvednessCurvatureDir(MeshType &m)
335 {
336  tri::RequirePerVertexQuality(m);
337  tri::RequirePerVertexCurvatureDir(m);
338  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
339  {
340  const ScalarType k1=(*vi).K1();
341  const ScalarType k2=(*vi).K2();
342 
343  (*vi).Q() = math::Sqrt((k1*k1+k2*k2)/2.0);
344  }
345 }
346 
347 
348 /*
349  * Absolute Curvature
350  *
351  * 2|H| if K >= 0
352  * |k1| + |k2| = <
353  * 2 * sqrt(|H|^2-K) otherwise
354  *
355  * defs and formulas taken from
356  *
357  * Improved curvature estimation for watershed segmentation of 3-dimensional meshes
358  * S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
359  * and from
360  * Optimizing 3D triangulations using discrete curvature analysis
361  * N Dyn, K Hormann, SJ Kim, D Levin - Mathematical Methods for Curves and Surfaces: Oslo, 2000
362  */
363 
364 static void VertexFromAbsoluteCurvature(MeshType &m)
365 {
366  tri::RequirePerVertexQuality(m);
367  tri::RequirePerVertexCurvature(m);
368  VertexIterator vi;
369  for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
370  {
371  if((*vi).Kg() >= 0)
372  (*vi).Q() = math::Abs( 2*(*vi).Kh() );
373  else
374  (*vi).Q() = 2*math::Sqrt(math::Abs( (*vi).Kh()*(*vi).Kh() - (*vi).Kg()));
375  }
376 }
377 
378 /*
379  * RMS Curvature = sqrt(4H^2-2K)
380  * def and formula taken from
381  *
382  * Improved curvature estimation for watershed segmentation of 3-dimensional meshes
383  * S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
384  */
385 static void VertexFromRMSCurvature(MeshType &m)
386 {
387  tri::RequirePerVertexQuality(m);
388  tri::RequirePerVertexCurvature(m);
389  VertexIterator vi;
390  for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
391  (*vi).Q() = math::Sqrt(math::Abs( 4*(*vi).Kh()*(*vi).Kh() - 2*(*vi).Kg()));
392 }
393 
394 /*
395  Saturate the Face quality so that for each vertex the gradient of the quality is lower than the given threshold value (in absolute value)
396  The saturation is done in a conservative way (quality is always decreased and never increased)
397 
398  Note: requires FF adjacency.
399  */
400 static void FaceSaturate(MeshType &m, FaceQualityType gradientThr=1.0)
401 {
402  tri::RequirePerFaceQuality(m);
403  tri::RequireFFAdjacency(m);
404  UpdateFlags<MeshType>::FaceClearV(m);
405  std::stack<FacePointer> st;
406 
407  st.push(&*m.face.begin());
408 
409  while(!st.empty())
410  {
411  FacePointer fc = st.top(); // the center
412  //printf("Stack size %i\n",st.size());
413  //printf("Pop elem %i %f\n",st.top() - &*m.vert.begin(), st.top()->Q());
414  st.pop();
415  fc->SetV();
416  std::vector<FacePointer> star;
417  typename std::vector<FacePointer>::iterator ffi;
418  for (int i=0;i<3;i++)
419  {
420  FacePointer fnext=fc->FFp(i);
421  if (fnext!=fc)star.push_back(fnext);
422  }
423  CoordType bary0=(fc->P(0)+fc->P(1)+fc->P(2))/3;
424  for(ffi=star.begin();ffi!=star.end();++ffi )
425  {
426  assert(fc!=(*ffi));
427  FaceQualityType &qi = (*ffi)->Q();
428  CoordType bary1=((*ffi)->P(0)+(*ffi)->P(1)+(*ffi)->P(2))/3;
429  FaceQualityType distGeom = Distance(bary0,bary1) / gradientThr;
430  // Main test if the quality varies more than the geometric displacement we have to lower something.
431  if( distGeom < fabs(qi - fc->Q()))
432  {
433  // center = 0 other=10 -> other =
434  // center = 10 other=0
435  if(fc->Q() > qi) // first case: the center of the star has to be lowered (and re-inserted in the queue).
436  {
437  //printf("Reinserting center %i \n",vc - &*m.vert.begin());
438  fc->Q() = qi+distGeom-(ScalarType)0.00001;
439  assert( distGeom > fabs(qi - fc->Q()));
440  st.push(fc);
441  break;
442  }
443  else
444  {
445  // second case: you have to lower qi, the vertex under examination.
446  assert( distGeom < fabs(qi - fc->Q()));
447  assert(fc->Q() < qi);
448  FaceQualityType newQi = fc->Q() + distGeom -(FaceQualityType)0.00001;
449  assert(newQi <= qi);
450  assert(fc->Q() < newQi);
451  assert( distGeom > fabs(newQi - fc->Q()) );
452 // printf("distGeom %f, qi %f, vc->Q() %f, fabs(qi - vc->Q()) %f\n",distGeom,qi,vc->Q(),fabs(qi - vc->Q()));
453  qi = newQi;
454  (*ffi)->ClearV();
455  }
456  }
457  if(!(*ffi)->IsV())
458  {
459  st.push( *ffi);
460 // printf("Reinserting side %i \n",*vvi - &*m.vert.begin());
461  (*ffi)->SetV();
462  }
463  }
464  }
465  }
466 
474 static void VertexSaturate(MeshType &m, ScalarType gradientThr=1.0)
475 {
476  tri::RequirePerVertexQuality(m);
477  tri::RequireVFAdjacency(m);
478 
480  std::stack<VertexPointer> st;
481 
482  st.push(&*m.vert.begin());
483 
484  while(!st.empty())
485  {
486  VertexPointer vc = st.top(); // the center
487  //printf("Stack size %i\n",st.size());
488  //printf("Pop elem %i %f\n",st.top() - &*m.vert.begin(), st.top()->Q());
489  st.pop();
490  vc->SetV();
491  std::vector<VertexPointer> star;
492  face::VVStarVF<FaceType>(vc,star);
493  for(auto vvi=star.begin();vvi!=star.end();++vvi )
494  {
495  ScalarType qi = (*vvi)->Q();
496  ScalarType distGeom = Distance((*vvi)->cP(),vc->cP()) / gradientThr;
497  // Main test if the quality varies more than the geometric displacement we have to lower something.
498  if( distGeom < fabs(qi - vc->Q()))
499  {
500  // center = 0 other=10 -> other =
501  // center = 10 other=0
502  if(vc->Q() > qi) // first case: the center of the star has to be lowered (and re-inserted in the queue).
503  {
504  //printf("Reinserting center %i \n",vc - &*m.vert.begin());
505  VertexQualityType delta=std::min(ScalarType(0.00001),ScalarType(distGeom/2.0));
506  vc->Q() = VertexQualityType(qi+distGeom-delta);
507  assert( distGeom > fabs(qi - vc->Q()));
508  st.push(vc);
509  break;
510  }
511  else
512  {
513  // second case: you have to lower qi, the vertex under examination.
514  assert( distGeom < fabs(qi - vc->Q()));
515  assert(vc->Q() < qi);
516  VertexQualityType delta=std::min(VertexQualityType(0.00001),VertexQualityType(distGeom/2.0));
517  VertexQualityType newQi = vc->Q() + distGeom -delta;
518  assert(newQi <= qi);
519  assert(vc->Q() < newQi);
520  assert( distGeom > fabs(newQi - vc->Q()) );
521 // printf("distGeom %f, qi %f, vc->Q() %f, fabs(qi - vc->Q()) %f\n",distGeom,qi,vc->Q(),fabs(qi - vc->Q()));
522  qi = newQi;
523  (*vvi)->ClearV();
524  }
525  }
526  if(!(*vvi)->IsV())
527  {
528  st.push( *vvi);
529 // printf("Reinserting side %i \n",*vvi - &*m.vert.begin());
530  (*vvi)->SetV();
531  }
532  }
533  }
534  }
535 
536 
537 }; //end class
538 } // end namespace
539 } // end namespace
540 #endif