diff --git a/vcg/complex/algorithms/curve_on_manifold.h b/vcg/complex/algorithms/curve_on_manifold.h index 1bc35938..2c8e4abe 100644 --- a/vcg/complex/algorithms/curve_on_manifold.h +++ b/vcg/complex/algorithms/curve_on_manifold.h @@ -176,7 +176,7 @@ public: /** - * @brief MarkFauxEdgeWithPolyLine marks the edges of basemesh as non-faux when they coincide with the polyline ones * + * @brief TagFaceEdgeSelWithPolyLine selects edges of basemesh when they coincide with the polyline ones * * @param poly * @return true if all the edges of the polyline are snapped onto the mesh. * @@ -642,7 +642,7 @@ bool TagFaceEdgeSelWithPolyLine(MeshType &poly,bool markFlag=true) tri::UpdateTopology::TestVertexEdge(poly); tri::Allocator::CompactEveryVector(poly); tri::UpdateTopology::TestVertexEdge(poly); - printf("Simplify %5i -> %5i (total len %5.2f)\n",startEn,poly.en,hist.Sum()); +// printf("Simplify %5i -> %5i (total len %5.2f)\n",startEn,poly.en,hist.Sum()); } void EvaluateHausdorffDistance(MeshType &poly, Distribution &dist) @@ -1007,7 +1007,7 @@ bool TagFaceEdgeSelWithPolyLine(MeshType &poly,bool markFlag=true) } } // tri::Allocator::CompactEveryVector(poly); - printf("Refine %i -> %i\n",startEdgeSize,poly.en);fflush(stdout); +// printf("Refine %i -> %i\n",startEdgeSize,poly.en);fflush(stdout); } /** @@ -1068,7 +1068,7 @@ bool TagFaceEdgeSelWithPolyLine(MeshType &poly,bool markFlag=true) { tri::RequireCompactness(poly); tri::UpdateTopology::VertexEdge(poly); - printf("SmoothProject: Selected vert num %i\n",tri::UpdateSelection::VertexCount(poly)); +// printf("SmoothProject: Selected vert num %i\n",tri::UpdateSelection::VertexCount(poly)); assert(poly.en>0 && base.fn>0); for(int k=0;k::TestVertexEdge(poly); int dupVertNum = Clean::RemoveDuplicateVertex(poly); if(dupVertNum) { - printf("****REMOVED %i Duplicated\n",dupVertNum); +// printf("****REMOVED %i Duplicated\n",dupVertNum); tri::Allocator::CompactEveryVector(poly); tri::UpdateTopology::VertexEdge(poly); } diff --git a/vcg/complex/algorithms/smooth.h b/vcg/complex/algorithms/smooth.h index 5c1fd691..9eb4ced0 100644 --- a/vcg/complex/algorithms/smooth.h +++ b/vcg/complex/algorithms/smooth.h @@ -217,26 +217,28 @@ class Smooth //if we are applying to a tetrahedral mesh: ForEachTetra(m, [&](TetraType &t) { - for (int i = 0; i < 4; ++i) - if (!t.IsB(i)) + for (int i = 0; i < 6; ++i) + { + VertexPointer v0, v1, vo0, vo1; + v0 = t.V(Tetra::VofE(i, 0)); + v1 = t.V(Tetra::VofE(i, 1)); + + if (cotangentFlag) { - VertexPointer v0, v1, v2; - v0 = t.V(Tetra::VofF(i, 0)); - v1 = t.V(Tetra::VofF(i, 1)); - v2 = t.V(Tetra::VofF(i, 2)); + vo0 = t.V(Tetra::VofE(5 - i, 0)); + vo1 = t.V(Tetra::VofE(5 - i, 1)); - TD[v0].sum += v1->P() * weight; - TD[v0].sum += v2->P() * weight; - TD[v0].cnt += 2 * weight; + ScalarType angle = Tetra::DihedralAngle(t, 5 - i); + ScalarType length = vcg::Distance(vo0->P(), vo1->P()); - TD[v1].sum += v0->P() * weight; - TD[v1].sum += v2->P() * weight; - TD[v1].cnt += 2 * weight; - - TD[v2].sum += v0->P() * weight; - TD[v2].sum += v1->P() * weight; - TD[v2].cnt += 2 * weight; + weight = (length / 6.) * (tan(M_PI / 2. - angle)); } + + TD[v0].sum += v1->cP() * weight; + TD[v1].sum += v0->cP() * weight; + TD[v0].cnt += weight; + TD[v1].cnt += weight; + } }); ForEachTetra(m, [&](TetraType &t) { @@ -258,28 +260,28 @@ class Smooth } }); - ForEachTetra(m, [&](TetraType &t) { - for (int i = 0; i < 4; ++i) - if (t.IsB(i)) - { - VertexPointer v0, v1, v2; - v0 = t.V(Tetra::VofF(i, 0)); - v1 = t.V(Tetra::VofF(i, 1)); - v2 = t.V(Tetra::VofF(i, 2)); +// ForEachTetra(m, [&](TetraType &t) { +// for (int i = 0; i < 4; ++i) +// if (t.IsB(i)) +// { +// VertexPointer v0, v1, v2; +// v0 = t.V(Tetra::VofF(i, 0)); +// v1 = t.V(Tetra::VofF(i, 1)); +// v2 = t.V(Tetra::VofF(i, 2)); - TD[v0].sum += v1->P() * weight; - TD[v0].sum += v2->P() * weight; - TD[v0].cnt += 2 * weight; +// TD[v0].sum += v1->P(); +// TD[v0].sum += v2->P(); +// TD[v0].cnt += 2; - TD[v1].sum += v0->P() * weight; - TD[v1].sum += v2->P() * weight; - TD[v1].cnt += 2 * weight; +// TD[v1].sum += v0->P(); +// TD[v1].sum += v2->P(); +// TD[v1].cnt += 2; - TD[v2].sum += v0->P() * weight; - TD[v2].sum += v1->P() * weight; - TD[v2].cnt += 2 * weight; - } - }); +// TD[v2].sum += v0->P(); +// TD[v2].sum += v1->P(); +// TD[v2].cnt += 2; +// } +// }); FaceIterator fi; for (fi = m.face.begin(); fi != m.face.end(); ++fi) diff --git a/vcg/complex/algorithms/stat.h b/vcg/complex/algorithms/stat.h index 64b9a2aa..676314a0 100644 --- a/vcg/complex/algorithms/stat.h +++ b/vcg/complex/algorithms/stat.h @@ -31,6 +31,7 @@ #include #include #include +#include namespace vcg { @@ -264,6 +265,17 @@ public: return area/ScalarType(2.0); } + static ScalarType ComputePolyMeshArea(MeshType & m) + { + ScalarType area=0; + + for(FaceIterator fi = m.face.begin(); fi != m.face.end(); ++fi) + if(!(*fi).IsD()) + area += PolyArea(*fi); + + return area; + } + static ScalarType ComputeBorderLength(MeshType & m) { RequireFFAdjacency(m); diff --git a/vcg/complex/algorithms/tetra_implicit_smooth.h b/vcg/complex/algorithms/tetra_implicit_smooth.h new file mode 100644 index 00000000..49b40a10 --- /dev/null +++ b/vcg/complex/algorithms/tetra_implicit_smooth.h @@ -0,0 +1,495 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ +#ifndef __VCG_IMPLICIT_TETRA_SMOOTHER +#define __VCG_IMPLICIT_TETRA_SMOOTHER + +#include +#include +#include +#include + +#define PENALTY 10000 + +namespace vcg +{ + +template +class ImplicitTetraSmoother +{ + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::TetraType TetraType; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType::ScalarType ScalarType; + typedef typename Eigen::Matrix MatrixXm; + + public: + struct FaceConstraint + { + int numF; + std::vector BarycentricW; + CoordType TargetPos; + + FaceConstraint() + { + numF = -1; + } + + FaceConstraint(int _numF, + const std::vector &_BarycentricW, + const CoordType &_TargetPos) + { + numF = _numF; + BarycentricW = std::vector(_BarycentricW.begin(), _BarycentricW.end()); + TargetPos = _TargetPos; + } + }; + + struct Parameter + { + //the amount of smoothness, useful only if we set the mass matrix + ScalarType lambda; + //the use of mass matrix to keep the mesh close to its original position + //(weighted per area distributed on vertices) + bool useMassMatrix; + //this bool is used to fix the border vertices of the mesh or not + bool fixBorder; + //this bool is used to set if cotangent weight is used, this flag to false means uniform laplacian + bool useCotWeight; + //use this weight for the laplacian when the cotangent one is not used + ScalarType lapWeight; + //the set of fixed vertices + std::vector FixedV; + //the set of faces for barycentric constraints + std::vector ConstrainedF; + //the degree of laplacian + int degree; + //this is to say if we smooth the positions or the quality + bool SmoothQ; + + Parameter() + { + degree = 2; + lambda = 0.05; + useMassMatrix = true; + fixBorder = true; + useCotWeight = false; + lapWeight = 1; + SmoothQ = false; + } + }; + + private: + static void InitSparse(const std::vector> &Index, + const std::vector &Values, + const int m, + const int n, + Eigen::SparseMatrix &X) + { + assert(Index.size() == Values.size()); + + std::vector> IJV; + IJV.reserve(Index.size()); + + for (size_t i = 0; i < Index.size(); i++) + { + int row = Index[i].first; + int col = Index[i].second; + ScalarType val = Values[i]; + + assert(row < m); + assert(col < n); + + IJV.push_back(Eigen::Triplet(row, col, val)); + } + X.resize(m, n); + X.setFromTriplets(IJV.begin(), IJV.end()); + } + + static void CollectHardConstraints(MeshType &mesh, const Parameter &SParam, + std::vector> &IndexC, + std::vector &WeightC, + bool SmoothQ = false) + { + std::vector To_Fix; + + //collect fixed vert + if (SParam.fixBorder) + { + //add penalization constra + for (size_t i = 0; i < mesh.vert.size(); i++) + { + if (!mesh.vert[i].IsB()) + continue; + To_Fix.push_back(i); + } + } + //add additional fixed vertices constraint + To_Fix.insert(To_Fix.end(), SParam.FixedV.begin(), SParam.FixedV.end()); + + //sort and make them unique + std::sort(To_Fix.begin(), To_Fix.end()); + typename std::vector::iterator it = std::unique(To_Fix.begin(), To_Fix.end()); + To_Fix.resize(std::distance(To_Fix.begin(), it)); + + for (size_t i = 0; i < To_Fix.size(); i++) + { + if (!SmoothQ) + { + for (int j = 0; j < 3; j++) + { + int IndexV = (To_Fix[i] * 3) + j; + IndexC.push_back(std::pair(IndexV, IndexV)); + WeightC.push_back((ScalarType)PENALTY); + } + } + else + { + int IndexV = To_Fix[i]; + IndexC.push_back(std::pair(IndexV, IndexV)); + WeightC.push_back((ScalarType)PENALTY); + } + } + } + + static void CollectBarycentricConstraints(MeshType &mesh, + const Parameter &SParam, + std::vector> &IndexC, + std::vector &WeightC, + std::vector &IndexRhs, + std::vector &ValueRhs) + { + ScalarType penalty; + int baseIndex = mesh.vert.size(); + for (size_t i = 0; i < SParam.ConstrainedF.size(); i++) + { + //get the index of the current constraint + int IndexConstraint = baseIndex + i; + + //add one hard constraint + int FaceN = SParam.ConstrainedF[i].numF; + assert(FaceN >= 0); + assert(FaceN < (int)mesh.face.size()); + assert(mesh.face[FaceN].VN() == (int)SParam.ConstrainedF[i].BarycentricW.size()); + penalty = ScalarType(1) - SParam.lapWeight; + assert(penalty > ScalarType(0) && penalty < ScalarType(1)); + + //then add all the weights to impose the constraint + for (int j = 0; j < mesh.face[FaceN].VN(); j++) + { + //get the current weight + ScalarType currW = SParam.ConstrainedF[i].BarycentricW[j]; + + //get the index of the current vertex + int FaceVert = vcg::tri::Index(mesh, mesh.face[FaceN].V(j)); + + //then add the constraints componentwise + for (int k = 0; k < 3; k++) + { + //multiply times 3 per component + int IndexV = (FaceVert * 3) + k; + + //get the index of the current constraint + int ComponentConstraint = (IndexConstraint * 3) + k; + IndexC.push_back(std::pair(ComponentConstraint, IndexV)); + + WeightC.push_back(currW * penalty); + + IndexC.push_back(std::pair(IndexV, ComponentConstraint)); + WeightC.push_back(currW * penalty); + + //this to avoid the 1 on diagonal last entry of mass matrix + IndexC.push_back(std::pair(ComponentConstraint, ComponentConstraint)); + WeightC.push_back(-1); + } + } + + for (int j = 0; j < 3; j++) + { + //get the index of the current constraint + int ComponentConstraint = (IndexConstraint * 3) + j; + + //get per component value + ScalarType ComponentV = SParam.ConstrainedF[i].TargetPos.V(j); + + //add the diagonal value + IndexRhs.push_back(ComponentConstraint); + ValueRhs.push_back(ComponentV * penalty); + } + } + } + + static void MassMatrixEntry(MeshType &m, + std::vector> &index, + std::vector &entry, + bool vertexCoord = true) + { + tri::RequireCompactness(m); + + typename MeshType::template PerVertexAttributeHandle h = + tri::Allocator::template GetPerVertexAttribute(m, "volume"); + for (int i = 0; i < m.vn; ++i) + h[i] = 0; + + ForEachTetra(m, [&](TetraType &t) { + ScalarType v = Tetra::ComputeVolume(t); + for (int i = 0; i < 4; ++i) + h[tri::Index(m, t.V(i))] += v; + }); + + ScalarType maxV = 0; + for (int i = 0; i < m.vn; ++i) + maxV = max(maxV, h[i]); + + for (int i = 0; i < m.vn; ++i) + { + int currI = i; + index.push_back(std::pair(currI, currI)); + entry.push_back(h[i] / maxV); + } + + tri::Allocator::template DeletePerVertexAttribute(m, h); + } + + static ScalarType ComputeCotangentWeight(TetraType &t, const int i) + { + //i is the edge in the tetra + tetra::Pos pp(&t, Tetra::FofE(i, 0), i, Tetra::VofE(i, 0)); + tetra::Pos pt(&t, Tetra::FofE(i, 0), i, Tetra::VofE(i, 0)); + + ScalarType weight = 0; + + do + { + CoordType po0 = t.V(Tetra::VofE(5 - pt.E(), 0))->cP(); + CoordType po1 = t.V(Tetra::VofE(5 - pt.E(), 1))->cP(); + + ScalarType length = vcg::Distance(po0, po1); + ScalarType cot = std::tan((M_PI / 2.) - Tetra::DihedralAngle(*pt.T(), 5 - pt.E())); + + weight = (length / 6.) * cot; + pt.FlipT(); + pt.FlipF(); + } while (pp != pt); + + return weight; + } + + static void GetLaplacianEntry(MeshType &mesh, + TetraType &t, + std::vector> &index, + std::vector &entry, + bool cotangent, + ScalarType weight = 1, + bool vertexCoord = true) + { + // if (cotangent) + // vcg::tri::MeshAssert::OnlyT(mesh); + //iterate on edges + for (int i = 0; i < 6; ++i) + { + weight = 1;//ComputeCotangentWeight(t, i); + + int indexV0 = Index(mesh, t.V(Tetra::VofE(i, 0))); + int indexV1 = Index(mesh, t.V(Tetra::VofE(i, 1))); + + for (int j = 0; j < 3; j++) + { + //multiply by 3 and add the component + int currI0 = (indexV0 * 3) + j; + int currI1 = (indexV1 * 3) + j; + + index.push_back(std::pair(currI0, currI0)); + entry.push_back(weight); + index.push_back(std::pair(currI0, currI1)); + entry.push_back(-weight); + + index.push_back(std::pair(currI1, currI1)); + entry.push_back(weight); + index.push_back(std::pair(currI1, currI0)); + entry.push_back(-weight); + } + } + } + + static void GetLaplacianMatrix(MeshType &mesh, + std::vector> &index, + std::vector &entry, + bool cotangent, + ScalarType weight = 1, + bool vertexCoord = true) + { + //store the index and the scalar for the sparse matrix + ForEachTetra(mesh, [&](TetraType &t) { + GetLaplacianEntry(mesh, t, index, entry, cotangent, weight); + }); + } + + public: + static void Compute(MeshType &mesh, Parameter &SParam) + { + //calculate the size of the system + int matr_size = mesh.vert.size() + SParam.ConstrainedF.size(); + + //the laplacian and the mass matrix + Eigen::SparseMatrix L, M, B; + + //initialize the mass matrix + std::vector> IndexM; + std::vector ValuesM; + + //add the entries for mass matrix + if (SParam.useMassMatrix) + MassMatrixEntry(mesh, IndexM, ValuesM, !SParam.SmoothQ); + + //then add entries for lagrange mult due to barycentric constraints + for (size_t i = 0; i < SParam.ConstrainedF.size(); i++) + { + int baseIndex = (mesh.vert.size() + i) * 3; + + if (SParam.SmoothQ) + baseIndex = (mesh.vert.size() + i); + + if (SParam.SmoothQ) + { + IndexM.push_back(std::pair(baseIndex, baseIndex)); + ValuesM.push_back(1); + } + else + { + for (int j = 0; j < 3; j++) + { + IndexM.push_back(std::pair(baseIndex + j, baseIndex + j)); + ValuesM.push_back(1); + } + } + } + //add the hard constraints + CollectHardConstraints(mesh, SParam, IndexM, ValuesM, SParam.SmoothQ); + + //initialize sparse mass matrix + if (!SParam.SmoothQ) + InitSparse(IndexM, ValuesM, matr_size * 3, matr_size * 3, M); + else + InitSparse(IndexM, ValuesM, matr_size, matr_size, M); + + //initialize the barycentric matrix + std::vector> IndexB; + std::vector ValuesB; + + std::vector IndexRhs; + std::vector ValuesRhs; + + //then also collect hard constraints + if (!SParam.SmoothQ) + { + CollectBarycentricConstraints(mesh, SParam, IndexB, ValuesB, IndexRhs, ValuesRhs); + //initialize sparse constraint matrix + InitSparse(IndexB, ValuesB, matr_size * 3, matr_size * 3, B); + } + else + InitSparse(IndexB, ValuesB, matr_size, matr_size, B); + + //get the entries for laplacian matrix + std::vector> IndexL; + std::vector ValuesL; + GetLaplacianMatrix(mesh, IndexL, ValuesL, SParam.useCotWeight, SParam.lapWeight, !SParam.SmoothQ); + + //initialize sparse laplacian matrix + if (!SParam.SmoothQ) + InitSparse(IndexL, ValuesL, matr_size * 3, matr_size * 3, L); + else + InitSparse(IndexL, ValuesL, matr_size, matr_size, L); + + for (int i = 0; i < (SParam.degree - 1); i++) + L = L * L; + + //then solve the system + Eigen::SparseMatrix S = (M + B + SParam.lambda * L); + + //SimplicialLDLT + Eigen::SimplicialCholesky> solver(S); + assert(solver.info() == Eigen::Success); + + MatrixXm V; + if (!SParam.SmoothQ) + V = MatrixXm(matr_size * 3, 1); + else + V = MatrixXm(matr_size, 1); + + //set the first part of the matrix with vertex values + if (!SParam.SmoothQ) + { + for (size_t i = 0; i < mesh.vert.size(); i++) + { + int index = i * 3; + V(index, 0) = mesh.vert[i].P().X(); + V(index + 1, 0) = mesh.vert[i].P().Y(); + V(index + 2, 0) = mesh.vert[i].P().Z(); + } + } + else + { + for (size_t i = 0; i < mesh.vert.size(); i++) + { + int index = i; + V(index, 0) = mesh.vert[i].Q(); + } + } + + //then set the second part by considering RHS gien by barycentric constraint + for (size_t i = 0; i < IndexRhs.size(); i++) + { + int index = IndexRhs[i]; + ScalarType val = ValuesRhs[i]; + V(index, 0) = val; + } + + //solve the system + V = solver.solve(M * V).eval(); + + //then copy back values + if (!SParam.SmoothQ) + { + for (size_t i = 0; i < mesh.vert.size(); i++) + { + int index = i * 3; + mesh.vert[i].P().X() = V(index, 0); + mesh.vert[i].P().Y() = V(index + 1, 0); + mesh.vert[i].P().Z() = V(index + 2, 0); + } + } + else + { + for (size_t i = 0; i < mesh.vert.size(); i++) + { + int index = i; + mesh.vert[i].Q() = V(index, 0); + } + } + } +}; + +} //end namespace vcg + +#endif diff --git a/vcg/complex/algorithms/update/flag.h b/vcg/complex/algorithms/update/flag.h index 1ff86976..c4d74ac4 100644 --- a/vcg/complex/algorithms/update/flag.h +++ b/vcg/complex/algorithms/update/flag.h @@ -40,466 +40,486 @@ class UpdateFlags { public: - typedef UpdateMeshType MeshType; - typedef typename MeshType::ScalarType ScalarType; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::VertexPointer VertexPointer; - typedef typename MeshType::VertexIterator VertexIterator; - typedef typename MeshType::EdgeType EdgeType; - typedef typename MeshType::EdgePointer EdgePointer; - typedef typename MeshType::EdgeIterator EdgeIterator; - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::FacePointer FacePointer; - typedef typename MeshType::FaceIterator FaceIterator; - typedef typename MeshType::TetraType TetraType; - typedef typename MeshType::TetraPointer TetraPointer; - typedef typename MeshType::TetraIterator TetraIterator; + typedef UpdateMeshType MeshType; + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::EdgeType EdgeType; + typedef typename MeshType::EdgePointer EdgePointer; + typedef typename MeshType::EdgeIterator EdgeIterator; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::TetraType TetraType; + typedef typename MeshType::TetraPointer TetraPointer; + typedef typename MeshType::TetraIterator TetraIterator; - /// \brief Reset all the mesh flags (vertexes edge faces) setting everithing to zero (the default value for flags) + /// \brief Reset all the mesh flags (vertexes edge faces) setting everithing to zero (the default value for flags) - static void Clear(MeshType &m) - { - if(HasPerVertexFlags(m) ) - for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) - (*vi).Flags() = 0; - if(HasPerEdgeFlags(m) ) - for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) - (*ei).Flags() = 0; - if(HasPerFaceFlags(m) ) - for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) - (*fi).Flags() = 0; - if(HasPerTetraFlags(m) ) - for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) - (*ti).Flags() = 0; - } + static void Clear(MeshType &m) + { + if(HasPerVertexFlags(m) ) + for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) + (*vi).Flags() = 0; + if(HasPerEdgeFlags(m) ) + for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) + (*ei).Flags() = 0; + if(HasPerFaceFlags(m) ) + for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) + (*fi).Flags() = 0; + if(HasPerTetraFlags(m) ) + for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) + (*ti).Flags() = 0; + } - static void VertexClear(MeshType &m, unsigned int FlagMask = 0xffffffff) - { - RequirePerVertexFlags(m); - int andMask = ~FlagMask; - for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) - if(!(*vi).IsD()) (*vi).Flags() &= andMask ; - } + static void VertexClear(MeshType &m, unsigned int FlagMask = 0xffffffff) + { + RequirePerVertexFlags(m); + int andMask = ~FlagMask; + for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) + if(!(*vi).IsD()) (*vi).Flags() &= andMask ; + } - static void EdgeClear(MeshType &m, unsigned int FlagMask = 0xffffffff) - { - RequirePerEdgeFlags(m); - int andMask = ~FlagMask; - for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) - if(!(*ei).IsD()) (*ei).Flags() &= andMask ; - } + static void EdgeClear(MeshType &m, unsigned int FlagMask = 0xffffffff) + { + RequirePerEdgeFlags(m); + int andMask = ~FlagMask; + for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) + if(!(*ei).IsD()) (*ei).Flags() &= andMask ; + } - static void FaceClear(MeshType &m, unsigned int FlagMask = 0xffffffff) - { - RequirePerFaceFlags(m); - int andMask = ~FlagMask; - for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) - if(!(*fi).IsD()) (*fi).Flags() &= andMask ; - } + static void FaceClear(MeshType &m, unsigned int FlagMask = 0xffffffff) + { + RequirePerFaceFlags(m); + int andMask = ~FlagMask; + for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) + if(!(*fi).IsD()) (*fi).Flags() &= andMask ; + } - static void TetraClear(MeshType &m, unsigned int FlagMask = 0xffffffff) - { - RequirePerTetraFlags(m); - int andMask = ~FlagMask; - for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) - if(!(*ti).IsD()) (*ti).Flags() &= andMask ; - } + static void TetraClear(MeshType &m, unsigned int FlagMask = 0xffffffff) + { + RequirePerTetraFlags(m); + int andMask = ~FlagMask; + for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) + if(!(*ti).IsD()) (*ti).Flags() &= andMask ; + } - static void VertexSet(MeshType &m, unsigned int FlagMask) - { - RequirePerVertexFlags(m); - for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) - if(!(*vi).IsD()) (*vi).Flags() |= FlagMask ; - } + static void VertexSet(MeshType &m, unsigned int FlagMask) + { + RequirePerVertexFlags(m); + for(VertexIterator vi=m.vert.begin(); vi!=m.vert.end(); ++vi) + if(!(*vi).IsD()) (*vi).Flags() |= FlagMask ; + } - static void EdgeSet(MeshType &m, unsigned int FlagMask) - { - RequirePerEdgeFlags(m); - for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) - if(!(*ei).IsD()) (*ei).Flags() |= FlagMask ; - } + static void EdgeSet(MeshType &m, unsigned int FlagMask) + { + RequirePerEdgeFlags(m); + for(EdgeIterator ei=m.edge.begin(); ei!=m.edge.end(); ++ei) + if(!(*ei).IsD()) (*ei).Flags() |= FlagMask ; + } - static void FaceSet(MeshType &m, unsigned int FlagMask) - { - RequirePerFaceFlags(m); - for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) - if(!(*fi).IsD()) (*fi).Flags() |= FlagMask ; - } + static void FaceSet(MeshType &m, unsigned int FlagMask) + { + RequirePerFaceFlags(m); + for(FaceIterator fi=m.face.begin(); fi!=m.face.end(); ++fi) + if(!(*fi).IsD()) (*fi).Flags() |= FlagMask ; + } - static void TetraSet(MeshType &m, unsigned int FlagMask) - { - RequirePerTetraFlags(m); - for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) - if(!(*ti).IsD()) (*ti).Flags() |= FlagMask ; - } + static void TetraSet(MeshType &m, unsigned int FlagMask) + { + RequirePerTetraFlags(m); + for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) + if(!(*ti).IsD()) (*ti).Flags() |= FlagMask ; + } - static void VertexClearV(MeshType &m) { VertexClear(m,VertexType::VISITED);} - static void VertexClearS(MeshType &m) { VertexClear(m,VertexType::SELECTED);} - static void VertexClearB(MeshType &m) { VertexClear(m,VertexType::BORDER);} - static void EdgeClearV(MeshType &m) { EdgeClear(m,EdgeType::VISITED);} - static void FaceClearV(MeshType &m) { FaceClear(m,FaceType::VISITED);} - static void FaceClearB(MeshType &m) { FaceClear(m,FaceType::BORDER012);} - static void FaceClearS(MeshType &m) {FaceClear(m,FaceType::SELECTED);} - static void FaceClearF(MeshType &m) { FaceClear(m,FaceType::FAUX012);} - static void FaceClearFaceEdgeS(MeshType &m) { FaceClear(m,FaceType::FACEEDGESEL012 ); } + static void VertexClearV(MeshType &m) { VertexClear(m,VertexType::VISITED);} + static void VertexClearS(MeshType &m) { VertexClear(m,VertexType::SELECTED);} + static void VertexClearB(MeshType &m) { VertexClear(m,VertexType::BORDER);} + static void EdgeClearV(MeshType &m) { EdgeClear(m,EdgeType::VISITED);} + static void FaceClearV(MeshType &m) { FaceClear(m,FaceType::VISITED);} + static void FaceClearB(MeshType &m) { FaceClear(m,FaceType::BORDER012);} + static void FaceClearS(MeshType &m) {FaceClear(m,FaceType::SELECTED);} + static void FaceClearF(MeshType &m) { FaceClear(m,FaceType::FAUX012);} + static void FaceClearFaceEdgeS(MeshType &m) { FaceClear(m,FaceType::FACEEDGESEL012 ); } - static void EdgeSetV(MeshType &m) { EdgeSet(m,EdgeType::VISITED);} - static void VertexSetV(MeshType &m) { VertexSet(m,VertexType::VISITED);} - static void VertexSetS(MeshType &m) { VertexSet(m,VertexType::SELECTED);} - static void VertexSetB(MeshType &m) { VertexSet(m,VertexType::BORDER);} - static void FaceSetV(MeshType &m) { FaceSet(m,FaceType::VISITED);} - static void FaceSetB(MeshType &m) { FaceSet(m,FaceType::BORDER);} - static void FaceSetF(MeshType &m) { FaceSet(m,FaceType::FAUX012);} - static void TetraClearV(MeshType &m) { TetraClear(m, TetraType::VISITED); } - static void TetraClearS(MeshType &m) { TetraClear(m, TetraType::SELECTED); } - static void TetraClearB(MeshType &m) { TetraClear(m, TetraType::BORDER0123); } - static void TetraSetV(MeshType &m) { TetraSet(m, TetraType::VISITED); } - static void TetraSetS(MeshType &m) { TetraSet(m, TetraType::SELECTED); } - static void TetraSetB(MeshType &m) { TetraSet(m, TetraType::BORDER0123); } - /// \brief Compute the border flags for the faces using the Face-Face Topology. - /** + static void EdgeSetV(MeshType &m) { EdgeSet(m,EdgeType::VISITED);} + static void VertexSetV(MeshType &m) { VertexSet(m,VertexType::VISITED);} + static void VertexSetS(MeshType &m) { VertexSet(m,VertexType::SELECTED);} + static void VertexSetB(MeshType &m) { VertexSet(m,VertexType::BORDER);} + static void FaceSetV(MeshType &m) { FaceSet(m,FaceType::VISITED);} + static void FaceSetB(MeshType &m) { FaceSet(m,FaceType::BORDER);} + static void FaceSetF(MeshType &m) { FaceSet(m,FaceType::FAUX012);} + static void TetraClearV(MeshType &m) { TetraClear(m, TetraType::VISITED); } + static void TetraClearS(MeshType &m) { TetraClear(m, TetraType::SELECTED); } + static void TetraClearB(MeshType &m) { TetraClear(m, TetraType::BORDER0123); } + static void TetraSetV(MeshType &m) { TetraSet(m, TetraType::VISITED); } + static void TetraSetS(MeshType &m) { TetraSet(m, TetraType::SELECTED); } + static void TetraSetB(MeshType &m) { TetraSet(m, TetraType::BORDER0123); } + /// \brief Compute the border flags for the faces using the Face-Face Topology. + /** \warning Obviously it assumes that the topology has been correctly computed (see: UpdateTopology::FaceFace ) */ - static void FaceBorderFromFF(MeshType &m) - { - RequirePerFaceFlags(m); - RequireFFAdjacency(m); + static void FaceBorderFromFF(MeshType &m) + { + RequirePerFaceFlags(m); + RequireFFAdjacency(m); - for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD()) - for(int j=0;jVN();++j) - { - if(face::IsBorder(*fi,j)) (*fi).SetB(j); - else (*fi).ClearB(j); - } - } + for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD()) + for(int j=0;jVN();++j) + { + if(face::IsBorder(*fi,j)) (*fi).SetB(j); + else (*fi).ClearB(j); + } + } - /// \brief Compute the border flags for the tetras using the Tetra-Tetra Topology. - /** + /// \brief Compute the border flags for the tetras using the Tetra-Tetra Topology. + /** \warning Obviously it assumes that the topology has been correctly computed (see: UpdateTopology::FaceFace ) */ - static void TetraBorderFromTT(MeshType &m) - { - RequirePerTetraFlags(m); - RequireTTAdjacency(m); - - for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) - if(!(*ti).IsD()) - for(int j = 0; j < 4; ++j) - { - if (tetrahedron::IsBorder(*ti,j)) (*ti).SetB(j); - else (*ti).ClearB(j); - } - } - - static void FaceBorderFromVF(MeshType &m) - { - RequirePerFaceFlags(m); - RequireVFAdjacency(m); - - FaceClearB(m); - int visitedBit=VertexType::NewBitFlag(); - - // Calcolo dei bordi - // per ogni vertice vi si cercano i vertici adiacenti che sono toccati da una faccia sola - // (o meglio da un numero dispari di facce) - - const int BORDERFLAG[3]={FaceType::BORDER0, FaceType::BORDER1, FaceType::BORDER2}; - - for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD()) - { - for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) - { - vfi.f->V1(vfi.z)->ClearUserBit(visitedBit); - vfi.f->V2(vfi.z)->ClearUserBit(visitedBit); - } - for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) - { - if(vfi.f->V1(vfi.z)->IsUserBit(visitedBit)) vfi.f->V1(vfi.z)->ClearUserBit(visitedBit); - else vfi.f->V1(vfi.z)->SetUserBit(visitedBit); - if(vfi.f->V2(vfi.z)->IsUserBit(visitedBit)) vfi.f->V2(vfi.z)->ClearUserBit(visitedBit); - else vfi.f->V2(vfi.z)->SetUserBit(visitedBit); - } - for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) - { - if(vfi.f->V(vfi.z)< vfi.f->V1(vfi.z) && vfi.f->V1(vfi.z)->IsUserBit(visitedBit)) - vfi.f->Flags() |= BORDERFLAG[vfi.z]; - if(vfi.f->V(vfi.z)< vfi.f->V2(vfi.z) && vfi.f->V2(vfi.z)->IsUserBit(visitedBit)) - vfi.f->Flags() |= BORDERFLAG[(vfi.z+2)%3]; - } - } - VertexType::DeleteBitFlag(visitedBit); - } - - - class EdgeSorter - { - public: - - VertexPointer v[2]; // Puntatore ai due vertici (Ordinati) - FacePointer f; // Puntatore alla faccia generatrice - int z; // Indice dell'edge nella faccia - - EdgeSorter() {} // Nothing to do - - - void Set( const FacePointer pf, const int nz ) + static void TetraBorderFromTT(MeshType &m) { - assert(pf!=0); - assert(nz>=0); - assert(nz<3); + RequirePerTetraFlags(m); + RequireTTAdjacency(m); - v[0] = pf->V(nz); - v[1] = pf->V((nz+1)%3); - assert(v[0] != v[1]); - - if( v[0] > v[1] ) std::swap(v[0],v[1]); - f = pf; - z = nz; + for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) + if(!(*ti).IsD()) + for(int j = 0; j < 4; ++j) + { + if (tetrahedron::IsBorder(*ti,j)) (*ti).SetB(j); + else (*ti).ClearB(j); + } } - inline bool operator < ( const EdgeSorter & pe ) const { - if( v[0]pe.v[0] ) return false; - else return v[1] < pe.v[1]; - } - - inline bool operator == ( const EdgeSorter & pe ) const + static void VertexBorderFromTT(MeshType &m) { - return v[0]==pe.v[0] && v[1]==pe.v[1]; + RequirePerVertexFlags(m); + RequireTTAdjacency(m); + + VertexClearB(m); + + for(TetraIterator ti=m.tetra.begin(); ti!=m.tetra.end(); ++ti) + if(!(*ti).IsD()) + for(int j = 0; j < 4; ++j) + { + if (tetrahedron::IsBorder(*ti,j)) + { + for (int i = 0; i < 3; ++i) + ti->V(Tetra::VofF(j, i))->SetB(); + } + } } - inline bool operator != ( const EdgeSorter & pe ) const + + + static void FaceBorderFromVF(MeshType &m) { - return v[0]!=pe.v[0] || v[1]!=pe.v[1]; + RequirePerFaceFlags(m); + RequireVFAdjacency(m); + + FaceClearB(m); + int visitedBit=VertexType::NewBitFlag(); + + // Calcolo dei bordi + // per ogni vertice vi si cercano i vertici adiacenti che sono toccati da una faccia sola + // (o meglio da un numero dispari di facce) + + const int BORDERFLAG[3]={FaceType::BORDER0, FaceType::BORDER1, FaceType::BORDER2}; + + for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD()) + { + for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) + { + vfi.f->V1(vfi.z)->ClearUserBit(visitedBit); + vfi.f->V2(vfi.z)->ClearUserBit(visitedBit); + } + for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) + { + if(vfi.f->V1(vfi.z)->IsUserBit(visitedBit)) vfi.f->V1(vfi.z)->ClearUserBit(visitedBit); + else vfi.f->V1(vfi.z)->SetUserBit(visitedBit); + if(vfi.f->V2(vfi.z)->IsUserBit(visitedBit)) vfi.f->V2(vfi.z)->ClearUserBit(visitedBit); + else vfi.f->V2(vfi.z)->SetUserBit(visitedBit); + } + for(face::VFIterator vfi(&*vi) ; !vfi.End(); ++vfi ) + { + if(vfi.f->V(vfi.z)< vfi.f->V1(vfi.z) && vfi.f->V1(vfi.z)->IsUserBit(visitedBit)) + vfi.f->Flags() |= BORDERFLAG[vfi.z]; + if(vfi.f->V(vfi.z)< vfi.f->V2(vfi.z) && vfi.f->V2(vfi.z)->IsUserBit(visitedBit)) + vfi.f->Flags() |= BORDERFLAG[(vfi.z+2)%3]; + } + } + VertexType::DeleteBitFlag(visitedBit); } - }; + + class EdgeSorter + { + public: + + VertexPointer v[2]; // Puntatore ai due vertici (Ordinati) + FacePointer f; // Puntatore alla faccia generatrice + int z; // Indice dell'edge nella faccia + + EdgeSorter() {} // Nothing to do - // versione minimale che non calcola i complex flag. - static void VertexBorderFromNone(MeshType &m) - { - RequirePerVertexFlags(m); - - std::vector e; - typename UpdateMeshType::FaceIterator pf; - typename std::vector::iterator p; - - if( m.fn == 0 ) - return; - - e.resize(m.fn*3); // Alloco il vettore ausiliario - p = e.begin(); - for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce - if( ! (*pf).IsD() ) - for(int j=0;j<3;++j) + void Set( const FacePointer pf, const int nz ) { - (*p).Set(&(*pf),j); - (*pf).ClearB(j); - ++p; - } - assert(p==e.end()); - sort(e.begin(), e.end()); // Lo ordino per vertici + assert(pf!=0); + assert(nz>=0); + assert(nz<3); - typename std::vector::iterator pe,ps; - for(ps = e.begin(), pe = e.begin(); pe < e.end(); ++pe) // Scansione vettore ausiliario + v[0] = pf->V(nz); + v[1] = pf->V((nz+1)%3); + assert(v[0] != v[1]); + + if( v[0] > v[1] ) std::swap(v[0],v[1]); + f = pf; + z = nz; + } + + inline bool operator < ( const EdgeSorter & pe ) const { + if( v[0]pe.v[0] ) return false; + else return v[1] < pe.v[1]; + } + + inline bool operator == ( const EdgeSorter & pe ) const + { + return v[0]==pe.v[0] && v[1]==pe.v[1]; + } + inline bool operator != ( const EdgeSorter & pe ) const + { + return v[0]!=pe.v[0] || v[1]!=pe.v[1]; + } + + }; + + + // versione minimale che non calcola i complex flag. + static void VertexBorderFromNone(MeshType &m) { - if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali - { - if(pe-ps==1) { - ps->v[0]->SetB(); - ps->v[1]->SetB(); - }/* else + RequirePerVertexFlags(m); + + std::vector e; + typename UpdateMeshType::FaceIterator pf; + typename std::vector::iterator p; + + if( m.fn == 0 ) + return; + + e.resize(m.fn*3); // Alloco il vettore ausiliario + p = e.begin(); + for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce + if( ! (*pf).IsD() ) + for(int j=0;j<3;++j) + { + (*p).Set(&(*pf),j); + (*pf).ClearB(j); + ++p; + } + assert(p==e.end()); + sort(e.begin(), e.end()); // Lo ordino per vertici + + typename std::vector::iterator pe,ps; + for(ps = e.begin(), pe = e.begin(); pe < e.end(); ++pe) // Scansione vettore ausiliario + { + if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali + { + if(pe-ps==1) { + ps->v[0]->SetB(); + ps->v[1]->SetB(); + }/* else if(pe-ps!=2) { // not twomanyfold! for(;ps!=pe;++ps) { ps->v[0]->SetB(); // Si settano border anche i complex. ps->v[1]->SetB(); } }*/ - ps = pe; - } - } - } - - /// Computes per-face border flags without requiring any kind of topology - /// It has a O(fn log fn) complexity. - static void FaceBorderFromNone(MeshType &m) - { - RequirePerFaceFlags(m); - - std::vector e; - typename UpdateMeshType::FaceIterator pf; - typename std::vector::iterator p; - - for(VertexIterator v=m.vert.begin();v!=m.vert.end();++v) - (*v).ClearB(); - - if( m.fn == 0 ) - return; - - FaceIterator fi; - int n_edges = 0; - for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD()) n_edges+=(*fi).VN(); - e.resize(n_edges); - - p = e.begin(); - for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce - if( ! (*pf).IsD() ) - for(int j=0;j<(*pf).VN();++j) - { - (*p).Set(&(*pf),j); - (*pf).ClearB(j); - ++p; + ps = pe; + } } - assert(p==e.end()); - sort(e.begin(), e.end()); // Lo ordino per vertici + } - typename std::vector::iterator pe,ps; - ps = e.begin();pe=e.begin(); - do + /// Computes per-face border flags without requiring any kind of topology + /// It has a O(fn log fn) complexity. + static void FaceBorderFromNone(MeshType &m) { - if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali - { - if(pe-ps==1) { - ps->f->SetB(ps->z); - } /*else + RequirePerFaceFlags(m); + + std::vector e; + typename UpdateMeshType::FaceIterator pf; + typename std::vector::iterator p; + + for(VertexIterator v=m.vert.begin();v!=m.vert.end();++v) + (*v).ClearB(); + + if( m.fn == 0 ) + return; + + FaceIterator fi; + int n_edges = 0; + for(fi = m.face.begin(); fi != m.face.end(); ++fi) if(! (*fi).IsD()) n_edges+=(*fi).VN(); + e.resize(n_edges); + + p = e.begin(); + for(pf=m.face.begin();pf!=m.face.end();++pf) // Lo riempio con i dati delle facce + if( ! (*pf).IsD() ) + for(int j=0;j<(*pf).VN();++j) + { + (*p).Set(&(*pf),j); + (*pf).ClearB(j); + ++p; + } + assert(p==e.end()); + sort(e.begin(), e.end()); // Lo ordino per vertici + + typename std::vector::iterator pe,ps; + ps = e.begin();pe=e.begin(); + do + { + if( pe==e.end() || *pe != *ps ) // Trovo blocco di edge uguali + { + if(pe-ps==1) { + ps->f->SetB(ps->z); + } /*else if(pe-ps!=2) { // Caso complex!! for(;ps!=pe;++ps) ps->f->SetB(ps->z); // Si settano border anche i complex. }*/ - ps = pe; - } - if(pe==e.end()) break; - ++pe; - } while(true); - // TRACE("found %i border (%i complex) on %i edges\n",nborder,ncomplex,ne); - } - - /// Compute the PerVertex Border flag deriving it from the face-face adjacency - static void VertexBorderFromFaceAdj(MeshType &m) -{ - RequirePerFaceFlags(m); - RequirePerVertexFlags(m); - RequireFFAdjacency(m); - // MeshAssert::FFAdjacencyIsInitialized(m); - - VertexClearB(m); - for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - { - - for(int z=0;z<(*fi).VN();++z) - if( face::IsBorder(*fi,z)) - { - (*fi).V0(z)->SetB(); - (*fi).V1(z)->SetB(); - } - } - } - - /// Compute the PerVertex Border flag deriving it from the border flag of faces - static void VertexBorderFromFaceBorder(MeshType &m) - { - RequirePerFaceFlags(m); - RequirePerVertexFlags(m); - VertexClearB(m); - for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - { - for(int z=0;z<(*fi).VN();++z) - if( (*fi).IsB(z) ) - { - (*fi).V(z)->SetB(); - (*fi).V((*fi).Next(z))->SetB(); - } - } - } - - /// Compute the PerVertex Border flag deriving it from the Edge-Edge adjacency (made for edgemeshes) - static void VertexBorderFromEdgeAdj(MeshType &m) - { - RequirePerVertexFlags(m); - RequireEEAdjacency(m); - - VertexClearB(m); - for (EdgeIterator ei=m.edge.begin();ei!=m.edge.end();++ei) - if (!ei->IsD()) - { - for (int z=0; z<2; ++z) - if (edge::IsEdgeBorder(*ei, z)) - { - ei->V(z)->SetB(); - } - } - } - - /// \brief Marks feature edges according to two signed dihedral angles. - /// Actually it uses the face_edge selection bit on faces, - /// we select the edges where the signed dihedral angle between the normal of two incident faces , - /// is outside the two given thresholds. - /// In this way all the edges that are almost planar are marked as non selected (e.g. edges to be ignored) - /// Note that it uses the signed dihedral angle convention (negative for concave edges and positive for convex ones); - /// - /// Optionally it can also mark as feature edges also the boundary edges. - /// - static void FaceEdgeSelSignedCrease(MeshType &m, float AngleRadNeg, float AngleRadPos, bool MarkBorderFlag = false ) - { - RequirePerFaceFlags(m); - RequireFFAdjacency(m); - //initially Nothing is faux (e.g all crease) - FaceClearFaceEdgeS(m); - // Then mark faux only if the signed angle is the range. - for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - for(int z=0;z<(*fi).VN();++z) - { - if(!face::IsBorder(*fi,z) ) - { - ScalarType angle = DihedralAngleRad(*fi,z); - if(angleAngleRadPos) - (*fi).SetFaceEdgeS(z); - } - else - { - if(MarkBorderFlag) (*fi).SetFaceEdgeS(z); - } - } + ps = pe; + } + if(pe==e.end()) break; + ++pe; + } while(true); + // TRACE("found %i border (%i complex) on %i edges\n",nborder,ncomplex,ne); } - } - /// \brief Selects feature edges according to Face adjacency. - /// - static void FaceEdgeSelBorder(MeshType &m) - { - RequirePerFaceFlags(m); - RequireFFAdjacency(m); - //initially Nothing is selected - FaceClearFaceEdgeS(m); - for (FaceIterator fi=m.face.begin(); fi!=m.face.end();++fi) - if (!fi->IsD()) - { - for (int z=0; z<(*fi).VN(); ++z) - { - if (face::IsBorder(*fi,z)) - fi->SetFaceEdgeS(z); - } - } - } + /// Compute the PerVertex Border flag deriving it from the face-face adjacency + static void VertexBorderFromFaceAdj(MeshType &m) + { + RequirePerFaceFlags(m); + RequirePerVertexFlags(m); + RequireFFAdjacency(m); + // MeshAssert::FFAdjacencyIsInitialized(m); + + VertexClearB(m); + for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + { + + for(int z=0;z<(*fi).VN();++z) + if( face::IsBorder(*fi,z)) + { + (*fi).V0(z)->SetB(); + (*fi).V1(z)->SetB(); + } + } + } + + /// Compute the PerVertex Border flag deriving it from the border flag of faces + static void VertexBorderFromFaceBorder(MeshType &m) + { + RequirePerFaceFlags(m); + RequirePerVertexFlags(m); + VertexClearB(m); + for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + { + for(int z=0;z<(*fi).VN();++z) + if( (*fi).IsB(z) ) + { + (*fi).V(z)->SetB(); + (*fi).V((*fi).Next(z))->SetB(); + } + } + } + + /// Compute the PerVertex Border flag deriving it from the Edge-Edge adjacency (made for edgemeshes) + static void VertexBorderFromEdgeAdj(MeshType &m) + { + RequirePerVertexFlags(m); + RequireEEAdjacency(m); + + VertexClearB(m); + for (EdgeIterator ei=m.edge.begin();ei!=m.edge.end();++ei) + if (!ei->IsD()) + { + for (int z=0; z<2; ++z) + if (edge::IsEdgeBorder(*ei, z)) + { + ei->V(z)->SetB(); + } + } + } + + /// \brief Marks feature edges according to two signed dihedral angles. + /// Actually it uses the face_edge selection bit on faces, + /// we select the edges where the signed dihedral angle between the normal of two incident faces , + /// is outside the two given thresholds. + /// In this way all the edges that are almost planar are marked as non selected (e.g. edges to be ignored) + /// Note that it uses the signed dihedral angle convention (negative for concave edges and positive for convex ones); + /// + /// Optionally it can also mark as feature edges also the boundary edges. + /// + static void FaceEdgeSelSignedCrease(MeshType &m, float AngleRadNeg, float AngleRadPos, bool MarkBorderFlag = false ) + { + RequirePerFaceFlags(m); + RequireFFAdjacency(m); + //initially Nothing is faux (e.g all crease) + FaceClearFaceEdgeS(m); + // Then mark faux only if the signed angle is the range. + for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) + { + for(int z=0;z<(*fi).VN();++z) + { + if(!face::IsBorder(*fi,z) ) + { + ScalarType angle = DihedralAngleRad(*fi,z); + if(angleAngleRadPos) + (*fi).SetFaceEdgeS(z); + } + else + { + if(MarkBorderFlag) (*fi).SetFaceEdgeS(z); + } + } + } + } + + /// \brief Selects feature edges according to Face adjacency. + /// + static void FaceEdgeSelBorder(MeshType &m) + { + RequirePerFaceFlags(m); + RequireFFAdjacency(m); + //initially Nothing is selected + FaceClearFaceEdgeS(m); + for (FaceIterator fi=m.face.begin(); fi!=m.face.end();++fi) + if (!fi->IsD()) + { + for (int z=0; z<(*fi).VN(); ++z) + { + if (face::IsBorder(*fi,z)) + fi->SetFaceEdgeS(z); + } + } + } + + /// \brief Marks feature edges according to a given angle + /// Actually it uses the face_edge selection bit on faces, + /// we select the edges where the dihedral angle between the normal of two incident faces is larger than , + /// the given thresholds. + /// In this way all the near planar edges are marked remains not selected (e.g. edges to be ignored) + static void FaceEdgeSelCrease(MeshType &m,float AngleRad) + { + FaceEdgeSelSignedCrease(m,-AngleRad,AngleRad); + } - /// \brief Marks feature edges according to a given angle - /// Actually it uses the face_edge selection bit on faces, - /// we select the edges where the dihedral angle between the normal of two incident faces is larger than , - /// the given thresholds. - /// In this way all the near planar edges are marked remains not selected (e.g. edges to be ignored) - static void FaceEdgeSelCrease(MeshType &m,float AngleRad) - { - FaceEdgeSelSignedCrease(m,-AngleRad,AngleRad); - } - }; // end class } // End namespace tri diff --git a/vcg/complex/base.h b/vcg/complex/base.h index 5032399a..7014f12e 100644 --- a/vcg/complex/base.h +++ b/vcg/complex/base.h @@ -413,6 +413,7 @@ public: vert.clear(); face.clear(); edge.clear(); + tetra.clear(); // textures.clear(); // normalmaps.clear(); vn = 0; diff --git a/vcg/simplex/tetrahedron/pos.h b/vcg/simplex/tetrahedron/pos.h index 216c811e..59b35ae7 100644 --- a/vcg/simplex/tetrahedron/pos.h +++ b/vcg/simplex/tetrahedron/pos.h @@ -165,7 +165,7 @@ public: return _e; } - /// Return the index of face as seen from the tetrahedron + /// Return the index of edge as seen from the tetrahedron inline const char & E() const { return _e; @@ -269,7 +269,7 @@ public: //get the current vertex VertexType *vcurr=T()->V(V()); - //get new tetrahedron according to faceto face topology + //get new tetrahedron according to tetra to tetra topology TetraType *nt=T()->TTp(F()); char nfa=T()->TTi(F()); if (nfa!=-1) diff --git a/vcg/space/polygon3.h b/vcg/space/polygon3.h index 2e02e6fc..069830c4 100644 --- a/vcg/space/polygon3.h +++ b/vcg/space/polygon3.h @@ -126,6 +126,9 @@ typename PolygonType::ScalarType PolyArea(const PolygonType &F) typedef typename PolygonType::CoordType CoordType; typedef typename PolygonType::ScalarType ScalarType; + if (F.VN() == 3) + return vcg::DoubleArea(F) / 2; + CoordType bary=PolyBarycenter(F); ScalarType Area=0; for (size_t i=0;i<(size_t)F.VN();i++) @@ -535,12 +538,12 @@ typename PolygonType::ScalarType PolyAspectRatio(const PolygonType &F, template typename PolygonType::ScalarType PolygonPointDistance(const PolygonType &F, const vcg::Point3 &pos, - vcg::Point3 &ClosestP) + vcg::Point3 &ClosestP, + typename PolygonType::ScalarType minD = std::numeric_limits::max()) { typedef typename PolygonType::ScalarType ScalarType; typedef typename PolygonType::CoordType CoordType; - ScalarType minD=std::numeric_limits::max(); CoordType bary=vcg::PolyBarycenter(F); for (size_t j=0;j +template class GlTetramesh:public GLW{ public: - typedef typename CONT_TETRA::value_type TetraType; - typedef typename TetraType::VertexType VertexType; + typedef typename MeshType::TetraType TetraType; + typedef typename TetraType::VertexType VertexType; typedef typename VertexType::ScalarType ScalarType; - typedef typename VertexType::CoordType Point3x; + typedef typename VertexType::CoordType CoordType; //subclass for clipping planes class ClipPlane { private: - Point3x D; - Point3x D0; - GLdouble eqn[4]; - vcg::Matrix44 TR; - Point3x pp0; - Point3x pp1; - Point3x pp2; - Point3x pp3; + CoordType D, D0; + ScalarType dist; + + GLdouble eqn[4]; + public: - bool active; - Point3x P; - ClipPlane (){active=false;} ~ClipPlane (){} - ClipPlane(Point3x p0, Point3x p1,Point3x p2) + ClipPlane(CoordType & p0, CoordType & p1, CoordType & p2) { - Point3x N=((p1-p0)^(p2-p0)).Normalize(); - N.Normalize(); - D=N; - D0=D; - P=(p0+p1+p2)/3.f; + CoordType N = ((p1-p0)^(p2-p0)).Normalize(); - Point3x v0=N; - Point3x v1=(P-p0); - v1.Normalize(); - Point3x v2=(v0^v1); - v2.Normalize(); - - v0=v0*2; - v1=v1*2; - v2=v2*2; - - pp0=-v1-v2; - pp1=-v1+v2; - pp2=v1+v2; - pp3=v1-v2; + D = N; + D0 = D; + dist = vcg::Norm((p0 + p1 + p2) / 3.f); } + //set normal of the clipping plane - void SetD(Point3x d) + void SetD(CoordType d) { - D=d; + D = d; + D0 = d; } //set the point of the clipping plane - void SetP(Point3x p) + void SetDist(ScalarType d) { - P=p; + dist = d; } - bool IsClipped(Point3x p) + bool IsClipped(CoordType p) { - return D.V(0) * p.X() + D.V(1) * p.Y() + D.V(2) * p.Z() - vcg::Norm(P) < 0; + return D.V(0) * p.X() + D.V(1) * p.Y() + D.V(2) * p.Z() - dist > 0; } void GlClip() @@ -132,55 +113,58 @@ public: void GlDraw() { - glPushMatrix(); - glPushAttrib(0xffffffff); - glDisable(GL_CLIP_PLANE0); + const ScalarType w = 50; + glColor4f(1., 1., 0., 0.3); + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + glDisable(GL_LIGHTING); - glEnable(GL_LIGHTING); - glEnable(GL_NORMALIZE); - - glTranslate(P); - glMultMatrix(TR); - glLineWidth(0.5); - glColor3d(0.7,0,0.7); - glBegin(GL_LINE_LOOP); - glVertex(pp0); - glVertex(pp1); - glVertex(pp2); - glVertex(pp3); + glBegin(GL_LINES); + glVertex3f(0, 0, 0); + glVertex3f(dist, 0, 0); + glEnd(); + glBegin(GL_TRIANGLES); + glVertex3f(dist, -w, w); + glVertex3f(dist, w, w); + glVertex3f(dist, -w, -w); + glVertex3f(dist, w, w); + glVertex3f(dist, w, -w); + glVertex3f(dist, -w, -w); glEnd(); - glPopAttrib(); - glPopMatrix(); + glDisable(GL_BLEND); + glEnable(GL_LIGHTING); } - void Transform(vcg::Matrix44 Tr) + void Transform(vcg::Matrix44 & tr) { - //thath's for casting in case of trackball using - //float to double and vice-versa - Point3f p=Point3f((float)D0.V(0),(float)D0.V(1),(float)D0.V(2)); - TR=Tr; - p=TR*p; - D=Point3x((ScalarType) p.V(0),(ScalarType) p.V(1),(ScalarType) p.V(2)); + D = (tr * D0).Normalize(); } - void Translate(float L) + void offsetDist(ScalarType off) { - Point3x D1=D*L; - P+=D1; + dist = (off < -dist) ? 0.001f : dist + off; } + bool IsActive() + { + return active; + } + bool switchActive() + { + return active ^= true; + } }; - GlTetramesh(CONT_TETRA * _t):tetra(_t){} + GlTetramesh(MeshType * m) : _m(m){} GlTetramesh( ) {} - CONT_TETRA * tetra; + MeshType * _m; ClipPlane section; private: - ScalarType shrink_factor = 0.95f; + ScalarType shrink_factor = 0.98f; public: @@ -191,7 +175,7 @@ public: } } - void AddClipSection(Point3x p0,Point3x p1,Point3x p2) + void AddClipSection(CoordType p0, CoordType p1, CoordType p2) { section=ClipPlane(p0,p1,p2); section.active=true; @@ -206,49 +190,34 @@ public: void Draw(){ switch (dm){ case DMNone: break; - case DMSmallTetra:_DrawSmallTetra();break; - case DMFlat:_DrawSurface();break; - case DMWire:_DrawSurface();break; - case DMHidden:_DrawSurface();break; - case DMFlatWire:_DrawFlatWire(); break; - case DMTransparent:break; + case DMSmallTetra: _DrawSmallTetra();break; + case DMFlat: _DrawSurface();break; + case DMWire: _DrawSurface();break; + case DMHidden: _DrawSurface();break; + case DMFlatWire: _DrawFlatWire(); break; + case DMTransparent: break; } } private: template void _DrawSmallTetra(){ - typename CONT_TETRA::iterator it; -// glPushAttrib(0xffff); -// glEnable(GL_COLOR_MATERIAL); -// glEnable(GL_NORMALIZE); -// glPolygonMode(GL_FRONT, GL_FILL); -// glLight(GL_LIGHT0, GL_DIFFUSE, vcg::Color4b::White); -// glEnable(GL_LIGHT0); -// glEnable(GL_LIGHTING); - /*glBegin(GL_TRIANGLES);*/ - for( it = tetra->begin(); it != tetra->end(); ++it) - if((!it->IsD())&&(!(it->IsS()))) //draw as normal + glEnable(GL_LIGHT0); + glEnable(GL_LIGHTING); + ForEachTetra(*_m, [&] (TetraType & t) { + if (!t.IsD()) { - _DrawSmallTetra(*it); + if (!t.IsS()) //draw as normal + _DrawSmallTetra(t); + else //draw in selected mode + _DrawSelectedTetra(t); } - else - if((!it->IsD())&&((it->IsS())))//draw in selection mode - { - _DrawSelectedTetra(*it); - } - //glEnd(); -// glPopAttrib(); - if (section.active) - { -// section.GlClip(); - section.GlDraw(); - } + }); } template - void _DrawFlatWire(){ + void _DrawFlatWire(){ glPushAttrib(0xffff); glEnable(GL_COLOR_MATERIAL); glEnable(GL_DEPTH); @@ -264,7 +233,7 @@ private: template void _DrawSurface(){ - typename CONT_TETRA::iterator it; + glPushAttrib(0xffff); glEnable(GL_COLOR_MATERIAL); @@ -281,10 +250,11 @@ private: glEnable(GL_NORMALIZE); glPolygonMode(GL_FRONT,GL_FILL); } - //glBegin(GL_TRIANGLES); - for( it = tetra->begin(); it != tetra->end(); ++it) - _DrawTetra((*it)); - //glEnd(); + + ForEachTetra(*_m, [&] (TetraType & t) { + _DrawTetra(t); + }); + glPopAttrib(); } @@ -355,14 +325,14 @@ private: } if (cm == CMPerTetra) vcg::glColor(t.C()); -// else -// if(cm == CMPerTetraF) -// { -// Color4b c; -// c = color_tetra(t); -// GLint ic[4]; ic[0] = c[0];ic[1] = c[1];ic[2] = c[2];ic[3] = c[3]; -// glMaterialiv(GL_FRONT,GL_DIFFUSE ,ic); -// } + // else + // if(cm == CMPerTetraF) + // { + // Color4b c; + // c = color_tetra(t); + // GLint ic[4]; ic[0] = c[0];ic[1] = c[1];ic[2] = c[2];ic[3] = c[3]; + // glMaterialiv(GL_FRONT,GL_DIFFUSE ,ic); + // } } template @@ -370,16 +340,16 @@ private: { if (cm!=CMNone) { -// if(cm == CMPerVertexF) -// { -// Color4b c; -// c = color_vertex(v); -// GLint ic[4]; ic[0] = c[0];ic[1] = c[1];ic[2] = c[2];ic[3] = c[3]; -// glMaterialiv(GL_FRONT,GL_DIFFUSE ,ic); -// } -// else - if(cm == CMPerVertex) - vcg::glColor(v.C()); + // if(cm == CMPerVertexF) + // { + // Color4b c; + // c = color_vertex(v); + // GLint ic[4]; ic[0] = c[0];ic[1] = c[1];ic[2] = c[2];ic[3] = c[3]; + // glMaterialiv(GL_FRONT,GL_DIFFUSE ,ic); + // } + // else + if(cm == CMPerVertex) + vcg::glColor(v.C()); } } @@ -424,24 +394,29 @@ private: template < ColorMode cm > void _DrawSmallTetra(TetraType &t) { - Point3x p[4], br; + CoordType p[4], br; br = Tetra::Barycenter(t); - if (section.IsClipped(br)) - return; - bool border = false; - bool clipBorder = false; - for (int i = 0; i < 4; ++i) - { - border = border || t.IsB(i); - Point3x br1 = Tetra::Barycenter(*t.TTp(i)); - clipBorder = clipBorder || section.IsClipped(br1); + if (section.active) + { + if (section.IsClipped(br)) + return; + bool border = false; + bool clipBorder = false; + for (int i = 0; i < 4; ++i) + { + border = border || t.IsB(i); + + CoordType br1 = Tetra::Barycenter(*t.TTp(i)); + clipBorder = clipBorder || section.IsClipped(br1); + } + + if (!border && !clipBorder) + return; } - if (!border && !clipBorder) - return; - for(int i = 0; i < 4; ++i) + // p[i] = t.V(i)->P(); p[i] = t.V(i)->P() * shrink_factor + br * (1 - shrink_factor); diff --git a/wrap/igl/miq_parametrization.h b/wrap/igl/miq_parametrization.h index 6ed928fa..2f52070d 100644 --- a/wrap/igl/miq_parametrization.h +++ b/wrap/igl/miq_parametrization.h @@ -106,7 +106,8 @@ public: for (int j=0;j<3;j++) { //if (!trimesh.face[i].IsB(j))continue; - if (!trimesh.face[i].IsCrease(j))continue; + //if (!trimesh.face[i].IsCrease(j))continue; + if (!trimesh.face[i].IsFaceEdgeS(j))continue; feature_lines.push_back(std::vector()); feature_lines.back().push_back(i); @@ -316,7 +317,8 @@ public: for (int j=0;jClearCrease(j); + //f0->ClearCrease(j); + f0->ClearFaceEdgeS(j); } @@ -326,12 +328,14 @@ public: FaceType *f0=&mesh.face[i]; FaceType *f1=f0->FFp(j); - if (f0==f1){f0->SetCrease(j);continue;} + //if (f0==f1){f0->SetCrease(j);continue;} + if (f0==f1){f0->SetFaceEdgeS(j);continue;} CoordType N0=f0->N(); CoordType N1=f1->N(); if ((N0*N1)>thr)continue; - f0->SetCrease(j); + //f0->SetCrease(j); + f0->SetFaceEdgeS(j); } } diff --git a/wrap/io_tetramesh/import.h b/wrap/io_tetramesh/import.h new file mode 100644 index 00000000..b2aaa7de --- /dev/null +++ b/wrap/io_tetramesh/import.h @@ -0,0 +1,152 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +/**************************************************************************** +****************************************************************************/ + +#ifndef __VCGLIB_TETRA_IMPORT +#define __VCGLIB_TETRA_IMPORT + +#include +#include + +#include + +namespace vcg { +namespace tetra { +namespace io { + +/** +This class encapsulate a filter for automatically importing meshes by guessing +the right filter according to the extension +*/ + +template +class Importer +{ +private: + enum KnownTypes { KT_UNKNOWN, KT_PLY, KT_MSH }; + +static int &LastType() +{ + static int lastType= KT_UNKNOWN; +return lastType; +} + +public: +enum ImporterError { + E_NOERROR =0 // No error =0 is the standard for ALL the imported files. +}; +// simple aux function that returns true if a given file has a given extesnion +static bool FileExtension(std::string filename, std::string extension) +{ + std::transform(filename.begin(), filename.end(), filename.begin(), ::tolower); + std::transform(extension.begin(), extension.end(), extension.begin(), ::tolower); + + std::string end=filename.substr(filename.length() - extension.length(), extension.length()); + + return end == extension; +} +// Open Mesh, returns 0 on success. +static int Open(OpenMeshType &m, const char *filename, CallBackPos *cb=0) +{ + int dummymask = 0; + return Open(m, filename, dummymask, cb); +} + +/// Open Mesh and fills the load mask (the load mask must be initialized first); returns 0 on success. +static int Open(OpenMeshType &m, const char *filename, int &loadmask, CallBackPos *cb=0) +{ + int err; + if(FileExtension(filename,"ply")) + { + err = tetra::io::ImporterPLY::Open(m, filename, loadmask, cb); + LastType()=KT_PLY; + } + else if(FileExtension(filename,"msh")) + { + err = tetra::io::ImporterMSH::Open(m, filename, cb); + LastType()=KT_MSH; + } + else + { + err=1; + LastType()=KT_UNKNOWN; + } + + return err; +} + +static bool ErrorCritical(int error) +{ + switch(LastType()) + { + case KT_PLY : return (error > 0); break; + case KT_MSH : return (error > 0); break; + } + + return true; +} + +static const char *ErrorMsg(int error) +{ + switch(LastType()) + { + // case KT_PLY : return ImporterPLY::ErrorMsg(error); break; + // case KT_STL : return ImporterSTL::ErrorMsg(error); break; + // case KT_OFF : return ImporterOFF::ErrorMsg(error); break; + // case KT_OBJ : return ImporterOBJ::ErrorMsg(error); break; + // case KT_VMI : return ImporterVMI::ErrorMsg(error); break; + } + return "Unknown type"; +} + +static bool LoadMask(const char * filename, int &mask) +{ + bool err; + + if(FileExtension(filename, "ply")) + { + err = tetra::io::ImporterPLY::LoadMask(filename, mask); + LastType() = KT_PLY; + } + else if(FileExtension(filename, "msh")) + { + mask = Mask::IOM_VERTCOORD | Mask::IOM_TETRAINDEX; + err = true; + LastType() = KT_MSH; + } + else + { + err = false; + LastType()=KT_UNKNOWN; + } + + return err; +} +}; // end class +} // end Namespace tetra +} // end Namespace io +} // end Namespace vcg + +#endif diff --git a/wrap/io_tetramesh/import_msh.h b/wrap/io_tetramesh/import_msh.h new file mode 100644 index 00000000..cc7cbf24 --- /dev/null +++ b/wrap/io_tetramesh/import_msh.h @@ -0,0 +1,402 @@ +#ifndef __VCGLIB_IMPORTTETMSH_H +#define __VCGLIB_IMPORTTETMSH_H + +#include + +namespace vcg +{ +namespace tetra +{ +namespace io +{ +template +class ImporterMSH +{ + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::TetraType TetraType; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::TetraIterator TetraIterator; + + enum ErrorCodes + { + INVALID_FORMAT = 1, + INVALID_VERSION, + NOT_IMPLEMENTED, + IO_ERROR + }; + + static inline void parseWhiteSpace(std::ifstream &fin) + { + //we don't want to consume non whitespace bytes, just peek it.. + char next = fin.peek(); + + while (next == '\n' || next == ' ' || next == '\t' || next == '\r') + { + fin.get(); + next = fin.peek(); + } + } + + static int parseNodes(MeshType &m, std::ifstream &fin, bool binary) + { + int numOfNodes; + fin >> numOfNodes; + if (numOfNodes < 0) + return INVALID_FORMAT; + + VertexIterator vi = vcg::tri::Allocator::AddVertices(m, numOfNodes); + + if (binary) + { + size_t lineBytes = (4 + 3 * 8); //int index + 3 * double coords + size_t bytes = numOfNodes * lineBytes; + char *data = new char[bytes]; + parseWhiteSpace(fin); + fin.read(data, bytes); + + for (int i = 0; i < numOfNodes; ++i) + { + int index = *reinterpret_cast(&data[i * lineBytes]) - 1; + if (index < 0) + return INVALID_FORMAT; + + m.vert[index].P().X() = *reinterpret_cast(&data[i * lineBytes + 4]); + m.vert[index].P().Y() = *reinterpret_cast(&data[i * lineBytes + 4 + 8]); + m.vert[index].P().Z() = *reinterpret_cast(&data[i * lineBytes + 4 + 2 * 8]); + } + delete[] data; + } + else + { + for (int i = 0; i < numOfNodes; ++i) + { + int index; + fin >> index; + --index; + + if (index < 0) + return INVALID_FORMAT; + + fin >> m.vert[index].P().X(); + fin >> m.vert[index].P().Y(); + fin >> m.vert[index].P().Z(); + } + } + return 0; + } + + static int parseElements(MeshType &m, std::ifstream &fin, bool binary) + { + int numOfElements; + fin >> numOfElements; + if (numOfElements < 0) + return INVALID_FORMAT; + + TetraIterator ti = vcg::tri::Allocator::AddTetras(m, numOfElements); + + if (binary) + { + parseWhiteSpace(fin); + size_t parsedElems = 0; + + while (parsedElems < numOfElements) + { + //MSH in binary format has a elem-header 3*4 bytes: {elems_type, numElems, tagsPerElem} + //followed by the list of elems under this header and eventually a new header and list. + int type, elements, tags; + fin.read((char *)&type, sizeof(int)); + fin.read((char *)&elements, sizeof(int)); + fin.read((char *)&tags, sizeof(int)); + + //check for tetra type + if (type != 4) + return NOT_IMPLEMENTED; + + //read tags and throw them + for (size_t j = 0; j < tags; ++j) + { + int tag; + fin.read((char *)&tag, sizeof(int)); + } + + //foreach element + for (int i = 0; i < elements; ++i) + { + int index; + fin.read((char *)&index, sizeof(int)); + --index; + + //check index validity + if (index < 0) + return INVALID_FORMAT; + + //read element nodes + TetraType * t = &m.tetra[index]; + for (int i = 0; i < 4; ++i) + { + int nodeIndex; + fin.read((char *)&nodeIndex, sizeof(int)); + --nodeIndex; + + if (nodeIndex < 0 || nodeIndex >= m.VN()) + return INVALID_FORMAT; + + t->V(i) = &m.vert[nodeIndex]; + } + ++parsedElems; + } + } + } + else + { + for (int i = 0; i < numOfElements; ++i) + { + int index, type, tags; + fin >> index >> type >> tags; + --index; + + //check for tetra type + if (type != 4) + return NOT_IMPLEMENTED; + //check index validity + if (index < 0) + return INVALID_FORMAT; + //read tags and throw them + for (size_t j = 0; j < tags; ++j) + { + int tag; + fin >> tag; + } + + TetraType *t = &m.tetra[index]; + for (int i = 0; i < 4; ++i) + { + int nodeIndex; + fin >> nodeIndex; + --nodeIndex; + + if (nodeIndex < 0 || nodeIndex > m.VN()) + return INVALID_FORMAT; + + t->V(i) = &m.vert[nodeIndex]; + } + } + } + + return 0; + } + + static int parseDataField(MeshType &m, std::ifstream &fin, bool binary) + { + int numString, numReal, numInteger; + + fin >> numString; + + std::string *strTags = new std::string[numString]; + for (int i = 0; i < numString; ++i) + { + parseWhiteSpace(fin); + fin >> strTags[i]; + } + + fin >> numReal; + + double *doubleTags = new double[numReal]; + for (int i = 0; i < numReal; ++i) + fin >> doubleTags[i]; + + fin >> numInteger; + + if (numString <= 0 || numInteger < 3) + return INVALID_FORMAT; + + int *integerTags = new int[numInteger]; + + for (int i = 0; i < numInteger; ++i) + fin >> integerTags[i]; + + std::string fieldName = strTags[0]; + int fieldComponents = integerTags[1]; + int fieldSize = integerTags[2]; + + double *fieldVec = new double[fieldComponents * fieldSize]; + + delete[] strTags; + delete[] doubleTags; + delete[] integerTags; + + if (binary) + { + size_t totalBytes = (4 + 8 * fieldComponents) * fieldSize; + char *data = new char[totalBytes]; + parseWhiteSpace(fin); + fin.read(data, totalBytes); + + for (int i = 0; i < fieldSize; ++i) + { + int index = *reinterpret_cast(&data[i * (4 + fieldComponents * 8)]); + --index; + + if (index < 0) + return INVALID_FORMAT; + + //values + int baseIndex = i * (4 + fieldComponents * 8) + 4; + + for (int j = 0; j < fieldComponents; ++j) + fieldVec[index * fieldComponents + j] = *reinterpret_cast(&data[baseIndex + j * 8]); + } + } + else + { + for (int i = 0; i < fieldSize; ++i) + { + int index; + fin >> index; + --index; + + if (index < 0) + return INVALID_FORMAT; + + for (int j = 0; j < fieldComponents; ++j) + fin >> fieldVec[index * fieldComponents + j]; + } + } + } + + static int parseNodeData(MeshType &m, std::ifstream &fin, bool binary) + { + return parseDataField(m, fin, binary); + } + + static int parseElementData(MeshType &m, std::ifstream &fin, bool binary) + { + return parseDataField(m, fin, binary); + } + + static int parseUnsupportedTag(std::ifstream &fin, std::string &tag) + { + std::cerr << "found unsupported tag" << std::endl; + std::string tagName = tag.substr(1, tag.size() - 1); + + std::string tagEnd = tag.substr(0, 1) + "End" + tagName; + + std::string buf; + while (buf != tagEnd && !fin.eof()) + fin >> buf; + } + + static int parseMshMesh(MeshType &m, std::string &filename) + { + std::ifstream fin(filename.c_str(), std::ios::in | std::ios::binary); + + if (!fin.is_open()) + return IO_ERROR; + + std::string lookAhead; + + fin >> lookAhead; + if (lookAhead != "$MeshFormat") + return INVALID_FORMAT; + + double version; + int type, dataSize; + + fin >> version >> type >> dataSize; + + if (version != 2.2) + return INVALID_VERSION; + + bool binary = (type == 1); + + if (dataSize != 8) + return INVALID_FORMAT; + + // Read endiannes info in binary header...it's a 1 used to detect endiannes. + if (binary) + { + int one; + parseWhiteSpace(fin); + fin.read(reinterpret_cast(&one), sizeof(int)); + if (one != 1) + { + std::cerr << "Warning: binary msh file " << filename + << " is saved with different endianness than this machine." + << std::endl; + throw NOT_IMPLEMENTED; + } + } + + lookAhead.clear(); + fin >> lookAhead; + if (lookAhead != "$EndMeshFormat") + return INVALID_FORMAT; + + while (!fin.eof()) + { + lookAhead.clear(); + fin >> lookAhead; + + if (lookAhead == "$Nodes") + { + int res = parseNodes(m, fin, binary); + if (res != 0) + return res; + + fin >> lookAhead; + if (lookAhead != "$EndNodes") + return INVALID_FORMAT; + } + else if (lookAhead == "$Elements") + { + int res = parseElements(m, fin, binary); + + if (res != 0) + return res; + + fin >> lookAhead; + if (lookAhead != "$EndElements") + return INVALID_FORMAT; + } + else if (lookAhead == "$NodeData") + { + parseNodeData(m, fin, binary); + fin >> lookAhead; + if (lookAhead != "$EndNodeData") + return INVALID_FORMAT; + } + else if (lookAhead == "$ElementData") + { + parseElementData(m, fin, binary); + fin >> lookAhead; + if (lookAhead != "$EndElementData") + return INVALID_FORMAT; + + } + else if (fin.eof()) + { + break; + } + else + { + parseUnsupportedTag(fin, lookAhead); + } + } + fin.close(); + return 0; + } + + public: + static int Open(MeshType &m, const char *filename, CallBackPos *cb = 0) + { + std::string name(filename); + return parseMshMesh(m, name); + } +}; +} // namespace io +} // namespace tetra +} // namespace vcg + +#endif diff --git a/wrap/io_tetramesh/import_ply.h b/wrap/io_tetramesh/import_ply.h index 325a23e1..af7e2fe9 100644 --- a/wrap/io_tetramesh/import_ply.h +++ b/wrap/io_tetramesh/import_ply.h @@ -141,16 +141,21 @@ public: static const PropDescriptor &VertDesc(int i) { - const static PropDescriptor pv[9]={ - {"vertex", "x", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p),0,0,0,0,0}, - {"vertex", "y", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p) + sizeof(ScalarType),0,0,0,0,0}, - {"vertex", "z", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p) + 2 * sizeof(ScalarType),0,0,0,0,0}, - {"vertex", "flags", ply::T_INT, ply::T_INT, offsetof(LoadPly_VertAux,flags),0,0,0,0,0}, - {"vertex", "quality", ply::T_FLOAT, ply::T_FLOAT, offsetof(LoadPly_VertAux,q),0,0,0,0,0}, - {"vertex", "red" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,r),0,0,0,0,0}, - {"vertex", "green", ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,g),0,0,0,0,0}, - {"vertex", "blue" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,b),0,0,0,0,0}, - {"vertex", "alpha" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,a),0,0,0,0,0}, + const static PropDescriptor pv[13]={ + /*00*/ {"vertex", "x", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p),0,0,0,0,0}, + /*01*/ {"vertex", "y", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p) + sizeof(ScalarType),0,0,0,0,0}, + /*02*/ {"vertex", "z", ply::T_FLOAT, PlyType(), offsetof(LoadPly_VertAux,p) + 2 * sizeof(ScalarType),0,0,0,0,0}, + /*03*/ {"vertex", "flags", ply::T_INT, ply::T_INT, offsetof(LoadPly_VertAux,flags),0,0,0,0,0}, + /*04*/ {"vertex", "quality", ply::T_FLOAT, ply::T_FLOAT, offsetof(LoadPly_VertAux,q),0,0,0,0,0}, + /*05*/ {"vertex", "red" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,r),0,0,0,0,0}, + /*06*/ {"vertex", "green", ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,g),0,0,0,0,0}, + /*07*/ {"vertex", "blue" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,b),0,0,0,0,0}, + /*08*/ {"vertex", "alpha" , ply::T_UCHAR, ply::T_UCHAR, offsetof(LoadPly_VertAux,a),0,0,0,0,0}, + // DOUBLE + /*09*/ {"vertex", "x", ply::T_DOUBLE, PlyType(),offsetof(LoadPly_VertAux,p),0,0,0,0,0 ,0}, + /*10*/ {"vertex", "y", ply::T_DOUBLE, PlyType(),offsetof(LoadPly_VertAux,p) + sizeof(ScalarType) ,0,0,0,0,0 ,0}, + /*11*/ {"vertex", "z", ply::T_DOUBLE, PlyType(),offsetof(LoadPly_VertAux,p) + 2*sizeof(ScalarType),0,0,0,0,0 ,0}, + /*12*/ {"vertex", "quality", ply::T_DOUBLE, PlyType(),offsetof(LoadPly_VertAux,q),0,0,0,0,0 ,0} }; return pv[i]; } @@ -224,18 +229,18 @@ public: pi.header = pf.GetHeader(); // Descrittori dati standard (vertex coord e faces) - if( pf.AddToRead(VertDesc(0)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return -1; } - if( pf.AddToRead(VertDesc(1)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return -1; } - if( pf.AddToRead(VertDesc(2)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return -1; } - if( pf.AddToRead(TetraDesc(0))== -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return -1; } + if( pf.AddToRead(VertDesc(0)) == -1 && pf.AddToRead(VertDesc(9)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return pi.status; } + if( pf.AddToRead(VertDesc(1)) == -1 && pf.AddToRead(VertDesc(10)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return pi.status; } + if( pf.AddToRead(VertDesc(2)) == -1 && pf.AddToRead(VertDesc(11)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return pi.status; } + if( pf.AddToRead(TetraDesc(0)) == -1 ) { pi.status = PlyInfo::E_NO_VERTEX; return pi.status; } // Descrittori facoltativi dei flags - if( pf.AddToRead(VertDesc(3))!=-1 ) + if( pf.AddToRead(VertDesc(3)) != -1 ) pi.mask |= vcg::tetra::io::Mask::IOM_VERTFLAGS; if( VertexType::HasQuality() ) { - if( pf.AddToRead(VertDesc(4))!=-1 || - pf.AddToRead(VertDesc(9))!=-1 ) + if( pf.AddToRead(VertDesc(4)) != -1 || + pf.AddToRead(VertDesc(9)) != -1 ) pi.mask |= vcg::tetra::io::Mask::IOM_VERTQUALITY; } if( VertexType::HasColor() )