← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.viewshed/rbbst.cpp

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.viewshed/rbbst.cpp
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.viewshed
 
5
 *
 
6
 * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
 
7
 *               Yi Zhuang - yzhuang@bowdoin.edu
 
8
 
 
9
 *               Ported to GRASS by William Richard -
 
10
 *               wkrichar@bowdoin.edu or willster3021@gmail.com
 
11
 *               Markus Metz: surface interpolation
 
12
 *
 
13
 * Date:         april 2011 
 
14
 * 
 
15
 * PURPOSE: To calculate the viewshed (the visible cells in the
 
16
 * raster) for the given viewpoint (observer) location.  The
 
17
 * visibility model is the following: Two points in the raster are
 
18
 * considered visible to each other if the cells where they belong are
 
19
 * visible to each other.  Two cells are visible to each other if the
 
20
 * line-of-sight that connects their centers does not intersect the
 
21
 * terrain. The terrain is NOT viewed as a tesselation of flat cells, 
 
22
 * i.e. if the line-of-sight does not pass through the cell center, 
 
23
 * elevation is determined using bilinear interpolation.
 
24
 * The viewshed algorithm is efficient both in
 
25
 * terms of CPU operations and I/O operations. It has worst-case
 
26
 * complexity O(n lg n) in the RAM model and O(sort(n)) in the
 
27
 * I/O-model.  For the algorithm and all the other details see the
 
28
 * paper: "Computing Visibility on * Terrains in External Memory" by
 
29
 * Herman Haverkort, Laura Toma and Yi Zhuang.
 
30
 *
 
31
 * COPYRIGHT: (C) 2008 by the GRASS Development Team
 
32
 *
 
33
 * This program is free software under the GNU General Public License
 
34
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 
35
 *
 
36
 *****************************************************************************/
 
37
 
 
38
/*
 
39
 
 
40
   A R/B BST. Always call initNILnode() before using the tree.
 
41
   Version 0.0.0
 
42
 
 
43
   Version 0.0.1
 
44
   Rewrote BST Deletion to improve efficiency
 
45
 
 
46
   Version 0.0.2
 
47
   Bug fixed in deletion.
 
48
   CLRS pseudocode forgot to make sure that x is not NIL before
 
49
   calling rbDeleteFixup(root,x).
 
50
 
 
51
   Version 0.0.3
 
52
   Some Cleanup. Separated the public portion and the 
 
53
   private porthion of the interface in the header
 
54
 
 
55
 
 
56
   =================================
 
57
   This is based on BST 1.0.4
 
58
   BST change log
 
59
   <---------------->
 
60
   find max is implemented in this version.
 
61
   Version 1.0.2
 
62
 
 
63
   Version 1.0.4 
 
64
   Major bug fix in deletion (when the node has two children, 
 
65
   one of them has a wrong parent pointer after the rotation in the deletion.)
 
66
   <----------------->
 
67
 */
 
68
 
 
69
 
 
70
#include <stdlib.h>
 
71
#include <stdio.h>
 
72
#include <assert.h>
 
73
#include <math.h>
 
74
#include "rbbst.h"
 
75
 
 
76
extern "C"
 
77
{
 
78
#include <grass/gis.h>
 
79
#include <grass/glocale.h>
 
80
}
 
81
 
 
82
 
 
83
 
 
84
TreeNode *NIL;
 
85
 
 
86
#define EPSILON 0.0000001
 
87
 
 
88
 
 
89
/*public:--------------------------------- */
 
90
RBTree *create_tree(TreeValue tv)
 
91
{
 
92
    init_nil_node();
 
93
    RBTree *rbt = (RBTree *) G_malloc(sizeof(RBTree));
 
94
    TreeNode *root = (TreeNode *) G_malloc(sizeof(TreeNode));
 
95
 
 
96
    rbt->root = root;
 
97
    rbt->root->value = tv;
 
98
    rbt->root->left = NIL;
 
99
    rbt->root->right = NIL;
 
100
    rbt->root->parent = NIL;
 
101
    rbt->root->color = RB_BLACK;
 
102
 
 
103
    return rbt;
 
104
}
 
