QGIS API Documentation  3.23.0-Master (7c4a6de034)
NormVecDecorator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  NormVecDecorator.cpp
3  --------------------
4  copyright : (C) 2004 by Marco Hugentobler
5  email : [email protected]
6  ***************************************************************************/
7 
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 <QApplication>
18 
19 #include "qgsfeedback.h"
20 #include "qgslogger.h"
21 #include "qgsfields.h"
22 #include "qgspoint.h"
23 
24 #include "NormVecDecorator.h"
25 
26 
28 {
29  //remove all the normals
30  if ( !mNormVec->isEmpty() )
31  {
32  for ( int i = 0; i < mNormVec->count(); i++ )
33  {
34  delete ( *mNormVec )[i];
35  }
36  }
37 
38  delete mNormVec;
39  delete mPointState;
40  delete mTIN;
41 }
42 
44 {
45  if ( mTIN )
46  {
47  int pointno;
48  pointno = mTIN->addPoint( p );
49 
50  if ( pointno == -100 )//a numerical error occurred
51  {
52 // QgsDebugMsg("warning, numerical error");
53  return -100;
54  }
55 
56  //recalculate the necessary normals
57  if ( alreadyestimated )
58  {
59  estimateFirstDerivative( pointno );
60  //update also the neighbours of the new point
61  const QList<int> list = mTIN->surroundingTriangles( pointno );
62  auto it = list.constBegin();//iterate through the list and analyze it
63  while ( it != list.constEnd() )
64  {
65  int point;
66  point = *it;
67  if ( point != -1 )
68  {
70  }
71  ++it;
72  ++it;
73  ++it;
74  ++it;
75  }
76  }
77  return pointno;
78  }
79 
80  return -1;
81 }
82 
83 bool NormVecDecorator::calcNormal( double x, double y, QgsPoint &result )
84 {
85  if ( !alreadyestimated )
86  {
88  alreadyestimated = true;
89  }
90 
91  if ( mInterpolator )
92  {
93  const bool b = mInterpolator->calcNormVec( x, y, result );
94  return b;
95  }
96  else
97  {
98  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
99  return false;
100  }
101 }
102 
103 bool NormVecDecorator::calcNormalForPoint( double x, double y, int pointIndex, Vector3D *result )
104 {
105  if ( !alreadyestimated )
106  {
108  alreadyestimated = true;
109  }
110 
111  if ( result )
112  {
113  int numberofbreaks = 0;//number of breaklines around the point
114  int ffirstbp = -1000; //numbers of the points related to the first breakline
115  int lfirstbp = -1000;
116  bool pointfound = false;//is set to true, if the triangle with the point in is found
117  int numberofruns = 0;//number of runs of the loop. This integer can be used to prevent endless loops
118  const int limit = 100000;//ater this number of iterations, the method is terminated
119 
120  result->setX( 0 );
121  result->setY( 0 );
122  result->setZ( 0 );
123 
124  const QList<int> vlist = surroundingTriangles( pointIndex );//get the value list
125  if ( vlist.empty() )//an error occurred in 'getSurroundingTriangles'
126  {
127  return false;
128  }
129 
130  if ( ( ( vlist.count() ) % 4 ) != 0 ) //number of items in vlist has to be a multiple of 4
131  {
132  QgsDebugMsg( QStringLiteral( "warning, wrong number of items in vlist" ) );
133  return false;
134  }
135 
136  auto it = vlist.constBegin();
137 
138  bool firstrun;
139 
140  while ( true )//endless loop to analyze vlist
141  {
142  numberofruns++;
143  if ( numberofruns > limit )
144  {
145  QgsDebugMsg( QStringLiteral( "warning, a probable endless loop is detected" ) );
146  return false;
147  }
148 
149  int p1, p2, p3, line;
150  firstrun = false;//flag which tells, if it is the first run with a breakline
151  p1 = ( *it );
152  ++it;
153  p2 = ( *it );
154  ++it;
155  p3 = ( *it );
156  ++it;
157  line = ( *it );
158 
159 
160  if ( numberofbreaks > 0 )
161  {
162 
163  if ( p1 != -1 && p2 != -1 && p3 != -1 )
164  {
165  if ( MathUtils::pointInsideTriangle( x, y, point( p1 ), point( p2 ), point( p3 ) ) )
166  {
167  pointfound = true;
168  }
169 
170  Vector3D addvec( 0, 0, 0 );
171  MathUtils::normalFromPoints( point( p1 ), point( p2 ), point( p3 ), &addvec );
172  result->setX( result->getX() + addvec.getX() );
173  result->setY( result->getY() + addvec.getY() );
174  result->setZ( result->getZ() + addvec.getZ() );
175  }
176  }
177 
178  if ( line == -10 )//we found a breakline
179  {
180 
181  if ( numberofbreaks == 0 )//it is the first breakline
182  {
183  firstrun = true;
184  ffirstbp = p2;//set the marks to recognize the breakline later
185  lfirstbp = p3;
186  }
187 
188  if ( p2 == ffirstbp && p3 == lfirstbp && !firstrun )//we are back at the first breakline
189  {
190  if ( !pointfound )//the point with coordinates x, y was in no triangle
191  {
192  QgsDebugMsg( QStringLiteral( "warning: point (x,y) was in no triangle" ) );
193  return false;
194  }
195  result->standardise();
196  break;
197  }
198 
199  if ( numberofbreaks > 0 && pointfound )//we found the second break line and the point is between the first and the second
200  {
201  result->standardise();
202  numberofbreaks++;//to make the distinction between endpoints and points on a breakline easier
203  break;
204  }
205 
206  result->setX( 0 );
207  result->setY( 0 );
208  result->setZ( 0 );
209  numberofbreaks++;
210  }
211 
212  ++it;
213  if ( it == vlist.constEnd() )//restart at the beginning of the loop
214  {
215  it = vlist.constBegin();
216  }
217  }
218  return true;
219  }
220  else
221  {
222  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
223  return false;
224  }
225 
226 }
227 
228 bool NormVecDecorator::calcPoint( double x, double y, QgsPoint &result )
229 {
230 
231  if ( !alreadyestimated )
232  {
234  alreadyestimated = true;
235  }
236 
237  if ( mInterpolator )
238  {
239  const bool b = mInterpolator->calcPoint( x, y, result );
240  return b;
241  }
242  else
243  {
244  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
245  return false;
246  }
247 }
248 
249 bool NormVecDecorator::getTriangle( double x, double y, QgsPoint &p1, Vector3D *v1, QgsPoint &p2, Vector3D *v2, QgsPoint &p3, Vector3D *v3 )
250 {
251  if ( v1 && v2 && v3 )
252  {
253  int nr1 = 0;
254  int nr2 = 0;
255  int nr3 = 0;
256 
257  if ( TriDecorator::triangleVertices( x, y, p1, nr1, p2, nr2, p3, nr3 ) )//everything alright
258  {
259  if ( ( *mNormVec )[ nr1 ] && ( *mNormVec )[ nr2 ] && ( *mNormVec )[ nr3 ] )
260  {
261  v1->setX( ( *mNormVec )[ nr1 ]->getX() );
262  v1->setY( ( *mNormVec )[nr1 ]->getY() );
263  v1->setZ( ( *mNormVec )[nr1 ]->getZ() );
264 
265  v2->setX( ( *mNormVec )[nr2 ]->getX() );
266  v2->setY( ( *mNormVec )[nr2 ]->getY() );
267  v2->setZ( ( *mNormVec )[nr2 ]->getZ() );
268 
269  v3->setX( ( *mNormVec )[nr3 ]->getX() );
270  v3->setY( ( *mNormVec )[nr3 ]->getY() );
271  v3->setZ( ( *mNormVec )[nr3 ]->getZ() );
272  }
273  else
274  {
275  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
276  return false;
277  }
278  return true;
279  }
280  else
281  {
282  return false;
283  }
284  }
285 
286  else
287  {
288  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
289  return false;
290  }
291 }
292 
294 {
295  if ( pointno >= 0 )
296  {
297  return mPointState->at( pointno );
298  }
299  else
300  {
301  QgsDebugMsg( QStringLiteral( "warning, number below 0" ) );
302  return mPointState->at( 0 );//just to avoid a compiler warning
303  }
304 }
305 
306 
307 bool NormVecDecorator::getTriangle( double x, double y, QgsPoint &p1, int &ptn1, Vector3D *v1, PointState *state1, QgsPoint &p2, int &ptn2, Vector3D *v2, PointState *state2, QgsPoint &p3, int &ptn3, Vector3D *v3, PointState *state3 )
308 {
309  if ( v1 && v2 && v3 && state1 && state2 && state3 )
310  {
311  if ( TriDecorator::triangleVertices( x, y, p1, ptn1, p2, ptn2, p3, ptn3 ) )//everything alright
312  {
313  v1->setX( ( *mNormVec )[( ptn1 )]->getX() );
314  v1->setY( ( *mNormVec )[( ptn1 )]->getY() );
315  v1->setZ( ( *mNormVec )[( ptn1 )]->getZ() );
316 
317  ( *state1 ) = ( *mPointState )[( ptn1 )];
318 
319  v2->setX( ( *mNormVec )[( ptn2 )]->getX() );
320  v2->setY( ( *mNormVec )[( ptn2 )]->getY() );
321  v2->setZ( ( *mNormVec )[( ptn2 )]->getZ() );
322 
323  ( *state2 ) = ( *mPointState )[( ptn2 )];
324 
325  v3->setX( ( *mNormVec )[( ptn3 )]->getX() );
326  v3->setY( ( *mNormVec )[( ptn3 )]->getY() );
327  v3->setZ( ( *mNormVec )[( ptn3 )]->getZ() );
328 
329  ( *state3 ) = ( *mPointState )[( ptn3 )];
330 
331  return true;
332  }
333  else
334  {
335  return false;
336  }
337 
338  }
339  else
340  {
341  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
342  return false;
343  }
344 }
345 
347 {
348  if ( pointno == -1 )
349  {
350  return false;
351  }
352 
353  Vector3D part;
354  Vector3D total;
355  total.setX( 0 );
356  total.setY( 0 );
357  total.setZ( 0 );
358  int numberofbreaks = 0;//number of counted breaklines
359  double weights = 0;//sum of the weights
360  double currentweight = 0;//current weight
361  PointState status;
362 
363  const QList<int> vlist = surroundingTriangles( pointno );//get the value list
364 
365  if ( vlist.empty() )
366  {
367  //something went wrong in getSurroundingTriangles, set the normal to (0,0,0)
368  if ( mNormVec->size() <= mNormVec->count() )//allocate more memory if necessary
369  {
370  QgsDebugMsg( QStringLiteral( "resizing mNormVec from %1 to %2" ).arg( mNormVec->size() ).arg( mNormVec->size() + 1 ) );
371  mNormVec->resize( mNormVec->size() + 1 );
372  }
373 
374  //todo:resize mNormVec if necessary
375 
376  if ( !( ( *mNormVec )[pointno] ) ) //insert a pointer to a Vector3D, if there is none at this position
377  {
378  Vector3D *vec = new Vector3D( total.getX(), total.getY(), total.getZ() );
379  mNormVec->insert( pointno, vec );
380  }
381  else
382  {
383  ( *mNormVec )[pointno]->setX( 0 );
384  ( *mNormVec )[pointno]->setY( 0 );
385  ( *mNormVec )[pointno]->setZ( 0 );
386  }
387  return false;
388  }
389 
390  if ( ( vlist.count() % 4 ) != 0 ) //number of items in vlist has to be a multiple of 4
391  {
392  QgsDebugMsg( QStringLiteral( "warning, wrong number of items in vlist" ) );
393  return false;
394  }
395 
396  auto it = vlist.constBegin();//iterate through the list and analyze it
397  while ( it != vlist.constEnd() )
398  {
399  int p1, p2, p3, flag;
400  part.setX( 0 );
401  part.setY( 0 );
402  part.setZ( 0 );
403 
404  currentweight = 0;
405 
406  p1 = ( *it );
407  ++it;
408  p2 = ( *it );
409  ++it;
410  p3 = ( *it );
411  ++it;
412  flag = ( *it );
413 
414  if ( flag == -10 )//we found a breakline.
415  {
416  numberofbreaks++;
417  }
418 
419  if ( p1 != -1 && p2 != -1 && p3 != -1 )//don't calculate normal, if a point is a virtual point
420  {
421  MathUtils::normalFromPoints( point( p1 ), point( p2 ), point( p3 ), &part );
422  const double dist1 = point( p3 )->distance3D( *point( p1 ) );
423  const double dist2 = point( p3 )->distance3D( *point( p2 ) );
424  //don't add the normal if the triangle is horizontal
425  if ( ( point( p1 )->z() != point( p2 )->z() ) || ( point( p1 )->z() != point( p3 )->z() ) )
426  {
427  currentweight = 1 / ( dist1 * dist1 * dist2 * dist2 );
428  total.setX( total.getX() + part.getX()*currentweight );
429  total.setY( total.getY() + part.getY()*currentweight );
430  total.setZ( total.getZ() + part.getZ()*currentweight );
431  weights += currentweight;
432  }
433  }
434  ++it;
435  }
436 
437  if ( total.getX() == 0 && total.getY() == 0 && total.getZ() == 0 )//we have a point surrounded by horizontal triangles
438  {
439  total.setZ( 1 );
440  }
441  else
442  {
443  total.setX( total.getX() / weights );
444  total.setY( total.getY() / weights );
445  total.setZ( total.getZ() / weights );
446  total.standardise();
447  }
448 
449 
450  if ( numberofbreaks == 0 )
451  {
452  status = Normal;
453  }
454  else if ( numberofbreaks == 1 )
455  {
456  status = EndPoint;
457  }
458  else
459  {
460  status = BreakLine;
461  }
462 
463  //insert the new calculated vector
464  if ( mNormVec->size() <= mNormVec->count() )//allocate more memory if necessary
465  {
466  mNormVec->resize( mNormVec->size() + 1 );
467  }
468 
469  if ( !( ( *mNormVec )[pointno] ) ) //insert a pointer to a Vector3D, if there is none at this position
470  {
471  Vector3D *vec = new Vector3D( total.getX(), total.getY(), total.getZ() );
472  mNormVec->insert( pointno, vec );
473  }
474  else
475  {
476  ( *mNormVec )[pointno]->setX( total.getX() );
477  ( *mNormVec )[pointno]->setY( total.getY() );
478  ( *mNormVec )[pointno]->setZ( total.getZ() );
479  }
480 
481  //insert the new status
482 
483  if ( pointno >= mPointState->size() )
484  {
485  mPointState->resize( mPointState->size() + 1 );
486  }
487 
488  ( *mPointState )[pointno] = status;
489 
490  return true;
491 }
492 
493 //weighted method of little
495 {
496  const int numberPoints = pointsCount();
497  for ( int i = 0; i < numberPoints; i++ )
498  {
499  if ( feedback )
500  {
501  feedback->setProgress( 100.0 * static_cast< double >( i ) / numberPoints );
502  }
504  }
505 
506  return true;
507 }
508 
510 {
511  if ( mTIN )
512  {
513  if ( alreadyestimated )
514  {
517  }
518  else
519  {
521  }
522  }
523  else
524  {
525  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
526  }
527 }
528 
530 {
531  if ( pointno >= 0 )
532  {
533  ( *mPointState )[pointno] = s;
534  }
535  else
536  {
537  QgsDebugMsg( QStringLiteral( "warning, pointno>0" ) );
538  }
539 }
540 
541 bool NormVecDecorator::swapEdge( double x, double y )
542 {
543  if ( mTIN )
544  {
545  bool b = false;
546  if ( alreadyestimated )
547  {
548  const QList<int> list = pointsAroundEdge( x, y );
549  if ( !list.empty() )
550  {
551  b = mTIN->swapEdge( x, y );
552  for ( const int i : list )
553  {
555  }
556  }
557  }
558  else
559  {
560  b = mTIN->swapEdge( x, y );
561  }
562  return b;
563  }
564  else
565  {
566  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
567  return false;
568  }
569 }
570 
572 {
573  if ( !mTIN )
574  {
575  return false;
576  }
577  return mTIN->saveTriangulation( sink, feedback );
578 }
579 
581 {
582  if ( !mTIN )
583  {
584  return QgsMesh();
585  }
586  return mTIN->triangulationToMesh( feedback );
587 }
588 
bool getTriangle(double x, double y, QgsPoint &p1, Vector3D *v1, QgsPoint &p2, Vector3D *v2, QgsPoint &p3, Vector3D *v3)
Finds out, in which triangle a point with coordinates x and y is and assigns the triangle points to p...
QVector< PointState > * mPointState
Vector who stores, it a point is not on a breakline, if it is a normal point of the breakline or if i...
bool alreadyestimated
Is true, if the normals already have been estimated.
bool calcNormal(double x, double y, QgsPoint &result) override
Calculates the normal at a point on the surface and assigns it to 'result'. Returns true in case of s...
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
~NormVecDecorator() override
QgsMesh triangulationToMesh(QgsFeedback *feedback=nullptr) const override
Returns a QgsMesh corresponding to the triangulation.
void eliminateHorizontalTriangles() override
Eliminates the horizontal triangles by swapping or by insertion of new points. If alreadyestimated is...
bool calcPoint(double x, double y, QgsPoint &result) override
Calculates x-, y and z-value of the point on the surface and assigns it to 'result'.
TriangleInterpolator * mInterpolator
Association with an interpolator object.
void setState(int pointno, PointState s)
Sets the state (BreakLine, Normal, EndPoint) of a point.
QVector< Vector3D * > * mNormVec
Vector that stores the normals for the points. If 'estimateFirstDerivatives()' was called and there i...
bool swapEdge(double x, double y) override
Swaps the edge which is closest to the point with x and y coordinates (if this is possible) and force...
PointState getState(int pointno) const
Returns the state of the point with the number 'pointno'.
bool calcNormalForPoint(double x, double y, int pointIndex, Vector3D *result)
Calculates the normal of a triangle-point for the point with coordinates x and y. This is needed,...
bool estimateFirstDerivatives(QgsFeedback *feedback=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
bool estimateFirstDerivative(int pointno)
Estimates the first derivative a point. Return true in case of success and false otherwise.
PointState
Enumeration for the state of a point. Normal means, that the point is not on a BreakLine,...
int addPoint(const QgsPoint &p) override
Adds a point to the triangulation.
An interface for objects which accept features via addFeature(s) methods.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
double distance3D(double x, double y, double z) const SIP_HOLDGIL
Returns the Cartesian 3D distance between this point and a specified x, y, z coordinate.
Definition: qgspoint.cpp:680
virtual QgsMesh triangulationToMesh(QgsFeedback *feedback=nullptr) const =0
Returns a QgsMesh corresponding to the triangulation.
virtual void eliminateHorizontalTriangles()=0
Eliminates the horizontal triangles by swapping.
virtual QList< int > surroundingTriangles(int pointno)=0
Returns a value list with the information of the triangles surrounding (counterclockwise) a point.
virtual bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const =0
Saves the triangulation features to a feature sink.
virtual int addPoint(const QgsPoint &point)=0
Adds a point to the triangulation.
virtual bool swapEdge(double x, double y)=0
Reads the content of a taff-file.
QgsPoint * point(int i) const override
Returns a pointer to the point with number i.
QgsTriangulation * mTIN
Association with a Triangulation object.
Definition: TriDecorator.h:65
QList< int > pointsAroundEdge(double x, double y) override
Returns a value list with the numbers of the four points, which would be affected by an edge swap.
QList< int > surroundingTriangles(int pointno) override
Returns a value list with the information of the triangles surrounding (counterclockwise) a point.
int pointsCount() const override
Returns the number of points.
bool triangleVertices(double x, double y, QgsPoint &p1, int &n1, QgsPoint &p2, int &n2, QgsPoint &p3, int &n3) override
Finds out in which triangle the point with coordinates x and y is and assigns the numbers of the vert...
virtual bool calcNormVec(double x, double y, QgsPoint &result)=0
Calculates the normal vector and assigns it to vec.
virtual bool calcPoint(double x, double y, QgsPoint &result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point.
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values.
Definition: Vector3D.h:36
void standardise()
Standardises the vector.
Definition: Vector3D.cpp:24
void setX(double x)
Sets the x-component of the vector.
Definition: Vector3D.h:106
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:96
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:91
void setY(double y)
Sets the y-component of the vector.
Definition: Vector3D.h:111
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:101
void setZ(double z)
Sets the z-component of the vector.
Definition: Vector3D.h:116
bool ANALYSIS_EXPORT pointInsideTriangle(double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3)
Returns true, if the point with coordinates x and y is inside (or at the edge) of the triangle p1,...
Definition: MathUtils.cpp:680
void ANALYSIS_EXPORT normalFromPoints(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec)
Calculates the normal vector of the plane through the points p1, p2 and p3 and assigns the result to ...
Definition: MathUtils.cpp:632
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
Mesh - vertices, edges and faces.