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 
59 
62 static void VertexConstant(MeshType &m, VertexQualityType q)
63 {
64  tri::RequirePerVertexQuality(m);
65  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
66  (*vi).Q()=q;
67 }
68 
71 static void VertexValence(UpdateMeshType &m)
72 {
73  tri::RequirePerVertexQuality(m);
74  VertexConstant(m,0);
75  for (size_t i=0;i<m.face.size();i++)
76  {
77  if (m.face[i].IsD())continue;
78 
79  for (int j=0;j<m.face[i].VN();j++)
80  {
81  VertexType *v=m.face[i].V(j);
82  v->Q()+=1;
83  }
84  }
85 }
86 
89 static void VertexClamp(MeshType &m,
90  VertexQualityType qmin,
91  VertexQualityType qmax)
92 {
93  tri::RequirePerVertexQuality(m);
94  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
95  (*vi).Q()=std::min(qmax, std::max(qmin,(*vi).Q()));
96 }
97 
100 static void VertexNormalize(MeshType &m, VertexQualityType qmin=0.0, VertexQualityType qmax=1.0)
101 {
102  tri::RequirePerVertexQuality(m);
103  ScalarType deltaRange = qmax-qmin;
104  std::pair<ScalarType,ScalarType> minmax = tri::Stat<MeshType>::ComputePerVertexQualityMinMax(m);
105  VertexIterator vi;
106  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
107  (*vi).Q() = qmin+deltaRange*((*vi).Q() - minmax.first)/(minmax.second - minmax.first);
108 }
109 
112 static void FaceNormalize(MeshType &m, FaceQualityType qmin=0.0, FaceQualityType qmax=1.0)
113 {
114  tri::RequirePerFaceQuality(m);
115  FaceQualityType deltaRange = qmax-qmin;
116  std::pair<FaceQualityType,FaceQualityType> minmax = tri::Stat<MeshType>::ComputePerFaceQualityMinMax(m);
117  for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi)
118  (*fi).Q() = qmin+deltaRange*((*fi).Q() - minmax.first)/(minmax.second - minmax.first);
119 }
120 
123 static void FaceConstant(MeshType &m, FaceQualityType q)
124 {
125  tri::RequirePerFaceQuality(m);
126  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
127  (*fi).Q()=q;
128 }
129 
132 static void FaceArea(MeshType &m)
133 {
134  tri::RequirePerFaceQuality(m);
135  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
136  (*fi).Q()=FaceQualityType(vcg::DoubleArea(*fi)/ScalarType(2.0));
137 }
138 
139 static void VertexFromFace( MeshType &m, bool areaWeighted=true)
140 {
141  tri::RequirePerFaceQuality(m);
142  tri::RequirePerVertexQuality(m);
143  SimpleTempData<typename MeshType::VertContainer, ScalarType> TQ(m.vert,0);
144  SimpleTempData<typename MeshType::VertContainer, ScalarType> TCnt(m.vert,0);
145 
146  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
147  if(!(*fi).IsD())
148  {
149  VertexQualityType weight=1.0;
150  if(areaWeighted) weight = vcg::DoubleArea(*fi);
151  for(int j=0;j<3;++j)
152  {
153  TQ[(*fi).V(j)]+=(*fi).Q()*weight;
154  TCnt[(*fi).V(j)]+=weight;
155  }
156  }
157 
158  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
159  if(!(*vi).IsD() && TCnt[*vi]>0 )
160  {
161  (*vi).Q() = TQ[*vi] / TCnt[*vi];
162  }
163 }
164 
165 template <class HandleScalar>
166 static void VertexFromAttributeHandle(MeshType &m, typename MeshType::template PerVertexAttributeHandle<HandleScalar> &h)
167 {
168  tri::RequirePerVertexQuality(m);
169  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
170  if(!(*vi).IsD())
171  (*vi).Q()=VertexQualityType(h[vi]);
172 }
173 
174 template <class HandleScalar>
175 static void FaceFromAttributeHandle(MeshType &m, typename MeshType::template PerFaceAttributeHandle<HandleScalar> &h)
176 {
177  tri::RequirePerFaceQuality(m);
178  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
179  (*fi).Q() =FaceQualityType(h[fi]);
180 }
181 
182 static void FaceFromVertex( MeshType &m)
183 {
184  tri::RequirePerFaceQuality(m);
185  tri::RequirePerVertexQuality(m);
186  for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
187  {
188  (*fi).Q() =0;
189  for (size_t i=0;i<(*fi).VN();i++)
190  (*fi).Q() += (*fi).V(i)->Q();
191  (*fi).Q()/=(FaceQualityType)(*fi).VN();
192  }
193 }
194 
195 static void VertexFromPlane(MeshType &m, const Plane3<ScalarType> &pl)
196 {
197  tri::RequirePerVertexQuality(m);
198  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
199  (*vi).Q() =SignedDistancePlanePoint(pl,(*vi).cP());
200 }
201 
202 static void VertexFromGaussianCurvatureHG(MeshType &m)
203 {
204  tri::RequirePerVertexQuality(m);
205  tri::RequirePerVertexCurvature(m);
206  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
207  (*vi).Q() = (*vi).Kg();
208 }
209 
210 static void VertexFromMeanCurvatureHG(MeshType &m)
211 {
212  tri::RequirePerVertexQuality(m);
213  tri::RequirePerVertexCurvature(m);
214  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
215  (*vi).Q() = (*vi).Kh();
216 }
217 
218 static void VertexFromGaussianCurvatureDir(MeshType &m)
219 {
220  tri::RequirePerVertexQuality(m);
221  tri::RequirePerVertexCurvatureDir(m);
222  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
223  (*vi).Q() = (*vi).K1()*(*vi).K2();
224 }
225 
226 static void VertexFromMeanCurvatureDir(MeshType &m)
227 {
228  tri::RequirePerVertexQuality(m);
229  tri::RequirePerVertexCurvatureDir(m);
230  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
231  (*vi).Q() = ((*vi).K1()+(*vi).K2())/2.0f;
232 }
233 static void VertexFromMinCurvatureDir(MeshType &m)
234 {
235  tri::RequirePerVertexQuality(m);
236  tri::RequirePerVertexCurvatureDir(m);
237  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
238  (*vi).Q() = (*vi).K1();
239 }
240 static void VertexFromMaxCurvatureDir(MeshType &m)
241 {
242  tri::RequirePerVertexQuality(m);
243  tri::RequirePerVertexCurvatureDir(m);
244  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
245  (*vi).Q() = (*vi).K2();
246 }
247 
259 static void VertexFromShapeIndexCurvatureDir(MeshType &m)
260 {
261  tri::RequirePerVertexQuality(m);
262  tri::RequirePerVertexCurvatureDir(m);
263  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
264  {
265  ScalarType k1=(*vi).K1();
266  ScalarType k2=(*vi).K2();
267  if(k1<k2) std::swap(k1,k2);
268  (*vi).Q() = (2.0/M_PI)*atan2(k1+k2,k1-k2);
269  }
270 }
281 static void VertexFromCurvednessCurvatureDir(MeshType &m)
282 {
283  tri::RequirePerVertexQuality(m);
284  tri::RequirePerVertexCurvatureDir(m);
285  for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
286  {
287  const ScalarType k1=(*vi).K1();
288  const ScalarType k2=(*vi).K2();
289 
290  (*vi).Q() = math::Sqrt((k1*k1+k2*k2)/2.0);
291  }
292 }
293 
294 
295 /*
296  * Absolute Curvature
297  *
298  * 2|H| if K >= 0
299  * |k1| + |k2| = <
300  * 2 * sqrt(|H|^2-K) otherwise
301  *
302  * defs and formulas taken from
303  *
304  * Improved curvature estimation for watershed segmentation of 3-dimensional meshes
305  * S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
306  * and from
307  * Optimizing 3D triangulations using discrete curvature analysis
308  * N Dyn, K Hormann, SJ Kim, D Levin - Mathematical Methods for Curves and Surfaces: Oslo, 2000
309  */
310 
311 static void VertexFromAbsoluteCurvature(MeshType &m)
312 {
313  tri::RequirePerVertexQuality(m);
314  tri::RequirePerVertexCurvature(m);
315  VertexIterator vi;
316  for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
317  {
318  if((*vi).Kg() >= 0)
319  (*vi).Q() = math::Abs( 2*(*vi).Kh() );
320  else
321  (*vi).Q() = 2*math::Sqrt(math::Abs( (*vi).Kh()*(*vi).Kh() - (*vi).Kg()));
322  }
323 }
324 
325 /*
326  * RMS Curvature = sqrt(4H^2-2K)
327  * def and formula taken from
328  *
329  * Improved curvature estimation for watershed segmentation of 3-dimensional meshes
330  * S Pulla, A Razdan, G Farin - Arizona State University, Tech. Rep, 2001
331  */
332 static void VertexFromRMSCurvature(MeshType &m)
333 {
334  tri::RequirePerVertexQuality(m);
335  tri::RequirePerVertexCurvature(m);
336  VertexIterator vi;
337  for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
338  (*vi).Q() = math::Sqrt(math::Abs( 4*(*vi).Kh()*(*vi).Kh() - 2*(*vi).Kg()));
339 }
340 
341 /*
342  Saturate the Face quality so that for each vertex the gradient of the quality is lower than the given threshold value (in absolute value)
343  The saturation is done in a conservative way (quality is always decreased and never increased)
344 
345  Note: requires FF adjacency.
346  */
347 static void FaceSaturate(MeshType &m, FaceQualityType gradientThr=1.0)
348 {
349  tri::RequirePerFaceQuality(m);
350  tri::RequireFFAdjacency(m);
351  UpdateFlags<MeshType>::FaceClearV(m);
352  std::stack<FacePointer> st;
353 
354  st.push(&*m.face.begin());
355 
356  while(!st.empty())
357  {
358  FacePointer fc = st.top(); // the center
359  //printf("Stack size %i\n",st.size());
360  //printf("Pop elem %i %f\n",st.top() - &*m.vert.begin(), st.top()->Q());
361  st.pop();
362  fc->SetV();
363  std::vector<FacePointer> star;
364  typename std::vector<FacePointer>::iterator ffi;
365  for (int i=0;i<3;i++)
366  {
367  FacePointer fnext=fc->FFp(i);
368  if (fnext!=fc)star.push_back(fnext);
369  }
370  CoordType bary0=(fc->P(0)+fc->P(1)+fc->P(2))/3;
371  for(ffi=star.begin();ffi!=star.end();++ffi )
372  {
373  assert(fc!=(*ffi));
374  FaceQualityType &qi = (*ffi)->Q();
375  CoordType bary1=((*ffi)->P(0)+(*ffi)->P(1)+(*ffi)->P(2))/3;
376  FaceQualityType distGeom = Distance(bary0,bary1) / gradientThr;
377  // Main test if the quality varies more than the geometric displacement we have to lower something.
378  if( distGeom < fabs(qi - fc->Q()))
379  {
380  // center = 0 other=10 -> other =
381  // center = 10 other=0
382  if(fc->Q() > qi) // first case: the center of the star has to be lowered (and re-inserted in the queue).
383  {
384  //printf("Reinserting center %i \n",vc - &*m.vert.begin());
385  fc->Q() = qi+distGeom-(ScalarType)0.00001;
386  assert( distGeom > fabs(qi - fc->Q()));
387  st.push(fc);
388  break;
389  }
390  else
391  {
392  // second case: you have to lower qi, the vertex under examination.
393  assert( distGeom < fabs(qi - fc->Q()));
394  assert(fc->Q() < qi);
395  FaceQualityType newQi = fc->Q() + distGeom -(FaceQualityType)0.00001;
396  assert(newQi <= qi);
397  assert(fc->Q() < newQi);
398  assert( distGeom > fabs(newQi - fc->Q()) );
399 // printf("distGeom %f, qi %f, vc->Q() %f, fabs(qi - vc->Q()) %f\n",distGeom,qi,vc->Q(),fabs(qi - vc->Q()));
400  qi = newQi;
401  (*ffi)->ClearV();
402  }
403  }
404  if(!(*ffi)->IsV())
405  {
406  st.push( *ffi);
407 // printf("Reinserting side %i \n",*vvi - &*m.vert.begin());
408  (*ffi)->SetV();
409  }
410  }
411  }
412  }
413 
421 static void VertexSaturate(MeshType &m, ScalarType gradientThr=1.0)
422 {
423  tri::RequirePerVertexQuality(m);
424  tri::RequireVFAdjacency(m);
425 
427  std::stack<VertexPointer> st;
428 
429  st.push(&*m.vert.begin());
430 
431  while(!st.empty())
432  {
433  VertexPointer vc = st.top(); // the center
434  //printf("Stack size %i\n",st.size());
435  //printf("Pop elem %i %f\n",st.top() - &*m.vert.begin(), st.top()->Q());
436  st.pop();
437  vc->SetV();
438  std::vector<VertexPointer> star;
439  face::VVStarVF<FaceType>(vc,star);
440  for(auto vvi=star.begin();vvi!=star.end();++vvi )
441  {
442  ScalarType qi = (*vvi)->Q();
443  ScalarType distGeom = Distance((*vvi)->cP(),vc->cP()) / gradientThr;
444  // Main test if the quality varies more than the geometric displacement we have to lower something.
445  if( distGeom < fabs(qi - vc->Q()))
446  {
447  // center = 0 other=10 -> other =
448  // center = 10 other=0
449  if(vc->Q() > qi) // first case: the center of the star has to be lowered (and re-inserted in the queue).
450  {
451  //printf("Reinserting center %i \n",vc - &*m.vert.begin());
452  VertexQualityType delta=std::min(ScalarType(0.00001),ScalarType(distGeom/2.0));
453  vc->Q() = VertexQualityType(qi+distGeom-delta);
454  assert( distGeom > fabs(qi - vc->Q()));
455  st.push(vc);
456  break;
457  }
458  else
459  {
460  // second case: you have to lower qi, the vertex under examination.
461  assert( distGeom < fabs(qi - vc->Q()));
462  assert(vc->Q() < qi);
463  VertexQualityType delta=std::min(VertexQualityType(0.00001),VertexQualityType(distGeom/2.0));
464  VertexQualityType newQi = vc->Q() + distGeom -delta;
465  assert(newQi <= qi);
466  assert(vc->Q() < newQi);
467  assert( distGeom > fabs(newQi - vc->Q()) );
468 // printf("distGeom %f, qi %f, vc->Q() %f, fabs(qi - vc->Q()) %f\n",distGeom,qi,vc->Q(),fabs(qi - vc->Q()));
469  qi = newQi;
470  (*vvi)->ClearV();
471  }
472  }
473  if(!(*vvi)->IsV())
474  {
475  st.push( *vvi);
476 // printf("Reinserting side %i \n",*vvi - &*m.vert.begin());
477  (*vvi)->SetV();
478  }
479  }
480  }
481  }
482 
483 
484 }; //end class
485 } // end namespace
486 } // end namespace
487 #endif