105
 
 
106
/*LT: not sure if this is correct */
 
107
int is_empty(RBTree * t)
 
108
{
 
109
    assert(t);
 
110
    return (t->root == NIL);
 
111
}
 
112
 
 
113
void delete_tree(RBTree * t)
 
114
{
 
115
    destroy_sub_tree(t->root);
 
116
    return;
 
117
}
 
118
 
 
119
void destroy_sub_tree(TreeNode * node)
 
120
{
 
121
    if (node == NIL)
 
122
        return;
 
123
 
 
124
    destroy_sub_tree(node->left);
 
125
    destroy_sub_tree(node->right);
 
126
    G_free(node);
 
127
    return;
 
128
}
 
129
 
 
130
void insert_into(RBTree * rbt, TreeValue value)
 
131
{
 
132
    insert_into_tree(&(rbt->root), value);
 
133
    return;
 
134
}
 
135
 
 
136
void delete_from(RBTree * rbt, double key)
 
137
{
 
138
    delete_from_tree(&(rbt->root), key);
 
139
    return;
 
140
}
 
141
 
 
142
TreeNode *search_for_node_with_key(RBTree * rbt, double key)
 
143
{
 
144
    return search_for_node(rbt->root, key);
 
145
}
 
146
 
 
147
/*------------The following is designed for viewshed's algorithm-------*/
 
148
double find_max_gradient_within_key(RBTree * rbt, double key, double angle, double gradient)
 
149
{
 
150
    return find_max_value_within_key(rbt->root, key, angle, gradient);
 
151
}
 
152
 
 
153
/*<--------------------------------->
 
154
   //Private below this line */
 
155
void init_nil_node()
 
156
{
 
157
    NIL = (TreeNode *) G_malloc(sizeof(TreeNode));
 
158
    NIL->color = RB_BLACK;
 
159
    NIL->value.angle[0] = 0;
 
160
    NIL->value.angle[1] = 0;
 
161
    NIL->value.angle[2] = 0;
 
162
    NIL->value.gradient[0] = SMALLEST_GRADIENT;
 
163
    NIL->value.gradient[1] = SMALLEST_GRADIENT;
 
164
    NIL->value.gradient[2] = SMALLEST_GRADIENT;
 
165
    NIL->value.maxGradient = SMALLEST_GRADIENT;
 
166
    NIL->value.key = 0;
 
167
 
 
168
    NIL->parent = NULL;
 
169
    NIL->left = NULL;
 
170
    NIL->right = NULL;
 
171
    return;
 
172
}
 
173
 
 
174
/*you can write change this compare function, depending on your TreeValue struct
 
175
   //compare function used by findMaxValue
 
176
   //-1: v1 < v2
 
177
   //0:  v1 = v2
 
178
   //2:  v1 > v2 */
 
179
char compare_values(TreeValue * v1, TreeValue * v2)
 
180
{
 
181
    if (v1->gradient[1] > v2->gradient[1])
 
182
        return 1;
 
183
    if (v1->gradient[1] < v2->gradient[1])
 
184
        return -1;
 
185
 
 
186
    return 0;
 
187
}
 
188
 
 
189
 
 
190
/*a function used to compare two doubles */
 
191
/* a < b : -1
 
192
 * a > b : 1
 
193
 * a == b : 0 */
 
194
char compare_double(double a, double b)
 
195
{
 
196
    return (a < b ? -1 : (a > b));
 
197
}
 
198
 
 
199
 
 
200
 
 
201
/*create a tree node */
 
202
TreeNode *create_tree_node(TreeValue value)
 
203
{
 
204
    TreeNode *ret;
 
205
 
 
206
    ret = (TreeNode *) G_malloc(sizeof(TreeNode));
 
207
 
 
208
    ret->color = RB_RED;
 
209
 
 
210
    ret->left = NIL;
 
211
    ret->right = NIL;
 
212
    ret->parent = NIL;
 
213
 
 
214
    ret->value = value;
 
215
    ret->value.maxGradient = SMALLEST_GRADIENT;
 
216
    return ret;
 
217
}
 
