https://github.com/ialhashim/topo-blend
Revision 38b5bd374cf0f331f148053866a9b44f60d13791 authored by jjcao1231@gmail.com on 08 January 2014, 23:53:28 UTC, committed by jjcao1231@gmail.com on 08 January 2014, 23:53:28 UTC
1 parent aa10e30
Tip revision: 38b5bd374cf0f331f148053866a9b44f60d13791 authored by jjcao1231@gmail.com on 08 January 2014, 23:53:28 UTC
Tip revision: 38b5bd3
ParametricSurface.cpp
// Geometric Tools, LLC
// Copyright (c) 1998-2012
// Distributed under the Boost Software License, Version 1.0.
// http://www.boost.org/LICENSE_1_0.txt
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
//
// File Version: 5.0.1 (2010/10/01)
#include "ParametricSurface.h"
#include "Matrix2.h"
namespace NURBS
{
//----------------------------------------------------------------------------
template <typename Real>
ParametricSurface<Real>::ParametricSurface (Real umin, Real umax,
Real vmin, Real vmax, bool rectangular)
{
assertion(umin < umax && vmin < vmax, "Invalid domain\n");
mUMin = umin;
mUMax = umax;
mVMin = vmin;
mVMax = vmax;
mRectangular = rectangular;
}
//----------------------------------------------------------------------------
template <typename Real>
Real ParametricSurface<Real>::GetUMin () const
{
return mUMin;
}
//----------------------------------------------------------------------------
template <typename Real>
Real ParametricSurface<Real>::GetUMax () const
{
return mUMax;
}
//----------------------------------------------------------------------------
template <typename Real>
Real ParametricSurface<Real>::GetVMin () const
{
return mVMin;
}
//----------------------------------------------------------------------------
template <typename Real>
Real ParametricSurface<Real>::GetVMax () const
{
return mVMax;
}
//----------------------------------------------------------------------------
template <typename Real>
bool ParametricSurface<Real>::IsRectangular () const
{
return mRectangular;
}
//----------------------------------------------------------------------------
template <typename Real>
void ParametricSurface<Real>::GetFrame (Real u, Real v,
Vector3& position, Vector3& tangent0,
Vector3& tangent1, Vector3& normal)
{
position = P(u, v);
tangent0 = PU(u, v);
tangent1 = PV(u, v);
tangent0.normalize(); // T0
tangent1.normalize(); // temporary T1 just to compute N
normal = cross(tangent0, tangent1).normalized(); // N
// The normalized first derivatives are not necessarily orthogonal.
// Recompute T1 so that {T0,T1,N} is an orthonormal set.
tangent1 = cross(normal,tangent0);
}
//----------------------------------------------------------------------------
template <typename Real>
void ParametricSurface<Real>::ComputePrincipalCurvatureInfo (Real u, Real v,
Real& curv0, Real& curv1, Vector3& dir0,
Vector3& dir1)
{
// Tangents: T0 = (x_u,y_u,z_u), T1 = (x_v,y_v,z_v)
// Normal: N = Cross(T0,T1)/Length(Cross(T0,T1))
// Metric Tensor: G = +- -+
// | Dot(T0,T0) Dot(T0,T1) |
// | Dot(T1,T0) Dot(T1,T1) |
// +- -+
//
// Curvature Tensor: B = +- -+
// | -Dot(N,T0_u) -Dot(N,T0_v) |
// | -Dot(N,T1_u) -Dot(N,T1_v) |
// +- -+
//
// Principal curvatures k are the generalized eigenvalues of
//
// Bw = kGw
//
// If k is a curvature and w=(a,b) is the corresponding solution to
// Bw = kGw, then the principal direction as a 3D vector is d = a*U+b*V.
//
// Let k1 and k2 be the principal curvatures. The mean curvature
// is (k1+k2)/2 and the Gaussian curvature is k1*k2.
// Compute derivatives.
Vector3 derU = PU(u,v);
Vector3 derV = PV(u,v);
Vector3 derUU = PUU(u,v);
Vector3 derUV = PUV(u,v);
Vector3 derVV = PVV(u,v);
// Compute the metric tensor.
Matrix2<Real> metricTensor;
metricTensor[0][0] = dot(derU,derU);
metricTensor[0][1] = dot(derU,derV);
metricTensor[1][0] = metricTensor[0][1];
metricTensor[1][1] = dot(derV,derV);
// Compute the curvature tensor.
Vector3 normal = cross(derU,derV).normalized();
Matrix2<Real> curvatureTensor;
curvatureTensor[0][0] = -dot(normal, derUU);
curvatureTensor[0][1] = -dot(normal, derUV);
curvatureTensor[1][0] = curvatureTensor[0][1];
curvatureTensor[1][1] = -dot(normal, derVV);
// Characteristic polynomial is 0 = det(B-kG) = c2*k^2+c1*k+c0.
Real c0 = curvatureTensor.Determinant();
Real c1 = ((Real)2)*curvatureTensor[0][1]* metricTensor[0][1] -
curvatureTensor[0][0]*metricTensor[1][1] -
curvatureTensor[1][1]*metricTensor[0][0];
Real c2 = metricTensor.Determinant();
// Principal curvatures are roots of characteristic polynomial.
Real temp = sqrt(abs(c1*c1 - ((Real)4)*c0*c2));
Real mult = ((Real)0.5)/c2;
curv0 = -mult*(c1+temp);
curv1 = mult*(-c1+temp);
// Principal directions are solutions to (B-kG)w = 0,
// w1 = (b12-k1*g12,-(b11-k1*g11)) OR (b22-k1*g22,-(b12-k1*g12)).
Real a0 = curvatureTensor[0][1] - curv0*metricTensor[0][1];
Real a1 = curv0*metricTensor[0][0] - curvatureTensor[0][0];
Real length = sqrt(a0*a0 + a1*a1);
if (length >= REAL_ZERO_TOLERANCE)
{
dir0 = a0*derU + a1*derV;
}
else
{
a0 = curvatureTensor[1][1] - curv0*metricTensor[1][1];
a1 = curv0*metricTensor[0][1] - curvatureTensor[0][1];
length = sqrt(a0*a0 + a1*a1);
if (length >= REAL_ZERO_TOLERANCE)
{
dir0 = a0*derU + a1*derV;
}
else
{
// Umbilic (surface is locally sphere, any direction principal).
dir0 = derU;
}
}
dir0.normalize();
// Second tangent is cross product of first tangent and normal.
dir1 = cross(dir0, normal);
}
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
// Explicit instantiation.
//----------------------------------------------------------------------------
//template
//class ParametricSurface<float>;
template
class ParametricSurface<double>;
//----------------------------------------------------------------------------
}
Computing file changes ...