Subversion Repositories gelsvn

Rev

Rev 220 | Rev 357 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 220 Rev 315
1 
#ifndef __QEM_H
 1 
#ifndef __QEM_H
2 
#define __QEM_H
 2 
#define __QEM_H
3 
 
 3 
 
4 
#include <cfloat>
 4 
#include <cfloat>
5 
#include "CGLA/Vec3d.h"
 5 
#include "CGLA/Vec3d.h"
6 
#include "CGLA/Mat3x3d.h"
 6 
#include "CGLA/Mat3x3d.h"
7 
 
 7 
 
8 
 
 8 
 
9 
namespace
 9 
namespace
10 
{
 10 
{
11 
	inline const CGLA::Mat3x3d direct_product(const CGLA::Vec3d& v0, const CGLA::Vec3d& v1)
 11 
	inline const CGLA::Mat3x3d direct_product(const CGLA::Vec3d& v0, const CGLA::Vec3d& v1)
12 
		{
 12 
		{
13 
			CGLA::Mat3x3d m;
 13 
			CGLA::Mat3x3d m;
14 
			for(int i=0;i<3;++i)
 14 
			for(int i=0;i<3;++i)
15 
				for(int j=0;j<3;++j)
 15 
				for(int j=0;j<3;++j)
16 
					m[i][j] = v0[i]*v1[j];
 16 
					m[i][j] = v0[i]*v1[j];
17 
			return m;
 17 
			return m;
18 
		}
 18 
		}
19 
}
 19 
}
20 
 
 20 
 
21 
namespace Geometry
 21 
namespace Geometry
22 
{
 22 
{
23 
	class QEM
 23 
	class QEM
24 
		{
 24 
		{
25 
			CGLA::Mat3x3d A;
 25 
			CGLA::Mat3x3d A;
26 
			CGLA::Vec3d   b;
 26 
			CGLA::Vec3d   b;
27 
			double   c;
 27 
			double   c;
28 
		public:
 28 
		public:
29
29
30 
			QEM(): A(0), b(0), c(0) {}
 30 
			QEM(): A(0), b(0), c(0) {}
31
31
32 
			QEM(const CGLA::Vec3d& p0, const CGLA::Vec3d& n0, double w=1.0f):
 32 
			QEM(const CGLA::Vec3d& p0, const CGLA::Vec3d& n0, double w=1.0f):
33 
				A(direct_product(n0,n0) * w),
33 
				A(direct_product(n0,n0) * w),
34 
				b(-n0*dot(n0,p0) * w),
34 
				b(-2*n0*dot(n0,p0) * w),
35 
				c(dot(p0,n0)*dot(p0,n0) * w) {}
 35 
				c(dot(p0,n0)*dot(p0,n0) * w) {}
36 
 
 36 
 
37
37
38 
			void operator+=(const QEM& q)
 38 
			void operator+=(const QEM& q)
39 
				{
 39 
				{
40 
					A += q.A;
 40 
					A += q.A;
41 
					b += q.b;
 41 
					b += q.b;
42 
					c += q.c;
 42 
					c += q.c;
43 
				}
 43 
				}
44
44
45 
			float error(const CGLA::Vec3d& p) const
 45 
			float error(const CGLA::Vec3d& p) const
46 
				{
 46 
				{
47 
					return dot(p,A*p) + 2*dot(b,p)+ c;
 47 
					return dot(p,A*p) + dot(b,p)+ c;
48 
				}
 48 
				}
49 
 
 49 
 
50 
			double determinant() const
 50 
			double determinant() const
51 
				{
 51 
				{
52 
					return CGLA::determinant(A);
 52 
					return CGLA::determinant(A);
53 
				}
 53 
				}
54
54
55 
			const CGLA::Vec3d grad(const CGLA::Vec3d& _p) const
 55 
			const CGLA::Vec3d grad(const CGLA::Vec3d& p) const
56 
				{
 56 
				{
57 
					CGLA::Vec3d p(_p);
 - 
 
58 
					CGLA::Vec3d g(2*A*p+2*b);
 57 
					return CGLA::Vec3d(2*A*p+b);
59 
					return g;
 - 
 
60 
				}
 58 
				}
61
59
62 
			CGLA::Vec3d opt_pos(double QEM_thresh = 0.005) const;
 60 
			CGLA::Vec3d opt_pos(double QEM_thresh = 0.005) const;
63
61
64 
		};
 62 
		};
65 
}
 63 
}
66 
 
 64 
 
67 
namespace GEO = Geometry;
 65 
namespace GEO = Geometry;
68 
 
 66 
 
69 
#endif
 67 
#endif
70 
 
 68