218
 
 
219
/*create node with its value set to the value given
 
220
   //and insert the node into the tree
 
221
   //rbInsertFixup may change the root pointer, so TreeNode** is passed in */
 
222
void insert_into_tree(TreeNode ** root, TreeValue value)
 
223
{
 
224
    TreeNode *curNode;
 
225
    TreeNode *nextNode;
 
226
 
 
227
    curNode = *root;
 
228
 
 
229
    if (compare_double(value.key, curNode->value.key) == -1) {
 
230
        nextNode = curNode->left;
 
231
    }
 
232
    else {
 
233
        nextNode = curNode->right;
 
234
    }
 
235
 
 
236
 
 
237
    while (nextNode != NIL) {
 
238
        curNode = nextNode;
 
239
 
 
240
        if (compare_double(value.key, curNode->value.key) == -1) {
 
241
            nextNode = curNode->left;
 
242
        }
 
243
        else {
 
244
            nextNode = curNode->right;
 
245
        }
 
246
    }
 
247
 
 
248
    /*create a new node 
 
249
       //and place it at the right place
 
250
       //created node is RED by default */
 
251
    nextNode = create_tree_node(value);
 
252
 
 
253
    nextNode->parent = curNode;
 
254
 
 
255
    if (compare_double(value.key, curNode->value.key) == -1) {
 
256
        curNode->left = nextNode;
 
257
    }
 
258
    else {
 
259
        curNode->right = nextNode;
 
260
    }
 
261
 
 
262
    TreeNode *inserted = nextNode;
 
263
 
 
264
    /*update augmented maxGradient */
 
265
    nextNode->value.maxGradient = find_value_min_value(&nextNode->value);
 
266
    while (nextNode->parent != NIL) {
 
267
        if (nextNode->parent->value.maxGradient < nextNode->value.maxGradient)
 
268
            nextNode->parent->value.maxGradient = nextNode->value.maxGradient;
 
269
 
 
270
        if (nextNode->parent->value.maxGradient > nextNode->value.maxGradient)
 
271
            break;
 
272
        nextNode = nextNode->parent;
 
273
    }
 
274
 
 
275
    /*fix rb tree after insertion */
 
276
    rb_insert_fixup(root, inserted);
 
277
 
 
278
    return;
 
279
}
 
280
 
 
281
void rb_insert_fixup(TreeNode ** root, TreeNode * z)
 
282
{
 
283
    /*see pseudocode on page 281 in CLRS */
 
284
    TreeNode *y;
 
285
 
 
286
    while (z->parent->color == RB_RED) {
 
287
        if (z->parent == z->parent->parent->left) {
 
288
            y = z->parent->parent->right;
 
289
            if (y->color == RB_RED) {   /*case 1 */
 
290
                z->parent->color = RB_BLACK;
 
291
                y->color = RB_BLACK;
 
292
                z->parent->parent->color = RB_RED;
 
293
                z = z->parent->parent;
 
294
            }
 
295
            else {
 
296
                if (z == z->parent->right) {    /*case 2 */
 
297
                    z = z->parent;
 
298
                    left_rotate(root, z);       /*convert case 2 to case 3 */
 
299
                }
 
300
                z->parent->color = RB_BLACK;    /*case 3 */
 
301
                z->parent->parent->color = RB_RED;
 
302
                right_rotate(root, z->parent->parent);
 
303
            }
 
304
 
 
305
        }
 
306
        else {                  /*(z->parent == z->parent->parent->right) */
 
307
            y = z->parent->parent->left;
 
308
            if (y->color == RB_RED) {   /*case 1 */
 
309
                z->parent->color = RB_BLACK;
 
310
                y->color = RB_BLACK;
 
311
                z->parent->parent->color = RB_RED;
 
312
                z = z->parent->parent;
 
313
            }
 
314
            else {
 
315
                if (z == z->parent->left) {     /*case 2 */
 
316
                    z = z->parent;
 
317
                    right_rotate(root, z);      /*convert case 2 to case 3 */
 
318
                }
 
319
                z->parent->color = RB_BLACK;    /*case 3 */
 
320
                z->parent->parent->color = RB_RED;
 
321
                left_rotate(root, z->parent->parent);
 
322
            }
 
323
        }
 
324
    }
 
325
    (*root)->color = RB_BLACK;
 
326
 
 
327
    return;
 
328
}
 
