1 // This is brl/bbas/bmsh3d/bmsh3d_face.cxx
2 //---------------------------------------------------------------------
3 #include <sstream>
4 #include <iostream>
5 #include <cstdio>
6 #include "bmsh3d_face.h"
7 //:
8 // \file
9 // \brief 3D mesh face.
10 //
11 // \author
12 // MingChing Chang Apr 22, 2005
13 //
14 // \verbatim
15 // Modifications
16 // <none>
17 // \endverbatim
18 //
19 //-------------------------------------------------------------------------
20
21 #include <cassert>
22 #ifdef _MSC_VER
23 # include "vcl_msvc_warnings.h"
24 #endif
25 #include "vgl/vgl_point_2d.h"
26 #include "vgl/vgl_distance.h"
27
28 #include "bmsh3d_triangle.h"
29 #include "bmsh3d_mesh.h"
30
31 //###############################################################
32 //###### Connectivity Query Functions ######
33 //###############################################################
34
get_incident_HEs(std::vector<bmsh3d_halfedge * > & incident_HEs) const35 void bmsh3d_face::get_incident_HEs (std::vector<bmsh3d_halfedge*>& incident_HEs) const
36 {
37 bmsh3d_halfedge* HE = halfedge_;
38 if (HE == nullptr)
39 return;
40 do {
41 assert (HE != nullptr);
42 incident_HEs.push_back (HE);
43 HE = HE->next();
44 }
45 while (HE != halfedge_ && HE != nullptr);
46 }
47
get_incident_Es(std::vector<bmsh3d_edge * > & incident_Es) const48 void bmsh3d_face::get_incident_Es (std::vector<bmsh3d_edge*>& incident_Es) const
49 {
50 bmsh3d_halfedge* HE = halfedge_;
51 do {
52 incident_Es.push_back (HE->edge());
53 HE = HE->next();
54 }
55 while (HE != halfedge_);
56 }
57
n_incident_Es() const58 unsigned int bmsh3d_face::n_incident_Es () const
59 {
60 unsigned int count = 0;
61 bmsh3d_halfedge* HE = halfedge_;
62 do {
63 count++;
64 HE = HE->next();
65 }
66 while (HE != halfedge_ && HE != nullptr);
67 return count;
68 }
69
is_E_incident(const bmsh3d_edge * inputE) const70 bool bmsh3d_face::is_E_incident (const bmsh3d_edge* inputE) const
71 {
72 bmsh3d_halfedge* HE = halfedge_;
73 // If the next is NULL, it is a loop curve.
74 if (HE->next() == nullptr) {
75 if (HE->edge() == inputE)
76 return true;
77 else
78 return false;
79 }
80 // Traverse through the circular list of halfedges,
81 // and find the vertex incident with both HE and nextHE
82 do {
83 if (HE->edge() == inputE)
84 return true;
85 HE = HE->next();
86 }
87 while (HE != halfedge_);
88 return false;
89 }
90
find_bnd_HE() const91 bmsh3d_halfedge* bmsh3d_face::find_bnd_HE () const
92 {
93 bmsh3d_halfedge* HE = halfedge_;
94 if (HE->next() == nullptr) {
95 if (HE->pair() == nullptr)
96 return HE;
97 else
98 return nullptr;
99 }
100 do { // Traverse through the circular list of halfedges.
101 if (HE->pair() == nullptr)
102 return HE;
103 HE = HE->next();
104 }
105 while (HE != halfedge_);
106 return nullptr;
107 }
108
is_V_incident_via_HE(const bmsh3d_vertex * inputV) const109 bool bmsh3d_face::is_V_incident_via_HE (const bmsh3d_vertex* inputV) const
110 {
111 bmsh3d_halfedge* HE = halfedge_;
112 // if the next is NULL, it is a loop curve.
113 // this will not happen for the fullshock mesh.
114 if (HE->next() == nullptr) {
115 if (HE->edge()->sV() == inputV || HE->edge()->eV() == inputV)
116 return true;
117 else
118 return false;
119 }
120 do { // traverse through the circular list of halfedges
121 if (HE->edge()->sV() == inputV || HE->edge()->eV() == inputV )
122 return true;
123 HE = HE->next();
124 }
125 while (HE != halfedge_);
126 return false;
127 }
128
129 //: loop through the halfedges and locate the inputV
get_next_V_via_HE(const bmsh3d_vertex * inputV) const130 bmsh3d_vertex* bmsh3d_face::get_next_V_via_HE (const bmsh3d_vertex* inputV) const
131 {
132 bmsh3d_halfedge* HE = halfedge_;
133 // if the next is NULL, it is a loop curve.
134 if (HE->next() == nullptr)
135 return nullptr;
136 do { // traverse through the circular list of halfedges,
137 bmsh3d_halfedge* nextHE = HE->next();
138 // find the vertex incident with both HE and nextHE
139 bmsh3d_vertex* vertex = incident_V_of_Es (HE->edge(), nextHE->edge());
140 if (vertex == inputV)
141 return nextHE->edge()->other_V (vertex);
142 HE = HE->next();
143 }
144 while (HE != halfedge_);
145 assert (0);
146 return nullptr;
147 }
148
149 //: Given a vertex V and an edge of this face incident to V, find the other edge of this face incident to V.
find_other_E(const bmsh3d_vertex * inputV,const bmsh3d_edge * inputE) const150 bmsh3d_edge* bmsh3d_face::find_other_E (const bmsh3d_vertex* inputV,
151 const bmsh3d_edge* inputE) const
152 {
153 bmsh3d_halfedge* HE = halfedge_;
154 do {
155 bmsh3d_edge* E = HE->edge();
156 if (E != inputE && E->is_V_incident(inputV))
157 return E;
158 HE = HE->next();
159 }
160 while (HE != halfedge_);
161 assert (0);
162 return nullptr;
163 }
164
165 //: Given a vertex V and a halfedge of this face incident to V, find the other halfedge of this face incident of V.
find_other_HE(const bmsh3d_vertex * inputV,const bmsh3d_halfedge * inputHE) const166 bmsh3d_halfedge* bmsh3d_face::find_other_HE (const bmsh3d_vertex* inputV,
167 const bmsh3d_halfedge* inputHE) const
168 {
169 bmsh3d_halfedge* HE = halfedge_;
170 do {
171 bmsh3d_edge* E = HE->edge();
172 if (HE != inputHE && E->is_V_incident(inputV))
173 return HE;
174 HE = HE->next();
175 }
176 while (HE != halfedge_);
177 assert (0);
178 return nullptr;
179 }
180
181 //: Given a vertex V and an edge of this face incident to V, find the next edge (following the circular halfedge list) of this face incident to V.
find_next_E(const bmsh3d_vertex * inputV,const bmsh3d_edge * inputE) const182 bmsh3d_edge* bmsh3d_face::find_next_E (const bmsh3d_vertex* inputV,
183 const bmsh3d_edge* inputE) const
184 {
185 if (halfedge_ == nullptr)
186 return nullptr;
187 if (halfedge_->next() == nullptr)
188 return nullptr;
189
190 // traverse through the circular list of halfedges
191 bmsh3d_halfedge* HE = halfedge_;
192 do {
193 bmsh3d_halfedge* nextHE = HE->next();
194 bmsh3d_edge* E = HE->edge();
195 bmsh3d_edge* nextE = nextHE->edge();
196 if (E == inputE) {
197 if (nextE->is_V_incident (inputV))
198 return nextE;
199 }
200 else if (nextE == inputE) {
201 if (E->is_V_incident (inputV))
202 return E;
203 }
204 HE = HE->next();
205 }
206 while (HE != halfedge_);
207 assert (0);
208 return nullptr;
209 }
210
211 //: Given a vertex V and a halfedge of this face incident to V, find the next halfedge (following the circular halfedge list) of this face incident to V.
212 // Note that there can exist multiple answers if inputV is multiply (>2) connected.
213 // This function returns the first qualify halfedge.
find_next_HE(const bmsh3d_vertex * inputV,const bmsh3d_halfedge * inputHE) const214 bmsh3d_halfedge* bmsh3d_face::find_next_HE(const bmsh3d_vertex* inputV,
215 const bmsh3d_halfedge* inputHE) const
216 {
217 if (halfedge_ == nullptr)
218 return nullptr;
219 if (halfedge_->next() == nullptr)
220 return nullptr;
221
222 // traverse through the circular list of halfedges
223 bmsh3d_halfedge* HE = halfedge_;
224 do {
225 bmsh3d_halfedge* nextHE = HE->next();
226 bmsh3d_edge* E = HE->edge();
227 bmsh3d_edge* nextE = nextHE->edge();
228 if (HE == inputHE) {
229 if (nextE->is_V_incident (inputV))
230 return nextHE;
231 }
232 else if (nextHE == inputHE) {
233 if (E->is_V_incident (inputV))
234 return HE;
235 }
236 HE = HE->next();
237 }
238 while (HE != halfedge_);
239 assert (0);
240 return nullptr;
241 }
242
243 //: get the two edges incident at this vertex and compute their angle
angle_at_V(const bmsh3d_vertex * inputV) const244 double bmsh3d_face::angle_at_V (const bmsh3d_vertex* inputV) const
245 {
246 bmsh3d_edge* E1 = nullptr;
247 bmsh3d_edge* E2 = nullptr;
248
249 // Get the two edges of the face incident at the inputV
250 bmsh3d_halfedge* HE = halfedge_;
251 do {
252 bmsh3d_edge* E = HE->edge();
253 if (E->is_V_incident (inputV)) {
254 if (E1 == nullptr)
255 E1 = E;
256 else if (E2 == nullptr) {
257 E2 = E;
258 }
259 else
260 assert (0);
261 }
262 HE = HE->next();
263 }
264 while (HE != halfedge_);
265
266 bmsh3d_vertex* v1 = E1->other_V (inputV);
267 bmsh3d_vertex* v2 = E2->other_V (inputV);
268 double a = vgl_distance (v1->pt(), v2->pt());
269 double b = E1->length();
270 double c = E2->length();
271 return std::acos ( (b*b + c*c - a*a)/(b*c*2) );
272 }
273
274 //: Return true if this face is incident to all given vertices.
n_incident_Vs_in_set(std::set<bmsh3d_vertex * > & vertices) const275 int bmsh3d_face::n_incident_Vs_in_set (std::set<bmsh3d_vertex*>& vertices) const
276 {
277 int n_inc_V = 0;
278 bmsh3d_halfedge* HE = halfedge_;
279 do {
280 bmsh3d_halfedge* nextHE = HE->next();
281 bmsh3d_vertex* V = Es_sharing_V (HE->edge(), nextHE->edge());
282 if (vertices.find (V) != vertices.end())
283 n_inc_V++;
284 HE = nextHE;
285 }
286 while (HE != halfedge_);
287 return n_inc_V;
288 }
289
290 //: Return true if this face is incident to all given vertices.
all_Vs_incident(std::vector<bmsh3d_vertex * > & vertices) const291 bool bmsh3d_face::all_Vs_incident (std::vector<bmsh3d_vertex*>& vertices) const
292 {
293 // Put all input vertices into the VSet.
294 std::set<bmsh3d_vertex*> VSet;
295 for (auto vertice : vertices)
296 VSet.insert (vertice);
297 assert (VSet.size() == vertices.size());
298
299 // Go through the halfedge loop and remove incident vertices from VSet.
300 bmsh3d_halfedge* HE = halfedge_;
301 do {
302 bmsh3d_edge* E = HE->edge();
303 VSet.erase (E->vertices(0));
304 VSet.erase (E->vertices(1));
305 HE = HE->next();
306 }
307 while (HE != halfedge_);
308
309 return VSet.empty();
310 }
311
get_ordered_Vs(std::vector<bmsh3d_vertex * > & vertices) const312 void bmsh3d_face::get_ordered_Vs (std::vector<bmsh3d_vertex*>& vertices) const
313 {
314 assert (vertices.size()==0);
315 if (vertices_.size() != 0)
316 _get_ordered_Vs_IFS(vertices);
317 else
318 _get_ordered_Vs_MHE(vertices);
319 assert(vertices.size() > 2);
320 }
321
_get_ordered_Vs_MHE(std::vector<bmsh3d_vertex * > & vertices) const322 void bmsh3d_face::_get_ordered_Vs_MHE(std::vector<bmsh3d_vertex*>& vertices) const
323 {
324 bmsh3d_halfedge* HE = halfedge_;
325 assert (HE->next());
326 // Traverse through the circular list of halfedges,
327 // and find the vertex incident with both HE and nextHE
328 do {
329 bmsh3d_halfedge* nextHE = HE->next();
330 bmsh3d_vertex* V = incident_V_of_Es (HE->edge(), nextHE->edge());
331 vertices.push_back (V);
332 HE = HE->next();
333 }
334 while (HE != halfedge_);
335 }
336
_get_ordered_Vs_IFS(std::vector<bmsh3d_vertex * > & vertices) const337 void bmsh3d_face::_get_ordered_Vs_IFS(std::vector<bmsh3d_vertex*>& vertices) const
338 {
339 for (auto vertice : vertices_)
340 vertices.push_back (vertice);
341 }
342
get_ordered_V_ids(std::vector<int> & vids) const343 void bmsh3d_face::get_ordered_V_ids (std::vector<int>& vids) const
344 {
345 assert (vids.size()==0);
346 if (vertices_.size() != 0)
347 _get_ordered_V_ids_IFS(vids);
348 else
349 _get_ordered_V_ids_MHE(vids);
350 assert (vids.size() > 2);
351 }
352
_get_ordered_V_ids_MHE(std::vector<int> & vids) const353 void bmsh3d_face::_get_ordered_V_ids_MHE(std::vector<int>& vids) const
354 {
355 bmsh3d_halfedge* HE = halfedge_;
356 assert (HE->next());
357 // Traverse through the circular list of halfedges,
358 // and find the vertex incident with both HE and nextHE
359 do {
360 bmsh3d_halfedge* nextHE = HE->next();
361 bmsh3d_vertex* V = incident_V_of_Es (HE->edge(), nextHE->edge());
362 vids.push_back (V->id());
363 HE = HE->next();
364 }
365 while (HE != halfedge_);
366 }
367
_get_ordered_V_ids_IFS(std::vector<int> & vids) const368 void bmsh3d_face::_get_ordered_V_ids_IFS(std::vector<int>& vids) const
369 {
370 for (auto vertice : vertices_)
371 vids.push_back (vertice->id());
372 }
373
374 //###############################################################
375 //###### Handle local list of incident vertices ######
376 //###############################################################
377
378 //: track IFS ordered vertices use the halfedge data structure.
379 // Store result in the vertices[] vector.
_ifs_track_ordered_vertices()380 void bmsh3d_face::_ifs_track_ordered_vertices()
381 {
382 // skip if the vertices_[] is non-empty.
383 if (vertices_.size() != 0)
384 return;
385 // the starting halfedge is the face's pointing halfedge
386 bmsh3d_halfedge* HE = halfedge_;
387 // if the next is NULL, it is a loop curve.
388 // this will not happen for the fullshock mesh.
389 if (HE->next() == nullptr) {
390 vertices_.push_back (HE->edge()->eV());
391 return;
392 }
393 // traverse through the circular list of halfedges,
394 // and find the vertex incident with both HE and nextHE
395 do {
396 bmsh3d_halfedge* nextHE = HE->next();
397 bmsh3d_vertex* V = incident_V_of_Es (HE->edge(), nextHE->edge());
398 vertices_.push_back (V);
399 HE = HE->next();
400 }
401 while (HE != halfedge_);
402 assert (vertices_.size() > 2);
403 }
404
405 // Test if the face's IFS structure is correct (repeated or wrong Vids).
_is_ifs_valid(bmsh3d_mesh * M)406 bool bmsh3d_face::_is_ifs_valid(bmsh3d_mesh* M)
407 {
408 // Check if there is repeated vertices.
409 std::set<bmsh3d_vertex*> Vset;
410 for (auto V : vertices_) {
411 Vset.insert (V);
412 }
413 if (Vset.size() != vertices_.size())
414 return false;
415
416 // Check if each V is inside M.
417 for (auto V : vertices_) {
418 if (! M->contains_V (V->id()))
419 return false;
420 }
421 return true;
422 }
423
_ifs_inside_box(const vgl_box_3d<double> & box) const424 bool bmsh3d_face::_ifs_inside_box(const vgl_box_3d<double>& box) const
425 {
426 for (auto V : vertices_) {
427 if (!box.contains (V->pt()))
428 return false;
429 }
430 return true;
431 }
432
_ifs_outside_box(const vgl_box_3d<double> & box) const433 bool bmsh3d_face::_ifs_outside_box(const vgl_box_3d<double>& box) const
434 {
435 for (auto V : vertices_) {
436 if (box.contains (V->pt()))
437 return false;
438 }
439 return true;
440 }
441
442 //###############################################################
443 //###### Geometry Query Functions ######
444 //###############################################################
445
446 //: if any vertex out of the box return false, otherwise return true.
447 // Use the halfedge_ to traverse the face.
448 //
is_inside_box(const vgl_box_3d<double> & box) const449 bool bmsh3d_face::is_inside_box (const vgl_box_3d<double>& box) const
450 {
451 bmsh3d_halfedge* HE = halfedge_;
452 do {
453 bmsh3d_halfedge* nextHE = HE->next();
454 bmsh3d_vertex* V = Es_sharing_V (HE->edge(), nextHE->edge());
455 if (!box.contains (V->pt()))
456 return false;
457 HE = nextHE;
458 }
459 while (HE != halfedge_);
460 return true;
461 }
462
463 //: if any vertex inside box return false, otherwise return true.
464 // Use the halfedge_ to traverse the face.
465 //
is_outside_box(const vgl_box_3d<double> & box) const466 bool bmsh3d_face::is_outside_box (const vgl_box_3d<double>& box) const
467 {
468 bmsh3d_halfedge* HE = halfedge_;
469 do {
470 bmsh3d_halfedge* nextHE = HE->next();
471 bmsh3d_vertex* V = Es_sharing_V (HE->edge(), nextHE->edge());
472 if (box.contains (V->pt()))
473 return false;
474 HE = nextHE;
475 }
476 while (HE != halfedge_);
477 return true;
478 }
479
compute_center_pt() const480 vgl_point_3d<double> bmsh3d_face::compute_center_pt () const
481 {
482 std::vector<bmsh3d_vertex*> vertices;
483 get_ordered_Vs (vertices);
484 return compute_cen (vertices);
485 }
486
compute_center_pt(const std::vector<bmsh3d_vertex * > & vertices) const487 vgl_point_3d<double> bmsh3d_face::compute_center_pt (const std::vector<bmsh3d_vertex*>& vertices) const
488 {
489 return compute_cen (vertices);
490 }
491
492
493 //: Compute face normal using the order of halfedges
494 // This function works for both convex and non-convex faces.
compute_normal()495 vgl_vector_3d<double> bmsh3d_face::compute_normal()
496 {
497 std::vector<bmsh3d_edge*> inc_edges;
498 this->get_incident_Es (inc_edges);
499 vgl_point_3d<double> centroid = this->compute_center_pt();
500 vgl_vector_3d<double> normal(0.0, 0.0, 0.0);
501
502 for (auto E : inc_edges) {
503 bmsh3d_halfedge* HE = E->get_HE_of_F (this);
504 auto* sV = (bmsh3d_vertex*)HE->s_vertex();
505 auto* eV = (bmsh3d_vertex*)HE->e_vertex();
506 vgl_point_3d<double> p0 = sV->pt();
507 vgl_point_3d<double> p1 = eV->pt();
508 vgl_vector_3d<double> v0 = p0 - centroid;
509 vgl_vector_3d<double> v1 = p1 - centroid;
510 normal += cross_product(v0,v1);
511 }
512 return normal;
513 }
514
515 //###############################################################
516 //###### Connectivity Modification Functions ######
517 //###############################################################
518
519 //: Connect a halfedge to this mesh face.
520 // Note that the link list is circular, but not necessarily geometrically ordered!
521 // Also be careful in the empty and starting cases.
_connect_HE_to_end(bmsh3d_halfedge * inputHE)522 void bmsh3d_face::_connect_HE_to_end(bmsh3d_halfedge* inputHE)
523 {
524 if (halfedge_ == nullptr) { // 1)
525 halfedge_ = inputHE;
526 return;
527 }
528 else if (halfedge_->next() == nullptr) { // 2) Only one halfedge there
529 halfedge_->set_next (inputHE);
530 inputHE->set_next (halfedge_);
531 return;
532 }
533 else { // 3) The general circular list case, add to the end!
534 bmsh3d_halfedge* HE = halfedge_;
535 while (HE->next() != halfedge_)
536 HE = HE->next();
537 HE->set_next (inputHE);
538 inputHE->set_next (halfedge_);
539 return;
540 }
541 }
542
543 //: remove the input halfedge from the face's halfedge list.
544 // (memory of inputHE not released)
545 // return true if success.
_remove_HE(bmsh3d_halfedge * inputHE)546 bool bmsh3d_face::_remove_HE(bmsh3d_halfedge* inputHE)
547 {
548 if (halfedge_ == nullptr) // 1)
549 return false;
550 else if (halfedge_->next() == nullptr) { // 2)
551 if (halfedge_ == inputHE) {
552 halfedge_ = nullptr;
553 return true;
554 }
555 else
556 return false;
557 }
558 else { // 3)
559 bmsh3d_halfedge* HE = halfedge_;
560 do {
561 bmsh3d_halfedge* nextHE = HE->next();
562 if (nextHE == inputHE) {
563 HE->set_next (nextHE->next()); // found. remove and fix pointers.
564 if (nextHE == halfedge_)
565 halfedge_ = HE; // fix the case if inputHE is halfedge_
566 return true;
567 }
568 HE = HE->next();
569 }
570 while (HE != halfedge_);
571 return false;
572 }
573 }
574
575 //: Create a halfedge and connect a mesh face and an edge.
connect_E_to_end(bmsh3d_edge * E)576 void bmsh3d_face::connect_E_to_end (bmsh3d_edge* E)
577 {
578 // The halfedge will be deleted when the sheet disconnect from the sheet.
579 auto* HE = new bmsh3d_halfedge (E, this);
580 // Handle the both-way connectivity of halfedge-face.
581 _connect_HE_to_end(HE);
582 // Handle the both-way connectivity of halfedge-edge.
583 E->_connect_HE_to_end(HE);
584 }
585
586 //: Disconnect a face and an edge (delete the halfedge).
disconnect_E(bmsh3d_halfedge * HE)587 void bmsh3d_face::disconnect_E (bmsh3d_halfedge* HE)
588 {
589 bmsh3d_edge* E = HE->edge();
590 // Disconnect HE from this face.
591 _remove_HE(HE);
592 // Disconnect HE from E.
593 E->_disconnect_HE(HE);
594 delete HE;
595 }
596
597 //: Sort the incident halfedges to form a circular list
_sort_HEs_circular()598 bool bmsh3d_face::_sort_HEs_circular()
599 {
600 if (halfedge_ == nullptr)
601 return false;
602
603 // put all halfedges into a vector
604 std::vector<bmsh3d_halfedge*> incident_HEs;
605 get_incident_HEs (incident_HEs);
606
607 if (incident_HEs.size() == 1)
608 return false; // Skip face with only one edge.
609 if (incident_HEs.size() == 2)
610 return false; // Skip face with only two edges.
611
612 // Now reset the halfedges in a correct order.
613 // The current halfedge_ is still the starting halfedge.
614 bmsh3d_halfedge* HE = halfedge_;
615 bmsh3d_edge* E = HE->edge();
616 bmsh3d_vertex* eV = E->eV();
617 do {
618 // Find the next halfedge incident to the eV from the vector
619 bmsh3d_halfedge* nextHE = _find_next_halfedge(HE, eV, incident_HEs);
620 HE->set_next (nextHE);
621
622 HE = HE->next();
623 E = HE->edge();
624 eV = E->other_V (eV);
625 }
626 while (HE != halfedge_);
627 return true;
628 }
629
630 //: disconnect all associated halfedges from their edges and delete them.
_discon_all_incident_Es()631 void bmsh3d_face::_discon_all_incident_Es ()
632 {
633 _delete_HE_chain (halfedge_);
634 }
635
636 //: reverse the orientation of chain of halfedges of this face.
_reverse_HE_chain()637 void bmsh3d_face::_reverse_HE_chain ()
638 {
639 if (halfedge_ == nullptr)
640 return;
641 if (halfedge_->next() == nullptr)
642 return;
643
644 std::vector<bmsh3d_halfedge*> HE_chain;
645 bmsh3d_halfedge* HE = halfedge_;
646 do {
647 HE_chain.push_back (HE);
648 HE = HE->next();
649 }
650 while (HE != halfedge_);
651
652 // Build the circular list of halfedges in reverse order.
653 HE = halfedge_;
654 for (int i=(int) HE_chain.size()-1; i>=0; i--) {
655 bmsh3d_halfedge* nextHE = HE_chain[i];
656 HE->set_next (nextHE);
657 HE = nextHE;
658 }
659 assert (HE == halfedge_);
660 }
661
set_orientation(bmsh3d_halfedge * new_start_he,bmsh3d_vertex * new_next_v)662 void bmsh3d_face::set_orientation(bmsh3d_halfedge* new_start_he,
663 bmsh3d_vertex* new_next_v)
664 {
665 if (halfedge_ == nullptr)
666 return;
667 if (halfedge_->next() == nullptr)
668 return;
669
670 // set the new_start_he (should check)
671 halfedge_ = new_start_he;
672
673 bmsh3d_halfedge* nextHE = new_start_he->next();
674 if (!nextHE->edge()->is_V_incident (new_next_v))
675 _reverse_HE_chain ();
676 }
677
678 //###############################################################
679 //###### Other functions ######
680 //###############################################################
681
getInfo(std::ostringstream & ostrm)682 void bmsh3d_face::getInfo (std::ostringstream& ostrm)
683 {
684 char s[1024];
685
686 bmsh3d_halfedge* HE = halfedge_;
687 bmsh3d_edge* e0 = HE->edge();
688 HE = HE->next();
689 bmsh3d_edge* e1 = HE->edge();
690 HE = HE->next();
691 bmsh3d_edge* e2 = HE->edge();
692 unsigned int c0 = e0->n_incident_Fs();
693 unsigned int c1 = e1->n_incident_Fs();
694 unsigned int c2 = e2->n_incident_Fs();
695
696 std::sprintf (s, "\n==============================\n"); ostrm<<s;
697 std::sprintf (s, "bmsh3d_face id: %d, type: %s (%u-%u-%u)\n",
698 id_, tri_get_topo_string().c_str(), c0, c1, c2); ostrm<<s;
699
700 // the incident edges via halfedges
701 int n_sides = n_incident_Es ();
702 std::sprintf (s, " %d incident edges in order: ", n_sides); ostrm<<s;
703
704 if (halfedge_ == nullptr) {
705 std::sprintf (s, "NONE "); ostrm<<s;
706 }
707 else if (halfedge_->next() == nullptr) {
708 std::sprintf (s, "%d ", (halfedge_)->edge()->id()); ostrm<<s;
709 }
710 else {
711 bmsh3d_halfedge* HE = halfedge_;
712 do {
713 std::sprintf (s, "%d ", HE->edge()->id()); ostrm<<s;
714 HE = HE->next();
715 }
716 while (HE != halfedge_);
717 }
718
719 // the incident vertices in order
720 std::sprintf (s, "\n %d incident vertices in order: ", n_sides); ostrm<<s;
721 if (halfedge_ == nullptr) {
722 std::sprintf (s, "NONE "); ostrm<<s;
723 }
724 else if (halfedge_->next() == nullptr) {
725 bmsh3d_halfedge* HE = halfedge_;
726 assert (HE->edge()->sV() == HE->edge()->eV());
727 std::sprintf (s, "%d ", HE->edge()->sV()->id()); ostrm<<s;
728 }
729 else {
730 bmsh3d_halfedge* HE = halfedge_;
731 do {
732 bmsh3d_halfedge* nextHE = HE->next();
733 bmsh3d_vertex* V = incident_V_of_Es (HE->edge(), nextHE->edge());
734
735 std::sprintf (s, "%d ", V->id()); ostrm<<s;
736 HE = HE->next();
737 }
738 while (HE != halfedge_);
739 }
740 }
741
742 //###############################################################
743 //###### For triangular face only ######
744 //###############################################################
745
tri_get_topo_type() const746 TRIFACE_TYPE bmsh3d_face::tri_get_topo_type () const
747 {
748 TRIFACE_TYPE type = BOGUS_TRIFACE;
749
750 bmsh3d_halfedge* HE = halfedge_;
751 bmsh3d_edge* e0 = HE->edge();
752 HE = HE->next();
753 bmsh3d_edge* e1 = HE->edge();
754 HE = HE->next();
755 bmsh3d_edge* e2 = HE->edge();
756
757 // Non-triangular face: more than 3 edges.
758 HE = HE->next();
759 if (HE->edge() != e0)
760 return TRIFACE_E4P;
761
762 unsigned int c0 = e0->n_incident_Fs();
763 unsigned int c1 = e1->n_incident_Fs();
764 unsigned int c2 = e2->n_incident_Fs();
765
766 if (c0 == 1) {
767 if (c1 == 1) {
768 if (c2 == 1)
769 return TRIFACE_111;
770 else if (c2 == 2)
771 return TRIFACE_112;
772 else { assert (c2 > 2);
773 return TRIFACE_113P; }
774 }
775 else if (c1 == 2) {
776 if (c2 == 1)
777 return TRIFACE_112;
778 else if (c2 == 2)
779 return TRIFACE_122;
780 else { assert (c2 > 2);
781 return TRIFACE_123P; }
782 }
783 else { assert (c1 > 2);
784 if (c2 == 1)
785 return TRIFACE_111;
786 else if (c2 == 2)
787 return TRIFACE_112;
788 else { assert (c2 > 2);
789 return TRIFACE_113P; }
790 }
791 }
792 else if (c0 == 2) {
793 if (c1 == 1) {
794 if (c2 == 1)
795 return TRIFACE_112;
796 else if (c2 == 2)
797 return TRIFACE_122;
798 else { assert (c2 > 2);
799 return TRIFACE_123P; }
800 }
801 else if (c1 == 2) {
802 if (c2 == 1)
803 return TRIFACE_122;
804 else if (c2 == 2)
805 return TRIFACE_222;
806 else { assert (c2 > 2);
807 return TRIFACE_223P; }
808 }
809 else { assert (c1 > 2);
810 if (c2 == 1)
811 return TRIFACE_123P;
812 else if (c2 == 2)
813 return TRIFACE_223P;
814 else { assert (c2 > 2);
815 return TRIFACE_23P3P; }
816 }
817 }
818 else { assert (c0 > 2);
819 if (c1 == 1) {
820 if (c2 == 1)
821 return TRIFACE_113P;
822 else if (c2 == 2)
823 return TRIFACE_123P;
824 else { assert (c2 > 2);
825 return TRIFACE_13P3P; }
826 }
827 else if (c1 == 2) {
828 if (c2 == 1)
829 return TRIFACE_123P;
830 else if (c2 == 2)
831 return TRIFACE_223P;
832 else { assert (c2 > 2);
833 return TRIFACE_23P3P; }
834 }
835 else { assert (c1 > 2);
836 if (c2 == 1)
837 return TRIFACE_13P3P;
838 else if (c2 == 2)
839 return TRIFACE_123P;
840 else { assert (c2 > 2);
841 return TRIFACE_3P3P3P; }
842 }
843 }
844
845 return type;
846 }
847
tri_get_topo_string() const848 std::string bmsh3d_face::tri_get_topo_string () const
849 {
850 TRIFACE_TYPE type = tri_get_topo_type();
851 switch (type) {
852 case BOGUS_TRIFACE: return "BOGUS_TRIFACE";
853 case TRIFACE_111: return "1_1_1";
854 case TRIFACE_112: return "1_1_2";
855 case TRIFACE_113P: return "1_1_3+";
856 case TRIFACE_122: return "1_2_2";
857 case TRIFACE_123P: return "1_2_3+";
858 case TRIFACE_13P3P: return "1_3+_3+";
859 case TRIFACE_222: return "2_2_2";
860 case TRIFACE_223P: return "2_2_3+";
861 case TRIFACE_23P3P: return "2_3+_3+";
862 case TRIFACE_3P3P3P: return "3+_3+_3+";
863 case TRIFACE_E4P: return "edges_4+";
864 default: assert(0); return "";
865 }
866 }
867
868 //###############################################################
869 //###### For the face of a 2-manifold triangular mesh only ######
870 //###############################################################
871
m2t_edge_against_vertex(bmsh3d_vertex * inputV)872 bmsh3d_edge* bmsh3d_face::m2t_edge_against_vertex (bmsh3d_vertex* inputV)
873 {
874 // loop through all incident edges, look for the one sharing inputV
875 bmsh3d_halfedge* HE = halfedge_;
876
877 // traverse through the circular list of halfedges,
878 do {
879 bmsh3d_edge* E = HE->edge();
880 if (!E->is_V_incident (inputV))
881 return E;
882
883 HE = HE->next();
884 }
885 while (HE != halfedge_);
886
887 return nullptr;
888 }
889
m2t_halfedge_against_vertex(bmsh3d_vertex * inputV)890 bmsh3d_halfedge* bmsh3d_face::m2t_halfedge_against_vertex (bmsh3d_vertex* inputV)
891 {
892 // loop through all incident edges, look for the one sharing inputV
893 bmsh3d_halfedge* HE = halfedge_;
894
895 // traverse through the circular list of halfedges,
896 do {
897 bmsh3d_edge* E = HE->edge();
898 if (!E->is_V_incident (inputV))
899 return HE;
900
901 HE = HE->next();
902 }
903 while (HE != halfedge_);
904
905 return nullptr;
906 }
907
908 //: for 2-manifold mesh, given input_face, find the neighboring face sharing the inputV.
909 // It is possible that nothing could be found, if the inputV is not correct.
m2t_nbr_face_against_vertex(bmsh3d_vertex * inputV)910 bmsh3d_face* bmsh3d_face::m2t_nbr_face_against_vertex (bmsh3d_vertex* inputV)
911 {
912 // loop through all incident edges, look for the one sharing inputV
913 bmsh3d_halfedge* HE = halfedge_;
914
915 // traverse through the circular list of halfedges,
916 do {
917 bmsh3d_edge* E = HE->edge();
918 if (!E->is_V_incident (inputV)) {
919 // found it, return the other face
920 // traverse through edge's list of halfedges.
921 // Only work for 2-manifold!!
922 if (E->halfedge()->face() != this)
923 return E->halfedge()->face();
924 else if (E->halfedge()->pair()) {
925 bmsh3d_halfedge* HE = E->halfedge()->pair();
926 assert (HE->face() != this);
927 return HE->face();
928 }
929 else {
930 return nullptr;
931 }
932 }
933 HE = HE->next();
934 }
935 while (HE != halfedge_);
936
937 return nullptr;
938 }
939
940 //: (redundant) for 2-manifold mesh, given input_face, find the neighboring face sharing incident v1 and v2
m2t_nbr_face_sharing_edge(bmsh3d_vertex * v1,bmsh3d_vertex * v2)941 bmsh3d_face* bmsh3d_face::m2t_nbr_face_sharing_edge (bmsh3d_vertex* v1, bmsh3d_vertex* v2)
942 {
943 // loop through all incident edges, look for the one sharing inputV
944 bmsh3d_halfedge* HE = halfedge_;
945
946 // traverse through the circular list of halfedges,
947 do {
948 bmsh3d_edge* E = HE->edge();
949 if (E->is_V_incident (v1) && E->is_V_incident (v2)) {
950 // found it, return the other face
951 // traverse through edge's list of halfedges.
952 // Only work for 2-manifold!!
953 if (E->halfedge()->face() != this)
954 return E->halfedge()->face();
955 else if (E->halfedge()->pair()) {
956 bmsh3d_halfedge* HE = E->halfedge()->pair();
957 return HE->face();
958 }
959 else {
960 return nullptr;
961 }
962 }
963 HE = HE->next();
964 }
965 while (HE != halfedge_);
966
967 return nullptr;
968 }
969
970 //: for triangular mesh only, given v1, v2, find v3
t_3rd_vertex(const bmsh3d_vertex * V1,const bmsh3d_vertex * V2) const971 bmsh3d_vertex* bmsh3d_face::t_3rd_vertex (const bmsh3d_vertex* V1, const bmsh3d_vertex* V2) const
972 {
973 // loop through all incident edges, look for the one sharing inputV
974 bmsh3d_halfedge* HE = halfedge_;
975
976 // traverse through the circular list of halfedges,
977 do {
978 bmsh3d_edge* E = HE->edge();
979 // some duplication, check both s_v and e_v.
980 if (E->sV() != V1 && E->sV() != V2)
981 return E->sV();
982 if (E->eV() != V1 && E->eV() != V2)
983 return E->eV();
984 HE = HE->next();
985 }
986 while (HE != halfedge_);
987
988 assert (0);
989 return nullptr;
990 }
991
t_vertex_against_edge(const bmsh3d_edge * E) const992 bmsh3d_vertex* bmsh3d_face::t_vertex_against_edge (const bmsh3d_edge* E) const
993 {
994 return t_3rd_vertex (E->sV(), E->eV());
995 }
996
997 //##########################################################
998 //###### Additional Functions ######
999 //##########################################################
1000
_find_prev_in_HE_chain(const bmsh3d_halfedge * inputHE)1001 bmsh3d_halfedge* _find_prev_in_HE_chain (const bmsh3d_halfedge* inputHE)
1002 {
1003 auto* HE = (bmsh3d_halfedge*) inputHE; // casting away const !!!
1004 while (HE->next() != inputHE)
1005 HE = HE->next();
1006 return HE;
1007 }
1008
1009 //: disconnect all associated halfedges from their edges from the given headHE.
_delete_HE_chain(bmsh3d_halfedge * & headHE)1010 void _delete_HE_chain (bmsh3d_halfedge* & headHE)
1011 {
1012 if (headHE == nullptr)
1013 return;
1014
1015 bmsh3d_halfedge* nextHE = headHE->next();
1016 if (nextHE) {
1017 do {
1018 bmsh3d_halfedge* he_to_remove = nextHE;
1019 nextHE = nextHE->next();
1020 bmsh3d_edge* E = he_to_remove->edge();
1021 E->_disconnect_HE (he_to_remove);
1022 delete he_to_remove;
1023 }
1024 while (nextHE != headHE);
1025 }
1026
1027 // Delete the headHE and set to NULL.
1028 bmsh3d_edge* E = headHE->edge();
1029 E->_disconnect_HE (headHE);
1030 delete headHE;
1031 headHE = nullptr;
1032 }
1033
1034 // Return: the set of incident edges that get disconnected.
1035 // Also set the headHE to be NULL after calling it.
_delete_HE_chain(bmsh3d_halfedge * & headHE,std::vector<bmsh3d_edge * > & incident_edge_list)1036 void _delete_HE_chain (bmsh3d_halfedge* & headHE,
1037 std::vector<bmsh3d_edge*>& incident_edge_list)
1038 {
1039 if (headHE == nullptr)
1040 return;
1041
1042 bmsh3d_halfedge* nextHE = headHE->next();
1043 if (nextHE) {
1044 do {
1045 bmsh3d_halfedge* he_to_remove = nextHE;
1046 nextHE = nextHE->next();
1047 bmsh3d_edge* E = he_to_remove->edge();
1048 E->_disconnect_HE (he_to_remove);
1049 delete he_to_remove;
1050 incident_edge_list.push_back (E);
1051 }
1052 while (nextHE != headHE);
1053 }
1054
1055 // Delete the headHE and set to NULL.
1056 bmsh3d_edge* E = headHE->edge();
1057 E->_disconnect_HE (headHE);
1058 delete headHE;
1059 headHE = nullptr;
1060 incident_edge_list.push_back (E);
1061 }
1062
1063 //: Given the face, current halfedge, and current eV, find the next halfedge in the input storage.
_find_next_halfedge(bmsh3d_halfedge * input_he,bmsh3d_vertex * eV,std::vector<bmsh3d_halfedge * > & incident_HEs)1064 bmsh3d_halfedge* _find_next_halfedge(bmsh3d_halfedge* input_he,
1065 bmsh3d_vertex* eV,
1066 std::vector<bmsh3d_halfedge*>& incident_HEs)
1067 {
1068 // Search for the next halfedge that's not the input_he
1069 for (auto HE : incident_HEs) {
1070 bmsh3d_edge* E = HE->edge();
1071 // skip if HE is the input_he
1072 if (HE == input_he)
1073 continue;
1074 if (E->is_V_incident (eV))
1075 return HE;
1076 }
1077 return nullptr;
1078 }
1079
1080 //: Assume the mesh face is planar and compute a 2D planar coordinate for it.
get_2d_coord(const std::vector<bmsh3d_vertex * > & vertices,vgl_vector_3d<double> & N,vgl_vector_3d<double> & AX,vgl_vector_3d<double> & AY)1081 void get_2d_coord (const std::vector<bmsh3d_vertex*>& vertices, vgl_vector_3d<double>& N,
1082 vgl_vector_3d<double>& AX, vgl_vector_3d<double>& AY)
1083 {
1084 assert (vertices.size() > 2);
1085 // The first vertex A is the origin.
1086 bmsh3d_vertex* A = vertices[0];
1087
1088 // The second vertex B identifies the x-axis.
1089 bmsh3d_vertex* B = vertices[1];
1090 // Make AB the unit vector in the 2D x-direction.
1091 AX = B->pt() - A->pt();
1092 AX = AX / AX.length();
1093
1094 // Use The 3rd vertex C to identify normal N.
1095 bmsh3d_vertex* C = vertices[2];
1096 vgl_vector_3d<double> AC = C->pt() - A->pt();
1097 N = cross_product (AX, AC);
1098 N = N / N.length();
1099
1100 // The unit vector in the 2D y-direction.
1101 AY = cross_product (N, AX);
1102 AY = AY / AY.length();
1103 }
1104
1105 //: Return ordered set of vertices in 2D (x,y) coord.
1106 // Assume mesh face is planar. First identify the plane normal N = AB * AC.
1107 // Use the first vertex A as origin (0,0),
1108 // the second vertex B to identify x-axis, B(d, 0).
1109 // the y-axis is AY = cross (N, AB).
1110 //
get_2d_polygon(const std::vector<bmsh3d_vertex * > & vertices,std::vector<double> & xs,std::vector<double> & ys)1111 void get_2d_polygon (const std::vector<bmsh3d_vertex*>& vertices,
1112 std::vector<double>& xs, std::vector<double>& ys)
1113 {
1114 assert (vertices.size() > 2);
1115 xs.resize (vertices.size());
1116 ys.resize (vertices.size());
1117
1118 // The first vertex A is the origin.
1119 bmsh3d_vertex* A = vertices[0];
1120 xs[0] = 0;
1121 ys[0] = 0;
1122
1123 // The second vertex B identifies the x-axis.
1124 bmsh3d_vertex* B = vertices[1];
1125 double dAB = vgl_distance (A->pt(), B->pt());
1126 vgl_vector_3d<double> AB = B->pt() - A->pt();
1127 // Make AB the unit vector in the 2D x-direction.
1128 AB = AB / dAB;
1129 xs[1] = dAB;
1130 ys[1] = 0;
1131
1132 // Use The 3rd vertex C to identify normal N.
1133 bmsh3d_vertex* C = vertices[2];
1134 vgl_vector_3d<double> AC = C->pt() - A->pt();
1135 vgl_vector_3d<double> N = cross_product (AB, AC);
1136 // Check that N is valid.
1137 // The unit vector in the 2D y-direction.
1138 vgl_vector_3d<double> AY = cross_product (N, AB);
1139 AY = AY / AY.length();
1140
1141 // Loop through all other vertices.
1142 for (unsigned int i=2; i<vertices.size(); i++) {
1143 bmsh3d_vertex* C = vertices[i];
1144 vgl_vector_3d<double> AC = C->pt() - A->pt();
1145 #if 0
1146 cx = d * cos theta = dot(AB, AC) / AB.
1147 cy = std::sqrt(AC^2 - cx^2)
1148 #endif // 0
1149 xs[i] = dot_product(AC, AB);
1150 ys[i] = dot_product(AC, AY);
1151 }
1152 }
1153
1154 //: Return the projected point in the local 2D (x,y) coord.
get_2d_proj_pt(vgl_point_3d<double> P,const vgl_point_3d<double> & A,const vgl_vector_3d<double> & AX,const vgl_vector_3d<double> & AY)1155 vgl_point_2d<double> get_2d_proj_pt (vgl_point_3d<double> P, const vgl_point_3d<double>& A,
1156 const vgl_vector_3d<double>& AX,
1157 const vgl_vector_3d<double>& AY)
1158 {
1159 double x = dot_product (P - A, AX);
1160 double y = dot_product (P - A, AY);
1161 return {x, y};
1162 }
1163
1164 //: determine the center point of the patch
compute_cen(const std::vector<bmsh3d_vertex * > & vertices)1165 vgl_point_3d<double> compute_cen (const std::vector<bmsh3d_vertex*>& vertices)
1166 {
1167 double x=0.0, y=0.0, z=0.0;
1168 assert (vertices.size() != 0);
1169
1170 for (auto vertice : vertices) {
1171 x += vertice->pt().x();
1172 y += vertice->pt().y();
1173 z += vertice->pt().z();
1174 }
1175
1176 x /= vertices.size();
1177 y /= vertices.size();
1178 z /= vertices.size();
1179
1180 return {x,y,z};
1181 }
1182
1183 //: compute the normal of the face
compute_normal_ifs(const std::vector<bmsh3d_vertex * > & vertices)1184 vgl_vector_3d<double> compute_normal_ifs (const std::vector<bmsh3d_vertex*>& vertices)
1185 {
1186 vgl_vector_3d<double> normal;
1187 // for P[0..n], compute cross product of (P[1]-P[0])*(P[2]-P[1]) ...
1188 assert (vertices.size() > 2);
1189 for (int i=0; i< ((int)vertices.size())-2; i++) {
1190 const bmsh3d_vertex* v0 = vertices[i];
1191 const bmsh3d_vertex* v1 = vertices[i+1];
1192 const bmsh3d_vertex* v2 = vertices[i+2];
1193 vgl_vector_3d<double> a = v1->pt() - v0->pt();
1194 vgl_vector_3d<double> b = v2->pt() - v1->pt();
1195 vgl_vector_3d<double> n = cross_product (a, b);
1196 normal += n;
1197 }
1198 return normal;
1199 }
1200
1201 //: Compute face normal using the given edge and starting node.
compute_normal(const vgl_point_3d<double> & C,const bmsh3d_edge * E,const bmsh3d_vertex * startV)1202 vgl_vector_3d<double> compute_normal (const vgl_point_3d<double>& C,
1203 const bmsh3d_edge* E,
1204 const bmsh3d_vertex* startV)
1205 {
1206 const bmsh3d_vertex* Es = startV;
1207 const bmsh3d_vertex* Ee = E->other_V (Es);
1208 return cross_product (Es->pt() - C, Ee->pt() - C);
1209 }
1210
1211 //: Return true if vertices represent a polygon (Vs > 3) or obtuse triangle.
is_tri_non_acute(const std::vector<bmsh3d_vertex * > & vertices)1212 bool is_tri_non_acute (const std::vector<bmsh3d_vertex*>& vertices)
1213 {
1214 if (vertices.size() > 3)
1215 return true;
1216 bmsh3d_vertex* V0 = vertices[0];
1217 bmsh3d_vertex* V1 = vertices[1];
1218 bmsh3d_vertex* V2 = vertices[2];
1219 return bmsh3d_is_tri_non_acute (V0->pt(), V1->pt(), V2->pt());
1220 }
1221
1222 //: Return true if F (triangle) is (1,1,3+) or (1,3+,3+) extraneous
1223 // For a general polygon, return true if F has only edges of 1 or 3 incidences.
1224 //
is_F_extraneous(bmsh3d_face * F)1225 bool is_F_extraneous (bmsh3d_face* F)
1226 {
1227 assert (F->halfedge() != nullptr);
1228 assert (F->halfedge()->next() != nullptr);
1229 bmsh3d_halfedge* HE = F->halfedge();
1230 do {
1231 bmsh3d_edge* E = HE->edge();
1232 int n = E->n_incident_Fs();
1233 if (n==2)
1234 return false;
1235 HE = HE->next();
1236 }
1237 while (HE != F->halfedge());
1238 return true;
1239 }
1240
get_F_sharing_Es(bmsh3d_edge * E1,bmsh3d_edge * E2)1241 bmsh3d_face* get_F_sharing_Es (bmsh3d_edge* E1, bmsh3d_edge* E2)
1242 {
1243 // Loop through all incident faces of E1 and find the face incident to E2.
1244 if (E1->halfedge() == nullptr) {
1245 return nullptr;
1246 }
1247 else if (E1->halfedge()->pair() == nullptr) {
1248 bmsh3d_face* F = E1->halfedge()->face();
1249 if (F->is_E_incident (E2))
1250 return F;
1251 }
1252 else {
1253 bmsh3d_halfedge* HE = E1->halfedge();
1254 do {
1255 bmsh3d_face* F = HE->face();
1256 if (F->is_E_incident (E2))
1257 return F;
1258 HE = HE->pair();
1259 }
1260 while (HE != E1->halfedge());
1261 }
1262 return nullptr;
1263 }
1264
1265 //###############################################################
1266 #if 0 // Old code from Frederic
1267 //###############################################################
1268
1269 //: compute the normalized N = V1 * V2
1270 void _cross_product (double V1x, double V1y, double V1z,
1271 double V2x, double V2y, double V2z,
1272 double& Nx, double& Ny, double& Nz)
1273 {
1274 // Vector product : NormTri = V1 x V2
1275 Nx = (V1y * V2z) - (V1z * V2y);
1276 Ny = (V1z * V2x) - (V1x * V2z);
1277 Nz = (V1x * V2y) - (V1y * V2x);
1278 }
1279
1280 #define VECTOR_LENGTH_EPSILON 1E-14
1281
1282 void _normalize_vector (double& Vx, double& Vy, double& Vz)
1283 {
1284 // Normalize the normal vector
1285 double dLength = std::sqrt (Vx*Vx + Vy*Vy + Vz*Vz);
1286 if (dLength < VECTOR_LENGTH_EPSILON) {
1287 vul_printf (std::cerr, "NUMERICAL ERROR in computing vector length.\n");
1288 }
1289 Vx = Vx / dLength;
1290 Vy = Vy / dLength;
1291 Vz = Vz / dLength;
1292 }
1293
1294 void _cross_product_normalized (double V1x, double V1y, double V1z,
1295 double V2x, double V2y, double V2z,
1296 double& Nx, double& Ny, double& Nz)
1297 {
1298 _cross_product (V1x, V1y, V1z, V2x, V2y, V2z, Nx, Ny, Nz);
1299 _normalize_vector (Nx, Ny, Nz);
1300 }
1301
1302 //--------------------------------------------------------------------
1303
1304 //:
1305 // Given an ordered quadruplet (A,B,C,D) with C, the "left" gene,
1306 // A and B, the edge pair, and D the next one the "right", find
1307 // the dihedral angle between the pair of triangular faces
1308 // (ABC) and (BCD), N1 = BC x BA, N2 = BD x BC ,
1309 // that is the angle the two faces make with respect to the
1310 // shared edge BC.
1311 //
1312 // GeneA: on the linkElm triangle
1313 // edgePatchElm: GeneB, GeneC
1314 // GeneD: the Gene to test
1315 //
1316 // \returns the angle
1317 double computeDotNormal (const bmsh3d_vertex* GeneC,
1318 const dbsk3d_fs_patch_elm* edgePatchElm,
1319 const bmsh3d_vertex* GeneD)
1320 {
1321 const bmsh3d_vertex* GeneA = edgePatchElm->genes(0);
1322 const bmsh3d_vertex* GeneB = edgePatchElm->genes(1);
1323
1324 // --- Compute oriented normal to 1st triangle: N1 = AB x AC ---
1325 double ABx = GeneB->pt().x() - GeneA->pt().x();
1326 double ABy = GeneB->pt().y() - GeneA->pt().y();
1327 double ABz = GeneB->pt().z() - GeneA->pt().z();
1328
1329 double ACx = GeneC->pt().x() - GeneA->pt().x();
1330 double ACy = GeneC->pt().y() - GeneA->pt().y();
1331 double ACz = GeneC->pt().z() - GeneA->pt().z();
1332
1333 // Vector product : BC x BA = N1
1334 double N1x, N1y, N1z;
1335 _cross_product_normalized (ABx, ABy, ABz, ACx, ACy, ACz, N1x, N1y, N1z);
1336
1337 // --- Compute oriented normal to 2nd triangle: N2 = BD x BC ---
1338 double fBDx = GeneD->pt().x() - GeneA->pt().x();
1339 double fBDy = GeneD->pt().y() - GeneA->pt().y();
1340 double fBDz = GeneD->pt().z() - GeneA->pt().z();
1341
1342 // Vector product : BD x BC = N2
1343 double N2x, N2y, N2z;
1344 _cross_product_normalized(fBDx, fBDy, fBDz, ABx, ABy, ABz, N2x, N2y, N2z);
1345
1346 // scalar product of N1 and N2
1347 double dDotProduct = (N1x * N2x) + (N1y * N2y) + (N1z * N2z);
1348
1349 // dot product of unit vector can not be greater than 1
1350 if (dDotProduct > 1.0)
1351 dDotProduct = 1.0;
1352 else if (dDotProduct < -1.0)
1353 dDotProduct = -1.0;
1354
1355 return dDotProduct;
1356 }
1357
1358 #endif // 0
1359