QGIS API Documentation  3.2.0-Bonn (bc43194)
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 "NormVecDecorator.h"
18 #include "qgsfeedback.h"
19 #include "qgslogger.h"
20 #include "qgsfields.h"
21 #include <QApplication>
22 
24 {
25  //remove all the normals
26  if ( !mNormVec->isEmpty() )
27  {
28  for ( int i = 0; i < mNormVec->count(); i++ )
29  {
30  delete ( *mNormVec )[i];
31  }
32  }
33 
34  delete mNormVec;
35  delete mPointState;
36 
37  if ( mTIN )
38  {
39  delete mTIN;
40  }
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->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 {
85  if ( !alreadyestimated )
86  {
88  alreadyestimated = true;
89  }
90 
91  if ( mInterpolator )
92  {
93  bool b = mInterpolator->calcNormVec( x, y, result );
94  return b;
95  }
96  else
97  {
98  QgsDebugMsg( "warning, null pointer" );
99  return false;
100  }
101 }
102 
103 bool NormVecDecorator::calcNormalForPoint( double x, double y, int point, 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  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( "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( "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( "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( "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  bool b = mInterpolator->calcPoint( x, y, result );
240  return b;
241  }
242  else
243  {
244  QgsDebugMsg( "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( "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( "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( "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( "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( QString( "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( "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  {
513  if ( alreadyestimated )
514  {
517  }
518  else
519  {
521  }
522  }
523  else
524  {
525  QgsDebugMsg( "warning, null pointer" );
526  }
527 }
528 
530 {
531  if ( pointno >= 0 )
532  {
533  ( *mPointState )[pointno] = s;
534  }
535  else
536  {
537  QgsDebugMsg( "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  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( "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.
PointState getState(int pointno) const
Returns the state of the point with the number &#39;pointno&#39;.
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
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:98
Base class for feedback objects to be used for cancelation of something running in a worker thread...
Definition: qgsfeedback.h:44
int addPoint(const QgsPoint &p) override
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...
virtual bool calcPoint(double x, double y, QgsPoint &result SIP_OUT)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
QVector< Vector3D * > * mNormVec
Vector that stores the normals for the points. If &#39;estimateFirstDerivatives()&#39; was called and there i...
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.
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...
virtual int addPoint(const QgsPoint &point)=0
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.
bool alreadyestimated
Is true, if the normals already have been estimated.
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:93
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...
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...
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:88
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
virtual void eliminateHorizontalTriangles()=0
Eliminates the horizontal triangles by swapping.
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...
virtual bool calcNormVec(double x, double y, Vector3D *result SIP_OUT)=0
Calculates the normal vector and assigns it to vec.
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:576