PolyMeshT.cc 15.7 KB
Newer Older
Jan Möbius's avatar
Jan Möbius committed
1
/* ========================================================================= *
2
3
 *                                                                           *
 *                               OpenMesh                                    *
Jan Möbius's avatar
Jan Möbius committed
4
 *           Copyright (c) 2001-2015, RWTH-Aachen University                 *
Jan Möbius's avatar
Typo    
Jan Möbius committed
5
 *           Department of Computer Graphics and Multimedia                  *
Jan Möbius's avatar
Jan Möbius committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
 *                          All rights reserved.                             *
 *                            www.openmesh.org                               *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * This file is part of OpenMesh.                                            *
 *---------------------------------------------------------------------------*
 *                                                                           *
 * Redistribution and use in source and binary forms, with or without        *
 * modification, are permitted provided that the following conditions        *
 * are met:                                                                  *
 *                                                                           *
 * 1. Redistributions of source code must retain the above copyright notice, *
 *    this list of conditions and the following disclaimer.                  *
 *                                                                           *
 * 2. Redistributions in binary form must reproduce the above copyright      *
 *    notice, this list of conditions and the following disclaimer in the    *
 *    documentation and/or other materials provided with the distribution.   *
 *                                                                           *
 * 3. Neither the name of the copyright holder nor the names of its          *
 *    contributors may be used to endorse or promote products derived from   *
 *    this software without specific prior written permission.               *
 *                                                                           *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       *
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A           *
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,  *
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       *
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        *
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      *
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        *
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
Jan Möbius's avatar
Jan Möbius committed
39
40
 *                                                                           *
 * ========================================================================= */
41

42

Jan Möbius's avatar
Jan Möbius committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58


//=============================================================================
//
//  CLASS PolyMeshT - IMPLEMENTATION
//
//=============================================================================


#define OPENMESH_POLYMESH_C


//== INCLUDES =================================================================

#include <OpenMesh/Core/Mesh/PolyMeshT.hh>
#include <OpenMesh/Core/Geometry/LoopSchemeMaskT.hh>
59
#include <OpenMesh/Core/Utils/GenProg.hh>
Jan Möbius's avatar
Jan Möbius committed
60
#include <OpenMesh/Core/Utils/vector_cast.hh>
61
#include <OpenMesh/Core/Utils/vector_traits.hh>
Jan Möbius's avatar
Jan Möbius committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <OpenMesh/Core/System/omstream.hh>
#include <vector>


//== NAMESPACES ===============================================================


