libsurf
Programmer's Documentation

dntriangle.h (r6227/r5385)
1
2 /* Copyright (C) 2015 David Eller <david@larosterna.com>
3  *
5  * Licensees holding valid commercial licenses may use this file in accordance
6  * with the terms contained in their respective non-exclusive license agreement.
7  * For further information contact david@larosterna.com .
8  *
9  * GNU General Public License Usage
10  * Alternatively, this file may be used under the terms of the GNU General
12  * appearing in the file gpl.txt included in the packaging of this file.
13  */
14
15 #ifndef SURF_DNTRIANGLE_H
16 #define SURF_DNTRIANGLE_H
17
18 #include <surf/surface.h>
19 #include <predicates/predicates.h>
20 #include "dnvertex.h"
21 #include "dnedge.h"
22
32 {
33  public:
34
36  DnTriangle(uint a, uint b, uint c) {
37  assert(a != NotFound);
38  assert(b != NotFound);
39  assert(c != NotFound);
40  assert(a != b);
41  assert(a != c);
42  assert(b != c);
43  order(a, b, c);
44  nbe[2] = nbe[1] = nbe[0] = NotFound;
45  }
46
48  void reconnect(uint a, uint b, uint c) {
49  assert(a != NotFound);
50  assert(b != NotFound);
51  assert(c != NotFound);
52  assert(a != b);
53  assert(a != c);
54  assert(b != c);
55  order(a, b, c);
56  nbe[2] = nbe[1] = nbe[0] = NotFound;
57  }
58
60  void computeSphere(const Surface &, const DnVertexArray & vtx, bool spt) {
61  assert(isValid());
62  const Vct3 & p1(vtx[vi[0]].eval());
63  const Vct3 & p2(vtx[vi[1]].eval());
64  const Vct3 & p3(vtx[vi[2]].eval());
65  pcs = sCircumCenter(vtx);
66  if (spt) {
67  Real r = norm(pcs - p1);
68  Vct3 tn( cross(p2-p1, p3-p1) );
69  normalize(tn);
70  pcs -= r*tn;
71  }
72  }
73
75  bool isValid() const {return (vi[0] != NotFound);}
76
78  bool hasDuplicates() const {
79  if (vi[0] == vi[1])
80  return true;
81  else if (vi[0] == vi[2])
82  return true;
83  else if (vi[1] == vi[2])
84  return true;
85  else
86  return false;
87  }
88
90  void invalidate() {memset(vi, (int) NotFound, 3*sizeof(uint));}
91
93  const uint *vertices() const {return vi;}
94
96  const uint *nbEdges() const {return nbe;}
97
99  uint *nbEdges() {return nbe;}
100
102  void reverse() {
103  assert(isValid());
104  std::swap(vi[1], vi[2]);
105  }
106
108  void itranslate(const Indices & repl) {
109  order(repl[vi[0]], repl[vi[1]], repl[vi[2]]);
110  }
111
113  uint opposedVertex(const DnEdge & e) const {
114  const uint s(e.source());
115  const uint t(e.target());
116  for (uint k=0; k<3; ++k) {
117  uint v = vi[k];
118  if (v != s and v != t)
119  return v;
120  }
121  return NotFound;
122  }
123
125  uint replaceVertex(uint vold, uint vnew) {
126  for (uint k=0; k<3; ++k) {
127  if (vi[k] == vold) {
128  vi[k] = vnew;
129  order(vi[0], vi[1], vi[2]);
130  return k;
131  }
132  }
133  return NotFound;
134  }
135
137  uint attachEdge(uint ei) {
138  if (ei == NotFound)
139  return NotFound;
140  else if (nbe[0] == ei)
141  return 0;
142  else if (nbe[1] == ei)
143  return 1;
144  else if (nbe[2] == ei)
145  return 2;
146  else if (nbe[0] == NotFound) {
147  nbe[0] = ei;
148  return 0;
149  } else if (nbe[1] == NotFound) {
150  nbe[1] = ei;
151  return 1;
152  } else if (nbe[2] == NotFound) {
153  nbe[2] = ei;
154  return 2;
155  } else
156  return NotFound;
157  }
158
160  uint detachEdge(uint ei) {
161  if (ei == NotFound)
162  return NotFound;
163  else if (nbe[0] == ei) {
164  nbe[0] = NotFound;
165  return 0;
166  } else if (nbe[1] == ei) {
167  nbe[1] = NotFound;
168  return 1;
169  } else if (nbe[2] == ei) {
170  nbe[2] = NotFound;
171  return 2;
172  } else
173  return NotFound;
174  }
175
177  uint replaceEdge(uint fold, uint fnew) {
178  if (nbe[0] == fold) {
179  nbe[0] = fnew;
180  return 0;
181  } else if (nbe[1] == fold) {
182  nbe[1] = fnew;
183  return 1;
184  } else if (nbe[2] == fold) {
185  nbe[2] = fnew;
186  return 2;
187  } else
188  return NotFound;
189  }
190
192  Vct3 normal(const DnVertexArray & vtx) const {
193  assert(isValid());
194  const Vct3 & p1(vtx[vi[0]].eval());
195  const Vct3 & p2(vtx[vi[1]].eval());
196  const Vct3 & p3(vtx[vi[2]].eval());
197  return cross(p2-p1, p3-p1);
198  }
199
201  Vct2 pCenter(const DnVertexArray & vtx) const {
202  assert(isValid());
203  const Vct2 & q1(vtx[vi[0]].parpos());
204  const Vct2 & q2(vtx[vi[1]].parpos());
205  const Vct2 & q3(vtx[vi[2]].parpos());
206  return (q1+q2+q3) / 3.0;
207  }
208
210  Vct3 sCenter(const DnVertexArray & vtx) const {
211  const Vct3 & p1(vtx[vi[0]].eval());
212  const Vct3 & p2(vtx[vi[1]].eval());
213  const Vct3 & p3(vtx[vi[2]].eval());
214  return (p1+p2+p3) / 3.0;
215  }
216
218  Vct2 pCircumCenter(const DnVertexArray & vtx) const {
219  assert(isValid());
220  const Vct2 & q1(vtx[vi[0]].parpos());
221  const Vct2 & q2(vtx[vi[1]].parpos());
222  const Vct2 & q3(vtx[vi[2]].parpos());
223
224  // compute edge-normal directions
225  Vct2 e1( q2-q1 );
226  Vct2 e2( q3-q2 );
227  Vct2 n1, n2;
228  if (e1[1] != 0.0) {
229  n1[0] = 1.0;
230  n1[1] = -e1[0]/e1[1];
231  } else {
232  assert(e1[0] != 0.0);
233  n1[0] = -e1[1]/e1[0];
234  n1[1] = 1.0;
235  }
236  if (e2[1] != 0.0) {
237  n2[0] = 1.0;
238  n2[1] = -e2[0]/e2[1];
239  } else {
240  assert(e2[0] != 0.0);
241  n2[0] = -e2[1]/e2[0];
242  n2[1] = 1.0;
243  }
244
245  // intersect lines through edge midpoints
246  Vct2 m1( 0.5*(q1+q2) );
247  Vct2 m2( 0.5*(q2+q3) );
248
249  Real a11 = n1[0];
250  Real a12 = n1[1];
251  Real a21 = -n2[0];
252  Real a22 = -n2[1];
253  Real det = a11*a22 - a12*a21;
254  assert(det != 0.0);
255  Real r1 = m2[0] - m1[0];
256  Real r2 = m2[1] - m1[1];
257
258  // line through m1 is m1 + s*n1
259  Real s = (r1*a22 - r2*a12) / det;
260
261  return m1 + s*n1;
262  }
263
265  Vct3 sCircumCenter(const DnVertexArray & vtx) const {
266  assert(isValid());
267  const Vct3 & p1(vtx[vi[0]].eval());
268  const Vct3 & p2(vtx[vi[1]].eval());
269  const Vct3 & p3(vtx[vi[2]].eval());
270
271  // triangle normal and edge-normal directions
272  Vct3 tn( cross(p2-p1,p3-p1) );
273  Vct3 e1( cross(tn, p2-p1) );
274  Vct3 e2( cross(tn, p3-p2) );
275
276  // edge midpoints
277  Vct3 m1( 0.5*(p1+p2) );
278  Vct3 m2( 0.5*(p2+p3) );
279
280  Real a11(0.0), a12(0.0), a21(0.0), a22(0.0);
281  Real r1(0.0), r2(0.0);
282  for (uint i=0; i<3; i++) {
283  a11 += sq(e1[i]);
284  a12 -= e1[i]*e2[i];
285  a22 += sq(e2[i]);
286  r1 -= (m1[i] - m2[i])*e1[i];
287  r2 += (m1[i] - m2[i])*e2[i];
288  }
289  a21 = a12;
290
291  // zero means lines are parallel
292  Real det = a11*a22 - a12*a21;
293  if (det == 0.0) {
294  return m1 + huge*e1;
295  }
296
297  Real s = (r1*a22 - r2*a12) / det;
298  return m1 + s*e1;
299  }
300
304  int isInside(const DnEdgeArray & edges, const DnVertexArray & vtx,
305  const Vct2 & p) const
306  {
307  assert(isValid());
308  double ort;
309  uint v1, v2;
310  for (uint k=0; k<3; ++k) {
311  assert(nbe[k] != NotFound);
312  const DnEdge & e(edges[nbe[k]]);
313  v1 = find(e.source());
314  v2 = find(e.target());
315  assert(v1 != NotFound and v2 != NotFound);
316  if (v1 == 1 and v2 == 0)
317  std::swap(v1,v2);
318  else if (v1 == 2 and v2 == 1)
319  std::swap(v1,v2);
320  else if (v1 == 0 and v2 == 2)
321  std::swap(v1,v2);
322  const Vct2 & q1(vtx[vi[v1]].parpos());
323  const Vct2 & q2(vtx[vi[v2]].parpos());
324  ort = jrsOrient2d(q1.pointer(), q2.pointer(), p.pointer());
325  if (ort < 0)
326  return -2;
327  else if (ort == 0) {
328  Vct2 tmp = q2-q1;
329  Real t = dot(p-q1, tmp) / dot(tmp,tmp);
330  if (t >= 0.0 and t <= 1.0)
331  return k;
332  else
333  return -2;
334  }
335  }
336  return -1;
337  }
338
342  int inCircle(const DnVertexArray & vtx, uint i) const {
343  const Vct2 & q1(vtx[vi[0]].parpos());
344  const Vct2 & q2(vtx[vi[1]].parpos());
345  const Vct2 & q3(vtx[vi[2]].parpos());
346  const Vct2 & pt(vtx[i].parpos());
347  double ict = jrsInCircle(q1.pointer(), q2.pointer(), q3.pointer(), pt.pointer());
348  if (ict > 0.0)
349  return 1;
350  else if (ict < 0.0)
351  return -1;
352  else
353  return 0;
354  }
355
357  int inSphere(const DnVertexArray & vtx, uint i) const {
358  assert(isValid());
359  const Vct3 & p1(vtx[vi[0]].eval());
360  const Vct3 & p2(vtx[vi[1]].eval());
361  const Vct3 & p3(vtx[vi[2]].eval());
362  const Vct3 & pt(vtx[i].eval());
363  double inside = jrsInSphere(p1.pointer(), p2.pointer(), p3.pointer(),
364  pcs.pointer(), pt.pointer());
365  if (inside > 0.0)
366  return 1;
367  else if (inside < 0.0)
368  return -1;
369  else
370  return 0;
371  }
372
374  Vct3 project(const DnVertexArray & vtx, const Vct3 & pt) const {
375  const Vct3 & p1( vtx[vi[0]].eval() );
376  const Vct3 & p2( vtx[vi[1]].eval() );
377  const Vct3 & p3( vtx[vi[2]].eval() );
378
379  Vct3 va, vb, vXi, vEta;
380  va = p2 - p1;
381  vb = p3 - p1;
382  vXi = va - vb*(dot(va,vb)/dot(vb,vb));
383  vEta = vb - va*(dot(va,vb)/dot(va,va));
384
385  Vct3 s;
386  s[0] = dot(pt-p1, vXi) / dot(vXi,vXi);
387  s[1] = dot(pt-p1, vEta) / dot(vEta,vEta);
388  s[2] = 1.0 - s[0] - s[1];
389  return s;
390  }
391
393  Vct2 sProject(const DnVertexArray & vtx, const Vct3 & pt) const {
394  const Vct3 & p1( vtx[vi[0]].eval() );
395  const Vct3 & p2( vtx[vi[1]].eval() );
396  const Vct3 & p3( vtx[vi[2]].eval() );
397
398  Vct3 va, vb, vXi, vEta;
399  va = p2 - p1;
400  vb = p3 - p1;
401  vXi = va - vb*(dot(va,vb)/dot(vb,vb));
402  vEta = vb - va*(dot(va,vb)/dot(va,va));
403
404  Real up, vp, wp;
405  up = dot(pt-p1, vXi) / dot(vXi,vXi);
406  vp = dot(pt-p1, vEta) / dot(vEta,vEta);
407  wp = 1.0 - up - vp;
408
409  const Vct2 & q1(vtx[vi[0]].parpos());
410  const Vct2 & q2(vtx[vi[1]].parpos());
411  const Vct2 & q3(vtx[vi[2]].parpos());
412  Vct2 qp( wp*q1 + up*q2 + vp*q3 );
413  qp[0] = std::min(1.0, std::max(0.0, qp[0]));
414  qp[1] = std::min(1.0, std::max(0.0, qp[1]));
415  return qp;
416  }
417
419  void pFixDirection(const DnVertexArray & vtx) {
420  const Vct2 & q1(vtx[vi[0]].parpos());
421  const Vct2 & q2(vtx[vi[1]].parpos());
422  const Vct2 & q3(vtx[vi[2]].parpos());
423  assert(norm(q1-q2) > 0);
424  assert(norm(q1-q3) > 0);
425  assert(norm(q3-q2) > 0);
426  double ot = jrsOrient2d(q1.pointer(), q2.pointer(), q3.pointer());
427  if (ot < 0)
428  reverse();
429  assert(ot != 0);
430  }
431
433  void sFixDirection(const DnVertexArray & vtx) {
434  const Vct3 & n1( vtx[vi[0]].normal() );
435  const Vct3 & n2( vtx[vi[1]].normal() );
436  const Vct3 & n3( vtx[vi[2]].normal() );
437  Vct3 sn( n1+n2+n3 );
438
439  const Vct3 & p1( vtx[vi[0]].eval() );
440  const Vct3 & p2( vtx[vi[1]].eval() );
441  const Vct3 & p3( vtx[vi[2]].eval() );
442  Vct3 tn( cross(p2-p1, p3-p1) );
443
444  Real s = dot(sn, tn);
445  if (s < 0)
446  reverse();
447  // assert(s != 0);
448  }
449
451  uint find(uint v) const {
452  for (uint k=0; k<3; ++k)
453  if (vi[k] == v)
454  return k;
455  return NotFound;
456  }
457
459  uint findEdge(uint e) const {
460  for (uint k=0; k<3; ++k)
461  if (nbe[k] == e)
462  return k;
463  return NotFound;
464  }
465
466  private:
467
469  void order(uint a, uint b, uint c) {
470  if (a < b and a < c) {
471  vi[0] = a;
472  vi[1] = b;
473  vi[2] = c;
474  } else if (b < a and b < c) {
475  vi[0] = b;
476  vi[1] = c;
477  vi[2] = a;
478  } else {
479  vi[0] = c;
480  vi[1] = a;
481  vi[2] = b;
482  }
483  }
484
486  bool touchesUBound(const DnVertexArray & vtx) const {
487  assert(isValid());
488  for (uint k=0; k<3; ++k) {
489  const Vct2 & q(vtx[vi[k]].parpos());
490  if (q[0] == 0.0 or q[0] == 1.0)
491  return true;
492  }
493  return false;
494  }
495
496  private:
497
499  Vct3 pcs;
500
502  uint vi[3];
503
505  uint nbe[3];
506 };
507
508 #endif
uint find(uint v) const
find index of vertex index v
Definition: dntriangle.h:451
int inSphere(const DnVertexArray &vtx, uint i) const
check if point is inside circumsphere
Definition: dntriangle.h:357
uint opposedVertex(const DnEdge &e) const
find vertex opposed to edge e
Definition: dntriangle.h:113
void invalidate()
mark triangle as invalid
Definition: dntriangle.h:90
uint attachEdge(uint ei)
Definition: dntriangle.h:137
void pFixDirection(const DnVertexArray &vtx)
make sure that the normal direction is correct (u,v space)
Definition: dntriangle.h:419
const uint * nbEdges() const
access neighbor edges
Definition: dntriangle.h:96
Vct3 sCenter(const DnVertexArray &vtx) const
center in real space
Definition: dntriangle.h:210
bool isValid() const
check if triangle is defined
Definition: dntriangle.h:75
uint * nbEdges()
access neighbor edges
Definition: dntriangle.h:99
Vct3 project(const DnVertexArray &vtx, const Vct3 &pt) const
compute projection (ut,vt,wt) of a 3D point on triangle
Definition: dntriangle.h:374
Vct2 pCenter(const DnVertexArray &vtx) const
center in parametric space
Definition: dntriangle.h:201
Edge in a Delaunay triangulation.
Definition: dnedge.h:28
bool hasDuplicates() const
check if triangle has duplicate vertices
Definition: dntriangle.h:78
Triangle in 3D Delaunay triangulation.
Definition: dntriangle.h:31
const uint * vertices() const
access vertices
Definition: dntriangle.h:93
uint source() const
access source vertex index
Definition: dnedge.h:58
int inCircle(const DnVertexArray &vtx, uint i) const
Check if point is in circumcircle (uv-plane).
Definition: dntriangle.h:342
DnTriangle(uint a, uint b, uint c)
create new triangle
Definition: dntriangle.h:36
void itranslate(const Indices &repl)
translate vertex indices
Definition: dntriangle.h:108
uint detachEdge(uint ei)
remove edge from neighbor list
Definition: dntriangle.h:160
uint replaceVertex(uint vold, uint vnew)
replace a single vertex index
Definition: dntriangle.h:125
bool touchesUBound(const DnVertexArray &vtx) const
decide if triangle touches a u-boundary with at vertex
Definition: dntriangle.h:486
Vct2 sProject(const DnVertexArray &vtx, const Vct3 &pt) const
compute projection (u,v) of a 3D point on mesh surface
Definition: dntriangle.h:393
int isInside(const DnEdgeArray &edges, const DnVertexArray &vtx, const Vct2 &p) const
Check if point is inside triangle.
Definition: dntriangle.h:304
Vct3 sCircumCenter(const DnVertexArray &vtx) const
compute spatial circumcenter
Definition: dntriangle.h:265
uint vi[3]
vertex indices
Definition: dntriangle.h:502
void reverse()
reverse normal vector
Definition: dntriangle.h:102
void order(uint a, uint b, uint c)
establish correct ordering
Definition: dntriangle.h:469
uint target() const
access target vertex index
Definition: dnedge.h:61
void reconnect(uint a, uint b, uint c)
change triangle vertices
Definition: dntriangle.h:48
Surface interface.
Definition: surface.h:37
void sFixDirection(const DnVertexArray &vtx)
make sure that the normal direction is correct (3-space)
Definition: dntriangle.h:433
Vct3 normal(const DnVertexArray &vtx) const
compute normal vector (not normalized)
Definition: dntriangle.h:192
uint findEdge(uint e) const
find index of edge index e
Definition: dntriangle.h:459
uint nbe[3]
neighbor edges
Definition: dntriangle.h:505
Vct2 pCircumCenter(const DnVertexArray &vtx) const
compute parametric circumcenter
Definition: dntriangle.h:218
uint replaceEdge(uint fold, uint fnew)
replace neighbor edge
Definition: dntriangle.h:177
Vct3 pcs
point on circumsphere
Definition: dntriangle.h:499
void computeSphere(const Surface &, const DnVertexArray &vtx, bool spt)
compute circumsphere
Definition: dntriangle.h:60
Generated on Mon Jan 24 2022 03:03:16 for libsurf by   1.8.5