1 /*
2 Szymon Rusinkiewicz
3 Princeton University
4
5 TriMesh_curvature.cc
6 Computation of per-vertex principal curvatures and directions.
7
8 Uses algorithm from
9 Rusinkiewicz, Szymon.
10 "Estimating Curvatures and Their Derivatives on Triangle Meshes,"
11 Proc. 3DPVT, 2004.
12 */
13
14 #include "TriMesh.h"
15 #include "TriMesh_algo.h"
16 #include "lineqn.h"
17 using namespace std;
18
19
20 // i+1 and i-1 modulo 3
21 // This way of computing it tends to be faster than using %
22 #define NEXT(i) ((i)<2 ? (i)+1 : (i)-2)
23 #define PREV(i) ((i)>0 ? (i)-1 : (i)+2)
24
25
26 namespace trimesh {
27
28 // Rotate a coordinate system to be perpendicular to the given normal
rot_coord_sys(const vec & old_u,const vec & old_v,const vec & new_norm,vec & new_u,vec & new_v)29 static void rot_coord_sys(const vec &old_u, const vec &old_v,
30 const vec &new_norm,
31 vec &new_u, vec &new_v)
32 {
33 new_u = old_u;
34 new_v = old_v;
35 vec old_norm = old_u CROSS old_v;
36 float ndot = old_norm DOT new_norm;
37 if (unlikely(ndot <= -1.0f)) {
38 new_u = -new_u;
39 new_v = -new_v;
40 return;
41 }
42 vec perp_old = new_norm - ndot * old_norm;
43 vec dperp = 1.0f / (1 + ndot) * (old_norm + new_norm);
44 new_u -= dperp * (new_u DOT perp_old);
45 new_v -= dperp * (new_v DOT perp_old);
46 }
47
48
49 // Reproject a curvature tensor from the basis spanned by old_u and old_v
50 // (which are assumed to be unit-length and perpendicular) to the
51 // new_u, new_v basis.
proj_curv(const vec & old_u,const vec & old_v,float old_ku,float old_kuv,float old_kv,const vec & new_u,const vec & new_v,float & new_ku,float & new_kuv,float & new_kv)52 void proj_curv(const vec &old_u, const vec &old_v,
53 float old_ku, float old_kuv, float old_kv,
54 const vec &new_u, const vec &new_v,
55 float &new_ku, float &new_kuv, float &new_kv)
56 {
57 vec r_new_u, r_new_v;
58 rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v);
59
60 float u1 = r_new_u DOT old_u;
61 float v1 = r_new_u DOT old_v;
62 float u2 = r_new_v DOT old_u;
63 float v2 = r_new_v DOT old_v;
64 new_ku = old_ku * u1*u1 + old_kuv * (2.0f * u1*v1) + old_kv * v1*v1;
65 new_kuv = old_ku * u1*u2 + old_kuv * (u1*v2 + u2*v1) + old_kv * v1*v2;
66 new_kv = old_ku * u2*u2 + old_kuv * (2.0f * u2*v2) + old_kv * v2*v2;
67 }
68
69
70 // Like the above, but for dcurv
proj_dcurv(const vec & old_u,const vec & old_v,const Vec<4> old_dcurv,const vec & new_u,const vec & new_v,Vec<4> & new_dcurv)71 void proj_dcurv(const vec &old_u, const vec &old_v,
72 const Vec<4> old_dcurv,
73 const vec &new_u, const vec &new_v,
74 Vec<4> &new_dcurv)
75 {
76 vec r_new_u, r_new_v;
77 rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v);
78
79 float u1 = r_new_u DOT old_u;
80 float v1 = r_new_u DOT old_v;
81 float u2 = r_new_v DOT old_u;
82 float v2 = r_new_v DOT old_v;
83
84 new_dcurv[0] = old_dcurv[0]*u1*u1*u1 +
85 old_dcurv[1]*3.0f*u1*u1*v1 +
86 old_dcurv[2]*3.0f*u1*v1*v1 +
87 old_dcurv[3]*v1*v1*v1;
88 new_dcurv[1] = old_dcurv[0]*u1*u1*u2 +
89 old_dcurv[1]*(u1*u1*v2 + 2.0f*u2*u1*v1) +
90 old_dcurv[2]*(u2*v1*v1 + 2.0f*u1*v1*v2) +
91 old_dcurv[3]*v1*v1*v2;
92 new_dcurv[2] = old_dcurv[0]*u1*u2*u2 +
93 old_dcurv[1]*(u2*u2*v1 + 2.0f*u1*u2*v2) +
94 old_dcurv[2]*(u1*v2*v2 + 2.0f*u2*v2*v1) +
95 old_dcurv[3]*v1*v2*v2;
96 new_dcurv[3] = old_dcurv[0]*u2*u2*u2 +
97 old_dcurv[1]*3.0f*u2*u2*v2 +
98 old_dcurv[2]*3.0f*u2*v2*v2 +
99 old_dcurv[3]*v2*v2*v2;
100 }
101
102
103 // Given a curvature tensor, find principal directions and curvatures
104 // Makes sure that pdir1 and pdir2 are perpendicular to normal
diagonalize_curv(const vec & old_u,const vec & old_v,float ku,float kuv,float kv,const vec & new_norm,vec & pdir1,vec & pdir2,float & k1,float & k2)105 void diagonalize_curv(const vec &old_u, const vec &old_v,
106 float ku, float kuv, float kv,
107 const vec &new_norm,
108 vec &pdir1, vec &pdir2, float &k1, float &k2)
109 {
110 vec r_old_u, r_old_v;
111 rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);
112
113 float c = 1, s = 0, tt = 0;
114 if (likely(kuv != 0.0f)) {
115 // Jacobi rotation to diagonalize
116 float h = 0.5f * (kv - ku) / kuv;
117 tt = (h < 0.0f) ?
118 1.0f / (h - sqrt(1.0f + h*h)) :
119 1.0f / (h + sqrt(1.0f + h*h));
120 c = 1.0f / sqrt(1.0f + tt*tt);
121 s = tt * c;
122 }
123
124 k1 = ku - tt * kuv;
125 k2 = kv + tt * kuv;
126
127 if (fabs(k1) >= fabs(k2)) {
128 pdir1 = c*r_old_u - s*r_old_v;
129 } else {
130 swap(k1, k2);
131 pdir1 = s*r_old_u + c*r_old_v;
132 }
133 pdir2 = new_norm CROSS pdir1;
134 }
135
136
137 // Compute principal curvatures and directions.
need_curvatures()138 void TriMesh::need_curvatures()
139 {
140 if (curv1.size() == vertices.size())
141 return;
142 need_faces();
143 need_normals();
144 need_pointareas();
145
146 dprintf("Computing curvatures... ");
147
148 // Resize the arrays we'll be using
149 int nv = vertices.size(), nf = faces.size();
150 curv1.clear(); curv1.resize(nv); curv2.clear(); curv2.resize(nv);
151 pdir1.clear(); pdir1.resize(nv); pdir2.clear(); pdir2.resize(nv);
152 vector<float> curv12(nv);
153
154 // Set up an initial coordinate system per vertex
155 for (int i = 0; i < nf; i++) {
156 pdir1[faces[i][0]] = vertices[faces[i][1]] -
157 vertices[faces[i][0]];
158 pdir1[faces[i][1]] = vertices[faces[i][2]] -
159 vertices[faces[i][1]];
160 pdir1[faces[i][2]] = vertices[faces[i][0]] -
161 vertices[faces[i][2]];
162 }
163 #pragma omp parallel for
164 for (int i = 0; i < nv; i++) {
165 pdir1[i] = pdir1[i] CROSS normals[i];
166 normalize(pdir1[i]);
167 pdir2[i] = normals[i] CROSS pdir1[i];
168 }
169
170 // Compute curvature per-face
171 #pragma omp parallel for
172 for (int i = 0; i < nf; i++) {
173 // Edges
174 vec e[3] = { vertices[faces[i][2]] - vertices[faces[i][1]],
175 vertices[faces[i][0]] - vertices[faces[i][2]],
176 vertices[faces[i][1]] - vertices[faces[i][0]] };
177
178 // N-T-B coordinate system per face
179 vec t = e[0];
180 normalize(t);
181 vec n = e[0] CROSS e[1];
182 vec b = n CROSS t;
183 normalize(b);
184
185 // Estimate curvature based on variation of normals
186 // along edges
187 float m[3] = { 0, 0, 0 };
188 float w[3][3] = { {0,0,0}, {0,0,0}, {0,0,0} };
189 for (int j = 0; j < 3; j++) {
190 float u = e[j] DOT t;
191 float v = e[j] DOT b;
192 w[0][0] += u*u;
193 w[0][1] += u*v;
194 //w[1][1] += v*v + u*u;
195 //w[1][2] += u*v;
196 w[2][2] += v*v;
197 vec dn = normals[faces[i][PREV(j)]] -
198 normals[faces[i][NEXT(j)]];
199 float dnu = dn DOT t;
200 float dnv = dn DOT b;
201 m[0] += dnu*u;
202 m[1] += dnu*v + dnv*u;
203 m[2] += dnv*v;
204 }
205 w[1][1] = w[0][0] + w[2][2];
206 w[1][2] = w[0][1];
207
208 // Least squares solution
209 float diag[3];
210 if (!ldltdc<float,3>(w, diag)) {
211 //dprintf("ldltdc failed!\n");
212 continue;
213 }
214 ldltsl<float,3>(w, diag, m, m);
215
216 // Push it back out to the vertices
217 for (int j = 0; j < 3; j++) {
218 int vj = faces[i][j];
219 float c1, c12, c2;
220 proj_curv(t, b, m[0], m[1], m[2],
221 pdir1[vj], pdir2[vj], c1, c12, c2);
222 float wt = cornerareas[i][j] / pointareas[vj];
223 #pragma omp atomic
224 curv1[vj] += wt * c1;
225 #pragma omp atomic
226 curv12[vj] += wt * c12;
227 #pragma omp atomic
228 curv2[vj] += wt * c2;
229 }
230 }
231
232 // Compute principal directions and curvatures at each vertex
233 #pragma omp parallel for
234 for (int i = 0; i < nv; i++) {
235 diagonalize_curv(pdir1[i], pdir2[i],
236 curv1[i], curv12[i], curv2[i],
237 normals[i], pdir1[i], pdir2[i],
238 curv1[i], curv2[i]);
239 }
240 dprintf("Done.\n");
241 }
242
243
244 // Compute derivatives of curvature.
need_dcurv()245 void TriMesh::need_dcurv()
246 {
247 if (dcurv.size() == vertices.size())
248 return;
249 need_curvatures();
250
251 dprintf("Computing dcurv... ");
252
253 // Resize the arrays we'll be using
254 int nv = vertices.size(), nf = faces.size();
255 dcurv.clear(); dcurv.resize(nv);
256
257 // Compute dcurv per-face
258 #pragma omp parallel for
259 for (int i = 0; i < nf; i++) {
260 // Edges
261 vec e[3] = { vertices[faces[i][2]] - vertices[faces[i][1]],
262 vertices[faces[i][0]] - vertices[faces[i][2]],
263 vertices[faces[i][1]] - vertices[faces[i][0]] };
264
265 // N-T-B coordinate system per face
266 vec t = e[0];
267 normalize(t);
268 vec n = e[0] CROSS e[1];
269 vec b = n CROSS t;
270 normalize(b);
271
272 // Project curvature tensor from each vertex into this
273 // face's coordinate system
274 vec fcurv[3];
275 for (int j = 0; j < 3; j++) {
276 int vj = faces[i][j];
277 proj_curv(pdir1[vj], pdir2[vj], curv1[vj], 0, curv2[vj],
278 t, b, fcurv[j][0], fcurv[j][1], fcurv[j][2]);
279
280 }
281
282 // Estimate dcurv based on variation of curvature along edges
283 float m[4] = { 0, 0, 0, 0 };
284 float w[4][4] = { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} };
285 for (int j = 0; j < 3; j++) {
286 // Variation of curvature along each edge
287 vec dfcurv = fcurv[PREV(j)] - fcurv[NEXT(j)];
288 float u = e[j] DOT t;
289 float v = e[j] DOT b;
290 float u2 = u*u, v2 = v*v, uv = u*v;
291 w[0][0] += u2;
292 w[0][1] += uv;
293 //w[1][1] += 2.0f*u2 + v2;
294 //w[1][2] += 2.0f*uv;
295 //w[2][2] += u2 + 2.0f*v2;
296 //w[2][3] += uv;
297 w[3][3] += v2;
298 m[0] += u*dfcurv[0];
299 m[1] += v*dfcurv[0] + 2.0f*u*dfcurv[1];
300 m[2] += 2.0f*v*dfcurv[1] + u*dfcurv[2];
301 m[3] += v*dfcurv[2];
302 }
303 w[1][1] = 2.0f * w[0][0] + w[3][3];
304 w[1][2] = 2.0f * w[0][1];
305 w[2][2] = w[0][0] + 2.0f * w[3][3];
306 w[2][3] = w[0][1];
307
308 // Least squares solution
309 float d[4];
310 if (!ldltdc<float,4>(w, d)) {
311 //dprintf("ldltdc failed!\n");
312 continue;
313 }
314 ldltsl<float,4>(w, d, m, m);
315 Vec<4> face_dcurv(m);
316
317 // Push it back out to each vertex
318 for (int j = 0; j < 3; j++) {
319 int vj = faces[i][j];
320 Vec<4> this_vert_dcurv;
321 proj_dcurv(t, b, face_dcurv,
322 pdir1[vj], pdir2[vj], this_vert_dcurv);
323 float wt = cornerareas[i][j] / pointareas[vj];
324 dcurv[vj] += wt * this_vert_dcurv;
325 }
326 }
327
328 dprintf("Done.\n");
329 }
330
331 }; // namespace trimesh
332