329
 
 
330
 
 
331
 
 
332
 
 
333
/*search for a node with the given key */
 
334
TreeNode *search_for_node(TreeNode * root, double key)
 
335
{
 
336
    TreeNode *curNode = root;
 
337
 
 
338
    while (curNode != NIL && compare_double(key, curNode->value.key) != 0) {
 
339
 
 
340
        if (compare_double(key, curNode->value.key) == -1) {
 
341
            curNode = curNode->left;
 
342
        }
 
343
        else {
 
344
            curNode = curNode->right;
 
345
        }
 
346
 
 
347
    }
 
348
 
 
349
    return curNode;
 
350
}
 
351
 
 
352
/*function used by treeSuccessor */
 
353
TreeNode *tree_minimum(TreeNode * x)
 
354
{
 
355
    while (x->left != NIL)
 
356
        x = x->left;
 
357
 
 
358
    return x;
 
359
}
 
360
 
 
361
/*function used by deletion */
 
362
TreeNode *tree_successor(TreeNode * x)
 
363
{
 
364
    if (x->right != NIL)
 
365
        return tree_minimum(x->right);
 
366
    TreeNode *y = x->parent;
 
367
 
 
368
    while (y != NIL && x == y->right) {
 
369
        x = y;
 
370
        if (y->parent == NIL)
 
371
            return y;
 
372
        y = y->parent;
 
373
    }
 
374
    return y;
 
375
}
 
376
 
 
377
 
 
378
/*delete the node out of the tree */
 
379
void delete_from_tree(TreeNode ** root, double key)
 
