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