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