DGtal 1.3.0
Loading...
Searching...
No Matches
Mesh.ih
1/**
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.
6 *
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.
11 *
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/>.
14 *
15 **/
16
17/**
18 * @file Mesh.ih
19 * @author Bertrand Kerautret (\c kerautre@loria.fr )
20 * LORIA (CNRS, UMR 7503), University of Nancy, France
21 *
22 * @date 2012/06/29
23 *
24 * Implementation of inline methods defined in Mesh.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <limits>
32#include <cstdlib>
33#include <DGtal/kernel/BasicPointPredicates.h>
34//////////////////////////////////////////////////////////////////////////////
35
36///////////////////////////////////////////////////////////////////////////////
37// IMPLEMENTATION of inline methods.
38///////////////////////////////////////////////////////////////////////////////
39
40///////////////////////////////////////////////////////////////////////////////
41// ----------------------- Standard services ------------------------------
42
43
44/**
45 * Constructor.
46 */
47template <typename TPoint>
48inline
49DGtal::Mesh<TPoint>::Mesh(bool saveFaceColor)
50{
51 mySaveFaceColor=saveFaceColor;
52 myDefaultColor = DGtal::Color::White;
53}
54
55/**
56 * Constructor.
57 */
58template <typename TPoint>
59inline
60DGtal::Mesh<TPoint>::Mesh(const DGtal::Color &aColor)
61{
62 mySaveFaceColor=false;
63 myDefaultColor = aColor;
64}
65
66/**
67 * Destructor.
68 */
69template <typename TPoint>
70inline
71DGtal::Mesh<TPoint>::~Mesh()
72{
73}
74
75
76template <typename TPoint>
77inline
78DGtal::Mesh<TPoint>::Mesh ( const Mesh & other ): myFaceList(other.myFaceList),
79 myVertexList(other.myVertexList),
80 myFaceColorList(other.myFaceColorList),
81 mySaveFaceColor(other.mySaveFaceColor),
82 myDefaultColor(other.myDefaultColor)
83{
84
85}
86
87template <typename TPoint>
88inline
89DGtal::Mesh<TPoint> &
90DGtal::Mesh<TPoint>::operator= ( const Mesh & other )
91{
92 myFaceList = other.myFaceList;
93 myVertexList = other.myVertexList;
94 myFaceColorList = other.myFaceColorList;
95 mySaveFaceColor = other.mySaveFaceColor;
96 myDefaultColor = other. myDefaultColor;
97 return *this;
98}
99
100
101///////////////////////////////////////////////////////////////////////////////
102// Interface - public :
103
104/**
105 * Writes/Displays the object on an output stream.
106 * @param out the output stream where the object is written.
107 */
108template <typename TPoint>
109inline
110void
111DGtal::Mesh<TPoint>::selfDisplay ( std::ostream & out ) const
112{
113 out << "[Mesh]";
114}
115
116/**
117 * Checks the validity/consistency of the object.
118 * @return 'true' if the object is valid, 'false' otherwise.
119 */
120template <typename TPoint>
121inline
122bool
123DGtal::Mesh<TPoint>::isValid() const
124{
125 return true;
126}
127
128
129
130
131///////////////////////////////////////////////////////////////////////////////
132// Implementation of inline functions //
133
134
135
136
137
138
139
140template<typename TPoint>
141inline
142DGtal::Mesh<TPoint>::Mesh(const VertexStorage &vertexSet)
143{
144 mySaveFaceColor=false;
145 for(int i =0; i< vertexSet.size(); i++)
146 {
147 myVertexList.push_back(vertexSet.at(i));
148 }
149
150}
151
152
153
154template<typename TPoint>
155inline
156void
157DGtal::Mesh<TPoint>::addVertex(const TPoint &point)
158{
159 myVertexList.push_back(point);
160}
161
162
163
164template<typename TPoint>
165inline
166void
167DGtal::Mesh<TPoint>::addTriangularFace(Index indexVertex1, Index indexVertex2,
168Index indexVertex3, const DGtal::Color &aColor)
169{
170 MeshFace aFace;
171 aFace.push_back(indexVertex1);
172 aFace.push_back(indexVertex2);
173 aFace.push_back(indexVertex3);
174 myFaceList.push_back(aFace);
175 if(mySaveFaceColor)
176 {
177 myFaceColorList.push_back(aColor);
178 }
179}
180
181
182
183
184template<typename TPoint>
185inline
186void
187DGtal::Mesh<TPoint>::addQuadFace(Index indexVertex1,Index indexVertex2,
188 Index indexVertex3, Index indexVertex4,
189 const DGtal::Color &aColor)
190{
191 MeshFace aFace;
192 aFace.push_back(indexVertex1);
193 aFace.push_back(indexVertex2);
194 aFace.push_back(indexVertex3);
195 aFace.push_back(indexVertex4);
196 myFaceList.push_back(aFace);
197 if(mySaveFaceColor)
198 {
199 myFaceColorList.push_back(aColor);
200 }
201}
202
203
204
205template<typename TPoint>
206inline
207void
208DGtal::Mesh<TPoint>::addFace(const MeshFace &aFace, const DGtal::Color &aColor){
209 myFaceList.push_back(aFace);
210 if(mySaveFaceColor)
211 {
212 myFaceColorList.push_back(aColor);
213 }
214}
215
216
217
218template<typename TPoint>
219inline
220void
221DGtal::Mesh<TPoint>::removeFaces(const std::vector<Index> &facesIndex){
222 DGtal::Mesh<TPoint> newMesh(true);
223 std::vector<unsigned int> indexVertexFaceCard(nbVertex());
224 std::vector<bool> indexFaceOK(nbFaces());
225 std::fill(indexVertexFaceCard.begin(), indexVertexFaceCard.end(), 0);
226 std::fill(indexFaceOK.begin(), indexFaceOK.end(), true);
227 for (unsigned int i = 0; i<facesIndex.size(); i++){
228 indexFaceOK[facesIndex[i]]=false;
229 }
230 // for each face remaining in the mesh we add +1 to each vertex used in a face
231 for(unsigned int i = 0; i < nbFaces(); i++){
232 if( indexFaceOK[i] ){
233 DGtal::Mesh<TPoint>::MeshFace aFace = getFace(i);
234 for (unsigned int j=0; j< aFace.size() ; j++) {
235 indexVertexFaceCard[aFace[j]] += 1;
236 }
237 }
238 }
239 // we remove all vertex with a face == 0 and compute the new vertex association:
240 std::vector<unsigned int> newVertexIndex;
241 unsigned int currentIndex=0;
242 for (unsigned int i=0; i< nbVertex(); i++) {
243 if (indexVertexFaceCard[i]!=0){
244 newMesh.addVertex(getVertex(i));
245 newVertexIndex.push_back(currentIndex);
246 currentIndex++;
247 }else{
248 newVertexIndex.push_back(0);
249 }
250 }
251 for (unsigned int i = 0; i < nbFaces(); i++) {
252 if(indexFaceOK[i]){
253 MeshFace aFace = getFace(i);
254 MeshFace aNewFace = aFace;
255 // translate the old face with new index:
256 for (unsigned int j=0; j< aFace.size() ; j++) {
257 aNewFace[j] = newVertexIndex[aFace[j]];
258 }
259 newMesh.addFace(aNewFace);
260 newMesh.setFaceColor(newMesh.nbFaces()-1, getFaceColor(i));
261 }
262 }
263 myFaceList = newMesh.myFaceList;
264 myVertexList = newMesh.myVertexList;
265 myFaceColorList = newMesh.myFaceColorList;
266}
267
268
269
270
271
272template<typename TPoint>
273inline
274const TPoint &
275DGtal::Mesh<TPoint>::getVertex(Index i) const
276{
277 return myVertexList.at(i);
278}
279
280
281
282template<typename TPoint>
283inline
284TPoint &
285DGtal::Mesh<TPoint>::getVertex(Index i)
286{
287 return myVertexList.at(i);
288}
289
290
291
292template<typename TPoint>
293inline
294const typename DGtal::Mesh<TPoint>::MeshFace &
295DGtal::Mesh<TPoint>::getFace(Index i) const
296{
297 return myFaceList.at(i);
298}
299
300
301template<typename TPoint>
302inline
303typename DGtal::Mesh<TPoint>::MeshFace &
304DGtal::Mesh<TPoint>::getFace(Index i)
305{
306 return myFaceList.at(i);
307}
308
309
310template<typename TPoint>
311inline
312typename DGtal::Mesh<TPoint>::RealPoint
313DGtal::Mesh<TPoint>::getFaceBarycenter(Index i) const
314{
315 DGtal::Mesh<TPoint>::RealPoint c;
316 MeshFace aFace = getFace(i);
317 for ( auto &j: aFace){
318 TPoint p = getVertex(j);
319 for (typename TPoint::Dimension k = 0; k < TPoint::dimension; k++){
320 c[k] += static_cast<typename RealPoint::Component>(p[k]) ;
321 }
322 }
323 return c/static_cast<typename RealPoint::Component>(aFace.size());
324}
325
326
327template<typename TPoint>
328inline
329typename DGtal::Mesh<TPoint>::Size
330DGtal::Mesh<TPoint>::nbFaces() const
331{
332 return myFaceList.size();
333}
334
335template<typename TPoint>
336inline
337typename DGtal::Mesh<TPoint>::Size
338DGtal::Mesh<TPoint>::nbVertex() const
339{
340 return myVertexList.size();
341}
342
343
344template<typename TPoint>
345inline
346const DGtal::Color &
347DGtal::Mesh<TPoint>::getFaceColor(Index i) const
348{
349 if(mySaveFaceColor)
350 {
351 return myFaceColorList.at(i);
352 }
353 else
354 {
355 return myDefaultColor;
356 }
357}
358
359template <typename TPoint>
360struct MeshBoundingBoxCompPoints
361{
362MeshBoundingBoxCompPoints(typename TPoint::Dimension d): myDim(d){};
363bool operator() (const TPoint &p1, const TPoint &p2){return p1[myDim]<p2[myDim];};
364typename TPoint::Dimension myDim;
365};
366
367template<typename TPoint>
368inline
369std::pair<TPoint, TPoint>
370DGtal::Mesh<TPoint>::getBoundingBox() const
371{
372 std::pair<TPoint, TPoint> theResult;
373 TPoint lowerBound, upperBound;
374 for(unsigned int i=0; i< TPoint::size(); i++)
375 {
376 const MeshBoundingBoxCompPoints<TPoint> cmp_points(i);
377 upperBound[i] = (*(std::max_element(vertexBegin(), vertexEnd(), cmp_points)))[i];
378 lowerBound[i] = (*(std::min_element(vertexBegin(), vertexEnd(), cmp_points)))[i];
379 }
380 theResult.first = lowerBound ;
381 theResult.second = upperBound ;
382 return theResult;
383}
384
385
386template<typename TPoint>
387inline
388void
389DGtal::Mesh<TPoint>::setFaceColor(const Index index,
390 const DGtal::Color &aColor)
391{
392 if (!mySaveFaceColor)
393 {
394 for(unsigned int i = 0; i<myFaceList.size(); i++)
395 {
396 myFaceColorList.push_back(myDefaultColor);
397 }
398 mySaveFaceColor=true;
399 }
400 myFaceColorList.at(index) = aColor;
401}
402
403
404template<typename TPoint>
405inline
406bool
407DGtal::Mesh<TPoint>::isStoringFaceColors() const
408{
409 return mySaveFaceColor;
410}
411
412
413
414
415template<typename TPoint>
416inline
417void
418DGtal::Mesh<TPoint>::invertVertexFaceOrder(){
419 for(unsigned int i=0; i<myFaceList.size(); i++)
420 {
421 auto aFace = myFaceList.at(i);
422 for(unsigned int j=0; j < aFace.size()/2; j++)
423 {
424 unsigned int tmp=aFace.at(j);
425 aFace.at(j)=aFace.at(aFace.size()-1-j);
426 aFace.at(aFace.size()-1-j)=tmp;
427 }
428 }
429}
430
431template<typename TPoint>
432inline
433void
434DGtal::Mesh<TPoint>::clearFaces(){
435 myFaceList.clear();
436}
437
438template<typename TPoint>
439inline
440void
441DGtal::Mesh<TPoint>::changeScale(const typename TPoint::Component aScale){
442 for(typename VertexStorage::iterator it = vertexBegin(); it != vertexEnd(); it++)
443 {
444 (*it) *= aScale;
445 }
446}
447
448
449template<typename TPoint>
450inline
451double
452DGtal::Mesh<TPoint>::subDivideTriangularFaces(const double minArea){
453 double maxArea = 0;
454 std::vector<Mesh<TPoint>::MeshFace> facesToAdd;
455 for(unsigned int i =0; i< nbFaces(); i++)
456 {
457 typename Mesh<TPoint>::MeshFace aFace = getFace(i);
458 if(aFace.size()==3)
459 {
460 TPoint p1 = getVertex(aFace[0]);
461 TPoint p2 = getVertex(aFace[1]);
462 TPoint p3 = getVertex(aFace[2]);
463 TPoint c = (p1+p2+p3)/3.0;
464 double a = ((p2-p1).crossProduct(p3-p1)).norm()/2.0;
465 if (a>maxArea)
466 {
467 maxArea = a;
468 }
469
470 if(a>=minArea)
471 {
472 addVertex(c);
473 MeshFace f1, f2, f3;
474 f1.push_back(aFace[0]);
475 f1.push_back(aFace[1]);
476 f1.push_back(nbVertex()-1);
477 facesToAdd.push_back(f1);
478
479 f2.push_back(aFace[1]);
480 f2.push_back(aFace[2]);
481 f2.push_back(nbVertex()-1);
482 facesToAdd.push_back(f2);
483
484 f3.push_back(aFace[2]);
485 f3.push_back(aFace[0]);
486 f3.push_back(nbVertex()-1);
487 facesToAdd.push_back(f3);
488 }
489 else
490 {
491 facesToAdd.push_back(aFace);
492 }
493
494 }
495 }
496 clearFaces();
497 for(unsigned i=0; i<facesToAdd.size(); i++)
498 {
499 addFace(facesToAdd[i]);
500 }
501 return maxArea;
502}
503
504
505
506template<typename TPoint>
507inline
508unsigned int
509DGtal::Mesh<TPoint>::quadToTriangularFaces(){
510 unsigned int nbQuadT=0;
511 std::vector<Mesh<TPoint>::MeshFace> facesToAdd;
512 for(unsigned int i =0; i< nbFaces(); i++)
513 {
514 typename Mesh<TPoint>::MeshFace aFace = getFace(i);
515 if(aFace.size()==4)
516 {
517 MeshFace f1, f2;
518 f1.push_back(aFace[0]);
519 f1.push_back(aFace[1]);
520 f1.push_back(aFace[2]);
521 facesToAdd.push_back(f1);
522
523 f2.push_back(aFace[2]);
524 f2.push_back(aFace[3]);
525 f2.push_back(aFace[0]);
526 facesToAdd.push_back(f2);
527 nbQuadT++;
528 }
529 else
530 {
531 facesToAdd.push_back(aFace);
532 }
533 }
534 clearFaces();
535 for(unsigned i=0; i<facesToAdd.size(); i++)
536 {
537 addFace(facesToAdd[i]);
538 }
539 return nbQuadT;
540}
541
542
543
544
545//------------------------------------------------------------------------------
546template<typename TPoint>
547inline
548std::string
549DGtal::Mesh<TPoint>::className() const
550{
551 return "Mesh";
552}
553
554
555
556
557
558
559template <typename TPoint>
560inline
561void
562DGtal::Mesh<TPoint>::createTubularMesh(DGtal::Mesh<TPoint> &aMesh, const std::vector<TPoint> &aSkeleton,
563 const double aRadius,
564 const double angleStep, const DGtal::Color &aMeshColor)
565{
566 std::vector<double> aVecR;
567 aVecR.push_back(aRadius);
568 DGtal::Mesh<TPoint>::createTubularMesh(aMesh, aSkeleton,
569 aVecR, angleStep, aMeshColor);
570
571}
572
573
574template <typename TPoint>
575inline
576void
577DGtal::Mesh<TPoint>::createTubularMesh(DGtal::Mesh<TPoint> &aMesh, const std::vector<TPoint> &aSkeleton,
578 const std::vector<double> &aVectRadius,
579 const double angleStep, const DGtal::Color &aMeshColor)
580{
581 auto nbVertexInitial = aMesh.nbVertex();
582 ASSERT(aVectRadius.size() > 0);
583 // Generating vertices..
584 for(auto i = 0; i< (int)aSkeleton.size(); i++)
585 {
586 TPoint vectDir;
587 TPoint uDir1, uDirPrec;
588 TPoint uDir2;
589 TPoint firstPoint;
590
591 if(i != (int)aSkeleton.size()-1)
592 {
593 vectDir = aSkeleton.at(i+1) - aSkeleton.at(i);
594 }
595 else
596 {
597 vectDir = aSkeleton.at(i) - aSkeleton.at(i-1);
598 }
599
600 double d = -vectDir[0]* aSkeleton.at(i)[0] - vectDir[1]*aSkeleton.at(i)[1]
601 - vectDir[2]*aSkeleton.at(i)[2];
602 TPoint pRefOrigin;
603 if(vectDir[0]!=0)
604 {
605 pRefOrigin [0]= -d/vectDir[0];
606 pRefOrigin [1]= 0.0;
607 pRefOrigin [2]= 0.0;
608 if(aSkeleton.at(i) == pRefOrigin ||
609 (vectDir[1]==0 && vectDir[2]==0))
610 {
611 pRefOrigin[1]=-1.0;
612 }
613
614 }
615 else if (vectDir[1]!=0)
616 {
617 pRefOrigin [0]= 0.0;
618 pRefOrigin [1]= -d/vectDir[1];
619 pRefOrigin [2]= 0.0;
620 if(aSkeleton.at(i) == pRefOrigin ||
621 (vectDir[0]==0 && vectDir[2]==0))
622 {
623 pRefOrigin[0]=-1.0;
624 }
625 }else if (vectDir[2]!=0)
626 {
627 pRefOrigin [0]= 0.0;
628 pRefOrigin [1]= 0.0;
629 pRefOrigin [2]= -d/vectDir[2];
630 if(aSkeleton.at(i) == pRefOrigin ||
631 (vectDir[0]==0 && vectDir[1]==0))
632 {
633 pRefOrigin[0]=-1.0;
634 }
635 }
636 uDir1=(pRefOrigin-aSkeleton.at(i))/((pRefOrigin-aSkeleton.at(i)).norm());
637 uDir2[0] = uDir1[1]*vectDir[2]-uDir1[2]*vectDir[1];
638 uDir2[1] = uDir1[2]*vectDir[0]-uDir1[0]*vectDir[2];
639 uDir2[2] = uDir1[0]*vectDir[1]-uDir1[1]*vectDir[0];
640 uDir2/=uDir2.norm();
641 for(double a = 0.0; a < 2.0*M_PI; a += angleStep)
642 {
643 TPoint vMove = aVectRadius.at(i%aVectRadius.size())*(uDir1*cos(a) + uDir2*sin(a));
644 aMesh.addVertex(vMove + aSkeleton[i]);
645 if(a==0)
646 {
647 firstPoint = vMove + aSkeleton[i]+vectDir;
648 }
649
650 }
651 }
652 unsigned int nbPtPerFaces = static_cast<unsigned int>((aMesh.nbVertex()-nbVertexInitial)/aSkeleton.size());
653
654 // Generating faces...
655 for(auto i = 0; i< (int)aSkeleton.size()-1; i++)
656 {
657 if (aSkeleton.at(i)==aSkeleton.at(i+1)){
658 trace.warning() << "Two skeleton points are identical, ignoring one point." << std::endl;
659 continue;
660 }
661 // Computing best shift between two consecutive ring points to generate tube face.
662 // (criteria defined by the minimal distance between 4 sampling points)
663 double minDistance = std::numeric_limits<double>::max();
664 TPoint ptRefRing1 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces);
665 TPoint ptRefRing2 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+nbPtPerFaces/4);
666 TPoint ptRefRing3 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+2*(nbPtPerFaces/4));
667 TPoint ptRefRing4 = aMesh.getVertex(nbVertexInitial+i*nbPtPerFaces+3*(nbPtPerFaces/4));
668
669 unsigned int shift = 0;
670 TPoint vectDir;
671 if(i != (int)aSkeleton.size()-1)
672 {
673 vectDir = aSkeleton.at(i+1) - aSkeleton.at(i);
674 }
675 else
676 {
677 vectDir = aSkeleton.at(i) - aSkeleton.at(i-1);
678 }
679
680 for(unsigned int k=0; k<nbPtPerFaces; k++)
681 {
682 TPoint pScan1 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+k);
683 TPoint pScan2 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
684 (nbPtPerFaces/4+k)%nbPtPerFaces);
685 TPoint pScan3 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
686 (2*(nbPtPerFaces/4)+k)%nbPtPerFaces);
687 TPoint pScan4 = aMesh.getVertex(nbVertexInitial+(i+1)*nbPtPerFaces+
688 (3*(nbPtPerFaces/4)+k)%nbPtPerFaces);
689 double distance = (ptRefRing1 - pScan1).norm()+(ptRefRing2 - pScan2).norm()+
690 (ptRefRing3 - pScan3).norm()+(ptRefRing4 - pScan4).norm();
691 if(distance<minDistance){
692 shift = k;
693 minDistance = distance;
694 }
695
696 }
697 for(unsigned int k=0; k<nbPtPerFaces; k++)
698 {
699 Mesh<TPoint>::MeshFace aFace;
700 aMesh.addQuadFace(nbVertexInitial+k+i*nbPtPerFaces,
701 nbVertexInitial+(shift+k)%nbPtPerFaces+nbPtPerFaces*(i+1),
702 nbVertexInitial+(shift+k+1)%nbPtPerFaces+nbPtPerFaces*(i+1),
703 nbVertexInitial+(k+1)%nbPtPerFaces+i*nbPtPerFaces,
704 aMeshColor);
705 }
706 }
707}
708
709
710
711
712
713template <typename TPoint>
714template <typename TValue>
715inline
716void
717DGtal::Mesh<TPoint>::createMeshFromHeightSequence(Mesh<TPoint> &aMesh, const std::vector<TValue> & anValueSequence,
718 const unsigned int lengthSequence,
719 double stepX, double stepY, double stepZ,
720 const DGtal::Color &aMeshColor ){
721 const auto nbVertexInitial = aMesh.nbVertex();
722 // Generating vertices..
723 int i = 0;
724 unsigned int posY = 0;
725 while(i+(int)lengthSequence-1 < (int)anValueSequence.size()){
726 for(unsigned int j = 0; j < lengthSequence; j++, i++){
727 aMesh.addVertex(TPoint(j*stepX, posY*stepY, stepZ*anValueSequence.at(i)));
728 }
729 posY++;
730 }
731 // Generating faces...
732 i = 0;
733 posY = 0;
734 while(i+(int)lengthSequence-1 < (int)anValueSequence.size() - (int)lengthSequence){
735 for(auto j = 0; j < (int)lengthSequence-1; j++, i++){
736 aMesh.addQuadFace(nbVertexInitial+i, nbVertexInitial+i+1,
737 nbVertexInitial+i+1+lengthSequence,
738 nbVertexInitial+i+lengthSequence,
739 aMeshColor);
740 }
741 i++;
742 posY++;
743 }
744}
745
746
747
748
749template <typename TPoint>
750inline
751std::ostream&
752DGtal::operator<< ( std::ostream & out,
753 const Mesh<TPoint> & object )
754{
755 object.selfDisplay( out );
756 return out;
757}
758
759
760
761
762
763// //
764///////////////////////////////////////////////////////////////////////////////
765
766
Generated on Sun Nov 27 2022 15:15:50 for DGtal by doxygen 1.9.5