mirror of
https://github.com/lucaspalomodevelop/meshlab.git
synced 2026-03-20 11:26:11 +00:00
updatedupdated filter_mls to float/double independence
This commit is contained in:
Paolo Cignoni cignoni
parent
e9b370122a
commit
db7ea0cf7c
@ -30,472 +30,472 @@ namespace GaelMls {
|
||||
template<typename _MeshType>
|
||||
void APSS<_MeshType>::setSphericalParameter(Scalar v)
|
||||
{
|
||||
mSphericalParameter = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mSphericalParameter = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename APSS<_MeshType>::Scalar APSS<_MeshType>::potential(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
|
||||
if (mStatus==ASS_SPHERE)
|
||||
{
|
||||
Scalar aux = vcg::Norm(lx - mCenter) - mRadius;
|
||||
if (uQuad<0.)
|
||||
aux = -aux;
|
||||
return aux;
|
||||
}
|
||||
else if (mStatus==ASS_PLANE)
|
||||
return vcg::Dot(lx,uLinear) + uConstant;
|
||||
else
|
||||
{
|
||||
return uConstant + vcg::Dot(lx,uLinear) + uQuad * vcg::SquaredNorm(lx);
|
||||
}
|
||||
if (mStatus==ASS_SPHERE)
|
||||
{
|
||||
Scalar aux = vcg::Norm(lx - mCenter) - mRadius;
|
||||
if (uQuad<0.)
|
||||
aux = -aux;
|
||||
return aux;
|
||||
}
|
||||
else if (mStatus==ASS_PLANE)
|
||||
return (lx*uLinear) + uConstant;
|
||||
else
|
||||
{
|
||||
return uConstant + (lx*uLinear) + uQuad * vcg::SquaredNorm(lx);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename APSS<_MeshType>::Scalar APSS<_MeshType>::approxMeanCurvature(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
|
||||
if (mStatus==ASS_SPHERE)
|
||||
return (uQuad>0.?1.0:-1.0)/mRadius;
|
||||
else
|
||||
return 0.;
|
||||
if (mStatus==ASS_SPHERE)
|
||||
return (uQuad>0.?1.0:-1.0)/mRadius;
|
||||
else
|
||||
return 0.;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename APSS<_MeshType>::VectorType APSS<_MeshType>::gradient(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_OK;
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return VectorType(0,0,0);
|
||||
}
|
||||
}
|
||||
if (errorMask)
|
||||
*errorMask = MLS_OK;
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return VectorType(0,0,0);
|
||||
}
|
||||
}
|
||||
|
||||
if (mGradientHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
VectorType grad;
|
||||
mlsGradient(x,grad);
|
||||
return grad;
|
||||
}
|
||||
else
|
||||
{
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
if (mStatus==ASS_PLANE)
|
||||
return VectorType(uLinear.X(), uLinear.Y(), uLinear.Z());
|
||||
else
|
||||
{
|
||||
LVector g = uLinear + lx * (Scalar(2) * uQuad);
|
||||
return VectorType(g.X(), g.Y(), g.Z());
|
||||
}
|
||||
}
|
||||
if (mGradientHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
VectorType grad;
|
||||
mlsGradient(x,grad);
|
||||
return grad;
|
||||
}
|
||||
else
|
||||
{
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
if (mStatus==ASS_PLANE)
|
||||
return VectorType(uLinear.X(), uLinear.Y(), uLinear.Z());
|
||||
else
|
||||
{
|
||||
LVector g = uLinear + lx * (Scalar(2) * uQuad);
|
||||
return VectorType(g.X(), g.Y(), g.Z());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename APSS<_MeshType>::MatrixType APSS<_MeshType>::hessian(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return MatrixType();
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!fit(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return MatrixType();
|
||||
}
|
||||
}
|
||||
|
||||
MatrixType hessian;
|
||||
if (Base::mHessianHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
mlsHessian(x, hessian);
|
||||
}
|
||||
else
|
||||
{
|
||||
// this is very approximate !!
|
||||
Scalar c = Scalar(2) * uQuad;
|
||||
for (int i=0; i<3; ++i)
|
||||
for (int j=0; j<3; ++j)
|
||||
{
|
||||
if (i==j)
|
||||
hessian[i][j] = c;
|
||||
else
|
||||
hessian[i][j] = 0;
|
||||
}
|
||||
}
|
||||
return hessian;
|
||||
MatrixType hessian;
|
||||
if (Base::mHessianHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
mlsHessian(x, hessian);
|
||||
}
|
||||
else
|
||||
{
|
||||
// this is very approximate !!
|
||||
Scalar c = Scalar(2) * uQuad;
|
||||
for (int i=0; i<3; ++i)
|
||||
for (int j=0; j<3; ++j)
|
||||
{
|
||||
if (i==j)
|
||||
hessian[i][j] = c;
|
||||
else
|
||||
hessian[i][j] = 0;
|
||||
}
|
||||
}
|
||||
return hessian;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename APSS<_MeshType>::VectorType APSS<_MeshType>::project(const VectorType& x, VectorType* pNormal, int* errorMask) const
|
||||
{
|
||||
int iterationCount = 0;
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
LVector position = lx;
|
||||
LVector normal;
|
||||
LVector previousPosition;
|
||||
LScalar delta2;
|
||||
LScalar epsilon2 = mAveragePointSpacing * mProjectionAccuracy;
|
||||
epsilon2 = epsilon2 * epsilon2;
|
||||
do {
|
||||
if (!fit(VectorType(position.X(), position.Y(), position.Z())))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
//std::cout << " proj failed\n";
|
||||
return x;
|
||||
}
|
||||
int iterationCount = 0;
|
||||
LVector lx(x.X(), x.Y(), x.Z());
|
||||
LVector position = lx;
|
||||
LVector normal;
|
||||
LVector previousPosition;
|
||||
LScalar delta2;
|
||||
LScalar epsilon2 = mAveragePointSpacing * mProjectionAccuracy;
|
||||
epsilon2 = epsilon2 * epsilon2;
|
||||
do {
|
||||
if (!fit(VectorType(position.X(), position.Y(), position.Z())))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
//std::cout << " proj failed\n";
|
||||
return x;
|
||||
}
|
||||
|
||||
previousPosition = position;
|
||||
// local projection
|
||||
if (mStatus==ASS_SPHERE)
|
||||
{
|
||||
normal = lx - mCenter;
|
||||
normal.Normalize();
|
||||
position = mCenter + normal * mRadius;
|
||||
previousPosition = position;
|
||||
// local projection
|
||||
if (mStatus==ASS_SPHERE)
|
||||
{
|
||||
normal = lx - mCenter;
|
||||
normal.Normalize();
|
||||
position = mCenter + normal * mRadius;
|
||||
|
||||
normal = uLinear + position * (LScalar(2) * uQuad);
|
||||
normal.Normalize();
|
||||
}
|
||||
else if (mStatus==ASS_PLANE)
|
||||
{
|
||||
normal = uLinear;
|
||||
position = lx - uLinear * (vcg::Dot(lx,uLinear) + uConstant);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Newton iterations
|
||||
LVector grad;
|
||||
LVector dir = uLinear+lx*(2.*uQuad);
|
||||
LScalar ilg = 1./vcg::Norm(dir);
|
||||
dir *= ilg;
|
||||
LScalar ad = uConstant + vcg::Dot(uLinear,lx) + uQuad * vcg::SquaredNorm(lx);
|
||||
LScalar delta = -ad*std::min<Scalar>(ilg,1.);
|
||||
LVector p = lx + dir*delta;
|
||||
for (int i=0 ; i<2 ; ++i)
|
||||
{
|
||||
grad = uLinear+p*(2.*uQuad);
|
||||
ilg = 1./vcg::Norm(grad);
|
||||
delta = -(uConstant + vcg::Dot(uLinear,p) + uQuad * vcg::SquaredNorm(p))*std::min<Scalar>(ilg,1.);
|
||||
p += dir*delta;
|
||||
}
|
||||
position = p;
|
||||
normal = uLinear + position * (LScalar(2) * uQuad);
|
||||
normal.Normalize();
|
||||
}
|
||||
else if (mStatus==ASS_PLANE)
|
||||
{
|
||||
normal = uLinear;
|
||||
position = lx - uLinear * ((lx*uLinear) + uConstant);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Newton iterations
|
||||
LVector grad;
|
||||
LVector dir = uLinear+lx*(2.*uQuad);
|
||||
LScalar ilg = 1./vcg::Norm(dir);
|
||||
dir *= ilg;
|
||||
LScalar ad = uConstant + (uLinear*lx) + uQuad * vcg::SquaredNorm(lx);
|
||||
LScalar delta = -ad*std::min<Scalar>(ilg,1.);
|
||||
LVector p = lx + dir*delta;
|
||||
for (int i=0 ; i<2 ; ++i)
|
||||
{
|
||||
grad = uLinear+p*(2.*uQuad);
|
||||
ilg = 1./vcg::Norm(grad);
|
||||
delta = -(uConstant + (uLinear*p) + uQuad * vcg::SquaredNorm(p))*std::min<Scalar>(ilg,1.);
|
||||
p += dir*delta;
|
||||
}
|
||||
position = p;
|
||||
|
||||
normal = uLinear + position * (Scalar(2) * uQuad);
|
||||
normal.Normalize();
|
||||
}
|
||||
normal = uLinear + position * (Scalar(2) * uQuad);
|
||||
normal.Normalize();
|
||||
}
|
||||
|
||||
delta2 = vcg::SquaredNorm(previousPosition - position);
|
||||
} while ( delta2>epsilon2 && ++iterationCount<mMaxNofProjectionIterations);
|
||||
delta2 = vcg::SquaredNorm(previousPosition - position);
|
||||
} while ( delta2>epsilon2 && ++iterationCount<mMaxNofProjectionIterations);
|
||||
|
||||
if (pNormal)
|
||||
{
|
||||
if (mGradientHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
VectorType grad;
|
||||
mlsGradient(vcg::Point3<Scalar>::Construct(position), grad);
|
||||
grad.Normalize();
|
||||
*pNormal = grad;
|
||||
}
|
||||
else
|
||||
{
|
||||
for (int k=0; k<3; ++k)
|
||||
(*pNormal)[k] = normal[k];
|
||||
}
|
||||
}
|
||||
if (pNormal)
|
||||
{
|
||||
if (mGradientHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
VectorType grad;
|
||||
mlsGradient(vcg::Point3<Scalar>::Construct(position), grad);
|
||||
grad.Normalize();
|
||||
*pNormal = grad;
|
||||
}
|
||||
else
|
||||
{
|
||||
for (int k=0; k<3; ++k)
|
||||
(*pNormal)[k] = normal[k];
|
||||
}
|
||||
}
|
||||
|
||||
if (iterationCount>=mMaxNofProjectionIterations && errorMask)
|
||||
*errorMask = MLS_TOO_MANY_ITERS;
|
||||
if (iterationCount>=mMaxNofProjectionIterations && errorMask)
|
||||
*errorMask = MLS_TOO_MANY_ITERS;
|
||||
|
||||
return VectorType(position.X(), position.Y(), position.Z());
|
||||
return VectorType(position.X(), position.Y(), position.Z());
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool APSS<_MeshType>::fit(const VectorType& x) const
|
||||
{
|
||||
Base::computeNeighborhood(x, true);
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
Base::computeNeighborhood(x, true);
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
|
||||
if (nofSamples==0)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
else if (nofSamples==1)
|
||||
{
|
||||
int id = mNeighborhood.index(0);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
if (nofSamples==0)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
else if (nofSamples==1)
|
||||
{
|
||||
int id = mNeighborhood.index(0);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
|
||||
uLinear = n;
|
||||
uConstant = -vcg::Dot(p, uLinear);
|
||||
uQuad = 0;
|
||||
mStatus = ASS_PLANE;
|
||||
return true;
|
||||
}
|
||||
uLinear = n;
|
||||
uConstant = -(p*uLinear);
|
||||
uQuad = 0;
|
||||
mStatus = ASS_PLANE;
|
||||
return true;
|
||||
}
|
||||
|
||||
LVector sumP; sumP.SetZero();
|
||||
LVector sumN; sumN.SetZero();
|
||||
LScalar sumDotPN = 0.;
|
||||
LScalar sumDotPP = 0.;
|
||||
LScalar sumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar w = mCachedWeights.at(i);
|
||||
LVector sumP; sumP.SetZero();
|
||||
LVector sumN; sumN.SetZero();
|
||||
LScalar sumDotPN = 0.;
|
||||
LScalar sumDotPP = 0.;
|
||||
LScalar sumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar w = mCachedWeights.at(i);
|
||||
|
||||
sumP += p * w;
|
||||
sumN += n * w;
|
||||
sumDotPN += w * vcg::Dot(n,p);
|
||||
sumDotPP += w * vcg::SquaredNorm(p);
|
||||
sumW += w;
|
||||
}
|
||||
sumP += p * w;
|
||||
sumN += n * w;
|
||||
sumDotPN += w * (n*p);
|
||||
sumDotPP += w * vcg::SquaredNorm(p);
|
||||
sumW += w;
|
||||
}
|
||||
|
||||
LScalar invSumW = Scalar(1)/sumW;
|
||||
LScalar aux4 = mSphericalParameter * LScalar(0.5) *
|
||||
(sumDotPN - invSumW*vcg::Dot(sumP,sumN))
|
||||
/(sumDotPP - invSumW*vcg::SquaredNorm(sumP));
|
||||
uLinear = (sumN-sumP*(Scalar(2)*aux4))*invSumW;
|
||||
uConstant = -invSumW*(Dot(uLinear,sumP) + sumDotPP*aux4);
|
||||
uQuad = aux4;
|
||||
LScalar invSumW = Scalar(1)/sumW;
|
||||
LScalar aux4 = mSphericalParameter * LScalar(0.5) *
|
||||
(sumDotPN - invSumW*(sumP*sumN))
|
||||
/(sumDotPP - invSumW*vcg::SquaredNorm(sumP));
|
||||
uLinear = (sumN-sumP*(Scalar(2)*aux4))*invSumW;
|
||||
uConstant = -invSumW*((uLinear*sumP) + sumDotPP*aux4);
|
||||
uQuad = aux4;
|
||||
|
||||
// finalize
|
||||
if (fabs(uQuad)>1e-7)
|
||||
{
|
||||
mStatus = ASS_SPHERE;
|
||||
LScalar b = 1./uQuad;
|
||||
mCenter = uLinear*(-0.5*b);
|
||||
mRadius = sqrt( vcg::SquaredNorm(mCenter) - b*uConstant );
|
||||
}
|
||||
else if (uQuad==0.)
|
||||
{
|
||||
mStatus = ASS_PLANE;
|
||||
LScalar s = LScalar(1)/vcg::Norm(uLinear);
|
||||
assert(!vcg::math::IsNAN(s) && "normal should not have zero len!");
|
||||
uLinear *= s;
|
||||
uConstant *= s;
|
||||
}
|
||||
else
|
||||
{
|
||||
mStatus = ASS_UNDETERMINED;
|
||||
// normalize the gradient
|
||||
LScalar f = 1./sqrt(vcg::SquaredNorm(uLinear) - Scalar(4)*uConstant*uQuad);
|
||||
uConstant *= f;
|
||||
uLinear *= f;
|
||||
uQuad *= f;
|
||||
}
|
||||
// finalize
|
||||
if (fabs(uQuad)>1e-7)
|
||||
{
|
||||
mStatus = ASS_SPHERE;
|
||||
LScalar b = 1./uQuad;
|
||||
mCenter = uLinear*(-0.5*b);
|
||||
mRadius = sqrt( vcg::SquaredNorm(mCenter) - b*uConstant );
|
||||
}
|
||||
else if (uQuad==0.)
|
||||
{
|
||||
mStatus = ASS_PLANE;
|
||||
LScalar s = LScalar(1)/vcg::Norm(uLinear);
|
||||
assert(!vcg::math::IsNAN(s) && "normal should not have zero len!");
|
||||
uLinear *= s;
|
||||
uConstant *= s;
|
||||
}
|
||||
else
|
||||
{
|
||||
mStatus = ASS_UNDETERMINED;
|
||||
// normalize the gradient
|
||||
LScalar f = 1./sqrt(vcg::SquaredNorm(uLinear) - Scalar(4)*uConstant*uQuad);
|
||||
uConstant *= f;
|
||||
uLinear *= f;
|
||||
uQuad *= f;
|
||||
}
|
||||
|
||||
// cache some values to be used by the mls gradient
|
||||
mCachedSumP = sumP;
|
||||
mCachedSumN = sumN;
|
||||
mCachedSumW = sumW;
|
||||
mCachedSumDotPP = sumDotPP;
|
||||
mCachedSumDotPN = sumDotPN;
|
||||
// cache some values to be used by the mls gradient
|
||||
mCachedSumP = sumP;
|
||||
mCachedSumN = sumN;
|
||||
mCachedSumW = sumW;
|
||||
mCachedSumDotPP = sumDotPP;
|
||||
mCachedSumDotPN = sumDotPN;
|
||||
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
return true;
|
||||
}
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
return true;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool APSS<_MeshType>::mlsGradient(const VectorType& x, VectorType& grad) const
|
||||
{
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
template<typename _MeshType>
|
||||
bool APSS<_MeshType>::mlsGradient(const VectorType& x, VectorType& grad) const
|
||||
{
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
|
||||
const LVector& sumP = mCachedSumP;
|
||||
const LVector& sumN = mCachedSumN;
|
||||
const LScalar& sumDotPN = mCachedSumDotPN;
|
||||
const LScalar& sumDotPP = mCachedSumDotPP;
|
||||
const LScalar& sumW = mCachedSumW;
|
||||
const LScalar invSumW = 1.f/sumW;
|
||||
const LVector& sumP = mCachedSumP;
|
||||
const LVector& sumN = mCachedSumN;
|
||||
const LScalar& sumDotPN = mCachedSumDotPN;
|
||||
const LScalar& sumDotPP = mCachedSumDotPP;
|
||||
const LScalar& sumW = mCachedSumW;
|
||||
const LScalar invSumW = 1.f/sumW;
|
||||
|
||||
const LScalar nume = sumDotPN - invSumW * vcg::Dot(sumP, sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * vcg::Dot(sumP, sumP);
|
||||
const LScalar nume = sumDotPN - invSumW * (sumP* sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * (sumP* sumP);
|
||||
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
LVector dSumP; dSumP.SetZero();
|
||||
LVector dSumN; dSumN.SetZero();
|
||||
LScalar dSumDotPN = 0.;
|
||||
LScalar dSumDotPP = 0.;
|
||||
LScalar dSumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar dw = mCachedWeightGradients.at(i)[k];
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
LVector dSumP; dSumP.SetZero();
|
||||
LVector dSumN; dSumN.SetZero();
|
||||
LScalar dSumDotPN = 0.;
|
||||
LScalar dSumDotPP = 0.;
|
||||
LScalar dSumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar dw = mCachedWeightGradients.at(i)[k];
|
||||
|
||||
dSumW += dw;
|
||||
dSumP += p*dw;
|
||||
dSumN += n*dw;
|
||||
dSumDotPN += dw * vcg::Dot(n,p);
|
||||
dSumDotPP += dw * vcg::SquaredNorm(p);
|
||||
}
|
||||
dSumW += dw;
|
||||
dSumP += p*dw;
|
||||
dSumN += n*dw;
|
||||
dSumDotPN += dw * (n*p);
|
||||
dSumDotPP += dw * vcg::SquaredNorm(p);
|
||||
}
|
||||
|
||||
mCachedGradSumP[k] = dSumP;
|
||||
mCachedGradSumN[k] = dSumN;
|
||||
mCachedGradSumDotPN[k] = dSumDotPN;
|
||||
mCachedGradSumDotPP[k] = dSumDotPP;
|
||||
mCachedGradSumW[k] = dSumW;
|
||||
mCachedGradSumP[k] = dSumP;
|
||||
mCachedGradSumN[k] = dSumN;
|
||||
mCachedGradSumDotPN[k] = dSumDotPN;
|
||||
mCachedGradSumDotPP[k] = dSumDotPP;
|
||||
mCachedGradSumW[k] = dSumW;
|
||||
|
||||
LScalar dVecU0;
|
||||
LVector dVecU13;
|
||||
LScalar dVecU4;
|
||||
LScalar dVecU0;
|
||||
LVector dVecU13;
|
||||
LScalar dVecU4;
|
||||
|
||||
LScalar dNume = dSumDotPN - invSumW*invSumW*( sumW*(vcg::Dot(dSumP,sumN) + vcg::Dot(sumP,dSumN))
|
||||
- dSumW*vcg::Dot(sumP,sumN));
|
||||
LScalar dDeno = dSumDotPP - invSumW*invSumW*( 2.*sumW*vcg::Dot(dSumP,sumP)
|
||||
- dSumW*vcg::Dot(sumP,sumP));
|
||||
LScalar dNume = dSumDotPN - invSumW*invSumW*( sumW*((dSumP*sumN) + (sumP*dSumN))
|
||||
- dSumW*(sumP*sumN));
|
||||
LScalar dDeno = dSumDotPP - invSumW*invSumW*( 2.*sumW*(dSumP*sumP)
|
||||
- dSumW*(sumP*sumP));
|
||||
|
||||
dVecU4 = mSphericalParameter * 0.5 * (deno * dNume - dDeno * nume)/(deno*deno);
|
||||
dVecU13 = ((dSumN - (dSumP*uQuad + sumP*dVecU4)*2.0) - uLinear * dSumW) * invSumW;
|
||||
dVecU0 = -invSumW*( vcg::Dot(dVecU13,sumP) + dVecU4*sumDotPP + vcg::Dot(uLinear,dSumP) + uQuad*dSumDotPP + dSumW*uConstant);
|
||||
dVecU4 = mSphericalParameter * 0.5 * (deno * dNume - dDeno * nume)/(deno*deno);
|
||||
dVecU13 = ((dSumN - (dSumP*uQuad + sumP*dVecU4)*2.0) - uLinear * dSumW) * invSumW;
|
||||
dVecU0 = -invSumW*( (dVecU13*sumP) + dVecU4*sumDotPP + (uLinear*dSumP) + uQuad*dSumDotPP + dSumW*uConstant);
|
||||
|
||||
grad[k] = dVecU0 + vcg::Dot(dVecU13,vcg::Point3<LScalar>::Construct(x)) + dVecU4*vcg::SquaredNorm(x) + uLinear[k] + 2.*x[k]*uQuad;
|
||||
grad[k] = dVecU0 + (dVecU13*vcg::Point3<LScalar>::Construct(x)) + dVecU4*vcg::SquaredNorm(x) + uLinear[k] + 2.*x[k]*uQuad;
|
||||
|
||||
mCachedGradDeno[k] = dDeno;
|
||||
mCachedGradNume[k] = dNume;
|
||||
mCachedGradUConstant[k] = dVecU0;
|
||||
mCachedGradULinear[k] = dVecU13;
|
||||
mCachedGradUQuad[k] = dVecU4;
|
||||
}
|
||||
mCachedGradDeno[k] = dDeno;
|
||||
mCachedGradNume[k] = dNume;
|
||||
mCachedGradUConstant[k] = dVecU0;
|
||||
mCachedGradULinear[k] = dVecU13;
|
||||
mCachedGradUQuad[k] = dVecU4;
|
||||
}
|
||||
|
||||
return true;
|
||||
return true;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool APSS<_MeshType>::mlsHessian(const VectorType& x, MatrixType& hessian) const
|
||||
{
|
||||
this->requestSecondDerivatives();
|
||||
this->requestSecondDerivatives();
|
||||
|
||||
// TODO call mlsGradient first
|
||||
VectorType grad;
|
||||
mlsGradient(x,grad);
|
||||
// TODO call mlsGradient first
|
||||
VectorType grad;
|
||||
mlsGradient(x,grad);
|
||||
|
||||
uint nofSamples = mNeighborhood.size();
|
||||
uint nofSamples = mNeighborhood.size();
|
||||
|
||||
const LVector& sumP = mCachedSumP;
|
||||
const LVector& sumN = mCachedSumN;
|
||||
const LScalar& sumDotPN = mCachedSumDotPN;
|
||||
const LScalar& sumDotPP = mCachedSumDotPP;
|
||||
const LScalar& sumW = mCachedSumW;
|
||||
const LScalar invSumW = 1.f/sumW;
|
||||
const LVector& sumP = mCachedSumP;
|
||||
const LVector& sumN = mCachedSumN;
|
||||
const LScalar& sumDotPN = mCachedSumDotPN;
|
||||
const LScalar& sumDotPP = mCachedSumDotPP;
|
||||
const LScalar& sumW = mCachedSumW;
|
||||
const LScalar invSumW = 1.f/sumW;
|
||||
|
||||
const LScalar nume = sumDotPN - invSumW * vcg::Dot(sumP, sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * vcg::Dot(sumP, sumP);
|
||||
const LScalar nume = sumDotPN - invSumW * (sumP* sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * (sumP* sumP);
|
||||
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
const LVector& dSumP = mCachedGradSumP[k];
|
||||
const LVector& dSumN = mCachedGradSumN[k];
|
||||
const LScalar& dSumDotPN = mCachedGradSumDotPN[k];
|
||||
const LScalar& dSumDotPP = mCachedGradSumDotPP[k];
|
||||
const LScalar& dSumW = mCachedGradSumW[k];
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
const LVector& dSumP = mCachedGradSumP[k];
|
||||
const LVector& dSumN = mCachedGradSumN[k];
|
||||
const LScalar& dSumDotPN = mCachedGradSumDotPN[k];
|
||||
const LScalar& dSumDotPP = mCachedGradSumDotPP[k];
|
||||
const LScalar& dSumW = mCachedGradSumW[k];
|
||||
|
||||
LScalar dVecU0 = mCachedGradUConstant[k];
|
||||
LVector dVecU13 = mCachedGradULinear[k];
|
||||
LScalar dVecU4 = mCachedGradUQuad[k];
|
||||
LScalar dVecU0 = mCachedGradUConstant[k];
|
||||
LVector dVecU13 = mCachedGradULinear[k];
|
||||
LScalar dVecU4 = mCachedGradUQuad[k];
|
||||
|
||||
LScalar dNume = mCachedGradNume[k];
|
||||
LScalar dDeno = mCachedGradDeno[k];
|
||||
LScalar dNume = mCachedGradNume[k];
|
||||
LScalar dDeno = mCachedGradDeno[k];
|
||||
|
||||
// second order derivatives
|
||||
for (uint j=0 ; j<3 ; ++j)
|
||||
{
|
||||
LVector d2SumP; d2SumP.SetZero();
|
||||
LVector d2SumN; d2SumN.SetZero();
|
||||
LScalar d2SumDotPN = 0.;
|
||||
LScalar d2SumDotPP = 0.;
|
||||
LScalar d2SumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar dw = mCachedWeightGradients.at(i)[j];
|
||||
LScalar d2w = ((x[k]-p[k]))*((x[j]-p[j])) * mCachedWeightSecondDerivatives.at(i);
|
||||
// second order derivatives
|
||||
for (uint j=0 ; j<3 ; ++j)
|
||||
{
|
||||
LVector d2SumP; d2SumP.SetZero();
|
||||
LVector d2SumN; d2SumN.SetZero();
|
||||
LScalar d2SumDotPN = 0.;
|
||||
LScalar d2SumDotPP = 0.;
|
||||
LScalar d2SumW = 0.;
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
LVector p = vcg::Point3<LScalar>::Construct(mPoints[id].cP());
|
||||
LVector n = vcg::Point3<LScalar>::Construct(mPoints[id].cN());
|
||||
LScalar dw = mCachedWeightGradients.at(i)[j];
|
||||
LScalar d2w = ((x[k]-p[k]))*((x[j]-p[j])) * mCachedWeightSecondDerivatives.at(i);
|
||||
|
||||
if (j==k)
|
||||
d2w += mCachedWeightDerivatives.at(i);
|
||||
if (j==k)
|
||||
d2w += mCachedWeightDerivatives.at(i);
|
||||
|
||||
d2SumW += d2w;
|
||||
d2SumP += p*d2w;
|
||||
d2SumN += n*d2w;
|
||||
d2SumDotPN += d2w * vcg::Dot(n,p);
|
||||
d2SumDotPP += d2w * vcg::Dot(p,p);
|
||||
}
|
||||
d2SumW += d2w;
|
||||
d2SumP += p*d2w;
|
||||
d2SumN += n*d2w;
|
||||
d2SumDotPN += d2w * (n*p);
|
||||
d2SumDotPP += d2w * (p*p);
|
||||
}
|
||||
|
||||
LScalar d2u0;
|
||||
LVector d2u13;
|
||||
LScalar d2u4;
|
||||
LScalar d2u0;
|
||||
LVector d2u13;
|
||||
LScalar d2u4;
|
||||
|
||||
LScalar d2Nume = d2SumDotPN - invSumW*invSumW*invSumW*invSumW*(
|
||||
- 2.*sumW*mCachedGradSumW[j]*( sumW*(vcg::Dot(dSumP,sumN)+vcg::Dot(sumP,dSumN)) - dSumW* vcg::Dot(sumP,sumN))
|
||||
+ sumW*sumW*( mCachedGradSumW[j]*(vcg::Dot(dSumP,sumN)+vcg::Dot(sumP,dSumN))
|
||||
+ sumW*(vcg::Dot(d2SumP,sumN) + vcg::Dot(sumP,d2SumN)
|
||||
+ vcg::Dot(mCachedGradSumP[j],dSumN) + vcg::Dot(dSumP,mCachedGradSumN[j]))
|
||||
- d2SumW*vcg::Dot(sumP,sumN)
|
||||
- dSumW*(vcg::Dot(mCachedGradSumP[j],sumN)+vcg::Dot(sumP,mCachedGradSumN[j])) ));
|
||||
LScalar d2Nume = d2SumDotPN - invSumW*invSumW*invSumW*invSumW*(
|
||||
- 2.*sumW*mCachedGradSumW[j]*( sumW*((dSumP*sumN)+(sumP*dSumN)) - dSumW* (sumP*sumN))
|
||||
+ sumW*sumW*( mCachedGradSumW[j]*((dSumP*sumN)+(sumP*dSumN))
|
||||
+ sumW*((d2SumP*sumN) + (sumP*d2SumN)
|
||||
+ (mCachedGradSumP[j]*dSumN) + (dSumP*mCachedGradSumN[j]))
|
||||
- d2SumW*(sumP*sumN)
|
||||
- dSumW*((mCachedGradSumP[j]*sumN)+(sumP*mCachedGradSumN[j])) ));
|
||||
|
||||
LScalar d2Deno = d2SumDotPP - invSumW*invSumW*invSumW*invSumW*(
|
||||
- 2.*sumW*mCachedGradSumW[j] * ( 2.*sumW*vcg::Dot(dSumP,sumP) - dSumW*vcg::Dot(sumP,sumP))
|
||||
+ sumW*sumW*( 2.*mCachedGradSumW[j]*(vcg::Dot(dSumP,sumP))
|
||||
+ 2.*sumW*(vcg::Dot(mCachedGradSumP[j],dSumP)+vcg::Dot(d2SumP,sumP))
|
||||
- d2SumW*vcg::Dot(sumP,sumP) - dSumW*(2.*vcg::Dot(mCachedGradSumP[j],sumP))) );
|
||||
LScalar d2Deno = d2SumDotPP - invSumW*invSumW*invSumW*invSumW*(
|
||||
- 2.*sumW*mCachedGradSumW[j] * ( 2.*sumW*(dSumP*sumP) - dSumW*(sumP*sumP))
|
||||
+ sumW*sumW*( 2.*mCachedGradSumW[j]*((dSumP*sumP))
|
||||
+ 2.*sumW*((mCachedGradSumP[j]*dSumP)+(d2SumP*sumP))
|
||||
- d2SumW*(sumP*sumP) - dSumW*(2.*(mCachedGradSumP[j]*sumP))) );
|
||||
|
||||
LScalar deno2 = deno*deno;
|
||||
d2u4 = mSphericalParameter * 0.5 * (deno2*(d2Nume*deno + mCachedGradDeno[j] * dNume
|
||||
- d2Deno*nume - dDeno * mCachedGradNume[j])
|
||||
- 2.*deno*mCachedGradDeno[j]*(deno * dNume - dDeno * nume))/(deno2*deno2);
|
||||
LScalar deno2 = deno*deno;
|
||||
d2u4 = mSphericalParameter * 0.5 * (deno2*(d2Nume*deno + mCachedGradDeno[j] * dNume
|
||||
- d2Deno*nume - dDeno * mCachedGradNume[j])
|
||||
- 2.*deno*mCachedGradDeno[j]*(deno * dNume - dDeno * nume))/(deno2*deno2);
|
||||
|
||||
d2u13 = ( -dVecU13 * mCachedGradSumW[j]
|
||||
+ (d2SumN - (dSumP*mCachedGradUQuad[j] + d2SumP*uQuad + sumP*d2u4 + mCachedGradSumP[j]*dVecU4)*2.0 )
|
||||
- uLinear*d2SumW - mCachedGradULinear[j]*dSumW ) * invSumW;
|
||||
d2u13 = ( -dVecU13 * mCachedGradSumW[j]
|
||||
+ (d2SumN - (dSumP*mCachedGradUQuad[j] + d2SumP*uQuad + sumP*d2u4 + mCachedGradSumP[j]*dVecU4)*2.0 )
|
||||
- uLinear*d2SumW - mCachedGradULinear[j]*dSumW ) * invSumW;
|
||||
|
||||
d2u0 = ( -dVecU0 * mCachedGradSumW[j]
|
||||
- ( vcg::Dot(dVecU13,mCachedGradSumP[j]) + vcg::Dot(d2u13,sumP)
|
||||
+ d2u4*sumDotPP + dVecU4*mCachedGradSumDotPP[j]
|
||||
+ vcg::Dot(uLinear,d2SumP) + vcg::Dot(mCachedGradULinear[j],dSumP)
|
||||
+ dSumDotPP*mCachedGradUQuad[j] + d2SumDotPP*uQuad
|
||||
+ d2SumW*uConstant + dSumW*mCachedGradUConstant[j]) ) * invSumW;
|
||||
d2u0 = ( -dVecU0 * mCachedGradSumW[j]
|
||||
- ( (dVecU13*mCachedGradSumP[j]) + (d2u13*sumP)
|
||||
+ d2u4*sumDotPP + dVecU4*mCachedGradSumDotPP[j]
|
||||
+ (uLinear*d2SumP) + (mCachedGradULinear[j]*dSumP)
|
||||
+ dSumDotPP*mCachedGradUQuad[j] + d2SumDotPP*uQuad
|
||||
+ d2SumW*uConstant + dSumW*mCachedGradUConstant[j]) ) * invSumW;
|
||||
|
||||
hessian[j][k] =
|
||||
dVecU13[j] + 2.*dVecU4*x[j]
|
||||
+ d2u0 + vcg::Dot(d2u13,vcg::Point3<LScalar>::Construct(x)) + d2u4*vcg::Dot(x,x)
|
||||
+ mCachedGradULinear[j][k] + (j==k ? 2.*uQuad : 0.) + 2.*x[k]*mCachedGradUQuad[j];
|
||||
hessian[j][k] =
|
||||
dVecU13[j] + 2.*dVecU4*x[j]
|
||||
+ d2u0 + (d2u13*vcg::Point3<LScalar>::Construct(x)) + d2u4*(x*x)
|
||||
+ mCachedGradULinear[j][k] + (j==k ? 2.*uQuad : 0.) + 2.*x[k]*mCachedGradUQuad[j];
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -26,138 +26,139 @@
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _Scalar>
|
||||
BallTree<_Scalar>::BallTree(const ConstDataWrapper<VectorType>& points, const ConstDataWrapper<Scalar>& radii)
|
||||
: mPoints(points), mRadii(radii), mRadiusScale(1.), mTreeIsUptodate(false)
|
||||
BallTree<_Scalar>::BallTree(const vcg::ConstDataWrapper<VectorType>& points, const vcg::ConstDataWrapper<Scalar>& radii)
|
||||
: mPoints(points), mRadii(radii), mRadiusScale(1.), mTreeIsUptodate(false)
|
||||
{
|
||||
mRootNode = 0;
|
||||
mMaxTreeDepth = 12;
|
||||
mTargetCellSize = 24;
|
||||
mRootNode = 0;
|
||||
mMaxTreeDepth = 12;
|
||||
mTargetCellSize = 24;
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::computeNeighbors(const VectorType& x, Neighborhood<Scalar>* pNei) const
|
||||
{
|
||||
if (!mTreeIsUptodate)
|
||||
const_cast<BallTree*>(this)->rebuild();
|
||||
if (!mTreeIsUptodate)
|
||||
const_cast<BallTree*>(this)->rebuild();
|
||||
|
||||
pNei->clear();
|
||||
mQueryPosition = x;
|
||||
queryNode(*mRootNode, pNei);
|
||||
pNei->clear();
|
||||
mQueryPosition = x;
|
||||
queryNode(*mRootNode, pNei);
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::queryNode(Node& node, Neighborhood<Scalar>* pNei) const
|
||||
{
|
||||
if (node.leaf)
|
||||
{
|
||||
for (unsigned int i=0 ; i<node.size ; ++i)
|
||||
{
|
||||
int id = node.indices[i];
|
||||
Scalar d2 = vcg::SquaredNorm(mQueryPosition - mPoints[id]);
|
||||
Scalar r = mRadiusScale * mRadii[id];
|
||||
if (d2<r*r)
|
||||
pNei->insert(id, d2);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (mQueryPosition[node.dim] - node.splitValue < 0)
|
||||
queryNode(*node.children[0], pNei);
|
||||
else
|
||||
queryNode(*node.children[1], pNei);
|
||||
}
|
||||
if (node.leaf)
|
||||
{
|
||||
for (unsigned int i=0 ; i<node.size ; ++i)
|
||||
{
|
||||
int id = node.indices[i];
|
||||
Scalar d2 = vcg::SquaredNorm(mQueryPosition - mPoints[id]);
|
||||
Scalar r = mRadiusScale * mRadii[id];
|
||||
if (d2<r*r)
|
||||
pNei->insert(id, d2);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (mQueryPosition[node.dim] - node.splitValue < 0)
|
||||
queryNode(*node.children[0], pNei);
|
||||
else
|
||||
queryNode(*node.children[1], pNei);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline vcg::Point3<Scalar> CwiseAdd(vcg::Point3<Scalar> const & p1, Scalar s)
|
||||
{
|
||||
return vcg::Point3<Scalar>(p1.X() + s, p1.Y() + s, p1.Z() + s);
|
||||
return vcg::Point3<Scalar>(p1.X() + s, p1.Y() + s, p1.Z() + s);
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::rebuild(void)
|
||||
{
|
||||
delete mRootNode;
|
||||
delete mRootNode;
|
||||
|
||||
mRootNode = new Node();
|
||||
IndexArray indices(mPoints.size());
|
||||
AxisAlignedBoxType aabb;
|
||||
aabb.Set(mPoints[0]);
|
||||
for (unsigned int i=0 ; i<mPoints.size() ; ++i)
|
||||
{
|
||||
indices[i] = i;
|
||||
aabb.Add(mPoints[i],mRadii[i]*mRadiusScale);
|
||||
mRootNode = new Node();
|
||||
IndexArray indices(mPoints.size());
|
||||
AxisAlignedBoxType aabb;
|
||||
aabb.Set(mPoints[0]);
|
||||
for (unsigned int i=0 ; i<mPoints.size() ; ++i)
|
||||
{
|
||||
indices[i] = i;
|
||||
aabb.Add(mPoints[i],mRadii[i]*mRadiusScale);
|
||||
// aabb.min = Min(aabb.min, CwiseAdd(mPoints[i], -mRadii[i]*mRadiusScale));
|
||||
// aabb.max = Max(aabb.max, CwiseAdd(mPoints[i], mRadii[i]*mRadiusScale));
|
||||
}
|
||||
buildNode(*mRootNode, indices, aabb, 0);
|
||||
}
|
||||
buildNode(*mRootNode, indices, aabb, 0);
|
||||
|
||||
mTreeIsUptodate = true;
|
||||
mTreeIsUptodate = true;
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::split(const IndexArray& indices, const AxisAlignedBoxType& aabbLeft, const AxisAlignedBoxType& aabbRight, IndexArray& iLeft, IndexArray& iRight)
|
||||
{
|
||||
for (std::vector<int>::const_iterator it=indices.begin(), end=indices.end() ; it!=end ; ++it)
|
||||
{
|
||||
unsigned int i = *it;
|
||||
if (vcg::PointFilledBoxDistance(mPoints[i], aabbLeft) < mRadii[i]*mRadiusScale)
|
||||
iLeft.push_back(i);
|
||||
for (std::vector<int>::const_iterator it=indices.begin(), end=indices.end() ; it!=end ; ++it)
|
||||
{
|
||||
unsigned int i = *it;
|
||||
if (vcg::PointFilledBoxDistance(mPoints[i], aabbLeft) < mRadii[i]*mRadiusScale)
|
||||
iLeft.push_back(i);
|
||||
|
||||
if (vcg::PointFilledBoxDistance(mPoints[i], aabbRight) < mRadii[i]*mRadiusScale)
|
||||
iRight.push_back(i);
|
||||
}
|
||||
if (vcg::PointFilledBoxDistance(mPoints[i], aabbRight) < mRadii[i]*mRadiusScale)
|
||||
iRight.push_back(i);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::buildNode(Node& node, std::vector<int>& indices, AxisAlignedBoxType aabb, int level)
|
||||
{
|
||||
Scalar avgradius = 0.;
|
||||
for (std::vector<int>::const_iterator it=indices.begin(), end=indices.end() ; it!=end ; ++it)
|
||||
avgradius += mRadii[*it];
|
||||
avgradius = mRadiusScale * avgradius / Scalar(indices.size());
|
||||
VectorType diag = aabb.max - aabb.min;
|
||||
if (int(indices.size())<mTargetCellSize
|
||||
|| avgradius*0.9 > std::max(std::max(diag.X(), diag.Y()), diag.Z())
|
||||
|| int(level)>=mMaxTreeDepth)
|
||||
{
|
||||
node.leaf = true;
|
||||
node.size = indices.size();
|
||||
node.indices = new unsigned int[node.size];
|
||||
for (unsigned int i=0 ; i<node.size ; ++i)
|
||||
node.indices[i] = indices[i];
|
||||
return;
|
||||
}
|
||||
unsigned int dim = vcg::MaxCoeffId(diag);
|
||||
node.dim = dim;
|
||||
node.splitValue = Scalar(0.5*(aabb.max[dim] + aabb.min[dim]));
|
||||
node.leaf = 0;
|
||||
Scalar avgradius = 0.;
|
||||
for (std::vector<int>::const_iterator it=indices.begin(), end=indices.end() ; it!=end ; ++it)
|
||||
avgradius += mRadii[*it];
|
||||
avgradius = mRadiusScale * avgradius / Scalar(indices.size());
|
||||
VectorType diag = aabb.max - aabb.min;
|
||||
if (int(indices.size())<mTargetCellSize
|
||||
|| avgradius*0.9 > std::max(std::max(diag.X(), diag.Y()), diag.Z())
|
||||
|| int(level)>=mMaxTreeDepth)
|
||||
{
|
||||
node.leaf = true;
|
||||
node.size = indices.size();
|
||||
node.indices = new unsigned int[node.size];
|
||||
for (unsigned int i=0 ; i<node.size ; ++i)
|
||||
node.indices[i] = indices[i];
|
||||
return;
|
||||
}
|
||||
|
||||
AxisAlignedBoxType aabbLeft=aabb, aabbRight=aabb;
|
||||
aabbLeft.max[dim] = node.splitValue;
|
||||
aabbRight.min[dim] = node.splitValue;
|
||||
unsigned int dim = diag.MaxCoeffId();
|
||||
node.dim = dim;
|
||||
node.splitValue = Scalar(0.5*(aabb.max[dim] + aabb.min[dim]));
|
||||
node.leaf = 0;
|
||||
|
||||
std::vector<int> iLeft, iRight;
|
||||
split(indices, aabbLeft, aabbRight, iLeft,iRight);
|
||||
AxisAlignedBoxType aabbLeft=aabb, aabbRight=aabb;
|
||||
aabbLeft.max[dim] = node.splitValue;
|
||||
aabbRight.min[dim] = node.splitValue;
|
||||
|
||||
// we don't need the index list anymore
|
||||
indices.clear();
|
||||
std::vector<int> iLeft, iRight;
|
||||
split(indices, aabbLeft, aabbRight, iLeft,iRight);
|
||||
|
||||
{
|
||||
// left child
|
||||
//mNodes.resize(mNodes.size()+1);
|
||||
Node* pChild = new Node();
|
||||
node.children[0] = pChild;
|
||||
buildNode(*pChild, iLeft, aabbLeft, level+1);
|
||||
}
|
||||
// we don't need the index list anymore
|
||||
indices.clear();
|
||||
|
||||
{
|
||||
// right child
|
||||
//mNodes.resize(mNodes.size()+1);
|
||||
Node* pChild = new Node();
|
||||
node.children[1] = pChild;
|
||||
buildNode(*pChild, iRight, aabbRight, level+1);
|
||||
}
|
||||
{
|
||||
// left child
|
||||
//mNodes.resize(mNodes.size()+1);
|
||||
Node* pChild = new Node();
|
||||
node.children[0] = pChild;
|
||||
buildNode(*pChild, iLeft, aabbLeft, level+1);
|
||||
}
|
||||
|
||||
{
|
||||
// right child
|
||||
//mNodes.resize(mNodes.size()+1);
|
||||
Node* pChild = new Node();
|
||||
node.children[1] = pChild;
|
||||
buildNode(*pChild, iRight, aabbRight, level+1);
|
||||
}
|
||||
}
|
||||
|
||||
template class BallTree<float>;
|
||||
|
||||
@ -33,85 +33,85 @@ namespace GaelMls {
|
||||
template<typename _Scalar>
|
||||
class Neighborhood
|
||||
{
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
|
||||
int index(int i) const { return mIndices.at(i); }
|
||||
Scalar squaredDistance(int i) const { return mSqDists.at(i); }
|
||||
int index(int i) const { return mIndices.at(i); }
|
||||
Scalar squaredDistance(int i) const { return mSqDists.at(i); }
|
||||
|
||||
void clear() { mIndices.clear(); mSqDists.clear(); }
|
||||
void resize(int size) { mIndices.resize(size); mSqDists.resize(size); }
|
||||
void reserve(int size) { mIndices.reserve(size); mSqDists.reserve(size); }
|
||||
int size() { return mIndices.size(); }
|
||||
void clear() { mIndices.clear(); mSqDists.clear(); }
|
||||
void resize(int size) { mIndices.resize(size); mSqDists.resize(size); }
|
||||
void reserve(int size) { mIndices.reserve(size); mSqDists.reserve(size); }
|
||||
int size() { return mIndices.size(); }
|
||||
|
||||
void insert(int id, Scalar d2) { mIndices.push_back(id); mSqDists.push_back(d2); }
|
||||
void insert(int id, Scalar d2) { mIndices.push_back(id); mSqDists.push_back(d2); }
|
||||
|
||||
protected:
|
||||
std::vector<int> mIndices;
|
||||
std::vector<Scalar> mSqDists;
|
||||
protected:
|
||||
std::vector<int> mIndices;
|
||||
std::vector<Scalar> mSqDists;
|
||||
};
|
||||
|
||||
template<typename _Scalar>
|
||||
class BallTree
|
||||
{
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
|
||||
BallTree(const ConstDataWrapper<VectorType>& points, const ConstDataWrapper<Scalar>& radii);
|
||||
BallTree(const vcg::ConstDataWrapper<VectorType>& points, const vcg::ConstDataWrapper<Scalar>& radii);
|
||||
|
||||
void computeNeighbors(const VectorType& x, Neighborhood<Scalar>* pNei) const;
|
||||
void computeNeighbors(const VectorType& x, Neighborhood<Scalar>* pNei) const;
|
||||
|
||||
void setRadiusScale(Scalar v) { mRadiusScale = v; mTreeIsUptodate = false; }
|
||||
void setRadiusScale(Scalar v) { mRadiusScale = v; mTreeIsUptodate = false; }
|
||||
|
||||
protected:
|
||||
protected:
|
||||
|
||||
struct Node
|
||||
{
|
||||
~Node()
|
||||
{
|
||||
if (!leaf)
|
||||
{
|
||||
delete children[0];
|
||||
delete children[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
delete[] indices;
|
||||
}
|
||||
}
|
||||
Scalar splitValue;
|
||||
unsigned char dim:2;
|
||||
unsigned char leaf:1;
|
||||
union {
|
||||
Node* children[2];
|
||||
struct {
|
||||
unsigned int* indices;
|
||||
unsigned int size;
|
||||
};
|
||||
};
|
||||
};
|
||||
struct Node
|
||||
{
|
||||
~Node()
|
||||
{
|
||||
if (!leaf)
|
||||
{
|
||||
delete children[0];
|
||||
delete children[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
delete[] indices;
|
||||
}
|
||||
}
|
||||
Scalar splitValue;
|
||||
unsigned char dim:2;
|
||||
unsigned char leaf:1;
|
||||
union {
|
||||
Node* children[2];
|
||||
struct {
|
||||
unsigned int* indices;
|
||||
unsigned int size;
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
typedef std::vector<int> IndexArray;
|
||||
typedef vcg::Box3<Scalar> AxisAlignedBoxType;
|
||||
typedef std::vector<int> IndexArray;
|
||||
typedef vcg::Box3<Scalar> AxisAlignedBoxType;
|
||||
|
||||
void rebuild();
|
||||
void split(const IndexArray& indices, const AxisAlignedBoxType& aabbLeft, const AxisAlignedBoxType& aabbRight,
|
||||
IndexArray& iLeft, IndexArray& iRight);
|
||||
void buildNode(Node& node, std::vector<int>& indices, AxisAlignedBoxType aabb, int level);
|
||||
void queryNode(Node& node, Neighborhood<Scalar>* pNei) const;
|
||||
void rebuild();
|
||||
void split(const IndexArray& indices, const AxisAlignedBoxType& aabbLeft, const AxisAlignedBoxType& aabbRight,
|
||||
IndexArray& iLeft, IndexArray& iRight);
|
||||
void buildNode(Node& node, std::vector<int>& indices, AxisAlignedBoxType aabb, int level);
|
||||
void queryNode(Node& node, Neighborhood<Scalar>* pNei) const;
|
||||
|
||||
protected:
|
||||
ConstDataWrapper<VectorType> mPoints;
|
||||
ConstDataWrapper<Scalar> mRadii;
|
||||
Scalar mRadiusScale;
|
||||
protected:
|
||||
vcg::ConstDataWrapper<VectorType> mPoints;
|
||||
vcg::ConstDataWrapper<Scalar> mRadii;
|
||||
Scalar mRadiusScale;
|
||||
|
||||
int mMaxTreeDepth;
|
||||
int mTargetCellSize;
|
||||
mutable bool mTreeIsUptodate;
|
||||
mutable VectorType mQueryPosition;
|
||||
int mMaxTreeDepth;
|
||||
int mTargetCellSize;
|
||||
mutable bool mTreeIsUptodate;
|
||||
mutable VectorType mQueryPosition;
|
||||
|
||||
Node* mRootNode;
|
||||
Node* mRootNode;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -38,163 +38,163 @@ namespace implicits
|
||||
template<typename Scalar>
|
||||
Scalar GaussCurvature(const Point3<Scalar>& gradient, const Matrix33<Scalar>& hessian)
|
||||
{
|
||||
Scalar l2 = gradient.SquaredNorm();
|
||||
Matrix33<Scalar> adjugate;
|
||||
adjugate[0][0] = hessian[1][1]*hessian[2][2] - hessian[1][2]*hessian[2][1];
|
||||
adjugate[1][0] = hessian[0][2]*hessian[2][1] - hessian[0][1]*hessian[2][2];
|
||||
adjugate[2][0] = hessian[0][1]*hessian[1][2] - hessian[0][2]*hessian[1][1];
|
||||
Scalar l2 = gradient.SquaredNorm();
|
||||
Matrix33<Scalar> adjugate;
|
||||
adjugate[0][0] = hessian[1][1]*hessian[2][2] - hessian[1][2]*hessian[2][1];
|
||||
adjugate[1][0] = hessian[0][2]*hessian[2][1] - hessian[0][1]*hessian[2][2];
|
||||
adjugate[2][0] = hessian[0][1]*hessian[1][2] - hessian[0][2]*hessian[1][1];
|
||||
|
||||
adjugate[0][1] = hessian[1][2]*hessian[2][0] - hessian[1][0]*hessian[2][2];
|
||||
adjugate[1][1] = hessian[0][0]*hessian[2][2] - hessian[0][2]*hessian[2][0];
|
||||
adjugate[2][1] = hessian[1][0]*hessian[0][2] - hessian[0][0]*hessian[1][2];
|
||||
adjugate[0][1] = hessian[1][2]*hessian[2][0] - hessian[1][0]*hessian[2][2];
|
||||
adjugate[1][1] = hessian[0][0]*hessian[2][2] - hessian[0][2]*hessian[2][0];
|
||||
adjugate[2][1] = hessian[1][0]*hessian[0][2] - hessian[0][0]*hessian[1][2];
|
||||
|
||||
adjugate[0][2] = hessian[1][0]*hessian[2][1] - hessian[1][1]*hessian[2][0];
|
||||
adjugate[1][2] = hessian[0][1]*hessian[2][0] - hessian[0][0]*hessian[2][1];
|
||||
adjugate[2][2] = hessian[0][0]*hessian[1][1] - hessian[0][1]*hessian[1][0];
|
||||
return (vcg::Dot(gradient, adjugate * gradient)) / (l2*l2);
|
||||
adjugate[0][2] = hessian[1][0]*hessian[2][1] - hessian[1][1]*hessian[2][0];
|
||||
adjugate[1][2] = hessian[0][1]*hessian[2][0] - hessian[0][0]*hessian[2][1];
|
||||
adjugate[2][2] = hessian[0][0]*hessian[1][1] - hessian[0][1]*hessian[1][0];
|
||||
return ((gradient* ( adjugate * gradient))) / (l2*l2);
|
||||
}
|
||||
|
||||
/** \returns the mean curvature directly from the gradient and hessian */
|
||||
template<typename Scalar>
|
||||
Scalar MeanCurvature(const Point3<Scalar>& gradient, const Matrix33<Scalar>& hessian)
|
||||
{
|
||||
Scalar l = gradient.Norm();
|
||||
return (l*l*hessian.Trace() - vcg::Dot(gradient, hessian * gradient)) / (2.*l*l*l);
|
||||
Scalar l = gradient.Norm();
|
||||
return (l*l*hessian.Trace() - (gradient* (hessian * gradient))) / (2.*l*l*l);
|
||||
}
|
||||
|
||||
/** This class computes the Weingarten map of a scalar field and provides
|
||||
* methods to extract curvatures from it.
|
||||
*
|
||||
* The Weingarten map is equal to the gradient of the normal vector:
|
||||
* \f$ W = \nabla \mathbf(n)
|
||||
* = \nabla \frac{\mathbf{g}}{\Vert \mathbf{g} \Vert}
|
||||
* = \frac{(I - n n^T) H}{\Vert \mathbf{g} \Vert} \f$
|
||||
* This matrix can also be seen as the projection of the hessian
|
||||
* matrix onto the tangent plane of normal n.
|
||||
*/
|
||||
* methods to extract curvatures from it.
|
||||
*
|
||||
* The Weingarten map is equal to the gradient of the normal vector:
|
||||
* \f$ W = \nabla \mathbf(n)
|
||||
* = \nabla \frac{\mathbf{g}}{\Vert \mathbf{g} \Vert}
|
||||
* = \frac{(I - n n^T) H}{\Vert \mathbf{g} \Vert} \f$
|
||||
* This matrix can also be seen as the projection of the hessian
|
||||
* matrix onto the tangent plane of normal n.
|
||||
*/
|
||||
template<typename Scalar> class WeingartenMap
|
||||
{
|
||||
public:
|
||||
typedef Point3<Scalar> VectorType;
|
||||
typedef Matrix33<Scalar> MatrixType;
|
||||
typedef Point3<Scalar> VectorType;
|
||||
typedef Matrix33<Scalar> MatrixType;
|
||||
|
||||
/** Default constructor computing the Weingarten map from the
|
||||
* first and second derivatives, namely the gradient vector
|
||||
* and hessian matrix of the scalar field.
|
||||
*/
|
||||
WeingartenMap(const VectorType& grad, const MatrixType& hess)
|
||||
{
|
||||
Scalar invL = 1.0/grad.Norm();
|
||||
assert(grad.Norm()>1e-8);
|
||||
m_normal = grad * invL;
|
||||
/** Default constructor computing the Weingarten map from the
|
||||
* first and second derivatives, namely the gradient vector
|
||||
* and hessian matrix of the scalar field.
|
||||
*/
|
||||
WeingartenMap(const VectorType& grad, const MatrixType& hess)
|
||||
{
|
||||
Scalar invL = 1.0/grad.Norm();
|
||||
assert(grad.Norm()>1e-8);
|
||||
m_normal = grad * invL;
|
||||
assert(!math::IsNAN(invL) && "gradient should not be zero!");
|
||||
|
||||
Matrix33<Scalar> I; I.SetIdentity();
|
||||
m_nnT.ExternalProduct(m_normal,m_normal);
|
||||
|
||||
m_w = (I-m_nnT) * hess * invL;
|
||||
Matrix33<Scalar> I; I.SetIdentity();
|
||||
m_nnT.ExternalProduct(m_normal,m_normal);
|
||||
|
||||
m_kgIsDirty = true;
|
||||
m_kmIsDirty = true;
|
||||
m_kpAreDirty = true;
|
||||
m_kdirAreDirty = true;
|
||||
}
|
||||
m_w = (I-m_nnT) * hess * invL;
|
||||
|
||||
/** \returns the Weingarten map matrix */
|
||||
const MatrixType& W() const { return m_w; }
|
||||
m_kgIsDirty = true;
|
||||
m_kmIsDirty = true;
|
||||
m_kpAreDirty = true;
|
||||
m_kdirAreDirty = true;
|
||||
}
|
||||
|
||||
/** \returns the Gauss curvature = k1 * k2 */
|
||||
Scalar GaussCurvature() const
|
||||
{
|
||||
if (m_kgIsDirty)
|
||||
{
|
||||
// we add nn^T to W such that the third eigenvalue becomes 1
|
||||
// then det(W) = k1 * k2 * 1 = Gauss curvature !
|
||||
m_kg = (m_w + m_nnT).Determinant();
|
||||
m_kgIsDirty = false;
|
||||
}
|
||||
return m_kg;
|
||||
}
|
||||
/** \returns the Weingarten map matrix */
|
||||
const MatrixType& W() const { return m_w; }
|
||||
|
||||
/** \returns the mean curvature = (k1 + k2)/2 */
|
||||
Scalar MeanCurvature() const
|
||||
{
|
||||
if (m_kmIsDirty)
|
||||
{
|
||||
// the third eigenvalue of W is 0, then tr(W) = k1 + k2 + 0 = 2 k mean !
|
||||
m_km = m_w.Trace();
|
||||
m_kmIsDirty = false;
|
||||
}
|
||||
return m_km;
|
||||
}
|
||||
/** \returns the Gauss curvature = k1 * k2 */
|
||||
Scalar GaussCurvature() const
|
||||
{
|
||||
if (m_kgIsDirty)
|
||||
{
|
||||
// we add nn^T to W such that the third eigenvalue becomes 1
|
||||
// then det(W) = k1 * k2 * 1 = Gauss curvature !
|
||||
m_kg = (m_w + m_nnT).Determinant();
|
||||
m_kgIsDirty = false;
|
||||
}
|
||||
return m_kg;
|
||||
}
|
||||
|
||||
/** \returns the first principal curvature */
|
||||
Scalar K1() const { updateKp(); return m_k1; }
|
||||
/** \returns the mean curvature = (k1 + k2)/2 */
|
||||
Scalar MeanCurvature() const
|
||||
{
|
||||
if (m_kmIsDirty)
|
||||
{
|
||||
// the third eigenvalue of W is 0, then tr(W) = k1 + k2 + 0 = 2 k mean !
|
||||
m_km = m_w.Trace();
|
||||
m_kmIsDirty = false;
|
||||
}
|
||||
return m_km;
|
||||
}
|
||||
|
||||
/** \returns the second principal curvature */
|
||||
Scalar K2() const { updateKp(); return m_k2; }
|
||||
/** \returns the first principal curvature */
|
||||
Scalar K1() const { updateKp(); return m_k1; }
|
||||
|
||||
/** \returns the direction of the first principal curvature */
|
||||
const VectorType& K1Dir() const { extractEigenvectors(); return m_k1dir; }
|
||||
/** \returns the second principal curvature */
|
||||
Scalar K2() const { updateKp(); return m_k2; }
|
||||
|
||||
/** \returns the direction of the second principal curvature */
|
||||
const VectorType& K2Dir() const { extractEigenvectors(); return m_k2dir; }
|
||||
/** \returns the direction of the first principal curvature */
|
||||
const VectorType& K1Dir() const { extractEigenvectors(); return m_k1dir; }
|
||||
|
||||
/** \returns the direction of the second principal curvature */
|
||||
const VectorType& K2Dir() const { extractEigenvectors(); return m_k2dir; }
|
||||
|
||||
protected:
|
||||
|
||||
// direct computation of principal curvatures if needed
|
||||
inline void updateKp() const
|
||||
{
|
||||
if (m_kpAreDirty)
|
||||
{
|
||||
Scalar delta = sqrt(MeanCurvature()*m_km - 4.0*GaussCurvature());
|
||||
m_k1 = 0.5*(m_km + delta);
|
||||
m_k2 = 0.5*(m_km - delta);
|
||||
if (fabs(m_k1)<fabs(m_k2))
|
||||
std::swap(m_k1,m_k2);
|
||||
m_kpAreDirty = false;
|
||||
}
|
||||
}
|
||||
// direct computation of principal curvatures if needed
|
||||
inline void updateKp() const
|
||||
{
|
||||
if (m_kpAreDirty)
|
||||
{
|
||||
Scalar delta = sqrt(MeanCurvature()*m_km - 4.0*GaussCurvature());
|
||||
m_k1 = 0.5*(m_km + delta);
|
||||
m_k2 = 0.5*(m_km - delta);
|
||||
if (fabs(m_k1)<fabs(m_k2))
|
||||
std::swap(m_k1,m_k2);
|
||||
m_kpAreDirty = false;
|
||||
}
|
||||
}
|
||||
|
||||
inline void extractEigenvectors() const
|
||||
{
|
||||
if (m_kdirAreDirty)
|
||||
{
|
||||
Eigen::Matrix<Scalar,3,3> copy;
|
||||
m_w.ToEigenMatrix(copy);
|
||||
inline void extractEigenvectors() const
|
||||
{
|
||||
if (m_kdirAreDirty)
|
||||
{
|
||||
Eigen::Matrix<Scalar,3,3> copy;
|
||||
m_w.ToEigenMatrix(copy);
|
||||
// MatrixType copy = m_w;
|
||||
// int mrot = 0;
|
||||
// VectorType evals;
|
||||
// MatrixType evecs;
|
||||
|
||||
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eig(copy);
|
||||
Eigen::Vector3f eval = eig.eigenvalues();
|
||||
Eigen::Matrix3f evec = eig.eigenvectors();
|
||||
eval = eval.cwiseAbs();
|
||||
int ind0,ind1,ind2;
|
||||
eval.minCoeff(&ind0);
|
||||
ind1=(ind0+1)%3;
|
||||
ind2=(ind0+2)%3;
|
||||
if(eval[ind1]>eval[ind2]) std::swap(ind1,ind2);
|
||||
Eigen::SelfAdjointEigenSolver<Eigen::Matrix<Scalar,3,3> > eig(copy);
|
||||
Eigen::Matrix<Scalar,3,1> eval = eig.eigenvalues();
|
||||
Eigen::Matrix<Scalar,3,3> evec = eig.eigenvectors();
|
||||
eval = eval.cwiseAbs();
|
||||
int ind0,ind1,ind2;
|
||||
eval.minCoeff(&ind0);
|
||||
ind1=(ind0+1)%3;
|
||||
ind2=(ind0+2)%3;
|
||||
if(eval[ind1]>eval[ind2]) std::swap(ind1,ind2);
|
||||
|
||||
// Jacobi(copy, evals, evecs, mrot);
|
||||
// VectorType evalsAbs(fabs(evals[0]),fabs(evals[0]),fabs(evals[0]));
|
||||
// SortEigenvaluesAndEigenvectors(evals,evecs,true);
|
||||
m_k1 = eval[ind1];
|
||||
m_k2 = eval[ind2];
|
||||
m_k1dir.FromEigenVector(evec.col(ind1));
|
||||
m_k2dir.FromEigenVector(evec.col(ind2));
|
||||
m_kdirAreDirty = false;
|
||||
}
|
||||
}
|
||||
m_k1 = eval[ind1];
|
||||
m_k2 = eval[ind2];
|
||||
m_k1dir.FromEigenVector(evec.col(ind1));
|
||||
m_k2dir.FromEigenVector(evec.col(ind2));
|
||||
m_kdirAreDirty = false;
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
VectorType m_normal;
|
||||
MatrixType m_nnT;
|
||||
MatrixType m_w;
|
||||
mutable VectorType m_k1dir, m_k2dir;
|
||||
mutable Scalar m_kg, m_km, m_k1, m_k2;
|
||||
mutable bool m_kgIsDirty, m_kmIsDirty, m_kpAreDirty, m_kdirAreDirty;
|
||||
VectorType m_normal;
|
||||
MatrixType m_nnT;
|
||||
MatrixType m_w;
|
||||
mutable VectorType m_k1dir, m_k2dir;
|
||||
mutable Scalar m_kg, m_km, m_k1, m_k2;
|
||||
mutable bool m_kgIsDirty, m_kmIsDirty, m_kpAreDirty, m_kdirAreDirty;
|
||||
};
|
||||
|
||||
} // namespace implicits
|
||||
|
||||
@ -36,240 +36,241 @@ template <class MeshType, class SurfaceType>
|
||||
class MlsWalker
|
||||
{
|
||||
private:
|
||||
typedef int VertexIndex;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef vcg::Point3f VectorType;
|
||||
typedef typename MeshType::VertexPointer VertexPointer;
|
||||
//typedef GaelMls::MlsSurface<MeshType> SurfaceType;
|
||||
typedef long long unsigned int Key;
|
||||
typedef std::map<Key,VertexIndex> MapType;
|
||||
typedef int VertexIndex;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::CoordType VectorType;
|
||||
typedef typename MeshType::VertexPointer VertexPointer;
|
||||
//typedef GaelMls::MlsSurface<MeshType> SurfaceType;
|
||||
typedef long long unsigned int Key;
|
||||
typedef std::map<Key,VertexIndex> MapType;
|
||||
|
||||
struct GridElement
|
||||
{
|
||||
VectorType position;
|
||||
ScalarType value;
|
||||
};
|
||||
struct GridElement
|
||||
{
|
||||
VectorType position;
|
||||
ScalarType value;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
inline bool IsFinite(T value)
|
||||
{
|
||||
return (value>=-std::numeric_limits<T>::max()) && (value<=std::numeric_limits<T>::max());
|
||||
}
|
||||
template <typename T>
|
||||
inline bool IsFinite(T value)
|
||||
{
|
||||
return (value>=-std::numeric_limits<T>::max()) && (value<=std::numeric_limits<T>::max());
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
int resolution;
|
||||
int resolution;
|
||||
|
||||
MlsWalker()
|
||||
{
|
||||
resolution = 150;
|
||||
mMaxBlockSize = 64;
|
||||
mIsoValue = 0;
|
||||
}
|
||||
MlsWalker()
|
||||
{
|
||||
resolution = 150;
|
||||
mMaxBlockSize = 64;
|
||||
mIsoValue = 0;
|
||||
}
|
||||
|
||||
template<class EXTRACTOR_TYPE>
|
||||
void BuildMesh(MeshType &mesh, SurfaceType &surface, EXTRACTOR_TYPE &extractor, CallBackPos *cb = 0)
|
||||
{
|
||||
mpSurface = &surface;
|
||||
const ScalarType invalidValue = SurfaceType::InvalidValue();
|
||||
template<class EXTRACTOR_TYPE>
|
||||
void BuildMesh(MeshType &mesh, SurfaceType &surface, EXTRACTOR_TYPE &extractor, CallBackPos *cb = 0)
|
||||
{
|
||||
mpSurface = &surface;
|
||||
const ScalarType invalidValue = SurfaceType::InvalidValue();
|
||||
|
||||
// precomputed offsets to access the cell corners
|
||||
const int offsets[8] = {
|
||||
0,
|
||||
1,
|
||||
1+mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize,
|
||||
1+mMaxBlockSize,
|
||||
1+mMaxBlockSize+mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize+mMaxBlockSize*mMaxBlockSize
|
||||
};
|
||||
// precomputed offsets to access the cell corners
|
||||
const int offsets[8] = {
|
||||
0,
|
||||
1,
|
||||
1+mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize,
|
||||
1+mMaxBlockSize,
|
||||
1+mMaxBlockSize+mMaxBlockSize*mMaxBlockSize,
|
||||
mMaxBlockSize+mMaxBlockSize*mMaxBlockSize
|
||||
};
|
||||
|
||||
mAABB = mpSurface->boundingBox();
|
||||
VectorType diag = mAABB.max - mAABB.min;
|
||||
mAABB.min -= diag * 0.1f;
|
||||
mAABB.max += diag * 0.1f;
|
||||
diag = mAABB.max - mAABB.min;
|
||||
mAABB = mpSurface->boundingBox();
|
||||
|
||||
if ( (diag[0]<=0.)
|
||||
|| (diag[1]<=0.)
|
||||
|| (diag[2]<=0.)
|
||||
|| (resolution==0))
|
||||
{
|
||||
return;
|
||||
}
|
||||
VectorType diag = mAABB.max - mAABB.min;
|
||||
mAABB.min -= diag * 0.1f;
|
||||
mAABB.max += diag * 0.1f;
|
||||
diag = mAABB.max - mAABB.min;
|
||||
|
||||
mCache = new GridElement[(mMaxBlockSize)*(mMaxBlockSize)*(mMaxBlockSize)];
|
||||
ScalarType step = vcg::MaxCoeff(diag)/ScalarType(resolution);
|
||||
if ( (diag[0]<=0.)
|
||||
|| (diag[1]<=0.)
|
||||
|| (diag[2]<=0.)
|
||||
|| (resolution==0))
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
unsigned int nofCells[3];
|
||||
unsigned int nofBlocks[3];
|
||||
mCache = new GridElement[(mMaxBlockSize)*(mMaxBlockSize)*(mMaxBlockSize)];
|
||||
ScalarType step = vcg::math::Max(diag[0],diag[1],diag[2])/ScalarType(resolution);
|
||||
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
nofCells[k] = int(diag[k]/step)+2;
|
||||
nofBlocks[k] = nofCells[k]/mMaxBlockSize + ( (nofCells[k]%mMaxBlockSize)==0 ? 0 : 1);
|
||||
}
|
||||
unsigned int nofCells[3];
|
||||
unsigned int nofBlocks[3];
|
||||
|
||||
_mesh = &mesh;
|
||||
_mesh->Clear();
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
nofCells[k] = int(diag[k]/step)+2;
|
||||
nofBlocks[k] = nofCells[k]/mMaxBlockSize + ( (nofCells[k]%mMaxBlockSize)==0 ? 0 : 1);
|
||||
}
|
||||
|
||||
int countSubSlice = 0;
|
||||
int totalSubSlices = nofBlocks[2] * nofBlocks[1] * nofCells[0];
|
||||
_mesh = &mesh;
|
||||
_mesh->Clear();
|
||||
|
||||
extractor.Initialize();
|
||||
// for each macro block
|
||||
vcg::Point3i bi; // block id
|
||||
for (bi[2]=0 ; bi[2]<nofBlocks[2] ; ++bi[2])
|
||||
for (bi[1]=0 ; bi[1]<nofBlocks[1] ; ++bi[1])
|
||||
for (bi[0]=0 ; bi[0]<nofBlocks[0] ; ++bi[0])
|
||||
{
|
||||
mBlockOrigin = bi * (mMaxBlockSize-1);
|
||||
int countSubSlice = 0;
|
||||
int totalSubSlices = nofBlocks[2] * nofBlocks[1] * nofCells[0];
|
||||
|
||||
// compute the size of the local grid
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
mGridSize[k] = std::min<int>(mMaxBlockSize, nofCells[k]-(mMaxBlockSize-1)*bi[k]);
|
||||
}
|
||||
VectorType origin = mAABB.min + VectorType(bi[0],bi[1],bi[2]) * (step * (mMaxBlockSize-1));
|
||||
extractor.Initialize();
|
||||
// for each macro block
|
||||
vcg::Point3i bi; // block id
|
||||
for (bi[2]=0 ; bi[2]<nofBlocks[2] ; ++bi[2])
|
||||
for (bi[1]=0 ; bi[1]<nofBlocks[1] ; ++bi[1])
|
||||
for (bi[0]=0 ; bi[0]<nofBlocks[0] ; ++bi[0])
|
||||
{
|
||||
mBlockOrigin = bi * (mMaxBlockSize-1);
|
||||
|
||||
// fill the grid
|
||||
vcg::Point3i ci; // local cell id
|
||||
// compute the size of the local grid
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
mGridSize[k] = std::min<int>(mMaxBlockSize, nofCells[k]-(mMaxBlockSize-1)*bi[k]);
|
||||
}
|
||||
VectorType origin = mAABB.min + VectorType(bi[0],bi[1],bi[2]) * (step * (mMaxBlockSize-1));
|
||||
|
||||
// for each corners...
|
||||
for (ci[0]=0 ; ci[0]<mGridSize[0] ; ++ci[0])
|
||||
{
|
||||
if (cb)
|
||||
cb((100*(++countSubSlice))/totalSubSlices, "Marching cube...");
|
||||
for (ci[1]=0 ; ci[1]<mGridSize[1] ; ++ci[1])
|
||||
for (ci[2]=0 ; ci[2]<mGridSize[2] ; ++ci[2])
|
||||
{
|
||||
GridElement& el = mCache[(ci[2]*mMaxBlockSize + ci[1])*mMaxBlockSize + ci[0]];
|
||||
el.position = origin + VectorType(ci[0],ci[1],ci[2]) * step;
|
||||
el.value = mpSurface->potential(el.position);
|
||||
if (!mpSurface->isInDomain(el.position))
|
||||
el.value = invalidValue;
|
||||
}
|
||||
}
|
||||
// fill the grid
|
||||
vcg::Point3i ci; // local cell id
|
||||
|
||||
// polygonize the grid (marching cube)
|
||||
// for each cell...
|
||||
for (ci[0]=0 ; ci[0]<mGridSize[0]-1 ; ++ci[0])
|
||||
for (ci[1]=0 ; ci[1]<mGridSize[1]-1 ; ++ci[1])
|
||||
for (ci[2]=0 ; ci[2]<mGridSize[2]-1 ; ++ci[2])
|
||||
{
|
||||
uint cellId = ci[0]+mMaxBlockSize*(ci[1]+mMaxBlockSize*ci[2]);
|
||||
// check if one corner is outside the surface definition domain
|
||||
bool out =false;
|
||||
for (int k=0; k<8 && (!out); ++k)
|
||||
out = out || (!IsFinite(mCache[cellId+offsets[k]].value))
|
||||
|| mCache[cellId+offsets[k]].value==invalidValue;
|
||||
// for each corners...
|
||||
for (ci[0]=0 ; ci[0]<mGridSize[0] ; ++ci[0])
|
||||
{
|
||||
if (cb)
|
||||
cb((100*(++countSubSlice))/totalSubSlices, "Marching cube...");
|
||||
for (ci[1]=0 ; ci[1]<mGridSize[1] ; ++ci[1])
|
||||
for (ci[2]=0 ; ci[2]<mGridSize[2] ; ++ci[2])
|
||||
{
|
||||
GridElement& el = mCache[(ci[2]*mMaxBlockSize + ci[1])*mMaxBlockSize + ci[0]];
|
||||
el.position = origin + VectorType(ci[0],ci[1],ci[2]) * step;
|
||||
el.value = mpSurface->potential(el.position);
|
||||
if (!mpSurface->isInDomain(el.position))
|
||||
el.value = invalidValue;
|
||||
}
|
||||
}
|
||||
|
||||
if (!out)
|
||||
{
|
||||
extractor.ProcessCell(ci+mBlockOrigin, ci+mBlockOrigin+vcg::Point3i(1,1,1));
|
||||
}
|
||||
}
|
||||
}
|
||||
extractor.Finalize();
|
||||
_mesh = NULL;
|
||||
delete[] mCache;
|
||||
};
|
||||
// polygonize the grid (marching cube)
|
||||
// for each cell...
|
||||
for (ci[0]=0 ; ci[0]<mGridSize[0]-1 ; ++ci[0])
|
||||
for (ci[1]=0 ; ci[1]<mGridSize[1]-1 ; ++ci[1])
|
||||
for (ci[2]=0 ; ci[2]<mGridSize[2]-1 ; ++ci[2])
|
||||
{
|
||||
uint cellId = ci[0]+mMaxBlockSize*(ci[1]+mMaxBlockSize*ci[2]);
|
||||
// check if one corner is outside the surface definition domain
|
||||
bool out =false;
|
||||
for (int k=0; k<8 && (!out); ++k)
|
||||
out = out || (!IsFinite(mCache[cellId+offsets[k]].value))
|
||||
|| mCache[cellId+offsets[k]].value==invalidValue;
|
||||
|
||||
int GetLocalCellId(const vcg::Point3i& p)
|
||||
{
|
||||
return p[0] + (p[1] + p[2]*mMaxBlockSize)*mMaxBlockSize;
|
||||
}
|
||||
if (!out)
|
||||
{
|
||||
extractor.ProcessCell(ci+mBlockOrigin, ci+mBlockOrigin+vcg::Point3i(1,1,1));
|
||||
}
|
||||
}
|
||||
}
|
||||
extractor.Finalize();
|
||||
_mesh = NULL;
|
||||
delete[] mCache;
|
||||
};
|
||||
|
||||
int GetLocalCellIdFromGlobal(vcg::Point3i p)
|
||||
{
|
||||
p -= mBlockOrigin;
|
||||
return GetLocalCellId(p);
|
||||
}
|
||||
int GetLocalCellId(const vcg::Point3i& p)
|
||||
{
|
||||
return p[0] + (p[1] + p[2]*mMaxBlockSize)*mMaxBlockSize;
|
||||
}
|
||||
|
||||
float V(int pi, int pj, int pk)
|
||||
{
|
||||
return mCache[GetLocalCellIdFromGlobal(vcg::Point3i(pi, pj, pk))].value;
|
||||
}
|
||||
int GetLocalCellIdFromGlobal(vcg::Point3i p)
|
||||
{
|
||||
p -= mBlockOrigin;
|
||||
return GetLocalCellId(p);
|
||||
}
|
||||
|
||||
int GetGlobalCellId(const vcg::Point3i &p)
|
||||
{
|
||||
return p[0] + p[1] * resolution + p[2] * resolution * resolution;
|
||||
}
|
||||
float V(int pi, int pj, int pk)
|
||||
{
|
||||
return mCache[GetLocalCellIdFromGlobal(vcg::Point3i(pi, pj, pk))].value;
|
||||
}
|
||||
|
||||
void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v, bool create)
|
||||
{
|
||||
int id1 = GetGlobalCellId(p1);
|
||||
int id2 = GetGlobalCellId(p2);
|
||||
if (id1>id2)
|
||||
std::swap(id1,id2);
|
||||
Key k = Key(id1) + (Key(id2)<<32);
|
||||
MapType::iterator it = mVertexMap.find(k);
|
||||
if (it!=mVertexMap.end())
|
||||
{
|
||||
// a vertex already exist
|
||||
v = &_mesh->vert[it->second];
|
||||
}
|
||||
else if (create)
|
||||
{
|
||||
// let's create a new vertex
|
||||
VertexIndex vi = (VertexIndex) _mesh->vert.size();
|
||||
Allocator<MeshType>::AddVertices( *_mesh, 1 );
|
||||
mVertexMap[k] = vi;
|
||||
v = &_mesh->vert[vi];
|
||||
id1 = GetLocalCellIdFromGlobal(p1);
|
||||
id2 = GetLocalCellIdFromGlobal(p2);
|
||||
// interpol along the edge
|
||||
ScalarType epsilon = ScalarType(1e-5);
|
||||
const GridElement& c1 = mCache[id1];
|
||||
const GridElement& c2 = mCache[id2];
|
||||
if (fabs(mIsoValue-c1.value) < epsilon)
|
||||
v->P() = c1.position;
|
||||
else if (fabs(mIsoValue-c2.value) < epsilon)
|
||||
v->P() = c2.position;
|
||||
else if (fabs(c1.value-c2.value) < epsilon)
|
||||
v->P() = (c1.position+c1.position)*0.5;
|
||||
else
|
||||
{
|
||||
ScalarType a = (mIsoValue - c1.value) / (c2.value - c1.value);
|
||||
v->P() = c1.position + (c2.position - c1.position) * a;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
v = 0;
|
||||
}
|
||||
}
|
||||
int GetGlobalCellId(const vcg::Point3i &p)
|
||||
{
|
||||
return p[0] + p[1] * resolution + p[2] * resolution * resolution;
|
||||
}
|
||||
|
||||
bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p0, p1, v, false);
|
||||
return v!=0;
|
||||
}
|
||||
void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v, bool create)
|
||||
{
|
||||
int id1 = GetGlobalCellId(p1);
|
||||
int id2 = GetGlobalCellId(p2);
|
||||
if (id1>id2)
|
||||
std::swap(id1,id2);
|
||||
Key k = Key(id1) + (Key(id2)<<32);
|
||||
MapType::iterator it = mVertexMap.find(k);
|
||||
if (it!=mVertexMap.end())
|
||||
{
|
||||
// a vertex already exist
|
||||
v = &_mesh->vert[it->second];
|
||||
}
|
||||
else if (create)
|
||||
{
|
||||
// let's create a new vertex
|
||||
VertexIndex vi = (VertexIndex) _mesh->vert.size();
|
||||
Allocator<MeshType>::AddVertices( *_mesh, 1 );
|
||||
mVertexMap[k] = vi;
|
||||
v = &_mesh->vert[vi];
|
||||
id1 = GetLocalCellIdFromGlobal(p1);
|
||||
id2 = GetLocalCellIdFromGlobal(p2);
|
||||
// interpol along the edge
|
||||
ScalarType epsilon = ScalarType(1e-5);
|
||||
const GridElement& c1 = mCache[id1];
|
||||
const GridElement& c2 = mCache[id2];
|
||||
if (fabs(mIsoValue-c1.value) < epsilon)
|
||||
v->P().Import(c1.position);
|
||||
else if (fabs(mIsoValue-c2.value) < epsilon)
|
||||
v->P().Import(c2.position);
|
||||
else if (fabs(c1.value-c2.value) < epsilon)
|
||||
v->P().Import( (c1.position+c1.position)*0.5);
|
||||
else
|
||||
{
|
||||
ScalarType a = (mIsoValue - c1.value) / (c2.value - c1.value);
|
||||
v->P().Import(c1.position + (c2.position - c1.position) * a);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
v = 0;
|
||||
}
|
||||
}
|
||||
|
||||
void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p0, p1, v, false);
|
||||
return v!=0;
|
||||
}
|
||||
|
||||
void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v)
|
||||
{
|
||||
GetIntercept(p1, p2, v, true);
|
||||
}
|
||||
|
||||
protected:
|
||||
Box3f mAABB;
|
||||
MapType mVertexMap;
|
||||
MeshType *_mesh;
|
||||
SurfaceType* mpSurface;
|
||||
GridElement* mCache;
|
||||
vcg::Point3i mBlockOrigin;
|
||||
vcg::Point3i mGridSize;
|
||||
int mMaxBlockSize;
|
||||
ScalarType mIsoValue;
|
||||
Box3m mAABB;
|
||||
MapType mVertexMap;
|
||||
MeshType *_mesh;
|
||||
SurfaceType* mpSurface;
|
||||
GridElement* mCache;
|
||||
vcg::Point3i mBlockOrigin;
|
||||
vcg::Point3i mGridSize;
|
||||
int mMaxBlockSize;
|
||||
ScalarType mIsoValue;
|
||||
|
||||
};
|
||||
} // end namespace
|
||||
|
||||
@ -52,114 +52,114 @@ using namespace vcg;
|
||||
// - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly
|
||||
|
||||
enum {
|
||||
CT_MEAN = 0,
|
||||
CT_GAUSS = 1,
|
||||
CT_K1 = 2,
|
||||
CT_K2 = 3,
|
||||
CT_APSS = 4
|
||||
CT_MEAN = 0,
|
||||
CT_GAUSS = 1,
|
||||
CT_K1 = 2,
|
||||
CT_K2 = 3,
|
||||
CT_APSS = 4
|
||||
};
|
||||
|
||||
MlsPlugin::MlsPlugin()
|
||||
{
|
||||
typeList
|
||||
<< FP_RIMLS_PROJECTION << FP_APSS_PROJECTION
|
||||
typeList
|
||||
<< FP_RIMLS_PROJECTION << FP_APSS_PROJECTION
|
||||
// << FP_RIMLS_AFRONT << FP_APSS_AFRONT
|
||||
<< FP_RIMLS_MCUBE << FP_APSS_MCUBE
|
||||
<< FP_RIMLS_COLORIZE << FP_APSS_COLORIZE
|
||||
<< FP_RADIUS_FROM_DENSITY
|
||||
<< FP_SELECT_SMALL_COMPONENTS;
|
||||
<< FP_RIMLS_MCUBE << FP_APSS_MCUBE
|
||||
<< FP_RIMLS_COLORIZE << FP_APSS_COLORIZE
|
||||
<< FP_RADIUS_FROM_DENSITY
|
||||
<< FP_SELECT_SMALL_COMPONENTS;
|
||||
|
||||
// initFilterList(this);
|
||||
foreach(FilterIDType tt , types())
|
||||
actionList << new QAction(filterName(tt), this);
|
||||
foreach(FilterIDType tt , types())
|
||||
actionList << new QAction(filterName(tt), this);
|
||||
}
|
||||
|
||||
// ST() must return the very short string describing each filtering action
|
||||
// (this string is used also to define the menu entry)
|
||||
QString MlsPlugin::filterName(FilterIDType filterId) const
|
||||
{
|
||||
switch(filterId) {
|
||||
case FP_APSS_PROJECTION : return QString("MLS projection (APSS)");
|
||||
case FP_RIMLS_PROJECTION : return QString("MLS projection (RIMLS)");
|
||||
case FP_APSS_AFRONT : return QString("MLS meshing/APSS Advancing Front");
|
||||
case FP_RIMLS_AFRONT : return QString("MLS meshing/RIMLS Advancing Front");
|
||||
case FP_APSS_MCUBE : return QString("Marching Cubes (APSS)");
|
||||
case FP_RIMLS_MCUBE : return QString("Marching Cubes (RIMLS)");
|
||||
case FP_APSS_COLORIZE : return QString("Colorize curvature (APSS)");
|
||||
case FP_RIMLS_COLORIZE : return QString("Colorize curvature (RIMLS)");
|
||||
case FP_RADIUS_FROM_DENSITY : return QString("Estimate radius from density");
|
||||
case FP_SELECT_SMALL_COMPONENTS : return QString("Small component selection");
|
||||
default : assert(0);
|
||||
}
|
||||
switch(filterId) {
|
||||
case FP_APSS_PROJECTION : return QString("MLS projection (APSS)");
|
||||
case FP_RIMLS_PROJECTION : return QString("MLS projection (RIMLS)");
|
||||
case FP_APSS_AFRONT : return QString("MLS meshing/APSS Advancing Front");
|
||||
case FP_RIMLS_AFRONT : return QString("MLS meshing/RIMLS Advancing Front");
|
||||
case FP_APSS_MCUBE : return QString("Marching Cubes (APSS)");
|
||||
case FP_RIMLS_MCUBE : return QString("Marching Cubes (RIMLS)");
|
||||
case FP_APSS_COLORIZE : return QString("Colorize curvature (APSS)");
|
||||
case FP_RIMLS_COLORIZE : return QString("Colorize curvature (RIMLS)");
|
||||
case FP_RADIUS_FROM_DENSITY : return QString("Estimate radius from density");
|
||||
case FP_SELECT_SMALL_COMPONENTS : return QString("Small component selection");
|
||||
default : assert(0);
|
||||
}
|
||||
return QString("Filter Unknown");
|
||||
}
|
||||
|
||||
MeshFilterInterface::FilterClass MlsPlugin::getClass(QAction *a)
|
||||
{
|
||||
int filterId = ID(a);
|
||||
int filterId = ID(a);
|
||||
|
||||
switch(filterId) {
|
||||
case FP_APSS_PROJECTION :
|
||||
case FP_RIMLS_PROJECTION : return FilterClass(MeshFilterInterface::PointSet + MeshFilterInterface::Smoothing);
|
||||
case FP_APSS_AFRONT :
|
||||
case FP_RIMLS_AFRONT :
|
||||
case FP_APSS_MCUBE :
|
||||
case FP_RIMLS_MCUBE : return FilterClass(MeshFilterInterface::PointSet | MeshFilterInterface::Remeshing);
|
||||
case FP_APSS_COLORIZE :
|
||||
case FP_RIMLS_COLORIZE : return FilterClass(MeshFilterInterface::PointSet | MeshFilterInterface::VertexColoring);
|
||||
case FP_RADIUS_FROM_DENSITY : return MeshFilterInterface::PointSet;
|
||||
case FP_SELECT_SMALL_COMPONENTS : return MeshFilterInterface::Selection;
|
||||
}
|
||||
assert(0);
|
||||
return MeshFilterInterface::Generic;
|
||||
switch(filterId) {
|
||||
case FP_APSS_PROJECTION :
|
||||
case FP_RIMLS_PROJECTION : return FilterClass(MeshFilterInterface::PointSet + MeshFilterInterface::Smoothing);
|
||||
case FP_APSS_AFRONT :
|
||||
case FP_RIMLS_AFRONT :
|
||||
case FP_APSS_MCUBE :
|
||||
case FP_RIMLS_MCUBE : return FilterClass(MeshFilterInterface::PointSet | MeshFilterInterface::Remeshing);
|
||||
case FP_APSS_COLORIZE :
|
||||
case FP_RIMLS_COLORIZE : return FilterClass(MeshFilterInterface::PointSet | MeshFilterInterface::VertexColoring);
|
||||
case FP_RADIUS_FROM_DENSITY : return MeshFilterInterface::PointSet;
|
||||
case FP_SELECT_SMALL_COMPONENTS : return MeshFilterInterface::Selection;
|
||||
}
|
||||
assert(0);
|
||||
return MeshFilterInterface::Generic;
|
||||
}
|
||||
|
||||
// Info() must return the longer string describing each filtering action
|
||||
// (this string is used in the About plugin dialog)
|
||||
QString MlsPlugin::filterInfo(FilterIDType filterId) const
|
||||
{
|
||||
QString str = "";
|
||||
if (filterId & _PROJECTION_)
|
||||
{
|
||||
str += "Project a mesh (or a point set) onto the MLS surface defined by itself or another point set.<br>";
|
||||
}
|
||||
QString str = "";
|
||||
if (filterId & _PROJECTION_)
|
||||
{
|
||||
str += "Project a mesh (or a point set) onto the MLS surface defined by itself or another point set.<br>";
|
||||
}
|
||||
|
||||
if (filterId & _MCUBE_)
|
||||
{
|
||||
str +=
|
||||
"Extract the iso-surface (as a mesh) of a MLS surface defined by the current point set (or mesh)"
|
||||
"using the marching cubes algorithm. The coarse extraction is followed by an accurate projection"
|
||||
"step onto the MLS, and an extra zero removal procedure.<br>";
|
||||
}
|
||||
if (filterId & _MCUBE_)
|
||||
{
|
||||
str +=
|
||||
"Extract the iso-surface (as a mesh) of a MLS surface defined by the current point set (or mesh)"
|
||||
"using the marching cubes algorithm. The coarse extraction is followed by an accurate projection"
|
||||
"step onto the MLS, and an extra zero removal procedure.<br>";
|
||||
}
|
||||
|
||||
if (filterId & _COLORIZE_)
|
||||
{
|
||||
str += "Colorize the vertices of a mesh or point set using the curfvature of the underlying surface.<br>";
|
||||
}
|
||||
if (filterId & _COLORIZE_)
|
||||
{
|
||||
str += "Colorize the vertices of a mesh or point set using the curfvature of the underlying surface.<br>";
|
||||
}
|
||||
|
||||
if (filterId & _APSS_)
|
||||
{
|
||||
str +=
|
||||
"<br>This is the <i>algebraic point set surfaces</i> (APSS) variant which is based on "
|
||||
"the local fitting of algebraic spheres. It requires points equipped with oriented normals. <br>"
|
||||
"For all the details about APSS see: <br> Guennebaud and Gross, 'Algebraic Point Set Surfaces', Siggraph 2007, and<br>"
|
||||
"Guennebaud et al., 'Dynamic Sampling and Rendering of APSS', Eurographics 2008";
|
||||
}
|
||||
if (filterId & _APSS_)
|
||||
{
|
||||
str +=
|
||||
"<br>This is the <i>algebraic point set surfaces</i> (APSS) variant which is based on "
|
||||
"the local fitting of algebraic spheres. It requires points equipped with oriented normals. <br>"
|
||||
"For all the details about APSS see: <br> Guennebaud and Gross, 'Algebraic Point Set Surfaces', Siggraph 2007, and<br>"
|
||||
"Guennebaud et al., 'Dynamic Sampling and Rendering of APSS', Eurographics 2008";
|
||||
}
|
||||
|
||||
if (filterId & _RIMLS_)
|
||||
{
|
||||
str +=
|
||||
"<br>This is the Robust Implicit MLS (RIMLS) variant which is an extension of "
|
||||
"Implicit MLS preserving sharp features using non linear regression. For more details see: <br>"
|
||||
"Oztireli, Guennebaud and Gross, 'Feature Preserving Point Set Surfaces based on Non-Linear Kernel Regression' Eurographics 2009.";
|
||||
}
|
||||
if (filterId & _RIMLS_)
|
||||
{
|
||||
str +=
|
||||
"<br>This is the Robust Implicit MLS (RIMLS) variant which is an extension of "
|
||||
"Implicit MLS preserving sharp features using non linear regression. For more details see: <br>"
|
||||
"Oztireli, Guennebaud and Gross, 'Feature Preserving Point Set Surfaces based on Non-Linear Kernel Regression' Eurographics 2009.";
|
||||
}
|
||||
|
||||
if (filterId == FP_RADIUS_FROM_DENSITY)
|
||||
str = "Estimate the local point spacing (aka radius) around each vertex using a basic estimate of the local density.";
|
||||
else if (filterId == FP_SELECT_SMALL_COMPONENTS)
|
||||
str = "Select the small disconnected components of a mesh.";
|
||||
if (filterId == FP_RADIUS_FROM_DENSITY)
|
||||
str = "Estimate the local point spacing (aka radius) around each vertex using a basic estimate of the local density.";
|
||||
else if (filterId == FP_SELECT_SMALL_COMPONENTS)
|
||||
str = "Select the small disconnected components of a mesh.";
|
||||
|
||||
return str;
|
||||
return str;
|
||||
}
|
||||
|
||||
|
||||
@ -172,155 +172,155 @@ return QString("Filter Unknown");
|
||||
// - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog)
|
||||
void MlsPlugin::initParameterSet(QAction* action, MeshDocument& md, RichParameterSet& parlst)
|
||||
{
|
||||
int id = ID(action);
|
||||
MeshModel *target = md.mm();
|
||||
int id = ID(action);
|
||||
MeshModel *target = md.mm();
|
||||
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
parlst.addParam(new RichFloat("NbFaceRatio",
|
||||
0.1f,
|
||||
"Small component ratio",
|
||||
"This ratio (between 0 and 1) defines the meaning of <i>small</i> as the threshold ratio between the number of faces"
|
||||
"of the largest component and the other ones. A larger value will select more components."));
|
||||
parlst.addParam(new RichBool( "NonClosedOnly",
|
||||
false,
|
||||
"Select only non closed components",
|
||||
""));
|
||||
return;
|
||||
}
|
||||
else if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
parlst.addParam(new RichInt("NbNeighbors",
|
||||
16,
|
||||
"Number of neighbors",
|
||||
"Number of neighbors used to estimate the local density. Larger values lead to smoother variations."));
|
||||
return;
|
||||
}
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
parlst.addParam(new RichFloat("NbFaceRatio",
|
||||
0.1f,
|
||||
"Small component ratio",
|
||||
"This ratio (between 0 and 1) defines the meaning of <i>small</i> as the threshold ratio between the number of faces"
|
||||
"of the largest component and the other ones. A larger value will select more components."));
|
||||
parlst.addParam(new RichBool( "NonClosedOnly",
|
||||
false,
|
||||
"Select only non closed components",
|
||||
""));
|
||||
return;
|
||||
}
|
||||
else if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
parlst.addParam(new RichInt("NbNeighbors",
|
||||
16,
|
||||
"Number of neighbors",
|
||||
"Number of neighbors used to estimate the local density. Larger values lead to smoother variations."));
|
||||
return;
|
||||
}
|
||||
|
||||
if ((id & _PROJECTION_))
|
||||
{
|
||||
parlst.addParam(new RichMesh( "ControlMesh", target,&md, "Point set",
|
||||
"The point set (or mesh) which defines the MLS surface."));
|
||||
parlst.addParam(new RichMesh( "ProxyMesh", target, &md, "Proxy Mesh",
|
||||
"The mesh that will be projected/resampled onto the MLS surface."));
|
||||
}
|
||||
if ((id & _PROJECTION_) || (id & _COLORIZE_))
|
||||
{
|
||||
parlst.addParam(new RichBool( "SelectionOnly",
|
||||
target->cm.sfn>0,
|
||||
"Selection only",
|
||||
"If checked, only selected vertices will be projected."));
|
||||
}
|
||||
if ((id & _PROJECTION_))
|
||||
{
|
||||
parlst.addParam(new RichMesh( "ControlMesh", target,&md, "Point set",
|
||||
"The point set (or mesh) which defines the MLS surface."));
|
||||
parlst.addParam(new RichMesh( "ProxyMesh", target, &md, "Proxy Mesh",
|
||||
"The mesh that will be projected/resampled onto the MLS surface."));
|
||||
}
|
||||
if ((id & _PROJECTION_) || (id & _COLORIZE_))
|
||||
{
|
||||
parlst.addParam(new RichBool( "SelectionOnly",
|
||||
target->cm.sfn>0,
|
||||
"Selection only",
|
||||
"If checked, only selected vertices will be projected."));
|
||||
}
|
||||
|
||||
if ( (id & _APSS_) || (id & _RIMLS_) )
|
||||
{
|
||||
parlst.addParam(new RichFloat("FilterScale",
|
||||
2.0,
|
||||
"MLS - Filter scale",
|
||||
"Scale of the spatial low pass filter.\n"
|
||||
"It is relative to the radius (local point spacing) of the vertices."));
|
||||
parlst.addParam(new RichFloat("ProjectionAccuracy",
|
||||
1e-4f,
|
||||
"Projection - Accuracy (adv)",
|
||||
"Threshold value used to stop the projections.\n"
|
||||
"This value is scaled by the mean point spacing to get the actual threshold."));
|
||||
parlst.addParam(new RichInt( "MaxProjectionIters",
|
||||
15,
|
||||
"Projection - Max iterations (adv)",
|
||||
"Max number of iterations for the projection."));
|
||||
}
|
||||
if ( (id & _APSS_) || (id & _RIMLS_) )
|
||||
{
|
||||
parlst.addParam(new RichFloat("FilterScale",
|
||||
2.0,
|
||||
"MLS - Filter scale",
|
||||
"Scale of the spatial low pass filter.\n"
|
||||
"It is relative to the radius (local point spacing) of the vertices."));
|
||||
parlst.addParam(new RichFloat("ProjectionAccuracy",
|
||||
1e-4f,
|
||||
"Projection - Accuracy (adv)",
|
||||
"Threshold value used to stop the projections.\n"
|
||||
"This value is scaled by the mean point spacing to get the actual threshold."));
|
||||
parlst.addParam(new RichInt( "MaxProjectionIters",
|
||||
15,
|
||||
"Projection - Max iterations (adv)",
|
||||
"Max number of iterations for the projection."));
|
||||
}
|
||||
|
||||
if (id & _APSS_)
|
||||
{
|
||||
parlst.addParam(new RichFloat("SphericalParameter",
|
||||
1,
|
||||
"MLS - Spherical parameter",
|
||||
"Control the curvature of the fitted spheres: 0 is equivalent to a pure plane fit,"
|
||||
"1 to a pure spherical fit, values between 0 and 1 gives intermediate results,"
|
||||
"while others real values might give interresting results, but take care with extreme"
|
||||
"settings !"));
|
||||
if (!(id & _COLORIZE_))
|
||||
parlst.addParam(new RichBool( "AccurateNormal",
|
||||
true,
|
||||
"Accurate normals",
|
||||
"If checked, use the accurate MLS gradient instead of the local approximation"
|
||||
"to compute the normals."));
|
||||
}
|
||||
if (id & _APSS_)
|
||||
{
|
||||
parlst.addParam(new RichFloat("SphericalParameter",
|
||||
1,
|
||||
"MLS - Spherical parameter",
|
||||
"Control the curvature of the fitted spheres: 0 is equivalent to a pure plane fit,"
|
||||
"1 to a pure spherical fit, values between 0 and 1 gives intermediate results,"
|
||||
"while others real values might give interresting results, but take care with extreme"
|
||||
"settings !"));
|
||||
if (!(id & _COLORIZE_))
|
||||
parlst.addParam(new RichBool( "AccurateNormal",
|
||||
true,
|
||||
"Accurate normals",
|
||||
"If checked, use the accurate MLS gradient instead of the local approximation"
|
||||
"to compute the normals."));
|
||||
}
|
||||
|
||||
if (id & _RIMLS_)
|
||||
{
|
||||
parlst.addParam(new RichFloat("SigmaN",
|
||||
0.75,
|
||||
"MLS - Sharpness",
|
||||
"Width of the filter used by the normal refitting weight."
|
||||
"This weight function is a Gaussian on the distance between two unit vectors:"
|
||||
"the current gradient and the input normal. Therefore, typical value range between 0.5 (sharp) to 2 (smooth)."));
|
||||
parlst.addParam(new RichInt( "MaxRefittingIters",
|
||||
3,
|
||||
"MLS - Max fitting iterations",
|
||||
"Max number of fitting iterations. (0 or 1 is equivalent to the standard IMLS)"));
|
||||
}
|
||||
if (id & _RIMLS_)
|
||||
{
|
||||
parlst.addParam(new RichFloat("SigmaN",
|
||||
0.75,
|
||||
"MLS - Sharpness",
|
||||
"Width of the filter used by the normal refitting weight."
|
||||
"This weight function is a Gaussian on the distance between two unit vectors:"
|
||||
"the current gradient and the input normal. Therefore, typical value range between 0.5 (sharp) to 2 (smooth)."));
|
||||
parlst.addParam(new RichInt( "MaxRefittingIters",
|
||||
3,
|
||||
"MLS - Max fitting iterations",
|
||||
"Max number of fitting iterations. (0 or 1 is equivalent to the standard IMLS)"));
|
||||
}
|
||||
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
parlst.addParam(new RichInt( "MaxSubdivisions",
|
||||
0,
|
||||
"Refinement - Max subdivisions",
|
||||
"Max number of subdivisions."));
|
||||
parlst.addParam(new RichFloat("ThAngleInDegree",
|
||||
2,
|
||||
"Refinement - Crease angle (degree)",
|
||||
"Threshold angle between two faces controlling the refinement."));
|
||||
}
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
parlst.addParam(new RichInt( "MaxSubdivisions",
|
||||
0,
|
||||
"Refinement - Max subdivisions",
|
||||
"Max number of subdivisions."));
|
||||
parlst.addParam(new RichFloat("ThAngleInDegree",
|
||||
2,
|
||||
"Refinement - Crease angle (degree)",
|
||||
"Threshold angle between two faces controlling the refinement."));
|
||||
}
|
||||
|
||||
if (id & _AFRONT_)
|
||||
{
|
||||
}
|
||||
if (id & _AFRONT_)
|
||||
{
|
||||
}
|
||||
|
||||
if ((id & _COLORIZE_))
|
||||
{
|
||||
QStringList lst;
|
||||
lst << "Mean" << "Gauss" << "K1" << "K2";
|
||||
if (id & _APSS_)
|
||||
lst << "ApproxMean";
|
||||
if ((id & _COLORIZE_))
|
||||
{
|
||||
QStringList lst;
|
||||
lst << "Mean" << "Gauss" << "K1" << "K2";
|
||||
if (id & _APSS_)
|
||||
lst << "ApproxMean";
|
||||
|
||||
parlst.addParam(new RichEnum("CurvatureType", CT_MEAN,
|
||||
lst,
|
||||
"Curvature type",
|
||||
QString("The type of the curvature to plot.")
|
||||
+ ((id & _APSS_) ? "<br>ApproxMean uses the radius of the fitted sphere as an approximation of the mean curvature." : "")));
|
||||
parlst.addParam(new RichEnum("CurvatureType", CT_MEAN,
|
||||
lst,
|
||||
"Curvature type",
|
||||
QString("The type of the curvature to plot.")
|
||||
+ ((id & _APSS_) ? "<br>ApproxMean uses the radius of the fitted sphere as an approximation of the mean curvature." : "")));
|
||||
// if ((id & _APSS_))
|
||||
// parlst.addParam(new RichBool( "ApproxCurvature",
|
||||
// false,
|
||||
// "Approx mean curvature",
|
||||
// "If checked, use the radius of the fitted sphere as an approximation of the mean curvature.");
|
||||
}
|
||||
}
|
||||
|
||||
if (id & _MCUBE_)
|
||||
{
|
||||
parlst.addParam(new RichInt( "Resolution",
|
||||
200,
|
||||
"Grid Resolution",
|
||||
"The resolution of the grid on which we run the marching cubes."
|
||||
"This marching cube is memory friendly, so you can safely set large values up to 1000 or even more."));
|
||||
}
|
||||
if (id & _MCUBE_)
|
||||
{
|
||||
parlst.addParam(new RichInt( "Resolution",
|
||||
200,
|
||||
"Grid Resolution",
|
||||
"The resolution of the grid on which we run the marching cubes."
|
||||
"This marching cube is memory friendly, so you can safely set large values up to 1000 or even more."));
|
||||
}
|
||||
}
|
||||
|
||||
int MlsPlugin::getRequirements(QAction *action)
|
||||
{
|
||||
return 0;
|
||||
return 0;
|
||||
}
|
||||
|
||||
/** Predicate functor for adaptive refinement according to crease angle.
|
||||
*
|
||||
*/
|
||||
*
|
||||
*/
|
||||
template <class MESH_TYPE, typename Scalar>
|
||||
struct EdgeAnglePredicate
|
||||
{
|
||||
Scalar thCosAngle;
|
||||
bool operator()(vcg::face::Pos<typename MESH_TYPE::FaceType> ep) const
|
||||
{
|
||||
// FIXME why does the following fails:
|
||||
Scalar thCosAngle;
|
||||
bool operator()(vcg::face::Pos<typename MESH_TYPE::FaceType> ep) const
|
||||
{
|
||||
// FIXME why does the following fails:
|
||||
// vcg::face::Pos<typename MESH_TYPE::FaceType> op = ep;
|
||||
// op.FlipF();
|
||||
// if (op.f)
|
||||
@ -328,288 +328,288 @@ struct EdgeAnglePredicate
|
||||
// else
|
||||
// return true;
|
||||
|
||||
return vcg::Dot(ep.F()->cN(), ep.FFlip()->cN()) < thCosAngle;
|
||||
}
|
||||
return (ep.F()->cN() * ep.FFlip()->cN()) < thCosAngle;
|
||||
}
|
||||
};
|
||||
|
||||
/** compute the normal of a face as the average of its vertices */
|
||||
template<typename MeshType>
|
||||
void UpdateFaceNormalFromVertex(MeshType& m)
|
||||
{
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename VertexType::NormalType NormalType;
|
||||
typedef typename VertexType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename VertexType::NormalType NormalType;
|
||||
typedef typename VertexType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
|
||||
for (FaceIterator f=m.face.begin(); f!=m.face.end(); ++f)
|
||||
{
|
||||
NormalType n;
|
||||
n.SetZero();
|
||||
for(int j=0; j<3; ++j)
|
||||
n += f->V(j)->cN();
|
||||
n.Normalize();
|
||||
f->N() = n;
|
||||
}
|
||||
for (FaceIterator f=m.face.begin(); f!=m.face.end(); ++f)
|
||||
{
|
||||
NormalType n;
|
||||
n.SetZero();
|
||||
for(int j=0; j<3; ++j)
|
||||
n += f->V(j)->cN();
|
||||
n.Normalize();
|
||||
f->N() = n;
|
||||
}
|
||||
}
|
||||
|
||||
bool MlsPlugin::applyFilter(QAction* filter, MeshDocument& md, RichParameterSet& par, vcg::CallBackPos* cb)
|
||||
{
|
||||
int id = ID(filter);
|
||||
int id = ID(filter);
|
||||
|
||||
if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
md.mm()->updateDataMask(MeshModel::MM_VERTRADIUS);
|
||||
APSS<CMeshO> mls(md.mm()->cm);
|
||||
mls.computeVertexRaddi(par.getInt("NbNeighbors"));
|
||||
return true;
|
||||
}
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
MeshModel* mesh = md.mm();
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
bool nonClosedOnly = par.getBool("NonClosedOnly");
|
||||
float ratio = par.getFloat("NbFaceRatio");
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, ratio, nonClosedOnly);
|
||||
return true;
|
||||
}
|
||||
if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
md.mm()->updateDataMask(MeshModel::MM_VERTRADIUS);
|
||||
APSS<CMeshO> mls(md.mm()->cm);
|
||||
mls.computeVertexRaddi(par.getInt("NbNeighbors"));
|
||||
return true;
|
||||
}
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
MeshModel* mesh = md.mm();
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
bool nonClosedOnly = par.getBool("NonClosedOnly");
|
||||
float ratio = par.getFloat("NbFaceRatio");
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, ratio, nonClosedOnly);
|
||||
return true;
|
||||
}
|
||||
|
||||
// we are doing some MLS based stuff
|
||||
{
|
||||
if(md.mm()->cm.fn > 0)
|
||||
{ // if we start from a mesh, and it has unreferenced vertices
|
||||
// normals are undefined on that vertices.
|
||||
int delvert=tri::Clean<CMeshO>::RemoveUnreferencedVertex(md.mm()->cm);
|
||||
if(delvert) Log( "Pre-MLS Cleaning: Removed %d unreferenced vertices",delvert);
|
||||
}
|
||||
tri::Allocator<CMeshO>::CompactVertexVector(md.mm()->cm);
|
||||
// we are doing some MLS based stuff
|
||||
{
|
||||
if(md.mm()->cm.fn > 0)
|
||||
{ // if we start from a mesh, and it has unreferenced vertices
|
||||
// normals are undefined on that vertices.
|
||||
int delvert=tri::Clean<CMeshO>::RemoveUnreferencedVertex(md.mm()->cm);
|
||||
if(delvert) Log( "Pre-MLS Cleaning: Removed %d unreferenced vertices",delvert);
|
||||
}
|
||||
tri::Allocator<CMeshO>::CompactVertexVector(md.mm()->cm);
|
||||
|
||||
// We require a per vertex radius so as a first thing check it
|
||||
if(!md.mm()->hasDataMask(MeshModel::MM_VERTRADIUS))
|
||||
{
|
||||
md.mm()->updateDataMask(MeshModel::MM_VERTRADIUS);
|
||||
APSS<CMeshO> mls(md.mm()->cm);
|
||||
mls.computeVertexRaddi();
|
||||
Log( "Mesh has no per vertex radius. Computed and added using default neighbourhood");
|
||||
}
|
||||
// We require a per vertex radius so as a first thing check it
|
||||
if(!md.mm()->hasDataMask(MeshModel::MM_VERTRADIUS))
|
||||
{
|
||||
md.mm()->updateDataMask(MeshModel::MM_VERTRADIUS);
|
||||
APSS<CMeshO> mls(md.mm()->cm);
|
||||
mls.computeVertexRaddi();
|
||||
Log( "Mesh has no per vertex radius. Computed and added using default neighbourhood");
|
||||
}
|
||||
|
||||
MeshModel* pPoints = 0;
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
if (par.getMesh("ControlMesh") == par.getMesh("ProxyMesh"))
|
||||
{
|
||||
// clone the control mesh
|
||||
MeshModel* ref = par.getMesh("ControlMesh");
|
||||
pPoints = md.addNewMesh("","TempMesh");
|
||||
pPoints->updateDataMask(ref);
|
||||
vcg::tri::Append<CMeshO,CMeshO>::Mesh(pPoints->cm, ref->cm); // the last true means "copy all vertices"
|
||||
vcg::tri::UpdateBounding<CMeshO>::Box(pPoints->cm);
|
||||
pPoints->cm.Tr = ref->cm.Tr;
|
||||
}
|
||||
else
|
||||
pPoints = par.getMesh("ControlMesh");
|
||||
}
|
||||
else // for curvature
|
||||
pPoints = md.mm();
|
||||
MeshModel* pPoints = 0;
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
if (par.getMesh("ControlMesh") == par.getMesh("ProxyMesh"))
|
||||
{
|
||||
// clone the control mesh
|
||||
MeshModel* ref = par.getMesh("ControlMesh");
|
||||
pPoints = md.addNewMesh("","TempMesh");
|
||||
pPoints->updateDataMask(ref);
|
||||
vcg::tri::Append<CMeshO,CMeshO>::Mesh(pPoints->cm, ref->cm); // the last true means "copy all vertices"
|
||||
vcg::tri::UpdateBounding<CMeshO>::Box(pPoints->cm);
|
||||
pPoints->cm.Tr = ref->cm.Tr;
|
||||
}
|
||||
else
|
||||
pPoints = par.getMesh("ControlMesh");
|
||||
}
|
||||
else // for curvature
|
||||
pPoints = md.mm();
|
||||
|
||||
// create the MLS surface
|
||||
cb(1, "Create the MLS data structures...");
|
||||
MlsSurface<CMeshO>* mls = 0;
|
||||
// create the MLS surface
|
||||
cb(1, "Create the MLS data structures...");
|
||||
MlsSurface<CMeshO>* mls = 0;
|
||||
|
||||
RIMLS<CMeshO>* rimls = 0;
|
||||
APSS<CMeshO>* apss = 0;
|
||||
RIMLS<CMeshO>* rimls = 0;
|
||||
APSS<CMeshO>* apss = 0;
|
||||
|
||||
if (id & _RIMLS_)
|
||||
mls = rimls = new RIMLS<CMeshO>(pPoints->cm);
|
||||
else if (id & _APSS_)
|
||||
mls = apss = new APSS<CMeshO>(pPoints->cm);
|
||||
else
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
if (id & _RIMLS_)
|
||||
mls = rimls = new RIMLS<CMeshO>(pPoints->cm);
|
||||
else if (id & _APSS_)
|
||||
mls = apss = new APSS<CMeshO>(pPoints->cm);
|
||||
else
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
mls->setFilterScale(par.getFloat("FilterScale"));
|
||||
mls->setMaxProjectionIters(par.getInt("MaxProjectionIters"));
|
||||
mls->setProjectionAccuracy(par.getFloat("ProjectionAccuracy"));
|
||||
mls->setFilterScale(par.getFloat("FilterScale"));
|
||||
mls->setMaxProjectionIters(par.getInt("MaxProjectionIters"));
|
||||
mls->setProjectionAccuracy(par.getFloat("ProjectionAccuracy"));
|
||||
|
||||
if (rimls)
|
||||
{
|
||||
rimls->setMaxRefittingIters(par.getInt("MaxRefittingIters"));
|
||||
//mls.setMinRefittingIters(par.getFloat("MinRefittingIters"));
|
||||
rimls->setSigmaN(par.getFloat("SigmaN"));
|
||||
}
|
||||
if (rimls)
|
||||
{
|
||||
rimls->setMaxRefittingIters(par.getInt("MaxRefittingIters"));
|
||||
//mls.setMinRefittingIters(par.getFloat("MinRefittingIters"));
|
||||
rimls->setSigmaN(par.getFloat("SigmaN"));
|
||||
}
|
||||
|
||||
if (apss)
|
||||
{
|
||||
apss->setSphericalParameter(par.getFloat("SphericalParameter"));
|
||||
if (!(id & _COLORIZE_))
|
||||
apss->setGradientHint(par.getBool("AccurateNormal") ? GaelMls::MLS_DERIVATIVE_ACCURATE : GaelMls::MLS_DERIVATIVE_APPROX);
|
||||
}
|
||||
if (apss)
|
||||
{
|
||||
apss->setSphericalParameter(par.getFloat("SphericalParameter"));
|
||||
if (!(id & _COLORIZE_))
|
||||
apss->setGradientHint(par.getBool("AccurateNormal") ? GaelMls::MLS_DERIVATIVE_ACCURATE : GaelMls::MLS_DERIVATIVE_APPROX);
|
||||
}
|
||||
|
||||
MeshModel * mesh = 0;
|
||||
MeshModel * mesh = 0;
|
||||
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
mesh = par.getMesh("ProxyMesh");
|
||||
bool selectionOnly = par.getBool("SelectionOnly");
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
mesh = par.getMesh("ProxyMesh");
|
||||
bool selectionOnly = par.getBool("SelectionOnly");
|
||||
|
||||
if (selectionOnly)
|
||||
vcg::tri::UpdateSelection<CMeshO>::VertexFromFaceStrict(mesh->cm);
|
||||
EdgeAnglePredicate<CMeshO,float> edgePred;
|
||||
edgePred.thCosAngle = cos(M_PI * par.getFloat("ThAngleInDegree")/180.);
|
||||
if (selectionOnly)
|
||||
vcg::tri::UpdateSelection<CMeshO>::VertexFromFaceStrict(mesh->cm);
|
||||
EdgeAnglePredicate<CMeshO,float> edgePred;
|
||||
edgePred.thCosAngle = cos(M_PI * par.getFloat("ThAngleInDegree")/180.);
|
||||
|
||||
int nbRefinements = par.getInt("MaxSubdivisions");
|
||||
for (int k=0; k<nbRefinements+1; ++k)
|
||||
{
|
||||
//UpdateFaceNormalFromVertex(m.cm);
|
||||
if (k!=0)
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
int nbRefinements = par.getInt("MaxSubdivisions");
|
||||
for (int k=0; k<nbRefinements+1; ++k)
|
||||
{
|
||||
//UpdateFaceNormalFromVertex(m.cm);
|
||||
if (k!=0)
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
|
||||
vcg::tri::UpdateNormal<CMeshO>::PerFace(mesh->cm);
|
||||
vcg::tri::UpdateNormal<CMeshO>::NormalizePerFace(mesh->cm);
|
||||
//vcg::RefineE<CMeshO,vcg::MidPoint<CMeshO> >(m.cm, vcg::MidPoint<CMeshO>(), edgePred, false, cb);
|
||||
vcg::tri::RefineOddEvenE<CMeshO, tri::OddPointLoop<CMeshO>, tri::EvenPointLoop<CMeshO> >
|
||||
(mesh->cm, tri::OddPointLoop<CMeshO>(mesh->cm), tri::EvenPointLoop<CMeshO>(), edgePred, selectionOnly, cb);
|
||||
}
|
||||
// project all vertices onto the MLS surface
|
||||
for (unsigned int i = 0; i< mesh->cm.vert.size(); i++)
|
||||
{
|
||||
cb(1+98*i/mesh->cm.vert.size(), "MLS projection...");
|
||||
vcg::tri::UpdateNormal<CMeshO>::PerFace(mesh->cm);
|
||||
vcg::tri::UpdateNormal<CMeshO>::NormalizePerFace(mesh->cm);
|
||||
//vcg::RefineE<CMeshO,vcg::MidPoint<CMeshO> >(m.cm, vcg::MidPoint<CMeshO>(), edgePred, false, cb);
|
||||
vcg::tri::RefineOddEvenE<CMeshO, tri::OddPointLoop<CMeshO>, tri::EvenPointLoop<CMeshO> >
|
||||
(mesh->cm, tri::OddPointLoop<CMeshO>(mesh->cm), tri::EvenPointLoop<CMeshO>(), edgePred, selectionOnly, cb);
|
||||
}
|
||||
// project all vertices onto the MLS surface
|
||||
for (unsigned int i = 0; i< mesh->cm.vert.size(); i++)
|
||||
{
|
||||
cb(1+98*i/mesh->cm.vert.size(), "MLS projection...");
|
||||
|
||||
if ( (!selectionOnly) || (mesh->cm.vert[i].IsS()) )
|
||||
mesh->cm.vert[i].P() = mls->project(mesh->cm.vert[i].P(), &mesh->cm.vert[i].N());
|
||||
}
|
||||
}
|
||||
if ( (!selectionOnly) || (mesh->cm.vert[i].IsS()) )
|
||||
mesh->cm.vert[i].P() = mls->project(mesh->cm.vert[i].P(), &mesh->cm.vert[i].N());
|
||||
}
|
||||
}
|
||||
|
||||
Log( "Successfully projected %i vertices", mesh->cm.vn);
|
||||
}
|
||||
else if (id & _COLORIZE_)
|
||||
{
|
||||
mesh = md.mm();
|
||||
mesh->updateDataMask(MeshModel::MM_VERTCOLOR);
|
||||
mesh->updateDataMask(MeshModel::MM_VERTQUALITY);
|
||||
mesh->updateDataMask(MeshModel::MM_VERTCURVDIR);
|
||||
Log( "Successfully projected %i vertices", mesh->cm.vn);
|
||||
}
|
||||
else if (id & _COLORIZE_)
|
||||
{
|
||||
mesh = md.mm();
|
||||
mesh->updateDataMask(MeshModel::MM_VERTCOLOR);
|
||||
mesh->updateDataMask(MeshModel::MM_VERTQUALITY);
|
||||
mesh->updateDataMask(MeshModel::MM_VERTCURVDIR);
|
||||
|
||||
bool selectionOnly = par.getBool("SelectionOnly");
|
||||
//bool approx = apss && par.getBool("ApproxCurvature");
|
||||
int ct = par.getEnum("CurvatureType");
|
||||
bool selectionOnly = par.getBool("SelectionOnly");
|
||||
//bool approx = apss && par.getBool("ApproxCurvature");
|
||||
int ct = par.getEnum("CurvatureType");
|
||||
|
||||
uint size = mesh->cm.vert.size();
|
||||
//std::vector<float> curvatures(size);
|
||||
float minc=1e9, maxc=-1e9, minabsc=1e9;
|
||||
vcg::Point3f grad;
|
||||
vcg::Matrix33f hess;
|
||||
uint size = mesh->cm.vert.size();
|
||||
//std::vector<float> curvatures(size);
|
||||
float minc=1e9, maxc=-1e9, minabsc=1e9;
|
||||
Point3m grad;
|
||||
Matrix33m hess;
|
||||
|
||||
// pass 1: computes curvatures and statistics
|
||||
for (unsigned int i = 0; i< size; i++)
|
||||
{
|
||||
cb(1+98*i/pPoints->cm.vert.size(), "MLS colorization...");
|
||||
// pass 1: computes curvatures and statistics
|
||||
for (unsigned int i = 0; i< size; i++)
|
||||
{
|
||||
cb(1+98*i/pPoints->cm.vert.size(), "MLS colorization...");
|
||||
|
||||
if ( (!selectionOnly) || (pPoints->cm.vert[i].IsS()) )
|
||||
{
|
||||
Point3f p = mls->project(mesh->cm.vert[i].P());
|
||||
float c = 0;
|
||||
if ( (!selectionOnly) || (pPoints->cm.vert[i].IsS()) )
|
||||
{
|
||||
Point3m p = mls->project(mesh->cm.vert[i].P());
|
||||
float c = 0;
|
||||
|
||||
if (ct==CT_APSS)
|
||||
c = apss->approxMeanCurvature(p);
|
||||
else
|
||||
{
|
||||
int errorMask;
|
||||
grad = mls->gradient(p, &errorMask);
|
||||
if (errorMask == MLS_OK && grad.Norm() > 1e-8)
|
||||
{
|
||||
hess = mls->hessian(p);
|
||||
implicits::WeingartenMap<float> W(grad,hess);
|
||||
if (ct==CT_APSS)
|
||||
c = apss->approxMeanCurvature(p);
|
||||
else
|
||||
{
|
||||
int errorMask;
|
||||
grad = mls->gradient(p, &errorMask);
|
||||
if (errorMask == MLS_OK && grad.Norm() > 1e-8)
|
||||
{
|
||||
hess = mls->hessian(p);
|
||||
implicits::WeingartenMap<CMeshO::ScalarType> W(grad,hess);
|
||||
|
||||
mesh->cm.vert[i].PD1() = W.K1Dir();
|
||||
mesh->cm.vert[i].PD2() = W.K2Dir();
|
||||
mesh->cm.vert[i].K1() = W.K1();
|
||||
mesh->cm.vert[i].K2() = W.K2();
|
||||
mesh->cm.vert[i].PD1() = W.K1Dir();
|
||||
mesh->cm.vert[i].PD2() = W.K2Dir();
|
||||
mesh->cm.vert[i].K1() = W.K1();
|
||||
mesh->cm.vert[i].K2() = W.K2();
|
||||
|
||||
switch(ct)
|
||||
{
|
||||
case CT_MEAN: c = W.MeanCurvature(); break;
|
||||
case CT_GAUSS: c = W.GaussCurvature(); break;
|
||||
case CT_K1: c = W.K1(); break;
|
||||
case CT_K2: c = W.K2(); break;
|
||||
default: assert(0 && "invalid curvature type");
|
||||
}
|
||||
}
|
||||
assert(!math::IsNAN(c) && "You should never try to compute Histogram with Invalid Floating points numbers (NaN)");
|
||||
}
|
||||
mesh->cm.vert[i].Q() = c;
|
||||
minc = std::min(c,minc);
|
||||
maxc = std::max(c,maxc);
|
||||
minabsc = std::min(fabsf(c),minabsc);
|
||||
}
|
||||
}
|
||||
// pass 2: convert the curvature to color
|
||||
cb(99, "Curvature to color...");
|
||||
float d = maxc-minc;
|
||||
minc += 0.05*d;
|
||||
maxc -= 0.05*d;
|
||||
switch(ct)
|
||||
{
|
||||
case CT_MEAN: c = W.MeanCurvature(); break;
|
||||
case CT_GAUSS: c = W.GaussCurvature(); break;
|
||||
case CT_K1: c = W.K1(); break;
|
||||
case CT_K2: c = W.K2(); break;
|
||||
default: assert(0 && "invalid curvature type");
|
||||
}
|
||||
}
|
||||
assert(!math::IsNAN(c) && "You should never try to compute Histogram with Invalid Floating points numbers (NaN)");
|
||||
}
|
||||
mesh->cm.vert[i].Q() = c;
|
||||
minc = std::min(c,minc);
|
||||
maxc = std::max(c,maxc);
|
||||
minabsc = std::min(fabsf(c),minabsc);
|
||||
}
|
||||
}
|
||||
// pass 2: convert the curvature to color
|
||||
cb(99, "Curvature to color...");
|
||||
float d = maxc-minc;
|
||||
minc += 0.05*d;
|
||||
maxc -= 0.05*d;
|
||||
|
||||
vcg::Histogramf H;
|
||||
vcg::Histogramf H;
|
||||
vcg::tri::Stat<CMeshO>::ComputePerVertexQualityHistogram(mesh->cm,H);
|
||||
vcg::tri::UpdateColor<CMeshO>::PerVertexQualityRamp(mesh->cm,H.Percentile(0.01f),H.Percentile(0.99f));
|
||||
}
|
||||
// else if (id & _AFRONT_)
|
||||
// {
|
||||
// // create a new mesh
|
||||
// mesh = md.addNewMesh("afront mesh");
|
||||
// vcg::tri::AdvancingMLS<CMeshO> afront(mesh->cm, *mls);
|
||||
// //afront.BuildMesh(cb);
|
||||
// afront._SeedFace();
|
||||
// for (int i=0; i<20120; ++i)
|
||||
// afront.AddFace();
|
||||
// Log(0, "Advancing front MLS meshing done.");
|
||||
// }
|
||||
else if (id & _MCUBE_)
|
||||
{
|
||||
// create a new mesh
|
||||
mesh = md.addNewMesh("","mc_mesh");
|
||||
}
|
||||
// else if (id & _AFRONT_)
|
||||
// {
|
||||
// // create a new mesh
|
||||
// mesh = md.addNewMesh("afront mesh");
|
||||
// vcg::tri::AdvancingMLS<CMeshO> afront(mesh->cm, *mls);
|
||||
// //afront.BuildMesh(cb);
|
||||
// afront._SeedFace();
|
||||
// for (int i=0; i<20120; ++i)
|
||||
// afront.AddFace();
|
||||
// Log(0, "Advancing front MLS meshing done.");
|
||||
// }
|
||||
else if (id & _MCUBE_)
|
||||
{
|
||||
// create a new mesh
|
||||
mesh = md.addNewMesh("","mc_mesh");
|
||||
|
||||
typedef vcg::tri::MlsWalker<CMeshO,MlsSurface<CMeshO> > MlsWalker;
|
||||
typedef vcg::tri::MarchingCubes<CMeshO, MlsWalker> MlsMarchingCubes;
|
||||
MlsWalker walker;
|
||||
walker.resolution = par.getInt("Resolution");
|
||||
typedef vcg::tri::MlsWalker<CMeshO,MlsSurface<CMeshO> > MlsWalker;
|
||||
typedef vcg::tri::MarchingCubes<CMeshO, MlsWalker> MlsMarchingCubes;
|
||||
MlsWalker walker;
|
||||
walker.resolution = par.getInt("Resolution");
|
||||
|
||||
// iso extraction
|
||||
MlsMarchingCubes mc(mesh->cm, walker);
|
||||
walker.BuildMesh<MlsMarchingCubes>(mesh->cm, *mls, mc, cb);
|
||||
// iso extraction
|
||||
MlsMarchingCubes mc(mesh->cm, walker);
|
||||
walker.BuildMesh<MlsMarchingCubes>(mesh->cm, *mls, mc, cb);
|
||||
|
||||
// accurate projection
|
||||
for (unsigned int i = 0; i< mesh->cm.vert.size(); i++)
|
||||
{
|
||||
cb(1+98*i/mesh->cm.vert.size(), "MLS projection...");
|
||||
mesh->cm.vert[i].P() = mls->project(mesh->cm.vert[i].P(), &mesh->cm.vert[i].N());
|
||||
}
|
||||
// accurate projection
|
||||
for (unsigned int i = 0; i< mesh->cm.vert.size(); i++)
|
||||
{
|
||||
cb(1+98*i/mesh->cm.vert.size(), "MLS projection...");
|
||||
mesh->cm.vert[i].P() = mls->project(mesh->cm.vert[i].P(), &mesh->cm.vert[i].N());
|
||||
}
|
||||
|
||||
// extra zero detection and removal
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO );
|
||||
// selection...
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, 0.1f);
|
||||
// deletion...
|
||||
vcg::tri::SmallComponent<CMeshO>::DeleteFaceVert(mesh->cm);
|
||||
mesh->clearDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
}
|
||||
// extra zero detection and removal
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO );
|
||||
// selection...
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, 0.1f);
|
||||
// deletion...
|
||||
vcg::tri::SmallComponent<CMeshO>::DeleteFaceVert(mesh->cm);
|
||||
mesh->clearDataMask(MeshModel::MM_FACEFACETOPO);
|
||||
}
|
||||
|
||||
Log( "Marching cubes MLS meshing done.");
|
||||
}
|
||||
Log( "Marching cubes MLS meshing done.");
|
||||
}
|
||||
|
||||
delete mls;
|
||||
if ( (id & _PROJECTION_) && par.getMesh("ControlMesh")!=pPoints)
|
||||
{
|
||||
delete mls;
|
||||
if ( (id & _PROJECTION_) && par.getMesh("ControlMesh")!=pPoints)
|
||||
{
|
||||
md.delMesh(pPoints);
|
||||
}
|
||||
}
|
||||
|
||||
if (mesh)
|
||||
mesh->UpdateBoxAndNormals();
|
||||
if (mesh)
|
||||
mesh->UpdateBoxAndNormals();
|
||||
|
||||
} // end MLS based stuff
|
||||
} // end MLS based stuff
|
||||
|
||||
return true;
|
||||
return true;
|
||||
}
|
||||
|
||||
MESHLAB_PLUGIN_NAME_EXPORTER(MlsPlugin)
|
||||
|
||||
@ -33,168 +33,168 @@
|
||||
namespace GaelMls {
|
||||
|
||||
enum {
|
||||
MLS_OK,
|
||||
MLS_TOO_FAR,
|
||||
MLS_TOO_MANY_ITERS,
|
||||
MLS_NOT_SUPPORTED,
|
||||
MLS_OK,
|
||||
MLS_TOO_FAR,
|
||||
MLS_TOO_MANY_ITERS,
|
||||
MLS_NOT_SUPPORTED,
|
||||
|
||||
MLS_DERIVATIVE_ACCURATE,
|
||||
MLS_DERIVATIVE_APPROX,
|
||||
MLS_DERIVATIVE_FINITEDIFF
|
||||
MLS_DERIVATIVE_ACCURATE,
|
||||
MLS_DERIVATIVE_APPROX,
|
||||
MLS_DERIVATIVE_FINITEDIFF
|
||||
};
|
||||
|
||||
template<typename MeshType>
|
||||
class MlsSurface
|
||||
{
|
||||
public:
|
||||
typedef typename MeshType::ScalarType Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
typedef vcg::Matrix33<Scalar> MatrixType;
|
||||
typedef typename MeshType::VertContainer PointsType;
|
||||
public:
|
||||
typedef typename MeshType::ScalarType Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
typedef vcg::Matrix33<Scalar> MatrixType;
|
||||
typedef typename MeshType::VertContainer PointsType;
|
||||
|
||||
MlsSurface(const MeshType& mesh)
|
||||
: mMesh(mesh), mPoints(mesh.vert)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
MlsSurface(const MeshType& mesh)
|
||||
: mMesh(mesh), mPoints(mesh.vert)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
|
||||
mAABB = mesh.bbox;
|
||||
mAABB = mesh.bbox;
|
||||
|
||||
// compute radii using a basic meshless density estimator
|
||||
if (!mPoints.RadiusEnabled)
|
||||
{
|
||||
const_cast<PointsType&>(mPoints).EnableRadius();
|
||||
computeVertexRaddi();
|
||||
}
|
||||
// compute radii using a basic meshless density estimator
|
||||
if (!mPoints.RadiusEnabled)
|
||||
{
|
||||
const_cast<PointsType&>(mPoints).EnableRadius();
|
||||
computeVertexRaddi();
|
||||
}
|
||||
|
||||
mFilterScale = 4.0;
|
||||
mMaxNofProjectionIterations = 20;
|
||||
mProjectionAccuracy = (Scalar)1e-4;
|
||||
mBallTree = 0;
|
||||
mGradientHint = MLS_DERIVATIVE_ACCURATE;
|
||||
mHessianHint = MLS_DERIVATIVE_ACCURATE;
|
||||
mFilterScale = 4.0;
|
||||
mMaxNofProjectionIterations = 20;
|
||||
mProjectionAccuracy = (Scalar)1e-4;
|
||||
mBallTree = 0;
|
||||
mGradientHint = MLS_DERIVATIVE_ACCURATE;
|
||||
mHessianHint = MLS_DERIVATIVE_ACCURATE;
|
||||
|
||||
mDomainMinNofNeighbors = 4;
|
||||
mDomainRadiusScale = 2.;
|
||||
mDomainNormalScale = 1.;
|
||||
}
|
||||
mDomainMinNofNeighbors = 4;
|
||||
mDomainRadiusScale = 2.;
|
||||
mDomainNormalScale = 1.;
|
||||
}
|
||||
|
||||
/** \returns the value of the reconstructed scalar field at point \a x */
|
||||
virtual Scalar potential(const VectorType& x, int* errorMask = 0) const = 0;
|
||||
/** \returns the value of the reconstructed scalar field at point \a x */
|
||||
virtual Scalar potential(const VectorType& x, int* errorMask = 0) const = 0;
|
||||
|
||||
/** \returns the gradient of the reconstructed scalar field at point \a x
|
||||
*
|
||||
* The method used to compute the gradient can be controlled with setGradientHint().
|
||||
*/
|
||||
virtual VectorType gradient(const VectorType& x, int* errorMask = 0) const = 0;
|
||||
/** \returns the gradient of the reconstructed scalar field at point \a x
|
||||
*
|
||||
* The method used to compute the gradient can be controlled with setGradientHint().
|
||||
*/
|
||||
virtual VectorType gradient(const VectorType& x, int* errorMask = 0) const = 0;
|
||||
|
||||
/** \returns the hessian matrix of the reconstructed scalar field at point \a x
|
||||
*
|
||||
* The method used to compute the hessian matrix can be controlled with setHessianHint().
|
||||
*/
|
||||
virtual MatrixType hessian(const VectorType& x, int* errorMask = 0) const
|
||||
{ if (errorMask) *errorMask = MLS_NOT_SUPPORTED; return MatrixType(); }
|
||||
/** \returns the hessian matrix of the reconstructed scalar field at point \a x
|
||||
*
|
||||
* The method used to compute the hessian matrix can be controlled with setHessianHint().
|
||||
*/
|
||||
virtual MatrixType hessian(const VectorType& x, int* errorMask = 0) const
|
||||
{ if (errorMask) *errorMask = MLS_NOT_SUPPORTED; return MatrixType(); }
|
||||
|
||||
/** \returns the projection of point x onto the MLS surface, and optionnaly returns the normal in \a pNormal */
|
||||
virtual VectorType project(const VectorType& x, VectorType* pNormal = 0, int* errorMask = 0) const = 0;
|
||||
/** \returns the projection of point x onto the MLS surface, and optionnaly returns the normal in \a pNormal */
|
||||
virtual VectorType project(const VectorType& x, VectorType* pNormal = 0, int* errorMask = 0) const = 0;
|
||||
|
||||
/** \returns whether \a x is inside the restricted surface definition domain */
|
||||
virtual bool isInDomain(const VectorType& x) const;
|
||||
/** \returns whether \a x is inside the restricted surface definition domain */
|
||||
virtual bool isInDomain(const VectorType& x) const;
|
||||
|
||||
/** \returns the mean curvature from the gradient vector and Hessian matrix.
|
||||
*/
|
||||
Scalar meanCurvature(const VectorType& gradient, const MatrixType& hessian) const;
|
||||
/** \returns the mean curvature from the gradient vector and Hessian matrix.
|
||||
*/
|
||||
Scalar meanCurvature(const VectorType& gradient, const MatrixType& hessian) const;
|
||||
|
||||
/** set the scale of the spatial filter */
|
||||
void setFilterScale(Scalar v);
|
||||
/** set the maximum number of iterations during the projection */
|
||||
void setMaxProjectionIters(int n);
|
||||
/** set the threshold factor to detect convergence of the iterations */
|
||||
void setProjectionAccuracy(Scalar v);
|
||||
/** set the scale of the spatial filter */
|
||||
void setFilterScale(Scalar v);
|
||||
/** set the maximum number of iterations during the projection */
|
||||
void setMaxProjectionIters(int n);
|
||||
/** set the threshold factor to detect convergence of the iterations */
|
||||
void setProjectionAccuracy(Scalar v);
|
||||
|
||||
/** set a hint on how to compute the gradient
|
||||
*
|
||||
* Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF
|
||||
*/
|
||||
void setGradientHint(int h);
|
||||
/** set a hint on how to compute the gradient
|
||||
*
|
||||
* Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF
|
||||
*/
|
||||
void setGradientHint(int h);
|
||||
|
||||
/** set a hint on how to compute the hessian matrix
|
||||
*
|
||||
* Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF
|
||||
*/
|
||||
void setHessianHint(int h);
|
||||
/** set a hint on how to compute the hessian matrix
|
||||
*
|
||||
* Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF
|
||||
*/
|
||||
void setHessianHint(int h);
|
||||
|
||||
inline const MeshType& mesh() const { return mMesh; }
|
||||
/** a shortcut for mesh().vert */
|
||||
inline const PointsType& points() const { return mPoints; }
|
||||
inline const MeshType& mesh() const { return mMesh; }
|
||||
/** a shortcut for mesh().vert */
|
||||
inline const PointsType& points() const { return mPoints; }
|
||||
|
||||
inline ConstDataWrapper<VectorType> positions() const
|
||||
{
|
||||
return ConstDataWrapper<VectorType>(&mPoints[0].P(), mPoints.size(),
|
||||
size_t(mPoints[1].P().V()) - size_t(mPoints[0].P().V()));
|
||||
}
|
||||
inline ConstDataWrapper<VectorType> normals() const
|
||||
{
|
||||
return ConstDataWrapper<VectorType>(&mPoints[0].N(), mPoints.size(),
|
||||
size_t(mPoints[1].N().V()) - size_t(mPoints[0].N().V()));
|
||||
}
|
||||
inline ConstDataWrapper<Scalar> radii() const
|
||||
{
|
||||
return ConstDataWrapper<Scalar>(&mPoints[0].R(), mPoints.size(),
|
||||
size_t(&mPoints[1].R()) - size_t(&mPoints[0].R()));
|
||||
}
|
||||
const vcg::Box3<Scalar>& boundingBox() const { return mAABB; }
|
||||
inline vcg::ConstDataWrapper<VectorType> positions() const
|
||||
{
|
||||
return vcg::ConstDataWrapper<VectorType>(&mPoints[0].P(), mPoints.size(),
|
||||
size_t(mPoints[1].P().V()) - size_t(mPoints[0].P().V()));
|
||||
}
|
||||
inline vcg::ConstDataWrapper<VectorType> normals() const
|
||||
{
|
||||
return vcg::ConstDataWrapper<VectorType>(&mPoints[0].N(), mPoints.size(),
|
||||
size_t(mPoints[1].N().V()) - size_t(mPoints[0].N().V()));
|
||||
}
|
||||
inline vcg::ConstDataWrapper<Scalar> radii() const
|
||||
{
|
||||
return vcg::ConstDataWrapper<Scalar>(&mPoints[0].R(), mPoints.size(),
|
||||
size_t(&mPoints[1].R()) - size_t(&mPoints[0].R()));
|
||||
}
|
||||
const vcg::Box3<Scalar>& boundingBox() const { return mAABB; }
|
||||
|
||||
static const Scalar InvalidValue() { return Scalar(12345679810.11121314151617); }
|
||||
static const Scalar InvalidValue() { return Scalar(12345679810.11121314151617); }
|
||||
|
||||
void computeVertexRaddi(const int nbNeighbors = 16);
|
||||
protected:
|
||||
void computeNeighborhood(const VectorType& x, bool computeDerivatives) const;
|
||||
void requestSecondDerivatives() const;
|
||||
void computeVertexRaddi(const int nbNeighbors = 16);
|
||||
protected:
|
||||
void computeNeighborhood(const VectorType& x, bool computeDerivatives) const;
|
||||
void requestSecondDerivatives() const;
|
||||
|
||||
struct PointToPointSqDist
|
||||
{
|
||||
inline bool operator()(const VectorType &a, const VectorType &b, Scalar& refD2, VectorType &q) const
|
||||
{
|
||||
struct PointToPointSqDist
|
||||
{
|
||||
inline bool operator()(const VectorType &a, const VectorType &b, Scalar& refD2, VectorType &q) const
|
||||
{
|
||||
// std::cout << a.X() << a.Y() << a.Z() << " - " << b.X() << b.Y() << b.Z() <<
|
||||
// " => " << vcg::Distance(a, b) << " < " << refD2 << "\n";
|
||||
Scalar d2 = vcg::SquaredDistance(a, b);
|
||||
if (d2>refD2)
|
||||
return false;
|
||||
Scalar d2 = vcg::SquaredDistance(a, b);
|
||||
if (d2>refD2)
|
||||
return false;
|
||||
|
||||
refD2 = d2;
|
||||
q = a;
|
||||
return true;
|
||||
}
|
||||
};
|
||||
refD2 = d2;
|
||||
q = a;
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
class DummyObjectMarker {};
|
||||
class DummyObjectMarker {};
|
||||
|
||||
protected:
|
||||
const MeshType& mMesh;
|
||||
const PointsType& mPoints;
|
||||
vcg::Box3<Scalar> mAABB;
|
||||
int mGradientHint;
|
||||
int mHessianHint;
|
||||
protected:
|
||||
const MeshType& mMesh;
|
||||
const PointsType& mPoints;
|
||||
vcg::Box3<Scalar> mAABB;
|
||||
int mGradientHint;
|
||||
int mHessianHint;
|
||||
|
||||
BallTree<Scalar>* mBallTree;
|
||||
BallTree<Scalar>* mBallTree;
|
||||
|
||||
int mMaxNofProjectionIterations;
|
||||
Scalar mFilterScale;
|
||||
Scalar mAveragePointSpacing;
|
||||
Scalar mProjectionAccuracy;
|
||||
int mMaxNofProjectionIterations;
|
||||
Scalar mFilterScale;
|
||||
Scalar mAveragePointSpacing;
|
||||
Scalar mProjectionAccuracy;
|
||||
|
||||
int mDomainMinNofNeighbors;
|
||||
float mDomainRadiusScale;
|
||||
float mDomainNormalScale;
|
||||
int mDomainMinNofNeighbors;
|
||||
float mDomainRadiusScale;
|
||||
float mDomainNormalScale;
|
||||
|
||||
// cached values:
|
||||
mutable bool mCachedQueryPointIsOK;
|
||||
mutable VectorType mCachedQueryPoint;
|
||||
mutable Neighborhood<Scalar> mNeighborhood;
|
||||
mutable std::vector<Scalar> mCachedWeights;
|
||||
mutable std::vector<Scalar> mCachedWeightDerivatives;
|
||||
mutable std::vector<VectorType> mCachedWeightGradients;
|
||||
mutable std::vector<Scalar> mCachedWeightSecondDerivatives;
|
||||
// cached values:
|
||||
mutable bool mCachedQueryPointIsOK;
|
||||
mutable VectorType mCachedQueryPoint;
|
||||
mutable Neighborhood<Scalar> mNeighborhood;
|
||||
mutable std::vector<Scalar> mCachedWeights;
|
||||
mutable std::vector<Scalar> mCachedWeightDerivatives;
|
||||
mutable std::vector<VectorType> mCachedWeightGradients;
|
||||
mutable std::vector<Scalar> mCachedWeightSecondDerivatives;
|
||||
};
|
||||
|
||||
} // namespace
|
||||
|
||||
@ -58,201 +58,202 @@ namespace GaelMls {
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setFilterScale(Scalar v)
|
||||
{
|
||||
mFilterScale = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
if (mBallTree)
|
||||
mBallTree->setRadiusScale(mFilterScale);
|
||||
mFilterScale = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
if (mBallTree)
|
||||
mBallTree->setRadiusScale(mFilterScale);
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setMaxProjectionIters(int n)
|
||||
{
|
||||
mMaxNofProjectionIterations = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mMaxNofProjectionIterations = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setProjectionAccuracy(Scalar v)
|
||||
{
|
||||
mProjectionAccuracy = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mProjectionAccuracy = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setGradientHint(int h)
|
||||
{
|
||||
mGradientHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mGradientHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setHessianHint(int h)
|
||||
{
|
||||
mHessianHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mHessianHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::computeVertexRaddi(const int nbNeighbors)
|
||||
{
|
||||
#if 0
|
||||
int nbNeighbors = 16;
|
||||
vcg::Octree<VectorType, Scalar> knn;
|
||||
knn.Set(mPoints.begin(), mPoints.end());
|
||||
std::vector<VectorType*> nearest_objects;
|
||||
std::vector<VectorType> nearest_points;
|
||||
std::vector<Scalar> sqDistances;
|
||||
mAveragePointSpacing = 0;
|
||||
for (uint i = 0; i< mPoints.size(); i++)
|
||||
{
|
||||
DummyObjectMarker dom;
|
||||
PointToPointSqDist dfunc;
|
||||
Scalar max_dist2 = 1e9;//std::numeric_limits<Scalar>::max();
|
||||
knn.GetKClosest(dfunc, dom, nbNeighbors, mPoints[i],
|
||||
max_dist2, nearest_objects, sqDistances, nearest_points);
|
||||
#if 0
|
||||
int nbNeighbors = 16;
|
||||
vcg::Octree<VectorType, Scalar> knn;
|
||||
knn.Set(mPoints.begin(), mPoints.end());
|
||||
std::vector<VectorType*> nearest_objects;
|
||||
std::vector<VectorType> nearest_points;
|
||||
std::vector<Scalar> sqDistances;
|
||||
mAveragePointSpacing = 0;
|
||||
for (uint i = 0; i< mPoints.size(); i++)
|
||||
{
|
||||
DummyObjectMarker dom;
|
||||
PointToPointSqDist dfunc;
|
||||
Scalar max_dist2 = 1e9;//std::numeric_limits<Scalar>::max();
|
||||
knn.GetKClosest(dfunc, dom, nbNeighbors, mPoints[i],
|
||||
max_dist2, nearest_objects, sqDistances, nearest_points);
|
||||
// for (int j=0; i<sqDistances.size(); ++j)
|
||||
// std::cout << sqDistances[j] << " ";
|
||||
// std::cout << "$\n";
|
||||
mRadii[i] = 2. * sqrt(sqDistances.at(0)/nearest_objects.size());
|
||||
mAveragePointSpacing += mRadii[i];
|
||||
}
|
||||
mAveragePointSpacing /= Scalar(mPoints.size());
|
||||
mRadii[i] = 2. * sqrt(sqDistances.at(0)/nearest_objects.size());
|
||||
mAveragePointSpacing += mRadii[i];
|
||||
}
|
||||
mAveragePointSpacing /= Scalar(mPoints.size());
|
||||
|
||||
#else
|
||||
#else
|
||||
|
||||
// int nbNeighbors = 16;
|
||||
|
||||
assert(mPoints.size()>=2);
|
||||
KdTree<Scalar> knn(positions());
|
||||
assert(mPoints.size()>=2);
|
||||
vcg::KdTree<Scalar> knn(positions());
|
||||
typename vcg::KdTree<Scalar>::PriorityQueue pq;
|
||||
// knn.setMaxNofNeighbors(nbNeighbors);
|
||||
mAveragePointSpacing = 0;
|
||||
for (size_t i = 0; i< mPoints.size(); i++)
|
||||
{
|
||||
knn.doQueryK(mPoints[i].cP(),nbNeighbors,pq);
|
||||
const_cast<PointsType&>(mPoints)[i].R() = 2. * sqrt(pq.getTopWeight()/Scalar(pq.getNofElements()));
|
||||
mAveragePointSpacing += mPoints[i].cR();
|
||||
}
|
||||
mAveragePointSpacing /= Scalar(mPoints.size());
|
||||
|
||||
knn.setMaxNofNeighbors(nbNeighbors);
|
||||
mAveragePointSpacing = 0;
|
||||
for (size_t i = 0; i< mPoints.size(); i++)
|
||||
{
|
||||
knn.doQueryK(mPoints[i].cP());
|
||||
const_cast<PointsType&>(mPoints)[i].R() = 2. * sqrt(knn.getNeighborSquaredDistance(0)/Scalar(knn.getNofFoundNeighbors()));
|
||||
mAveragePointSpacing += mPoints[i].cR();
|
||||
}
|
||||
mAveragePointSpacing /= Scalar(mPoints.size());
|
||||
|
||||
#endif
|
||||
#endif
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::computeNeighborhood(const VectorType& x, bool computeDerivatives) const
|
||||
{
|
||||
if (!mBallTree)
|
||||
{
|
||||
const_cast<BallTree<Scalar>*&>(mBallTree) = new BallTree<Scalar>(positions(), radii());
|
||||
const_cast<BallTree<Scalar>*>(mBallTree)->setRadiusScale(mFilterScale);
|
||||
}
|
||||
mBallTree->computeNeighbors(x, &mNeighborhood);
|
||||
size_t nofSamples = mNeighborhood.size();
|
||||
if (!mBallTree)
|
||||
{
|
||||
const_cast<BallTree<Scalar>*&>(mBallTree) = new BallTree<Scalar>(positions(), radii());
|
||||
const_cast<BallTree<Scalar>*>(mBallTree)->setRadiusScale(mFilterScale);
|
||||
}
|
||||
mBallTree->computeNeighbors(x, &mNeighborhood);
|
||||
size_t nofSamples = mNeighborhood.size();
|
||||
|
||||
// compute spatial weights and partial derivatives
|
||||
mCachedWeights.resize(nofSamples);
|
||||
if (computeDerivatives)
|
||||
{
|
||||
mCachedWeightDerivatives.resize(nofSamples);
|
||||
mCachedWeightGradients.resize(nofSamples);
|
||||
}
|
||||
else
|
||||
mCachedWeightGradients.clear();
|
||||
// compute spatial weights and partial derivatives
|
||||
mCachedWeights.resize(nofSamples);
|
||||
if (computeDerivatives)
|
||||
{
|
||||
mCachedWeightDerivatives.resize(nofSamples);
|
||||
mCachedWeightGradients.resize(nofSamples);
|
||||
}
|
||||
else
|
||||
mCachedWeightGradients.clear();
|
||||
|
||||
for (size_t i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar s = 1./(mPoints[id].cR()*mFilterScale);
|
||||
s = s*s;
|
||||
Scalar w = Scalar(1) - mNeighborhood.squaredDistance(i) * s;
|
||||
if (w<0)
|
||||
w = 0;
|
||||
Scalar aux = w;
|
||||
w = w * w;
|
||||
w = w * w;
|
||||
mCachedWeights[i] = w;
|
||||
for (size_t i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar s = 1./(mPoints[id].cR()*mFilterScale);
|
||||
s = s*s;
|
||||
Scalar w = Scalar(1) - mNeighborhood.squaredDistance(i) * s;
|
||||
if (w<0)
|
||||
w = 0;
|
||||
Scalar aux = w;
|
||||
w = w * w;
|
||||
w = w * w;
|
||||
mCachedWeights[i] = w;
|
||||
|
||||
if (computeDerivatives)
|
||||
{
|
||||
mCachedWeightDerivatives[i] = (-2. * s) * (4. * aux * aux * aux);
|
||||
mCachedWeightGradients[i] = (x - mPoints[id].cP()) * mCachedWeightDerivatives[i];
|
||||
}
|
||||
}
|
||||
if (computeDerivatives)
|
||||
{
|
||||
mCachedWeightDerivatives[i] = (-2. * s) * (4. * aux * aux * aux);
|
||||
mCachedWeightGradients[i] = (x - mPoints[id].cP()) * mCachedWeightDerivatives[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::requestSecondDerivatives() const
|
||||
{
|
||||
//if (!mSecondDerivativeUptodate)
|
||||
{
|
||||
size_t nofSamples = mNeighborhood.size();
|
||||
if (nofSamples>mCachedWeightSecondDerivatives.size())
|
||||
mCachedWeightSecondDerivatives.resize(nofSamples+10);
|
||||
//if (!mSecondDerivativeUptodate)
|
||||
{
|
||||
size_t nofSamples = mNeighborhood.size();
|
||||
if (nofSamples>mCachedWeightSecondDerivatives.size())
|
||||
mCachedWeightSecondDerivatives.resize(nofSamples+10);
|
||||
|
||||
{
|
||||
for (size_t i=0 ; i<nofSamples ; ++i)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar s = 1./(mPoints[id].cR()*mFilterScale);
|
||||
s = s*s;
|
||||
Scalar x2 = s * mNeighborhood.squaredDistance(i);
|
||||
x2 = 1.0 - x2;
|
||||
if (x2<0)
|
||||
x2 = 0.;
|
||||
mCachedWeightSecondDerivatives[i] = (4.0*s*s) * (12.0 * x2 * x2);
|
||||
}
|
||||
}
|
||||
//mSecondDerivativeUptodate = true;
|
||||
}
|
||||
{
|
||||
for (size_t i=0 ; i<nofSamples ; ++i)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar s = 1./(mPoints[id].cR()*mFilterScale);
|
||||
s = s*s;
|
||||
Scalar x2 = s * mNeighborhood.squaredDistance(i);
|
||||
x2 = 1.0 - x2;
|
||||
if (x2<0)
|
||||
x2 = 0.;
|
||||
mCachedWeightSecondDerivatives[i] = (4.0*s*s) * (12.0 * x2 * x2);
|
||||
}
|
||||
}
|
||||
//mSecondDerivativeUptodate = true;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename MlsSurface<_MeshType>::Scalar
|
||||
MlsSurface<_MeshType>::meanCurvature(const VectorType& gradient, const MatrixType& hessian) const
|
||||
{
|
||||
Scalar gl = gradient.Norm();
|
||||
return (gl*gl*hessian.Trace() - vcg::Dot(gradient, VectorType(hessian * gradient))) / (2.*gl*gl*gl);
|
||||
Scalar gl = gradient.Norm();
|
||||
// return (gl*gl*hessian.Trace() - vcg::Dot(gradient, VectorType(hessian * gradient))) / (2.*gl*gl*gl);
|
||||
return (gl*gl*hessian.Trace() - (gradient * VectorType(hessian * gradient))) / (2.*gl*gl*gl);
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool MlsSurface<_MeshType>::isInDomain(const VectorType& x) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
computeNeighborhood(x, false);
|
||||
}
|
||||
int nb = mNeighborhood.size();
|
||||
if (nb<mDomainMinNofNeighbors)
|
||||
return false;
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
computeNeighborhood(x, false);
|
||||
}
|
||||
int nb = mNeighborhood.size();
|
||||
if (nb<mDomainMinNofNeighbors)
|
||||
return false;
|
||||
|
||||
int i=0;
|
||||
bool out = true;
|
||||
bool hasNormal = true;
|
||||
if ((mDomainNormalScale==1.f) || (!hasNormal))
|
||||
{
|
||||
while (out && i<nb)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar rs2 = mPoints[id].cR() * mDomainRadiusScale;
|
||||
rs2 = rs2*rs2;
|
||||
out = mNeighborhood.squaredDistance(i) > rs2;
|
||||
++i;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
Scalar s = 1./(mDomainNormalScale*mDomainNormalScale) - 1.f;
|
||||
while (out && i<nb)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar rs2 = mPoints[id].cR() * mDomainRadiusScale;
|
||||
rs2 = rs2*rs2;
|
||||
Scalar dn = mPoints[id].cN().dot(x-mPoints[id].cP());
|
||||
out = (mNeighborhood.squaredDistance(i) + s*dn*dn) > rs2;
|
||||
++i;
|
||||
}
|
||||
}
|
||||
return !out;
|
||||
int i=0;
|
||||
bool out = true;
|
||||
bool hasNormal = true;
|
||||
if ((mDomainNormalScale==1.f) || (!hasNormal))
|
||||
{
|
||||
while (out && i<nb)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar rs2 = mPoints[id].cR() * mDomainRadiusScale;
|
||||
rs2 = rs2*rs2;
|
||||
out = mNeighborhood.squaredDistance(i) > rs2;
|
||||
++i;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
Scalar s = 1./(mDomainNormalScale*mDomainNormalScale) - 1.f;
|
||||
while (out && i<nb)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
Scalar rs2 = mPoints[id].cR() * mDomainRadiusScale;
|
||||
rs2 = rs2*rs2;
|
||||
Scalar dn = mPoints[id].cN().dot(x-mPoints[id].cP());
|
||||
out = (mNeighborhood.squaredDistance(i) + s*dn*dn) > rs2;
|
||||
++i;
|
||||
}
|
||||
}
|
||||
return !out;
|
||||
}
|
||||
|
||||
// template class MlsSurface<float>;
|
||||
|
||||
@ -43,263 +43,263 @@ namespace GaelMls {
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setSigmaR(Scalar v)
|
||||
{
|
||||
mSigmaR = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mSigmaR = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setSigmaN(Scalar v)
|
||||
{
|
||||
mSigmaN = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mSigmaN = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setRefittingThreshold(Scalar v)
|
||||
{
|
||||
mRefittingThreshold = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mRefittingThreshold = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setMinRefittingIters(int n)
|
||||
{
|
||||
mMinRefittingIters = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mMinRefittingIters = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setMaxRefittingIters(int n)
|
||||
{
|
||||
mMaxRefittingIters = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
mMaxRefittingIters = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename RIMLS<_MeshType>::Scalar RIMLS<_MeshType>::potential(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return Base::InvalidValue();
|
||||
}
|
||||
}
|
||||
|
||||
return mCachedPotential;
|
||||
return mCachedPotential;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename RIMLS<_MeshType>::VectorType RIMLS<_MeshType>::gradient(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return VectorType(0,0,0);
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return VectorType(0,0,0);
|
||||
}
|
||||
}
|
||||
|
||||
return mCachedGradient;
|
||||
return mCachedGradient;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename RIMLS<_MeshType>::MatrixType RIMLS<_MeshType>::hessian(const VectorType& x, int* errorMask) const
|
||||
{
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return MatrixType();
|
||||
}
|
||||
}
|
||||
if ((!mCachedQueryPointIsOK) || mCachedQueryPoint!=x)
|
||||
{
|
||||
if (!computePotentialAndGradient(x))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
return MatrixType();
|
||||
}
|
||||
}
|
||||
|
||||
MatrixType hessian;
|
||||
mlsHessian(x, hessian);
|
||||
return hessian;
|
||||
MatrixType hessian;
|
||||
mlsHessian(x, hessian);
|
||||
return hessian;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
typename RIMLS<_MeshType>::VectorType RIMLS<_MeshType>::project(const VectorType& x, VectorType* pNormal, int* errorMask) const
|
||||
{
|
||||
int iterationCount = 0;
|
||||
VectorType position = x;
|
||||
VectorType normal;
|
||||
Scalar delta;
|
||||
Scalar epsilon = mAveragePointSpacing * mProjectionAccuracy;
|
||||
do {
|
||||
if (!computePotentialAndGradient(position))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
//std::cerr << " proj failed\n";
|
||||
return x;
|
||||
}
|
||||
int iterationCount = 0;
|
||||
VectorType position = x;
|
||||
VectorType normal;
|
||||
Scalar delta;
|
||||
Scalar epsilon = mAveragePointSpacing * mProjectionAccuracy;
|
||||
do {
|
||||
if (!computePotentialAndGradient(position))
|
||||
{
|
||||
if (errorMask)
|
||||
*errorMask = MLS_TOO_FAR;
|
||||
//std::cerr << " proj failed\n";
|
||||
return x;
|
||||
}
|
||||
|
||||
normal = mCachedGradient;
|
||||
normal.Normalize();
|
||||
delta = mCachedPotential;
|
||||
position = position - normal*delta;
|
||||
} while ( fabs(delta)>epsilon && ++iterationCount<mMaxNofProjectionIterations);
|
||||
normal = mCachedGradient;
|
||||
normal.Normalize();
|
||||
delta = mCachedPotential;
|
||||
position = position - normal*delta;
|
||||
} while ( fabs(delta)>epsilon && ++iterationCount<mMaxNofProjectionIterations);
|
||||
|
||||
if (iterationCount>=mMaxNofProjectionIterations && errorMask)
|
||||
*errorMask = MLS_TOO_MANY_ITERS;
|
||||
if (iterationCount>=mMaxNofProjectionIterations && errorMask)
|
||||
*errorMask = MLS_TOO_MANY_ITERS;
|
||||
|
||||
if (pNormal)
|
||||
*pNormal = normal;
|
||||
if (pNormal)
|
||||
*pNormal = normal;
|
||||
|
||||
return position;
|
||||
return position;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool RIMLS<_MeshType>::computePotentialAndGradient(const VectorType& x) const
|
||||
{
|
||||
Base::computeNeighborhood(x, true);
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
Base::computeNeighborhood(x, true);
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
|
||||
if (nofSamples<1)
|
||||
{
|
||||
mCachedGradient.SetZero();
|
||||
mCachedQueryPoint = x;
|
||||
mCachedPotential = 1e9;
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
if (nofSamples<1)
|
||||
{
|
||||
mCachedGradient.SetZero();
|
||||
mCachedQueryPoint = x;
|
||||
mCachedPotential = 1e9;
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
|
||||
if (mCachedRefittingWeights.size()<nofSamples)
|
||||
mCachedRefittingWeights.resize(nofSamples+5);
|
||||
if (mCachedRefittingWeights.size()<nofSamples)
|
||||
mCachedRefittingWeights.resize(nofSamples+5);
|
||||
|
||||
VectorType source = x;
|
||||
VectorType grad; grad.SetZero();
|
||||
VectorType previousGrad;
|
||||
VectorType sumN; sumN.SetZero();
|
||||
Scalar potential = 0.;
|
||||
Scalar invSigma2 = Scalar(1) / (mSigmaN*mSigmaN);
|
||||
Scalar invSigmaR2 = 0;
|
||||
if (mSigmaR>0)
|
||||
invSigmaR2 = Scalar(1) / (mSigmaR*mSigmaR);
|
||||
VectorType sumGradWeight;
|
||||
VectorType sumGradWeightPotential;
|
||||
Scalar sumW;
|
||||
VectorType source = x;
|
||||
VectorType grad; grad.SetZero();
|
||||
VectorType previousGrad;
|
||||
VectorType sumN; sumN.SetZero();
|
||||
Scalar potential = 0.;
|
||||
Scalar invSigma2 = Scalar(1) / (mSigmaN*mSigmaN);
|
||||
Scalar invSigmaR2 = 0;
|
||||
if (mSigmaR>0)
|
||||
invSigmaR2 = Scalar(1) / (mSigmaR*mSigmaR);
|
||||
VectorType sumGradWeight;
|
||||
VectorType sumGradWeightPotential;
|
||||
Scalar sumW;
|
||||
|
||||
int iterationCount = 0;
|
||||
do
|
||||
{
|
||||
previousGrad = grad;
|
||||
sumGradWeight.SetZero();
|
||||
sumGradWeightPotential.SetZero();
|
||||
sumN.SetZero();
|
||||
potential = 0.;
|
||||
sumW = 0.;
|
||||
int iterationCount = 0;
|
||||
do
|
||||
{
|
||||
previousGrad = grad;
|
||||
sumGradWeight.SetZero();
|
||||
sumGradWeightPotential.SetZero();
|
||||
sumN.SetZero();
|
||||
potential = 0.;
|
||||
sumW = 0.;
|
||||
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
VectorType diff = source - mPoints[id].cP();
|
||||
VectorType normal = mPoints[id].cN();
|
||||
Scalar f = Dot(diff, normal);
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
VectorType diff = source - mPoints[id].cP();
|
||||
VectorType normal = mPoints[id].cN();
|
||||
Scalar f = diff *normal;
|
||||
|
||||
Scalar refittingWeight = 1;
|
||||
if (iterationCount > 0)
|
||||
{
|
||||
refittingWeight = exp(-vcg::SquaredNorm(normal - previousGrad) * invSigma2);
|
||||
Scalar refittingWeight = 1;
|
||||
if (iterationCount > 0)
|
||||
{
|
||||
refittingWeight = exp(-vcg::SquaredNorm(normal - previousGrad) * invSigma2);
|
||||
// if (mSigmaR>0)
|
||||
// {
|
||||
// Scalar residual = (ddotn - potentialPrev) / mRadii.at(id);
|
||||
// refittingWeight *= exp(-residual*residual * invSigmaR2);
|
||||
// }
|
||||
}
|
||||
mCachedRefittingWeights.at(i) = refittingWeight;
|
||||
Scalar w = mCachedWeights.at(i) * refittingWeight;
|
||||
VectorType gw = mCachedWeightGradients.at(i) * refittingWeight;
|
||||
}
|
||||
mCachedRefittingWeights.at(i) = refittingWeight;
|
||||
Scalar w = mCachedWeights.at(i) * refittingWeight;
|
||||
VectorType gw = mCachedWeightGradients.at(i) * refittingWeight;
|
||||
|
||||
sumGradWeight += gw;
|
||||
sumGradWeightPotential += gw * f;
|
||||
sumN += normal * w;
|
||||
potential += w * f;
|
||||
sumW += w;
|
||||
}
|
||||
sumGradWeight += gw;
|
||||
sumGradWeightPotential += gw * f;
|
||||
sumN += normal * w;
|
||||
potential += w * f;
|
||||
sumW += w;
|
||||
}
|
||||
|
||||
if(sumW==0.)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
if(sumW==0.)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
potential /= sumW;
|
||||
grad = (-sumGradWeight*potential + sumGradWeightPotential + sumN) * (1./sumW);
|
||||
potential /= sumW;
|
||||
grad = (-sumGradWeight*potential + sumGradWeightPotential + sumN) * (1./sumW);
|
||||
|
||||
iterationCount++;
|
||||
iterationCount++;
|
||||
|
||||
} while ( (iterationCount < mMinRefittingIters)
|
||||
|| ( vcg::SquaredNorm(grad - previousGrad) > mRefittingThreshold && iterationCount < mMaxRefittingIters) );
|
||||
} while ( (iterationCount < mMinRefittingIters)
|
||||
|| ( vcg::SquaredNorm(grad - previousGrad) > mRefittingThreshold && iterationCount < mMaxRefittingIters) );
|
||||
|
||||
mCachedGradient = grad;
|
||||
mCachedPotential = potential;
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
mCachedGradient = grad;
|
||||
mCachedPotential = potential;
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
|
||||
mCachedSumGradWeight = sumGradWeight;
|
||||
mCachedSumN = sumN;
|
||||
mCachedSumW = sumW;
|
||||
mCachedSumGradPotential = sumGradWeightPotential;
|
||||
mCachedSumGradWeight = sumGradWeight;
|
||||
mCachedSumN = sumN;
|
||||
mCachedSumW = sumW;
|
||||
mCachedSumGradPotential = sumGradWeightPotential;
|
||||
|
||||
return true;
|
||||
return true;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool RIMLS<_MeshType>::mlsHessian(const VectorType& x, MatrixType& hessian) const
|
||||
{
|
||||
this->requestSecondDerivatives();
|
||||
// at this point we assume computePotentialAndGradient has been called first
|
||||
this->requestSecondDerivatives();
|
||||
// at this point we assume computePotentialAndGradient has been called first
|
||||
|
||||
uint nofSamples = mNeighborhood.size();
|
||||
uint nofSamples = mNeighborhood.size();
|
||||
|
||||
const VectorType& sumGradWeight = mCachedSumGradWeight;
|
||||
const VectorType& sumGradWeightPotential = mCachedSumGradPotential ;
|
||||
const VectorType& sumN = mCachedSumN;
|
||||
const Scalar& sumW = mCachedSumW;
|
||||
const Scalar invW = 1.f/sumW;
|
||||
const VectorType& sumGradWeight = mCachedSumGradWeight;
|
||||
const VectorType& sumGradWeightPotential = mCachedSumGradPotential ;
|
||||
const VectorType& sumN = mCachedSumN;
|
||||
const Scalar& sumW = mCachedSumW;
|
||||
const Scalar invW = 1.f/sumW;
|
||||
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
VectorType sumDGradWeight; sumDGradWeight.SetZero();
|
||||
VectorType sumDWeightNormal; sumDWeightNormal.SetZero();
|
||||
VectorType sumGradWeightNk; sumGradWeightNk.SetZero();
|
||||
VectorType sumDGradWeightPotential; sumDGradWeightPotential.SetZero();
|
||||
for (uint k=0 ; k<3 ; ++k)
|
||||
{
|
||||
VectorType sumDGradWeight; sumDGradWeight.SetZero();
|
||||
VectorType sumDWeightNormal; sumDWeightNormal.SetZero();
|
||||
VectorType sumGradWeightNk; sumGradWeightNk.SetZero();
|
||||
VectorType sumDGradWeightPotential; sumDGradWeightPotential.SetZero();
|
||||
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
VectorType p = mPoints[id].cP();
|
||||
VectorType diff = x - p;
|
||||
Scalar f = Dot(diff, mPoints[id].cN());
|
||||
for (unsigned int i=0; i<nofSamples; i++)
|
||||
{
|
||||
int id = mNeighborhood.index(i);
|
||||
VectorType p = mPoints[id].cP();
|
||||
VectorType diff = x - p;
|
||||
Scalar f = (diff * mPoints[id].cN());
|
||||
|
||||
VectorType gradW = mCachedWeightGradients.at(i) * mCachedRefittingWeights.at(i);
|
||||
VectorType dGradW = (x-p) * ( mCachedWeightSecondDerivatives.at(i) * (x[k]-p[k]) * mCachedRefittingWeights.at(i));
|
||||
dGradW[k] += mCachedWeightDerivatives.at(i);
|
||||
VectorType gradW = mCachedWeightGradients.at(i) * mCachedRefittingWeights.at(i);
|
||||
VectorType dGradW = (x-p) * ( mCachedWeightSecondDerivatives.at(i) * (x[k]-p[k]) * mCachedRefittingWeights.at(i));
|
||||
dGradW[k] += mCachedWeightDerivatives.at(i);
|
||||
|
||||
sumDGradWeight += dGradW;
|
||||
sumDWeightNormal += mPoints[id].cN() * gradW[k];
|
||||
sumGradWeightNk += gradW * mPoints[id].cN()[k];
|
||||
sumDGradWeightPotential += dGradW * f;
|
||||
}
|
||||
sumDGradWeight += dGradW;
|
||||
sumDWeightNormal += mPoints[id].cN() * gradW[k];
|
||||
sumGradWeightNk += gradW * mPoints[id].cN()[k];
|
||||
sumDGradWeightPotential += dGradW * f;
|
||||
}
|
||||
|
||||
VectorType dGrad = (
|
||||
sumDWeightNormal + sumGradWeightNk + sumDGradWeightPotential
|
||||
- sumDGradWeight * mCachedPotential
|
||||
- sumGradWeight * mCachedGradient[k]
|
||||
- mCachedGradient * sumGradWeight[k] ) * invW;
|
||||
VectorType dGrad = (
|
||||
sumDWeightNormal + sumGradWeightNk + sumDGradWeightPotential
|
||||
- sumDGradWeight * mCachedPotential
|
||||
- sumGradWeight * mCachedGradient[k]
|
||||
- mCachedGradient * sumGradWeight[k] ) * invW;
|
||||
|
||||
hessian.SetColumn(k,dGrad);
|
||||
}
|
||||
hessian.SetColumn(k,dGrad);
|
||||
}
|
||||
|
||||
return true;
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user