PolyMeshT.cc 15.5 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,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,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, 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, 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, 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
//=============================================================================