gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
curvature.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <stdint.h>
9 #include "curvature.h"
10 
11 #define tolerance 0.1e-20
12 
13 static bool Inv4x4ColumnMajor(double m[16], double inv[16], double *det)
14 {
15  int i;
16 
17  inv[0] = m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] +
18  m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10];
19  inv[4] = -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] -
20  m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10];
21  inv[8] = m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] +
22  m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[9];
23  inv[12] = -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] -
24  m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[9];
25  inv[1] = -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] -
26  m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10];
27  inv[5] = m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] +
28  m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10];
29  inv[9] = -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] -
30  m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[9];
31  inv[13] = m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] +
32  m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[9];
33  inv[2] = m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] +
34  m[5] * m[3] * m[14] + m[13] * m[2] * m[7] - m[13] * m[3] * m[6];
35  inv[6] = -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] -
36  m[4] * m[3] * m[14] - m[12] * m[2] * m[7] + m[12] * m[3] * m[6];
37  inv[10] = m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] +
38  m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5];
39  inv[14] = -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] -
40  m[4] * m[2] * m[13] - m[12] * m[1] * m[6] + m[12] * m[2] * m[5];
41  inv[3] = -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -
42  m[5] * m[3] * m[10] - m[9] * m[2] * m[7] + m[9] * m[3] * m[6];
43  inv[7] = m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] +
44  m[4] * m[3] * m[10] + m[8] * m[2] * m[7] - m[8] * m[3] * m[6];
45  inv[11] = -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] -
46  m[4] * m[3] * m[9] - m[8] * m[1] * m[7] + m[8] * m[3] * m[5];
47  inv[15] = m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] +
48  m[4] * m[2] * m[9] + m[8] * m[1] * m[6] - m[8] * m[2] * m[5];
49 
50  *det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
51 
52  if(*det == 0) return false;
53 
54  double invDet = 1.0 / *det;
55 
56  for(i = 0; i < 16; i++) inv[i] *= invDet;
57 
58  return true;
59 }
60 
61 static void solveEig(double A, double B, double C, double D, double *lambda1,
62  double *v1x, double *v1y, double *lambda2, double *v2x,
63  double *v2y)
64 {
65  if(B * C <= tolerance) {
66  *lambda1 = A;
67  *v1x = 1;
68  *v1y = 0;
69  *lambda2 = D;
70  *v2x = 0;
71  *v2y = 1;
72  return;
73  }
74 
75  double tr = A + D;
76  double det = A * D - B * C;
77  double S = sqrt((tr * tr / 4) - det);
78  *lambda1 = tr / 2 + S;
79  *lambda2 = tr / 2 - S;
80 
81  double temp = ((A - D) * (A - D) / 4) + B * C;
82 
83  double SS = sqrt(temp > 0 ? temp : 0.0);
84  if(A - D < 0) {
85  *v1x = C;
86  *v1y = -(A - D) / 2 + SS;
87  *v2x = +(A - D) / 2 - SS;
88  *v2y = B;
89  }
90  else {
91  *v2x = C;
92  *v2y = -(A - D) / 2 - SS;
93  *v1x = +(A - D) / 2 + SS;
94  *v1y = B;
95  }
96 
97  double n1 = sqrt((*v1x) * (*v1x) + (*v1y) * (*v1y));
98  *v1x /= n1;
99  *v1y /= n1;
100  double n2 = sqrt((*v2x) * (*v2x) + (*v2y) * (*v2y));
101  *v2x /= n2;
102  *v2y /= n2;
103 }
104 
105 static inline int node2trianglescmp(const void *p0, const void *p1)
106 {
107  const std::size_t *e0 = (const std::size_t *)p0;
108  const std::size_t *e1 = (const std::size_t *)p1;
109 
110  if(e0[0] < e1[0]) return -1;
111  if(e0[0] > e1[0]) return 1;
112  return 0;
113 }
114 
115 static inline double normalize(double *n)
116 {
117  double d = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
118  if(d != 0.0) {
119  n[0] /= d;
120  n[1] /= d;
121  n[2] /= d;
122  }
123  return fabs(d);
124 }
125 
126 static inline void makevector(const double *a, const double *b, double *ba)
127 {
128  ba[0] = b[0] - a[0];
129  ba[1] = b[1] - a[1];
130  ba[2] = b[2] - a[2];
131 }
132 
133 static inline double dotprod(double *a, double *b)
134 {
135  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
136 }
137 
138 static inline void crossprod(double *a, double *b, double *n)
139 {
140  n[0] = a[1] * b[2] - a[2] * b[1];
141  n[1] = -(a[0] * b[2] - a[2] * b[0]);
142  n[2] = a[0] * b[1] - a[1] * b[0];
143 }
144 
145 static inline void unitNormal2Triangle(const double *x1, const double *x2,
146  const double *x3, double *n, double *s)
147 {
148  double a[3] = {x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2]};
149  double b[3] = {x3[0] - x1[0], x3[1] - x1[1], x3[2] - x1[2]};
150 
151  crossprod(a, b, n);
152  *s = 2 * normalize(n);
153 }
154 
155 static void computeLocalFrame(double *n, double *u, double *v)
156 {
157  if(fabs(n[0]) > fabs(n[1]) && fabs(n[0]) > fabs(n[2])) {
158  u[0] = 0;
159  u[1] = 0;
160  u[2] = 1;
161  }
162  else {
163  u[0] = 1;
164  u[1] = 0;
165  u[2] = 0;
166  }
167  crossprod(n, u, v);
168  normalize(v);
169  crossprod(v, n, u);
170 }
171 
172 bool CurvatureRusinkiewicz(
173  const std::vector<int> &triangles, const std::vector<SPoint3> &nodes,
174  std::vector<std::pair<SVector3, SVector3> > &nodalCurvatures)
175 {
176  std::size_t nTriangles = triangles.size() / 3;
177  std::size_t nVertices = nodes.size();
178 
179  nodalCurvatures.resize(nVertices);
180 
181  // for(int i = 0; i < nodes.size(); ++i){
182  // nodalCurvatures[i] = std::make_pair(0.0,0.0);
183  // }
184 
185  std::vector<std::size_t> node2tri(3 * 2 * nTriangles);
186 
187  // first compute node normals and node-to-triangle connectivity
188 
189  std::vector<double> nodeNormals(3 * nVertices, 0.);
190  std::size_t counter = 0;
191  double n[3], surf;
192 
193  for(std::size_t i = 0; i < nTriangles; i++) {
194  node2tri[counter++] = triangles[3 * i + 0];
195  node2tri[counter++] = i;
196  node2tri[counter++] = triangles[3 * i + 1];
197  node2tri[counter++] = i;
198  node2tri[counter++] = triangles[3 * i + 2];
199  node2tri[counter++] = i;
200  unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
201  nodes[triangles[3 * i + 1]].data(),
202  nodes[triangles[3 * i + 2]].data(), n, &surf);
203  double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
204  double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
205  double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
206  for(std::size_t i1 = 0; i1 < 3; i1++) {
207  n0[i1] += n[i1];
208  n1[i1] += n[i1];
209  n2[i1] += n[i1];
210  }
211  }
212  for(std::size_t i = 0; i < nVertices; i++) normalize(&nodeNormals[3 * i]);
213 
214  qsort(&node2tri[0], 3 * nTriangles, 2 * sizeof(std::size_t),
215  node2trianglescmp);
216 
217  // compute the second fundamental tensor on each triangle using least squares
218 
219  std::vector<double> CURV(4 * nTriangles);
220 
221  for(std::size_t i = 0; i < nTriangles; i++) {
222  double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
223  double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
224  double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
225 
226  double e0[3], e1[3], e2[3];
227  makevector(nodes[triangles[3 * i + 2]].data(),
228  nodes[triangles[3 * i + 1]].data(), e0);
229  makevector(nodes[triangles[3 * i + 0]].data(),
230  nodes[triangles[3 * i + 2]].data(), e1);
231  makevector(nodes[triangles[3 * i + 1]].data(),
232  nodes[triangles[3 * i + 0]].data(), e2);
233  unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
234  nodes[triangles[3 * i + 1]].data(),
235  nodes[triangles[3 * i + 2]].data(), n, &surf);
236  double u[3], v[3];
237  makevector(nodes[triangles[3 * i + 0]].data(),
238  nodes[triangles[3 * i + 1]].data(), u);
239  normalize(u);
240  crossprod(n, u, v);
241 
242  double sys[6][4], rhs[6], temp[3], invA[16], A[16], B[4];
243  sys[0][0] = sys[1][2] = dotprod(e0, u);
244  sys[0][1] = sys[1][3] = dotprod(e0, v);
245  sys[0][2] = sys[0][3] = sys[1][0] = sys[1][1] = 0;
246  makevector(n2, n1, temp);
247  rhs[0] = dotprod(temp, u);
248  rhs[1] = dotprod(temp, v);
249 
250  sys[2][0] = sys[3][2] = dotprod(e1, u);
251  sys[2][1] = sys[3][3] = dotprod(e1, v);
252  sys[2][2] = sys[2][3] = sys[3][0] = sys[3][1] = 0;
253  makevector(n0, n2, temp);
254  rhs[2] = dotprod(temp, u);
255  rhs[3] = dotprod(temp, v);
256 
257  sys[4][0] = sys[5][2] = dotprod(e2, u);
258  sys[4][1] = sys[5][3] = dotprod(e2, v);
259  sys[4][2] = sys[4][3] = sys[5][0] = sys[5][1] = 0;
260  makevector(n1, n0, temp);
261  rhs[4] = dotprod(temp, u);
262  rhs[5] = dotprod(temp, v);
263 
264  for(std::size_t i1 = 0; i1 < 4; i1++) {
265  B[i1] = 0.0;
266  for(std::size_t i3 = 0; i3 < 6; i3++) { B[i1] += sys[i3][i1] * rhs[i3]; }
267  for(std::size_t i2 = 0; i2 < 4; i2++) {
268  A[i1 + 4 * i2] = 0.0;
269  for(std::size_t i3 = 0; i3 < 6; i3++) {
270  A[i1 + 4 * i2] += sys[i3][i2] * sys[i3][i1];
271  }
272  }
273  }
274  double det;
275  Inv4x4ColumnMajor(A, invA, &det);
276  for(std::size_t i1 = 0; i1 < 4; i1++) {
277  CURV[4 * i + i1] = 0.0;
278  for(std::size_t i2 = 0; i2 < 4; i2++) {
279  CURV[4 * i + i1] += invA[i1 + 4 * i2] * B[i2];
280  }
281  }
282  CURV[4 * i + 1] = CURV[4 * i + 2] =
283  0.5 * (CURV[4 * i + 1] + CURV[4 * i + 2]);
284  }
285 
286  uint64_t count = 0;
287  double uP[3] = {0., 0., 0.}, vP[3] = {0., 0., 0.};
288  double A = 0., B = 0., D = 0.;
289  uint64_t iTriangle;
290  uint64_t ind = 0;
291  for(uint64_t iVert = 0; iVert < nVertices; ++iVert) {
292  count = 0;
293  A = 0.0;
294  B = 0.0;
295  D = 0.0;
296  while(node2tri[2 * ind + 0] == iVert) {
297  iTriangle = node2tri[2 * ind + 1];
298  computeLocalFrame(&nodeNormals[3 * iVert], uP, vP);
299  // computing each curvature around a vertex
300  unitNormal2Triangle(nodes[triangles[3 * iTriangle + 0]].data(),
301  nodes[triangles[3 * iTriangle + 1]].data(),
302  nodes[triangles[3 * iTriangle + 2]].data(), n, &surf);
303  double uF[3], vF[3];
304  makevector(nodes[triangles[3 * iTriangle + 0]].data(),
305  nodes[triangles[3 * iTriangle + 1]].data(), uF);
306  normalize(uF);
307  crossprod(n, uF, vF);
308  double *c = &CURV[4 * iTriangle];
309  double UP[3] = {dotprod(uP, uF), dotprod(uP, vF), 0};
310  normalize(UP);
311  double VP[3] = {dotprod(vP, uF), dotprod(vP, vF), 0};
312  normalize(VP);
313  A += (UP[0] * UP[0] * c[0] + 2 * UP[0] * UP[1] * c[1] +
314  UP[1] * UP[1] * c[3]);
315  D += (VP[0] * VP[0] * c[0] + 2 * VP[0] * VP[1] * c[1] +
316  VP[1] * VP[1] * c[3]);
317  B += (VP[0] * UP[0] * c[0] + (VP[1] * UP[0] + VP[0] * UP[1]) * c[1] +
318  VP[1] * UP[1] * c[3]);
319  ++count;
320  ++ind;
321  if(ind >= 3 * nTriangles) break;
322  }
323 
324  A /= (double)count;
325  B /= (double)count;
326  D /= (double)count;
327  double lambda1, lambda2, v1x, v1y, v2x, v2y;
328  solveEig(A, B, B, D, &lambda1, &v1x, &v1y, &lambda2, &v2x, &v2y);
329  SVector3 cMax(fabs(lambda1) * (v1x * uP[0] + v1y * vP[0]),
330  fabs(lambda1) * (v1x * uP[1] + v1y * vP[1]),
331  fabs(lambda1) * (v1x * uP[2] + v1y * vP[2]));
332  SVector3 cMin(fabs(lambda2) * (v2x * uP[0] + v2y * vP[0]),
333  fabs(lambda2) * (v2x * uP[1] + v2y * vP[1]),
334  fabs(lambda2) * (v2x * uP[2] + v2y * vP[2]));
335  nodalCurvatures[iVert] = std::make_pair(cMax, cMin);
336  } // for iVert
337 
338  return true;
339 }
340 
341 bool CurvatureRusinkiewicz(
342  const std::vector<int> &triangles, const std::vector<SPoint3> &nodes,
343  std::vector<std::pair<SVector3, SVector3> > &nodalCurvatures,
344  std::vector<double> &nodeNormals)
345 {
346  std::size_t nTriangles = triangles.size() / 3;
347  std::size_t nVertices = nodes.size();
348 
349  nodalCurvatures.resize(nVertices);
350 
351  // for(int i = 0; i < nodes.size(); ++i){
352  // nodalCurvatures[i] = std::make_pair(0.0,0.0);
353  // }
354 
355  std::vector<std::size_t> node2tri(3 * 2 * nTriangles);
356 
357  // first compute node normals and node-to-triangle connectivity
358 
359  // std::vector<double> nodeNormals(3 * nVertices, 0.);
360  std::size_t counter = 0;
361  double n[3], surf;
362 
363  for(std::size_t i = 0; i < nTriangles; i++) {
364  node2tri[counter++] = triangles[3 * i + 0];
365  node2tri[counter++] = i;
366  node2tri[counter++] = triangles[3 * i + 1];
367  node2tri[counter++] = i;
368  node2tri[counter++] = triangles[3 * i + 2];
369  node2tri[counter++] = i;
370  unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
371  nodes[triangles[3 * i + 1]].data(),
372  nodes[triangles[3 * i + 2]].data(), n, &surf);
373  double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
374  double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
375  double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
376  for(std::size_t i1 = 0; i1 < 3; i1++) {
377  n0[i1] += n[i1];
378  n1[i1] += n[i1];
379  n2[i1] += n[i1];
380  }
381  }
382  for(std::size_t i = 0; i < nVertices; i++) normalize(&nodeNormals[3 * i]);
383 
384  qsort(&node2tri[0], 3 * nTriangles, 2 * sizeof(std::size_t),
385  node2trianglescmp);
386 
387  // compute the second fundamental tensor on each triangle using least squares
388 
389  std::vector<double> CURV(4 * nTriangles);
390 
391  for(std::size_t i = 0; i < nTriangles; i++) {
392  double *n0 = &nodeNormals[3 * triangles[3 * i + 0]];
393  double *n1 = &nodeNormals[3 * triangles[3 * i + 1]];
394  double *n2 = &nodeNormals[3 * triangles[3 * i + 2]];
395 
396  double e0[3], e1[3], e2[3];
397  makevector(nodes[triangles[3 * i + 2]].data(),
398  nodes[triangles[3 * i + 1]].data(), e0);
399  makevector(nodes[triangles[3 * i + 0]].data(),
400  nodes[triangles[3 * i + 2]].data(), e1);
401  makevector(nodes[triangles[3 * i + 1]].data(),
402  nodes[triangles[3 * i + 0]].data(), e2);
403  unitNormal2Triangle(nodes[triangles[3 * i + 0]].data(),
404  nodes[triangles[3 * i + 1]].data(),
405  nodes[triangles[3 * i + 2]].data(), n, &surf);
406  double u[3], v[3];
407  makevector(nodes[triangles[3 * i + 0]].data(),
408  nodes[triangles[3 * i + 1]].data(), u);
409  normalize(u);
410  crossprod(n, u, v);
411 
412  double sys[6][4], rhs[6], temp[3], invA[16], A[16], B[4];
413  sys[0][0] = sys[1][2] = dotprod(e0, u);
414  sys[0][1] = sys[1][3] = dotprod(e0, v);
415  sys[0][2] = sys[0][3] = sys[1][0] = sys[1][1] = 0;
416  makevector(n2, n1, temp);
417  rhs[0] = dotprod(temp, u);
418  rhs[1] = dotprod(temp, v);
419 
420  sys[2][0] = sys[3][2] = dotprod(e1, u);
421  sys[2][1] = sys[3][3] = dotprod(e1, v);
422  sys[2][2] = sys[2][3] = sys[3][0] = sys[3][1] = 0;
423  makevector(n0, n2, temp);
424  rhs[2] = dotprod(temp, u);
425  rhs[3] = dotprod(temp, v);
426 
427  sys[4][0] = sys[5][2] = dotprod(e2, u);
428  sys[4][1] = sys[5][3] = dotprod(e2, v);
429  sys[4][2] = sys[4][3] = sys[5][0] = sys[5][1] = 0;
430  makevector(n1, n0, temp);
431  rhs[4] = dotprod(temp, u);
432  rhs[5] = dotprod(temp, v);
433 
434  for(std::size_t i1 = 0; i1 < 4; i1++) {
435  B[i1] = 0.0;
436  for(std::size_t i3 = 0; i3 < 6; i3++) { B[i1] += sys[i3][i1] * rhs[i3]; }
437  for(std::size_t i2 = 0; i2 < 4; i2++) {
438  A[i1 + 4 * i2] = 0.0;
439  for(std::size_t i3 = 0; i3 < 6; i3++) {
440  A[i1 + 4 * i2] += sys[i3][i2] * sys[i3][i1];
441  }
442  }
443  }
444  double det;
445  Inv4x4ColumnMajor(A, invA, &det);
446  for(std::size_t i1 = 0; i1 < 4; i1++) {
447  CURV[4 * i + i1] = 0.0;
448  for(std::size_t i2 = 0; i2 < 4; i2++) {
449  CURV[4 * i + i1] += invA[i1 + 4 * i2] * B[i2];
450  }
451  }
452  CURV[4 * i + 1] = CURV[4 * i + 2] =
453  0.5 * (CURV[4 * i + 1] + CURV[4 * i + 2]);
454  }
455 
456  uint64_t count = 0;
457  double uP[3] = {0., 0., 0.}, vP[3] = {0., 0., 0.};
458  double A = 0., B = 0., D = 0.;
459  uint64_t iTriangle;
460  uint64_t ind = 0;
461  for(uint64_t iVert = 0; iVert < nVertices; ++iVert) {
462  count = 0;
463  A = 0.0;
464  B = 0.0;
465  D = 0.0;
466  while(node2tri[2 * ind + 0] == iVert) {
467  iTriangle = node2tri[2 * ind + 1];
468  computeLocalFrame(&nodeNormals[3 * iVert], uP, vP);
469  // computing each curvature around a vertex
470  unitNormal2Triangle(nodes[triangles[3 * iTriangle + 0]].data(),
471  nodes[triangles[3 * iTriangle + 1]].data(),
472  nodes[triangles[3 * iTriangle + 2]].data(), n, &surf);
473  double uF[3], vF[3];
474  makevector(nodes[triangles[3 * iTriangle + 0]].data(),
475  nodes[triangles[3 * iTriangle + 1]].data(), uF);
476  normalize(uF);
477  crossprod(n, uF, vF);
478  double *c = &CURV[4 * iTriangle];
479  double UP[3] = {dotprod(uP, uF), dotprod(uP, vF), 0};
480  normalize(UP);
481  double VP[3] = {dotprod(vP, uF), dotprod(vP, vF), 0};
482  normalize(VP);
483  A += (UP[0] * UP[0] * c[0] + 2 * UP[0] * UP[1] * c[1] +
484  UP[1] * UP[1] * c[3]);
485  D += (VP[0] * VP[0] * c[0] + 2 * VP[0] * VP[1] * c[1] +
486  VP[1] * VP[1] * c[3]);
487  B += (VP[0] * UP[0] * c[0] + (VP[1] * UP[0] + VP[0] * UP[1]) * c[1] +
488  VP[1] * UP[1] * c[3]);
489  ++count;
490  ++ind;
491  if(ind >= 3 * nTriangles) break;
492  }
493 
494  A /= (double)count;
495  B /= (double)count;
496  D /= (double)count;
497  double lambda1, lambda2, v1x, v1y, v2x, v2y;
498  solveEig(A, B, B, D, &lambda1, &v1x, &v1y, &lambda2, &v2x, &v2y);
499  SVector3 cMax(fabs(lambda1) * (v1x * uP[0] + v1y * vP[0]),
500  fabs(lambda1) * (v1x * uP[1] + v1y * vP[1]),
501  fabs(lambda1) * (v1x * uP[2] + v1y * vP[2]));
502  SVector3 cMin(fabs(lambda2) * (v2x * uP[0] + v2y * vP[0]),
503  fabs(lambda2) * (v2x * uP[1] + v2y * vP[1]),
504  fabs(lambda2) * (v2x * uP[2] + v2y * vP[2]));
505  nodalCurvatures[iVert] = std::make_pair(cMax, cMin);
506  } // for iVert
507 
508  return true;
509 }
makevector
static void makevector(const double *a, const double *b, double *ba)
Definition: curvature.cpp:126
normalize
static double normalize(double *n)
Definition: curvature.cpp:115
D
#define D
Definition: DefaultOptions.h:24
unitNormal2Triangle
static void unitNormal2Triangle(const double *x1, const double *x2, const double *x3, double *n, double *s)
Definition: curvature.cpp:145
computeLocalFrame
static void computeLocalFrame(double *n, double *u, double *v)
Definition: curvature.cpp:155
curvature.h
crossprod
static void crossprod(double *a, double *b, double *n)
Definition: curvature.cpp:138
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
SVector3
Definition: SVector3.h:16
node2trianglescmp
static int node2trianglescmp(const void *p0, const void *p1)
Definition: curvature.cpp:105
dotprod
static double dotprod(double *a, double *b)
Definition: curvature.cpp:133
tolerance
#define tolerance
Definition: curvature.cpp:11
S
#define S
Definition: DefaultOptions.h:21
Inv4x4ColumnMajor
static bool Inv4x4ColumnMajor(double m[16], double inv[16], double *det)
Definition: curvature.cpp:13
solveEig
static void solveEig(double A, double B, double C, double D, double *lambda1, double *v1x, double *v1y, double *lambda2, double *v2x, double *v2y)
Definition: curvature.cpp:61
CurvatureRusinkiewicz
bool CurvatureRusinkiewicz(const std::vector< int > &triangles, const std::vector< SPoint3 > &nodes, std::vector< std::pair< SVector3, SVector3 > > &nodalCurvatures)
Definition: curvature.cpp:172