mirror of
https://github.com/lucaspalomodevelop/meshlab.git
synced 2026-03-20 03:16:10 +00:00
renamed mlsfilters to the standard naming scheme filter_mls
This commit is contained in:
Paolo Cignoni cignoni
parent
27e420778c
commit
83d559efe6
@ -1,11 +0,0 @@
|
||||
|
||||
#include <qglobal.h>
|
||||
#include "apss.h"
|
||||
#include "apss.tpp"
|
||||
#include <meshlab/meshmodel.h>
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template class APSS<CMeshO>;
|
||||
|
||||
}
|
||||
@ -1,120 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef APSS_H
|
||||
#define APSS_H
|
||||
|
||||
#include "mlssurface.h"
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _MeshType>
|
||||
class APSS : public MlsSurface<_MeshType>
|
||||
{
|
||||
typedef MlsSurface<_MeshType> Base;
|
||||
|
||||
typedef typename Base::Scalar Scalar;
|
||||
typedef typename Base::VectorType VectorType;
|
||||
typedef typename Base::MatrixType MatrixType;
|
||||
typedef _MeshType MeshType;
|
||||
using Base::mCachedQueryPointIsOK;
|
||||
using Base::mCachedQueryPoint;
|
||||
using Base::mNeighborhood;
|
||||
using Base::mCachedWeights;
|
||||
using Base::mCachedWeightDerivatives;
|
||||
using Base::mCachedWeightGradients;
|
||||
using Base::mCachedWeightSecondDerivatives;
|
||||
using Base::mBallTree;
|
||||
using Base::mPoints;
|
||||
using Base::mFilterScale;
|
||||
using Base::mMaxNofProjectionIterations;
|
||||
using Base::mAveragePointSpacing;
|
||||
using Base::mProjectionAccuracy;
|
||||
using Base::mGradientHint;
|
||||
|
||||
enum Status {ASS_SPHERE, ASS_PLANE, ASS_UNDETERMINED};
|
||||
|
||||
public:
|
||||
|
||||
APSS(const MeshType& m)
|
||||
: Base(m)
|
||||
{
|
||||
mSphericalParameter = 1;
|
||||
}
|
||||
|
||||
virtual Scalar potential(const VectorType& x, int* errorMask = 0) const;
|
||||
virtual VectorType gradient(const VectorType& x, int* errorMask = 0) const;
|
||||
virtual MatrixType hessian(const VectorType& x, int* errorMask) const;
|
||||
virtual VectorType project(const VectorType& x, VectorType* pNormal = 0, int* errorMask = 0) const;
|
||||
|
||||
/** \returns the approximation of the mean curvature obtained from the radius of the fitted sphere */
|
||||
virtual Scalar approxMeanCurvature(const VectorType& x, int* errorMask = 0) const;
|
||||
|
||||
void setSphericalParameter(Scalar v);
|
||||
|
||||
protected:
|
||||
bool fit(const VectorType& x) const;
|
||||
bool mlsGradient(const VectorType& x, VectorType& grad) const;
|
||||
bool mlsHessian(const VectorType& x, MatrixType& hessian) const;
|
||||
|
||||
protected:
|
||||
Scalar mSphericalParameter;
|
||||
|
||||
// use double precision anyway
|
||||
typedef double LScalar;
|
||||
typedef vcg::Point3<LScalar> LVector;
|
||||
|
||||
// cached algebraic sphere coefficients
|
||||
mutable LScalar uConstant;
|
||||
mutable LVector uLinear;
|
||||
mutable LScalar uQuad;
|
||||
|
||||
mutable LVector mCenter;
|
||||
mutable LScalar mRadius;
|
||||
mutable Status mStatus;
|
||||
|
||||
mutable LVector mCachedSumP;
|
||||
mutable LVector mCachedSumN;
|
||||
mutable LScalar mCachedSumDotPP;
|
||||
mutable LScalar mCachedSumDotPN;
|
||||
mutable LScalar mCachedSumW;
|
||||
|
||||
mutable LVector mCachedGradSumP[3];
|
||||
mutable LVector mCachedGradSumN[3];
|
||||
mutable LScalar mCachedGradSumDotPN[3];
|
||||
mutable LScalar mCachedGradSumDotPP[3];
|
||||
mutable LScalar mCachedGradSumW[3];
|
||||
|
||||
mutable LScalar mCachedGradNume[3];
|
||||
mutable LScalar mCachedGradDeno[3];
|
||||
|
||||
mutable LScalar mCachedGradUConstant[3];
|
||||
mutable LVector mCachedGradULinear[3];
|
||||
mutable LScalar mCachedGradUQuad[3];
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
//#include "apss.tpp"
|
||||
|
||||
#endif
|
||||
@ -1,499 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include "apss.h"
|
||||
#include "kdtree.h"
|
||||
#include "mlsutils.h"
|
||||
#include <iostream>
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _MeshType>
|
||||
void APSS<_MeshType>::setSphericalParameter(Scalar v)
|
||||
{
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
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 (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 ((!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());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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 * (Scalar(2) * uQuad);
|
||||
normal.Normalize();
|
||||
}
|
||||
|
||||
delta2 = vcg::SquaredNorm(previousPosition - position);
|
||||
} while ( delta2>epsilon2 && ++iterationCount<mMaxNofProjectionIterations);
|
||||
|
||||
if (pNormal)
|
||||
{
|
||||
if (mGradientHint==MLS_DERIVATIVE_ACCURATE)
|
||||
{
|
||||
VectorType grad;
|
||||
mlsGradient(vcg::Point3Cast<Scalar>(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;
|
||||
|
||||
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();
|
||||
|
||||
if (nofSamples==0)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
else if (nofSamples==1)
|
||||
{
|
||||
int id = mNeighborhood.index(0);
|
||||
LVector p = vcg::Point3Cast<LScalar>(mPoints[id].cP());
|
||||
LVector n = vcg::Point3Cast<LScalar>(mPoints[id].cN());
|
||||
|
||||
uLinear = n;
|
||||
uConstant = -vcg::Dot(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::Point3Cast<LScalar>(mPoints[id].cP());
|
||||
LVector n = vcg::Point3Cast<LScalar>(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;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
// 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);
|
||||
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;
|
||||
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
return true;
|
||||
}
|
||||
|
||||
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 LScalar nume = sumDotPN - invSumW * vcg::Dot(sumP, sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * vcg::Dot(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::Point3Cast<LScalar>(mPoints[id].cP());
|
||||
LVector n = vcg::Point3Cast<LScalar>(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);
|
||||
}
|
||||
|
||||
mCachedGradSumP[k] = dSumP;
|
||||
mCachedGradSumN[k] = dSumN;
|
||||
mCachedGradSumDotPN[k] = dSumDotPN;
|
||||
mCachedGradSumDotPP[k] = dSumDotPP;
|
||||
mCachedGradSumW[k] = dSumW;
|
||||
|
||||
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));
|
||||
|
||||
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);
|
||||
|
||||
grad[k] = dVecU0 + vcg::Dot(dVecU13,vcg::Point3Cast<LScalar>(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;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool APSS<_MeshType>::mlsHessian(const VectorType& x, MatrixType& hessian) const
|
||||
{
|
||||
this->requestSecondDerivatives();
|
||||
|
||||
// TODO call mlsGradient first
|
||||
VectorType grad;
|
||||
mlsGradient(x,grad);
|
||||
|
||||
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 LScalar nume = sumDotPN - invSumW * vcg::Dot(sumP, sumN);
|
||||
const LScalar deno = sumDotPP - invSumW * vcg::Dot(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];
|
||||
|
||||
LScalar dVecU0 = mCachedGradUConstant[k];
|
||||
LVector dVecU13 = mCachedGradULinear[k];
|
||||
LScalar dVecU4 = mCachedGradUQuad[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::Point3Cast<LScalar>(mPoints[id].cP());
|
||||
LVector n = vcg::Point3Cast<LScalar>(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);
|
||||
|
||||
d2SumW += d2w;
|
||||
d2SumP += p*d2w;
|
||||
d2SumN += n*d2w;
|
||||
d2SumDotPN += d2w * vcg::Dot(n,p);
|
||||
d2SumDotPP += d2w * vcg::Dot(p,p);
|
||||
}
|
||||
|
||||
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 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 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;
|
||||
|
||||
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;
|
||||
|
||||
hessian[j][k] =
|
||||
dVecU13[j] + 2.*dVecU4*x[j]
|
||||
+ d2u0 + vcg::Dot(d2u13,vcg::Point3Cast<LScalar>(x)) + d2u4*vcg::Dot(x,x)
|
||||
+ mCachedGradULinear[j][k] + (j==k ? 2.*uQuad : 0.) + 2.*x[k]*mCachedGradUQuad[j];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,159 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include "balltree.h"
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _Scalar>
|
||||
BallTree<_Scalar>::BallTree(const ConstDataWrapper<VectorType>& points, const ConstDataWrapper<Scalar>& radii)
|
||||
: mPoints(points), mRadii(radii), mRadiusScale(1.), mTreeIsUptodate(false)
|
||||
{
|
||||
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();
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _Scalar>
|
||||
void BallTree<_Scalar>::rebuild(void)
|
||||
{
|
||||
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.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);
|
||||
|
||||
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::Distance(mPoints[i], aabbLeft) < mRadii[i]*mRadiusScale)
|
||||
iLeft.push_back(i);
|
||||
|
||||
if (vcg::Distance(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 = 0.5*(aabb.max[dim] + aabb.min[dim]);
|
||||
node.leaf = 0;
|
||||
|
||||
AxisAlignedBoxType aabbLeft=aabb, aabbRight=aabb;
|
||||
aabbLeft.max[dim] = node.splitValue;
|
||||
aabbRight.min[dim] = node.splitValue;
|
||||
|
||||
std::vector<int> iLeft, iRight;
|
||||
split(indices, aabbLeft, aabbRight, iLeft,iRight);
|
||||
|
||||
// we don't need the index list anymore
|
||||
indices.clear();
|
||||
|
||||
{
|
||||
// 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>;
|
||||
template class BallTree<double>;
|
||||
|
||||
}
|
||||
@ -1,119 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef BALLTREE_H
|
||||
#define BALLTREE_H
|
||||
|
||||
#include <vcg/space/point3.h>
|
||||
#include <vcg/space/box3.h>
|
||||
#include "mlsutils.h"
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _Scalar>
|
||||
class Neighborhood
|
||||
{
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
|
||||
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 insert(int id, Scalar d2) { mIndices.push_back(id); mSqDists.push_back(d2); }
|
||||
|
||||
protected:
|
||||
std::vector<int> mIndices;
|
||||
std::vector<Scalar> mSqDists;
|
||||
};
|
||||
|
||||
template<typename _Scalar>
|
||||
class BallTree
|
||||
{
|
||||
public:
|
||||
typedef _Scalar Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
|
||||
BallTree(const ConstDataWrapper<VectorType>& points, const ConstDataWrapper<Scalar>& radii);
|
||||
|
||||
void computeNeighbors(const VectorType& x, Neighborhood<Scalar>* pNei) const;
|
||||
|
||||
void setRadiusScale(Scalar v) { mRadiusScale = v; mTreeIsUptodate = false; }
|
||||
|
||||
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;
|
||||
};
|
||||
};
|
||||
};
|
||||
|
||||
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;
|
||||
|
||||
protected:
|
||||
ConstDataWrapper<VectorType> mPoints;
|
||||
ConstDataWrapper<Scalar> mRadii;
|
||||
Scalar mRadiusScale;
|
||||
|
||||
int mMaxTreeDepth;
|
||||
int mTargetCellSize;
|
||||
mutable bool mTreeIsUptodate;
|
||||
mutable VectorType mQueryPosition;
|
||||
|
||||
Node* mRootNode;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -1,187 +0,0 @@
|
||||
/****************************************************************************
|
||||
* VCGLib o o *
|
||||
* Visual and Computer Graphics Library o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2004 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef __IMPLICITS
|
||||
#define __IMPLICITS
|
||||
|
||||
#include <vcg/space/point2.h>
|
||||
#include <vcg/space/point3.h>
|
||||
#include <vcg/math/lin_algebra.h>
|
||||
|
||||
namespace vcg
|
||||
{
|
||||
namespace implicits
|
||||
{
|
||||
|
||||
/** \returns the Gauss curvature directly from the gradient and hessian */
|
||||
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];
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
/** \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);
|
||||
}
|
||||
|
||||
/** This class computes the Weingarten map of a scalar field and provides
|
||||
* methods to extract curvatures from it.
|
||||
*
|
||||
* The Weingarten map is equal to 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;
|
||||
|
||||
/** 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();
|
||||
m_normal = grad * invL;
|
||||
|
||||
Matrix33<Scalar> I; I.SetIdentity();
|
||||
m_nnT.ExternalProduct(m_normal,m_normal);
|
||||
|
||||
m_w = (I-m_nnT) * hess * invL;
|
||||
|
||||
m_kgIsDirty = true;
|
||||
m_kmIsDirty = true;
|
||||
m_kpAreDirty = true;
|
||||
m_kdirAreDirty = true;
|
||||
}
|
||||
|
||||
/** \returns the Weingarten map matrix */
|
||||
const MatrixType& W() const { return m_w; }
|
||||
|
||||
/** \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 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 first principal curvature */
|
||||
Scalar K1() const { updateKp(); return m_k1; }
|
||||
|
||||
/** \returns the second principal curvature */
|
||||
Scalar K2() const { updateKp(); return m_k2; }
|
||||
|
||||
/** \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;
|
||||
}
|
||||
}
|
||||
|
||||
inline void extractEigenvectors() const
|
||||
{
|
||||
if (m_kdirAreDirty)
|
||||
{
|
||||
MatrixType copy = m_w;
|
||||
int mrot = 0;
|
||||
VectorType evals;
|
||||
MatrixType evecs;
|
||||
Jacobi(copy, evals, evecs, mrot);
|
||||
VectorType evalsAbs(fabs(evals[0]),fabs(evals[0]),fabs(evals[0]));
|
||||
SortEigenvaluesAndEigenvectors(evals,evecs,true);
|
||||
m_k1 = evals[0];
|
||||
m_k2 = evals[1];
|
||||
m_k1dir = evecs[0];
|
||||
m_k2dir = evecs[1];
|
||||
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;
|
||||
};
|
||||
|
||||
} // namespace implicits
|
||||
} // namespace vcg
|
||||
|
||||
#endif //__IMPLICITS
|
||||
@ -1,211 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include "kdtree.h"
|
||||
#include <limits>
|
||||
#include <iostream>
|
||||
|
||||
template<typename Scalar>
|
||||
KdTree<Scalar>::KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell, unsigned int maxDepth)
|
||||
: mPoints(points.size())
|
||||
{
|
||||
// compute the AABB of the input
|
||||
mPoints[0] = points[0];
|
||||
mAABB.Set(mPoints[0]);
|
||||
for (unsigned int i=1 ; i<mPoints.size() ; ++i)
|
||||
{
|
||||
mPoints[i] = points[i];
|
||||
mAABB.Add(mPoints[i]);
|
||||
}
|
||||
|
||||
mNodes.reserve(4*mPoints.size()/nofPointsPerCell);
|
||||
mNodes.resize(1);
|
||||
mNodes.back().leaf = 0;
|
||||
createTree(0, 0, mPoints.size(), 1, nofPointsPerCell, maxDepth);
|
||||
}
|
||||
|
||||
template<typename Scalar>
|
||||
KdTree<Scalar>::~KdTree()
|
||||
{
|
||||
}
|
||||
|
||||
template<typename Scalar>
|
||||
void KdTree<Scalar>::setMaxNofNeighbors(unsigned int k)
|
||||
{
|
||||
mNeighborQueue.setMaxSize(k);
|
||||
}
|
||||
|
||||
/** Performs the kNN query.
|
||||
*
|
||||
* This algorithm uses the simple distance to the split plane to prune nodes.
|
||||
* A more elaborated approach consists to track the closest corner of the cell
|
||||
* relatively to the current query point. This strategy allows to save about 5%
|
||||
* of the leaves. However, in practice the slight overhead due to this tracking
|
||||
* reduces the overall performance.
|
||||
*
|
||||
* This algorithm also use a simple stack while a priority queue using the squared
|
||||
* distances to the cells as a priority values allows to save about 10% of the leaves.
|
||||
* But, again, priority queue insertions and deletions are quite involved, and therefore
|
||||
* a simple stack is by far much faster.
|
||||
*/
|
||||
template<typename Scalar>
|
||||
void KdTree<Scalar>::doQueryK(const VectorType& queryPoint)
|
||||
{
|
||||
mNeighborQueue.init();
|
||||
mNeighborQueue.insert(0xffffffff, std::numeric_limits<float>::max());
|
||||
|
||||
mNodeStack[0].nodeId = 0;
|
||||
mNodeStack[0].sq = 0.f;
|
||||
unsigned int count = 1;
|
||||
|
||||
while (count)
|
||||
{
|
||||
QueryNode& qnode = mNodeStack[count-1];
|
||||
Node& node = mNodes[qnode.nodeId];
|
||||
|
||||
if (qnode.sq < mNeighborQueue.getTopWeight())
|
||||
{
|
||||
if (node.leaf)
|
||||
{
|
||||
--count; // pop
|
||||
unsigned int end = node.start+node.size;
|
||||
for (unsigned int i=node.start ; i<end ; ++i)
|
||||
mNeighborQueue.insert(i, vcg::SquaredNorm(queryPoint - mPoints[i]));
|
||||
}
|
||||
else
|
||||
{
|
||||
// replace the stack top by the farthest and push the closest
|
||||
float new_off = queryPoint[node.dim] - node.splitValue;
|
||||
if (new_off < 0.)
|
||||
{
|
||||
mNodeStack[count].nodeId = node.firstChildId;
|
||||
qnode.nodeId = node.firstChildId+1;
|
||||
}
|
||||
else
|
||||
{
|
||||
mNodeStack[count].nodeId = node.firstChildId+1;
|
||||
qnode.nodeId = node.firstChildId;
|
||||
}
|
||||
mNodeStack[count].sq = qnode.sq;
|
||||
qnode.sq = new_off*new_off;
|
||||
++count;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// pop
|
||||
--count;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Scalar>
|
||||
unsigned int KdTree<Scalar>::split(int start, int end, unsigned int dim, float splitValue)
|
||||
{
|
||||
int l(start), r(end-1);
|
||||
for ( ; l<r ; ++l, --r)
|
||||
{
|
||||
while (l < end && mPoints[l][dim] < splitValue)
|
||||
l++;
|
||||
while (r >= start && mPoints[r][dim] >= splitValue)
|
||||
r--;
|
||||
if (l > r)
|
||||
break;
|
||||
std::swap(mPoints[l],mPoints[r]);
|
||||
}
|
||||
return (mPoints[l][dim] < splitValue ? l+1 : l);
|
||||
}
|
||||
|
||||
/** recursively builds the kdtree
|
||||
*
|
||||
* The heuristic is the following:
|
||||
* - if the number of points in the node is lower than targetCellsize then make a leaf
|
||||
* - else compute the AABB of the points of the node and split it at the middle of
|
||||
* the largest AABB dimension.
|
||||
*
|
||||
* This strategy might look not optimal because it does not explicitly prune empty space,
|
||||
* unlike more advanced SAH-like techniques used for RT. On the other hand it leads to a shorter tree,
|
||||
* faster to traverse and our experience shown that in the special case of kNN queries,
|
||||
* this strategy is indeed more efficient (and much faster to build). Moreover, for volume data
|
||||
* (e.g., fluid simulation) pruning the empty space is useless.
|
||||
*
|
||||
* Actually, storing at each node the exact AABB (we therefore have a binary BVH) allows
|
||||
* to prune only about 10% of the leaves, but the overhead of this pruning (ball/ABBB intersection)
|
||||
* is more expensive than the gain it provides and the memory consumption is x4 higher !
|
||||
*/
|
||||
template<typename Scalar>
|
||||
void KdTree<Scalar>::createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellSize, unsigned int targetMaxDepth)
|
||||
{
|
||||
Node& node = mNodes[nodeId];
|
||||
AxisAlignedBoxType aabb;
|
||||
aabb.Set(mPoints[start]);
|
||||
for (unsigned int i=start+1 ; i<end ; ++i)
|
||||
aabb.Add(mPoints[i]);
|
||||
|
||||
VectorType diag = aabb.max - aabb.min;
|
||||
unsigned int dim = vcg::MaxCoeffId(diag);
|
||||
node.dim = dim;
|
||||
node.splitValue = 0.5*(aabb.max[dim] + aabb.min[dim]);
|
||||
|
||||
unsigned int midId = split(start, end, dim, node.splitValue);
|
||||
|
||||
node.firstChildId = mNodes.size();
|
||||
mNodes.resize(mNodes.size()+2);
|
||||
|
||||
{
|
||||
// left child
|
||||
unsigned int childId = mNodes[nodeId].firstChildId;
|
||||
Node& child = mNodes[childId];
|
||||
if (midId-start <= targetCellSize || level>=targetMaxDepth)
|
||||
{
|
||||
child.leaf = 1;
|
||||
child.start = start;
|
||||
child.size = midId-start;
|
||||
}
|
||||
else
|
||||
{
|
||||
child.leaf = 0;
|
||||
createTree(childId, start, midId, level+1, targetCellSize, targetMaxDepth);
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
// right child
|
||||
unsigned int childId = mNodes[nodeId].firstChildId+1;
|
||||
Node& child = mNodes[childId];
|
||||
if (end-midId <= targetCellSize || level>=targetMaxDepth)
|
||||
{
|
||||
child.leaf = 1;
|
||||
child.start = midId;
|
||||
child.size = end-midId;
|
||||
}
|
||||
else
|
||||
{
|
||||
child.leaf = 0;
|
||||
createTree(childId, midId, end, level+1, targetCellSize, targetMaxDepth);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template class KdTree<float>;
|
||||
template class KdTree<double>;
|
||||
@ -1,103 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef KDTREE_H
|
||||
#define KDTREE_H
|
||||
|
||||
#include <vcg/space/point3.h>
|
||||
#include <vcg/space/box3.h>
|
||||
#include "mlsutils.h"
|
||||
#include "priorityqueue.h"
|
||||
#include <vector>
|
||||
|
||||
template<typename _Scalar>
|
||||
class KdTree
|
||||
{
|
||||
public:
|
||||
|
||||
typedef _Scalar Scalar;
|
||||
typedef vcg::Point3<Scalar> VectorType;
|
||||
typedef vcg::Box3<Scalar> AxisAlignedBoxType;
|
||||
|
||||
struct Node
|
||||
{
|
||||
union {
|
||||
struct {
|
||||
Scalar splitValue;
|
||||
unsigned int firstChildId:24;
|
||||
unsigned int dim:2;
|
||||
unsigned int leaf:1;
|
||||
};
|
||||
struct {
|
||||
unsigned int start;
|
||||
unsigned short size;
|
||||
};
|
||||
};
|
||||
};
|
||||
typedef std::vector<Node> NodeList;
|
||||
|
||||
inline const NodeList& _getNodes(void) { return mNodes; }
|
||||
inline const std::vector<VectorType>& _getPoints(void) { return mPoints; }
|
||||
|
||||
void setMaxNofNeighbors(unsigned int k);
|
||||
inline int getNofFoundNeighbors(void) { return mNeighborQueue.getNofElements(); }
|
||||
inline const VectorType& getNeighbor(int i) { return mPoints[ mNeighborQueue.getIndex(i) ]; }
|
||||
inline unsigned int getNeighborId(int i) { return mNeighborQueue.getIndex(i); }
|
||||
inline float getNeighborSquaredDistance(int i) { return mNeighborQueue.getWeight(i); }
|
||||
|
||||
public:
|
||||
|
||||
KdTree(const ConstDataWrapper<VectorType>& points, unsigned int nofPointsPerCell = 16, unsigned int maxDepth = 64);
|
||||
|
||||
~KdTree();
|
||||
|
||||
void doQueryK(const VectorType& p);
|
||||
|
||||
protected:
|
||||
|
||||
// element of the stack
|
||||
struct QueryNode
|
||||
{
|
||||
QueryNode() {}
|
||||
QueryNode(unsigned int id) : nodeId(id) {}
|
||||
unsigned int nodeId; // id of the next node
|
||||
Scalar sq; // squared distance to the next node
|
||||
};
|
||||
|
||||
// used to build the tree: split the subset [start..end[ according to dim and splitValue,
|
||||
// and returns the index of the first element of the second subset
|
||||
unsigned int split(int start, int end, unsigned int dim, float splitValue);
|
||||
|
||||
void createTree(unsigned int nodeId, unsigned int start, unsigned int end, unsigned int level, unsigned int targetCellsize, unsigned int targetMaxDepth);
|
||||
|
||||
protected:
|
||||
|
||||
AxisAlignedBoxType mAABB;
|
||||
NodeList mNodes;
|
||||
std::vector<VectorType> mPoints;
|
||||
HeapMaxPriorityQueue<int,Scalar> mNeighborQueue;
|
||||
QueryNode mNodeStack[64];
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -1,18 +0,0 @@
|
||||
include(../../shared.pri)
|
||||
|
||||
# CONFIG += debug
|
||||
|
||||
HEADERS = rimls.h mlsplugin.h \
|
||||
../../meshlab/meshmodel.h
|
||||
|
||||
SOURCES = balltree.cpp kdtree.cpp mlsplugin.cpp \
|
||||
apss.cpp rimls.cpp \
|
||||
../../meshlab/filterparameter.cpp \
|
||||
../../meshlab/meshmodel.cpp \
|
||||
$$GLEWCODE
|
||||
|
||||
TARGET = mlsfilters
|
||||
|
||||
CONFIG += opengl
|
||||
CONFIG += warn_off
|
||||
|
||||
@ -1,280 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef MLSMARCHINGCUBES_H
|
||||
#define MLSMARCHINGCUBES_H
|
||||
|
||||
#include <vcg/space/point3.h>
|
||||
#include <vcg/space/box3.h>
|
||||
#include <map>
|
||||
#include "mlssurface.h"
|
||||
#include "mlsutils.h"
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
|
||||
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;
|
||||
|
||||
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());
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
int resolution;
|
||||
|
||||
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();
|
||||
|
||||
// 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.1;
|
||||
mAABB.max += diag * 0.1;
|
||||
diag = mAABB.max - mAABB.min;
|
||||
|
||||
if ( (diag[0]<=0.)
|
||||
|| (diag[1]<=0.)
|
||||
|| (diag[2]<=0.)
|
||||
|| (resolution==0))
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
mCache = new GridElement[(mMaxBlockSize)*(mMaxBlockSize)*(mMaxBlockSize)];
|
||||
ScalarType step = vcg::MaxCoeff(diag)/ScalarType(resolution);
|
||||
|
||||
unsigned int nofCells[3];
|
||||
unsigned int nofBlocks[3];
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
_mesh = &mesh;
|
||||
_mesh->Clear();
|
||||
|
||||
int countSubSlice = 0;
|
||||
int totalSubSlices = nofBlocks[2] * nofBlocks[1] * nofCells[0];
|
||||
|
||||
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);
|
||||
|
||||
// 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));
|
||||
|
||||
// fill the grid
|
||||
vcg::Point3i ci; // local cell id
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
|
||||
if (!out)
|
||||
{
|
||||
extractor.ProcessCell(ci+mBlockOrigin, ci+mBlockOrigin+vcg::Point3i(1,1,1));
|
||||
}
|
||||
}
|
||||
}
|
||||
extractor.Finalize();
|
||||
_mesh = NULL;
|
||||
delete[] mCache;
|
||||
};
|
||||
|
||||
int GetLocalCellId(const vcg::Point3i& p)
|
||||
{
|
||||
return p[0] + (p[1] + p[2]*mMaxBlockSize)*mMaxBlockSize;
|
||||
}
|
||||
|
||||
int GetLocalCellIdFromGlobal(vcg::Point3i p)
|
||||
{
|
||||
p -= mBlockOrigin;
|
||||
return GetLocalCellId(p);
|
||||
}
|
||||
|
||||
float V(int pi, int pj, int pk)
|
||||
{
|
||||
return mCache[GetLocalCellIdFromGlobal(vcg::Point3i(pi, pj, pk))].value;
|
||||
}
|
||||
|
||||
int GetGlobalCellId(const vcg::Point3i &p)
|
||||
{
|
||||
return p[0] + p[1] * resolution + p[2] * resolution * resolution;
|
||||
}
|
||||
|
||||
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 = 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;
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
};
|
||||
} // end namespace
|
||||
} // end namespace
|
||||
|
||||
#endif
|
||||
|
||||
@ -1,605 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include <QtGui>
|
||||
|
||||
#define NDEBUG
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <time.h>
|
||||
#include <iostream>
|
||||
|
||||
#include <meshlab/meshmodel.h>
|
||||
#include <meshlab/interfaces.h>
|
||||
|
||||
#include <vcg/complex/trimesh/clean.h>
|
||||
#include <vcg/complex/trimesh/refine.h>
|
||||
#include <vcg/complex/trimesh/update/selection.h>
|
||||
#include <vcg/complex/trimesh/update/normal.h>
|
||||
#include <vcg/complex/trimesh/update/bounding.h>
|
||||
#include <vcg/complex/trimesh/append.h>
|
||||
#include <vcg/complex/trimesh/create/advancing_front.h>
|
||||
#include <vcg/complex/trimesh/create/marching_cubes.h>
|
||||
#include <vcg/complex/trimesh/update/quality.h>
|
||||
#include "mlsmarchingcube.h"
|
||||
|
||||
#include "mlsplugin.h"
|
||||
#include "rimls.h"
|
||||
#include "apss.h"
|
||||
#include "mlsutils.h"
|
||||
#include "implicits.h"
|
||||
#include "../meshfilter/refine_loop.h"
|
||||
|
||||
#include "smallcomponentselection.h"
|
||||
|
||||
using namespace GaelMls;
|
||||
using namespace vcg;
|
||||
|
||||
// Constructor usually performs only two simple tasks of filling the two lists
|
||||
// - typeList: with all the possible id of the filtering actions
|
||||
// - 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
|
||||
};
|
||||
|
||||
MlsPlugin::MlsPlugin()
|
||||
{
|
||||
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;
|
||||
|
||||
// initFilterList(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)
|
||||
const QString MlsPlugin::filterName(FilterIDType filterId)
|
||||
{
|
||||
switch(filterId) {
|
||||
case FP_APSS_PROJECTION : return QString("MLS projection (APSS)");
|
||||
case FP_RIMLS_PROJECTION : return QString("MLS projection (#####)");
|
||||
case FP_APSS_AFRONT : return QString("MLS meshing/APSS Advancing Front");
|
||||
case FP_RIMLS_AFRONT : return QString("MLS meshing/##### Advancing Front");
|
||||
case FP_APSS_MCUBE : return QString("Marching Cubes (APSS)");
|
||||
case FP_RIMLS_MCUBE : return QString("Marching Cubes (#####)");
|
||||
case FP_APSS_COLORIZE : return QString("Colorize curvature (APSS)");
|
||||
case FP_RIMLS_COLORIZE : return QString("Colorize curvature (#####)");
|
||||
case FP_RADIUS_FROM_DENSITY : return QString("Estimate radius from density");
|
||||
case FP_SELECT_SMALL_COMPONENTS : return QString("Small component selection");
|
||||
default : assert(0);
|
||||
}
|
||||
}
|
||||
|
||||
const MeshFilterInterface::FilterClass MlsPlugin::getClass(QAction *a)
|
||||
{
|
||||
int id = ID(a);
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
return MeshFilterInterface::Selection;
|
||||
}
|
||||
|
||||
return MeshFilterInterface::PointSet;
|
||||
}
|
||||
|
||||
// Info() must return the longer string describing each filtering action
|
||||
// (this string is used in the About plugin dialog)
|
||||
const QString MlsPlugin::filterInfo(FilterIDType filterId)
|
||||
{
|
||||
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 & _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>"
|
||||
"See [Guennebaud and Gross, Algebraic Point Set Surfaces, Siggraph 2007]"
|
||||
"and [Guennebaud et al., Dynamic Sampling and Rendering of APSS, Eurographics 2008]"
|
||||
"for all the details about APSS.";
|
||||
}
|
||||
|
||||
if (filterId & _RIMLS_)
|
||||
{
|
||||
str +=
|
||||
"<br>This is the ##### MLS variant.";
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
// This function define the needed parameters for each filter. Return true if the filter has some parameters
|
||||
// it is called every time, so you can set the default value of parameters according to the mesh
|
||||
// For each parmeter you need to define,
|
||||
// - the name of the parameter,
|
||||
// - the string shown in the dialog
|
||||
// - the default value
|
||||
// - 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, FilterParameterSet& parlst)
|
||||
//void ExtraSamplePlugin::initParList(QAction *action, MeshModel &m, FilterParameterSet &parlst)
|
||||
{
|
||||
int id = ID(action);
|
||||
MeshModel *target = md.mm();
|
||||
|
||||
if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
parlst.addFloat("NbFaceRatio",
|
||||
0.1,
|
||||
"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.addBool( "NonClosedOnly",
|
||||
false,
|
||||
"Select only non closed components",
|
||||
"");
|
||||
return;
|
||||
}
|
||||
else if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
parlst.addInt("NbNeighbors",
|
||||
16,
|
||||
"Number of neighbors",
|
||||
"Number of neighbors used to estimate the local density. Larger values lead to smoother variations.");
|
||||
return;
|
||||
}
|
||||
|
||||
if ((id & _PROJECTION_) || (id & _COLORIZE_))
|
||||
{
|
||||
parlst.addMesh( "ControlMesh", target, "Point set",
|
||||
"The point set (or mesh) which defines the MLS surface.");
|
||||
parlst.addMesh( "ProxyMesh", target, "Proxy Mesh",
|
||||
"The mesh that will be projected/resampled onto the MLS surface.");
|
||||
parlst.addBool( "SelectionOnly",
|
||||
target->cm.sfn>0,
|
||||
"Selection only",
|
||||
"If checked, only selected vertices will be projected.");
|
||||
}
|
||||
|
||||
if ( (id & _APSS_) || (id & _RIMLS_) )
|
||||
{
|
||||
parlst.addFloat("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.addFloat("ProjectionAccuracy",
|
||||
1e-4,
|
||||
"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.addInt( "MaxProjectionIters",
|
||||
15,
|
||||
"Projection - Max iterations (adv)",
|
||||
"Max number of iterations for the projection.");
|
||||
}
|
||||
|
||||
if (id & _APSS_)
|
||||
{
|
||||
parlst.addFloat("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.addBool( "AccurateNormal",
|
||||
true,
|
||||
"Accurate normals",
|
||||
"If checked, use the accurate MLS gradient instead of the local approximation"
|
||||
"to compute the normals.");
|
||||
}
|
||||
|
||||
if (id & _RIMLS_)
|
||||
{
|
||||
parlst.addFloat("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.addInt( "MaxRefittingIters",
|
||||
3,
|
||||
"MLS - Max fitting iterations",
|
||||
"Max number of fitting iterations. (0 or 1 is equivalent to the standard IMLS)");
|
||||
}
|
||||
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
parlst.addInt( "MaxSubdivisions",
|
||||
0,
|
||||
"Refinement - Max subdivisions",
|
||||
"Max number of subdivisions.");
|
||||
parlst.addFloat("ThAngleInDegree",
|
||||
2,
|
||||
"Refinement - Crease angle (degree)",
|
||||
"Threshold angle between two faces controlling the refinement.");
|
||||
}
|
||||
|
||||
if (id & _AFRONT_)
|
||||
{
|
||||
}
|
||||
|
||||
if ((id & _COLORIZE_))
|
||||
{
|
||||
QStringList lst;
|
||||
lst << "Mean" << "Gauss" << "K1" << "K2";
|
||||
if (id & _APSS_)
|
||||
lst << "ApproxMean";
|
||||
|
||||
parlst.addEnum("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.addBool( "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.addInt( "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.");
|
||||
}
|
||||
}
|
||||
|
||||
const int MlsPlugin::getRequirements(QAction *action)
|
||||
{
|
||||
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:
|
||||
// vcg::face::Pos<typename MESH_TYPE::FaceType> op = ep;
|
||||
// op.FlipF();
|
||||
// if (op.f)
|
||||
// return vcg::Dot(ep.f->cN(), op.f->cN()) < thCosAngle;
|
||||
// else
|
||||
// return true;
|
||||
|
||||
return vcg::Dot(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;
|
||||
|
||||
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, FilterParameterSet& par, vcg::CallBackPos* cb)
|
||||
{
|
||||
int id = ID(filter);
|
||||
|
||||
if (id == FP_RADIUS_FROM_DENSITY)
|
||||
{
|
||||
CMeshO::VertContainer& points = md.mm()->cm.vert;
|
||||
if (!points.RadiusEnabled)
|
||||
{
|
||||
points.EnableRadius();
|
||||
}
|
||||
int nbNeighbors = par.getInt("NbNeighbors");
|
||||
|
||||
assert(points.size()>=2);
|
||||
KdTree<float> knn(ConstDataWrapper<vcg::Point3f>(&points[0].cP(), points.size(),
|
||||
size_t(points[1].cP().V()) - size_t(points[0].cP().V())));
|
||||
|
||||
knn.setMaxNofNeighbors(nbNeighbors);
|
||||
for (size_t i = 0; i< points.size(); i++)
|
||||
{
|
||||
knn.doQueryK(points[i].cP());
|
||||
points[i].R() = 2. * sqrt(knn.getNeighborSquaredDistance(0)/float(knn.getNofFoundNeighbors()));
|
||||
}
|
||||
}
|
||||
else if (id == FP_SELECT_SMALL_COMPONENTS)
|
||||
{
|
||||
MeshModel* mesh = md.mm();
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO | MeshModel::MM_FACEFLAGBORDER);
|
||||
bool nonClosedOnly = par.getBool("NonClosedOnly");
|
||||
float ratio = par.getFloat("NbFaceRatio");
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, ratio, nonClosedOnly);
|
||||
}
|
||||
else
|
||||
{
|
||||
// we are doing some MLS based stuff
|
||||
|
||||
MeshModel* pPoints = 0;
|
||||
if (id & _PROJECTION_)
|
||||
{
|
||||
if (par.getMesh("ControlMesh") == par.getMesh("ProxyMesh"))
|
||||
{
|
||||
// clone the control mesh
|
||||
MeshModel* ref = par.getMesh("ControlMesh");
|
||||
pPoints = new MeshModel();
|
||||
pPoints->updateDataMask(ref->currentDataMask);
|
||||
vcg::tri::Append<CMeshO,CMeshO>::Mesh(pPoints->cm, ref->cm, false);
|
||||
vcg::tri::UpdateBounding<CMeshO>::Box(pPoints->cm);
|
||||
pPoints->cm.Tr = ref->cm.Tr;
|
||||
}
|
||||
else
|
||||
pPoints = par.getMesh("ControlMesh");
|
||||
}
|
||||
else
|
||||
pPoints = md.mm();
|
||||
|
||||
// create the MLS surface
|
||||
cb(1, "Create the MLS data structures...");
|
||||
MlsSurface<CMeshO>* mls = 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);
|
||||
}
|
||||
|
||||
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 (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;
|
||||
|
||||
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.);
|
||||
|
||||
int nbRefinements = par.getInt("MaxSubdivisions");
|
||||
for (int k=0; k<nbRefinements+1; ++k)
|
||||
{
|
||||
//UpdateFaceNormalFromVertex(m.cm);
|
||||
if (k!=0)
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO | MeshModel::MM_FACEFLAGBORDER);
|
||||
|
||||
vcg::tri::UpdateNormals<CMeshO>::PerFace(mesh->cm);
|
||||
vcg::tri::UpdateNormals<CMeshO>::NormalizeFace(mesh->cm);
|
||||
//vcg::RefineE<CMeshO,vcg::MidPoint<CMeshO> >(m.cm, vcg::MidPoint<CMeshO>(), edgePred, false, cb);
|
||||
vcg::RefineOddEvenE<CMeshO, vcg::OddPointLoop<CMeshO>, vcg::EvenPointLoop<CMeshO> >
|
||||
(mesh->cm, vcg::OddPointLoop<CMeshO>(), vcg::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());
|
||||
}
|
||||
}
|
||||
|
||||
Log(0, "Successfully projected %i vertices", mesh->cm.vn);
|
||||
}
|
||||
else if (id & _COLORIZE_)
|
||||
{
|
||||
mesh = par.getMesh("ProxyMesh");
|
||||
|
||||
bool selectionOnly = par.getBool("SelectionOnly");
|
||||
//bool approx = apss && par.getBool("ApproxCurvature");
|
||||
int ct = par.getEnum("CurvatureType");
|
||||
|
||||
int size = mesh->cm.vert.size();
|
||||
//std::vector<float> curvatures(size);
|
||||
float minc=1e9, maxc=-1e9, minabsc=1e9;
|
||||
vcg::Point3f grad;
|
||||
vcg::Matrix33f 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...");
|
||||
|
||||
if ( (!selectionOnly) || (pPoints->cm.vert[i].IsS()) )
|
||||
{
|
||||
Point3f p = mls->project(mesh->cm.vert[i].P());
|
||||
float c = 0;
|
||||
|
||||
if (ct==CT_APSS)
|
||||
c = apss->approxMeanCurvature(p);
|
||||
else
|
||||
{
|
||||
grad = mls->gradient(p);
|
||||
hess = mls->hessian(p);
|
||||
implicits::WeingartenMap<float> W(grad,hess);
|
||||
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");
|
||||
}
|
||||
}
|
||||
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::tri::Stat<CMeshO>::ComputePerVertexQualityHistogram(mesh->cm,H);
|
||||
vcg::tri::UpdateColor<CMeshO>::VertexQualityRamp(mesh->cm,H.Percentile(0.01f),H.Percentile(0.99f));
|
||||
// for (unsigned int i = 0; i< size; i++)
|
||||
// {
|
||||
// if ( (!selectionOnly) || (mesh->cm.vert[i].IsS()) )
|
||||
// {
|
||||
// // float c = curvatures[i];
|
||||
// mesh->cm.vert[i].C().ColorRamp(minc,maxc, mesh->cm.vert[i].Q());
|
||||
// }
|
||||
// }
|
||||
}
|
||||
// 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");
|
||||
|
||||
// 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());
|
||||
}
|
||||
|
||||
// extra zero detection and removal
|
||||
{
|
||||
mesh->updateDataMask(MeshModel::MM_FACEFACETOPO | MeshModel::MM_FACEFLAGBORDER);
|
||||
// selection...
|
||||
vcg::tri::SmallComponent<CMeshO>::Select(mesh->cm, 0.1);
|
||||
// deletion...
|
||||
vcg::tri::SmallComponent<CMeshO>::DeleteFaceVert(mesh->cm);
|
||||
mesh->clearDataMask(MeshModel::MM_FACEFACETOPO | MeshModel::MM_FACEFLAGBORDER);
|
||||
}
|
||||
|
||||
Log(0, "Marching cubes MLS meshing done.");
|
||||
}
|
||||
|
||||
delete mls;
|
||||
if ( (id & _PROJECTION_) && par.getMesh("ControlMesh")!=pPoints)
|
||||
{
|
||||
delete pPoints;
|
||||
}
|
||||
|
||||
if (mesh)
|
||||
{
|
||||
cb(99, "Update face normals...");
|
||||
vcg::tri::UpdateNormals<CMeshO>::PerFace(mesh->cm);
|
||||
cb(100, "Update box...");
|
||||
vcg::tri::UpdateBounding<CMeshO>::Box(mesh->cm);
|
||||
}
|
||||
|
||||
} // end MLS based stuff
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
Q_EXPORT_PLUGIN(MlsPlugin)
|
||||
@ -1,76 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef MLSPLUGIN_H
|
||||
#define MLSPLUGIN_H
|
||||
|
||||
#include <QObject>
|
||||
|
||||
#include <meshlab/meshmodel.h>
|
||||
#include <meshlab/interfaces.h>
|
||||
|
||||
class MlsPlugin : public QObject, public MeshFilterInterface
|
||||
{
|
||||
Q_OBJECT
|
||||
Q_INTERFACES(MeshFilterInterface)
|
||||
|
||||
public:
|
||||
|
||||
enum {
|
||||
_RIMLS_ = 0x1,
|
||||
_APSS_ = 0x2,
|
||||
_PROJECTION_ = 0x1000,
|
||||
_AFRONT_ = 0x2000,
|
||||
_MCUBE_ = 0x4000,
|
||||
_COLORIZE_ = 0x8000,
|
||||
|
||||
FP_RIMLS_PROJECTION = _RIMLS_ | _PROJECTION_,
|
||||
FP_APSS_PROJECTION = _APSS_ | _PROJECTION_,
|
||||
|
||||
FP_RIMLS_AFRONT = _RIMLS_ | _AFRONT_,
|
||||
FP_APSS_AFRONT = _APSS_ | _AFRONT_,
|
||||
|
||||
FP_RIMLS_MCUBE = _RIMLS_ | _MCUBE_,
|
||||
FP_APSS_MCUBE = _APSS_ | _MCUBE_,
|
||||
|
||||
FP_RIMLS_COLORIZE = _RIMLS_ | _COLORIZE_,
|
||||
FP_APSS_COLORIZE = _APSS_ | _COLORIZE_,
|
||||
|
||||
FP_RADIUS_FROM_DENSITY = 0x10000,
|
||||
FP_SELECT_SMALL_COMPONENTS = 0x20000
|
||||
};
|
||||
|
||||
MlsPlugin();
|
||||
|
||||
virtual const QString filterName(FilterIDType filter);
|
||||
virtual const QString filterInfo(FilterIDType filter);
|
||||
const FilterClass getClass(QAction *a);
|
||||
virtual bool autoDialog(QAction *) {return true;}
|
||||
virtual void initParameterSet(QAction *,MeshDocument &md, FilterParameterSet &parent);
|
||||
virtual const int getRequirements(QAction *action);
|
||||
virtual bool applyFilter(QAction *filter, MeshDocument &m, FilterParameterSet &parent, vcg::CallBackPos *cb) ;
|
||||
virtual bool applyFilter(QAction *filter, MeshModel &m, FilterParameterSet &parent, vcg::CallBackPos *cb)
|
||||
{assert(0); return false;}
|
||||
};
|
||||
|
||||
#endif
|
||||
@ -1,203 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef MLSSURFACE_H
|
||||
#define MLSSURFACE_H
|
||||
|
||||
#include "balltree.h"
|
||||
#include <vcg/space/box3.h>
|
||||
#include <vcg/math/matrix33.h>
|
||||
#include <iostream>
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
enum {
|
||||
MLS_OK,
|
||||
MLS_TOO_FAR,
|
||||
MLS_TOO_MANY_ITERS,
|
||||
MLS_NOT_SUPPORTED,
|
||||
|
||||
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;
|
||||
|
||||
MlsSurface(const MeshType& mesh)
|
||||
: mMesh(mesh), mPoints(mesh.vert)
|
||||
{
|
||||
mCachedQueryPointIsOK = false;
|
||||
|
||||
mAABB = mesh.bbox;
|
||||
|
||||
// compute radii using a basic meshless density estimator
|
||||
if (!mPoints.RadiusEnabled)
|
||||
{
|
||||
const_cast<PointsType&>(mPoints).EnableRadius();
|
||||
computeVertexRaddi();
|
||||
}
|
||||
|
||||
mFilterScale = 4.0;
|
||||
mMaxNofProjectionIterations = 20;
|
||||
mProjectionAccuracy = 1e-4;
|
||||
mBallTree = 0;
|
||||
mGradientHint = MLS_DERIVATIVE_ACCURATE;
|
||||
mHessianHint = MLS_DERIVATIVE_ACCURATE;
|
||||
|
||||
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 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 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 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 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);
|
||||
|
||||
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].cP(), mPoints.size(),
|
||||
size_t(mPoints[1].cP().V()) - size_t(mPoints[0].cP().V()));
|
||||
}
|
||||
inline ConstDataWrapper<VectorType> normals() const
|
||||
{
|
||||
return ConstDataWrapper<VectorType>(&mPoints[0].cN(), mPoints.size(),
|
||||
size_t(mPoints[1].cN().V()) - size_t(mPoints[0].cN().V()));
|
||||
}
|
||||
inline ConstDataWrapper<Scalar> radii() const
|
||||
{
|
||||
return ConstDataWrapper<Scalar>(&mPoints[0].cR(), mPoints.size(),
|
||||
size_t(&mPoints[1].cR()) - size_t(&mPoints[0].cR()));
|
||||
}
|
||||
const vcg::Box3<Scalar>& boundingBox() const { return mAABB; }
|
||||
|
||||
static const Scalar InvalidValue() { return 12345679810.11121314151617; }
|
||||
|
||||
protected:
|
||||
void computeNeighborhood(const VectorType& x, bool computeDerivatives) const;
|
||||
void computeVertexRaddi();
|
||||
void requestSecondDerivatives() 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;
|
||||
|
||||
refD2 = d2;
|
||||
q = a;
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
class DummyObjectMarker {};
|
||||
|
||||
protected:
|
||||
const MeshType& mMesh;
|
||||
const PointsType& mPoints;
|
||||
vcg::Box3<Scalar> mAABB;
|
||||
int mGradientHint;
|
||||
int mHessianHint;
|
||||
|
||||
BallTree<Scalar>* mBallTree;
|
||||
|
||||
int mMaxNofProjectionIterations;
|
||||
Scalar mFilterScale;
|
||||
Scalar mAveragePointSpacing;
|
||||
Scalar mProjectionAccuracy;
|
||||
|
||||
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;
|
||||
};
|
||||
|
||||
} // namespace
|
||||
|
||||
#include "mlssurface.tpp"
|
||||
|
||||
#endif // MLSSURFACE_H
|
||||
@ -1,262 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include "mlssurface.h"
|
||||
#include "kdtree.h"
|
||||
#include "mlsutils.h"
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <vcg/space/index/octree.h>
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
// template<typename _Scalar>
|
||||
// MlsSurface<_Scalar>::MlsSurface(const MeshModel& m)
|
||||
// : mMesh(m)
|
||||
// {
|
||||
// mCachedQueryPointIsOK = false;
|
||||
//
|
||||
// mPoints.resize(m.cm.vert.size());
|
||||
// mNormals.resize(m.cm.vert.size());
|
||||
// mRadii.resize(m.cm.vert.size());
|
||||
// mAABB.Set(mMesh.cm.vert[0].cP());
|
||||
// for (uint i = 0; i< m.cm.vert.size(); i++)
|
||||
// {
|
||||
// mPoints[i] = /*vcg::vector_cast<Scalar>*/(mMesh.cm.vert[i].cP());
|
||||
// mNormals[i] = /*vcg::vector_cast<Scalar>*/(mMesh.cm.vert[i].cN());
|
||||
// mAABB.Add(mMesh.cm.vert[i].cP());
|
||||
// }
|
||||
//
|
||||
// // compute radii using a basic meshless density estimator
|
||||
// computeVertexRaddi();
|
||||
//
|
||||
// mFilterScale = 4.0;
|
||||
// mMaxNofProjectionIterations = 20;
|
||||
// mProjectionAccuracy = 1e-4;
|
||||
// mBallTree = 0;
|
||||
// }
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setFilterScale(Scalar v)
|
||||
{
|
||||
mFilterScale = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
if (mBallTree)
|
||||
mBallTree->setRadiusScale(mFilterScale);
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setMaxProjectionIters(int n)
|
||||
{
|
||||
mMaxNofProjectionIterations = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setProjectionAccuracy(Scalar v)
|
||||
{
|
||||
mProjectionAccuracy = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setGradientHint(int h)
|
||||
{
|
||||
mGradientHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::setHessianHint(int h)
|
||||
{
|
||||
mHessianHint = h;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void MlsSurface<_MeshType>::computeVertexRaddi()
|
||||
{
|
||||
#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());
|
||||
|
||||
#else
|
||||
|
||||
int nbNeighbors = 16;
|
||||
|
||||
assert(mPoints.size()>=2);
|
||||
KdTree<Scalar> knn(positions());
|
||||
|
||||
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
|
||||
}
|
||||
|
||||
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();
|
||||
|
||||
// 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;
|
||||
|
||||
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);
|
||||
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
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>;
|
||||
// template class MlsSurface<double>;
|
||||
|
||||
}
|
||||
@ -1,140 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef VCGADDONS_H
|
||||
#define VCGADDONS_H
|
||||
|
||||
template<typename _DataType>
|
||||
class ConstDataWrapper
|
||||
{
|
||||
public:
|
||||
typedef _DataType DataType;
|
||||
inline ConstDataWrapper()
|
||||
: mpData(0), mStride(0), mSize(0)
|
||||
{}
|
||||
inline ConstDataWrapper(const DataType* pData, int size, int stride = sizeof(DataType))
|
||||
: mpData(reinterpret_cast<const unsigned char*>(pData)), mStride(stride), mSize(size)
|
||||
{}
|
||||
inline const DataType& operator[] (int i) const
|
||||
{
|
||||
return *reinterpret_cast<const DataType*>(mpData + i*mStride);
|
||||
}
|
||||
inline int size() const { return mSize; }
|
||||
protected:
|
||||
const unsigned char* mpData;
|
||||
int mStride;
|
||||
int mSize;
|
||||
};
|
||||
|
||||
namespace vcg {
|
||||
|
||||
template <typename Scalar>
|
||||
inline Point3<Scalar> CwiseMul(Point3<Scalar> const & p1, Point3<Scalar> const & p2)
|
||||
{
|
||||
return Point3<Scalar>(p1.X()*p2.X(), p1.Y()*p2.Y(), p1.Z()*p2.Z());
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Point3<Scalar> Min(Point3<Scalar> const & p1, Point3<Scalar> const & p2)
|
||||
{
|
||||
return Point3<Scalar>(std::min(p1.X(), p2.X()), std::min(p1.Y(), p2.Y()), std::min(p1.Z(), p2.Z()));
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Point3<Scalar> Max(Point3<Scalar> const & p1, Point3<Scalar> const & p2)
|
||||
{
|
||||
return Point3<Scalar>(std::max(p1.X(), p2.X()), std::max(p1.Y(), p2.Y()), std::max(p1.Z(), p2.Z()));
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Scalar MaxCoeff(Point3<Scalar> const & p)
|
||||
{
|
||||
return std::max(std::max(p.X(), p.Y()), p.Z());
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Scalar MinCoeff(Point3<Scalar> const & p)
|
||||
{
|
||||
return std::min(std::min(p.X(), p.Y()), p.Z());
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Scalar Dot(Point3<Scalar> const & p1, Point3<Scalar> const & p2)
|
||||
{
|
||||
return p1.X() * p2.X() + p1.Y() * p2.Y() + p1.Z() * p2.Z();
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Point3<Scalar> Cross(Point3<Scalar> const & p1, Point3<Scalar> const & p2)
|
||||
{
|
||||
return p1 ^ p2;
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline Point3<Scalar> CwiseAdd(Point3<Scalar> const & p1, Scalar s)
|
||||
{
|
||||
return Point3<Scalar>(p1.X() + s, p1.Y() + s, p1.Z() + s);
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline int MaxCoeffId(Point3<Scalar> const & p)
|
||||
{
|
||||
if (p.X()>p.Y())
|
||||
return p.X()>p.Z() ? 0 : 2;
|
||||
else
|
||||
return p.Y()>p.Z() ? 1 : 2;
|
||||
}
|
||||
|
||||
template <typename Scalar>
|
||||
inline int MinCoeffId(Point3<Scalar> const & p)
|
||||
{
|
||||
if (p.X()<p.Y())
|
||||
return p.X()<p.Z() ? 0 : 2;
|
||||
else
|
||||
return p.Y()<p.Z() ? 1 : 2;
|
||||
}
|
||||
|
||||
template <typename ToType, typename Scalar>
|
||||
inline Point3<ToType> Point3Cast(const Point3<Scalar>& p)
|
||||
{
|
||||
return Point3<ToType>(p.X(), p.Y(), p.Z());
|
||||
}
|
||||
|
||||
template<class Scalar>
|
||||
Scalar Distance(const Point3<Scalar> &p, const Box3<Scalar> &bbox)
|
||||
{
|
||||
Scalar dist2 = 0.;
|
||||
Scalar aux;
|
||||
for (int k=0 ; k<3 ; ++k)
|
||||
{
|
||||
if ( (aux = (p[k]-bbox.min[k]))<0. )
|
||||
dist2 += aux*aux;
|
||||
else if ( (aux = (bbox.max[k]-p[k]))<0. )
|
||||
dist2 += aux*aux;
|
||||
}
|
||||
return sqrt(dist2);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -1,122 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef _PriorityQueue_h_
|
||||
#define _PriorityQueue_h_
|
||||
|
||||
/** Implements a bounded-size max priority queue using a heap
|
||||
*/
|
||||
template <typename Index, typename Weight>
|
||||
class HeapMaxPriorityQueue
|
||||
{
|
||||
struct Element
|
||||
{
|
||||
Weight weight;
|
||||
Index index;
|
||||
};
|
||||
|
||||
public:
|
||||
|
||||
HeapMaxPriorityQueue(void)
|
||||
{
|
||||
mElements = 0;
|
||||
mMaxSize = 0;
|
||||
}
|
||||
|
||||
inline void setMaxSize(int maxSize)
|
||||
{
|
||||
if (mMaxSize!=maxSize)
|
||||
{
|
||||
mMaxSize = maxSize;
|
||||
delete[] mElements;
|
||||
mElements = new Element[mMaxSize];
|
||||
mpOffsetedElements = (mElements-1);
|
||||
}
|
||||
init();
|
||||
}
|
||||
|
||||
inline void init() { mCount = 0; }
|
||||
|
||||
inline bool isFull() const { return mCount == mMaxSize; }
|
||||
|
||||
/** returns number of elements inserted in queue
|
||||
*/
|
||||
inline int getNofElements() const { return mCount; }
|
||||
|
||||
inline Weight getWeight(int i) const { return mElements[i].weight; }
|
||||
inline Index getIndex(int i) const { return mElements[i].index; }
|
||||
|
||||
inline Weight getTopWeight() const { return mElements[0].weight; }
|
||||
|
||||
inline void insert(Index index, Weight weight)
|
||||
{
|
||||
if (mCount==mMaxSize)
|
||||
{
|
||||
if (weight<mElements[0].weight)
|
||||
{
|
||||
register int j, k;
|
||||
j = 1;
|
||||
k = 2;
|
||||
while (k <= mMaxSize)
|
||||
{
|
||||
Element* z = &(mpOffsetedElements[k]);
|
||||
if ((k < mMaxSize) && (z->weight < mpOffsetedElements[k+1].weight))
|
||||
z = &(mpOffsetedElements[++k]);
|
||||
|
||||
if(weight >= z->weight)
|
||||
break;
|
||||
mpOffsetedElements[j] = *z;
|
||||
j = k;
|
||||
k = 2 * j;
|
||||
}
|
||||
mpOffsetedElements[j].weight = weight;
|
||||
mpOffsetedElements[j].index = index;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
int i, j;
|
||||
i = ++mCount;
|
||||
while (i >= 2)
|
||||
{
|
||||
j = i >> 1;
|
||||
Element& y = mpOffsetedElements[j];
|
||||
if(weight <= y.weight)
|
||||
break;
|
||||
mpOffsetedElements[i] = y;
|
||||
i = j;
|
||||
}
|
||||
mpOffsetedElements[i].index = index;
|
||||
mpOffsetedElements[i].weight = weight;
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
|
||||
int mCount;
|
||||
int mMaxSize;
|
||||
Element* mElements;
|
||||
Element* mpOffsetedElements;
|
||||
};
|
||||
|
||||
#endif
|
||||
@ -1,11 +0,0 @@
|
||||
|
||||
#include <qglobal.h>
|
||||
#include "rimls.h"
|
||||
#include "rimls.tpp"
|
||||
#include <meshlab/meshmodel.h>
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template class RIMLS<CMeshO>;
|
||||
|
||||
}
|
||||
@ -1,104 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef RIMLS_H
|
||||
#define RIMLS_H
|
||||
|
||||
#include "mlssurface.h"
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
template<typename _MeshType>
|
||||
class RIMLS : public MlsSurface<_MeshType>
|
||||
{
|
||||
typedef _MeshType MeshType;
|
||||
typedef MlsSurface<_MeshType> Base;
|
||||
|
||||
typedef typename Base::Scalar Scalar;
|
||||
typedef typename Base::VectorType VectorType;
|
||||
typedef typename Base::MatrixType MatrixType;
|
||||
using Base::mCachedQueryPointIsOK;
|
||||
using Base::mCachedQueryPoint;
|
||||
using Base::mNeighborhood;
|
||||
using Base::mCachedWeights;
|
||||
using Base::mCachedWeightDerivatives;
|
||||
using Base::mCachedWeightGradients;
|
||||
using Base::mCachedWeightSecondDerivatives;
|
||||
using Base::mBallTree;
|
||||
using Base::mPoints;
|
||||
using Base::mFilterScale;
|
||||
using Base::mMaxNofProjectionIterations;
|
||||
using Base::mAveragePointSpacing;
|
||||
using Base::mProjectionAccuracy;
|
||||
|
||||
public:
|
||||
|
||||
RIMLS(const MeshType& points)
|
||||
: Base(points)
|
||||
{
|
||||
mSigmaR = 0;
|
||||
mSigmaN = 0.8;
|
||||
mRefittingThreshold = 1e-3;
|
||||
mMinRefittingIters = 1;
|
||||
mMaxRefittingIters = 3;
|
||||
}
|
||||
|
||||
virtual Scalar potential(const VectorType& x, int* errorMask = 0) const;
|
||||
virtual VectorType gradient(const VectorType& x, int* errorMask = 0) const;
|
||||
virtual MatrixType hessian(const VectorType& x, int* errorMask = 0) const;
|
||||
virtual VectorType project(const VectorType& x, VectorType* pNormal = 0, int* errorMask = 0) const;
|
||||
|
||||
void setSigmaR(Scalar v);
|
||||
void setSigmaN(Scalar v);
|
||||
void setRefittingThreshold(Scalar v);
|
||||
void setMinRefittingIters(int n);
|
||||
void setMaxRefittingIters(int n);
|
||||
|
||||
protected:
|
||||
bool computePotentialAndGradient(const VectorType& x) const;
|
||||
bool mlsHessian(const VectorType& x, MatrixType& hessian) const;
|
||||
|
||||
protected:
|
||||
|
||||
int mMinRefittingIters;
|
||||
int mMaxRefittingIters;
|
||||
Scalar mRefittingThreshold;
|
||||
Scalar mSigmaN;
|
||||
Scalar mSigmaR;
|
||||
|
||||
// cached values:
|
||||
mutable VectorType mCachedGradient;
|
||||
mutable Scalar mCachedPotential;
|
||||
|
||||
mutable Scalar mCachedSumW;
|
||||
mutable std::vector<Scalar> mCachedRefittingWeights;;
|
||||
mutable VectorType mCachedSumN;
|
||||
mutable VectorType mCachedSumGradWeight;
|
||||
mutable VectorType mCachedSumGradPotential;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
//#include "rimls.tpp"
|
||||
|
||||
#endif
|
||||
@ -1,307 +0,0 @@
|
||||
/****************************************************************************
|
||||
* MeshLab o o *
|
||||
* A versatile mesh processing toolbox o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2005 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#include "rimls.h"
|
||||
#include "kdtree.h"
|
||||
#include "mlsutils.h"
|
||||
#include <iostream>
|
||||
|
||||
using namespace vcg;
|
||||
|
||||
namespace GaelMls {
|
||||
|
||||
// template<typename _Scalar>
|
||||
// RIMLS<_Scalar>::RIMLS(const MeshModel& m)
|
||||
// : Base(m)
|
||||
// {
|
||||
// mSigmaR = 0;
|
||||
// mSigmaN = 0.8;
|
||||
// mRefittingThreshold = 1e-3;
|
||||
// mMinRefittingIters = 1;
|
||||
// mMaxRefittingIters = 3;
|
||||
// }
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setSigmaR(Scalar v)
|
||||
{
|
||||
mSigmaR = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setSigmaN(Scalar v)
|
||||
{
|
||||
mSigmaN = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setRefittingThreshold(Scalar v)
|
||||
{
|
||||
mRefittingThreshold = v;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setMinRefittingIters(int n)
|
||||
{
|
||||
mMinRefittingIters = n;
|
||||
mCachedQueryPointIsOK = false;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
void RIMLS<_MeshType>::setMaxRefittingIters(int n)
|
||||
{
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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 (pNormal)
|
||||
*pNormal = normal;
|
||||
|
||||
return position;
|
||||
}
|
||||
|
||||
template<typename _MeshType>
|
||||
bool RIMLS<_MeshType>::computePotentialAndGradient(const VectorType& x) const
|
||||
{
|
||||
Base::computeNeighborhood(x, true);
|
||||
unsigned int nofSamples = mNeighborhood.size();
|
||||
|
||||
if (nofSamples<1)
|
||||
{
|
||||
mCachedGradient.SetZero();
|
||||
mCachedQueryPoint = x;
|
||||
mCachedPotential = 1e9;
|
||||
mCachedQueryPointIsOK = false;
|
||||
return false;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
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);
|
||||
|
||||
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;
|
||||
|
||||
sumGradWeight += gw;
|
||||
sumGradWeightPotential += gw * f;
|
||||
sumN += normal * w;
|
||||
potential += w * f;
|
||||
sumW += w;
|
||||
}
|
||||
|
||||
if(sumW==0.)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
potential /= sumW;
|
||||
grad = (-sumGradWeight*potential + sumGradWeightPotential + sumN) * (1./sumW);
|
||||
|
||||
iterationCount++;
|
||||
|
||||
} while ( (iterationCount < mMinRefittingIters)
|
||||
|| ( vcg::SquaredNorm(grad - previousGrad) > mRefittingThreshold && iterationCount < mMaxRefittingIters) );
|
||||
|
||||
mCachedGradient = grad;
|
||||
mCachedPotential = potential;
|
||||
mCachedQueryPoint = x;
|
||||
mCachedQueryPointIsOK = true;
|
||||
|
||||
mCachedSumGradWeight = sumGradWeight;
|
||||
mCachedSumN = sumN;
|
||||
mCachedSumW = sumW;
|
||||
mCachedSumGradPotential = sumGradWeightPotential;
|
||||
|
||||
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
|
||||
|
||||
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;
|
||||
|
||||
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());
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
VectorType dGrad = (
|
||||
sumDWeightNormal + sumGradWeightNk + sumDGradWeightPotential
|
||||
- sumDGradWeight * mCachedPotential
|
||||
- sumGradWeight * mCachedGradient[k]
|
||||
- mCachedGradient * sumGradWeight[k] ) * invW;
|
||||
|
||||
hessian.SetColumn(k,dGrad);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -1,150 +0,0 @@
|
||||
/****************************************************************************
|
||||
* VCGLib o o *
|
||||
* Visual and Computer Graphics Library o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2004 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
|
||||
#ifndef SMALL_COMPONENT_SELECTION_H
|
||||
#define SMALL_COMPONENT_SELECTION_H
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
|
||||
template<class _MeshType>
|
||||
class SmallComponent
|
||||
{
|
||||
|
||||
public:
|
||||
typedef _MeshType MeshType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::VertexPointer VertexPointer;
|
||||
typedef typename MeshType::VertexIterator VertexIterator;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::FacePointer FacePointer;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
|
||||
static int Select(MeshType &m, float nbFaceRatio = 0.1, bool nonClosedOnly = false)
|
||||
{
|
||||
assert(tri::HasFFAdjacency(m) && "The small component selection procedure requires face to face adjacency.");
|
||||
|
||||
// the different components as a list of face pointer
|
||||
std::vector< std::vector<FacePointer> > components;
|
||||
|
||||
for(int faceSeed = 0; faceSeed<m.face.size(); )
|
||||
{
|
||||
// find the first not selected face
|
||||
bool foundSeed = false;
|
||||
while (faceSeed<m.face.size())
|
||||
{
|
||||
CFaceO& f = m.face[faceSeed];
|
||||
if (!f.IsS())
|
||||
{
|
||||
if (nonClosedOnly)
|
||||
{
|
||||
for (int k=0; k<3; ++k)
|
||||
if (f.IsB(k))
|
||||
{
|
||||
foundSeed = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
else
|
||||
foundSeed = true;
|
||||
if (foundSeed)
|
||||
break;
|
||||
}
|
||||
++faceSeed;
|
||||
}
|
||||
if (!foundSeed) // no more seed, stop
|
||||
break;
|
||||
|
||||
// expand the region from this face...
|
||||
components.resize(components.size()+1);
|
||||
std::vector<FacePointer> activefaces;
|
||||
activefaces.push_back(&m.face[faceSeed]);
|
||||
while (!activefaces.empty())
|
||||
{
|
||||
FacePointer f = activefaces.back();
|
||||
activefaces.pop_back();
|
||||
if (f->IsS())
|
||||
continue;
|
||||
|
||||
f->SetS();
|
||||
components.back().push_back(f);
|
||||
for (int k=0; k<3; ++k)
|
||||
{
|
||||
if (f->IsB(k))
|
||||
continue;
|
||||
FacePointer of = f->FFp(k);
|
||||
if (!of->IsS())
|
||||
activefaces.push_back(of);
|
||||
}
|
||||
}
|
||||
++faceSeed;
|
||||
}
|
||||
UpdateSelection<MeshType>::ClearFace(m);
|
||||
|
||||
// now the segmentation is done, let's compute the absolute face count threshold
|
||||
int total_selected = 0;
|
||||
int maxComponent = 0;
|
||||
for (int i=0; i<components.size(); ++i)
|
||||
{
|
||||
//std::cout << "Component " << i << " -> " << components[i].size() << "\n";
|
||||
total_selected += components[i].size();
|
||||
maxComponent = std::max<int>(maxComponent,components[i].size());
|
||||
}
|
||||
int remaining = m.face.size() - total_selected;
|
||||
int th = std::max(maxComponent,remaining) * nbFaceRatio;
|
||||
|
||||
int selCount = 0;
|
||||
for (int i=0; i<components.size(); ++i)
|
||||
{
|
||||
if (components[i].size()<th)
|
||||
{
|
||||
selCount += components[i].size();
|
||||
for (int j=0; j<components[i].size(); ++j)
|
||||
components[i][j]->SetS();
|
||||
}
|
||||
}
|
||||
return selCount;
|
||||
}
|
||||
|
||||
static void DeleteFaceVert(MeshType &m)
|
||||
{
|
||||
typename MeshType::FaceIterator fi;
|
||||
typename MeshType::VertexIterator vi;
|
||||
UpdateSelection<MeshType>::ClearVertex(m);
|
||||
UpdateSelection<MeshType>::VertexFromFaceStrict(m);
|
||||
|
||||
for(fi=m.face.begin();fi!=m.face.end();++fi)
|
||||
if (!(*fi).IsD() && (*fi).IsS() )
|
||||
tri::Allocator<CMeshO>::DeleteFace(m,*fi);
|
||||
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
|
||||
if (!(*vi).IsD() && (*vi).IsS() )
|
||||
tri::Allocator<CMeshO>::DeleteVertex(m,*vi);
|
||||
}
|
||||
|
||||
}; // end class
|
||||
|
||||
} // End namespace
|
||||
} // End namespace
|
||||
|
||||
|
||||
#endif
|
||||
Loading…
x
Reference in New Issue
Block a user