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:--------------------------------- */
create_tree(TreeValue tv)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 */
is_empty(RBTree * t)107 int is_empty(RBTree * t)
108 {
109     assert(t);
110     return (t->root == NIL);
111 }
112 
delete_tree(RBTree * t)113 void delete_tree(RBTree * t)
114 {
115     destroy_sub_tree(t->root);
116     return;
117 }
118 
destroy_sub_tree(TreeNode * node)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 
insert_into(RBTree * rbt,TreeValue value)130 void insert_into(RBTree * rbt, TreeValue value)
131 {
132     insert_into_tree(&(rbt->root), value);
133     return;
134 }
135 
delete_from(RBTree * rbt,double key)136 void delete_from(RBTree * rbt, double key)
137 {
138     delete_from_tree(&(rbt->root), key);
139     return;
140 }
141 
search_for_node_with_key(RBTree * rbt,double key)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-------*/
find_max_gradient_within_key(RBTree * rbt,double key,double angle,double gradient)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 */
init_nil_node()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 */
compare_values(TreeValue * v1,TreeValue * 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 */
compare_double(double a,double b)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 */
create_tree_node(TreeValue value)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 */
insert_into_tree(TreeNode ** root,TreeValue value)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 
rb_insert_fixup(TreeNode ** root,TreeNode * z)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 */
search_for_node(TreeNode * root,double 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 */
tree_minimum(TreeNode * x)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 */
tree_successor(TreeNode * x)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 */
delete_from_tree(TreeNode ** root,double key)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 */
rb_delete_fixup(TreeNode ** root,TreeNode * x)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 */
find_value_min_value(TreeValue * v)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 */
find_max_value(TreeNode * root)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 should 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 */
find_max_value_within_key(TreeNode * root,double maxKey,double angle,double gradient)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 
left_rotate(TreeNode ** root,TreeNode * x)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 
right_rotate(TreeNode ** root,TreeNode * y)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 }
852