2 * This program is free software: you can redistribute it and/or modify
3 * it under the terms of the GNU Lesser General Public License as
4 * published by the Free Software Foundation, either version 3 of the
5 * License, or (at your option) any later version.
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
12 * You should have received a copy of the GNU General Public License
13 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 * @author Bertrand Kerautret (\c kerautre@loria.fr )
20 * LORIA (CNRS, UMR 7503), University of Nancy, France
24 * Implementation of inline methods defined in Mesh.h
26 * This file is part of the DGtal library.
30//////////////////////////////////////////////////////////////////////////////
34#include <DGtal/kernel/BasicPointPredicates.h>
35//////////////////////////////////////////////////////////////////////////////
37///////////////////////////////////////////////////////////////////////////////
38// IMPLEMENTATION of inline methods.
39///////////////////////////////////////////////////////////////////////////////
41///////////////////////////////////////////////////////////////////////////////
42// ----------------------- Standard services ------------------------------
48template <typename TPoint>
50DGtal::Mesh<TPoint>::Mesh(bool saveFaceColor)
52 mySaveFaceColor=saveFaceColor;
53 myDefaultColor = DGtal::Color::White;
59template <typename TPoint>
61DGtal::Mesh<TPoint>::Mesh(const DGtal::Color &aColor)
63 mySaveFaceColor=false;
64 myDefaultColor = aColor;
70template <typename TPoint>
72DGtal::Mesh<TPoint>::~Mesh()
77template <typename TPoint>
79DGtal::Mesh<TPoint>::Mesh ( const Mesh & other ): myFaceList(other.myFaceList),
80 myVertexList(other.myVertexList),
81 myFaceColorList(other.myFaceColorList),
82 mySaveFaceColor(other.mySaveFaceColor),
83 myDefaultColor(other.myDefaultColor)
88template <typename TPoint>
91DGtal::Mesh<TPoint>::operator= ( const Mesh & other )
93 myFaceList = other.myFaceList;
94 myVertexList = other.myVertexList;
95 myFaceColorList = other.myFaceColorList;
96 mySaveFaceColor = other.mySaveFaceColor;
97 myDefaultColor = other. myDefaultColor;
102///////////////////////////////////////////////////////////////////////////////
103// Interface - public :
106 * Writes/Displays the object on an output stream.
107 * @param out the output stream where the object is written.
109template <typename TPoint>
112DGtal::Mesh<TPoint>::selfDisplay ( std::ostream & out ) const
118 * Checks the validity/consistency of the object.
119 * @return 'true' if the object is valid, 'false' otherwise.
121template <typename TPoint>
124DGtal::Mesh<TPoint>::isValid() const
132///////////////////////////////////////////////////////////////////////////////
133// Implementation of inline functions //
141template<typename TPoint>
143DGtal::Mesh<TPoint>::Mesh(const VertexStorage &vertexSet)
145 mySaveFaceColor=false;
146 for(int i =0; i< vertexSet.size(); i++)
148 myVertexList.push_back(vertexSet.at(i));
155template<typename TPoint>
158DGtal::Mesh<TPoint>::addVertex(const TPoint &point)
160 myVertexList.push_back(point);
165template<typename TPoint>
168DGtal::Mesh<TPoint>::addTriangularFace(Index indexVertex1, Index indexVertex2,
169 Index indexVertex3, const DGtal::Color &aColor)
172 aFace.push_back(indexVertex1);
173 aFace.push_back(indexVertex2);
174 aFace.push_back(indexVertex3);
175 myFaceList.push_back(aFace);
178 myFaceColorList.push_back(aColor);
185template<typename TPoint>
188DGtal::Mesh<TPoint>::addQuadFace(Index indexVertex1,Index indexVertex2,
189 Index indexVertex3, Index indexVertex4,
190 const DGtal::Color &aColor)
193 aFace.push_back(indexVertex1);
194 aFace.push_back(indexVertex2);
195 aFace.push_back(indexVertex3);
196 aFace.push_back(indexVertex4);
197 myFaceList.push_back(aFace);
200 myFaceColorList.push_back(aColor);
206template<typename TPoint>
209DGtal::Mesh<TPoint>::addFace(const MeshFace &aFace, const DGtal::Color &aColor){
210 myFaceList.push_back(aFace);
213 myFaceColorList.push_back(aColor);
219template<typename TPoint>
222DGtal::Mesh<TPoint>::removeFaces(const std::vector<Index> &facesIndex){
223 DGtal::Mesh<TPoint> newMesh(true);
225 std::vector<unsigned int> indexVertexFaceCard(nbVertex());
226 std::vector<bool> indexFaceOK(nbFaces());
227 std::fill(indexVertexFaceCard.begin(), indexVertexFaceCard.end(), 0);
228 std::fill(indexFaceOK.begin(), indexFaceOK.end(), true);
229 for (unsigned int i = 0; i<facesIndex.size(); i++){
230 indexFaceOK[facesIndex[i]]=false;
232 // for each face remaining in the mesh we add +1 to each vertex used in a face
233 for(unsigned int i = 0; i < nbFaces(); i++){
234 if( indexFaceOK[i] ){
235 DGtal::Mesh<TPoint>::MeshFace aFace = getFace(i);
236 for (unsigned int j=0; j< aFace.size() ; j++) {
237 indexVertexFaceCard[aFace[j]] += 1;
241 // we remove all vertex with a face == 0 and compute the new vertex association:
242 std::vector<unsigned int> newVertexIndex;
243 unsigned int currentIndex=0;
244 for (unsigned int i=0; i< nbVertex(); i++) {
245 if (indexVertexFaceCard[i]!=0){
246 newMesh.addVertex(getVertex(i));
247 newVertexIndex.push_back(currentIndex);
250 newVertexIndex.push_back(0);
253 for (unsigned int i = 0; i < nbFaces(); i++) {
255 MeshFace aFace = getFace(i);
256 MeshFace aNewFace = aFace;
257 // translate the old face with new index:
258 for (unsigned int j=0; j< aFace.size() ; j++) {
259 aNewFace[j] = newVertexIndex[aFace[j]];
261 newMesh.addFace(aNewFace);
262 newMesh.setFaceColor(newMesh.nbFaces()-1, getFaceColor(i));
265 myFaceList = newMesh.myFaceList;
266 myVertexList = newMesh.myVertexList;
267 myFaceColorList = newMesh.myFaceColorList;
274template<typename TPoint>
277DGtal::Mesh<TPoint>::getVertex(Index i) const
279 return myVertexList.at(i);
284template<typename TPoint>
287DGtal::Mesh<TPoint>::getVertex(Index i)
289 return myVertexList.at(i);
294template<typename TPoint>
296const typename DGtal::Mesh<TPoint>::MeshFace &
297DGtal::Mesh<TPoint>::getFace(Index i) const
299 return myFaceList.at(i);
303template<typename TPoint>
305typename DGtal::Mesh<TPoint>::MeshFace &
306DGtal::Mesh<TPoint>::getFace(Index i)
308 return myFaceList.at(i);
312template<typename TPoint>
314typename DGtal::Mesh<TPoint>::RealPoint
315DGtal::Mesh<TPoint>::getFaceBarycenter(Index i) const
317 DGtal::Mesh<TPoint>::RealPoint c;
318 MeshFace aFace = getFace(i);
319 for ( auto &j: aFace){
320 TPoint p = getVertex(j);
321 for (typename TPoint::Dimension k = 0; k < TPoint::dimension; k++){
322 c[k] += static_cast<typename RealPoint::Component>(p[k]) ;
325 return c/static_cast<typename RealPoint::Component>(aFace.size());
329template<typename TPoint>
331typename DGtal::Mesh<TPoint>::Size
332DGtal::Mesh<TPoint>::nbFaces() const
334 return myFaceList.size();
337template<typename TPoint>
339typename DGtal::Mesh<TPoint>::Size
340DGtal::Mesh<TPoint>::nbVertex() const
342 return myVertexList.size();
346template<typename TPoint>
349DGtal::Mesh<TPoint>::getFaceColor(Index i) const
353 return myFaceColorList.at(i);
357 return myDefaultColor;
361template <typename TPoint>
362struct MeshBoundingBoxCompPoints
364 MeshBoundingBoxCompPoints(typename TPoint::Dimension d): myDim(d){};
365 bool operator() (const TPoint &p1, const TPoint &p2){return p1[myDim]<p2[myDim];};
366 typename TPoint::Dimension myDim;
369template<typename TPoint>
371std::pair<TPoint, TPoint>
372DGtal::Mesh<TPoint>::getBoundingBox() const
374 std::pair<TPoint, TPoint> theResult;
375 TPoint lowerBound, upperBound;
376 for(unsigned int i=0; i< TPoint::size(); i++)
378 const MeshBoundingBoxCompPoints<TPoint> cmp_points(i);
379 upperBound[i] = (*(std::max_element(vertexBegin(), vertexEnd(), cmp_points)))[i];
380 lowerBound[i] = (*(std::min_element(vertexBegin(), vertexEnd(), cmp_points)))[i];
382 theResult.first = lowerBound ;
383 theResult.second = upperBound ;
388template<typename TPoint>
391DGtal::Mesh<TPoint>::setFaceColor(const Index index,
392 const DGtal::Color &aColor)
394 if (!mySaveFaceColor)
396 for(unsigned int i = 0; i<myFaceList.size(); i++)
398 myFaceColorList.push_back(myDefaultColor);
400 mySaveFaceColor=true;
402 myFaceColorList.at(index) = aColor;
406template<typename TPoint>
409DGtal::Mesh<TPoint>::isStoringFaceColors() const
411 return mySaveFaceColor;
417template<typename TPoint>
420DGtal::Mesh<TPoint>::invertVertexFaceOrder(){
421 for(unsigned int i=0; i<myFaceList.size(); i++)
423 auto aFace = myFaceList.at(i);
424 for(unsigned int j=0; j < aFace.size()/2; j++)
426 const auto tmp=aFace.at(j);
427 aFace.at(j)=aFace.at(aFace.size()-1-j);
428 aFace.at(aFace.size()-1-j)=tmp;
433template<typename TPoint>
436DGtal::Mesh<TPoint>::clearFaces(){
441template<typename TPoint>
444DGtal::Mesh<TPoint>::clearVertices(){
445 myVertexList.clear();
448template <typename TPoint>
450DGtal::Mesh<TPoint>::removeIsolatedVertices(){
451 typedef typename Mesh<TPoint>::Index MIndex;
452 DGtal::Mesh<TPoint>::VertexStorage vSt;
453 std::vector<bool> vertexUsed (nbVertex(), false);
454 for ( MIndex f = 0; f< nbFaces(); f++ )
456 auto face = getFace(f);
457 for (MIndex i = 0; i<face.size(); i++)
459 vertexUsed[face[i]] = true;
462 std::vector<MIndex> translateIndexId;
463 MIndex currentIndex = 0;
464 MIndex nbV = nbVertex();
465 for(MIndex i = 0; i < nbV; i++ )
469 translateIndexId.push_back(currentIndex);
470 vSt.push_back(myVertexList[i]);
475 translateIndexId.push_back(0);
479 for ( MIndex f = 0; f< nbFaces(); f++ )
481 auto &face = getFace(f);
482 for (MIndex i = 0; i<face.size(); i++)
484 face[i]=translateIndexId[face[i]];
489template<typename TPoint>
492DGtal::Mesh<TPoint>::rescale(const typename TPoint::Component aScale){
493 for(typename VertexStorage::iterator it = vertexBegin(); it != vertexEnd(); it++)
500template<typename TPoint>
503DGtal::Mesh<TPoint>::subDivideTriangularFaces(const double minArea){
505 std::vector<Mesh<TPoint>::MeshFace> facesToAdd;
506 for(unsigned int i =0; i< nbFaces(); i++)
508 typename Mesh<TPoint>::MeshFace aFace = getFace(i);
511 TPoint p1 = getVertex(aFace[0]);
512 TPoint p2 = getVertex(aFace[1]);
513 TPoint p3 = getVertex(aFace[2]);
514 TPoint c = (p1+p2+p3)/3.0;
515 double a = ((p2-p1).crossProduct(p3-p1)).norm()/2.0;
525 f1.push_back(aFace[0]);
526 f1.push_back(aFace[1]);
527 f1.push_back(nbVertex()-1);
528 facesToAdd.push_back(f1);
530 f2.push_back(aFace[1]);
531 f2.push_back(aFace[2]);
532 f2.push_back(nbVertex()-1);
533 facesToAdd.push_back(f2);
535 f3.push_back(aFace[2]);
536 f3.push_back(aFace[0]);
537 f3.push_back(nbVertex()-1);
538 facesToAdd.push_back(f3);
542 facesToAdd.push_back(aFace);
548 for(unsigned i=0; i<facesToAdd.size(); i++)
550 addFace(facesToAdd[i]);
557template<typename TPoint>
560DGtal::Mesh<TPoint>::quadToTriangularFaces(){
561 unsigned int nbQuadT=0;
562 std::vector<Mesh<TPoint>::MeshFace> facesToAdd;
563 for(unsigned int i =0; i< nbFaces(); i++)
565 typename Mesh<TPoint>::MeshFace aFace = getFace(i);
569 f1.push_back(aFace[0]);
570 f1.push_back(aFace[1]);
571 f1.push_back(aFace[2]);
572 facesToAdd.push_back(f1);
574 f2.push_back(aFace[2]);
575 f2.push_back(aFace[3]);
576 f2.push_back(aFace[0]);
577 facesToAdd.push_back(f2);
582 facesToAdd.push_back(aFace);
586 for(unsigned i=0; i<facesToAdd.size(); i++)
588 addFace(facesToAdd[i]);
596//------------------------------------------------------------------------------
597template<typename TPoint>
600DGtal::Mesh<TPoint>::className() const
610template <typename TPoint>
613DGtal::Mesh<TPoint>::createTubularMesh(DGtal::Mesh<TPoint> &aMesh, const std::vector<TPoint> &aSkeleton,
614 const double aRadius,
615 const double angleStep, const DGtal::Color &aMeshColor)
617 std::vector<double> aVecR;
618 aVecR.push_back(aRadius);
619 DGtal::Mesh<TPoint>::createTubularMesh(aMesh, aSkeleton,
620 aVecR, angleStep, aMeshColor);
625template <typename TPoint>
628DGtal::Mesh<TPoint>::createTubularMesh(DGtal::Mesh<TPoint> &aMesh, const std::vector<TPoint> &aSkeleton,
629 const std::vector<double> &aVectRadius,
630 const double angleStep, const DGtal::Color &aMeshColor)
632 auto nbVertexInitial = aMesh.nbVertex();
633 ASSERT(aVectRadius.size() > 0);
634 // Generating vertices..
635 for(auto i = 0; i< (int)aSkeleton.size(); i++)
638 TPoint uDir1, uDirPrec;
642 if(i != (int)aSkeleton.size()-1)
644 vectDir = aSkeleton.at(i+1) - aSkeleton.at(i);
648 vectDir = aSkeleton.at(i) - aSkeleton.at(i-1);
651 double d = -vectDir[0]* aSkeleton.at(i)[0] - vectDir[1]*aSkeleton.at(i)[1]
652 - vectDir[2]*aSkeleton.at(i)[2];
656 pRefOrigin [0]= -d/vectDir[0];
659 if(aSkeleton.at(i) == pRefOrigin ||
660 (vectDir[1]==0 && vectDir[2]==0))
666 else if (vectDir[1]!=0)
669 pRefOrigin [1]= -d/vectDir[1];
671 if(aSkeleton.at(i) == pRefOrigin ||
672 (vectDir[0]==0 && vectDir[2]==0))
676 }else if (vectDir[2]!=0)
680 pRefOrigin [2]= -d/vectDir[2];
681 if(aSkeleton.at(i) == pRefOrigin ||
682 (vectDir[0]==0 && vectDir[1]==0))
687 uDir1=(pRefOrigin-aSkeleton.at(i))/((pRefOrigin-aSkeleton.at(i)).norm());
688 uDir2[0] = uDir1[1]*vectDir[2]-uDir1[2]*vectDir[1];
689 uDir2[1] = uDir1[2]*vectDir[0]-uDir1[0]*vectDir[2];
690 uDir2[2] = uDir1[0]*vectDir[1]-uDir1[1]*vectDir[0];
692 for(double a = 0.0; a < 2.0*M_PI; a += angleStep)
694 TPoint vMove = aVectRadius.at(i%aVectRadius.size())*(uDir1*cos(a) + uDir2*sin(a));
695 aMesh.addVertex(vMove + aSkeleton[i]);
698 firstPoint = vMove + aSkeleton[i]+vectDir;
703 unsigned int nbPtPerFaces = static_cast<unsigned int>((aMesh.nbVertex()-nbVertexInitial)/aSkeleton.size());
705 // Generating faces...
706 for(auto i = 0; i< (int)aSkeleton.size()-1; i++)
708 if (aSkeleton.at(i)==aSkeleton.at(i+1)){
709 trace.warning() << "Two skeleton points are identical, ignoring one point." << std::endl;
712 // Computing best shift between two consecutive ring points to generate tube face.
713 // (criteria defined by the minimal distance between 4 sampling points)
714 double minDistance = std::numeric_limits<double>::max();
715 TPoint ptRefRing1 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces);
716 TPoint ptRefRing2 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+nbPtPerFaces/4);
717 TPoint ptRefRing3 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+2*(nbPtPerFaces/4));
718 TPoint ptRefRing4 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+3*(nbPtPerFaces/4));
720 unsigned int shift = 0;
722 if(i != (int)aSkeleton.size()-1)
724 vectDir = aSkeleton.at(i+1) - aSkeleton.at(i);
728 vectDir = aSkeleton.at(i) - aSkeleton.at(i-1);
731 for(unsigned int k=0; k<nbPtPerFaces; k++)
733 TPoint pScan1 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+k);
734 TPoint pScan2 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
735 (nbPtPerFaces/4+k)%nbPtPerFaces);
736 TPoint pScan3 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
737 (2*(nbPtPerFaces/4)+k)%nbPtPerFaces);
738 TPoint pScan4 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
739 (3*(nbPtPerFaces/4)+k)%nbPtPerFaces);
740 double distance = (ptRefRing1 - pScan1).norm()+(ptRefRing2 - pScan2).norm()+
741 (ptRefRing3 - pScan3).norm()+(ptRefRing4 - pScan4).norm();
742 if(distance<minDistance){
744 minDistance = distance;
748 for(unsigned int k=0; k<nbPtPerFaces; k++)
750 Mesh<TPoint>::MeshFace aFace;
751 aMesh.addQuadFace(nbVertexInitial+k+i*nbPtPerFaces,
752 nbVertexInitial+(shift+k)%nbPtPerFaces+nbPtPerFaces*(i+1),
753 nbVertexInitial+(shift+k+1)%nbPtPerFaces+nbPtPerFaces*(i+1),
754 nbVertexInitial+(k+1)%nbPtPerFaces+i*nbPtPerFaces,
764template <typename TPoint>
765template <typename TValue>
768DGtal::Mesh<TPoint>::createMeshFromHeightSequence(Mesh<TPoint> &aMesh, const std::vector<TValue> & anValueSequence,
769 const unsigned int lengthSequence,
770 double stepX, double stepY, double stepZ,
771 const DGtal::Color &aMeshColor ){
772 const auto nbVertexInitial = aMesh.nbVertex();
773 // Generating vertices..
775 unsigned int posY = 0;
776 while(i+(int)lengthSequence-1 < (int)anValueSequence.size()){
777 for(unsigned int j = 0; j < lengthSequence; j++, i++){
778 aMesh.addVertex(TPoint(j*stepX, posY*stepY, stepZ*anValueSequence.at(i)));
782 // Generating faces...
785 while(i+(int)lengthSequence-1 < (int)anValueSequence.size() - (int)lengthSequence){
786 for(auto j = 0; j < (int)lengthSequence-1; j++, i++){
787 aMesh.addQuadFace(nbVertexInitial+i, nbVertexInitial+i+1,
788 nbVertexInitial+i+1+lengthSequence,
789 nbVertexInitial+i+lengthSequence,
800template <typename TPoint>
803DGtal::operator<< ( std::ostream & out,
804 const Mesh<TPoint> & object )
806 object.selfDisplay( out );
815///////////////////////////////////////////////////////////////////////////////