← Back to branch summary

~ubuntu-branches/ubuntu/wily/qgis/wily

« back to all changes in this revision

Viewing changes to src/analysis/interpolation/NormVecDecorator.cc

  • Committer: Bazaar Package Importer
  • Author(s): Johan Van de Wauw
  • Date: 2010-07-11 20:23:24 UTC
  • mfrom: (3.1.4 squeeze)
  • Revision ID: james.westby@ubuntu.com-20100711202324-5ktghxa7hracohmr
Tags: 1.4.0+12730-3ubuntu1
https://launchpad.net/bugs/540941
* Merge from Debian unstable (LP: #540941).
* Fix compilation issues with QT 4.7
* Add build-depends on libqt4-webkit-dev 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/analysis/interpolation/NormVecDecorator.cc
 
1
/***************************************************************************
 
2
                          NormVecDecorator.cc  -  description
 
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 "NormVecDecorator.h"
 
18
#include "qgslogger.h"
 
19
#include <QApplication>
 
20
#include <QProgressDialog>
 
21
 
 
22
NormVecDecorator::~NormVecDecorator()
 
23
{
 
24
  //remove all the normals
 
25
  if ( mNormVec->count() > 0 )
 
26
  {
 
27
    for ( int i = 0; i < mNormVec->count(); i++ )
 
28
    {
 
29
      delete( *mNormVec )[i];
 
30
    }
 
31
  }
 
32
 
 
33
  delete mNormVec;
 
34
  delete mPointState;
 
35
 
 
36
  if ( mTIN )
 
37
  {
 
38
    delete mTIN;
 
39
  }
 
40
}
 
41
 
 
42
int NormVecDecorator::addPoint( Point3D* p )
 
43
{
 
44
  if ( mTIN )
 
45
  {
 
46
    int pointno;
 
47
    pointno = mTIN->addPoint( p );
 
48
 
 
49
    if ( pointno == -100 )//a numerical error occured
 
50
    {
 
51
// QgsDebugMsg("warning, numerical error");
 
52
      return -100;
 
53
    }
 
54
 
 
55
    //recalculate the necessary normals
 
56
    if ( alreadyestimated )
 
57
    {
 
58
      estimateFirstDerivative( pointno );
 
59
      //update also the neighbours of the new point
 
60
      QList<int>* list = mTIN->getSurroundingTriangles( pointno );
 
61
      QList<int>::iterator it = list->begin();//iterate through the list and analize it
 
62
      while ( it != list->end() )
 
63
      {
 
64
        int point;
 
65
        point = ( *it );
 
66
        if ( point != -1 )
 
67
        {
 
68
          estimateFirstDerivative( point );
 
69
        }
 
70
        it++;
 
71
        it++;
 
72
        it++;
 
73
        it++;
 
74
      }
 
75
      delete list;
 
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 == false )
 
86
  {
 
87
    estimateFirstDerivatives();
 
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 == false )
 
106
  {
 
107
    estimateFirstDerivatives();
 
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
    QList<int>* vlist = getSurroundingTriangles( point );//get the value list
 
125
 
 
126
    if ( !vlist )//an error occured in 'getSurroundingTriangles'
 
127
    {
 
128
      return false;
 
129
    }
 
130
 
 
131
    if ((( vlist->count() ) % 4 ) != 0 )//number of items in vlist has to be a multiple of 4
 
132
    {
 
133
      QgsDebugMsg( "warning, wrong number of items in vlist" );
 
134
      return false;
 
135
    }
 
136
 
 
137
    QList<int>::iterator it = vlist->begin();
 
138
 
 
139
    bool firstrun;
 
140
 
 
141
    while ( true )//endless loop to analyze vlist
 
142
    {
 
143
      numberofruns++;
 
144
      if ( numberofruns > limit )
 
145
      {
 
146
        QgsDebugMsg( "warning, a probable endless loop is detected" );
 
147
        return false;
 
148
      }
 
149
 
 
150
      int p1, p2, p3, line;
 
151
      firstrun = false;//flag which tells, if it is the first run with a breakline
 
152
      p1 = ( *it );
 
153
      ++it;
 
154
      p2 = ( *it );
 
155
      ++it;
 
156
      p3 = ( *it );
 
157
      ++it;
 
158
      line = ( *it );
 
159
 
 
160
 
 
161
      if ( numberofbreaks > 0 )
 
162
      {
 
163
 
 
164
        if ( p1 != -1 && p2 != -1 && p2 != -1 )
 
165
        {
 
166
          if ( MathUtils::pointInsideTriangle( x, y, getPoint( p1 ), getPoint( p2 ), getPoint( p3 ) ) )
 
167
          {
 
168
            pointfound = true;
 
169
          }
 
170
 
 
171
          Vector3D addvec( 0, 0, 0 );
 
172
          MathUtils::normalFromPoints( getPoint( p1 ), getPoint( p2 ), getPoint( p3 ), &addvec );
 
173
          result->setX( result->getX() + addvec.getX() );
 
174
          result->setY( result->getY() + addvec.getY() );
 
175
          result->setZ( result->getZ() + addvec.getZ() );
 
176
        }
 
177
      }
 
178
 
 
179
      if ( line == -10 )//we found a breakline
 
180
      {
 
181
 
 
182
        if ( numberofbreaks == 0 )//it is the first breakline
 
183
        {
 
184
          firstrun = true;
 
185
          ffirstbp = p2;//set the marks to recognize the breakline later
 
186
          lfirstbp = p3;
 
187
        }
 
188
 
 
189
        if ( p2 == ffirstbp && p3 == lfirstbp && firstrun == false )//we are back at the first breakline
 
190
        {
 
191
          if ( pointfound == false )//the point with coordinates x, y was in no triangle
 
192
          {
 
193
            QgsDebugMsg( "warning: point (x,y) was in no triangle" );
 
194
            return false;
 
195
          }
 
196
          result->standardise();
 
197
          break;
 
198
        }
 
199
 
 
200
        if ( numberofbreaks > 0 && pointfound == true )//we found the second break line and the point is between the first and the second
 
201
        {
 
202
          result->standardise();
 
203
          numberofbreaks++;//to make the distinction between endpoints and points on a breakline easier
 
204
          break;
 
205
        }
 
206
 
 
207
        result->setX( 0 );
 
208
        result->setY( 0 );
 
209
        result->setZ( 0 );
 
210
        numberofbreaks++;
 
211
      }
 
212
 
 
213
      ++it;
 
214
      if ( it == vlist->end() )//restart at the beginning of the loop
 
215
        {it = vlist->begin();}
 
216
 
 
217
 
 
218
    }
 
219
 
 
220
    delete vlist;
 
221
 
 
222
    return true;
 
223
 
 
224
  }
 
225
 
 
226
  else
 
227
  {
 
228
    QgsDebugMsg( "warning, null pointer" );
 
229
    return false;
 
230
  }
 
231
 
 
232
}
 
233
 
 
234
bool NormVecDecorator::calcPoint( double x, double y, Point3D* result )
 
235
{
 
236
 
 
237
  if ( alreadyestimated == false )
 
238
  {
 
239
    estimateFirstDerivatives();
 
240
    alreadyestimated = true;
 
241
  }
 
242
 
 
243
  if ( mInterpolator )
 
244
  {
 
245
    bool b = mInterpolator->calcPoint( x, y, result );
 
246
    return b;
 
247
  }
 
248
  else
 
249
  {
 
250
    QgsDebugMsg( "warning, null pointer" );
 
251
    return false;
 
252
  }
 
253
}
 
254
 
 
255
bool NormVecDecorator::getTriangle( double x, double y, Point3D* p1, Vector3D* v1, Point3D* p2, Vector3D* v2, Point3D* p3, Vector3D* v3 )
 
256
{
 
257
  if ( p1 && p2 && p3 && v1 && v2 && v3 )
 
258
  {
 
259
    int* nr1 = new int();
 
260
    int* nr2 = new int();
 
261
    int* nr3 = new int();
 
262
 
 
263
    if ( TriDecorator::getTriangle( x, y, p1, nr1, p2, nr2, p3, nr3 ) )//everything allright
 
264
    {
 
265
      if (( *mNormVec )[( *nr1 )] && ( *mNormVec )[( *nr2 )] && ( *mNormVec )[( *nr3 )] )
 
266
      {
 
267
        v1->setX(( *mNormVec )[( *nr1 )]->getX() );
 
268
        v1->setY(( *mNormVec )[( *nr1 )]->getY() );
 
269
        v1->setZ(( *mNormVec )[( *nr1 )]->getZ() );
 
270
 
 
271
        v2->setX(( *mNormVec )[( *nr2 )]->getX() );
 
272
        v2->setY(( *mNormVec )[( *nr2 )]->getY() );
 
273
        v2->setZ(( *mNormVec )[( *nr2 )]->getZ() );
 
274
 
 
275
        v3->setX(( *mNormVec )[( *nr3 )]->getX() );
 
276
        v3->setY(( *mNormVec )[( *nr3 )]->getY() );
 
277
        v3->setZ(( *mNormVec )[( *nr3 )]->getZ() );
 
278
      }
 
279
      else
 
280
      {
 
281
        QgsDebugMsg( "warning, null pointer" );
 
282
        delete nr1;
 
283
        delete nr2;
 
284
        delete nr3;
 
285
        return false;
 
286
      }
 
287
 
 
288
      delete nr1;
 
289
      delete nr2;
 
290
      delete nr3;
 
291
      return true;
 
292
 
 
293
    }
 
294
 
 
295
    else
 
296
    {
 
297
      delete nr1;
 
298
      delete nr2;
 
299
      delete nr3;
 
300
      return false;
 
301
    }
 
302
  }
 
303
 
 
304
  else
 
305
  {
 
306
    QgsDebugMsg( "warning, null pointer" );
 
307
    return false;
 
308
  }
 
309
}
 
310
 
 
311
NormVecDecorator::pointState NormVecDecorator::getState( int pointno ) const
 
312
{
 
313
  if ( pointno >= 0 )
 
314
  {
 
315
    return mPointState->at( pointno );
 
316
  }
 
317
  else
 
318
  {
 
319
    QgsDebugMsg( "warning, number below 0" );
 
320
    return mPointState->at( 0 );//just to avoid a compiler warning
 
321
  }
 
322
}
 
323
 
 
324
 
 
325
bool NormVecDecorator::getTriangle( double x, double y, Point3D* p1, int* ptn1, Vector3D* v1, pointState* state1, Point3D* p2, int* ptn2, Vector3D* v2, pointState* state2, Point3D* p3, int* ptn3, Vector3D* v3, pointState* state3 )
 
326
{
 
327
  if ( p1 && p2 && p3 && v1 && v2 && v3 && ptn1 && ptn2 && ptn3 && state1 && state2 && state3 )
 
328
  {
 
329
    if ( TriDecorator::getTriangle( x, y, p1, ptn1, p2, ptn2, p3, ptn3 ) )//everything allright
 
330
    {
 
331
      v1->setX(( *mNormVec )[( *ptn1 )]->getX() );
 
332
      v1->setY(( *mNormVec )[( *ptn1 )]->getY() );
 
333
      v1->setZ(( *mNormVec )[( *ptn1 )]->getZ() );
 
334
 
 
335
      ( *state1 ) = ( *mPointState )[( *ptn1 )];
 
336
 
 
337
      v2->setX(( *mNormVec )[( *ptn2 )]->getX() );
 
338
      v2->setY(( *mNormVec )[( *ptn2 )]->getY() );
 
339
      v2->setZ(( *mNormVec )[( *ptn2 )]->getZ() );
 
340
 
 
341
      ( *state2 ) = ( *mPointState )[( *ptn2 )];
 
342
 
 
343
      v3->setX(( *mNormVec )[( *ptn3 )]->getX() );
 
344
      v3->setY(( *mNormVec )[( *ptn3 )]->getY() );
 
345
      v3->setZ(( *mNormVec )[( *ptn3 )]->getZ() );
 
346
 
 
347
      ( *state3 ) = ( *mPointState )[( *ptn3 )];
 
348
 
 
349
      return true;
 
350
    }
 
351
    else
 
352
    {
 
353
      QgsDebugMsg( "warning, getTriangle returned false" );
 
354
      return false;
 
355
    }
 
356
 
 
357
  }
 
358
  else
 
359
  {
 
360
    QgsDebugMsg( "warning, null pointer" );
 
361
    return false;
 
362
  }
 
363
}
 
364
 
 
365
bool NormVecDecorator::estimateFirstDerivative( int pointno )
 
366
{
 
367
  if ( pointno == -1 )
 
368
  {
 
369
    return false;
 
370
  }
 
371
 
 
372
  Vector3D part;
 
373
  Vector3D total;
 
374
  total.setX( 0 );
 
375
  total.setY( 0 );
 
376
  total.setZ( 0 );
 
377
  int numberofbreaks = 0;//number of counted breaklines
 
378
  double weights = 0;//sum of the weights
 
379
  double currentweight = 0;//current weight
 
380
  pointState status;
 
381
 
 
382
  QList<int>* vlist = getSurroundingTriangles( pointno );//get the value list
 
383
 
 
384
  if ( !vlist )
 
385
  {
 
386
    //something went wrong in getSurroundingTriangles, set the normal to (0,0,0)
 
387
    if ( mNormVec->size() <= mNormVec->count() )//allocate more memory if necessary
 
388
    {
 
389
      QgsDebugMsg( QString( "resizing mNormVec from %1 to %2" ).arg( mNormVec->size() ).arg( mNormVec->size() + 1 ) );
 
390
      mNormVec->resize( mNormVec->size() + 1 );
 
391
    }
 
392
 
 
393
    //todo:resize mNormVec if necessary
 
394
 
 
395
    if ( !(( *mNormVec )[pointno] ) )//insert a pointer to a Vector3D, if there is none at this position
 
396
    {
 
397
      Vector3D* vec = new Vector3D( total.getX(), total.getY(), total.getZ() );
 
398
      mNormVec->insert( pointno, vec );
 
399
    }
 
400
    else
 
401
    {
 
402
      ( *mNormVec )[pointno]->setX( 0 );
 
403
      ( *mNormVec )[pointno]->setY( 0 );
 
404
      ( *mNormVec )[pointno]->setZ( 0 );
 
405
    }
 
406
    return false;
 
407
  }
 
408
 
 
409
  if (( vlist->count() % 4 ) != 0 ) //number of items in vlist has to be a multiple of 4
 
410
  {
 
411
    QgsDebugMsg( "warning, wrong number of items in vlist" );
 
412
    return false;
 
413
  }
 
414
 
 
415
  QList<int>::iterator it = vlist->begin();//iterate through the list and analize it
 
416
  while ( it != vlist->end() )
 
417
  {
 
418
    int p1, p2, p3, flag;
 
419
    part.setX( 0 );
 
420
    part.setY( 0 );
 
421
    part.setZ( 0 );
 
422
 
 
423
    currentweight = 0;
 
424
 
 
425
    p1 = ( *it );
 
426
    ++it;
 
427
    p2 = ( *it );
 
428
    ++it;
 
429
    p3 = ( *it );
 
430
    ++it;
 
431
    flag = ( *it );
 
432
 
 
433
    if ( flag == -10 )//we found a breakline.
 
434
    {
 
435
      numberofbreaks++;
 
436
    }
 
437
 
 
438
    if ( p1 != -1 && p2 != -1 && p3 != -1 )//don't calculate normal, if a point is a virtual point
 
439
    {
 
440
      MathUtils::normalFromPoints( getPoint( p1 ), getPoint( p2 ), getPoint( p3 ), &part );
 
441
      double dist1 = getPoint( p3 )->dist3D( getPoint( p1 ) );
 
442
      double dist2 = getPoint( p3 )->dist3D( getPoint( p2 ) );
 
443
      //don't add the normal if the triangle is horizontal
 
444
      if (( getPoint( p1 )->getZ() != getPoint( p2 )->getZ() ) || ( getPoint( p1 )->getZ() != getPoint( p3 )->getZ() ) )
 
445
      {
 
446
        currentweight = 1 / ( dist1 * dist1 * dist2 * dist2 );
 
447
        total.setX( total.getX() + part.getX()*currentweight );
 
448
        total.setY( total.getY() + part.getY()*currentweight );
 
449
        total.setZ( total.getZ() + part.getZ()*currentweight );
 
450
        weights += currentweight;
 
451
      }
 
452
    }
 
453
    ++it;
 
454
  }
 
455
 
 
456
  if ( total.getX() == 0 && total.getY() == 0 && total.getZ() == 0 )//we have a point surrounded by horizontal triangles
 
457
  {
 
458
    total.setZ( 1 );
 
459
  }
 
460
  else
 
461
  {
 
462
    total.setX( total.getX() / weights );
 
463
    total.setY( total.getY() / weights );
 
464
    total.setZ( total.getZ() / weights );
 
465
    total.standardise();
 
466
  }
 
467
 
 
468
 
 
469
  if ( numberofbreaks == 0 )
 
470
  {
 
471
    status = NORMAL;
 
472
  }
 
473
  else if ( numberofbreaks == 1 )
 
474
  {
 
475
    status = ENDPOINT;
 
476
  }
 
477
  else
 
478
  {
 
479
    status = BREAKLINE;
 
480
  }
 
481
 
 
482
  delete vlist;
 
483
 
 
484
  //insert the new calculated vector
 
485
  if ( mNormVec->size() <= mNormVec->count() )//allocate more memory if necessary
 
486
  {
 
487
    mNormVec->resize( mNormVec->size() + 1 );
 
488
  }
 
489
 
 
490
  if ( !(( *mNormVec )[pointno] ) )//insert a pointer to a Vector3D, if there is none at this position
 
491
  {
 
492
    Vector3D* vec = new Vector3D( total.getX(), total.getY(), total.getZ() );
 
493
    mNormVec->insert( pointno, vec );
 
494
  }
 
495
  else
 
496
  {
 
497
    ( *mNormVec )[pointno]->setX( total.getX() );
 
498
    ( *mNormVec )[pointno]->setY( total.getY() );
 
499
    ( *mNormVec )[pointno]->setZ( total.getZ() );
 
500
  }
 
501
 
 
502
  //insert the new status
 
503
 
 
504
  if ( pointno >= mPointState->size() )
 
505
  {
 
506
    mPointState->resize( mPointState->size() + 1 );
 
507
  }
 
508
 
 
509
  ( *mPointState )[pointno] = status;
 
510
 
 
511
  return true;
 
512
}
 
513
 
 
514
//weighted method of little
 
515
bool NormVecDecorator::estimateFirstDerivatives( QProgressDialog* d )
 
516
{
 
517
  if ( d )
 
518
  {
 
519
    d->setMinimum( 0 );
 
520
    d->setMaximum( getNumberOfPoints() );
 
521
    d->setCancelButton( 0 ); //we cannot cancel derivative estimation
 
522
    d->show();
 
523
  }
 
524
 
 
525
  for ( int i = 0; i < getNumberOfPoints(); i++ )
 
526
  {
 
527
    if ( d )
 
528
    {
 
529
      d->setValue( i );
 
530
    }
 
531
    estimateFirstDerivative( i );
 
532
  }
 
533
 
 
534
  if ( d )
 
535
  {
 
536
    d->setValue( getNumberOfPoints() );
 
537
  }
 
538
  return true;
 
539
}
 
540
 
 
541
void NormVecDecorator::eliminateHorizontalTriangles()
 
542
{
 
543
  if ( mTIN )
 
544
  {
 
545
    if ( alreadyestimated )
 
546
    {
 
547
      mTIN->eliminateHorizontalTriangles();
 
548
      estimateFirstDerivatives();
 
549
    }
 
550
    else
 
551
    {
 
552
      mTIN->eliminateHorizontalTriangles();
 
553
    }
 
554
  }
 
555
  else
 
556
  {
 
557
    QgsDebugMsg( "warning, null pointer" );
 
558
  }
 
559
}
 
560
 
 
561
void NormVecDecorator::setState( int pointno, pointState s )
 
562
{
 
563
  if ( pointno >= 0 )
 
564
  {
 
565
    ( *mPointState )[pointno] = s;
 
566
  }
 
567
  else
 
568
  {
 
569
    QgsDebugMsg( "warning, pointno>0" );
 
570
  }
 
571
}
 
572
 
 
573
bool NormVecDecorator::swapEdge( double x, double y )
 
574
{
 
575
  if ( mTIN )
 
576
  {
 
577
    bool b = false;
 
578
    if ( alreadyestimated )
 
579
    {
 
580
      QList<int>* list = getPointsAroundEdge( x, y );
 
581
      if ( list )
 
582
      {
 
583
        b = mTIN->swapEdge( x, y );
 
584
        QList<int>::iterator it;
 
585
        for ( it = list->begin(); it != list->end(); ++it )
 
586
        {
 
587
          estimateFirstDerivative(( *it ) );
 
588
        }
 
589
        delete list;
 
590
      }
 
591
    }
 
592
    else
 
593
    {
 
594
      b = mTIN->swapEdge( x, y );
 
595
    }
 
596
    return b;
 
597
  }
 
598
  else
 
599
  {
 
600
    QgsDebugMsg( "warning, null pointer" );
 
601
    return false;
 
602
  }
 
603
}
 
604
 
 
605
bool NormVecDecorator::saveAsShapefile( const QString& fileName ) const
 
606
{
 
607
  if ( !mTIN )
 
608
  {
 
609
    return false;
 
610
  }
 
611
  return mTIN->saveAsShapefile( fileName );
 
612
}
 
613
 
 
614