380
{
 
381
    double tmpMax;
 
382
    TreeNode *z;
 
383
    TreeNode *x;
 
384
    TreeNode *y;
 
385
    TreeNode *toFix;
 
386
 
 
387
    z = search_for_node(*root, key);
 
388
 
 
389
    if (z == NIL) {
 
390
        /*node to delete is not found */
 
391
        G_fatal_error(_("Attempt to delete node with key=%f failed"), key);
 
392
    }
 
393
 
 
394
    /*1-3 */
 
395
    if (z->left == NIL || z->right == NIL)
 
396
        y = z;
 
397
    else
 
398
        y = tree_successor(z);
 
399
        
 
400
    if (y == NIL) {
 
401
        G_fatal_error(_("Successor node not found. Deletion fails."));
 
402
    }
 
403
 
 
404
    /*4-6 */
 
405
    if (y->left != NIL)
 
406
        x = y->left;
 
407
    else
 
408
        x = y->right;
 
409
 
 
410
    /*7 */
 
411
    x->parent = y->parent;
 
412
 
 
413
    /*8-12 */
 
414
    if (y->parent == NIL) {
 
415
        *root = x;
 
416
 
 
417
        toFix = *root;          /*augmentation to be fixed */
 
418
    }
 
419
    else {
 
420
        if (y == y->parent->left)
 
421
            y->parent->left = x;
 
422
        else
 
423
            y->parent->right = x;
 
424
 
 
425
        toFix = y->parent;      /*augmentation to be fixed */
 
426
    }
 
427
 
 
428
    /*fix augmentation for removing y */
 
429
    TreeNode *curNode = y;
 
430
    double left, right;
 
431
 
 
432
    while (curNode->parent != NIL) {
 
433
        if (curNode->parent->value.maxGradient == find_value_min_value(&y->value)) {
 
434
            left = find_max_value(curNode->parent->left);
 
435
            right = find_max_value(curNode->parent->right);
 
436
 
 
437
            if (left > right)
 
438
                curNode->parent->value.maxGradient = left;
 
439
            else
 
440
                curNode->parent->value.maxGradient = right;
 
441
 
 
442
            if (find_value_min_value(&curNode->parent->value) >
 
443
                curNode->parent->value.maxGradient)
 
444
                curNode->parent->value.maxGradient =
 
445
                    find_value_min_value(&curNode->parent->value);
 
446
        }
 
447
        else {
 
448
            break;
 
449
        }
 
450
        curNode = curNode->parent;
 
451
    }
 
452
 
 
453
 
 
454
    /*fix augmentation for x */
 
455
    tmpMax =
 
456
        toFix->left->value.maxGradient >
 
457
        toFix->right->value.maxGradient ? toFix->left->value.
 
458
        maxGradient : toFix->right->value.maxGradient;
 
459
    if (tmpMax > find_value_min_value(&toFix->value))
 
460
        toFix->value.maxGradient = tmpMax;
 
461
    else
 
462
        toFix->value.maxGradient = find_value_min_value(&toFix->value);
 
463
 
 
464
    /*13-15 */
 
465
    if (y != NIL && y != z) {
 
466
        double zGradient = find_value_min_value(&z->value);
 
467
 
 
468
        z->value.key = y->value.key;
 
469
        z->value.gradient[0] = y->value.gradient[0];
 
470
        z->value.gradient[1] = y->value.gradient[1];
 
471
        z->value.gradient[2] = y->value.gradient[2];
 
472
        z->value.angle[0] = y->value.angle[0];
 
473
        z->value.angle[1] = y->value.angle[1];
 
474
        z->value.angle[2] = y->value.angle[2];
 
475
 
 
476
 
 
477
        toFix = z;
 
478
        /*fix augmentation */
 
479
        tmpMax =
 
480
            toFix->left->value.maxGradient >
 
481
            toFix->right->value.maxGradient ? toFix->left->value.
 
482
            maxGradient : toFix->right->value.maxGradient;
 
483
        if (tmpMax > find_value_min_value(&toFix->value))
 
484
            toFix->value.maxGradient = tmpMax;
 
485
        else
 
486
            toFix->value.maxGradient = find_value_min_value(&toFix->value);
 
487
 
 
488
        while (z->parent != NIL) {
 
489
            if (z->parent->value.maxGradient == zGradient) {
 
490
                if (find_value_min_value(&z->parent->value) != zGradient &&
 
491
                    (!(z->parent->left->value.maxGradient == zGradient &&
 
492
                       z->parent->right->value.maxGradient == zGradient))) {
 
493
 
 
494
                    left = find_max_value(z->parent->left);
 
495
                    right = find_max_value(z->parent->right);
 
496
 
 
497
                    if (left > right)
 
498
                        z->parent->value.maxGradient = left;
 
499
                    else
 
500
                        z->parent->value.maxGradient = right;
 
501
 
 
502
                    if (find_value_min_value(&z->parent->value) >
 
503
                        z->parent->value.maxGradient)
 
504
                        z->parent->value.maxGradient =
 
505
                            find_value_min_value(&z->parent->value);
 
506
 
 
507
                }
 
508
 
 
509
            }
 
510
            else {
 
511
                if (z->value.maxGradient > z->parent->value.maxGradient)
 
512
                    z->parent->value.maxGradient = z->value.maxGradient;
 
513
            }
 
514
            z = z->parent;
 
515
        }
 
516
 
 
517
    }
 
518
 
 
519
    /*16-17 */
 
520
    if (y->color == RB_BLACK && x != NIL)
 
521
        rb_delete_fixup(root, x);
 
522
 
 
523
    /*18 */
 
524
    G_free(y);
 
525
 
 
526
    return;
 
527
}
 
528
 
 
529
/*fix the rb tree after deletion */
 
530
void rb_delete_fixup(TreeNode ** root, TreeNode * x)
 
