From c321c42e8731e1bfe8fb249738f75c2edc9b2baa Mon Sep 17 00:00:00 2001 From: alemuntoni Date: Thu, 14 Oct 2021 17:09:37 +0200 Subject: [PATCH] filter mls uses a custom per vertex attribute for the radius --- src/meshlabplugins/filter_mls/mlsplugin.cpp | 12 +--- src/meshlabplugins/filter_mls/mlsplugin.h | 2 - src/meshlabplugins/filter_mls/mlssurface.h | 22 ++++--- src/meshlabplugins/filter_mls/mlssurface.tpp | 65 +++++++++++++------- src/vcglib | 2 +- 5 files changed, 61 insertions(+), 42 deletions(-) diff --git a/src/meshlabplugins/filter_mls/mlsplugin.cpp b/src/meshlabplugins/filter_mls/mlsplugin.cpp index 0c91ab28c..6d86b1966 100644 --- a/src/meshlabplugins/filter_mls/mlsplugin.cpp +++ b/src/meshlabplugins/filter_mls/mlsplugin.cpp @@ -316,9 +316,7 @@ std::map MlsPlugin::applyFilter( computeColorize(md, par, mls, pPoints, cb); break; case FP_RADIUS_FROM_DENSITY: { - md.mm()->updateDataMask(MeshModel::MM_VERTRADIUS); - APSS mls(md.mm()->cm); - mls.computeVertexRaddi(par.getInt("NbNeighbors")); + GaelMls::computeVertexRadius(md.mm()->cm, par.getInt("NbNeighbors")); break; } case FP_SELECT_SMALL_COMPONENTS: @@ -466,13 +464,7 @@ void MlsPlugin::initMLS(MeshDocument& md) } 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"); - } + GaelMls::computeVertexRadius(md.mm()->cm); } MeshModel* MlsPlugin::getProjectionPointsMesh(MeshDocument& md, const RichParameterList& params) diff --git a/src/meshlabplugins/filter_mls/mlsplugin.h b/src/meshlabplugins/filter_mls/mlsplugin.h index 14db8407f..180fd1399 100644 --- a/src/meshlabplugins/filter_mls/mlsplugin.h +++ b/src/meshlabplugins/filter_mls/mlsplugin.h @@ -71,8 +71,6 @@ private: void addColorizeParameters(RichParameterList& parlst, bool apss); void addMarchingCubesParameters(RichParameterList& parlst); - - void initMLS(MeshDocument& md); MeshModel* getProjectionPointsMesh(MeshDocument& md, const RichParameterList& params); GaelMls::MlsSurface* createMlsRimls(MeshModel* pPoints, const RichParameterList& par); diff --git a/src/meshlabplugins/filter_mls/mlssurface.h b/src/meshlabplugins/filter_mls/mlssurface.h index e13c4c101..f5082fdb1 100644 --- a/src/meshlabplugins/filter_mls/mlssurface.h +++ b/src/meshlabplugins/filter_mls/mlssurface.h @@ -29,9 +29,13 @@ #include #include #include +#include namespace GaelMls { +template +void computeVertexRadius(MeshType& m, int nNeighbors = 16); + enum { MLS_OK, MLS_TOO_FAR, @@ -58,11 +62,9 @@ public: mAABB = mesh.bbox; - // compute radii using a basic meshless density estimator - if (!mMesh.vert.RadiusEnabled) { - const_cast(mMesh.vert).EnableRadius(); - computeVertexRaddi(); - } + typename MeshType::template ConstPerVertexAttributeHandle h; + h = vcg::tri::Allocator::template FindPerVertexAttribute(mMesh, "radius"); + assert(vcg::tri::Allocator::IsValidHandle(mMesh, h)); mFilterScale = 4.0; mMaxNofProjectionIterations = 20; @@ -149,16 +151,20 @@ public: } inline vcg::ConstDataWrapper radii() const { + typename MeshType::template ConstPerVertexAttributeHandle h; + h = vcg::tri::Allocator::template FindPerVertexAttribute(mMesh, "radius"); + assert(vcg::tri::Allocator<_MeshType>::template IsValidHandle(mMesh, h)); + return vcg::ConstDataWrapper( - &mMesh.vert[0].R(), + &h[0], mMesh.vert.size(), - size_t(&mMesh.vert[1].R()) - size_t(&mMesh.vert[0].R())); + size_t(&h[1]) - size_t(&h[0])); } const vcg::Box3& boundingBox() const { return mAABB; } static const Scalar InvalidValue() { return Scalar(12345679810.11121314151617); } - void computeVertexRaddi(const int nbNeighbors = 16); + //void computeVertexRaddi(const int nbNeighbors = 16); protected: void computeNeighborhood(const VectorType& x, bool computeDerivatives) const; diff --git a/src/meshlabplugins/filter_mls/mlssurface.tpp b/src/meshlabplugins/filter_mls/mlssurface.tpp index b0613c59e..805769dbe 100644 --- a/src/meshlabplugins/filter_mls/mlssurface.tpp +++ b/src/meshlabplugins/filter_mls/mlssurface.tpp @@ -26,9 +26,36 @@ #include #include #include +#include +#include namespace GaelMls { +template +void computeVertexRadius(_MeshType& mesh, int nNeighbors) +{ + typedef typename _MeshType::ScalarType Scalar; + if (!vcg::tri::HasPerVertexAttribute(mesh, "radius")) { + vcg::tri::Allocator<_MeshType>::template AddPerVertexAttribute(mesh, "radius"); + } + + typename _MeshType::template PerVertexAttributeHandle h; + h = vcg::tri::Allocator<_MeshType>::template FindPerVertexAttribute(mesh, "radius"); + assert(vcg::tri::Allocator<_MeshType>::template IsValidHandle(mesh, h)); + + auto positions = vcg::ConstDataWrapper>( + &mesh.vert[0].P(), + mesh.vert.size(), + size_t(mesh.vert[1].P().V()) - size_t(mesh.vert[0].P().V())); + + vcg::KdTree knn(positions); + typename vcg::KdTree::PriorityQueue pq; + for (size_t i = 0; i < mesh.vert.size(); i++) { + knn.doQueryK(mesh.vert[i].cP(), nNeighbors, pq); + h[i] = 2. * sqrt(pq.getTopWeight() / Scalarm(pq.getNofElements())); + } +} + template void MlsSurface<_MeshType>::setFilterScale(Scalar v) { @@ -66,23 +93,6 @@ void MlsSurface<_MeshType>::setHessianHint(int h) mCachedQueryPointIsOK = false; } -template -void MlsSurface<_MeshType>::computeVertexRaddi(const int nbNeighbors) -{ - assert(mMesh.vert.size() >= 2); - vcg::KdTree knn(positions()); - typename vcg::KdTree::PriorityQueue pq; - // knn.setMaxNofNeighbors(nbNeighbors); - mAveragePointSpacing = 0; - for (size_t i = 0; i < mMesh.vert.size(); i++) { - knn.doQueryK(mMesh.vert[i].cP(), nbNeighbors, pq); - const_cast(mMesh.vert)[i].R() = - 2. * sqrt(pq.getTopWeight() / Scalar(pq.getNofElements())); - mAveragePointSpacing += mMesh.vert[i].cR(); - } - mAveragePointSpacing /= Scalar(mMesh.vert.size()); -} - template void MlsSurface<_MeshType>::computeNeighborhood(const VectorType& x, bool computeDerivatives) const { @@ -102,9 +112,13 @@ void MlsSurface<_MeshType>::computeNeighborhood(const VectorType& x, bool comput else mCachedWeightGradients.clear(); + typename _MeshType::template ConstPerVertexAttributeHandle h; + h = vcg::tri::Allocator<_MeshType>::template FindPerVertexAttribute(mMesh, "radius"); + assert(vcg::tri::Allocator<_MeshType>::template IsValidHandle(mMesh, h)); + for (size_t i = 0; i < nofSamples; i++) { int id = mNeighborhood.index(i); - Scalar s = 1. / (mMesh.vert[id].cR() * mFilterScale); + Scalar s = 1. / (h[id] * mFilterScale); s = s * s; Scalar w = Scalar(1) - mNeighborhood.squaredDistance(i) * s; if (w < 0) @@ -124,6 +138,10 @@ void MlsSurface<_MeshType>::computeNeighborhood(const VectorType& x, bool comput template void MlsSurface<_MeshType>::requestSecondDerivatives() const { + typename _MeshType::template ConstPerVertexAttributeHandle h; + h = vcg::tri::Allocator<_MeshType>::template FindPerVertexAttribute(mMesh, "radius"); + assert(vcg::tri::Allocator<_MeshType>::template IsValidHandle(mMesh, h)); + // if (!mSecondDerivativeUptodate) { size_t nofSamples = mNeighborhood.size(); @@ -133,7 +151,7 @@ void MlsSurface<_MeshType>::requestSecondDerivatives() const { for (size_t i = 0; i < nofSamples; ++i) { int id = mNeighborhood.index(i); - Scalar s = 1. / (mMesh.vert[id].cR() * mFilterScale); + Scalar s = 1. / (h[id] * mFilterScale); s = s * s; Scalar x2 = s * mNeighborhood.squaredDistance(i); x2 = 1.0 - x2; @@ -160,6 +178,11 @@ MlsSurface<_MeshType>::meanCurvature(const VectorType& gradient, const MatrixTyp template bool MlsSurface<_MeshType>::isInDomain(const VectorType& x) const { + typename _MeshType::template ConstPerVertexAttributeHandle h; + h = vcg::tri::Allocator<_MeshType>::template FindPerVertexAttribute(mMesh, "radius"); + assert(vcg::tri::Allocator<_MeshType>::template IsValidHandle(mMesh, h)); + + if ((!mCachedQueryPointIsOK) || mCachedQueryPoint != x) { computeNeighborhood(x, false); } @@ -173,7 +196,7 @@ bool MlsSurface<_MeshType>::isInDomain(const VectorType& x) const if ((mDomainNormalScale == 1.f) || (!hasNormal)) { while (out && i < nb) { int id = mNeighborhood.index(i); - Scalar rs2 = mMesh.vert[id].cR() * mDomainRadiusScale; + Scalar rs2 = h[id] * mDomainRadiusScale; rs2 = rs2 * rs2; out = mNeighborhood.squaredDistance(i) > rs2; ++i; @@ -183,7 +206,7 @@ bool MlsSurface<_MeshType>::isInDomain(const VectorType& x) const Scalar s = 1. / (mDomainNormalScale * mDomainNormalScale) - 1.f; while (out && i < nb) { int id = mNeighborhood.index(i); - Scalar rs2 = mMesh.vert[id].cR() * mDomainRadiusScale; + Scalar rs2 = h[id] * mDomainRadiusScale; rs2 = rs2 * rs2; Scalar dn = mMesh.vert[id].cN().dot(x - mMesh.vert[id].cP()); out = (mNeighborhood.squaredDistance(i) + s * dn * dn) > rs2; diff --git a/src/vcglib b/src/vcglib index fa6bab759..3bb6cfc71 160000 --- a/src/vcglib +++ b/src/vcglib @@ -1 +1 @@ -Subproject commit fa6bab75912840168c7814d9157bf11bc5be4224 +Subproject commit 3bb6cfc71aa3081c4f8a30bce45df3bfcd0ce547