1 #include "hxt_curvature.h"
2 #include "hxt_tools.h"
3 #include <stdio.h>
4
5 // #include <time.h> // for commented timings
6
7
8 #define tolerance 0.1e-20
9
hxtInv4x4ColumnMajor(double m[16],double invOut[16],double * det)10 static HXTStatus hxtInv4x4ColumnMajor(double m[16], double invOut[16], double *det)
11 {
12 double inv[16];
13 int i;
14
15 inv[ 0] = m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] + m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10];
16 inv[ 4] = -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] - m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10];
17 inv[ 8] = m[4] * m[ 9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] + m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[ 9];
18 inv[12] = -m[4] * m[ 9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] - m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[ 9];
19 inv[ 1] = -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] - m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10];
20 inv[ 5] = m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] + m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10];
21 inv[ 9] = -m[0] * m[ 9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] - m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[ 9];
22 inv[13] = m[0] * m[ 9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] + m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[ 9];
23 inv[ 2] = m[1] * m[ 6] * m[15] - m[1] * m[ 7] * m[14] - m[5] * m[2] * m[15] + m[5] * m[3] * m[14] + m[13] * m[2] * m[ 7] - m[13] * m[3] * m[ 6];
24 inv[ 6] = -m[0] * m[ 6] * m[15] + m[0] * m[ 7] * m[14] + m[4] * m[2] * m[15] - m[4] * m[3] * m[14] - m[12] * m[2] * m[ 7] + m[12] * m[3] * m[ 6];
25 inv[10] = m[0] * m[ 5] * m[15] - m[0] * m[ 7] * m[13] - m[4] * m[1] * m[15] + m[4] * m[3] * m[13] + m[12] * m[1] * m[ 7] - m[12] * m[3] * m[ 5];
26 inv[14] = -m[0] * m[ 5] * m[14] + m[0] * m[ 6] * m[13] + m[4] * m[1] * m[14] - m[4] * m[2] * m[13] - m[12] * m[1] * m[ 6] + m[12] * m[2] * m[ 5];
27 inv[ 3] = -m[1] * m[ 6] * m[11] + m[1] * m[ 7] * m[10] + m[5] * m[2] * m[11] - m[5] * m[3] * m[10] - m[ 9] * m[2] * m[ 7] + m[ 9] * m[3] * m[ 6];
28 inv[ 7] = m[0] * m[ 6] * m[11] - m[0] * m[ 7] * m[10] - m[4] * m[2] * m[11] + m[4] * m[3] * m[10] + m[ 8] * m[2] * m[ 7] - m[ 8] * m[3] * m[ 6];
29 inv[11] = -m[0] * m[ 5] * m[11] + m[0] * m[ 7] * m[ 9] + m[4] * m[1] * m[11] - m[4] * m[3] * m[ 9] - m[ 8] * m[1] * m[ 7] + m[ 8] * m[3] * m[ 5];
30 inv[15] = m[0] * m[ 5] * m[10] - m[0] * m[ 6] * m[ 9] - m[4] * m[1] * m[10] + m[4] * m[2] * m[ 9] + m[ 8] * m[1] * m[ 6] - m[ 8] * m[2] * m[ 5];
31
32 *det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
33
34 if(*det == 0)
35 return HXT_STATUS_ERROR;
36
37 double invDet = 1.0 / *det;
38
39 for(i = 0; i < 16; i++)
40 invOut[i] = inv[i] * invDet;
41
42 return HXT_STATUS_OK;
43 }
44
solveEig(double A,double B,double C,double D,double * lambda1,double * v1x,double * v1y,double * lambda2,double * v2x,double * v2y)45 static void solveEig(double A, double B, double C, double D,
46 double* lambda1, double *v1x, double*v1y,
47 double* lambda2, double *v2x, double*v2y )
48 {
49 if(B*C <= tolerance ) {
50 *lambda1 = A; *v1x = 1; *v1y = 0;
51 *lambda2 = D; *v2x = 0; *v2y = 1;
52 return;
53 }
54
55 double tr = A + D;
56 double det = A * D - B * C;
57 double S = sqrt( (tr*tr/4) - det );
58 *lambda1 = tr/2 + S;
59 *lambda2 = tr/2 - S;
60
61 double temp = ((A-D)*(A-D)/4) + B * C;
62
63 double SS = sqrt( temp > 0 ? temp: 0.0);
64 if( A - D < 0 ) {
65 *v1x = C;
66 *v1y = - (A-D)/2 + SS;
67 *v2x = + (A-D)/2 - SS;
68 *v2y = B;
69 } else {
70 *v2x = C;
71 *v2y = - (A-D)/2 - SS;
72 *v1x = + (A-D)/2 + SS;
73 *v1y = B;
74 }
75
76 double n1 = sqrt((*v1x)*(*v1x)+(*v1y)*(*v1y));
77 *v1x /= n1; *v1y /= n1;
78 double n2 = sqrt((*v2x)*(*v2x)+(*v2y)*(*v2y));
79 *v2x /= n2; *v2y /= n2;
80 }
81
node2trianglescmp(const void * p0,const void * p1)82 static inline int node2trianglescmp(const void *p0, const void *p1)
83 {
84 const uint64_t *e0 = (const uint64_t*)p0;
85 const uint64_t *e1 = (const uint64_t*)p1;
86
87 if (e0[0] < e1[0]) return -1;
88 if (e0[0] > e1[0]) return 1;
89 return 0;
90 }
91
92
normalize(double * n)93 static inline double normalize (double *n){
94 double d = sqrt (n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
95 if (d != 0.0){
96 n[0] /= d;
97 n[1] /= d;
98 n[2] /= d;
99 }
100 return fabs(d);
101 }
102
makevector(double * a,double * b,double * ba)103 static inline void makevector (double *a, double *b, double *ba){
104 ba[0] = b[0]-a[0];
105 ba[1] = b[1]-a[1];
106 ba[2] = b[2]-a[2];
107 }
108
dotprod(double * a,double * b)109 static inline double dotprod (double *a, double *b){
110 return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
111 }
112
crossprod(double * a,double * b,double * n)113 static inline void crossprod (double *a, double *b, double *n){
114 n[0] = a[1]*b[2] - a[2]*b[1];
115 n[1] = - (a[0]*b[2] - a[2]*b[0]);
116 n[2] = a[0]*b[1] - a[1]*b[0];
117 }
118
unitNormal2Triangle(double * x1,double * x2,double * x3,double * n,double * s)119 static inline void unitNormal2Triangle (double *x1, double *x2, double *x3, double *n, double *s) {
120
121 double a[3] = {x2[0]-x1[0],x2[1]-x1[1],x2[2]-x1[2]};
122 double b[3] = {x3[0]-x1[0],x3[1]-x1[1],x3[2]-x1[2]};
123
124
125 crossprod (a,b,n);
126 *s = 2*normalize(n);
127 // printf("%g %g %g = %g %g %g\n",(x1[0]+x2[0]+x3[0])/3.,
128 // (x1[1]+x2[1]+x3[1])/3.,(x1[2]+x2[2]+x3[2])/3.,
129 // n[0],n[1],n[2]);
130 }
131
computeLocalFrame(double * n,double * u,double * v)132 void computeLocalFrame (double *n, double *u, double *v) {
133 if (fabs(n[0]) > fabs(n[1]) && fabs(n[0]) > fabs(n[2])){
134 u[0] = 0;u[1] = 0; u[2] = 1;
135 }
136 else {
137 u[0] = 1;u[1] = 0; u[2] = 0;
138 }
139 crossprod (n,u,v);
140 normalize (v);
141 crossprod (v,n,u);
142 // printf("%g %g %g\n",v[0],v[1],v[2]);
143 // printf("%g %g %g\n",u[0],u[1],u[2]);
144 // printf("%g %g %g\n",v[0],v[1],v[2]);
145 }
146
147
saveNodalField(HXTMesh * mesh,double * v,int ncomp,const char * fn)148 void saveNodalField (HXTMesh *mesh, double *v, int ncomp, const char *fn){
149 FILE *f = fopen (fn , "w");
150 uint64_t nTriangles = mesh->triangles.num;
151 fprintf(f,"View \"Vector\" {\n");
152 for (uint64_t i = 0; i<nTriangles; i++){
153 double *p0 = mesh->vertices.coord + 4*mesh->triangles.node[3*i+0];
154 double *p1 = mesh->vertices.coord + 4*mesh->triangles.node[3*i+1];
155 double *p2 = mesh->vertices.coord + 4*mesh->triangles.node[3*i+2];
156 double *v0 = &v[abs(ncomp)*mesh->triangles.node[3*i+0]];
157 double *v1 = &v[abs(ncomp)*mesh->triangles.node[3*i+1]];
158 double *v2 = &v[abs(ncomp)*mesh->triangles.node[3*i+2]];
159 if (ncomp == 3)
160 fprintf(f,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
161 p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],
162 v0[0],v0[1],v0[2],v1[0],v1[1],v1[2],v2[0],v2[1],v2[2]);
163 else if (ncomp == 1)
164 fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
165 p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],
166 v0[0],v1[0],v2[0]);
167 else if (ncomp == 6)
168 fprintf(f,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
169 p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],
170 v0[0],v0[1],v0[2],v1[0],v1[1],v1[2],v2[0],v2[1],v2[2]);
171 else if (ncomp == -6)
172 fprintf(f,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
173 p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],
174 v0[3],v0[4],v0[5],v1[3],v1[4],v1[5],v2[3],v2[4],v2[5]);
175 }
176 fprintf(f,"};\n");
177
178 fclose (f);
179 }
180
hxtCurvatureRusinkiewicz(HXTMesh * mesh,double ** nodalCurvatures,double ** crossField,HXTEdges * edges,int debug)181 HXTStatus hxtCurvatureRusinkiewicz (HXTMesh *mesh, double **nodalCurvatures, double **crossField, HXTEdges* edges, int debug)
182 {
183 // clock_t T1 = clock();
184
185 uint64_t nTriangles = mesh->triangles.num;
186 // uint64_t nEdgesBdry = mesh->lines.num;
187 uint64_t nVertices = mesh->vertices.num;
188
189 HXT_CHECK(hxtMalloc(nodalCurvatures,6*nVertices*sizeof(double)));
190
191 // FIXME: not all nodalCurvatures are initialized by the following algorithm,
192 // => it result in an undefined behavior when computing cross-field (line 374+)
193 // To at least have a deterministic result, I decided to initialize nodalCurvature to 0
194 memset(*nodalCurvatures, 0, 6*nVertices*sizeof(double));
195
196 uint64_t *node2tri;
197 HXT_CHECK(hxtMalloc(&node2tri,3*2*nTriangles*sizeof(uint64_t)));
198
199 // --> FIRST COMPUTE NODE NORMALS and node2triangle connectivity
200 double *nodeNormals;
201 HXT_CHECK(hxtMalloc(&nodeNormals,3*nVertices*sizeof(double)));
202 for (uint64_t i = 0; i<3*nVertices; i++)nodeNormals[i] = 0.0;
203 uint64_t counter = 0;
204
205 double n[3],surf;
206
207 for (uint64_t i = 0; i<nTriangles; i++){
208 node2tri[counter++] = mesh->triangles.node[3*i+0];
209 node2tri[counter++] = i;
210 node2tri[counter++] = mesh->triangles.node[3*i+1];
211 node2tri[counter++] = i;
212 node2tri[counter++] = mesh->triangles.node[3*i+2];
213 node2tri[counter++] = i;
214 unitNormal2Triangle ( mesh->vertices.coord + 4*mesh->triangles.node[3*i+0],
215 mesh->vertices.coord + 4*mesh->triangles.node[3*i+1],
216 mesh->vertices.coord + 4*mesh->triangles.node[3*i+2], n, &surf);
217 double *n0 = &nodeNormals[3*mesh->triangles.node[3*i+0]];
218 double *n1 = &nodeNormals[3*mesh->triangles.node[3*i+1]];
219 double *n2 = &nodeNormals[3*mesh->triangles.node[3*i+2]];
220 for (uint64_t i1 = 0; i1 < 3; i1++) {
221 n0[i1]+=n[i1];
222 n1[i1]+=n[i1];
223 n2[i1]+=n[i1];
224 }
225 }
226 for (uint64_t i = 0; i<nVertices; i++)normalize(&nodeNormals[3*i]);
227 if (debug) saveNodalField (mesh,nodeNormals, 3, "normals.pos");
228
229 qsort(node2tri,3*nTriangles,2*sizeof(uint64_t),node2trianglescmp);
230
231
232 // --> COMPUTE THE SECOND FUNDAMENTAL TENSOR ON EACH TRIANGLE USING LEAST SQUARES
233
234 double *CURV;
235 HXT_CHECK(hxtMalloc(&CURV,4*nTriangles*sizeof(double)));
236
237 for (uint64_t i = 0; i<nTriangles; i++){
238
239 double *n0 = &nodeNormals[3*mesh->triangles.node[3*i+0]];
240 double *n1 = &nodeNormals[3*mesh->triangles.node[3*i+1]];
241 double *n2 = &nodeNormals[3*mesh->triangles.node[3*i+2]];
242
243 double e0[3],e1[3],e2[3];
244 makevector (mesh->vertices.coord + 4*mesh->triangles.node[3*i+2],
245 mesh->vertices.coord + 4*mesh->triangles.node[3*i+1],e0);
246 makevector (mesh->vertices.coord + 4*mesh->triangles.node[3*i+0],
247 mesh->vertices.coord + 4*mesh->triangles.node[3*i+2],e1);
248 makevector (mesh->vertices.coord + 4*mesh->triangles.node[3*i+1],
249 mesh->vertices.coord + 4*mesh->triangles.node[3*i+0],e2);
250
251 unitNormal2Triangle ( mesh->vertices.coord + 4*mesh->triangles.node[3*i+0],
252 mesh->vertices.coord + 4*mesh->triangles.node[3*i+1],
253 mesh->vertices.coord + 4*mesh->triangles.node[3*i+2], n, &surf);
254 double u[3],v[3];
255 makevector (mesh->vertices.coord + 4*mesh->triangles.node[3*i+0],
256 mesh->vertices.coord + 4*mesh->triangles.node[3*i+1],u);
257 normalize(u);
258 crossprod(n,u,v);
259
260 double sys[6][4], rhs[6], temp[3], invA[16], A[16], B[4];
261 sys[0][0] = sys[1][2] = dotprod(e0,u);
262 sys[0][1] = sys[1][3] = dotprod(e0,v);
263 sys[0][2] = sys[0][3] = sys[1][0] = sys[1][1] = 0;
264 makevector (n2,n1,temp);
265 rhs[0] = dotprod(temp,u);
266 rhs[1] = dotprod(temp,v);
267
268 sys[2][0] = sys[3][2] = dotprod(e1,u);
269 sys[2][1] = sys[3][3] = dotprod(e1,v);
270 sys[2][2] = sys[2][3] = sys[3][0] = sys[3][1] = 0;
271 makevector (n0,n2,temp);
272 rhs[2] = dotprod(temp,u);
273 rhs[3] = dotprod(temp,v);
274
275 sys[4][0] = sys[5][2] = dotprod(e2,u);
276 sys[4][1] = sys[5][3] = dotprod(e2,v);
277 sys[4][2] = sys[4][3] = sys[5][0] = sys[5][1] = 0;
278 makevector (n1,n0,temp);
279 rhs[4] = dotprod(temp,u);
280 rhs[5] = dotprod(temp,v);
281
282 for (uint64_t i1=0;i1<4;i1++){
283 B[i1] = 0.0;
284 for (uint64_t i3=0;i3<6;i3++){
285 B[i1] += sys[i3][i1] * rhs[i3];
286 }
287 for (uint64_t i2=0;i2<4;i2++){
288 A[i1+4*i2] = 0.0;
289 for (uint64_t i3=0;i3<6;i3++){
290 A[i1+4*i2] += sys[i3][i2] * sys[i3][i1];
291 }
292 }
293 }
294 double det;
295 hxtInv4x4ColumnMajor(A, invA,&det);
296 for (uint64_t i1=0;i1<4;i1++){
297 CURV[4*i+i1] = 0.0;
298 for (uint64_t i2=0;i2<4;i2++){
299 CURV[4*i+i1] += invA[i1+4*i2] * B[i2];
300 }
301 }
302 CURV[4*i+1] = CURV[4*i+2] = 0.5* (CURV[4*i+1] + CURV[4*i+2]);
303 }
304
305 // Get vertex curvatures by averaging triangle curvatures
306 uint64_t currentVertex = nVertices + 1;
307 uint64_t count = 0;
308 double uP[3],vP[3], A, B, D;
309 for (uint64_t i = 0; i<6*nTriangles; i+=2){
310 uint64_t iVertex = node2tri[i];
311 uint64_t iTriangle = node2tri[i+1];
312 // printf("%d %d %d %d\n",iVertex,iTriangle,nVertices,currentVertex);
313 if (currentVertex != iVertex){
314 // compute the real stuff
315 if (currentVertex != nVertices + 1){
316 // printf("%g %g %g %d \n",A,B,D,count);
317 A /= (double) count;
318 B /= (double) count;
319 D /= (double) count;
320 double lambda1, lambda2, v1x, v1y, v2x, v2y;
321 solveEig(A, B, B, D,
322 & lambda1, & v1x, &v1y,
323 & lambda2, & v2x, & v2y );
324 (*nodalCurvatures) [6 * currentVertex + 0] = fabs(lambda1) * (v1x * uP[0] + v1y * vP[0]);
325 (*nodalCurvatures) [6 * currentVertex + 1] = fabs(lambda1) * (v1x * uP[1] + v1y * vP[1]);
326 (*nodalCurvatures) [6 * currentVertex + 2] = fabs(lambda1) * (v1x * uP[2] + v1y * vP[2]);
327 (*nodalCurvatures) [6 * currentVertex + 3] = fabs(lambda2) * (v2x * uP[0] + v2y * vP[0]);
328 (*nodalCurvatures) [6 * currentVertex + 4] = fabs(lambda2) * (v2x * uP[1] + v2y * vP[1]);
329 (*nodalCurvatures) [6 * currentVertex + 5] = fabs(lambda2) * (v2x * uP[2] + v2y * vP[2]);
330 }
331
332 count = 0;
333 A = 0.0;
334 B = 0.0;
335 D = 0.0;
336 computeLocalFrame (&nodeNormals[3*iVertex], uP, vP);
337 currentVertex = iVertex;
338 }
339 unitNormal2Triangle ( mesh->vertices.coord + 4*mesh->triangles.node[3*iTriangle+0],
340 mesh->vertices.coord + 4*mesh->triangles.node[3*iTriangle+1],
341 mesh->vertices.coord + 4*mesh->triangles.node[3*iTriangle+2], n, &surf);
342 double uF[3],vF[3];
343 makevector (mesh->vertices.coord + 4*mesh->triangles.node[3*iTriangle+0],
344 mesh->vertices.coord + 4*mesh->triangles.node[3*iTriangle+1],uF);
345 normalize(uF);
346 crossprod(n,uF,vF);
347 double *c = &CURV[4*iTriangle];
348 double UP[3] = {dotprod (uP,uF),dotprod (uP,vF),0};
349 normalize(UP);
350 double VP[3] = {dotprod (vP,uF),dotprod (vP,vF),0};
351 normalize(VP);
352 // printf("C[%d] = %g %g %g %g V= %g %g %g %g %g %g\n",iTriangle,c[0],c[1],c[2],c[3],uF[0],uF[1],uF[2],vF[0],vF[1],vF[2]);
353 A += (UP[0]*UP[0]*c[0] + 2*UP[0]*UP[1]*c[1] + UP[1]*UP[1]*c[3]) ;
354 D += (VP[0]*VP[0]*c[0] + 2*VP[0]*VP[1]*c[1] + VP[1]*VP[1]*c[3]) ;
355 B += (VP[0]*UP[0]*c[0] + (VP[1]*UP[0]+VP[0]*UP[1])*c[1] + VP[1]*UP[1]*c[3]) ;
356 count++;
357 }
358
359 if (debug)saveNodalField (mesh,*nodalCurvatures, 6, "curvaturesMax.pos");
360 if (debug)saveNodalField (mesh,*nodalCurvatures, -6, "curvaturesMin.pos");
361
362 //-----------------------------------------------------------------------------
363 // C R O S S F I E L D
364 //-----------------------------------------------------------------------------
365
366 // printf("coucou1 %d\n",edges->numEdges);
367 HXT_CHECK(hxtMalloc(crossField,3*edges->numEdges*sizeof(double)));
368 for (uint64_t i = 0; i<edges->numEdges; i++){
369 uint32_t v0 = edges->node[2*i ];
370 uint32_t v1 = edges->node[2*i+1];
371 double dir0[3] = {(*nodalCurvatures) [6 * v0 + 0],
372 (*nodalCurvatures) [6 * v0 + 1],
373 (*nodalCurvatures) [6 * v0 + 2]};
374 double dir1[3] = {(*nodalCurvatures) [6 * v0 + 3],
375 (*nodalCurvatures) [6 * v0 + 4],
376 (*nodalCurvatures) [6 * v0 + 5]};
377 double x1 = normalize(dir0);
378 double x2 = normalize(dir1);
379 double *n = &(nodeNormals [3 * v0]);
380 double dir2[3];
381 crossprod(dir0,dir1,dir2);
382 if (dotprod(n,dir2) < 0){
383 dir1[0] = -dir1[0];
384 dir1[1] = -dir1[1];
385 dir1[2] = -dir1[2];
386 }
387
388 // printf ("vertex %d C %g %g %g\n",v0,dir0[0],dir0[1],dir0[2]);
389 double ed[3];
390 makevector (mesh->vertices.coord + 4*v0,
391 mesh->vertices.coord + 4*v1, ed);
392 normalize(ed);
393 double SIN = dotprod(ed,dir0);
394 double COS = dotprod(ed,dir1);
395 double ANG = atan2(SIN,COS);
396 // printf ("%g %g %g\n",SIN,COS,ANG);
397 (*crossField)[2*i ] = cos(4*ANG);
398 (*crossField)[2*i+1] = sin(4*ANG);
399 if (x2 > x1){
400 double temp = x1;
401 x1 = x2;x2 = temp;
402 }
403 (*crossField)[2*edges->numEdges+i] = x1 > 1.e-4 ? x1/x2 : 1.e-12;
404 // printf("%g %g %g\n",x1,x2,(*crossField)[2*edges->numEdges+i]);
405 }
406 // printf("coucou2\n");
407
408 //-----------------------------------------------------------------------------
409 // E N D C R O S S F I E L D
410 //-----------------------------------------------------------------------------
411
412 HXT_CHECK(hxtFree(&node2tri));
413 HXT_CHECK(hxtFree(&nodeNormals));
414 HXT_CHECK(hxtFree(&CURV));
415
416 // clock_t T2 = clock();
417
418
419 // HXT_INFO ("Curvature has been computed for %u vertices in %8.4f seconds",nVertices,(double)(T2-T1)/CLOCKS_PER_SEC);
420
421 return HXT_STATUS_OK;
422
423 }
424