1 /**************************************************************************** 2 * MeshLab o o * 3 * A versatile mesh processing toolbox o o * 4 * _ O _ * 5 * Copyright(C) 2005 \/)\/ * 6 * Visual Computing Lab /\/| * 7 * ISTI - Italian National Research Council | * 8 * \ * 9 * All rights reserved. * 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 * This program is distributed in the hope that it will be useful, * 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * 20 * for more details. * 21 * * 22 ****************************************************************************/ 23 24 #ifndef MLSSURFACE_H 25 #define MLSSURFACE_H 26 27 #include "balltree.h" 28 #include <vcg/space/box3.h> 29 #include <vcg/math/matrix33.h> 30 #include <Eigen/Dense> 31 #include <iostream> 32 33 namespace GaelMls { 34 35 enum { 36 MLS_OK, 37 MLS_TOO_FAR, 38 MLS_TOO_MANY_ITERS, 39 MLS_NOT_SUPPORTED, 40 41 MLS_DERIVATIVE_ACCURATE, 42 MLS_DERIVATIVE_APPROX, 43 MLS_DERIVATIVE_FINITEDIFF 44 }; 45 46 template<typename MeshType> 47 class MlsSurface 48 { 49 public: 50 typedef typename MeshType::ScalarType Scalar; 51 typedef vcg::Point3<Scalar> VectorType; 52 typedef vcg::Matrix33<Scalar> MatrixType; 53 typedef typename MeshType::VertContainer PointsType; 54 MlsSurface(const MeshType & mesh)55 MlsSurface(const MeshType& mesh) 56 : mMesh(mesh), mPoints(mesh.vert) 57 { 58 mCachedQueryPointIsOK = false; 59 60 mAABB = mesh.bbox; 61 62 // compute radii using a basic meshless density estimator 63 if (!mPoints.RadiusEnabled) 64 { 65 const_cast<PointsType&>(mPoints).EnableRadius(); 66 computeVertexRaddi(); 67 } 68 69 mFilterScale = 4.0; 70 mMaxNofProjectionIterations = 20; 71 mProjectionAccuracy = (Scalar)1e-4; 72 mBallTree = 0; 73 mGradientHint = MLS_DERIVATIVE_ACCURATE; 74 mHessianHint = MLS_DERIVATIVE_ACCURATE; 75 76 mDomainMinNofNeighbors = 4; 77 mDomainRadiusScale = 2.; 78 mDomainNormalScale = 1.; 79 } 80 ~MlsSurface()81 virtual ~MlsSurface() {} 82 83 /** \returns the value of the reconstructed scalar field at point \a x */ 84 virtual Scalar potential(const VectorType& x, int* errorMask = 0) const = 0; 85 86 /** \returns the gradient of the reconstructed scalar field at point \a x 87 * 88 * The method used to compute the gradient can be controlled with setGradientHint(). 89 */ 90 virtual VectorType gradient(const VectorType& x, int* errorMask = 0) const = 0; 91 92 /** \returns the hessian matrix of the reconstructed scalar field at point \a x 93 * 94 * The method used to compute the hessian matrix can be controlled with setHessianHint(). 95 */ 96 virtual MatrixType hessian(const VectorType& x, int* errorMask = 0) const 97 { if (errorMask) *errorMask = MLS_NOT_SUPPORTED; return MatrixType(); } 98 99 /** \returns the projection of point x onto the MLS surface, and optionally returns the normal in \a pNormal */ 100 virtual VectorType project(const VectorType& x, VectorType* pNormal = 0, int* errorMask = 0) const = 0; 101 102 /** \returns whether \a x is inside the restricted surface definition domain */ 103 virtual bool isInDomain(const VectorType& x) const; 104 105 /** \returns the mean curvature from the gradient vector and Hessian matrix. 106 */ 107 Scalar meanCurvature(const VectorType& gradient, const MatrixType& hessian) const; 108 109 /** set the scale of the spatial filter */ 110 void setFilterScale(Scalar v); 111 /** set the maximum number of iterations during the projection */ 112 void setMaxProjectionIters(int n); 113 /** set the threshold factor to detect convergence of the iterations */ 114 void setProjectionAccuracy(Scalar v); 115 116 /** set a hint on how to compute the gradient 117 * 118 * Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF 119 */ 120 void setGradientHint(int h); 121 122 /** set a hint on how to compute the hessian matrix 123 * 124 * Possible values are MLS_DERIVATIVE_ACCURATE, MLS_DERIVATIVE_APPROX, MLS_DERIVATIVE_FINITEDIFF 125 */ 126 void setHessianHint(int h); 127 mesh()128 inline const MeshType& mesh() const { return mMesh; } 129 /** a shortcut for mesh().vert */ points()130 inline const PointsType& points() const { return mPoints; } 131 positions()132 inline vcg::ConstDataWrapper<VectorType> positions() const 133 { 134 return vcg::ConstDataWrapper<VectorType>(&mPoints[0].P(), mPoints.size(), 135 size_t(mPoints[1].P().V()) - size_t(mPoints[0].P().V())); 136 } normals()137 inline vcg::ConstDataWrapper<VectorType> normals() const 138 { 139 return vcg::ConstDataWrapper<VectorType>(&mPoints[0].N(), mPoints.size(), 140 size_t(mPoints[1].N().V()) - size_t(mPoints[0].N().V())); 141 } radii()142 inline vcg::ConstDataWrapper<Scalar> radii() const 143 { 144 return vcg::ConstDataWrapper<Scalar>(&mPoints[0].R(), mPoints.size(), 145 size_t(&mPoints[1].R()) - size_t(&mPoints[0].R())); 146 } boundingBox()147 const vcg::Box3<Scalar>& boundingBox() const { return mAABB; } 148 InvalidValue()149 static const Scalar InvalidValue() { return Scalar(12345679810.11121314151617); } 150 151 void computeVertexRaddi(const int nbNeighbors = 16); 152 protected: 153 void computeNeighborhood(const VectorType& x, bool computeDerivatives) const; 154 void requestSecondDerivatives() const; 155 156 struct PointToPointSqDist 157 { operatorPointToPointSqDist158 inline bool operator()(const VectorType &a, const VectorType &b, Scalar& refD2, VectorType &q) const 159 { 160 // std::cout << a.X() << a.Y() << a.Z() << " - " << b.X() << b.Y() << b.Z() << 161 // " => " << vcg::Distance(a, b) << " < " << refD2 << "\n"; 162 Scalar d2 = vcg::SquaredDistance(a, b); 163 if (d2>refD2) 164 return false; 165 166 refD2 = d2; 167 q = a; 168 return true; 169 } 170 }; 171 172 class DummyObjectMarker {}; 173 174 protected: 175 const MeshType& mMesh; 176 const PointsType& mPoints; 177 vcg::Box3<Scalar> mAABB; 178 int mGradientHint; 179 int mHessianHint; 180 181 BallTree<Scalar>* mBallTree; 182 183 int mMaxNofProjectionIterations; 184 Scalar mFilterScale; 185 Scalar mAveragePointSpacing; 186 Scalar mProjectionAccuracy; 187 188 int mDomainMinNofNeighbors; 189 float mDomainRadiusScale; 190 float mDomainNormalScale; 191 192 // cached values: 193 mutable bool mCachedQueryPointIsOK; 194 mutable VectorType mCachedQueryPoint; 195 mutable Neighborhood<Scalar> mNeighborhood; 196 mutable std::vector<Scalar> mCachedWeights; 197 mutable std::vector<Scalar> mCachedWeightDerivatives; 198 mutable std::vector<VectorType> mCachedWeightGradients; 199 mutable std::vector<Scalar> mCachedWeightSecondDerivatives; 200 }; 201 202 } // namespace 203 204 #include "mlssurface.tpp" 205 206 #endif // MLSSURFACE_H 207