QGIS API Documentation  3.23.0-Master (dd0cd13a00)
qgstopologicalmesh.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstopologicalmesh.cpp - QgsTopologicalMesh
3 
4  ---------------------
5  begin : 18.6.2021
6  copyright : (C) 2021 by Vincent Cloarec
7  email : vcloarec at gmail dot com
8  ***************************************************************************
9  * *
10  * This program is free software; you can redistribute it and/or modify *
11  * it under the terms of the GNU General Public License as published by *
12  * the Free Software Foundation; either version 2 of the License, or *
13  * (at your option) any later version. *
14  * *
15  ***************************************************************************/
16 
17 #include "qgis.h"
18 #include "qgstopologicalmesh.h"
19 #include "qgsmesheditor.h"
20 #include "qgsmessagelog.h"
21 #include "qgsgeometryutils.h"
22 #include "qgslinestring.h"
23 
24 #include <poly2tri.h>
25 #include <QSet>
26 #include <QQueue>
27 
28 static int vertexPositionInFace( int vertexIndex, const QgsMeshFace &face )
29 {
30  return face.indexOf( vertexIndex );
31 }
32 
33 static int vertexPositionInFace( const QgsMesh &mesh, int vertexIndex, int faceIndex )
34 {
35  if ( faceIndex < 0 || faceIndex >= mesh.faceCount() )
36  return -1;
37 
38  return vertexPositionInFace( vertexIndex, mesh.face( faceIndex ) );
39 }
40 
41 static double crossProduct( int centralVertex, int vertex1, int vertex2, const QgsMesh &mesh )
42 {
43  QgsMeshVertex vc = mesh.vertices.at( centralVertex );
44  QgsMeshVertex v1 = mesh.vertices.at( vertex1 );
45  QgsMeshVertex v2 = mesh.vertices.at( vertex2 );
46 
47  double ux1 = v1.x() - vc.x();
48  double uy1 = v1.y() - vc.y();
49  double vx1 = v2.x() - vc.x();
50  double vy1 = v2.y() - vc.y();
51 
52  return ux1 * vy1 - uy1 * vx1;
53 }
54 
55 
56 QgsMeshVertexCirculator::QgsMeshVertexCirculator( const QgsTopologicalMesh &topologicalMesh, int vertexIndex )
57  : mFaces( topologicalMesh.mMesh->faces )
58  , mFacesNeighborhood( topologicalMesh.mFacesNeighborhood )
59  , mVertexIndex( vertexIndex )
60 {
61  if ( vertexIndex >= 0 && vertexIndex < topologicalMesh.mMesh->vertexCount() )
62  {
63  mCurrentFace = topologicalMesh.mVertexToFace[vertexIndex];
64  mIsValid = vertexPositionInFace( *topologicalMesh.mesh(), vertexIndex, mCurrentFace ) != -1;
65  }
66  else
67  {
68  mIsValid = false;
69  }
70 
71  if ( mIsValid )
72  mLastValidFace = mCurrentFace;
73 }
74 
75 QgsMeshVertexCirculator::QgsMeshVertexCirculator( const QgsTopologicalMesh::TopologicalFaces &topologicalFaces, int faceIndex, int vertexIndex )
76  : mFaces( topologicalFaces.mFaces )
77  , mFacesNeighborhood( topologicalFaces.mFacesNeighborhood )
78  , mVertexIndex( vertexIndex )
79 {
80  const QgsMeshFace &face = topologicalFaces.mFaces.at( faceIndex );
81  mIsValid = vertexPositionInFace( vertexIndex, face ) != -1;
82 
83  mCurrentFace = faceIndex;
84  mLastValidFace = mCurrentFace;
85 }
86 
88  : mFaces( topologicalFaces.mFaces )
89  , mFacesNeighborhood( topologicalFaces.mFacesNeighborhood )
90  , mVertexIndex( vertexIndex )
91 {
92  if ( topologicalFaces.mVerticesToFace.contains( vertexIndex ) )
93  mCurrentFace = topologicalFaces.mVerticesToFace.values( vertexIndex ).first();
94  mLastValidFace = mCurrentFace;
95  mIsValid = mCurrentFace != -1;
96 }
97 
99 {
100  if ( mCurrentFace == -1 )
101  mCurrentFace = mLastValidFace;
102  else
103  {
104  int currentPos = positionInCurrentFace();
105  Q_ASSERT( currentPos != -1 );
106 
107  const QgsMeshFace &currentFace = mFaces.at( mCurrentFace );
108  int faceSize = currentFace.size();
109  mLastValidFace = mCurrentFace;
110  mCurrentFace = mFacesNeighborhood.at( mCurrentFace ).at( ( currentPos + faceSize - 1 ) % currentFace.count() );
111  }
112 
113  return mCurrentFace;
114 }
115 
117 {
118  if ( mCurrentFace == -1 )
119  mCurrentFace = mLastValidFace;
120  else
121  {
122  int currentPos = positionInCurrentFace();
123  Q_ASSERT( currentPos != -1 );
124 
125  const QgsMeshFace &currentFace = mFaces.at( mCurrentFace );
126  int faceSize = currentFace.size();
127  mLastValidFace = mCurrentFace;
128  mCurrentFace = mFacesNeighborhood.at( mCurrentFace ).at( ( currentPos ) % faceSize );
129  }
130 
131  return mCurrentFace;
132 }
133 
135 {
136  return mCurrentFace;
137 }
138 
140 {
141  if ( mCurrentFace != -1 )
142  return mFaces.at( mCurrentFace );
143  else
144  return QgsMeshFace();
145 }
146 
148 {
149  if ( !isValid() )
150  return false;
151 
152  if ( mCurrentFace == -1 )
153  mCurrentFace = mLastValidFace;
154 
155  int firstFace = mCurrentFace;
156 
157  while ( turnClockwise() != -1 && currentFaceIndex() != firstFace ) {}
158  if ( mCurrentFace == firstFace )
159  return false; // a complete turn, so the vertex is not a boundary vertex, something went wrong
160 
161  return ( turnCounterClockwise() != -1 );
162 }
163 
165 {
166  if ( !isValid() )
167  return false;
168 
169  if ( mCurrentFace == -1 )
170  mCurrentFace = mLastValidFace;
171 
172  int firstFace = mCurrentFace;
173 
174  while ( turnCounterClockwise() != -1 && currentFaceIndex() != firstFace ) {}
175  if ( mCurrentFace == firstFace )
176  return false; // a complete turn, so the vertex is not a boundary vertex, something went wrong
177 
178  return ( turnClockwise() != -1 );
179 }
180 
182 {
183  if ( mCurrentFace == -1 )
184  return -1;
185 
186  const QgsMeshFace &face = currentFace();
187 
188  if ( face.isEmpty() )
189  return -1;
190 
191  int vertexPosition = vertexPositionInFace( mVertexIndex, currentFace() );
192 
193  if ( vertexPosition == -1 )
194  return -1;
195 
196  return face.at( ( vertexPosition + 1 ) % face.count() );
197 }
198 
200 {
201  if ( mCurrentFace == -1 )
202  return -1;
203 
204  const QgsMeshFace &face = currentFace();
205 
206  if ( face.isEmpty() )
207  return -1;
208 
209  int vertexPosition = vertexPositionInFace( mVertexIndex, currentFace() );
210 
211  if ( vertexPosition == -1 )
212  return -1;
213 
214  return face.at( ( vertexPosition - 1 + face.count() ) % face.count() );
215 }
216 
218 {
219  return mIsValid;
220 }
221 
223 {
224  QList<int> ret;
225  if ( !isValid() )
226  return ret;
227 
228  if ( mCurrentFace != -1 )
229  ret.append( mCurrentFace );
230  while ( turnCounterClockwise() != ret.first() && currentFaceIndex() != -1 )
231  ret.append( currentFaceIndex() );
232 
233 
234  if ( currentFaceIndex() == -1 ) //we encounter a boundary, restart with other direction
235  {
236  ret.clear();
237  if ( turnClockwise() == -1 )
238  return ret;
239  ret.append( currentFaceIndex() );
240  while ( turnClockwise() != -1 )
241  {
242  ret.append( currentFaceIndex() );
243  }
244  }
245 
246  return ret;
247 }
248 
250 {
251  if ( mDegree != -1 )
252  return mDegree;
253 
254  mDegree = 0;
255  if ( !mIsValid )
256  return mDegree;
257 
259  int firstFace = currentFaceIndex();
260 
261  mDegree = 2; //if the vertex is not free, the vertex is linked with at least 2 other vertex
262 
263  while ( turnClockwise() != firstFace && currentFaceIndex() != -1 )
264  ++mDegree;
265 
266  return mDegree;
267 }
268 
269 
270 int QgsMeshVertexCirculator::positionInCurrentFace() const
271 {
272  if ( mCurrentFace < 0 || mCurrentFace > mFaces.count() )
273  return -1;
274 
275  return vertexPositionInFace( mVertexIndex, mFaces.at( mCurrentFace ) );
276 }
277 
279 {
280  Changes changes;
281  changes.mFacesToAdd = topologicalFaces.mFaces;
282  changes.mAddedFacesFirstIndex = mMesh->faceCount();
283 
284  changes.mFacesNeighborhoodToAdd.resize( changes.mFacesToAdd.count() );
285  for ( int i = 0; i < changes.mFacesToAdd.count(); ++i )
286  changes.mFacesNeighborhoodToAdd[i] = QVector<int>( changes.mFacesToAdd.at( i ).count(), -1 );
287 
288  for ( int boundary : topologicalFaces.mBoundaries )
289  {
290  //if the boundary is a free vertex in the destination mesh, no need to check
291  if ( mVertexToFace.at( boundary ) == -1 )
292  continue;
293 
294  const QList<int> &linkedFaces = topologicalFaces.mVerticesToFace.values( boundary );
295  for ( int linkedFace : linkedFaces )
296  {
297  QgsMeshVertexCirculator newFacesCirculator( topologicalFaces, linkedFace, boundary );
298  //search for face boundary on clockwise side of new faces
299  newFacesCirculator.goBoundaryClockwise();
300  int oppositeVertexForNewFace = newFacesCirculator.oppositeVertexClockwise();
301  if ( mVertexToFace.at( oppositeVertexForNewFace ) == -1 )
302  continue;
303 
304  QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
305  meshCirculator.goBoundaryCounterClockwise();
306  int oppositeVertexForMeshFace = meshCirculator.oppositeVertexCounterClockwise();
307 
308  const QgsMeshFace &newFaceBoundary = newFacesCirculator.currentFace();
309  int boundaryPositionInNewFace = vertexPositionInFace( boundary, newFaceBoundary );
310 
311  if ( oppositeVertexForMeshFace != oppositeVertexForNewFace )
312  {
313  changes.mFacesNeighborhoodToAdd[newFacesCirculator.currentFaceIndex()][boundaryPositionInNewFace] = -1 ;
314  }
315  else
316  {
317  const QgsMeshFace &meshFaceBoundary = meshCirculator.currentFace();
318  int boundaryPositionInMeshFace = vertexPositionInFace( meshCirculator.oppositeVertexCounterClockwise(), meshFaceBoundary );
319 
320  changes.mNeighborhoodChanges.append( std::array<int, 4>(
321  {
322  meshCirculator.currentFaceIndex(),
323  boundaryPositionInMeshFace,
324  -1,
325  changes.addedFaceIndexInMesh( newFacesCirculator.currentFaceIndex() )
326  } ) );
327 
328  changes.mFacesNeighborhoodToAdd[newFacesCirculator.currentFaceIndex()][boundaryPositionInNewFace] = meshCirculator.currentFaceIndex();
329  }
330  }
331  }
332 
333  for ( int f = 0; f < changes.mFacesToAdd.count(); ++f )
334  for ( int n = 0; n < changes.mFacesToAdd.at( f ).count(); ++n )
335  if ( changes.mFacesNeighborhoodToAdd.at( f ).at( n ) == -1 )
336  changes.mFacesNeighborhoodToAdd[f][n] = changes.addedFaceIndexInMesh( topologicalFaces.mFacesNeighborhood.at( f ).at( n ) );
337 
338  const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
339  for ( const int vtc : verticesToFaceToChange )
340  if ( mVertexToFace.at( vtc ) == -1 )
341  changes.mVerticesToFaceChanges.append( {vtc,
342  mVertexToFace.at( vtc ),
343  changes.addedFaceIndexInMesh( topologicalFaces.mVerticesToFace.values( vtc ).first() ) } );
344 
345  applyChanges( changes );
346 
347  return changes;
348 }
349 
351 {
352  int initialVerticesCount = mMesh->vertices.count();
353  if ( !changes.mVerticesToAdd.empty() )
354  {
355  int newSize = mMesh->vertices.count() + changes.mVerticesToAdd.count();
356  mMesh->vertices.resize( newSize );
357  mVertexToFace.resize( newSize );
358  }
359 
360  if ( !changes.mFacesToAdd.empty() )
361  {
362  int newSize = mMesh->faceCount() + changes.mFacesToAdd.count();
363  mMesh->faces.resize( newSize );
364  mFacesNeighborhood.resize( newSize );
365  }
366 
367  for ( int i = 0; i < changes.mFacesToRemove.count(); ++i )
368  {
369  mMesh->faces[changes.removedFaceIndexInMesh( i )] = QgsMeshFace();
370  mFacesNeighborhood[changes.removedFaceIndexInMesh( i )] = FaceNeighbors();//changes.facesNeighborhoodToRemove[i];
371  }
372 
373  for ( int i = 0; i < changes.mVerticesToRemoveIndexes.count(); ++i )
374  {
375  int vertexIndex = changes.mVerticesToRemoveIndexes.at( i );
376  if ( mVertexToFace.at( vertexIndex ) == -1 )
377  dereferenceAsFreeVertex( vertexIndex );
378  mMesh->vertices[vertexIndex] = QgsMeshVertex();
379  mVertexToFace[vertexIndex] = -1;
380  }
381 
382  for ( int i = 0; i < changes.mVerticesToAdd.count(); ++i )
383  {
384  mMesh->vertices[initialVerticesCount + i] = changes.mVerticesToAdd.at( i );
385  mVertexToFace[initialVerticesCount + i] = changes.mVertexToFaceToAdd.at( i );
386  if ( changes.mVertexToFaceToAdd.at( i ) == -1 )
387  referenceAsFreeVertex( initialVerticesCount + i );
388  }
389 
390  for ( int i = 0; i < changes.mFacesToAdd.count(); ++i )
391  {
392  mMesh->faces[changes.addedFaceIndexInMesh( i )] = changes.mFacesToAdd.at( i );
393  mFacesNeighborhood[changes.addedFaceIndexInMesh( i )] = changes.mFacesNeighborhoodToAdd.at( i );
394  }
395 
396  for ( const std::array<int, 4> neigborChange : std::as_const( changes.mNeighborhoodChanges ) )
397  {
398  const int faceIndex = neigborChange.at( 0 );
399  const int positionInFace = neigborChange.at( 1 );
400  const int valueToApply = neigborChange.at( 3 );
401  mFacesNeighborhood[faceIndex][positionInFace] = valueToApply;
402  }
403 
404  for ( const std::array<int, 3> vertexToFaceChange : std::as_const( changes.mVerticesToFaceChanges ) )
405  {
406  int vertexIndex = vertexToFaceChange.at( 0 );
407  mVertexToFace[vertexToFaceChange.at( 0 )] = vertexToFaceChange.at( 2 );
408 
409  if ( vertexToFaceChange.at( 2 ) == -1 &&
410  vertexToFaceChange.at( 1 ) != -1 &&
411  !mMesh->vertices.at( vertexIndex ).isEmpty() )
412  referenceAsFreeVertex( vertexIndex );
413 
414  if ( vertexToFaceChange.at( 1 ) == -1 && vertexToFaceChange.at( 2 ) != -1 )
415  dereferenceAsFreeVertex( vertexIndex );
416  }
417 
418  for ( int i = 0; i < changes.mChangeCoordinateVerticesIndexes.count(); ++i )
419  {
420  int vertexIndex = changes.mChangeCoordinateVerticesIndexes.at( i );
421  if ( !changes.mNewZValues.isEmpty() )
422  mMesh->vertices[vertexIndex].setZ( changes.mNewZValues.at( i ) );
423  if ( !changes.mNewXYValues.isEmpty() )
424  {
425  const QgsPointXY &pt = changes.mNewXYValues.at( i );
426  mMesh->vertices[vertexIndex].setX( pt.x() );
427  mMesh->vertices[vertexIndex].setY( pt.y() );
428  }
429  }
430 }
431 
433 {
434  for ( const std::array<int, 4> neigborChange : std::as_const( changes.mNeighborhoodChanges ) )
435  {
436  const int faceIndex = neigborChange.at( 0 );
437  const int positionInFace = neigborChange.at( 1 );
438  const int valueToApply = neigborChange.at( 2 );
439  mFacesNeighborhood[faceIndex][positionInFace] = valueToApply;
440  }
441 
442  for ( int i = 0; i < changes.mFacesToRemove.count(); ++i )
443  {
444  mMesh->faces[changes.removedFaceIndexInMesh( i )] = changes.mFacesToRemove.at( i );
445  mFacesNeighborhood[changes.removedFaceIndexInMesh( i )] = changes.mFacesNeighborhoodToRemove.at( i );
446  }
447 
448  for ( int i = 0; i < changes.mVerticesToRemoveIndexes.count(); ++i )
449  {
450  int vertexIndex = changes.mVerticesToRemoveIndexes.at( i );
451  mMesh->vertices[vertexIndex] = changes.mRemovedVertices.at( i );
452  mVertexToFace[vertexIndex] = changes.mVerticesToFaceRemoved.at( i );
453  if ( mVertexToFace.at( vertexIndex ) == -1 )
454  referenceAsFreeVertex( vertexIndex );
455  }
456 
457  int verticesToFaceChangesCount = changes.mVerticesToFaceChanges.count();
458  for ( int i = 0; i < verticesToFaceChangesCount; ++i )
459  {
460  const std::array<int, 3> vertexToFaceChange = changes.mVerticesToFaceChanges.at( verticesToFaceChangesCount - i - 1 );
461  int vertexIndex = vertexToFaceChange.at( 0 );
462  mVertexToFace[vertexIndex] = vertexToFaceChange.at( 1 );
463 
464  if ( vertexToFaceChange.at( 2 ) == -1 && vertexToFaceChange.at( 1 ) != -1 )
465  dereferenceAsFreeVertex( vertexIndex );
466 
467  if ( vertexToFaceChange.at( 1 ) == -1 &&
468  vertexToFaceChange.at( 2 ) != -1 &&
469  !mMesh->vertex( vertexIndex ).isEmpty() )
470  referenceAsFreeVertex( vertexIndex );
471  }
472 
473  if ( !changes.mFacesToAdd.empty() )
474  {
475  int newSize = mMesh->faceCount() - changes.mFacesToAdd.count();
476  mMesh->faces.resize( newSize );
477  mFacesNeighborhood.resize( newSize );
478  }
479 
480  if ( !changes.mVerticesToAdd.isEmpty() )
481  {
482  int newSize = mMesh->vertexCount() - changes.mVerticesToAdd.count();
483 
484  for ( int i = newSize; i < mMesh->vertexCount(); ++i )
485  if ( mVertexToFace.at( i ) == -1 )
486  dereferenceAsFreeVertex( i );
487 
488  mMesh->vertices.resize( newSize );
489  mVertexToFace.resize( newSize );
490  }
491 
492  for ( int i = 0; i < changes.mChangeCoordinateVerticesIndexes.count(); ++i )
493  {
494  int vertexIndex = changes.mChangeCoordinateVerticesIndexes.at( i );
495  if ( !changes.mOldZValues.isEmpty() )
496  mMesh->vertices[vertexIndex].setZ( changes.mOldZValues.at( i ) );
497  if ( !changes.mOldXYValues.isEmpty() )
498  {
499  const QgsPointXY &pt = changes.mOldXYValues.at( i );
500  mMesh->vertices[vertexIndex].setX( pt.x() );
501  mMesh->vertices[vertexIndex].setY( pt.y() );
502  }
503  }
504 }
505 
507 {
508  return QgsMeshVertexCirculator( *this, vertexIndex );
509 }
510 
511 QSet<int> QgsTopologicalMesh::concernedFacesBy( const QList<int> faceIndexes ) const
512 {
513  QSet<int> faces;
514  for ( const int faceIndex : faceIndexes )
515  {
516  const QgsMeshFace &face = mMesh->face( faceIndex );
517  for ( int i = 0; i < face.count(); ++i )
518  faces.unite( qgis::listToSet( facesAroundVertex( face.at( i ) ) ) );
519  }
520  return faces;
521 }
522 
523 void QgsTopologicalMesh::dereferenceAsFreeVertex( int vertexIndex )
524 {
525  mFreeVertices.remove( vertexIndex );
526 }
527 
528 void QgsTopologicalMesh::referenceAsFreeVertex( int vertexIndex )
529 {
530  // QSet used to retrieve free vertex without going through all the vertices container.
531  // But maybe later we could use more sophisticated reference (spatial index?), to retrieve free vertex in an extent
532  mFreeVertices.insert( vertexIndex );
533 }
534 
536 {
537  for ( int faceIndex = 0 ; faceIndex < mMesh->faces.count( ); ++faceIndex )
538  {
539  const QgsMeshFace &face = mMesh->faces.at( faceIndex );
540  const FaceNeighbors &neighborhood = mFacesNeighborhood.at( faceIndex );
541  if ( face.count() != neighborhood.count() )
543  for ( int i = 0; i < face.count(); ++i )
544  {
545  int vertexIndex = face.at( i );
546  // check if each vertices is linked to a face (not free vertex)
547  if ( mVertexToFace.at( vertexIndex ) == -1 )
549 
550  int neighborIndex = neighborhood.at( i );
551  if ( neighborIndex != -1 )
552  {
553  const QgsMeshFace &neighborFace = mMesh->faces.at( neighborIndex );
554  if ( neighborFace.isEmpty() )
556  int neighborSize = neighborFace.size();
557  const FaceNeighbors &neighborhoodOfNeighbor = mFacesNeighborhood.at( neighborIndex );
558  int posInNeighbor = vertexPositionInFace( *mMesh, vertexIndex, neighborIndex );
559  if ( neighborhoodOfNeighbor.isEmpty() || neighborhoodOfNeighbor.at( ( posInNeighbor + neighborSize - 1 ) % neighborSize ) != faceIndex )
561  }
562  }
563  }
564 
565  for ( int vertexIndex = 0; vertexIndex < mMesh->vertexCount(); ++vertexIndex )
566  {
567  if ( mVertexToFace.at( vertexIndex ) != -1 )
568  {
569  if ( !mMesh->face( mVertexToFace.at( vertexIndex ) ).contains( vertexIndex ) )
571 
572  if ( facesAroundVertex( vertexIndex ).count() == 0 )
574  }
575  }
576 
577  return QgsMeshEditingError();
578 }
579 
580 QgsMeshEditingError QgsTopologicalMesh::checkTopology( const QgsMesh &mesh, int maxVerticesPerFace )
581 {
582  QgsMesh temp = mesh;
583  QgsMeshEditingError error;
584  createTopologicalMesh( &temp, maxVerticesPerFace, error );
585  return error;
586 }
587 
589 {
590  return mMesh;
591 }
592 
593 int QgsTopologicalMesh::firstFaceLinked( int vertexIndex ) const
594 {
595  if ( vertexIndex < 0 || vertexIndex >= mMesh->vertexCount() )
596  return -1;
597  return mVertexToFace.at( vertexIndex );
598 }
599 
600 bool QgsTopologicalMesh::isVertexOnBoundary( int vertexIndex ) const
601 {
602  QgsMeshVertexCirculator circulator = vertexCirculator( vertexIndex );
603 
604  if ( circulator.isValid() )
605  return circulator.goBoundaryClockwise();
606 
607  return false;
608 }
609 
610 bool QgsTopologicalMesh::isVertexFree( int vertexIndex ) const
611 {
612  if ( vertexIndex < 0 || vertexIndex >= mMesh->vertexCount() )
613  return false;
614 
615  if ( mMesh->vertices.at( vertexIndex ).isEmpty() )
616  return false;
617 
618  return mVertexToFace.at( vertexIndex ) == -1;
619 }
620 
622 {
623 #if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
624  return mFreeVertices.values();
625 #else
626  return QList<int>( mFreeVertices.begin(), mFreeVertices.end() );
627 #endif
628 }
629 
631 {
632  // First check if the face is convex and put it counter clockwise
633  // If the index are not well ordered (edges intersect), invalid face --> return false
634  int faceSize = face.count();
635  if ( faceSize < 3 )
637 
638  int direction = 0;
639  for ( int i = 0; i < faceSize; ++i )
640  {
641  int iv0 = face[i];
642  int iv1 = face[( i + 1 ) % faceSize];
643  int iv2 = face[( i + 2 ) % faceSize];
644 
645  if ( iv0 < 0 || iv0 >= mesh->vertexCount() )
647 
648  if ( iv1 < 0 || iv1 >= mesh->vertexCount() )
650 
651  if ( iv2 < 0 || iv2 >= mesh->vertexCount() )
653 
654  const QgsMeshVertex &v0 = mesh->vertices.at( iv0 ) ;
655  const QgsMeshVertex &v1 = mesh->vertices.at( iv1 ) ;
656  const QgsMeshVertex &v2 = mesh->vertices.at( iv2 ) ;
657 
658  if ( v0.isEmpty() )
660 
661  if ( v1.isEmpty() )
663 
664  if ( v2.isEmpty() )
666 
667  double crossProd = crossProduct( iv1, iv0, iv2, *mesh ); //if cross product>0, we have two edges clockwise
668  if ( direction != 0 && crossProd * direction < 0 ) // We have a convex face or a (partially) flat face
670  else if ( crossProd == 0 )
672  else if ( direction == 0 && crossProd != 0 )
673  direction = crossProd / std::fabs( crossProd );
674  }
675 
676  if ( direction > 0 )// clockwise --> reverse the order of the index;
677  {
678  for ( int i = 0; i < faceSize / 2; ++i )
679  {
680  int temp = face[i];
681  face[i] = face.at( faceSize - i - 1 );
682  face[faceSize - i - 1] = temp;
683  }
684  }
685 
687 }
688 
690 {
691  QVector<int> oldToNewIndex( mMesh->vertices.count(), -1 );
692  int verticesTotalCount = mMesh->vertexCount();
693  int oldIndex = 0;
694  int newIndex = 0;
695  while ( oldIndex < verticesTotalCount )
696  {
697  if ( mMesh->vertex( oldIndex ).isEmpty() )
698  {
699  oldIndex++;
700  }
701  else
702  {
703  oldToNewIndex[oldIndex] = newIndex;
704  if ( oldIndex != newIndex )
705  mMesh->vertices[newIndex] = mMesh->vertices[oldIndex];
706  oldToNewIndex[oldIndex] = newIndex;
707  oldIndex++;
708  newIndex++;
709  }
710  }
711 
712  mMesh->vertices.resize( newIndex );
713 
714  oldIndex = 0;
715  newIndex = 0;
716  int facesTotalCount = mMesh->faceCount();
717  while ( oldIndex < facesTotalCount )
718  {
719  if ( mMesh->face( oldIndex ).isEmpty() )
720  oldIndex++;
721  else
722  {
723  if ( oldIndex != newIndex )
724  mMesh->faces[newIndex] = mMesh->faces[oldIndex];
725  QgsMeshFace &face = mMesh->faces[newIndex];
726  for ( int i = 0; i < face.count(); ++i )
727  face[i] = oldToNewIndex[face.at( i )];
728  newIndex++;
729  oldIndex++;
730  }
731  }
732 
733  mMesh->faces.resize( newIndex );
734 
735  mVertexToFace.clear();
736  mFacesNeighborhood.clear();
737 }
738 
740 {
741  QVector<int> oldToNewVerticesIndexes;
742  if ( !renumberVertices( oldToNewVerticesIndexes ) )
743  return false;
744 
745 
746  QVector<int> oldToNewFacesIndexes;
747 
748  if ( !renumberFaces( oldToNewFacesIndexes ) )
749  return false;
750 
751  // If we are here, we can apply the renumbering
752 
753  QVector<QgsMeshVertex> tempVertices( mMesh->vertices.count() );
754  for ( int i = 0; i < oldToNewVerticesIndexes.count(); ++i )
755  {
756  tempVertices[oldToNewVerticesIndexes.at( i )] = mMesh->vertex( i );
757  }
758  mMesh->vertices = tempVertices;
759 
760  QVector<QgsMeshFace> tempFaces( mMesh->faces.count() );
761  for ( int i = 0; i < oldToNewFacesIndexes.count(); ++i )
762  {
763  tempFaces[oldToNewFacesIndexes.at( i )] = mMesh->face( i );
764 
765  QgsMeshFace &face = tempFaces[oldToNewFacesIndexes.at( i )];
766 
767  for ( int fi = 0; fi < face.count(); ++fi )
768  {
769  face[fi] = oldToNewVerticesIndexes.at( face.at( fi ) );
770  }
771  }
772 
773  mMesh->faces = tempFaces;
774 
775  return true;
776 
777 }
778 
779 bool QgsTopologicalMesh::renumberVertices( QVector<int> &oldToNewIndex ) const
780 {
781  std::vector<QgsMeshVertexCirculator> circulators;
782  circulators.reserve( mMesh->vertices.count() );
783  int minDegree = std::numeric_limits<int>::max();
784  int minDegreeVertex = -1;
785 
786  QQueue<int> queue;
787  QSet<int> nonThreadedVertex;
788  oldToNewIndex = QVector<int> ( mMesh->vertexCount(), -1 );
789  for ( int i = 0; i < mMesh->vertexCount(); ++i )
790  {
791  circulators.emplace_back( *this, i );
792  const QgsMeshVertexCirculator &circulator = circulators.back();
793  if ( circulators.back().degree() < minDegree )
794  {
795  minDegreeVertex = circulators.size() - 1;
796  minDegree = circulator.degree();
797  }
798  nonThreadedVertex.insert( i );
799  }
800 
801  auto sortedNeighbor = [ = ]( QList<int> &neighbors, int index )
802  {
803  const QgsMeshVertexCirculator &circ = circulators.at( index );
804 
805  if ( !circ.isValid() )
806  return;
807 
809  neighbors.append( circ.oppositeVertexCounterClockwise() );
810 
811  int firsrFace = circ.currentFaceIndex();
812  do
813  {
814  int neighborIndex = circ.oppositeVertexClockwise();
815  int degree = circulators.at( neighborIndex ).degree();
816  QList<int>::Iterator it = neighbors.begin();
817  while ( it != neighbors.end() )
818  {
819  if ( degree <= circulators.at( *it ).degree() )
820  {
821  neighbors.insert( it, neighborIndex );
822  break;
823  }
824  it++;
825  }
826  if ( it == neighbors.end() )
827  neighbors.append( neighborIndex );
828  }
829  while ( circ.turnClockwise() != firsrFace && circ.currentFaceIndex() != -1 );
830  };
831 
832  int newIndex = 0;
833  int currentVertex = minDegreeVertex;
834  nonThreadedVertex.remove( minDegreeVertex );
835 
836  while ( newIndex < mMesh->vertexCount() )
837  {
838  if ( oldToNewIndex[currentVertex] == -1 )
839  oldToNewIndex[currentVertex] = newIndex++;
840 
841  if ( circulators.at( currentVertex ).isValid() )
842  {
843  QList<int> neighbors;
844  sortedNeighbor( neighbors, currentVertex );
845 
846  for ( const int i : std::as_const( neighbors ) )
847  if ( oldToNewIndex.at( i ) == -1 && nonThreadedVertex.contains( i ) )
848  {
849  queue.enqueue( i );
850  nonThreadedVertex.remove( i );
851  }
852  }
853 
854  if ( queue.isEmpty() )
855  {
856  if ( nonThreadedVertex.isEmpty() && newIndex < mMesh->vertexCount() )
857  return false;
858 
859  const QList<int> remainingVertex = qgis::setToList( nonThreadedVertex );
860  int minRemainingDegree = std::numeric_limits<int>::max();
861  int minRemainingVertex = -1;
862  for ( const int i : remainingVertex )
863  {
864  int degree = circulators.at( i ).degree();
865  if ( degree < minRemainingDegree )
866  {
867  minRemainingDegree = degree;
868  minRemainingVertex = i;
869  }
870  }
871  currentVertex = minRemainingVertex;
872  nonThreadedVertex.remove( currentVertex );
873  }
874  else
875  {
876  currentVertex = queue.dequeue();
877  }
878  }
879 
880  return true;
881 }
882 
883 bool QgsTopologicalMesh::renumberFaces( QVector<int> &oldToNewIndex ) const
884 {
885  QQueue<int> queue;
886  QSet<int> nonThreadedFaces;
887 
888  oldToNewIndex = QVector<int>( mMesh->faceCount(), -1 );
889 
890  QVector<int> faceDegrees( mMesh->faceCount(), 0 );
891 
892  int minDegree = std::numeric_limits<int>::max();
893  int minDegreeFace = -1;
894  for ( int faceIndex = 0; faceIndex < mMesh->faceCount(); ++faceIndex )
895  {
896  const FaceNeighbors &neighbors = mFacesNeighborhood.at( faceIndex );
897 
898  int degree = 0;
899  for ( int n = 0; n < neighbors.size(); ++n )
900  {
901  if ( neighbors.at( n ) != -1 )
902  degree++;
903  }
904 
905  if ( degree < minDegree )
906  {
907  minDegree = degree;
908  minDegreeFace = faceIndex;
909  }
910 
911  faceDegrees[faceIndex] = degree;
912  nonThreadedFaces.insert( faceIndex );
913  }
914 
915  int newIndex = 0;
916  int currentFace = minDegreeFace;
917  nonThreadedFaces.remove( minDegreeFace );
918 
919  auto sortedNeighbor = [ = ]( QList<int> &neighbors, int index )
920  {
921  const FaceNeighbors &neighborhood = mFacesNeighborhood.at( index );
922 
923  for ( int i = 0; i < neighborhood.count(); ++i )
924  {
925  int neighborIndex = neighborhood.at( i );
926  if ( neighborIndex == -1 )
927  continue;
928 
929  int degree = faceDegrees.at( neighborIndex );
930  if ( neighbors.isEmpty() )
931  neighbors.append( neighborIndex );
932  else
933  {
934  QList<int>::Iterator it = neighbors.begin();
935  while ( it != neighbors.end() )
936  {
937  if ( degree <= faceDegrees.at( *it ) )
938  {
939  neighbors.insert( it, neighborIndex );
940  break;
941  }
942  it++;
943  }
944  if ( it == neighbors.end() )
945  neighbors.append( neighborIndex );
946  }
947  }
948  };
949 
950  while ( newIndex < mMesh->faceCount() )
951  {
952  if ( oldToNewIndex[currentFace] == -1 )
953  oldToNewIndex[currentFace] = newIndex++;
954 
955  QList<int> neighbors;
956  sortedNeighbor( neighbors, currentFace );
957 
958  for ( const int i : std::as_const( neighbors ) )
959  if ( oldToNewIndex.at( i ) == -1 && nonThreadedFaces.contains( i ) )
960  {
961  queue.enqueue( i );
962  nonThreadedFaces.remove( i );
963  }
964 
965  if ( queue.isEmpty() )
966  {
967  if ( nonThreadedFaces.isEmpty() && newIndex < mMesh->faceCount() )
968  return false;
969 
970  const QList<int> remainingFace = qgis::setToList( nonThreadedFaces );
971  int minRemainingDegree = std::numeric_limits<int>::max();
972  int minRemainingFace = -1;
973  for ( const int i : remainingFace )
974  {
975  int degree = faceDegrees.at( i );
976  if ( degree < minRemainingDegree )
977  {
978  minRemainingDegree = degree;
979  minRemainingFace = i;
980  }
981  }
982  currentFace = minRemainingFace;
983  nonThreadedFaces.remove( currentFace );
984  }
985  else
986  {
987  currentFace = queue.dequeue();
988  }
989 
990  }
991 
992  return true;
993 }
994 
995 QVector<QgsMeshFace> QgsTopologicalMesh::Changes::addedFaces() const
996 {
997  return mFacesToAdd;
998 }
999 
1000 QVector<QgsMeshFace> QgsTopologicalMesh::Changes::removedFaces() const
1001 {
1002  return mFacesToRemove;
1003 }
1004 
1006 {
1007  return mFaceIndexesToRemove;
1008 }
1009 
1010 QVector<QgsMeshVertex> QgsTopologicalMesh::Changes::addedVertices() const
1011 {
1012  return mVerticesToAdd;
1013 }
1014 
1016 {
1017  return mChangeCoordinateVerticesIndexes;
1018 }
1019 
1021 {
1022  return mNewZValues;
1023 }
1024 
1026 {
1027  return mNewXYValues;
1028 }
1029 
1031 {
1032  return mOldXYValues;
1033 }
1034 
1036 {
1037  return mNativeFacesIndexesGeometryChanged;
1038 }
1039 
1041 {
1042  return ( mFaceIndexesToRemove.isEmpty() &&
1043  mFacesToAdd.isEmpty() &&
1044  mFacesNeighborhoodToAdd.isEmpty() &&
1045  mFacesToRemove.isEmpty() &&
1046  mFacesNeighborhoodToRemove.isEmpty() &&
1047  mNeighborhoodChanges.isEmpty() &&
1048  mVerticesToAdd.isEmpty() &&
1049  mVertexToFaceToAdd.isEmpty() &&
1050  mVerticesToRemoveIndexes.isEmpty() &&
1051  mRemovedVertices.isEmpty() &&
1052  mVerticesToFaceRemoved.isEmpty() &&
1053  mVerticesToFaceChanges.isEmpty() &&
1054  mChangeCoordinateVerticesIndexes.isEmpty() &&
1055  mNewZValues.isEmpty() &&
1056  mOldZValues.isEmpty() &&
1057  mNewXYValues.isEmpty() &&
1058  mOldXYValues.isEmpty() &&
1059  mNativeFacesIndexesGeometryChanged.isEmpty() );
1060 }
1061 
1063 {
1064  return mVerticesToRemoveIndexes;
1065 }
1066 
1067 int QgsTopologicalMesh::Changes::addedFaceIndexInMesh( int internalIndex ) const
1068 {
1069  if ( internalIndex == -1 )
1070  return -1;
1071 
1072  return internalIndex + mAddedFacesFirstIndex;
1073 }
1074 
1075 int QgsTopologicalMesh::Changes::removedFaceIndexInMesh( int internalIndex ) const
1076 {
1077  if ( internalIndex == -1 )
1078  return -1;
1079 
1080  return mFaceIndexesToRemove.at( internalIndex );
1081 }
1082 
1084 {
1085  mAddedFacesFirstIndex = 0;
1086  mFaceIndexesToRemove.clear();
1087  mFacesToAdd.clear();
1088  mFacesNeighborhoodToAdd.clear();
1089  mFacesToRemove.clear();
1090  mFacesNeighborhoodToRemove.clear();
1091  mNeighborhoodChanges.clear();
1092 
1093  mVerticesToAdd.clear();
1094  mVertexToFaceToAdd.clear();
1095  mVerticesToRemoveIndexes.clear();
1096  mRemovedVertices.clear();
1097  mVerticesToFaceRemoved.clear();
1098  mVerticesToFaceChanges.clear();
1099 
1100  mChangeCoordinateVerticesIndexes.clear();
1101  mNewZValues.clear();
1102  mOldZValues.clear();
1103  mNewXYValues.clear();
1104  mOldXYValues.clear();
1105  mNativeFacesIndexesGeometryChanged.clear();
1106 }
1107 
1109 {
1110  Changes changes;
1111  changes.mVerticesToAdd.append( vertex );
1112  changes.mVertexToFaceToAdd.append( -1 );
1113 
1114  mMesh->vertices.append( vertex );
1115  mVertexToFace.append( -1 );
1116  referenceAsFreeVertex( mMesh->vertices.count() - 1 );
1117 
1118  return changes;
1119 }
1120 
1121 // Returns the orientation of the polygon formed by mesh vertices, <0 counter clockwise; >0 clockwise
1122 static double vertexPolygonOrientation( const QgsMesh &mesh, const QList<int> &vertexIndexes )
1123 {
1124  if ( vertexIndexes.count() < 3 )
1125  return 0.0;
1126 
1127  int hullDomainVertexPos = -1;
1128  double xMin = std::numeric_limits<double>::max();
1129  double yMin = std::numeric_limits<double>::max();
1130  for ( int i = 0; i < vertexIndexes.count(); ++i )
1131  {
1132  const QgsMeshVertex &vertex = mesh.vertices.at( vertexIndexes.at( i ) );
1133  if ( xMin >= vertex.x() && yMin > vertex.y() )
1134  {
1135  hullDomainVertexPos = i;
1136  xMin = vertex.x();
1137  yMin = vertex.y();
1138  }
1139  }
1140 
1141  if ( hullDomainVertexPos >= 0 )
1142  {
1143  int iv1 = vertexIndexes.at( ( hullDomainVertexPos - 1 + vertexIndexes.count() ) % vertexIndexes.count() );
1144  int iv2 = vertexIndexes.at( ( hullDomainVertexPos + 1 ) % vertexIndexes.count() );
1145  int ivc = vertexIndexes.at( ( hullDomainVertexPos ) );
1146  double cp = crossProduct( ivc, iv1, iv2, mesh );
1147  return cp;
1148  }
1149 
1150  return 0.0;
1151 }
1152 
1154 {
1155  if ( vertexIndex >= mVertexToFace.count() )
1156  return Changes();
1157 
1158  if ( mVertexToFace.at( vertexIndex ) == -1 ) //it is a free vertex
1159  {
1160  Changes changes;
1161  changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
1162  changes.mVerticesToRemoveIndexes.append( vertexIndex );
1163  changes.mVerticesToFaceRemoved.append( -1 );
1164  dereferenceAsFreeVertex( vertexIndex );
1165  mMesh->vertices[vertexIndex] = QgsMeshVertex();
1166  return changes;
1167  }
1168 
1169  //search concerned faces
1170  QgsMeshVertexCirculator circulator = vertexCirculator( vertexIndex );
1171  circulator.goBoundaryClockwise();
1172  QList<int> boundariesVertexIndex;
1173  QList<int> associateFaceToBoundaries;
1174  QList<int> removedFacesIndexes;
1175  QSet<int> boundaryInGlobalMesh;
1176 
1177  do
1178  {
1179  removedFacesIndexes.append( circulator.currentFaceIndex() );
1180  boundariesVertexIndex.append( circulator.oppositeVertexClockwise() );
1181  Q_ASSERT( !mMesh->vertices.at( boundariesVertexIndex.last() ).isEmpty() );
1182  const QgsMeshFace &currentFace = circulator.currentFace();
1183  associateFaceToBoundaries.append( mFacesNeighborhood.at( circulator.currentFaceIndex() ).at(
1184  vertexPositionInFace( boundariesVertexIndex.last(), currentFace ) ) );
1185 
1186  if ( currentFace.count() > 3 ) // quad or more, need other vertices
1187  {
1188  int posInface = vertexPositionInFace( vertexIndex, currentFace );
1189  for ( int i = 2; i < currentFace.count() - 1; ++i )
1190  {
1191  boundariesVertexIndex.append( currentFace.at( ( posInface + i ) % currentFace.count() ) );
1192  Q_ASSERT( !mMesh->vertices.at( boundariesVertexIndex.last() ).isEmpty() );
1193  associateFaceToBoundaries.append( mFacesNeighborhood.at( circulator.currentFaceIndex() ).at(
1194  vertexPositionInFace( boundariesVertexIndex.last(), currentFace ) ) );
1195  }
1196  }
1197  }
1198  while ( circulator.turnCounterClockwise() != -1 && circulator.currentFaceIndex() != removedFacesIndexes.first() );
1199 
1200  bool boundaryFill = false;
1201  if ( circulator.currentFaceIndex() == -1 ) //we are on boundary of the mesh
1202  {
1203  boundaryFill = true;
1204  //we need to add last vertex/boundary faces that was not added because we are on mesh boundary
1205  circulator.goBoundaryCounterClockwise();
1206  int lastVertexIndex = circulator.oppositeVertexCounterClockwise();
1207  boundariesVertexIndex.append( lastVertexIndex );
1208 
1209  // but we can be on the case where the last vertex share an edge with the first point,
1210  // in the case, the associate face on boundarie is not -1, but the face on the other side of the edge
1211  QgsMeshVertexCirculator boundaryCirculator = vertexCirculator( lastVertexIndex );
1212  boundaryCirculator.goBoundaryCounterClockwise();
1213  if ( boundaryCirculator.oppositeVertexCounterClockwise() == boundariesVertexIndex.first() )
1214  {
1215  associateFaceToBoundaries.append( boundaryCirculator.currentFaceIndex() );
1216  boundaryFill = false; //we are not a boundary fill anymore
1217  }
1218  else
1219  associateFaceToBoundaries.append( -1 );
1220 
1221  for ( const int index : std::as_const( boundariesVertexIndex ) )
1222  {
1223  if ( isVertexOnBoundary( index ) )
1224  boundaryInGlobalMesh.insert( index );
1225  }
1226  }
1227 
1228  int currentVertexToFace = mVertexToFace.at( vertexIndex );
1229  // here, we use the method removeFaces that effectivly removes and then constructs changes
1230  Changes changes = removeFaces( removedFacesIndexes );
1231 
1232  QList<QList<int>> holes;
1233  QList<QList<int>> associateMeshFacesToHoles;
1234 
1235  bool cancelOperation = false;
1236 
1237  if ( boundaryFill )
1238  {
1239  // The hole is not a closed polygon, we need to close it, but the closing segment can intersect another segments/vertices.
1240  // In this case we consider as many polygons as necessary.
1241 
1242  int startPos = 0;
1243  int finalPos = boundariesVertexIndex.count() - 1;
1244  QList<int> uncoveredVertex;
1245 
1246  QList<int> partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
1247  QList<int> associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
1248  while ( startPos < finalPos && !partToCheck.isEmpty() )
1249  {
1250  // check if we intersect an edge between first and second
1251  int secondPos = partToCheck.count() - 1;
1252  const QgsPoint &closingSegmentExtremety1 = mMesh->vertex( partToCheck.at( 0 ) );
1253  const QgsPoint &closingSegmentExtremety2 = mMesh->vertex( partToCheck.last() );
1254  bool isEdgeIntersect = false;
1255  for ( int i = 1; i < secondPos - 1; ++i )
1256  {
1257  const QgsPoint &p1 = mMesh->vertex( partToCheck.at( i ) );
1258  const QgsPoint &p2 = mMesh->vertex( partToCheck.at( i + 1 ) );
1259  bool isLineIntersection;
1260  QgsPoint intersectPoint;
1261  isEdgeIntersect = QgsGeometryUtils::segmentIntersection( closingSegmentExtremety1, closingSegmentExtremety2, p1, p2, intersectPoint, isLineIntersection, 1e-8, true );
1262  if ( isEdgeIntersect )
1263  break;
1264  }
1265 
1266  int index = partToCheck.at( 0 );
1267  if ( boundaryInGlobalMesh.contains( index ) && index != boundariesVertexIndex.at( 0 ) )
1268  {
1269  cancelOperation = true;
1270  break;
1271  }
1272 
1273  // if uncovered vertex is a boundary vertex in the global mesh (except first that is always a boundary in the global mesh)
1274  // This operation will leads to a unique shared vertex that is not allowed, you have to cancel the operation
1275  if ( isEdgeIntersect || vertexPolygonOrientation( *mMesh, partToCheck ) >= 0 )
1276  {
1277  partToCheck.removeLast();
1278  associateFacePart.removeAt( associateFacePart.count() - 2 );
1279  if ( partToCheck.count() == 1 )
1280  {
1281  uncoveredVertex.append( index );
1282  startPos = startPos + 1;
1283  partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
1284  associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
1285  }
1286  }
1287  else
1288  {
1289  // store the well defined hole
1290  holes.append( partToCheck );
1291  associateMeshFacesToHoles.append( associateFacePart );
1292 
1293  startPos = startPos + partToCheck.count() - 1;
1294  uncoveredVertex.append( partToCheck.at( 0 ) );
1295  partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
1296  associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
1297  }
1298  }
1299  }
1300  else
1301  {
1302  holes.append( boundariesVertexIndex );
1303  associateMeshFacesToHoles.append( associateFaceToBoundaries );
1304  }
1305 
1306  if ( cancelOperation )
1307  {
1308  reverseChanges( changes );
1309  return Changes();
1310  }
1311 
1312  Q_ASSERT( holes.count() == associateMeshFacesToHoles.count() );
1313 
1314  changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
1315  changes.mVerticesToRemoveIndexes.append( vertexIndex );
1316  changes.mVerticesToFaceRemoved.append( currentVertexToFace );
1317  // these changes contain information that will lead to reference the removed vertex as free vertex when reverse/reapply
1318  dereferenceAsFreeVertex( vertexIndex );
1319  mMesh->vertices[vertexIndex] = QgsMeshVertex();
1320  mVertexToFace[vertexIndex] = -1;
1321 
1322  int oldFacesCount = mMesh->faceCount();
1323  for ( int h = 0; h < holes.count(); ++h )
1324  {
1325  const QList<int> &holeVertices = holes.at( h );
1326  const QList<int> &associateMeshFacesToHole = associateMeshFacesToHoles.at( h );
1327  QHash<p2t::Point *, int> mapPoly2TriPointToVertex;
1328  std::vector<p2t::Point *> holeToFill( holeVertices.count() );
1329  try
1330  {
1331  for ( int i = 0; i < holeVertices.count(); ++i )
1332  {
1333  const QgsMeshVertex &vertex = mMesh->vertex( holeVertices.at( i ) );
1334  holeToFill[i] = new p2t::Point( vertex.x(), vertex.y() );
1335  mapPoly2TriPointToVertex.insert( holeToFill[i], holeVertices.at( i ) );
1336  }
1337 
1338  std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( holeToFill ) );
1339 
1340  cdt->Triangulate();
1341  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
1342  QVector<QgsMeshFace> newFaces( triangles.size() );
1343  for ( size_t i = 0; i < triangles.size(); ++i )
1344  {
1345  QgsMeshFace &face = newFaces[i];
1346  face.resize( 3 );
1347  for ( int j = 0; j < 3; j++ )
1348  {
1349  int vertInd = mapPoly2TriPointToVertex.value( triangles.at( i )->GetPoint( j ), -1 );
1350  if ( vertInd == -1 )
1351  throw std::exception();
1352  Q_ASSERT( !mMesh->vertices.at( vertInd ).isEmpty() );
1353  face[j] = vertInd;
1354  }
1355  }
1356 
1357  QgsMeshEditingError error;
1358  QgsTopologicalMesh::TopologicalFaces topologicalFaces = createNewTopologicalFaces( newFaces, false, error );
1360  throw std::exception();
1361  int newFaceIndexStartIndex = mMesh->faceCount();
1362  QgsTopologicalMesh::Changes addChanges;
1363  addChanges.mAddedFacesFirstIndex = newFaceIndexStartIndex;
1364  addChanges.mFacesToAdd = topologicalFaces.meshFaces();
1365  addChanges.mFacesNeighborhoodToAdd = topologicalFaces.mFacesNeighborhood;
1366 
1367  // vertices to face changes
1368  const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
1369  for ( const int vtc : verticesToFaceToChange )
1370  if ( mVertexToFace.at( vtc ) == -1 )
1371  addChanges.mVerticesToFaceChanges.append( {vtc,
1372  mVertexToFace.at( vtc ),
1373  addChanges.addedFaceIndexInMesh( topologicalFaces.mVerticesToFace.values( vtc ).first() ) } );
1374 
1375 
1376  // reindex neighborhood for new faces
1377  for ( int i = 0; i < topologicalFaces.mFaces.count(); ++i )
1378  {
1379  FaceNeighbors &faceNeighbors = addChanges.mFacesNeighborhoodToAdd[i];
1380  faceNeighbors = topologicalFaces.mFacesNeighborhood.at( i );
1381  for ( int n = 0; n < faceNeighbors.count(); ++n )
1382  {
1383  if ( faceNeighbors.at( n ) != -1 )
1384  faceNeighbors[n] += newFaceIndexStartIndex; //reindex internal neighborhood
1385  }
1386  }
1387 
1388  // link neighborhood for boundaries of each side
1389  for ( int i = 0 ; i < holeVertices.count(); ++i )
1390  {
1391  int vertexHoleIndex = holeVertices.at( i );
1392  int meshFaceBoundaryIndex = associateMeshFacesToHole.at( i );
1393 
1394  const QgsMeshVertexCirculator circulator = QgsMeshVertexCirculator( topologicalFaces, vertexHoleIndex );
1395  circulator.goBoundaryClockwise();
1396  int newFaceBoundaryLocalIndex = circulator.currentFaceIndex();
1397  int newFaceBoundaryIndexInMesh = circulator.currentFaceIndex() + newFaceIndexStartIndex;
1398  const QgsMeshFace &newFace = circulator.currentFace();
1399  int positionInNewFaces = vertexPositionInFace( vertexHoleIndex, newFace );
1400 
1401  if ( meshFaceBoundaryIndex != -1 )
1402  {
1403  const QgsMeshFace meshFace = mMesh->face( meshFaceBoundaryIndex );
1404  int positionInMeshFaceBoundary = vertexPositionInFace( *mMesh, vertexHoleIndex, meshFaceBoundaryIndex );
1405  positionInMeshFaceBoundary = ( positionInMeshFaceBoundary - 1 + meshFace.count() ) % meshFace.count(); //take the position just before
1406 
1407  addChanges.mNeighborhoodChanges.append( {meshFaceBoundaryIndex, positionInMeshFaceBoundary, -1, newFaceBoundaryIndexInMesh} );
1408  }
1409 
1410  addChanges.mFacesNeighborhoodToAdd[newFaceBoundaryLocalIndex][positionInNewFaces] = meshFaceBoundaryIndex;
1411  }
1412 
1413  applyChanges( addChanges );
1414 
1415  changes.mFacesToAdd.append( addChanges.mFacesToAdd );
1416  changes.mFacesNeighborhoodToAdd.append( addChanges.mFacesNeighborhoodToAdd );
1417  //for each neighborhood change, check if a corresponding change already exist and merge them, if not just append
1418  for ( const std::array<int, 4> &neighborChangeToAdd : std::as_const( addChanges.mNeighborhoodChanges ) )
1419  {
1420  bool merged = false;
1421  for ( std::array<int, 4> &existingNeighborChange : changes.mNeighborhoodChanges )
1422  {
1423  if ( existingNeighborChange.at( 0 ) == neighborChangeToAdd.at( 0 ) &&
1424  existingNeighborChange.at( 1 ) == neighborChangeToAdd.at( 1 ) )
1425  {
1426  merged = true;
1427  Q_ASSERT( existingNeighborChange.at( 3 ) == neighborChangeToAdd.at( 2 ) );
1428  existingNeighborChange[3] = neighborChangeToAdd.at( 3 );
1429  }
1430  }
1431  if ( !merged )
1432  changes.mNeighborhoodChanges.append( neighborChangeToAdd );
1433  }
1434  //for each vertex to face change, check if a corresponding change already exist and merge them, if not just append
1435  for ( const std::array<int, 3> &verticesToFaceToAdd : std::as_const( addChanges.mVerticesToFaceChanges ) )
1436  {
1437  bool merged = false;
1438  for ( std::array<int, 3> &existingVerticesToFace : changes.mVerticesToFaceChanges )
1439  {
1440  if ( existingVerticesToFace.at( 0 ) == verticesToFaceToAdd.at( 0 ) )
1441  {
1442  merged = true;
1443  Q_ASSERT( existingVerticesToFace.at( 2 ) == verticesToFaceToAdd.at( 1 ) );
1444  existingVerticesToFace[2] = verticesToFaceToAdd.at( 2 );
1445  }
1446  }
1447  if ( !merged )
1448  changes.mVerticesToFaceChanges.append( verticesToFaceToAdd );
1449  }
1450 
1451  qDeleteAll( holeToFill );
1452  }
1453  catch ( ... )
1454  {
1455  qDeleteAll( holeToFill );
1456  QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping holeā€¦" ), QObject::tr( "Mesh Editing" ) );
1457  }
1458  }
1459  changes.mAddedFacesFirstIndex = oldFacesCount;
1460 
1461 
1462 
1463  changes.mAddedFacesFirstIndex = oldFacesCount;
1464 
1465  return changes;
1466 }
1467 
1469 {
1470  QSet<int> facesIndex;
1471  //Search for associated faces
1472  for ( int vertexIndex : vertices )
1473  facesIndex.unite( qgis::listToSet( facesAroundVertex( vertexIndex ) ) );
1474 
1475  // remove the faces
1476  Changes changes = removeFaces( facesIndex.values() );
1477 
1478  // removes the vertices
1479  for ( int vertexIndex : vertices )
1480  {
1481  int currentVertexToFace = mVertexToFace.at( vertexIndex );
1482  // here, we use the method removeFaces that effectivly removes and then constructs changes
1483  changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
1484  changes.mVerticesToRemoveIndexes.append( vertexIndex );
1485  changes.mVerticesToFaceRemoved.append( currentVertexToFace );
1486  // these changes contain information that will lead to reference the removed vertex as free vertex when reverse/reapply
1487  dereferenceAsFreeVertex( vertexIndex );
1488  mMesh->vertices[vertexIndex] = QgsMeshVertex();
1489  mVertexToFace[vertexIndex] = -1;
1490  }
1491 
1492  return changes;
1493 }
1494 
1496 {
1497  QList<int> boundariesToCheckClockwiseInNewFaces = topologicFaces.mBoundaries;
1498  QList<std::array<int, 2>> boundariesToCheckCounterClockwiseInNewFaces; //couple boundary / associate face in topologicFaces.mVerticesToFace
1499  QList<int> uniqueSharedVertexBoundary;
1500 
1501 
1502  // Go through the boundary and search for opposite boundary vertex clockwise in new faces and compare
1503  // with boundary opposite vertices on the mesh
1504  // If, in the mesh, the opposite vertex counter clockwise is not the same , another check will be done
1505  // later with counter clockwise in new faces
1506  // If, in the mesh, the opposite vertex clockwise is the same, this is an error
1507  while ( !boundariesToCheckClockwiseInNewFaces.isEmpty() )
1508  {
1509  int boundary = boundariesToCheckClockwiseInNewFaces.takeLast();
1510 
1511  const QList<int> &linkedFaces = topologicFaces.mVerticesToFace.values( boundary );
1512 
1513  for ( int const linkedFace : linkedFaces )
1514  {
1515 
1516  //if the boundary is a free vertex in the destination mesh, no need to check
1517  if ( mVertexToFace.at( boundary ) == -1 )
1518  continue;
1519 
1520  QgsMeshVertexCirculator newFacescirculator( topologicFaces, linkedFace, boundary );
1521  QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
1522 
1523  if ( !newFacescirculator.isValid() )
1525 
1526  //Get the opposite vertex on the clockwise side with new faces block
1527  if ( !newFacescirculator.goBoundaryClockwise() )
1529 
1530  int oppositeVertexInNewFaces = newFacescirculator.oppositeVertexClockwise();
1531 
1532  if ( !meshCirculator.goBoundaryCounterClockwise() )
1534 
1535  int oppositeVertexCCWInMesh = meshCirculator.oppositeVertexCounterClockwise();
1536 
1537  if ( oppositeVertexCCWInMesh == oppositeVertexInNewFaces ) //this boundary is OK, continue wit next one
1538  continue;
1539  else
1540  {
1541  //to avoid manifold face that could pass through the check, compare not only the boundary edges but also with the opposite internal edge in new face
1542  const QgsMeshFace &newFaceOnBoundary = newFacescirculator.currentFace();
1543  int faceSize = newFaceOnBoundary.size();
1544  int posInNewFace = vertexPositionInFace( boundary, newFaceOnBoundary );
1545  int previousVertexIndex = ( posInNewFace + faceSize - 1 ) % faceSize;
1546  if ( newFaceOnBoundary.at( previousVertexIndex ) == oppositeVertexCCWInMesh )
1548  }
1549 
1550  meshCirculator.goBoundaryClockwise();
1551 
1552  int oppositeVertexCWInMesh = meshCirculator.oppositeVertexClockwise();
1553 
1554  if ( oppositeVertexCWInMesh == oppositeVertexInNewFaces )
1556 
1557  //if we are here we need more checks
1558  boundariesToCheckCounterClockwiseInNewFaces.append( {boundary, linkedFace} );
1559  }
1560  }
1561 
1562  // Check now with opposite boundary vertex counter clockwise in new faces
1563  while ( !boundariesToCheckCounterClockwiseInNewFaces.isEmpty() )
1564  {
1565  std::array<int, 2> boundaryLinkedface = boundariesToCheckCounterClockwiseInNewFaces.takeLast();
1566  int boundary = boundaryLinkedface.at( 0 );
1567  int linkedFace = boundaryLinkedface.at( 1 );
1568 
1569  QgsMeshVertexCirculator newFacescirculator( topologicFaces, linkedFace, boundary );
1570  QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
1571 
1572  if ( !newFacescirculator.goBoundaryCounterClockwise() )
1574 
1575  int oppositeVertexInNewFaces = newFacescirculator.oppositeVertexCounterClockwise();
1576 
1577  if ( !meshCirculator.goBoundaryClockwise() )
1579 
1580  int oppositeVertexCWInMesh = meshCirculator.oppositeVertexClockwise();
1581 
1582  if ( oppositeVertexCWInMesh == oppositeVertexInNewFaces ) //this boundary is OK, continue with next one
1583  continue;
1584 
1585  meshCirculator.goBoundaryCounterClockwise();
1586 
1587  int oppositeVertexCCWInMesh = meshCirculator.oppositeVertexCounterClockwise();
1588 
1589  if ( oppositeVertexCCWInMesh == oppositeVertexInNewFaces )
1591 
1592  uniqueSharedVertexBoundary.append( boundary );
1593  }
1594 
1595  if ( !uniqueSharedVertexBoundary.isEmpty() )
1596  return QgsMeshEditingError( Qgis::MeshEditingErrorType::UniqueSharedVertex, uniqueSharedVertexBoundary.first() );
1597 
1598  // Check if internal vertices of new faces block are free in the mesh
1599  QSet<int> boundaryVertices = qgis::listToSet( topologicFaces.mBoundaries );
1600  for ( const QgsMeshFace &newFace : std::as_const( topologicFaces.mFaces ) )
1601  {
1602  for ( const int vertexIndex : newFace )
1603  {
1604  if ( boundaryVertices.contains( vertexIndex ) )
1605  continue;
1606  if ( mVertexToFace.at( vertexIndex ) != -1 )
1608  }
1609  }
1610 
1611  return QgsMeshEditingError();
1612 }
1613 
1615 {
1616  mFaces.clear();
1617  mFacesNeighborhood.clear();
1618  mVerticesToFace.clear();
1619  mBoundaries.clear();
1620 }
1621 
1622 QVector<QgsTopologicalMesh::FaceNeighbors> QgsTopologicalMesh::TopologicalFaces::facesNeighborhood() const
1623 {
1624  return mFacesNeighborhood;
1625 }
1626 
1628 {
1629  if ( mVerticesToFace.contains( vertexIndex ) )
1630  return mVerticesToFace.values( vertexIndex ).at( 0 );
1631 
1632  return -1;
1633 }
1634 
1636 {
1637  QgsTopologicalMesh topologicMesh;
1638  topologicMesh.mMesh = mesh;
1639  topologicMesh.mVertexToFace = QVector<int>( mesh->vertexCount(), -1 );
1640  topologicMesh.mMaximumVerticesPerFace = maxVerticesPerFace;
1642 
1643  for ( int i = 0; i < mesh->faceCount(); ++i )
1644  {
1645  if ( mesh->face( i ).isEmpty() )
1646  continue;
1647  if ( maxVerticesPerFace != 0 && mesh->face( i ).count() > maxVerticesPerFace )
1648  {
1650  break;
1651  }
1652 
1653  error = counterClockwiseFaces( mesh->faces[i], mesh );
1655  {
1657  error.elementIndex = i;
1658  break;
1659  }
1660  }
1661 
1663  {
1664  TopologicalFaces subMesh = createTopologicalFaces( mesh->faces, &topologicMesh.mVertexToFace, error, false );
1665  topologicMesh.mFacesNeighborhood = subMesh.mFacesNeighborhood;
1666 
1667  for ( int i = 0; i < topologicMesh.mMesh->vertexCount(); ++i )
1668  {
1669  if ( topologicMesh.mVertexToFace.at( i ) == -1 )
1670  topologicMesh.mFreeVertices.insert( i );
1671  }
1672  }
1673 
1674  return topologicMesh;
1675 }
1676 
1677 QgsTopologicalMesh::TopologicalFaces QgsTopologicalMesh::createNewTopologicalFaces( const QVector<QgsMeshFace> &faces, bool uniqueSharedVertexAllowed, QgsMeshEditingError &error )
1678 {
1679  return createTopologicalFaces( faces, nullptr, error, uniqueSharedVertexAllowed );
1680 }
1681 
1682 
1683 QgsTopologicalMesh::TopologicalFaces QgsTopologicalMesh::createTopologicalFaces(
1684  const QVector<QgsMeshFace> &faces,
1685  QVector<int> *globalVertexToFace,
1686  QgsMeshEditingError &error,
1687  bool allowUniqueSharedVertex )
1688 {
1689  int facesCount = faces.count();
1690  QVector<FaceNeighbors> faceTopologies;
1691  QMultiHash<int, int> verticesToFace;
1692 
1693  error = QgsMeshEditingError();
1694  TopologicalFaces ret;
1695 
1696  // Contains for each vertex a map (opposite vertex # edge) --> face index
1697  // when turning counter clockwise, if v1 first vertex and v2 second one, [v1][v2]--> neighbor face
1698  QMap<int, QMap<int, int>> verticesToNeighbor;
1699 
1700  for ( int faceIndex = 0; faceIndex < facesCount; ++faceIndex )
1701  {
1702  const QgsMeshFace &face = faces.at( faceIndex );
1703  int faceSize = face.count();
1704  //Fill vertices to neighbor faces
1705  for ( int i = 0; i < faceSize; ++i )
1706  {
1707  int v1 = face[i % faceSize];
1708  int v2 = face[( i + 1 ) % faceSize];
1709  if ( verticesToNeighbor[v2].contains( v1 ) )
1710  {
1712  return ret;
1713  }
1714  else
1715  verticesToNeighbor[v2].insert( v1, faceIndex );
1716  }
1717  }
1718 
1719  faceTopologies = QVector<FaceNeighbors>( faces.count() );
1720 
1721  QSet<int> boundaryVertices;
1722 
1723  for ( int faceIndex = 0; faceIndex < facesCount; ++faceIndex )
1724  {
1725  const QgsMeshFace &face = faces.at( faceIndex );
1726  int faceSize = face.size();
1727  FaceNeighbors &faceTopology = faceTopologies[faceIndex];
1728  faceTopology.resize( faceSize );
1729 
1730  for ( int i = 0; i < faceSize; ++i )
1731  {
1732  int v1 = face.at( i );
1733  int v2 = face.at( ( i + 1 ) % faceSize );
1734 
1735  if ( globalVertexToFace )
1736  {
1737  if ( ( *globalVertexToFace )[v1] == -1 )
1738  ( *globalVertexToFace )[v1] = faceIndex ;
1739  }
1740  else
1741  {
1742  if ( allowUniqueSharedVertex || !verticesToFace.contains( v1 ) )
1743  verticesToFace.insert( v1, faceIndex ) ;
1744  }
1745 
1746  QMap<int, int> &edges = verticesToNeighbor[v1];
1747  if ( edges.contains( v2 ) )
1748  faceTopology[i] = edges.value( v2 );
1749  else
1750  {
1751  faceTopology[i] = -1;
1752 
1753  if ( !allowUniqueSharedVertex )
1754  {
1755  if ( boundaryVertices.contains( v1 ) )
1756  {
1757  error = QgsMeshEditingError( Qgis::MeshEditingErrorType::UniqueSharedVertex, v1 ); // if a vertices is more than one time in the boundary, that means faces share only one vertices
1758  return ret;
1759  }
1760  }
1761  boundaryVertices.insert( v1 );
1762  }
1763  }
1764  }
1765 
1766  ret.mFaces = faces;
1767  ret.mFacesNeighborhood = faceTopologies;
1768  ret.mBoundaries = boundaryVertices.values();
1769  ret.mVerticesToFace = verticesToFace;
1770  return ret;
1771 }
1772 
1773 QVector<int> QgsTopologicalMesh::neighborsOfFace( int faceIndex ) const
1774 {
1775  return mFacesNeighborhood.at( faceIndex );
1776 }
1777 
1778 QList<int> QgsTopologicalMesh::facesAroundVertex( int vertexIndex ) const
1779 {
1780  QgsMeshVertexCirculator circ = vertexCirculator( vertexIndex );
1781 
1782  return circ.facesAround();
1783 }
1784 
1786 {
1787  QSet<int> removedFaces = qgis::listToSet( facesIndexes );
1788  QSet<int> concernedFaces = concernedFacesBy( facesIndexes );
1789 
1790  for ( const int f : std::as_const( removedFaces ) )
1791  concernedFaces.remove( f );
1792 
1793  QVector<QgsMeshFace> remainingFaces;
1794  remainingFaces.reserve( concernedFaces.count() );
1795  for ( const int f : std::as_const( concernedFaces ) )
1796  remainingFaces.append( mMesh->face( f ) );
1797 
1798  QgsMeshEditingError error;
1799  createTopologicalFaces( remainingFaces, nullptr, error, false );
1800 
1801  return error;
1802 }
1803 
1804 QgsTopologicalMesh::Changes QgsTopologicalMesh::removeFaces( const QList<int> facesIndexesToRemove )
1805 {
1806  Changes changes;
1807  changes.mFaceIndexesToRemove = facesIndexesToRemove;
1808  changes.mFacesToRemove.resize( facesIndexesToRemove.count() );
1809  changes.mFacesNeighborhoodToRemove.resize( facesIndexesToRemove.count() );
1810 
1811  QSet<int> indexSet = qgis::listToSet( facesIndexesToRemove );
1812  QSet<int> threatedVertex;
1813 
1814  for ( int i = 0; i < facesIndexesToRemove.count(); ++i )
1815  {
1816  const int faceIndex = facesIndexesToRemove.at( i );
1817  const QgsMeshFace &face = mMesh->face( faceIndex );
1818  changes.mFacesToRemove[i] = face;
1819  const FaceNeighbors &neighborhood = mFacesNeighborhood.at( faceIndex );
1820  changes.mFacesNeighborhoodToRemove[i] = neighborhood;
1821  for ( int j = 0; j < face.count(); ++j )
1822  {
1823  //change the neighborhood for each neighbor face
1824  int neighborIndex = neighborhood.at( j );
1825  if ( neighborIndex != -1 && !indexSet.contains( neighborIndex ) )
1826  {
1827  int positionInNeighbor = mFacesNeighborhood.at( neighborIndex ).indexOf( faceIndex );
1828  changes.mNeighborhoodChanges.append( {neighborIndex, positionInNeighbor, faceIndex, -1} );
1829  }
1830 
1831  //change vertexToFace
1832  int vertexIndex = face.at( j );
1833  if ( !threatedVertex.contains( vertexIndex ) && indexSet.contains( mVertexToFace.at( vertexIndex ) ) )
1834  {
1835  int oldValue = mVertexToFace.at( vertexIndex );
1836  //look for another face linked to this vertex
1837  int refValue = -1;
1838  if ( neighborIndex != -1 && !indexSet.contains( neighborIndex ) ) //if exist, simpler to take it
1839  refValue = neighborIndex;
1840  else
1841  {
1842  QList<int> aroundFaces = facesAroundVertex( vertexIndex );
1843  aroundFaces.removeOne( faceIndex );
1844  if ( !aroundFaces.isEmpty() )
1845  {
1846  while ( !aroundFaces.isEmpty() && refValue == -1 )
1847  {
1848  if ( !indexSet.contains( aroundFaces.first() ) )
1849  refValue = aroundFaces.first();
1850  else
1851  aroundFaces.removeFirst();
1852  }
1853  }
1854  }
1855  changes.mVerticesToFaceChanges.append( {vertexIndex, oldValue, refValue} );
1856  threatedVertex.insert( vertexIndex );
1857  }
1858  }
1859  }
1860 
1861  applyChanges( changes );
1862 
1863  return changes;
1864 }
1865 
1866 bool QgsTopologicalMesh::eitherSideFacesAndVertices( int vertexIndex1,
1867  int vertexIndex2,
1868  int &face1,
1869  int &face2,
1870  int &neighborVertex1InFace1,
1871  int &neighborVertex1InFace2,
1872  int &neighborVertex2inFace1,
1873  int &neighborVertex2inFace2 ) const
1874 {
1875  QgsMeshVertexCirculator circulator1 = vertexCirculator( vertexIndex1 );
1876  QgsMeshVertexCirculator circulator2 = vertexCirculator( vertexIndex2 );
1877 
1878  circulator1.goBoundaryClockwise();
1879  int firstFace1 = circulator1.currentFaceIndex();
1880  circulator2.goBoundaryClockwise();
1881  int firstFace2 = circulator2.currentFaceIndex();
1882 
1883  if ( circulator1.oppositeVertexCounterClockwise() != vertexIndex2 )
1884  while ( circulator1.turnCounterClockwise() != -1 &&
1885  circulator1.currentFaceIndex() != firstFace1 &&
1886  circulator1.oppositeVertexCounterClockwise() != vertexIndex2 ) {}
1887 
1888  if ( circulator2.oppositeVertexCounterClockwise() != vertexIndex1 )
1889  while ( circulator2.turnCounterClockwise() != -1 &&
1890  circulator2.currentFaceIndex() != firstFace2 &&
1891  circulator2.oppositeVertexCounterClockwise() != vertexIndex1 ) {}
1892 
1893  if ( circulator1.oppositeVertexCounterClockwise() != vertexIndex2
1894  || circulator2.oppositeVertexCounterClockwise() != vertexIndex1 )
1895  return false;
1896 
1897  face1 = circulator1.currentFaceIndex();
1898  face2 = circulator2.currentFaceIndex();
1899 
1900  neighborVertex1InFace1 = circulator1.oppositeVertexClockwise();
1901  circulator1.turnCounterClockwise();
1902  neighborVertex1InFace2 = circulator1.oppositeVertexCounterClockwise();
1903 
1904  neighborVertex2inFace2 = circulator2.oppositeVertexClockwise();
1905  circulator2.turnCounterClockwise();
1906  neighborVertex2inFace1 = circulator2.oppositeVertexCounterClockwise();
1907 
1908  return true;
1909 }
1910 
1911 bool QgsTopologicalMesh::edgeCanBeFlipped( int vertexIndex1, int vertexIndex2 ) const
1912 {
1913  int faceIndex1;
1914  int faceIndex2;
1915  int oppositeVertexFace1;
1916  int oppositeVertexFace2;
1917  int supposedOppositeVertexFace1;
1918  int supposedoppositeVertexFace2;
1919 
1920  bool result = eitherSideFacesAndVertices(
1921  vertexIndex1,
1922  vertexIndex2,
1923  faceIndex1,
1924  faceIndex2,
1925  oppositeVertexFace1,
1926  supposedoppositeVertexFace2,
1927  supposedOppositeVertexFace1,
1928  oppositeVertexFace2 );
1929 
1930  if ( !result ||
1931  faceIndex1 < 0 ||
1932  faceIndex2 < 0 ||
1933  oppositeVertexFace1 < 0 ||
1934  oppositeVertexFace2 < 0 ||
1935  supposedOppositeVertexFace1 != oppositeVertexFace1 ||
1936  supposedoppositeVertexFace2 != oppositeVertexFace2 )
1937  return false;
1938 
1939  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
1940  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
1941 
1942 
1943  if ( face1.count() != 3 || face2.count() != 3 )
1944  return false;
1945 
1946  double crossProduct1 = crossProduct( vertexIndex1, oppositeVertexFace1, oppositeVertexFace2, *mMesh );
1947  double crossProduct2 = crossProduct( vertexIndex2, oppositeVertexFace1, oppositeVertexFace2, *mMesh );
1948 
1949  return crossProduct1 * crossProduct2 < 0;
1950 }
1951 
1952 QgsTopologicalMesh::Changes QgsTopologicalMesh::flipEdge( int vertexIndex1, int vertexIndex2 )
1953 {
1954  int faceIndex1;
1955  int faceIndex2;
1956  int oppositeVertexFace1;
1957  int oppositeVertexFace2;
1958  int supposedOppositeVertexFace1;
1959  int supposedoppositeVertexFace2;
1960 
1961  bool result = eitherSideFacesAndVertices(
1962  vertexIndex1,
1963  vertexIndex2,
1964  faceIndex1,
1965  faceIndex2,
1966  oppositeVertexFace1,
1967  supposedoppositeVertexFace2,
1968  supposedOppositeVertexFace1,
1969  oppositeVertexFace2 );
1970 
1971  if ( !result ||
1972  faceIndex1 < 0 ||
1973  faceIndex2 < 0 ||
1974  oppositeVertexFace1 < 0 ||
1975  oppositeVertexFace2 < 0 ||
1976  supposedOppositeVertexFace1 != oppositeVertexFace1 ||
1977  supposedoppositeVertexFace2 != oppositeVertexFace2 )
1978  return Changes();
1979 
1980 
1981  Changes changes;
1982 
1983  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
1984  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
1985 
1986  Q_ASSERT( face1.count() == 3 );
1987  Q_ASSERT( face2.count() == 3 );
1988 
1989  int pos1 = vertexPositionInFace( vertexIndex1, face1 );
1990  int pos2 = vertexPositionInFace( vertexIndex2, face2 );
1991 
1992  int neighborFace1 = mFacesNeighborhood.at( faceIndex1 ).at( pos1 );
1993  int posInNeighbor1 = vertexPositionInFace( *mMesh, oppositeVertexFace1, neighborFace1 );
1994  int neighborFace2 = mFacesNeighborhood.at( faceIndex1 ).at( ( pos1 + 1 ) % 3 );
1995  int posInNeighbor2 = vertexPositionInFace( *mMesh, vertexIndex2, neighborFace2 );
1996  int neighborFace3 = mFacesNeighborhood.at( faceIndex2 ).at( pos2 );
1997  int posInNeighbor3 = vertexPositionInFace( *mMesh, oppositeVertexFace2, neighborFace3 );
1998  int neighborFace4 = mFacesNeighborhood.at( faceIndex2 ).at( ( pos2 + 1 ) % 3 );
1999  int posInNeighbor4 = vertexPositionInFace( *mMesh, vertexIndex1, neighborFace4 );
2000 
2001  changes.mFaceIndexesToRemove.append( faceIndex1 );
2002  changes.mFaceIndexesToRemove.append( faceIndex2 );
2003  changes.mFacesToRemove.append( face1 );
2004  changes.mFacesToRemove.append( face2 );
2005  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex1 ) );
2006  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex2 ) );
2007  int startIndex = mMesh->faceCount();
2008  changes.mAddedFacesFirstIndex = startIndex;
2009  changes.mFacesToAdd.append( {oppositeVertexFace1, oppositeVertexFace2, vertexIndex1} );
2010  changes.mFacesToAdd.append( {oppositeVertexFace2, oppositeVertexFace1, vertexIndex2} );
2011  changes.mFacesNeighborhoodToAdd.append( {startIndex + 1,
2012  mFacesNeighborhood.at( faceIndex2 ).at( ( pos2 + 1 ) % 3 ),
2013  mFacesNeighborhood.at( faceIndex1 ).at( pos1 )} );
2014  changes.mFacesNeighborhoodToAdd.append( {startIndex,
2015  mFacesNeighborhood.at( faceIndex1 ).at( ( pos1 + 1 ) % 3 ),
2016  mFacesNeighborhood.at( faceIndex2 ).at( pos2 )} );
2017 
2018  if ( neighborFace1 >= 0 )
2019  changes.mNeighborhoodChanges.append( {neighborFace1, posInNeighbor1, faceIndex1, startIndex} );
2020  if ( neighborFace2 >= 0 )
2021  changes.mNeighborhoodChanges.append( {neighborFace2, posInNeighbor2, faceIndex1, startIndex + 1} );
2022  if ( neighborFace3 >= 0 )
2023  changes.mNeighborhoodChanges.append( {neighborFace3, posInNeighbor3, faceIndex2, startIndex + 1} );
2024  if ( neighborFace4 >= 0 )
2025  changes.mNeighborhoodChanges.append( {neighborFace4, posInNeighbor4, faceIndex2, startIndex} );
2026 
2027 
2028  if ( mVertexToFace.at( vertexIndex1 ) == faceIndex1 || mVertexToFace.at( vertexIndex1 ) == faceIndex2 )
2029  changes.mVerticesToFaceChanges.append( {vertexIndex1, mVertexToFace.at( vertexIndex1 ), startIndex} );
2030  if ( mVertexToFace.at( vertexIndex2 ) == faceIndex1 || mVertexToFace.at( vertexIndex2 ) == faceIndex2 )
2031  changes.mVerticesToFaceChanges.append( {vertexIndex2, mVertexToFace.at( vertexIndex2 ), startIndex + 1} );
2032 
2033  if ( mVertexToFace.at( oppositeVertexFace1 ) == faceIndex1 )
2034  changes.mVerticesToFaceChanges.append( {oppositeVertexFace1, faceIndex1, startIndex} );
2035 
2036  if ( mVertexToFace.at( oppositeVertexFace2 ) == faceIndex2 )
2037  changes.mVerticesToFaceChanges.append( {oppositeVertexFace2, faceIndex2, startIndex + 1} );
2038 
2039  applyChanges( changes );
2040 
2041  return changes;
2042 }
2043 
2044 bool QgsTopologicalMesh::canBeMerged( int vertexIndex1, int vertexIndex2 ) const
2045 {
2046  int faceIndex1;
2047  int faceIndex2;
2048  int neighborVertex1InFace1;
2049  int neighborVertex1InFace2;
2050  int neighborVertex2inFace1;
2051  int neighborVertex2inFace2;
2052 
2053  bool result = eitherSideFacesAndVertices(
2054  vertexIndex1,
2055  vertexIndex2,
2056  faceIndex1,
2057  faceIndex2,
2058  neighborVertex1InFace1,
2059  neighborVertex1InFace2,
2060  neighborVertex2inFace1,
2061  neighborVertex2inFace2 );
2062 
2063  if ( !result ||
2064  faceIndex1 < 0 ||
2065  faceIndex2 < 0 )
2066  return false;
2067 
2068  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
2069  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
2070 
2071  if ( face1.count() + face2.count() - 2 > mMaximumVerticesPerFace )
2072  return false;
2073 
2074  QgsMeshVertex v1 = mMesh->vertices.at( vertexIndex1 );
2075  QgsMeshVertex v2 = mMesh->vertices.at( vertexIndex2 );
2076  QgsMeshVertex nv11 = mMesh->vertices.at( neighborVertex1InFace1 );
2077  QgsMeshVertex nv12 = mMesh->vertices.at( neighborVertex1InFace2 );
2078  QgsMeshVertex nv21 = mMesh->vertices.at( neighborVertex2inFace1 );
2079  QgsMeshVertex nv22 = mMesh->vertices.at( neighborVertex2inFace2 );
2080 
2081  double crossProduct1 = crossProduct( vertexIndex1, neighborVertex1InFace1, neighborVertex1InFace2, *mMesh );
2082  double crossProduct2 = crossProduct( vertexIndex2, neighborVertex2inFace1, neighborVertex2inFace2, *mMesh );
2083 
2084  return crossProduct1 * crossProduct2 < 0;
2085 }
2086 
2087 QgsTopologicalMesh::Changes QgsTopologicalMesh::merge( int vertexIndex1, int vertexIndex2 )
2088 {
2089  int faceIndex1;
2090  int faceIndex2;
2091  int neighborVertex1InFace1;
2092  int neighborVertex1InFace2;
2093  int neighborVertex2inFace1;
2094  int neighborVertex2inFace2;
2095 
2096  bool result = eitherSideFacesAndVertices(
2097  vertexIndex1,
2098  vertexIndex2,
2099  faceIndex1,
2100  faceIndex2,
2101  neighborVertex1InFace1,
2102  neighborVertex1InFace2,
2103  neighborVertex2inFace1,
2104  neighborVertex2inFace2 );
2105 
2106  if ( !result ||
2107  faceIndex1 < 0 ||
2108  faceIndex2 < 0 )
2109  return Changes();
2110 
2111  Changes changes;
2112 
2113  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
2114  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
2115  int faceSize1 = face1.count();
2116  int faceSize2 = face2.count();
2117 
2118  int pos1 = vertexPositionInFace( vertexIndex1, face1 );
2119  int pos2 = vertexPositionInFace( vertexIndex2, face2 );
2120 
2121  changes.mFaceIndexesToRemove.append( faceIndex1 );
2122  changes.mFaceIndexesToRemove.append( faceIndex2 );
2123  changes.mFacesToRemove.append( face1 );
2124  changes.mFacesToRemove.append( face2 );
2125  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex1 ) );
2126  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex2 ) );
2127  int startIndex = mMesh->faceCount();
2128  changes.mAddedFacesFirstIndex = startIndex;
2129 
2130  QgsMeshFace newface;
2131  FaceNeighbors newNeighborhood;
2132  for ( int i = 0; i < faceSize1 - 1; ++i )
2133  {
2134  int currentPos = ( pos1 + i ) % faceSize1;
2135  newface.append( face1.at( currentPos ) ); //add vertex of face1
2136 
2137  int currentNeighbor = mFacesNeighborhood.at( faceIndex1 ).at( currentPos );
2138  newNeighborhood.append( currentNeighbor );
2139 
2140  if ( currentNeighbor != -1 )
2141  {
2142  int currentPosInNeighbor = vertexPositionInFace( *mMesh, face1.at( ( currentPos + 1 ) % faceSize1 ), currentNeighbor );
2143  changes.mNeighborhoodChanges.append( {currentNeighbor, currentPosInNeighbor, faceIndex1, startIndex} );
2144  }
2145  }
2146  for ( int i = 0; i < faceSize2 - 1; ++i )
2147  {
2148  int currentPos = ( pos2 + i ) % faceSize2;
2149  newface.append( face2.at( currentPos ) ); //add vertex of face2
2150 
2151  int currentNeighbor = mFacesNeighborhood.at( faceIndex2 ).at( currentPos );
2152  newNeighborhood.append( currentNeighbor );
2153 
2154  if ( currentNeighbor != -1 )
2155  {
2156  int currentPosInNeighbor = vertexPositionInFace( *mMesh, face2.at( ( currentPos + 1 ) % faceSize2 ), currentNeighbor );
2157  changes.mNeighborhoodChanges.append( {currentNeighbor, currentPosInNeighbor, faceIndex2, startIndex} );
2158  }
2159  }
2160 
2161  for ( int i = 0; i < faceSize1; ++i )
2162  if ( mVertexToFace.at( face1.at( i ) ) == faceIndex1 )
2163  changes.mVerticesToFaceChanges.append( {face1.at( i ), faceIndex1, startIndex} );
2164 
2165  for ( int i = 0; i < faceSize2; ++i )
2166  if ( mVertexToFace.at( face2.at( i ) ) == faceIndex2 )
2167  changes.mVerticesToFaceChanges.append( {face2.at( i ), faceIndex2, startIndex} );
2168 
2169  changes.mFacesToAdd.append( newface );
2170  changes.mFacesNeighborhoodToAdd.append( newNeighborhood );
2171 
2172  applyChanges( changes );
2173 
2174  return changes;
2175 }
2176 
2177 bool QgsTopologicalMesh::canBeSplit( int faceIndex ) const
2178 {
2179  const QgsMeshFace face = mMesh->face( faceIndex );
2180 
2181  return face.count() == 4;
2182 }
2183 
2185 {
2186  //search for the spliited angle (greater angle)
2187  const QgsMeshFace &face = mMesh->face( faceIndex );
2188  int faceSize = face.count();
2189 
2190  Q_ASSERT( faceSize == 4 );
2191 
2192  double maxAngle = 0;
2193  int splitVertexPos = -1;
2194  for ( int i = 0; i < faceSize; ++i )
2195  {
2196  QgsVector vect1( mMesh->vertex( face.at( i ) ) - mMesh->vertex( face.at( ( i + 1 ) % faceSize ) ) );
2197  QgsVector vect2( mMesh->vertex( face.at( ( i + 2 ) % faceSize ) ) - mMesh->vertex( face.at( ( i + 1 ) % faceSize ) ) );
2198 
2199  double angle = std::abs( vect1.angle( vect2 ) );
2200  angle = std::min( angle, 2.0 * M_PI - angle );
2201  if ( angle > maxAngle )
2202  {
2203  maxAngle = angle;
2204  splitVertexPos = ( i + 1 ) % faceSize;
2205  }
2206  }
2207 
2208  Changes changes;
2209 
2210  const QgsMeshFace newFace1 = {face.at( splitVertexPos ),
2211  face.at( ( splitVertexPos + 1 ) % faceSize ),
2212  face.at( ( splitVertexPos + 2 ) % faceSize )
2213  };
2214 
2215  const QgsMeshFace newFace2 = {face.at( splitVertexPos ),
2216  face.at( ( splitVertexPos + 2 ) % faceSize ),
2217  face.at( ( splitVertexPos + 3 ) % faceSize )
2218  };
2219 
2220  QVector<int> neighborIndex( faceSize );
2221  QVector<int> posInNeighbor( faceSize );
2222 
2223  for ( int i = 0; i < faceSize; ++i )
2224  {
2225  neighborIndex[i] = mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + i ) % faceSize );
2226  posInNeighbor[i] = vertexPositionInFace( *mMesh, face.at( ( splitVertexPos + i + 1 ) % faceSize ), neighborIndex[i] );
2227  }
2228 
2229  changes.mFaceIndexesToRemove.append( faceIndex );
2230  changes.mFacesToRemove.append( face );
2231  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex ) );
2232  int startIndex = mMesh->faceCount();
2233  changes.mAddedFacesFirstIndex = startIndex;
2234  changes.mFacesToAdd.append( newFace1 );
2235  changes.mFacesToAdd.append( newFace2 );
2236 
2237  changes.mFacesNeighborhoodToAdd.append( {mFacesNeighborhood.at( faceIndex ).at( splitVertexPos ),
2238  mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 1 ) % faceSize ),
2239  startIndex + 1} );
2240  changes.mFacesNeighborhoodToAdd.append( {startIndex,
2241  mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 2 ) % faceSize ),
2242  mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 3 ) % faceSize )} );
2243 
2244  for ( int i = 0; i < faceSize; ++i )
2245  {
2246  if ( neighborIndex[i] >= 0 )
2247  changes.mNeighborhoodChanges.append( {neighborIndex[i], posInNeighbor[i], faceIndex, startIndex + int( i / 2 )} );
2248 
2249  int vertexIndex = face.at( ( splitVertexPos + i ) % faceSize );
2250  if ( mVertexToFace.at( vertexIndex ) == faceIndex )
2251  changes.mVerticesToFaceChanges.append( {vertexIndex, faceIndex, startIndex + int( i / 2 )} );
2252  }
2253 
2254  applyChanges( changes );
2255 
2256  return changes;
2257 }
2258 
2259 
2261 {
2262  Changes changes;
2263  changes.mVerticesToAdd.append( vertex );
2264  changes.mVertexToFaceToAdd.append( -1 );
2265 
2266  mMesh->vertices.append( vertex );
2267  mVertexToFace.append( -1 );
2268  changes.mAddedFacesFirstIndex = mMesh->faceCount();
2269 
2270  const QgsMeshFace includingFace = mMesh->face( includingFaceIndex );
2271  const FaceNeighbors includingFaceNeighborhood = mFacesNeighborhood.at( includingFaceIndex );
2272  int includingFaceSize = includingFace.count();
2273 
2274  for ( int i = 0; i < includingFaceSize; ++i )
2275  {
2276  // add a new face
2277  QgsMeshFace face( 3 );
2278  face[0] = mMesh->vertexCount() - 1;
2279  face[1] = includingFace.at( i );
2280  face[2] = includingFace.at( ( i + 1 ) % includingFaceSize );
2281  mMesh->faces.append( face );
2282  changes.mFacesToAdd.append( face );
2283 
2284  int currentVertexIndex = includingFace.at( i );
2285  if ( mVertexToFace.at( currentVertexIndex ) == includingFaceIndex )
2286  {
2287  int newFaceIndex = mMesh->faceCount() - 1;
2288  mVertexToFace[currentVertexIndex] = newFaceIndex;
2289  changes.mVerticesToFaceChanges.append( {currentVertexIndex, includingFaceIndex, newFaceIndex} );
2290  }
2291 
2292  int includingFaceNeighbor = includingFaceNeighborhood.at( i );
2293  FaceNeighbors neighbors(
2294  {
2295  changes.mAddedFacesFirstIndex + ( i + includingFaceSize - 1 ) % includingFaceSize,
2296  includingFaceNeighbor,
2297  changes.mAddedFacesFirstIndex + ( i + includingFaceSize + 1 ) % includingFaceSize
2298  } );
2299  mFacesNeighborhood.append( neighbors );
2300  changes.mFacesNeighborhoodToAdd.append( neighbors );
2301 
2302  if ( includingFaceNeighbor != -1 )
2303  {
2304  int indexInNeighbor = vertexPositionInFace( *mMesh, includingFace.at( ( i + 1 ) % includingFaceSize ), includingFaceNeighbor );
2305  int oldValue = mFacesNeighborhood[includingFaceNeighbor][indexInNeighbor];
2306  mFacesNeighborhood[includingFaceNeighbor][indexInNeighbor] = changes.mAddedFacesFirstIndex + i;
2307  changes.mNeighborhoodChanges.append( {includingFaceNeighbor, indexInNeighbor, oldValue, changes.mAddedFacesFirstIndex + i} );
2308  }
2309  }
2310 
2311  changes.mFacesToRemove.append( includingFace );
2312  changes.mFaceIndexesToRemove.append( includingFaceIndex );
2313  changes.mFacesNeighborhoodToRemove.append( includingFaceNeighborhood );
2314 
2315  mFacesNeighborhood[includingFaceIndex] = FaceNeighbors();
2316  mMesh->faces[includingFaceIndex] = QgsMeshFace();
2317  mVertexToFace[mVertexToFace.count() - 1] = mMesh->faceCount() - 1;
2318  changes.mVertexToFaceToAdd[changes.mVertexToFaceToAdd.count() - 1] = mMesh->faceCount() - 1 ;
2319 
2320  return changes;
2321 }
2322 
2324 {
2325  const QgsMeshFace face1 = mMesh->face( faceIndex );
2326 
2327  Changes changes;
2328  changes.mVerticesToAdd.append( vertexToInsert );
2329  changes.mAddedFacesFirstIndex = mMesh->faceCount();
2330 
2331  // triangulate the first face
2332  int newVertexPositionInFace1 = position + 1;
2333 
2334  auto triangulate = [this, &changes]( int removedFaceIndex, const QgsMeshVertex & newVertex, int newVertexPosition, QVector<int> &edgeFacesIndexes )->bool
2335  {
2336  const QgsMeshFace &initialFace = mMesh->face( removedFaceIndex );
2337  changes.mFacesToRemove.append( initialFace );
2338  changes.mFaceIndexesToRemove.append( removedFaceIndex );
2339  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( removedFaceIndex ) );
2340  const int addedVertexIndex = mMesh->vertexCount();
2341 
2342  int faceStartGlobalIndex = mMesh->faceCount() + changes.mFacesToAdd.count();
2343  int localStartIndex = changes.mFacesToAdd.count();
2344 
2345  QVector<int> newBoundary = initialFace;
2346  newBoundary.insert( newVertexPosition, addedVertexIndex );
2347 
2348  try
2349  {
2350  QHash<p2t::Point *, int> mapPoly2TriPointToVertex;
2351  std::vector<p2t::Point *> faceToFill( newBoundary.count() );
2352  for ( int i = 0; i < newBoundary.count(); ++i )
2353  {
2354  QgsMeshVertex vert;
2355 
2356  if ( newBoundary.at( i ) == addedVertexIndex )
2357  vert = newVertex;
2358  else
2359  vert = mMesh->vertex( newBoundary.at( i ) );
2360 
2361  faceToFill[i] = new p2t::Point( vert.x(), vert.y() );
2362  mapPoly2TriPointToVertex.insert( faceToFill[i], newBoundary.at( i ) );
2363  }
2364 
2365  std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( faceToFill ) );
2366  cdt->Triangulate();
2367  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
2368  QVector<QgsMeshFace> newFaces( triangles.size() );
2369  for ( size_t i = 0; i < triangles.size(); ++i )
2370  {
2371  QgsMeshFace &face = newFaces[i];
2372  face.resize( 3 );
2373  for ( int j = 0; j < 3; j++ )
2374  {
2375  int vertInd = mapPoly2TriPointToVertex.value( triangles.at( i )->GetPoint( j ), -1 );
2376  if ( vertInd == -1 )
2377  throw std::exception();
2378  face[j] = vertInd;
2379  }
2380  }
2381 
2382  QgsMeshEditingError error;
2383  QgsTopologicalMesh::TopologicalFaces topologicalFaces = createNewTopologicalFaces( newFaces, false, error );
2385  throw std::exception();
2386 
2387  changes.mFacesToAdd.append( topologicalFaces.meshFaces() );
2388  changes.mFacesNeighborhoodToAdd.append( topologicalFaces.mFacesNeighborhood );
2389 
2390  // vertices to face changes
2391  const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
2392  for ( const int vtc : verticesToFaceToChange )
2393  if ( vtc != addedVertexIndex && mVertexToFace.at( vtc ) == removedFaceIndex )
2394  changes.mVerticesToFaceChanges.append(
2395  {
2396  vtc,
2397  removedFaceIndex,
2398  topologicalFaces.vertexToFace( vtc ) + faceStartGlobalIndex
2399  } );
2400 
2401  // reindex neighborhood for new faces
2402  for ( int i = 0; i < topologicalFaces.mFaces.count(); ++i )
2403  {
2404  FaceNeighbors &faceNeighbors = changes.mFacesNeighborhoodToAdd[i + localStartIndex];
2405  for ( int n = 0; n < faceNeighbors.count(); ++n )
2406  {
2407  if ( faceNeighbors.at( n ) != -1 )
2408  faceNeighbors[n] += faceStartGlobalIndex; //reindex internal neighborhood
2409  }
2410  }
2411 
2412  edgeFacesIndexes.resize( 2 );
2413  // link neighborhood for boundaries of each side
2414  for ( int i = 0 ; i < newBoundary.count(); ++i )
2415  {
2416  int vertexIndex = newBoundary.at( i );
2417  QgsMeshVertexCirculator circulator = QgsMeshVertexCirculator( topologicalFaces, vertexIndex );
2418  circulator.goBoundaryClockwise();
2419  int newFaceBoundaryLocalIndex = localStartIndex + circulator.currentFaceIndex();
2420 
2421  int newFaceBoundaryIndexInMesh = faceStartGlobalIndex;
2422 
2423  int meshFaceBoundaryIndex;
2424  if ( i == newVertexPosition )
2425  {
2426  meshFaceBoundaryIndex = -1; //face that are on the opposite side of the edge, filled later
2427  edgeFacesIndexes[0] = newFaceBoundaryLocalIndex;
2428  }
2429  else if ( i == ( newVertexPosition + newBoundary.count() - 1 ) % newBoundary.count() )
2430  {
2431  meshFaceBoundaryIndex = -1; //face that are on the opposite side of the edge, filled later
2432  edgeFacesIndexes[1] = newFaceBoundaryLocalIndex;
2433  }
2434  else
2435  meshFaceBoundaryIndex = mFacesNeighborhood.at( removedFaceIndex ).at( vertexPositionInFace( vertexIndex, initialFace ) );
2436 
2437  const QgsMeshFace &newFace = circulator.currentFace();
2438  int positionInNewFaces = vertexPositionInFace( vertexIndex, newFace );
2439 
2440  if ( meshFaceBoundaryIndex != -1 )
2441  {
2442  const QgsMeshFace meshFace = mMesh->face( meshFaceBoundaryIndex );
2443  int positionInMeshFaceBoundary = vertexPositionInFace( *mMesh, vertexIndex, meshFaceBoundaryIndex );
2444  positionInMeshFaceBoundary = ( positionInMeshFaceBoundary - 1 + meshFace.count() ) % meshFace.count(); //take the position just before
2445 
2446  changes.mNeighborhoodChanges.append( {meshFaceBoundaryIndex,
2447  positionInMeshFaceBoundary,
2448  removedFaceIndex,
2449  newFaceBoundaryIndexInMesh +
2450  circulator.currentFaceIndex()} );
2451  }
2452 
2453  changes.mFacesNeighborhoodToAdd[newFaceBoundaryLocalIndex][positionInNewFaces] = meshFaceBoundaryIndex;
2454  }
2455 
2456  qDeleteAll( faceToFill );
2457  }
2458  catch ( ... )
2459  {
2460  return false;
2461  }
2462 
2463  return true;
2464  };
2465 
2466  QVector<int> edgeFacesIndexes;
2467  if ( !triangulate( faceIndex, vertexToInsert, newVertexPositionInFace1, edgeFacesIndexes ) )
2468  return Changes();
2469 
2470  changes.mVertexToFaceToAdd.append( edgeFacesIndexes.at( 0 ) + changes.mAddedFacesFirstIndex );
2471 
2472  int addedVertexIndex = mMesh->vertexCount();
2473 
2474  //if exist, triangulate the second face if exists
2475  int face2Index = mFacesNeighborhood.at( faceIndex ).at( position );
2476  if ( face2Index != -1 )
2477  {
2478  const QgsMeshFace &face2 = mMesh->face( face2Index );
2479  int vertexPositionInFace2 = vertexPositionInFace( face1.at( position ), face2 );
2480  QVector<int> edgeFacesIndexesFace2;
2481  if ( !triangulate( face2Index, vertexToInsert, vertexPositionInFace2, edgeFacesIndexesFace2 ) )
2482  return Changes();
2483 
2484  //link neighborhood with other side
2485  const QgsMeshFace &firstFaceSide1 = changes.mFacesToAdd.at( edgeFacesIndexes.at( 0 ) );
2486  int pos1InFaceSide1 = vertexPositionInFace( addedVertexIndex, firstFaceSide1 );
2487 
2488  const QgsMeshFace &secondFaceSide1 = changes.mFacesToAdd.at( edgeFacesIndexes.at( 1 ) );
2489  int pos2InFaceSide1 = vertexPositionInFace( addedVertexIndex, secondFaceSide1 );
2490  pos2InFaceSide1 = ( pos2InFaceSide1 + secondFaceSide1.size() - 1 ) % secondFaceSide1.size();
2491 
2492  const QgsMeshFace &firstFaceSide2 = changes.mFacesToAdd.at( edgeFacesIndexesFace2.at( 0 ) );
2493  int pos1InFaceSide2 = vertexPositionInFace( addedVertexIndex, firstFaceSide2 );
2494 
2495  const QgsMeshFace &secondFaceSide2 = changes.mFacesToAdd.at( edgeFacesIndexesFace2.at( 1 ) );
2496  int pos2InFaceSide2 = vertexPositionInFace( addedVertexIndex, secondFaceSide2 );
2497  pos2InFaceSide2 = ( pos2InFaceSide2 + secondFaceSide1.size() - 1 ) % secondFaceSide1.size();
2498 
2499  changes.mFacesNeighborhoodToAdd[edgeFacesIndexes.at( 0 )][pos1InFaceSide1] = edgeFacesIndexesFace2.at( 1 ) + changes.mAddedFacesFirstIndex;
2500  changes.mFacesNeighborhoodToAdd[edgeFacesIndexes.at( 1 )][pos2InFaceSide1] = edgeFacesIndexesFace2.at( 0 ) + changes.mAddedFacesFirstIndex;
2501  changes.mFacesNeighborhoodToAdd[edgeFacesIndexesFace2.at( 0 )][pos1InFaceSide2] = edgeFacesIndexes.at( 1 ) + changes.mAddedFacesFirstIndex;
2502  changes.mFacesNeighborhoodToAdd[edgeFacesIndexesFace2.at( 1 )][pos2InFaceSide2] = edgeFacesIndexes.at( 0 ) + changes.mAddedFacesFirstIndex;
2503  }
2504 
2505  applyChanges( changes );
2506  return changes;
2507 }
2508 
2509 QgsTopologicalMesh::Changes QgsTopologicalMesh::changeZValue( const QList<int> &verticesIndexes, const QList<double> &newValues )
2510 {
2511  Q_ASSERT( verticesIndexes.count() == newValues.count() );
2512  Changes changes;
2513  changes.mChangeCoordinateVerticesIndexes.reserve( verticesIndexes.count() );
2514  changes.mNewZValues.reserve( verticesIndexes.count() );
2515  changes.mOldZValues.reserve( verticesIndexes.count() );
2516  for ( int i = 0; i < verticesIndexes.count(); ++i )
2517  {
2518  changes.mChangeCoordinateVerticesIndexes.append( verticesIndexes.at( i ) );
2519  changes.mOldZValues.append( mMesh->vertices.at( verticesIndexes.at( i ) ).z() );
2520  changes.mNewZValues.append( newValues.at( i ) );
2521  }
2522 
2523  applyChanges( changes );
2524 
2525  return changes;
2526 }
2527 
2528 QgsTopologicalMesh::Changes QgsTopologicalMesh::changeXYValue( const QList<int> &verticesIndexes, const QList<QgsPointXY> &newValues )
2529 {
2530  Q_ASSERT( verticesIndexes.count() == newValues.count() );
2531  Changes changes;
2532  changes.mChangeCoordinateVerticesIndexes.reserve( verticesIndexes.count() );
2533  changes.mNewXYValues.reserve( verticesIndexes.count() );
2534  changes.mOldXYValues.reserve( verticesIndexes.count() );
2535  QSet<int> concernedFace;
2536  for ( int i = 0; i < verticesIndexes.count(); ++i )
2537  {
2538  changes.mChangeCoordinateVerticesIndexes.append( verticesIndexes.at( i ) );
2539  changes.mOldXYValues.append( mMesh->vertices.at( verticesIndexes.at( i ) ) );
2540  changes.mNewXYValues.append( newValues.at( i ) );
2541  concernedFace.unite( qgis::listToSet( facesAroundVertex( verticesIndexes.at( i ) ) ) );
2542  }
2543 
2544  changes.mNativeFacesIndexesGeometryChanged = concernedFace.values();
2545 
2546  applyChanges( changes );
2547 
2548  return changes;
2549 }
@ InvalidFace
An error occurs due to an invalid face (for example, vertex indexes are unordered)
@ UniqueSharedVertex
A least two faces share only one vertices.
@ ManifoldFace
ManifoldFace.
@ InvalidVertex
An error occurs due to an invalid vertex (for example, vertex index is out of range the available ver...
@ FlatFace
A flat face is present.
static bool segmentIntersection(const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &q1, const QgsPoint &q2, QgsPoint &intersectionPoint, bool &isIntersection, double tolerance=1e-8, bool acceptImproperIntersection=false) SIP_HOLDGIL
Compute the intersection between two segments.
Class that represents an error during mesh editing.
Definition: qgsmesheditor.h:43
Qgis::MeshEditingErrorType errorType
Definition: qgsmesheditor.h:52
Convenient class that turn around a vertex and provide information about faces and vertices.
bool isValid() const
Returns whether the vertex circulator is valid.
int turnClockwise() const
Turns counter clockwise around the vertex and returns the new current face, -1 if the circulator pass...
bool goBoundaryCounterClockwise() const
Sets the circulator on the boundary face turning counter clockwise, return false is there isn't bound...
int oppositeVertexCounterClockwise() const
Returns the opposite vertex of the current face and on the edge on the side turning counter clockwise...
int turnCounterClockwise() const
Turns counter clockwise around the vertex and returns the new current face, -1 if the circulator pass...
int currentFaceIndex() const
Returns the current face index, -1 if the circulator has passed a boundary or circulator is invalid.
bool goBoundaryClockwise() const
Sets the circulator on the boundary face turning clockwise, return false is there isn't boundary face...
QgsMeshFace currentFace() const
Returns the current face, empty face if the circulator pass a boundary or circulator is invalid.
QgsMeshVertexCirculator(const QgsTopologicalMesh &topologicalMesh, int vertexIndex)
Constructor with topologicalMesh and vertexIndex.
int oppositeVertexClockwise() const
Returns the opposite vertex of the current face and on the edge on the side turning clockwise.
int degree() const
Returns the degree of the vertex, that is the count of other vertices linked.
QList< int > facesAround() const
Returns all the faces indexes around the vertex.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
A class to represent a 2D point.
Definition: qgspointxy.h:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
bool isEmpty() const override SIP_HOLDGIL
Returns true if the geometry is empty.
Definition: qgspoint.cpp:767
Q_GADGET double x
Definition: qgspoint.h:52
double y
Definition: qgspoint.h:53
Class that contains topological differences between two states of a topological mesh,...
QList< int > mChangeCoordinateVerticesIndexes
void clearChanges()
Clears all changes.
QVector< FaceNeighbors > mFacesNeighborhoodToRemove
QVector< QgsMeshFace > removedFaces() const
Returns the faces that are removed with this changes.
QList< QgsPointXY > mNewXYValues
QList< QgsMeshVertex > mRemovedVertices
QVector< QgsMeshVertex > addedVertices() const
Returns the added vertices with this changes.
QList< std::array< int, 4 > > mNeighborhoodChanges
bool isEmpty() const
Returns whether changes are empty, that there is nothing to change.
QList< int > changedCoordinatesVerticesIndexes() const
Returns the indexes of vertices that have changed coordinates.
QList< int > mNativeFacesIndexesGeometryChanged
QVector< QgsMeshFace > mFacesToAdd
QList< int > removedFaceIndexes() const
Returns the indexes of the faces that are removed with this changes.
QVector< FaceNeighbors > mFacesNeighborhoodToAdd
QList< std::array< int, 3 > > mVerticesToFaceChanges
QList< QgsPointXY > mOldXYValues
QList< double > newVerticesZValues() const
Returns the new Z values of vertices that have changed their coordinates.
QVector< QgsMeshVertex > mVerticesToAdd
QVector< QgsMeshFace > addedFaces() const
Returns the face that are added with this changes.
QList< QgsPointXY > oldVerticesXYValues() const
Returns the old (X,Y) values of vertices that have changed their coordinates.
QList< QgsPointXY > newVerticesXYValues() const
Returns the new (X,Y) values of vertices that have changed their coordinates.
QVector< QgsMeshFace > mFacesToRemove
QList< int > nativeFacesIndexesGeometryChanged() const
Returns a list of the native face indexes that have a geometry changed.
QList< int > verticesToRemoveIndexes() const
Returns the indexes of vertices to remove.
Class that contains independent faces an topological information about this faces.
int vertexToFace(int vertexIndex) const
Returns a face linked to the vertices with index vertexIndex.
QVector< FaceNeighbors > facesNeighborhood() const
Returns the face neighborhood of the faces, indexing is local.
void clear()
Clears all data contained in the instance.
QVector< QgsMeshFace > meshFaces() const
Returns faces.
Class that wraps a QgsMesh to ensure the consistency of the mesh during editing and help to access to...
Changes changeZValue(const QList< int > &verticesIndexes, const QList< double > &newValues)
Changes the Z values of the vertices with indexes in vertices indexes with the values in newValues.
static QgsTopologicalMesh createTopologicalMesh(QgsMesh *mesh, int maxVerticesPerFace, QgsMeshEditingError &error)
Creates a topologicaly consistent mesh with mesh, this static method modifies mesh to be topological ...
bool isVertexFree(int vertexIndex) const
Returns whether the vertex is a free vertex.
static QgsMeshEditingError counterClockwiseFaces(QgsMeshFace &face, QgsMesh *mesh)
Checks the topology of the face and sets it counter clockwise if necessary.
Changes removeVertexFillHole(int vertexIndex)
Removes the vertex with index vertexIndex.
static QgsMeshEditingError checkTopology(const QgsMesh &mesh, int maxVerticesPerFace)
Checks the topology of the mesh mesh, if error occurs, this mesh can't be edited.
friend class QgsMeshVertexCirculator
void applyChanges(const Changes &changes)
Applies the changes.
int firstFaceLinked(int vertexIndex) const
Returns the index of the first face linked, returns -1 if it is a free vertex or out of range index.
QgsMeshEditingError facesCanBeRemoved(const QList< int > facesIndexes)
Returns whether faces with index in faceIndexes can be removed/ The method an error object with type ...
QgsMeshEditingError checkConsistency() const
Checks the consistency of the topological mesh and return false if there is a consistency issue.
Changes removeVertices(const QList< int > &vertices)
Removes all the vertices with index in the list vertices If vertices in linked with faces,...
Changes changeXYValue(const QList< int > &verticesIndexes, const QList< QgsPointXY > &newValues)
Changes the (X,Y) values of the vertices with indexes in vertices indexes with the values in newValue...
void reindex()
Reindexes faces and vertices, after this operation, the topological mesh can't be edited anymore and ...
QVector< int > neighborsOfFace(int faceIndex) const
Returns the indexes of neighbor faces of the face with index faceIndex.
QgsMeshEditingError facesCanBeAdded(const TopologicalFaces &topologicalFaces) const
Returns whether the faces can be added to the mesh.
bool renumber()
Renumbers the indexes of vertices and faces using the Reverse CutHill McKee Algorithm.
Changes flipEdge(int vertexIndex1, int vertexIndex2)
Flips edge (vertexIndex1, vertexIndex2) The method returns a instance of the class QgsTopologicalMesh...
QVector< int > FaceNeighbors
void reverseChanges(const Changes &changes)
Reverses the changes.
Changes addFaces(const TopologicalFaces &topologicFaces)
Adds faces topologicFaces to the topologic mesh.
Changes merge(int vertexIndex1, int vertexIndex2)
Merges faces separated by vertices with indexes vertexIndex1 and vertexIndex2 The method returns a in...
QList< int > freeVerticesIndexes() const
Returns a list of vertices are not linked to any faces.
bool edgeCanBeFlipped(int vertexIndex1, int vertexIndex2) const
Returns true if the edge can be flipped (only available for edge shared by two faces with 3 vertices)
Changes addVertexInFace(int faceIndex, const QgsMeshVertex &vertex)
Adds a vertex in the face with index faceIndex.
Changes removeFaces(const QList< int > facesIndexes)
Removes faces with index in faceIndexes.
bool canBeMerged(int vertexIndex1, int vertexIndex2) const
Returns true if faces separated by vertices with indexes vertexIndex1 and vertexIndex2 can be merged.
QList< int > facesAroundVertex(int vertexIndex) const
Returns the indexes of faces that are around the vertex with index vertexIndex.
bool canBeSplit(int faceIndex) const
Returns true if face with index faceIndex can be split.
Changes addFreeVertex(const QgsMeshVertex &vertex)
Adds a free vertex in the face, that is a vertex tha tis not included or linked with any faces.
Changes insertVertexInFacesEdge(int faceIndex, int position, const QgsMeshVertex &vertex)
Inserts a vertex in the edge of face with index faceIndex at position .
QgsMesh * mesh() const
Returns a pointer to the wrapped mesh.
bool isVertexOnBoundary(int vertexIndex) const
Returns whether the vertex is on a boundary.
Changes splitFace(int faceIndex)
Splits face with index faceIndex The method returns a instance of the class QgsTopologicalMesh::Chang...
static TopologicalFaces createNewTopologicalFaces(const QVector< QgsMeshFace > &faces, bool uniqueSharedVertexAllowed, QgsMeshEditingError &error)
Creates new topological faces that are not yet included in the mesh.
QgsMeshVertexCirculator vertexCirculator(int vertexIndex) const
Returns a vertex circulator linked to this mesh around the vertex with index vertexIndex.
A class to represent a vector.
Definition: qgsvector.h:30
double angle() const SIP_HOLDGIL
Returns the angle of the vector in radians.
Definition: qgsvector.h:172
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
QVector< int > QgsMeshFace
List of vertex indexes.
QgsPoint QgsMeshVertex
xyz coords of vertex
Mesh - vertices, edges and faces.
int vertexCount() const
Returns number of vertices.
QVector< QgsMeshVertex > vertices
QgsMeshFace face(int index) const
Returns a face at the index.
QVector< QgsMeshFace > faces
int faceCount() const
Returns number of faces.
QgsMeshVertex vertex(int index) const
Returns a vertex at the index.