Subversion Repositories gelsvn

Rev

Rev 632 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 632 Rev 646
Line 4... Line 4...
4 
 * For license and list of authors, see ../../doc/intro.pdf
4 
 * For license and list of authors, see ../../doc/intro.pdf
5 
 * ----------------------------------------------------------------------- */
 5 
 * ----------------------------------------------------------------------- */
6 
 
 6 
 
7 
#include "curvature.h"
 7 
#include "curvature.h"
8 
 
 8 
 
9 
#include "../GLGraphics/gel_glut.h"
 9 
#include "../GLGraphics/gel_glu.h"
10 
 
 10 
 
11 
#include <iostream>
 11 
#include <iostream>
12 
#include "../CGLA/eigensolution.h"
 12 
#include "../CGLA/eigensolution.h"
13 
#include "../CGLA/Vec2d.h"
 13 
#include "../CGLA/Vec2d.h"
14 
#include "../CGLA/Vec3d.h"
 14 
#include "../CGLA/Vec3d.h"
Line 246... Line 246...
246 
        x = LinearLSSolve(PMat,b);
 246 
        x = LinearLSSolve(PMat,b);
247 
        if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;
 247 
        if(isnan(x[0])) cout << __LINE__ << " " << PMat << b << endl ;
248 
 
 248 
 
249 
        // Finally compute the shape tensor from the coefficients
249 
        // Finally compute the shape tensor from the coefficients
250 
        // using the first and second fundamental forms.
250 
        // using the first and second fundamental forms.
251 
        curv_tensor =Mat2x2d(x[0],x[1],x[1],x[2]);
 251 
        curv_tensor = - Mat2x2d(x[0],x[1],x[1],x[2]);
252 
 
 252 
 
253 
    }
 253 
    }
254 
 
 254 
 
255 
    void curvature_tensors_from_edges(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)
 255 
    void curvature_tensors_from_edges(const Manifold& m, VertexAttributeVector<Mat3x3d>& curvature_tensors)
256 
    {
 256 
    {
Line 303... Line 303...
303 
 
 303 
 
304 
 
 304 
 
305 
    void curvature_paraboloids( const Manifold& m,
305 
    void curvature_paraboloids( const Manifold& m,
306 
                                VertexAttributeVector<Vec3d>& min_curv_direction,
306 
                                VertexAttributeVector<Vec3d>& min_curv_direction,
307 
                                VertexAttributeVector<Vec3d>& max_curv_direction,
 307 
                                VertexAttributeVector<Vec3d>& max_curv_direction,
308 
                                VertexAttributeVector<double>& curvature)
 308 
                                VertexAttributeVector<Vec2d>& curvature)
309 
    {
 309 
    {
310 
 
 310 
 
311 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 311 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
312 
            Mat2x2d tensor;
 312 
            Mat2x2d tensor;
313 
            Mat3x3d frame;
 313 
            Mat3x3d frame;
Line 320... Line 320...
320 
                cout << tensor << Q << L << endl;
 320 
                cout << tensor << Q << L << endl;
321 
 
 321 
 
322 
            int max_idx = 0;
 322 
            int max_idx = 0;
323 
            int min_idx = 1;
 323 
            int min_idx = 1;
324 
 
 324 
 
325 
            if(fabs(L[max_idx][max_idx])<fabs(L[min_idx][min_idx])) swap(max_idx, min_idx);
 325 
            if(L[max_idx][max_idx]<L[min_idx][min_idx]) swap(max_idx, min_idx);
326 
 
 326 
 
327 
            Mat3x3d frame_t = transpose(frame);
 327 
            Mat3x3d frame_t = transpose(frame);
328 
 
 328 
 
329 
            max_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0));
 329 
            max_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[max_idx][0], Q[max_idx][1], 0));
330 
 
 330 
 
331 
            min_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0));
 331 
            min_curv_direction[*v] = cond_normalize(frame_t * Vec3d(Q[min_idx][0], Q[min_idx][1], 0));
332 
 
 332 
 
- 
 
 333 
            curvature[*v][0] = L[min_idx][min_idx];
333 
            curvature[*v] = L[0][0]*L[1][1];
 334 
            curvature[*v][1] = L[max_idx][max_idx];
334 
        }
 335 
        }
335 
    }
 336 
    }
336 
 
 337 
 
337 
 
 338 
 