531
{
 
532
    TreeNode *w;
 
533
 
 
534
    while (x != *root && x->color == RB_BLACK) {
 
535
        if (x == x->parent->left) {
 
536
            w = x->parent->right;
 
537
            if (w->color == RB_RED) {
 
538
                w->color = RB_BLACK;
 
539
                x->parent->color = RB_RED;
 
540
                left_rotate(root, x->parent);
 
541
                w = x->parent->right;
 
542
            }
 
543
 
 
544
            if (w == NIL) {
 
545
                x = x->parent;
 
546
                continue;
 
547
            }
 
548
 
 
549
            if (w->left->color == RB_BLACK && w->right->color == RB_BLACK) {
 
550
                w->color = RB_RED;
 
551
                x = x->parent;
 
552
            }
 
553
            else {
 
554
                if (w->right->color == RB_BLACK) {
 
555
                    w->left->color = RB_BLACK;
 
556
                    w->color = RB_RED;
 
557
                    right_rotate(root, w);
 
558
                    w = x->parent->right;
 
559
                }
 
560
 
 
561
                w->color = x->parent->color;
 
562
                x->parent->color = RB_BLACK;
 
563
                w->right->color = RB_BLACK;
 
564
                left_rotate(root, x->parent);
 
565
                x = *root;
 
566
            }
 
567
 
 
568
        }
 
569
        else {                  /*(x==x->parent->right) */
 
570
            w = x->parent->left;
 
571
            if (w->color == RB_RED) {
 
572
                w->color = RB_BLACK;
 
573
                x->parent->color = RB_RED;
 
574
                right_rotate(root, x->parent);
 
575
                w = x->parent->left;
 
576
            }
 
577
 
 
578
            if (w == NIL) {
 
579
                x = x->parent;
 
580
                continue;
 
581
            }
 
582
 
 
583
            if (w->right->color == RB_BLACK && w->left->color == RB_BLACK) {
 
584
                w->color = RB_RED;
 
585
                x = x->parent;
 
586
            }
 
587
            else {
 
588
                if (w->left->color == RB_BLACK) {
 
589
                    w->right->color = RB_BLACK;
 
590
                    w->color = RB_RED;
 
591
                    left_rotate(root, w);
 
592
                    w = x->parent->left;
 
593
                }
 
594
 
 
595
                w->color = x->parent->color;
 
596
                x->parent->color = RB_BLACK;
 
597
                w->left->color = RB_BLACK;
 
598
                right_rotate(root, x->parent);
 
599
                x = *root;
 
600
            }
 
601
 
 
602
        }
 
603
    }
 
604
    x->color = RB_BLACK;
 
605
 
 
606
    return;
 
607
}
 
608
 
 
609
/*find the min value in the given tree value */
 
610
double find_value_min_value(TreeValue * v)
 
611
{
 
612
    if (v->gradient[0] < v->gradient[1]) {
 
613
        if (v->gradient[0] < v->gradient[2])
 
614
            return v->gradient[0];
 
615
        else
 
616
            return v->gradient[2];
 
617
    }
 
618
    else {
 
619
        if (v->gradient[1] < v->gradient[2])
 
620
            return v->gradient[1];
 
621
        else
 
622
            return v->gradient[2];
 
623
    }
 
624
    return v->gradient[0];
 
625
}
 
626
 
 
627
/*find the max value in the given tree
 
628
   //you need to provide a compare function to compare the nodes */
 
629
double find_max_value(TreeNode * root)
 
630
{
 
631
    if (!root)
 
632
        return SMALLEST_GRADIENT;
 
633
    assert(root);
 
634
    /*assert(root->value.maxGradient != SMALLEST_GRADIENT);
 
635
       //LT: this shoudl be fixed
 
636
       //if (root->value.maxGradient != SMALLEST_GRADIENT) */
 
637
    return root->value.maxGradient;
 
638
}
 
639
 
 
640
 
 
641
 
 
642
/* find max within the max key */
 
643
double find_max_value_within_key(TreeNode * root, double maxKey, double angle, double gradient)
 