namespace OpenMesh {

//== IMPLEMENTATION ==========================================================

template <class Kernel>
uint PolyMeshT<Kernel>::find_feature_edges(Scalar _angle_tresh)
{
  assert(Kernel::has_edge_status());//this function needs edge status property
  uint n_feature_edges = 0;
  for (EdgeIter e_it = Kernel::edges_begin(); e_it != Kernel::edges_end(); ++e_it)
  {
Jan Möbius's avatar
OM3 fix    
Jan Möbius committed
80
    if (fabs(calc_dihedral_angle(*e_it)) > _angle_tresh)
Jan Möbius's avatar
Jan Möbius committed
81
    {//note: could be optimized by comparing cos(dih_angle) vs. cos(_angle_tresh)
Matthias Möller's avatar
Matthias Möller committed
82
      this->status(*e_it).set_feature(true);
Jan Möbius's avatar
Jan Möbius committed
83
84
85
86
      n_feature_edges++;
    }
    else
    {
Matthias Möller's avatar
Matthias Möller committed
87
      this->status(*e_it).set_feature(false);
Jan Möbius's avatar
Jan Möbius committed
88
89
90
91
92
93
94
95
96
    }
  }
  return n_feature_edges;
}

//-----------------------------------------------------------------------------

template <class Kernel>
typename PolyMeshT<Kernel>::Normal
97
98
99
PolyMeshT<Kernel>::calc_face_normal(FaceHandle _fh) const
{
  return calc_face_normal_impl(_fh, typename GenProg::IF<
100
    vector_traits<PolyMeshT<Kernel>::Point>::size_ == 3,
101
102
103
104
105
106
107
108
    PointIs3DTag,
    PointIsNot3DTag
  >::Result());
}

template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::calc_face_normal_impl(FaceHandle _fh, PointIs3DTag) const
Jan Möbius's avatar
Jan Möbius committed
109
{
110
  assert(this->halfedge_handle(_fh).is_valid());
111
  ConstFaceVertexIter fv_it(this->cfv_iter(_fh));
112
  
113
114
115
116
117
  // Safeguard for 1-gons
  if (!(++fv_it).is_valid()) return Normal(0, 0, 0);

  // Safeguard for 2-gons
  if (!(++fv_it).is_valid()) return Normal(0, 0, 0);
118
119

  // use Newell's Method to compute the surface normal
120
  Normal n(0,0,0);
121
  for(fv_it = this->cfv_iter(_fh); fv_it.is_valid(); ++fv_it)
122
  {
123
124
125
126
127
128
129
130
131
132
133
    // next vertex
    ConstFaceVertexIter fv_itn = fv_it;
    ++fv_itn;

    if (!fv_itn.is_valid())
      fv_itn = this->cfv_iter(_fh);

    // http://www.opengl.org/wiki/Calculating_a_Surface_Normal
    const Point a = this->point(*fv_it) - this->point(*fv_itn);
    const Point b = this->point(*fv_it) + this->point(*fv_itn);

Jan Möbius's avatar
Jan Möbius committed
134
135
136

    // Due to traits, the value types of normals and points can be different.
    // Therefore we cast them here.
137
138
139
    n[0] += static_cast<typename vector_traits<Normal>::value_type>(a[1] * b[2]);
    n[1] += static_cast<typename vector_traits<Normal>::value_type>(a[2] * b[0]);
    n[2] += static_cast<typename vector_traits<Normal>::value_type>(a[0] * b[1]);
140
  }
Jan Möbius's avatar
Jan Möbius committed
141

142
  const typename vector_traits<Normal>::value_type length = norm(n);
143
144
145
146
  
  // The expression ((n *= (1.0/norm)),n) is used because the OpenSG
  // vector class does not return self after component-wise
  // self-multiplication with a scalar!!!
147
148
  return (length != typename vector_traits<Normal>::value_type(0))
          ? ((n *= (typename vector_traits<Normal>::value_type(1)/length)), n)
149
          : Normal(0, 0, 0);
Jan Möbius's avatar
Jan Möbius committed
150
151
}

152
153
154
155
156
template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::calc_face_normal_impl(FaceHandle, PointIsNot3DTag) const
{
  // Dummy fallback implementation
157
158
159
160
161
162
163
  // Returns just an initialized all 0 normal
  // This function is only used if we don't hate a matching implementation
  // for normal computation with the current vector type defined in the mesh traits

  assert(false);

  Normal normal;
164
  vectorize(normal,static_cast<typename Normal::value_type>(0.0));
165
  return normal;
166
}
Jan Möbius's avatar
Jan Möbius committed
167

168
//-----------------------------------------------------------------------------
Jan Möbius's avatar
Jan Möbius committed
169
170
171
172
173
174
175

template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::
calc_face_normal(const Point& _p0,
     const Point& _p1,
     const Point& _p2) const
176
177
{
  return calc_face_normal_impl(_p0, _p1, _p2, typename GenProg::IF<
178
    vector_traits<PolyMeshT<Kernel>::Point>::size_ == 3,
179
180
181
182
183
184
185
186
187
188
189
190
    PointIs3DTag,
    PointIsNot3DTag
  >::Result());
}

template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::
calc_face_normal_impl(const Point& _p0,
     const Point& _p1,
     const Point& _p2,
     PointIs3DTag) const
Jan Möbius's avatar
Jan Möbius committed
191
192
193
194
195
{
#if 1
  // The OpenSG <Vector>::operator -= () does not support the type Point
  // as rhs. Therefore use vector_cast at this point!!!
  // Note! OpenSG distinguishes between Normal and Point!!!
196
197
  Normal p1p0(vector_cast<Normal>(_p0));  p1p0 -= vector_cast<Normal>(_p1);
  Normal p1p2(vector_cast<Normal>(_p2));  p1p2 -= vector_cast<Normal>(_p1);
Jan Möbius's avatar
Jan Möbius committed
198
199

  Normal n    = cross(p1p2, p1p0);
200
  typename vector_traits<Normal>::value_type length = norm(n);
Jan Möbius's avatar
Jan Möbius committed
201
202
203
204

  // The expression ((n *= (1.0/norm)),n) is used because the OpenSG
  // vector class does not return self after component-wise
  // self-multiplication with a scalar!!!
205
206
207
  return (length != typename vector_traits<Normal>::value_type(0))
          ? ((n *= (typename vector_traits<Normal>::value_type(1)/length)),n)
          : Normal(0,0,0);
Jan Möbius's avatar
Jan Möbius committed
208
209
210
211
212
#else
  Point p1p0 = _p0;  p1p0 -= _p1;
  Point p1p2 = _p2;  p1p2 -= _p1;

  Normal n = vector_cast<Normal>(cross(p1p2, p1p0));
213
  typename vector_traits<Normal>::value_type length = norm(n);
Jan Möbius's avatar
Jan Möbius committed
214

215
  return (length != 0.0) ? n *= (1.0/length) : Normal(0,0,0);
Jan Möbius's avatar
Jan Möbius committed
216
217
218
#endif
}

219
220
template <class Kernel>
typename PolyMeshT<Kernel>::Normal
221
PolyMeshT<Kernel>::calc_face_normal_impl(const Point&, const Point&, const Point&, PointIsNot3DTag) const
222
{
223
224
225
226
227
228
229
230
231

  // Dummy fallback implementation
  // Returns just an initialized all 0 normal
  // This function is only used if we don't hate a matching implementation
  // for normal computation with the current vector type defined in the mesh traits

  assert(false);

  Normal normal;
232
  vectorize(normal,static_cast<typename Normal::value_type>(0.0));
233
  return normal;
234
235
}

Jan Möbius's avatar
Jan Möbius committed
236
237
238
//-----------------------------------------------------------------------------

template <class Kernel>
239
typename PolyMeshT<Kernel>::Point
Jan Möbius's avatar
Jan Möbius committed
240
PolyMeshT<Kernel>::
241
calc_face_centroid(FaceHandle _fh) const
Jan Möbius's avatar
Jan Möbius committed
242
{
243
  Point _pt;
244
  vectorize(_pt, static_cast<typename Normal::value_type>(0.0));
Jan Möbius's avatar
Jan Möbius committed
245
  Scalar valence = 0.0;
Jan Möbius's avatar
Jan Möbius committed
246
  for (ConstFaceVertexIter cfv_it = this->cfv_iter(_fh); cfv_it.is_valid(); ++cfv_it, valence += 1.0)
Jan Möbius's avatar
Jan Möbius committed
247
  {
Jan Möbius's avatar
Jan Möbius committed
248
    _pt += this->point(*cfv_it);
Jan Möbius's avatar
Jan Möbius committed
249
250
  }
  _pt /= valence;
251
  return _pt;
Jan Möbius's avatar
Jan Möbius committed
252
253
254
255
256
257
258
259
260
}
//-----------------------------------------------------------------------------


template <class Kernel>
void
PolyMeshT<Kernel>::
update_normals()
{
261
  // Face normals are required to compute the vertex and the halfedge normals
262
  if (Kernel::has_face_normals() ) {
263
264
265
266
267
    update_face_normals();

    if (Kernel::has_vertex_normals() ) update_vertex_normals();
    if (Kernel::has_halfedge_normals()) update_halfedge_normals();
  }
Jan Möbius's avatar
Jan Möbius committed
268
269
270
271
272
273
274
275
276
277
278
}


//-----------------------------------------------------------------------------


template <class Kernel>
void
PolyMeshT<Kernel>::
update_face_normals()
{
279
  FaceIter f_it(Kernel::faces_sbegin()), f_end(Kernel::faces_end());
Jan Möbius's avatar
Jan Möbius committed
280
281

  for (; f_it != f_end; ++f_it)
Jan Möbius's avatar
Jan Möbius committed
282
    this->set_normal(*f_it, calc_face_normal(*f_it));
Jan Möbius's avatar
Jan Möbius committed
283
284
285
286
287
288
}


//-----------------------------------------------------------------------------


289
290
291
292
293
294
295
296
template <class Kernel>
void
PolyMeshT<Kernel>::
update_halfedge_normals(const double _feature_angle)
{
  HalfedgeIter h_it(Kernel::halfedges_begin()), h_end(Kernel::halfedges_end());

  for (; h_it != h_end; ++h_it)
297
    this->set_normal(*h_it, calc_halfedge_normal(*h_it, _feature_angle));
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
}


//-----------------------------------------------------------------------------


template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::
calc_halfedge_normal(HalfedgeHandle _heh, const double _feature_angle) const
{
  if(Kernel::is_boundary(_heh))
    return Normal(0,0,0);
  else
  {
    std::vector<FaceHandle> fhs; fhs.reserve(10);

    HalfedgeHandle heh = _heh;

    // collect CW face-handles
    do
    {
      fhs.push_back(Kernel::face_handle(heh));

      heh = Kernel::next_halfedge_handle(heh);
      heh = Kernel::opposite_halfedge_handle(heh);
    }
    while(heh != _heh && !Kernel::is_boundary(heh) && !is_estimated_feature_edge(heh, _feature_angle));

    // collect CCW face-handles
    if(heh != _heh && !is_estimated_feature_edge(_heh, _feature_angle))
    {
      heh = Kernel::opposite_halfedge_handle(_heh);

332
333
334
      if ( !Kernel::is_boundary(heh) ) {
        do
        {
335

336
337
338
339
340
341
          fhs.push_back(Kernel::face_handle(heh));

          heh = Kernel::prev_halfedge_handle(heh);
          heh = Kernel::opposite_halfedge_handle(heh);
        }
        while(!Kernel::is_boundary(heh) && !is_estimated_feature_edge(heh, _feature_angle));
342
343
344
345
346
347
348
      }
    }

    Normal n(0,0,0);
    for(unsigned int i=0; i<fhs.size(); ++i)
      n += Kernel::normal(fhs[i]);

349
    return normalize(n);
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
  }
}


//-----------------------------------------------------------------------------


template <class Kernel>
bool
PolyMeshT<Kernel>::
is_estimated_feature_edge(HalfedgeHandle _heh, const double _feature_angle) const
{
  EdgeHandle eh = Kernel::edge_handle(_heh);

  if(Kernel::has_edge_status())
  {
    if(Kernel::status(eh).feature())
      return true;
  }

  if(Kernel::is_boundary(eh))
    return false;

  // compute angle between faces
  FaceHandle fh0 = Kernel::face_handle(_heh);
  FaceHandle fh1 = Kernel::face_handle(Kernel::opposite_halfedge_handle(_heh));

  Normal fn0 = Kernel::normal(fh0);
  Normal fn1 = Kernel::normal(fh1);

  // dihedral angle above angle threshold
381
  return ( dot(fn0,fn1) < cos(_feature_angle) );
382
383
384
385
386
387
}


//-----------------------------------------------------------------------------


Jan Möbius's avatar
Jan Möbius committed
388
389
390
391
392
393
394
395
template <class Kernel>
typename PolyMeshT<Kernel>::Normal
PolyMeshT<Kernel>::
calc_vertex_normal(VertexHandle _vh) const
{
  Normal n;
  calc_vertex_normal_fast(_vh,n);

396
397
  Scalar length = norm(n);
  if (length != 0.0) n *= (Scalar(1.0)/length);
Jan Möbius's avatar
Jan Möbius committed
398
399
400
401
402
403
404
405
406

  return n;
}

//-----------------------------------------------------------------------------
template <class Kernel>
void PolyMeshT<Kernel>::
calc_vertex_normal_fast(VertexHandle _vh, Normal& _n) const
{
407
  vectorize(_n, static_cast<typename Normal::value_type>(0.0) );
Jan Möbius's avatar
   
Jan Möbius committed
408
  for (ConstVertexFaceIter vf_it = this->cvf_iter(_vh); vf_it.is_valid(); ++vf_it)
409
    _n += this->normal(*vf_it);
Jan Möbius's avatar
Jan Möbius committed
410
411
412
413
414
415
416
}

//-----------------------------------------------------------------------------
template <class Kernel>
void PolyMeshT<Kernel>::
calc_vertex_normal_correct(VertexHandle _vh, Normal& _n) const
{
417
  vectorize(_n, static_cast<typename Normal::value_type>(0.0) );
Jan Möbius's avatar
   
Jan Möbius committed
418
419
  ConstVertexIHalfedgeIter cvih_it = this->cvih_iter(_vh);
  if (! cvih_it.is_valid() )
Jan Möbius's avatar
Jan Möbius committed
420
421
422
423
  {//don't crash on isolated vertices
    return;
  }
  Normal in_he_vec;
Jan Möbius's avatar
   
Jan Möbius committed
424
425
  calc_edge_vector(*cvih_it, in_he_vec);
  for ( ; cvih_it.is_valid(); ++cvih_it)
Jan Möbius's avatar
Jan Möbius committed
426
  {//calculates the sector normal defined by cvih_it and adds it to _n
Jan Möbius's avatar
   
Jan Möbius committed
427
    if (this->is_boundary(*cvih_it))
Jan Möbius's avatar
Jan Möbius committed
428
429
430
    {
      continue;
    }
Jan Möbius's avatar
   
Jan Möbius committed
431
    HalfedgeHandle out_heh(this->next_halfedge_handle(*cvih_it));
Jan Möbius's avatar
Jan Möbius committed
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
    Normal out_he_vec;
    calc_edge_vector(out_heh, out_he_vec);
    _n += cross(in_he_vec, out_he_vec);//sector area is taken into account
    in_he_vec = out_he_vec;
    in_he_vec *= -1;//change the orientation
  }
}

//-----------------------------------------------------------------------------
template <class Kernel>
void PolyMeshT<Kernel>::
calc_vertex_normal_loop(VertexHandle _vh, Normal& _n) const
{
  static const LoopSchemeMaskDouble& loop_scheme_mask__ =
                  LoopSchemeMaskDoubleSingleton::Instance();

  Normal t_v(0.0,0.0,0.0), t_w(0.0,0.0,0.0);
Jan Möbius's avatar
   
Jan Möbius committed
449
  unsigned int vh_val = this->valence(_vh);
Jan Möbius's avatar
Jan Möbius committed
450
  unsigned int i = 0;
Jan Möbius's avatar
   
Jan Möbius committed
451
  for (ConstVertexOHalfedgeIter cvoh_it = this->cvoh_iter(_vh); cvoh_it.is_valid(); ++cvoh_it, ++i)
Jan Möbius's avatar
Jan Möbius committed
452
  {
Jan Möbius's avatar
   
Jan Möbius committed
453
    VertexHandle r1_v( this->to_vertex_handle(*cvoh_it) );
454
455
    t_v += (typename vector_traits<Point>::value_type)(loop_scheme_mask__.tang0_weight(vh_val, i))*this->point(r1_v);
    t_w += (typename vector_traits<Point>::value_type)(loop_scheme_mask__.tang1_weight(vh_val, i))*this->point(r1_v);
Jan Möbius's avatar
Jan Möbius committed
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
  }
  _n = cross(t_w, t_v);//hack: should be cross(t_v, t_w), but then the normals are reversed?
}

//-----------------------------------------------------------------------------


template <class Kernel>
void
PolyMeshT<Kernel>::
update_vertex_normals()
{
  VertexIter  v_it(Kernel::vertices_begin()), v_end(Kernel::vertices_end());

  for (; v_it!=v_end; ++v_it)
471
    this->set_normal(*v_it, calc_vertex_normal(*v_it));
Jan Möbius's avatar
Jan Möbius committed
472
473
474
475
476
}

//=============================================================================
} // namespace OpenMesh
//=============================================================================