338 
    void curvature_from_tensors(const Manifold& m,
 339 
    void curvature_from_tensors(const Manifold& m,
339 
                                const VertexAttributeVector<Mat3x3d>& curvature_tensors,
 340 
                                const VertexAttributeVector<Mat3x3d>& curvature_tensors,
340 
                                VertexAttributeVector<Vec3d>& min_curv_direction,
 341 
                                VertexAttributeVector<Vec3d>& min_curv_direction,
341 
                                VertexAttributeVector<Vec3d>& max_curv_direction,
 342 
                                VertexAttributeVector<Vec3d>& max_curv_direction,
342 
                                VertexAttributeVector<double>& curvature)
 343 
                                VertexAttributeVector<Vec2d>& curvature)
343 
    {
 344 
    {
344 
        assert(curvature_tensors.size() == m.allocated_vertices());
 345 
        assert(curvature_tensors.size() == m.allocated_vertices());
345 
 
 346 
 
346 
        double max_val = -1e30;
 347 
        double max_val = -1e30;
347 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 348 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
Line 351... Line 352...
351 
            Vec3d dmin, dmax;
 352 
            Vec3d dmin, dmax;
352 
            if(s == 0)
 353 
            if(s == 0)
353 
            {
 354 
            {
354 
                Vec3d n(normal(m, *v));
 355 
                Vec3d n(normal(m, *v));
355 
                orthogonal(n, dmin, dmax);
 356 
                orthogonal(n, dmin, dmax);
356 
                curvature[*v] = 0;
 357 
                curvature[*v] = Vec2d(0);
357 
                cout << " rank 0 " << endl;
 358 
                cout << " rank 0 " << endl;
358 
            }
 359 
            }
359 
            else if(s == 1)
 360 
            else if(s == 1)
360 
            {
 361 
            {
361 
                Vec3d n(normal(m, *v));
 362 
                Vec3d n(normal(m, *v));
362 
                dmin = normalize(Q[0]);
 363 
                dmin = normalize(Q[0]);
363 
                dmax = cross(n, dmin);
 364 
                dmax = cross(n, dmin);
364 
                curvature[*v] = 0;
 365 
                curvature[*v] = Vec2d(0);
365 
                cout << " rank 1 " << endl;
 366 
                cout << " rank 1 " << endl;
366 
            }
 367 
            }
367 
            else
 368 
            else
368 
            {
 369 
            {
369 
                /*				Vec3d l(fabs(L[0][0]), fabs(L[1][1]), fabs(L[2][2]));
370 
                Vec2d l(fabs(L[0][0]), fabs(L[1][1]));
370
371 
 
371 
                int z_idx=2;
- 
 
372 
                if(s==3)
- 
 
373 
                {
- 
 
374 
                if(l[0] < l[1])
- 
 
375 
                z_idx = l[0]<l[2] ? 0 : 2;
- 
 
376 
                else
- 
 
377 
                z_idx = l[1]<l[2] ? 1 : 2;
- 
 
378 
                }
- 
 
379 
                int max_idx = (z_idx + 1) % 3;
- 
 
380 
                int min_idx = (z_idx + 2) % 3;
- 
 
381
- 
 
382 
                if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);
- 
 
383 
                */
 - 
 
384 
                int max_idx = 0;
 372 
                int max_idx = 0;
385 
                int min_idx = 1;
 373 
                int min_idx = 1;
- 
 
 374 
 
- 
 
 375 
                if(l[max_idx] < l[min_idx]) swap(max_idx, min_idx);
- 
 
 376
386 
                // Yes - the biggest eigenvalue corresponds to the min direction
377 
                // Yes - the biggest eigenvalue corresponds to the min direction
387 
                // and vice versa.
378 
                // and vice versa.
388 
                dmin = normalize(Q[max_idx]);
 379 
                dmin = normalize(Q[max_idx]);
389 
                dmax = normalize(Q[min_idx]);
 380 
                dmax = normalize(Q[min_idx]);
390 
 
 381 
 
391 
                curvature[*v] = L[max_idx][max_idx]*L[min_idx][min_idx];
 382 
                curvature[*v][0] = L[min_idx][min_idx];
- 
 
 383 
                curvature[*v][1] = L[max_idx][max_idx];
392 
 
 384 
 
393 
            }
 385 
            }
394 
            min_curv_direction[*v] = dmin;
 386 
            min_curv_direction[*v] = dmin;
395 
            max_curv_direction[*v] = dmax;
 387 
            max_curv_direction[*v] = dmax;
396 
            max_val = max(fabs(curvature[*v]), max_val);
 388 
            max_val = max(fabs(curvature[*v][1]), max_val);
397 
 
 389 
 
398 
        }
 390 
        }
399 
        //scal = 1.0/max_val;
391 
        //scal = 1.0/max_val;
400 
    }
 392 
    }
401 
}
 393 
}