644
{
 
645
    TreeNode *keyNode = search_for_node(root, maxKey);
 
646
 
 
647
    if (keyNode == NIL) {
 
648
        /*fprintf(stderr, "key node not found. error occurred!\n");
 
649
           //there is no point in the structure with key < maxKey */
 
650
        /*node not found */
 
651
        G_fatal_error(_("Attempt to find node with key=%f failed"), maxKey);
 
652
        return SMALLEST_GRADIENT;
 
653
    }
 
654
 
 
655
    TreeNode *currNode = keyNode;
 
656
    double max = SMALLEST_GRADIENT;
 
657
    double tmpMax;
 
658
    double curr_gradient;
 
659
 
 
660
    while (currNode->parent != NIL) {
 
661
        if (currNode == currNode->parent->right) {      /*its the right node of its parent; */
 
662
            tmpMax = find_max_value(currNode->parent->left);
 
663
            if (tmpMax > max)
 
664
                max = tmpMax;
 
665
            if (find_value_min_value(&currNode->parent->value) > max)
 
666
                max = find_value_min_value(&currNode->parent->value);
 
667
        }
 
668
        currNode = currNode->parent;
 
669
    }
 
670
 
 
671
    if (max > gradient)
 
672
        return max;
 
673
 
 
674
    /* traverse all nodes with smaller distance */
 
675
    max = SMALLEST_GRADIENT;
 
676
    currNode = keyNode;
 
677
    while (currNode != NIL) {
 
678
        int checkme = (currNode->value.angle[0] <= angle &&
 
679
                      currNode->value.angle[2] >= angle);
 
680
                      
 
681
        if (!checkme && currNode->value.key > 0) {
 
682
            G_warning(_("Angles outside angle %.4f"), angle);
 
683
            G_warning(_("ENTER angle %.4f"), currNode->value.angle[0]);
 
684
            G_warning(_("CENTER angle %.4f"), currNode->value.angle[1]);
 
685
            G_warning(_("EXIT angle %.4f"), currNode->value.angle[2]);
 
686
            G_warning(_("ENTER gradient %.4f"), currNode->value.gradient[0]);
 
687
            G_warning(_("CENTER gradient %.4f"), currNode->value.gradient[1]);
 
688
            G_warning(_("EXIT gradient %.4f"), currNode->value.gradient[2]);
 
689
        }
 
690
        
 
691
        if (currNode->value.key > maxKey) {
 
692
            G_fatal_error(_("current dist too large %.4f > %.4f"), currNode->value.key, maxKey);
 
693
        }
 
694
            
 
695
            
 
696
        if (checkme && currNode != keyNode) {
 
697
            if (angle < currNode->value.angle[1]) {
 
698
                curr_gradient = currNode->value.gradient[1] +
 
699
                  (currNode->value.gradient[0] - currNode->value.gradient[1]) *
 
700
                  (currNode->value.angle[1] - angle) /
 
701
                  (currNode->value.angle[1] - currNode->value.angle[0]);
 
702
            }
 
703
            else if (angle > currNode->value.angle[1]) {
 
704
                curr_gradient = currNode->value.gradient[1] +
 
705
                  (currNode->value.gradient[2] - currNode->value.gradient[1]) *
 
706
                  (angle - currNode->value.angle[1]) /
 
707
                  (currNode->value.angle[2] - currNode->value.angle[1]);
 
708
            }
 
709
            else /* angle == currNode->value.angle[1] */
 
710
                curr_gradient = currNode->value.gradient[1];
 
711
 
 
712
            if (curr_gradient > max)
 
713
                max = curr_gradient;
 
714
                
 
715
            if (max > gradient)
 
716
                return max;
 
717
        }
 
718
        /* get next smaller key */
 
719
        if (currNode->left != NIL) {
 
720
            currNode = currNode->left;
 
721
            while (currNode->right != NIL)
 
722
                currNode = currNode->right;
 
723
        }
 
724
        else {
 
725
            /* at smallest item in this branch, go back up */
 
726
            TreeNode *lastNode;
 
727
            
 
728
            do {
 
729
                lastNode = currNode;
 
730
                currNode = currNode->parent;
 
731
            } while (currNode != NIL && lastNode == currNode->left);
 
732
        }
 
733
    }
 
734
    return max;
 
735
 
 
736
    /* old code assuming flat cells */
 
737
    while (keyNode->parent != NIL) {
 
738
        if (keyNode == keyNode->parent->right) {        /*its the right node of its parent; */
 
739
            tmpMax = find_max_value(keyNode->parent->left);
 
740
            if (tmpMax > max)
 
741
                max = tmpMax;
 
742
            if (keyNode->parent->value.gradient[1] > max)
 
743
                max = keyNode->parent->value.gradient[1];
 
744
        }
 
745
        keyNode = keyNode->parent;
 
746
    }
 
747
 
 
748
    return max;
 
749
}
 
