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 *
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
58  {
59  estimateFirstDerivative( pointno );
60  //update also the neighbours of the new point
61  const QList<int> list = mTIN->getSurroundingTriangles( 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  {
69  estimateFirstDerivative( point );
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, Vector3D *result )
84 {
86  {
89  }
90
91  if ( mInterpolator )
92  {
93  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 point, Vector3D *result )
104 {
106  {
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  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 = getSurroundingTriangles( point );//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, getPoint( p1 ), getPoint( p2 ), getPoint( p3 ) ) )
166  {
167  pointfound = true;
168  }
169
170  Vector3D addvec( 0, 0, 0 );
171  MathUtils::normalFromPoints( getPoint( p1 ), getPoint( p2 ), getPoint( 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
232  {
235  }
236
237  if ( mInterpolator )
238  {
239  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::getTriangle( x, y, p1, nr1, p2, nr2, p3, nr3 ) )//everything allright
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::getTriangle( x, y, p1, ptn1, p2, ptn2, p3, ptn3 ) )//everything allright
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 = getSurroundingTriangles( 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( getPoint( p1 ), getPoint( p2 ), getPoint( p3 ), &part );
422  double dist1 = getPoint( p3 )->distance3D( *getPoint( p1 ) );
423  double dist2 = getPoint( p3 )->distance3D( *getPoint( p2 ) );
424  //don't add the normal if the triangle is horizontal
425  if ( ( getPoint( p1 )->z() != getPoint( p2 )->z() ) || ( getPoint( p1 )->z() != getPoint( 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  int numberPoints = getNumberOfPoints();
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  {
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;
547  {
548  QList<int> *list = getPointsAroundEdge( x, y );
549  if ( list )
550  {
551  b = mTIN->swapEdge( x, y );
552  QList<int>::iterator it;
553  for ( it = list->begin(); it != list->end(); ++it )
554  {
555  estimateFirstDerivative( ( *it ) );
556  }
557  delete list;
558  }
559  }
560  else
561  {
562  b = mTIN->swapEdge( x, y );
563  }
564  return b;
565  }
566  else
567  {
568  QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
569  return false;
570  }
571 }
572
574 {
575  if ( !mTIN )
576  {
577  return false;
578  }
579  return mTIN->saveTriangulation( sink, feedback );
580 }
581
TriangleInterpolator * mInterpolator
Association with an interpolator object.
void setX(double x)
Sets the x-component of the vector.
Definition: Vector3D.h:103
bool calcNormal(double x, double y, Vector3D *result) override
Calculates the normal at a point on the surface and assigns it to &#39;result&#39;. Returns true in case of s...
virtual bool swapEdge(double x, double y)=0
Reads the content of a taff-file.
Triangulation * mTIN
Association with a Triangulation object.
Definition: TriDecorator.h:67
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 &#39;result&#39;.
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QList< int > getSurroundingTriangles(int pointno) override
Returns a pointer to a value list with the information of the triangles surrounding (counterclockwise...
void standardise()
Standardises the vector.
Definition: Vector3D.cpp:24
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
void setState(int pointno, PointState s)
Sets the state (BreakLine, Normal, EndPoint) of a point.
An interface for objects which accept features via addFeature(s) methods.
QList< int > * getPointsAroundEdge(double x, double y) override
Returns a value list with the numbers of the four points, which would be affected by an edge swap...
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...
bool estimateFirstDerivative(int pointno)
Estimates the first derivative a point. Return true in case of success and false otherwise.
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
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:44
Adds a point to the triangulation.
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values...
Definition: Vector3D.h:33
virtual bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const =0
Saves the triangulation features to a feature sink.
void setZ(double z)
Sets the z-component of the vector.
Definition: Vector3D.h:113
bool getTriangle(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...
QVector< Vector3D * > * mNormVec
Vector that stores the normals for the points. If &#39;estimateFirstDerivatives()&#39; was called and there i...
PointState getState(int pointno) const
Returns the state of the point with the number &#39;pointno&#39;.
int getNumberOfPoints() const override
Returns the number of points.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
bool saveTriangulation(QgsFeatureSink *sink, QgsFeedback *feedback=nullptr) const override
Saves the triangulation features to a feature sink.
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:93
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...
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:88
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...
Adds a point to the triangulation.
PointState
Enumeration for the state of a point. Normal means, that the point is not on a BreakLine, BreakLine means that the point is on a breakline (but not an end point of it) and EndPoint means, that it is an endpoint of a breakline.
Is true, if the normals already have been estimated.
void setY(double y)
Sets the y-component of the vector.
Definition: Vector3D.h:108
virtual QList< int > getSurroundingTriangles(int pointno)=0
Returns a pointer to a value list with the information of the triangles surrounding (counterclockwise...
virtual bool calcNormVec(double x, double y, Vector3D *result)=0
Calculates the normal vector and assigns it to vec.
void eliminateHorizontalTriangles() override
Eliminates the horizontal triangles by swapping or by insertion of new points. If alreadyestimated is...
~NormVecDecorator() override
bool calcNormalForPoint(double x, double y, int point, Vector3D *result)
Calculates the normal of a triangle-point for the point with coordinates x and y. This is needed...
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
double distance3D(double x, double y, double z) const
Returns the 3D distance between this point and a specified x, y, z coordinate.
Definition: qgspoint.cpp:588
virtual void eliminateHorizontalTriangles()=0
Eliminates the horizontal triangles by swapping.
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:98
QgsPoint * getPoint(int i) const override
Returns a pointer to the point with number i. Any virtual points must have the number -1...
bool estimateFirstDerivatives(QgsFeedback *feedback=nullptr)
This method adds the functionality of estimating normals at the data points. Return true in the case ...
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...