750
 
 
751
 
 
752
void left_rotate(TreeNode ** root, TreeNode * x)
 
753
{
 
754
    TreeNode *y;
 
755
 
 
756
    y = x->right;
 
757
 
 
758
    /*maintain augmentation */
 
759
    double tmpMax;
 
760
 
 
761
    /*fix x */
 
762
    tmpMax = x->left->value.maxGradient > y->left->value.maxGradient ?
 
763
        x->left->value.maxGradient : y->left->value.maxGradient;
 
764
 
 
765
    if (tmpMax > find_value_min_value(&x->value))
 
766
        x->value.maxGradient = tmpMax;
 
767
    else
 
768
        x->value.maxGradient = find_value_min_value(&x->value);
 
769
 
 
770
 
 
771
    /*fix y */
 
772
    tmpMax = x->value.maxGradient > y->right->value.maxGradient ?
 
773
        x->value.maxGradient : y->right->value.maxGradient;
 
774
 
 
775
    if (tmpMax > find_value_min_value(&y->value))
 
776
        y->value.maxGradient = tmpMax;
 
777
    else
 
778
        y->value.maxGradient = find_value_min_value(&y->value);
 
779
 
 
780
    /*left rotation
 
781
       //see pseudocode on page 278 in CLRS */
 
782
 
 
783
    x->right = y->left;         /*turn y's left subtree into x's right subtree */
 
784
    y->left->parent = x;
 
785
 
 
786
    y->parent = x->parent;      /*link x's parent to y */
 
787
 
 
788
    if (x->parent == NIL) {
 
789
        *root = y;
 
790
    }
 
791
    else {
 
792
        if (x == x->parent->left)
 
793
            x->parent->left = y;
 
794
        else
 
795
            x->parent->right = y;
 
796
    }
 
797
 
 
798
    y->left = x;
 
799
    x->parent = y;
 
800
 
 
801
    return;
 
802
}
 
803
 
 
804
void right_rotate(TreeNode ** root, TreeNode * y)
 
805
{
 
806
    TreeNode *x;
 
807
 
 
808
    x = y->left;
 
809
 
 
810
    /*maintain augmentation
 
811
       //fix y */
 
812
    double tmpMax;
 
813
 
 
814
    tmpMax = x->right->value.maxGradient > y->right->value.maxGradient ?
 
815
        x->right->value.maxGradient : y->right->value.maxGradient;
 
816
 
 
817
    if (tmpMax > find_value_min_value(&y->value))
 
818
        y->value.maxGradient = tmpMax;
 
819
    else
 
820
        y->value.maxGradient = find_value_min_value(&y->value);
 
821
 
 
822
    /*fix x */
 
823
    tmpMax = x->left->value.maxGradient > y->value.maxGradient ?
 
824
        x->left->value.maxGradient : y->value.maxGradient;
 
825
 
 
826
    if (tmpMax > find_value_min_value(&x->value))
 
827
        x->value.maxGradient = tmpMax;
 
828
    else
 
829
        x->value.maxGradient = find_value_min_value(&x->value);
 
830
 
 
831
    /*ratation */
 
832
    y->left = x->right;
 
833
    x->right->parent = y;
 
834
 
 
835
    x->parent = y->parent;
 
836
 
 
837
    if (y->parent == NIL) {
 
838
        *root = x;
 
839
    }
 
840
    else {
 
841
        if (y->parent->left == y)
 
842
            y->parent->left = x;
 
843
        else
 
844
            y->parent->right = x;
 
845
    }
 
846
 
 
847
    x->right = y;
 
848
    y->parent = x;
 
849
 
 
850
    return;
 
851
}