00001
00102 #ifndef __PUSHRELABEL_NON_LINEAR_GRID_3D__
00103 #define __PUSHRELABEL_NON_LINEAR_GRID_3D__
00104
00105
00106
00107 #include <algorithm>
00108 #include <deque>
00109 #include <iterator>
00110
00111 #ifndef WITHOUT_LIMITS
00112 #include <limits>
00113 #endif
00114
00115 #include <list>
00116 #include <vector>
00117
00118 #include "array.h"
00119 #include "geom.h"
00120 #include "msg_stream.h"
00121
00122
00123
00124 namespace Graph_cut{
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00140
00144 template<typename Real>
00145 class PushRelabel_non_linear_grid_3D{
00146
00147 public:
00148
00150
00151
00152 typedef Real real_type;
00153 typedef Array_3D<bool>::size_type size_type;
00155
00157 PushRelabel_non_linear_grid_3D(
00158 const Array_3D<real_type>& cost_edge_cap,
00159 const Array_3D<bool>& dom,
00160 const Array_2D<real_type>& x_derivative_edge_cap,
00161 const Array_2D<real_type>& y_derivative_edge_cap,
00162 const real_type eps = 1e-3);
00163
00165 real_type max_flow();
00166
00168 void min_cut(Array_2D<size_type>* const cut) const;
00169
00170
00171
00172 bool DBG_LOGGING;
00173
00175
00176
00177 typedef Geometry::Vec3<unsigned int> position_type;
00178 typedef unsigned int index_type;
00179 typedef unsigned int label_type;
00181
00183
00184 typedef char direction_type;
00185
00186 static const direction_type X_MIN_DIR = 1;
00187 static const direction_type X_MAX_DIR = 2;
00188
00189 static const direction_type Y_MIN_DIR = 3;
00190 static const direction_type Y_MAX_DIR = 4;
00191
00192 static const direction_type Z_MIN_x_min_DIR = 5;
00193 static const direction_type Z_MIN_x_max_DIR = 6;
00194 static const direction_type Z_MIN_y_min_DIR = 7;
00195 static const direction_type Z_MIN_y_max_DIR = 8;
00196
00197 static const direction_type Z_MAX_x_min_DIR = 9;
00198 static const direction_type Z_MAX_x_max_DIR = 10;
00199 static const direction_type Z_MAX_y_min_DIR = 11;
00200 static const direction_type Z_MAX_y_max_DIR = 12;
00202
00204
00205 typedef char location_type;
00206
00207 static const location_type CENTER_NODE = 100;
00208 static const location_type Z_MAX_x_min_NODE = 101;
00209 static const location_type Z_MAX_x_max_NODE = 102;
00210 static const location_type Z_MAX_y_min_NODE = 103;
00211 static const location_type Z_MAX_y_max_NODE = 104;
00213
00215 struct node_type{
00216 position_type pos;
00217 location_type loc;
00218
00219 inline node_type(){}
00220 inline node_type(const position_type& p,const location_type l):pos(p),loc(l){}
00221 inline bool operator==(const node_type& n) const {
00222 return (pos==n.pos)&&(loc==n.loc);
00223 }
00224 inline bool operator!=(const node_type& n) const {
00225 return !((pos==n.pos)&&(loc==n.loc));
00226 }
00227 };
00228
00230 const real_type epsilon;
00231
00233 size_type nb_nodes;
00234
00236 size_type nb_relabel;
00237
00239 Array_3D<real_type> cost_edge_capacity;
00240
00243 const Array_3D<bool>& domain;
00244
00246
00247
00248 Array_2D<real_type> x_major_derivative_edge_capacity;
00249 Array_2D<real_type> x_minor_derivative_edge_capacity;
00250
00251 Array_2D<real_type> y_major_derivative_edge_capacity;
00252 Array_2D<real_type> y_minor_derivative_edge_capacity;
00254
00256
00257
00258 Array_3D<index_type> label_CENTER;
00259 Array_3D<index_type> label_Z_MAX_x_min;
00260 Array_3D<index_type> label_Z_MAX_x_max;
00261 Array_3D<index_type> label_Z_MAX_y_min;
00262 Array_3D<index_type> label_Z_MAX_y_max;
00264
00266
00267
00269 Array_3D<real_type> cost_edge_Z_MAX_x_min_flow;
00270 Array_3D<real_type> cost_edge_Z_MAX_x_max_flow;
00271 Array_3D<real_type> cost_edge_Z_MAX_y_min_flow;
00272 Array_3D<real_type> cost_edge_Z_MAX_y_max_flow;
00273
00274 Array_3D<real_type> cost_edge_Z_MAX_MAX_x_min_flow;
00275 Array_3D<real_type> cost_edge_Z_MAX_MAX_x_max_flow;
00276 Array_3D<real_type> cost_edge_Z_MAX_MAX_y_min_flow;
00277 Array_3D<real_type> cost_edge_Z_MAX_MAX_y_max_flow;
00279
00281
00282
00284 Array_3D<real_type> x_major_derivative_edge_flow;
00285 Array_3D<real_type> y_major_derivative_edge_flow;
00286
00287 Array_3D<real_type> x_minor_derivative_edge_flow;
00288 Array_3D<real_type> y_minor_derivative_edge_flow;
00290
00292 Array_2D<size_type> min_bound;
00293
00295 Array_2D<size_type> max_bound;
00296
00298 std::deque<index_type> active_node_index;
00299
00301 const node_type source;
00302
00304 const node_type target;
00305
00307
00308
00309 const size_type x_size,y_size,z_size;
00311
00312
00314 inline void to_index(const node_type& p,
00315 size_type* const i) const;
00316
00318 inline void from_index(node_type* const p,
00319 const size_type i) const;
00320
00322
00323
00324 inline const label_type& label(const node_type& node) const;
00325 inline label_type& label(const node_type& node);
00327
00329 inline void neighbors(const node_type& node,
00330 node_type* const result,
00331 size_type* const n_res) const;
00332
00334 inline void adjacent_nodes(const node_type& node,
00335 node_type* const result,
00336 size_type* const n_res) const;
00337
00339 inline real_type edge_capacity(const node_type& from,
00340 const direction_type dir) const;
00341
00343 inline real_type edge_flow(const node_type& from,
00344 const direction_type dir) const;
00345
00346
00348 inline real_type edge_residual(const node_type& from,
00349 const direction_type dir) const;
00350
00352 inline real_type node_excess(const node_type& node) const;
00353
00354
00356 inline void add_flow(const node_type& from,
00357 const direction_type dir,
00358 const real_type flow);
00359
00361 inline direction_type direction(const node_type& from,
00362 const node_type& to) const;
00363
00365 inline label_type relabel(const node_type& node);
00366
00368 void global_relabel();
00369
00371 inline bool discharge(const node_type& node);
00372
00374 inline void add_active_node(const node_type& node);
00375
00377 inline void get_next_node(node_type* const node,
00378 bool* const no_more_node);
00379
00381 inline bool is_active(const node_type& node) const;
00382
00384 inline bool is_admissible(const node_type& from,
00385 const node_type& to) const;
00386
00387
00388 };
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00433 template<typename Real>
00434 PushRelabel_non_linear_grid_3D<Real>::
00435 PushRelabel_non_linear_grid_3D(
00436 const Array_3D<real_type>& cost_edge_cap,
00437 const Array_3D<bool>& dom,
00438 const Array_2D<real_type>& x_derivative_edge_cap,
00439 const Array_2D<real_type>& y_derivative_edge_cap,
00440 const real_type eps)
00441 :epsilon(eps),
00442 nb_nodes(2),
00443 nb_relabel(0),
00444
00445 domain(dom),
00446 DBG_LOGGING(false),
00447
00448 x_major_derivative_edge_capacity(x_derivative_edge_cap),
00449 x_minor_derivative_edge_capacity(x_derivative_edge_cap),
00450
00451 y_major_derivative_edge_capacity(y_derivative_edge_cap),
00452 y_minor_derivative_edge_capacity(y_derivative_edge_cap),
00453
00454 cost_edge_Z_MAX_x_min_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00455 cost_edge_Z_MAX_x_max_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00456 cost_edge_Z_MAX_y_min_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00457 cost_edge_Z_MAX_y_max_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00458
00459 cost_edge_Z_MAX_MAX_x_min_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00460 cost_edge_Z_MAX_MAX_x_max_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00461 cost_edge_Z_MAX_MAX_y_min_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00462 cost_edge_Z_MAX_MAX_y_max_flow(dom.x_size(),dom.y_size(),dom.z_size(),0),
00463
00464 x_major_derivative_edge_flow(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00465 x_minor_derivative_edge_flow(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00466 y_major_derivative_edge_flow(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00467 y_minor_derivative_edge_flow(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00468
00469 label_CENTER(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00470 label_Z_MAX_x_min(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00471 label_Z_MAX_x_max(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00472 label_Z_MAX_y_min(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00473 label_Z_MAX_y_max(dom.x_size(),dom.y_size(),dom.z_size()+1,0),
00474
00475 min_bound(dom.x_size(),dom.y_size()),
00476 max_bound(dom.x_size(),dom.y_size()),
00477
00478 x_size(dom.x_size()),
00479 y_size(dom.y_size()),
00480 z_size(dom.z_size()),
00481
00482 #ifndef WITHOUT_LIMITS
00483 source(position_type(std::numeric_limits<size_type>::max(),
00484 std::numeric_limits<size_type>::max(),
00485 std::numeric_limits<size_type>::max()),CENTER_NODE),
00486 target(position_type(std::numeric_limits<size_type>::max()-1,
00487 std::numeric_limits<size_type>::max()-1,
00488 std::numeric_limits<size_type>::max()-1),CENTER_NODE)
00489 #else
00490 source(position_type(static_cast<size_type>(4294967295),
00491 static_cast<size_type>(4294967295),
00492 static_cast<size_type>(4294967295)),CENTER_NODE),
00493 target(position_type(static_cast<size_type>(4294967294),
00494 static_cast<size_type>(4294967294),
00495 static_cast<size_type>(4294967294)),CENTER_NODE)
00496 #endif
00497 {
00498
00499
00500 for(size_type x=0;x<x_size;x++){
00501 for(size_type y=0;y<y_size;y++){
00502
00503 size_type top = 0;
00504 bool found = false;
00505
00506 for(top=0;(top<z_size)&&(!found);top++){
00507 found = domain(x,y,top);
00508 }
00509 min_bound(x,y) = top-1;
00510
00511 add_active_node(node_type(position_type(x,y,min_bound(x,y)),CENTER_NODE));
00512
00513 size_type bottom = 0;
00514 found = false;
00515 for(bottom=0;(bottom<z_size)&&(!found);bottom++){
00516 found = domain(x,y,z_size-1-bottom);
00517 }
00518 max_bound(x,y) = z_size - bottom + 1;
00519
00520
00521
00522 }
00523 }
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538 for(size_type x=0;x<x_size;x++){
00539 for(size_type y=0;y<y_size;y++){
00540
00541 x_major_derivative_edge_capacity(x,y) *= 1.0;
00542 x_minor_derivative_edge_capacity(x,y) *= 1.0;
00543
00544 y_major_derivative_edge_capacity(x,y) *= 1.0;
00545 y_minor_derivative_edge_capacity(x,y) *= 1.0;
00546
00547
00548
00549
00550
00551
00552
00553
00554 }
00555 }
00556
00557 real_type max_minor = 0;
00558 for(size_type x=0;x<x_size;x++){
00559 for(size_type y=0;y<y_size;y++){
00560 if (max_minor < x_minor_derivative_edge_capacity(x,y)){
00561 max_minor = x_minor_derivative_edge_capacity(x,y);
00562 }
00563 if (max_minor < y_minor_derivative_edge_capacity(x,y)){
00564 max_minor = y_minor_derivative_edge_capacity(x,y);
00565 }
00566 }
00567 }
00568
00569 cost_edge_capacity.resize(x_size,y_size,z_size);
00570
00571 for(size_type x=0;x<x_size;x++){
00572 for(size_type y=0;y<y_size;y++){
00573 for(size_type z=0;z<z_size;z++){
00574 cost_edge_capacity(x,y,z) = cost_edge_cap(x,y,z)/4 + max_minor;
00575
00576 if (domain(x,y,z)){
00577 nb_nodes += 12;
00578 }
00579 }
00580 }
00581 }
00582
00583
00584 }
00585
00586
00587 template<typename Real>
00588 void
00589 PushRelabel_non_linear_grid_3D<Real>::
00590 to_index(const node_type& p,
00591 index_type* const i) const{
00592
00593 if (p==source) {
00594 #ifndef WITHOUT_LIMITS
00595 *i = std::numeric_limits<index_type>::max();
00596 #else
00597 *i = static_cast<size_type>(4294967295);
00598 #endif
00599 }
00600 else if (p==target){
00601 #ifndef WITHOUT_LIMITS
00602 *i = std::numeric_limits<index_type>::max()-1;
00603 #else
00604 *i = static_cast<size_type>(4294967294);
00605 #endif
00606 }
00607 else{
00608
00609 *i = ((p.pos.x() * y_size * z_size) + (p.pos.y() * z_size) + p.pos.z()) * 5;
00610 switch(p.loc){
00611 case CENTER_NODE:
00612 break;
00613
00614 case Z_MAX_x_min_NODE:
00615 *i += 1;
00616 break;
00617
00618 case Z_MAX_x_max_NODE:
00619 *i += 2;
00620 break;
00621
00622 case Z_MAX_y_min_NODE:
00623 *i += 3;
00624 break;
00625
00626 case Z_MAX_y_max_NODE:
00627 *i += 4;
00628 break;
00629
00630 default:
00631 Message::error<<"to_index: unknown value of p.loc"<<Message::done;
00632 }
00633
00634 }
00635 }
00636
00637
00638 template<typename Real>
00639 void
00640 PushRelabel_non_linear_grid_3D<Real>::
00641 from_index(node_type* const p,
00642 const index_type i) const{
00643
00644 #ifndef WITHOUT_LIMITS
00645 if (i==std::numeric_limits<index_type>::max()){
00646 #else
00647 if (i==static_cast<size_type>(4294967295)){
00648 #endif
00649 *p = source;
00650 }
00651 #ifndef WITHOUT_LIMITS
00652 else if (i==std::numeric_limits<index_type>::max()-1){
00653 #else
00654 else if (i==static_cast<size_type>(4294967294)){
00655 #endif
00656 *p = target;
00657 }
00658 else{
00659 switch(i%5){
00660 case 0:
00661 p->loc = CENTER_NODE;
00662 break;
00663 case 1:
00664 p->loc = Z_MAX_x_min_NODE;
00665 break;
00666 case 2:
00667 p->loc = Z_MAX_x_max_NODE;
00668 break;
00669 case 3:
00670 p->loc = Z_MAX_y_min_NODE;
00671 break;
00672 case 4:
00673 p->loc = Z_MAX_y_max_NODE;
00674 break;
00675 }
00676
00677 const index_type j = i/5;
00678 p->pos.x() = j / (y_size * z_size);
00679 p->pos.y() = (j / z_size) - (p->pos.x() * y_size);
00680 p->pos.z() = j % z_size;
00681 }
00682 }
00683
00684
00685 template<typename Real>
00686 const typename PushRelabel_non_linear_grid_3D<Real>::label_type&
00687 PushRelabel_non_linear_grid_3D<Real>::label(const node_type& node) const{
00688
00689 switch(node.loc){
00690 case CENTER_NODE: return label_CENTER[node.pos];
00691 case Z_MAX_x_min_NODE: return label_Z_MAX_x_min[node.pos];
00692 case Z_MAX_x_max_NODE: return label_Z_MAX_x_max[node.pos];
00693 case Z_MAX_y_min_NODE: return label_Z_MAX_y_min[node.pos];
00694 case Z_MAX_y_max_NODE: return label_Z_MAX_y_max[node.pos];
00695 }
00696
00697 Message::error<<"label : non expected exit"<<Message::done;
00698 return label_CENTER[node.pos];
00699 }
00700
00701
00702 template<typename Real>
00703 typename PushRelabel_non_linear_grid_3D<Real>::label_type&
00704 PushRelabel_non_linear_grid_3D<Real>::label(const node_type& node){
00705
00706 switch(node.loc){
00707 case CENTER_NODE: return label_CENTER[node.pos];
00708 case Z_MAX_x_min_NODE: return label_Z_MAX_x_min[node.pos];
00709 case Z_MAX_x_max_NODE: return label_Z_MAX_x_max[node.pos];
00710 case Z_MAX_y_min_NODE: return label_Z_MAX_y_min[node.pos];
00711 case Z_MAX_y_max_NODE: return label_Z_MAX_y_max[node.pos];
00712 }
00713
00714 Message::error<<"label : non expected exit"<<Message::done;
00715 return label_CENTER[node.pos];
00716 }
00717
00718
00719
00723 template<typename Real>
00724 typename PushRelabel_non_linear_grid_3D<Real>::real_type
00725 PushRelabel_non_linear_grid_3D<Real>::
00726 edge_capacity(const node_type& from,
00727 const direction_type dir) const{
00728
00729 switch(from.loc){
00730
00731 case CENTER_NODE:
00732
00733 switch(dir){
00734 case X_MIN_DIR:
00735 return x_major_derivative_edge_capacity(from.pos.x()-1,from.pos.y());
00736
00737 case X_MAX_DIR:
00738 return x_major_derivative_edge_capacity(from.pos.x(),from.pos.y());
00739
00740 case Y_MIN_DIR:
00741 return y_major_derivative_edge_capacity(from.pos.x(),from.pos.y()-1);
00742
00743 case Y_MAX_DIR:
00744 return y_major_derivative_edge_capacity(from.pos.x(),from.pos.y());
00745
00746 case Z_MIN_x_min_DIR:
00747 case Z_MIN_x_max_DIR:
00748 case Z_MIN_y_min_DIR:
00749 case Z_MIN_y_max_DIR:
00750 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z()-1);
00751
00752 case Z_MAX_x_min_DIR:
00753 case Z_MAX_x_max_DIR:
00754 case Z_MAX_y_min_DIR:
00755 case Z_MAX_y_max_DIR:
00756 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00757 }
00758
00759 Message::error<<"edge_capacity : bad end 1"<<Message::done;
00760
00761
00762 case Z_MAX_x_min_NODE:
00763
00764 switch(dir){
00765 case X_MIN_DIR:
00766 return x_minor_derivative_edge_capacity(from.pos.x()-1,from.pos.y());
00767
00768 case X_MAX_DIR:
00769 Message::error<<"edge_capacity : bad end 2.1"<<Message::done;
00770
00771 case Y_MIN_DIR:
00772 Message::error<<"edge_capacity : bad end 2.2"<<Message::done;
00773
00774 case Y_MAX_DIR:
00775 Message::error<<"edge_capacity : bad end 2.3"<<Message::done;
00776
00777 case Z_MIN_x_min_DIR:
00778 Message::error<<"edge_capacity : bad end 2.4"<<Message::done;
00779
00780 case Z_MIN_x_max_DIR:
00781 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00782
00783 case Z_MIN_y_min_DIR:
00784 Message::error<<"edge_capacity : bad end 2.5"<<Message::done;
00785
00786 case Z_MIN_y_max_DIR:
00787 Message::error<<"edge_capacity : bad end 2.6"<<Message::done;
00788
00789 case Z_MAX_x_min_DIR:
00790 Message::error<<"edge_capacity : bad end 2.7"<<Message::done;
00791
00792 case Z_MAX_x_max_DIR:
00793 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00794
00795 case Z_MAX_y_min_DIR:
00796 Message::error<<"edge_capacity : bad end 2.8"<<Message::done;
00797
00798 case Z_MAX_y_max_DIR:
00799 Message::error<<"edge_capacity : bad end 2.9"<<Message::done;
00800 }
00801
00802 Message::error<<"edge_capacity : bad end 2"<<Message::done;
00803
00804
00805 case Z_MAX_x_max_NODE:
00806
00807 switch(dir){
00808 case X_MIN_DIR:
00809 Message::error<<"edge_capacity : bad end 3.1"<<Message::done;
00810
00811 case X_MAX_DIR:
00812 return x_minor_derivative_edge_capacity(from.pos.x(),from.pos.y());
00813
00814 case Y_MIN_DIR:
00815 Message::error<<"edge_capacity : bad end 3.2"<<Message::done;
00816
00817 case Y_MAX_DIR:
00818 Message::error<<"edge_capacity : bad end 3.3"<<Message::done;
00819
00820 case Z_MIN_x_min_DIR:
00821 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00822
00823 case Z_MIN_x_max_DIR:
00824 Message::error<<"edge_capacity : bad end 3.4"<<Message::done;
00825
00826 case Z_MIN_y_min_DIR:
00827 Message::error<<"edge_capacity : bad end 3.5"<<Message::done;
00828
00829 case Z_MIN_y_max_DIR:
00830 Message::error<<"edge_capacity : bad end 3.6"<<Message::done;
00831
00832 case Z_MAX_x_min_DIR:
00833 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00834
00835 case Z_MAX_x_max_DIR:
00836 Message::error<<"edge_capacity : bad end 3.7"<<Message::done;
00837
00838 case Z_MAX_y_min_DIR:
00839 Message::error<<"edge_capacity : bad end 3.8"<<Message::done;;
00840
00841 case Z_MAX_y_max_DIR:
00842 Message::error<<"edge_capacity : bad end 3.9"<<Message::done;
00843 }
00844
00845 Message::error<<"edge_capacity : bad end 3"<<Message::done;
00846 return 0;
00847
00848
00849 case Z_MAX_y_min_NODE:
00850
00851 switch(dir){
00852 case X_MIN_DIR:
00853 Message::error<<"edge_capacity : bad end 4.1"<<Message::done;
00854
00855 case X_MAX_DIR:
00856 Message::error<<"edge_capacity : bad end 4.2"<<Message::done;
00857
00858 case Y_MIN_DIR:
00859 return y_minor_derivative_edge_capacity(from.pos.x(),from.pos.y()-1);
00860
00861 case Y_MAX_DIR:
00862 Message::error<<"edge_capacity : bad end 4.3"<<Message::done;
00863
00864 case Z_MIN_x_min_DIR:
00865 Message::error<<"edge_capacity : bad end 4.4"<<Message::done;
00866
00867 case Z_MIN_x_max_DIR:
00868 Message::error<<"edge_capacity : bad end 4.5"<<Message::done;
00869
00870 case Z_MIN_y_min_DIR:
00871 Message::error<<"edge_capacity : bad end 4.6"<<Message::done;
00872
00873 case Z_MIN_y_max_DIR:
00874 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00875
00876 case Z_MAX_x_min_DIR:
00877 Message::error<<"edge_capacity : bad end 4.7"<<Message::done;
00878
00879 case Z_MAX_x_max_DIR:
00880 Message::error<<"edge_capacity : bad end 4.8"<<Message::done;
00881
00882 case Z_MAX_y_min_DIR:
00883 Message::error<<"edge_capacity : bad end 4.9"<<Message::done;
00884
00885 case Z_MAX_y_max_DIR:
00886 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00887 }
00888
00889 Message::error<<"edge_capacity : bad end 4"<<Message::done;
00890
00891
00892 case Z_MAX_y_max_NODE:
00893
00894 switch(dir){
00895 case X_MIN_DIR:
00896 Message::error<<"edge_capacity : bad end 5.1"<<Message::done;
00897
00898 case X_MAX_DIR:
00899 Message::error<<"edge_capacity : bad end 5.2"<<Message::done;
00900
00901 case Y_MIN_DIR:
00902 Message::error<<"edge_capacity : bad end 5.3"<<Message::done;
00903
00904 case Y_MAX_DIR:
00905 return y_minor_derivative_edge_capacity(from.pos.x(),from.pos.y());
00906
00907 case Z_MIN_x_min_DIR:
00908 Message::error<<"edge_capacity : bad end 5.4"<<Message::done;
00909
00910 case Z_MIN_x_max_DIR:
00911 Message::error<<"edge_capacity : bad end 5.5"<<Message::done;
00912
00913 case Z_MIN_y_min_DIR:
00914 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00915
00916 case Z_MIN_y_max_DIR:
00917 Message::error<<"edge_capacity : bad end 5.6"<<Message::done;
00918
00919 case Z_MAX_x_min_DIR:
00920 Message::error<<"edge_capacity : bad end 5.7"<<Message::done;;
00921
00922 case Z_MAX_x_max_DIR:
00923 Message::error<<"edge_capacity : bad end 5.8"<<Message::done;
00924
00925 case Z_MAX_y_min_DIR:
00926 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z());
00927
00928 case Z_MAX_y_max_DIR:
00929 Message::error<<"edge_capacity : bad end 5.9"<<Message::done;
00930 }
00931
00932 Message::error<<"edge_capacity : bad end 5"<<Message::done;
00933 return 0;
00934 }
00935
00936 }
00937
00938 template<typename Real>
00939 typename PushRelabel_non_linear_grid_3D<Real>::real_type
00940 PushRelabel_non_linear_grid_3D<Real>::
00941 edge_flow(const node_type& from,
00942 const direction_type dir) const{
00943
00944 switch(from.loc){
00945
00946 case CENTER_NODE:
00947
00948 switch(dir){
00949 case X_MIN_DIR:
00950 return -x_major_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z());
00951
00952 case X_MAX_DIR:
00953 return x_major_derivative_edge_flow[from.pos];
00954
00955 case Y_MIN_DIR:
00956 return -y_major_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z());
00957
00958 case Y_MAX_DIR:
00959 return y_major_derivative_edge_flow[from.pos];
00960
00961 case Z_MIN_x_min_DIR:
00962 return -cost_edge_Z_MAX_MAX_x_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
00963
00964 case Z_MIN_x_max_DIR:
00965 return -cost_edge_Z_MAX_MAX_x_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
00966
00967 case Z_MIN_y_min_DIR:
00968 return -cost_edge_Z_MAX_MAX_y_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
00969
00970 case Z_MIN_y_max_DIR:
00971 return -cost_edge_Z_MAX_MAX_y_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
00972
00973 case Z_MAX_x_min_DIR:
00974 return cost_edge_Z_MAX_x_min_flow[from.pos];
00975
00976 case Z_MAX_x_max_DIR:
00977 return cost_edge_Z_MAX_x_max_flow[from.pos];
00978
00979 case Z_MAX_y_min_DIR:
00980 return cost_edge_Z_MAX_y_min_flow[from.pos];
00981
00982 case Z_MAX_y_max_DIR:
00983 return cost_edge_Z_MAX_y_max_flow[from.pos];
00984 }
00985
00986 Message::error<<"edge_flow : bad end 1"<<Message::done;
00987 return 0;
00988
00989
00990 case Z_MAX_x_min_NODE:
00991
00992 switch(dir){
00993 case X_MIN_DIR:
00994 return -x_minor_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z());
00995
00996 case X_MAX_DIR:
00997 Message::error<<"edge_flow : bad end 2.1"<<Message::done;
00998
00999 case Y_MIN_DIR:
01000 Message::error<<"edge_flow : bad end 2.2"<<Message::done;
01001
01002 case Y_MAX_DIR:
01003 Message::error<<"edge_flow : bad end 2.3"<<Message::done;
01004
01005 case Z_MIN_x_min_DIR:
01006 Message::error<<"edge_flow : bad end 2.4"<<Message::done;
01007
01008 case Z_MIN_x_max_DIR:
01009 return -cost_edge_Z_MAX_x_min_flow[from.pos];
01010
01011 case Z_MIN_y_min_DIR:
01012 Message::error<<"edge_flow : bad end 2.5"<<Message::done;
01013
01014 case Z_MIN_y_max_DIR:
01015 Message::error<<"edge_flow : bad end 2.6"<<Message::done;
01016
01017 case Z_MAX_x_min_DIR:
01018 Message::error<<"edge_flow : bad end 2.7"<<Message::done;
01019
01020 case Z_MAX_x_max_DIR:
01021 return cost_edge_Z_MAX_MAX_x_min_flow[from.pos];
01022
01023 case Z_MAX_y_min_DIR:
01024 Message::error<<"edge_flow : bad end 2.8"<<Message::done;
01025
01026 case Z_MAX_y_max_DIR:
01027 Message::error<<"edge_flow : bad end 2.9"<<Message::done;
01028 }
01029
01030 Message::error<<"edge_flow : bad end 2"<<Message::done;
01031 return 0;
01032
01033
01034 case Z_MAX_x_max_NODE:
01035
01036 switch(dir){
01037 case X_MIN_DIR:
01038 Message::error<<"edge_flow : bad end 3.1"<<Message::done;
01039
01040 case X_MAX_DIR:
01041 return x_minor_derivative_edge_flow[from.pos];
01042
01043 case Y_MIN_DIR:
01044 Message::error<<"edge_flow : bad end 3.2"<<Message::done;
01045
01046 case Y_MAX_DIR:
01047 Message::error<<"edge_flow : bad end 3.3"<<Message::done;
01048
01049 case Z_MIN_x_min_DIR:
01050 return -cost_edge_Z_MAX_x_max_flow[from.pos];
01051
01052 case Z_MIN_x_max_DIR:
01053 Message::error<<"edge_flow : bad end 3.4"<<Message::done;
01054
01055 case Z_MIN_y_min_DIR:
01056 Message::error<<"edge_flow : bad end 3.5"<<Message::done;
01057
01058 case Z_MIN_y_max_DIR:
01059 Message::error<<"edge_flow : bad end 3.6"<<Message::done;
01060
01061 case Z_MAX_x_min_DIR:
01062 return cost_edge_Z_MAX_MAX_x_max_flow[from.pos];
01063
01064 case Z_MAX_x_max_DIR:
01065 Message::error<<"edge_flow : bad end 3.7"<<Message::done;
01066
01067 case Z_MAX_y_min_DIR:
01068 Message::error<<"edge_flow : bad end 3.8"<<Message::done;
01069
01070 case Z_MAX_y_max_DIR:
01071 Message::error<<"edge_flow : bad end 3.9"<<Message::done;
01072 }
01073
01074 Message::error<<"edge_flow : bad end 3"<<Message::done;
01075 return 0;
01076
01077
01078 case Z_MAX_y_min_NODE:
01079
01080 switch(dir){
01081 case X_MIN_DIR:
01082 Message::error<<"edge_flow : bad end 4.1"<<Message::done;
01083
01084 case X_MAX_DIR:
01085 Message::error<<"edge_flow : bad end 4.2"<<Message::done;
01086
01087 case Y_MIN_DIR:
01088 return -y_minor_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z());
01089
01090
01091 case Y_MAX_DIR:
01092 Message::error<<"edge_flow : bad end 4.3"<<Message::done;
01093
01094 case Z_MIN_x_min_DIR:
01095 Message::error<<"edge_flow : bad end 4.4"<<Message::done;
01096
01097 case Z_MIN_x_max_DIR:
01098 Message::error<<"edge_flow : bad end 4.5"<<Message::done;
01099
01100 case Z_MIN_y_min_DIR:
01101 Message::error<<"edge_flow : bad end 4.6"<<Message::done;
01102
01103 case Z_MIN_y_max_DIR:
01104 return -cost_edge_Z_MAX_y_min_flow[from.pos];
01105
01106 case Z_MAX_x_min_DIR:
01107 Message::error<<"edge_flow : bad end 4.7"<<Message::done;
01108
01109 case Z_MAX_x_max_DIR:
01110 Message::error<<"edge_flow : bad end 4.8"<<Message::done;
01111
01112 case Z_MAX_y_min_DIR:
01113 Message::error<<"edge_flow : bad end 4.9"<<Message::done;
01114
01115 case Z_MAX_y_max_DIR:
01116 return cost_edge_Z_MAX_MAX_y_min_flow[from.pos];
01117
01118 }
01119
01120 Message::error<<"edge_flow : bad end 4"<<Message::done;
01121 return 0;
01122
01123
01124 case Z_MAX_y_max_NODE:
01125
01126 switch(dir){
01127 case X_MIN_DIR:
01128 Message::error<<"edge_flow : bad end 5.1"<<Message::done;
01129
01130 case X_MAX_DIR:
01131 Message::error<<"edge_flow : bad end 5.2"<<Message::done;
01132
01133 case Y_MIN_DIR:
01134 Message::error<<"edge_flow : bad end 5.3"<<Message::done;
01135
01136 case Y_MAX_DIR:
01137 return y_minor_derivative_edge_flow[from.pos];
01138
01139 case Z_MIN_x_min_DIR:
01140 Message::error<<"edge_flow : bad end 5.4"<<Message::done;
01141
01142 case Z_MIN_x_max_DIR:
01143 Message::error<<"edge_flow : bad end 5.5"<<Message::done;
01144
01145 case Z_MIN_y_min_DIR:
01146 return -cost_edge_Z_MAX_y_max_flow[from.pos];
01147
01148 case Z_MIN_y_max_DIR:
01149 Message::error<<"edge_flow : bad end 5.6"<<Message::done;
01150
01151 case Z_MAX_x_min_DIR:
01152 Message::error<<"edge_flow : bad end 5.7"<<Message::done;
01153
01154 case Z_MAX_x_max_DIR:
01155 Message::error<<"edge_flow : bad end 5.8"<<Message::done;
01156
01157 case Z_MAX_y_min_DIR:
01158 return cost_edge_Z_MAX_MAX_y_max_flow[from.pos];
01159
01160 case Z_MAX_y_max_DIR:
01161 Message::error<<"edge_flow : bad end 5.9"<<Message::done;
01162 }
01163
01164 Message::error<<"edge_flow : bad end 5"<<Message::done;
01165 return 0;
01166 }
01167
01168 Message::error<<"edge_flow : bad end 0"<<Message::done;
01169 return 0;
01170 }
01171
01172
01173 template<typename Real>
01174 typename PushRelabel_non_linear_grid_3D<Real>::real_type
01175 PushRelabel_non_linear_grid_3D<Real>::
01176 edge_residual(const node_type& from,
01177 const direction_type dir) const{
01178
01179 switch(from.loc){
01180
01181 case CENTER_NODE:
01182
01183 switch(dir){
01184 case X_MIN_DIR:
01185 return x_major_derivative_edge_capacity(from.pos.x()-1,from.pos.y())
01186 + x_major_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z());
01187
01188 case X_MAX_DIR:
01189 return x_major_derivative_edge_capacity(from.pos.x(),from.pos.y())
01190 - x_major_derivative_edge_flow[from.pos];
01191
01192 case Y_MIN_DIR:
01193 return y_major_derivative_edge_capacity(from.pos.x(),from.pos.y()-1)
01194 + y_major_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z());
01195
01196 case Y_MAX_DIR:
01197 return y_major_derivative_edge_capacity(from.pos.x(),from.pos.y())
01198 - y_major_derivative_edge_flow[from.pos];
01199
01200 case Z_MIN_x_min_DIR:
01201 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z()-1)
01202 + cost_edge_Z_MAX_MAX_x_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
01203
01204 case Z_MIN_x_max_DIR:
01205 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z()-1)
01206 + cost_edge_Z_MAX_MAX_x_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
01207
01208 case Z_MIN_y_min_DIR:
01209 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z()-1)
01210 + cost_edge_Z_MAX_MAX_y_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
01211
01212 case Z_MIN_y_max_DIR:
01213 return cost_edge_capacity(from.pos.x(),from.pos.y(),from.pos.z()-1)
01214 + cost_edge_Z_MAX_MAX_y_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1);
01215
01216 case Z_MAX_x_min_DIR:
01217 return cost_edge_capacity[from.pos]
01218 - cost_edge_Z_MAX_x_min_flow[from.pos];
01219
01220 case Z_MAX_x_max_DIR:
01221 return cost_edge_capacity[from.pos]
01222 - cost_edge_Z_MAX_x_max_flow[from.pos];
01223
01224 case Z_MAX_y_min_DIR:
01225 return cost_edge_capacity[from.pos]
01226 - cost_edge_Z_MAX_y_min_flow[from.pos];
01227
01228 case Z_MAX_y_max_DIR:
01229 return cost_edge_capacity[from.pos]
01230 - cost_edge_Z_MAX_y_max_flow[from.pos];
01231 }
01232
01233 Message::error<<"edge_residual : bad end 1"<<Message::done;
01234 return 0;
01235
01236
01237 case Z_MAX_x_min_NODE:
01238
01239 switch(dir){
01240 case X_MIN_DIR:
01241 return x_minor_derivative_edge_capacity(from.pos.x()-1,from.pos.y())
01242 + x_minor_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z());
01243
01244 case X_MAX_DIR:
01245 Message::error<<"edge_residual : bad end 2.1"<<Message::done;
01246
01247 case Y_MIN_DIR:
01248 Message::error<<"edge_residual : bad end 2.2"<<Message::done;
01249
01250 case Y_MAX_DIR:
01251 Message::error<<"edge_residual : bad end 2.3"<<Message::done;
01252
01253 case Z_MIN_x_min_DIR:
01254 Message::error<<"edge_residual : bad end 2.4"<<Message::done;
01255
01256 case Z_MIN_x_max_DIR:
01257 return cost_edge_capacity[from.pos]
01258 + cost_edge_Z_MAX_x_min_flow[from.pos];
01259
01260 case Z_MIN_y_min_DIR:
01261 Message::error<<"edge_residual : bad end 2.5"<<Message::done;
01262
01263 case Z_MIN_y_max_DIR:
01264 Message::error<<"edge_residual : bad end 2.6"<<Message::done;
01265
01266 case Z_MAX_x_min_DIR:
01267 Message::error<<"edge_residual : bad end 2.7"<<Message::done;
01268
01269 case Z_MAX_x_max_DIR:
01270 return cost_edge_capacity[from.pos]
01271 - cost_edge_Z_MAX_MAX_x_min_flow[from.pos];
01272
01273 case Z_MAX_y_min_DIR:
01274 Message::error<<"edge_residual : bad end 2.8"<<Message::done;
01275
01276 case Z_MAX_y_max_DIR:
01277 Message::error<<"edge_residual : bad end 2.9"<<Message::done;
01278 }
01279
01280 Message::error<<"edge_residual : bad end 2"<<Message::done;
01281 return 0;
01282
01283
01284 case Z_MAX_x_max_NODE:
01285
01286 switch(dir){
01287 case X_MIN_DIR:
01288 Message::error<<"edge_residual : bad end 3.1"<<Message::done;
01289
01290 case X_MAX_DIR:
01291 return x_minor_derivative_edge_capacity(from.pos.x(),from.pos.y())
01292 - x_minor_derivative_edge_flow[from.pos];
01293
01294 case Y_MIN_DIR:
01295 Message::error<<"edge_residual : bad end 3.2"<<Message::done;
01296
01297 case Y_MAX_DIR:
01298 Message::error<<"edge_residual : bad end 3.3"<<Message::done;
01299
01300 case Z_MIN_x_min_DIR:
01301 return cost_edge_capacity[from.pos]
01302 + cost_edge_Z_MAX_x_max_flow[from.pos];
01303
01304 case Z_MIN_x_max_DIR:
01305 Message::error<<"edge_residual : bad end 3.4"<<Message::done;
01306
01307 case Z_MIN_y_min_DIR:
01308 Message::error<<"edge_residual : bad end 3.5"<<Message::done;
01309
01310 case Z_MIN_y_max_DIR:
01311 Message::error<<"edge_residual : bad end 3.6"<<Message::done;
01312
01313 case Z_MAX_x_min_DIR:
01314 return cost_edge_capacity[from.pos]
01315 - cost_edge_Z_MAX_MAX_x_max_flow[from.pos];
01316
01317 case Z_MAX_x_max_DIR:
01318 Message::error<<"edge_residual : bad end 3.7"<<Message::done;
01319
01320 case Z_MAX_y_min_DIR:
01321 Message::error<<"edge_residual : bad end 3.8"<<Message::done;
01322
01323 case Z_MAX_y_max_DIR:
01324 Message::error<<"edge_residual : bad end 3.9"<<Message::done;
01325 }
01326
01327 Message::error<<"edge_residual : bad end 3"<<Message::done;
01328 return 0;
01329
01330
01331 case Z_MAX_y_min_NODE:
01332
01333 switch(dir){
01334 case X_MIN_DIR:
01335 Message::error<<"edge_residual : bad end 4.1"<<Message::done;
01336
01337 case X_MAX_DIR:
01338 Message::error<<"edge_residual : bad end 4.2"<<Message::done;
01339
01340 case Y_MIN_DIR:
01341 return y_minor_derivative_edge_capacity(from.pos.x(),from.pos.y()-1)
01342 + y_minor_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z());
01343
01344
01345 case Y_MAX_DIR:
01346 Message::error<<"edge_residual : bad end 4.3"<<Message::done;
01347
01348 case Z_MIN_x_min_DIR:
01349 Message::error<<"edge_residual : bad end 4.4"<<Message::done;
01350
01351 case Z_MIN_x_max_DIR:
01352 Message::error<<"edge_residual : bad end 4.5"<<Message::done;
01353
01354 case Z_MIN_y_min_DIR:
01355 Message::error<<"edge_residual : bad end 4.6"<<Message::done;
01356
01357 case Z_MIN_y_max_DIR:
01358 return cost_edge_capacity[from.pos]
01359 + cost_edge_Z_MAX_y_min_flow[from.pos];
01360
01361 case Z_MAX_x_min_DIR:
01362 Message::error<<"edge_residual : bad end 4.7"<<Message::done;
01363
01364 case Z_MAX_x_max_DIR:
01365 Message::error<<"edge_residual : bad end 4.8"<<Message::done;
01366
01367 case Z_MAX_y_min_DIR:
01368 Message::error<<"edge_residual : bad end 4.9"<<Message::done;
01369
01370 case Z_MAX_y_max_DIR:
01371 return cost_edge_capacity[from.pos]
01372 - cost_edge_Z_MAX_MAX_y_min_flow[from.pos];
01373
01374 }
01375
01376 Message::error<<"edge_residual : bad end 4"<<Message::done;
01377 return 0;
01378
01379
01380 case Z_MAX_y_max_NODE:
01381
01382 switch(dir){
01383 case X_MIN_DIR:
01384 Message::error<<"edge_residual : bad end 5.1"<<Message::done;
01385
01386 case X_MAX_DIR:
01387 Message::error<<"edge_residual : bad end 5.2"<<Message::done;
01388
01389 case Y_MIN_DIR:
01390 Message::error<<"edge_residual : bad end 5.3"<<Message::done;
01391
01392 case Y_MAX_DIR:
01393 return y_minor_derivative_edge_capacity(from.pos.x(),from.pos.y())
01394 - y_minor_derivative_edge_flow[from.pos];
01395
01396 case Z_MIN_x_min_DIR:
01397 Message::error<<"edge_residual : bad end 5.4"<<Message::done;
01398
01399 case Z_MIN_x_max_DIR:
01400 Message::error<<"edge_residual : bad end 5.5"<<Message::done;
01401
01402 case Z_MIN_y_min_DIR:
01403 return cost_edge_capacity[from.pos]
01404 + cost_edge_Z_MAX_y_max_flow[from.pos];
01405
01406 case Z_MIN_y_max_DIR:
01407 Message::error<<"edge_residual : bad end 5.6"<<Message::done;
01408
01409 case Z_MAX_x_min_DIR:
01410 Message::error<<"edge_residual : bad end 5.7"<<Message::done;
01411
01412 case Z_MAX_x_max_DIR:
01413 Message::error<<"edge_residual : bad end 5.8"<<Message::done;
01414
01415 case Z_MAX_y_min_DIR:
01416 return cost_edge_capacity[from.pos]
01417 - cost_edge_Z_MAX_MAX_y_max_flow[from.pos];
01418
01419 case Z_MAX_y_max_DIR:
01420 Message::error<<"edge_residual : bad end 5.9"<<Message::done;
01421 }
01422
01423 Message::error<<"edge_residual : bad end 5"<<Message::done;
01424 return 0;
01425 }
01426
01427 Message::error<<"edge_residual : bad end 0"<<Message::done;
01428 return 0;
01429 }
01430
01431
01432
01436 template<typename Real>
01437 void
01438 PushRelabel_non_linear_grid_3D<Real>::
01439 neighbors(const node_type& node,
01440 node_type* const result,
01441 size_type* const n_res) const{
01442
01443 *n_res = 0;
01444
01445
01446
01447 switch(node.loc){
01448
01449 case CENTER_NODE:
01450 if ((node.pos.x()>0)
01451 &&(domain(node.pos.x()-1,node.pos.y(),node.pos.z())||
01452 ((node.pos.z()>0)&&
01453 domain(node.pos.x()-1,node.pos.y(),node.pos.z()-1)))){
01454 result[*n_res].pos.x() = node.pos.x()-1;
01455 result[*n_res].pos.y() = node.pos.y();
01456 result[*n_res].pos.z() = node.pos.z();
01457 result[*n_res].loc = CENTER_NODE;
01458 (*n_res)++;
01459 }
01460 break;
01461
01462 case Z_MAX_x_min_NODE:
01463 if ((node.pos.x()>0)
01464 &&(domain(node.pos.x()-1,node.pos.y(),node.pos.z()))){
01465 result[*n_res].pos.x() = node.pos.x()-1;
01466 result[*n_res].pos.y() = node.pos.y();
01467 result[*n_res].pos.z() = node.pos.z();
01468 result[*n_res].loc = Z_MAX_x_max_NODE;
01469 (*n_res)++;
01470 }
01471 break;
01472
01473 case Z_MAX_x_max_NODE:
01474 case Z_MAX_y_min_NODE:
01475 case Z_MAX_y_max_NODE:
01476 break;
01477 }
01478
01479
01480
01481 switch(node.loc){
01482
01483 case CENTER_NODE:
01484 if ((node.pos.y()>0)
01485 &&(domain(node.pos.x(),node.pos.y()-1,node.pos.z())||
01486 ((node.pos.z()>0)&&
01487 domain(node.pos.x(),node.pos.y()-1,node.pos.z()-1)))){
01488 result[*n_res].pos.x() = node.pos.x();
01489 result[*n_res].pos.y() = node.pos.y()-1;
01490 result[*n_res].pos.z() = node.pos.z();
01491 result[*n_res].loc = CENTER_NODE;
01492 (*n_res)++;
01493 }
01494 break;
01495
01496 case Z_MAX_y_min_NODE:
01497 if ((node.pos.y()>0)
01498 &&(domain(node.pos.x(),node.pos.y()-1,node.pos.z()))){
01499 result[*n_res].pos.x() = node.pos.x();
01500 result[*n_res].pos.y() = node.pos.y()-1;
01501 result[*n_res].pos.z() = node.pos.z();
01502 result[*n_res].loc = Z_MAX_y_max_NODE;
01503 (*n_res)++;
01504 }
01505 break;
01506
01507 case Z_MAX_x_max_NODE:
01508 case Z_MAX_x_min_NODE:
01509 case Z_MAX_y_max_NODE:
01510 break;
01511 }
01512
01513
01514
01515 switch(node.loc){
01516
01517 case CENTER_NODE:
01518
01519 if ((node.pos.z()>0)
01520 &&(domain(node.pos.x(),node.pos.y(),node.pos.z()-1))){
01521
01522 result[*n_res].pos.x() = node.pos.x();
01523 result[*n_res].pos.y() = node.pos.y();
01524 result[*n_res].pos.z() = node.pos.z()-1;
01525 result[*n_res].loc = Z_MAX_x_min_NODE;
01526 (*n_res)++;
01527
01528 result[*n_res].pos.x() = node.pos.x();
01529 result[*n_res].pos.y() = node.pos.y();
01530 result[*n_res].pos.z() = node.pos.z()-1;
01531 result[*n_res].loc = Z_MAX_x_max_NODE;
01532 (*n_res)++;
01533
01534 result[*n_res].pos.x() = node.pos.x();
01535 result[*n_res].pos.y() = node.pos.y();
01536 result[*n_res].pos.z() = node.pos.z()-1;
01537 result[*n_res].loc = Z_MAX_y_min_NODE;
01538 (*n_res)++;
01539
01540 result[*n_res].pos.x() = node.pos.x();
01541 result[*n_res].pos.y() = node.pos.y();
01542 result[*n_res].pos.z() = node.pos.z()-1;
01543 result[*n_res].loc = Z_MAX_y_max_NODE;
01544 (*n_res)++;
01545 }
01546 break;
01547
01548 case Z_MAX_x_min_NODE:
01549 result[*n_res].pos.x() = node.pos.x();
01550 result[*n_res].pos.y() = node.pos.y();
01551 result[*n_res].pos.z() = node.pos.z();
01552 result[*n_res].loc = CENTER_NODE;
01553 (*n_res)++;
01554 break;
01555
01556 case Z_MAX_x_max_NODE:
01557 result[*n_res].pos.x() = node.pos.x();
01558 result[*n_res].pos.y() = node.pos.y();
01559 result[*n_res].pos.z() = node.pos.z();
01560 result[*n_res].loc = CENTER_NODE;
01561 (*n_res)++;
01562 break;
01563
01564 case Z_MAX_y_min_NODE:
01565 result[*n_res].pos.x() = node.pos.x();
01566 result[*n_res].pos.y() = node.pos.y();
01567 result[*n_res].pos.z() = node.pos.z();
01568 result[*n_res].loc = CENTER_NODE;
01569 (*n_res)++;
01570 break;
01571
01572 case Z_MAX_y_max_NODE:
01573 result[*n_res].pos.x() = node.pos.x();
01574 result[*n_res].pos.y() = node.pos.y();
01575 result[*n_res].pos.z() = node.pos.z();
01576 result[*n_res].loc = CENTER_NODE;
01577 (*n_res)++;
01578 break;
01579 }
01580
01581
01582
01583 switch(node.loc){
01584
01585 case CENTER_NODE:
01586
01587 if ((node.pos.x()<x_size-1)
01588 &&(domain(node.pos.x()+1,node.pos.y(),node.pos.z())||
01589 ((node.pos.z()>0)&&
01590 domain(node.pos.x()+1,node.pos.y(),node.pos.z()-1)))){
01591
01592 result[*n_res].pos.x() = node.pos.x()+1;
01593 result[*n_res].pos.y() = node.pos.y();
01594 result[*n_res].pos.z() = node.pos.z();
01595 result[*n_res].loc = CENTER_NODE;
01596 (*n_res)++;
01597 }
01598 break;
01599
01600 case Z_MAX_x_max_NODE:
01601
01602 if ((node.pos.x()<x_size-1)
01603 &&(domain(node.pos.x()+1,node.pos.y(),node.pos.z()))){
01604 result[*n_res].pos.x() = node.pos.x()+1;
01605 result[*n_res].pos.y() = node.pos.y();
01606 result[*n_res].pos.z() = node.pos.z();
01607 result[*n_res].loc = Z_MAX_x_min_NODE;
01608 (*n_res)++;
01609 }
01610 break;
01611
01612 case Z_MAX_x_min_NODE:
01613 case Z_MAX_y_min_NODE:
01614 case Z_MAX_y_max_NODE:
01615 break;
01616 }
01617
01618
01619
01620 switch(node.loc){
01621
01622 case CENTER_NODE:
01623
01624 if ((node.pos.y()<y_size-1)
01625 &&(domain(node.pos.x(),node.pos.y()+1,node.pos.z())||
01626 ((node.pos.z()>0)&&
01627 domain(node.pos.x(),node.pos.y()+1,node.pos.z()-1)))){
01628
01629 result[*n_res].pos.x() = node.pos.x();
01630 result[*n_res].pos.y() = node.pos.y()+1;
01631 result[*n_res].pos.z() = node.pos.z();
01632 result[*n_res].loc = CENTER_NODE;
01633 (*n_res)++;
01634 }
01635 break;
01636
01637 case Z_MAX_y_max_NODE:
01638 if ((node.pos.y()<y_size-1)
01639 &&(domain(node.pos.x(),node.pos.y()+1,node.pos.z()))){
01640
01641 result[*n_res].pos.x() = node.pos.x();
01642 result[*n_res].pos.y() = node.pos.y()+1;
01643 result[*n_res].pos.z() = node.pos.z();
01644 result[*n_res].loc = Z_MAX_y_min_NODE;
01645 (*n_res)++;
01646 }
01647 break;
01648
01649 case Z_MAX_x_min_NODE:
01650 case Z_MAX_x_max_NODE:
01651 case Z_MAX_y_min_NODE:
01652 break;
01653 }
01654
01655
01656
01657
01658 switch(node.loc){
01659 case CENTER_NODE:
01660 if ((node.pos.z()<z_size)
01661 &&(domain[node.pos])){
01662
01663 result[*n_res].pos.x() = node.pos.x();
01664 result[*n_res].pos.y() = node.pos.y();
01665 result[*n_res].pos.z() = node.pos.z();
01666 result[*n_res].loc = Z_MAX_x_min_NODE;
01667 (*n_res)++;
01668
01669 result[*n_res].pos.x() = node.pos.x();
01670 result[*n_res].pos.y() = node.pos.y();
01671 result[*n_res].pos.z() = node.pos.z();
01672 result[*n_res].loc = Z_MAX_x_max_NODE;
01673 (*n_res)++;
01674
01675 result[*n_res].pos.x() = node.pos.x();
01676 result[*n_res].pos.y() = node.pos.y();
01677 result[*n_res].pos.z() = node.pos.z();
01678 result[*n_res].loc = Z_MAX_y_min_NODE;
01679 (*n_res)++;
01680
01681 result[*n_res].pos.x() = node.pos.x();
01682 result[*n_res].pos.y() = node.pos.y();
01683 result[*n_res].pos.z() = node.pos.z();
01684 result[*n_res].loc = Z_MAX_y_max_NODE;
01685 (*n_res)++;
01686 }
01687 break;
01688
01689 case Z_MAX_x_min_NODE:
01690 if ((node.pos.z()<z_size)
01691 &&(domain[node.pos])){
01692
01693 result[*n_res].pos.x() = node.pos.x();
01694 result[*n_res].pos.y() = node.pos.y();
01695 result[*n_res].pos.z() = node.pos.z()+1;
01696 result[*n_res].loc = CENTER_NODE;
01697 (*n_res)++;
01698 }
01699 break;
01700
01701 case Z_MAX_x_max_NODE:
01702 if ((node.pos.z()<z_size)
01703 &&(domain[node.pos])){
01704
01705 result[*n_res].pos.x() = node.pos.x();
01706 result[*n_res].pos.y() = node.pos.y();
01707 result[*n_res].pos.z() = node.pos.z()+1;
01708 result[*n_res].loc = CENTER_NODE;
01709 (*n_res)++;
01710 }
01711 break;
01712
01713 case Z_MAX_y_min_NODE:
01714 if ((node.pos.z()<z_size)
01715 &&(domain[node.pos])){
01716
01717 result[*n_res].pos.x() = node.pos.x();
01718 result[*n_res].pos.y() = node.pos.y();
01719 result[*n_res].pos.z() = node.pos.z()+1;
01720 result[*n_res].loc = CENTER_NODE;
01721 (*n_res)++;
01722 }
01723 break;
01724
01725 case Z_MAX_y_max_NODE:
01726 if ((node.pos.z()<z_size)
01727 &&(domain[node.pos])){
01728
01729 result[*n_res].pos.x() = node.pos.x();
01730 result[*n_res].pos.y() = node.pos.y();
01731 result[*n_res].pos.z() = node.pos.z()+1;
01732 result[*n_res].loc = CENTER_NODE;
01733 (*n_res)++;
01734 }
01735 break;
01736 }
01737 }
01738
01739
01740
01741
01742
01746 template<typename Real>
01747 void
01748 PushRelabel_non_linear_grid_3D<Real>::
01749 adjacent_nodes(const node_type& node,
01750 node_type* const result,
01751 size_type* const n_res) const{
01752
01753 *n_res = 0;
01754
01755
01756
01757 switch(node.loc){
01758
01759 case CENTER_NODE:
01760 if ((node.pos.x()>0)
01761 &&(domain(node.pos.x()-1,node.pos.y(),node.pos.z())||
01762 ((node.pos.z()>0)&&
01763 domain(node.pos.x()-1,node.pos.y(),node.pos.z()-1)))){
01764
01765 if(edge_residual(node,X_MIN_DIR)>epsilon){
01766 result[*n_res].pos.x() = node.pos.x()-1;
01767 result[*n_res].pos.y() = node.pos.y();
01768 result[*n_res].pos.z() = node.pos.z();
01769 result[*n_res].loc = CENTER_NODE;
01770 (*n_res)++;
01771 }
01772 }
01773 break;
01774
01775 case Z_MAX_x_min_NODE:
01776 if ((node.pos.x()>0)
01777 &&(domain(node.pos.x()-1,node.pos.y(),node.pos.z()))){
01778 if(edge_residual(node,X_MIN_DIR)>epsilon){
01779 result[*n_res].pos.x() = node.pos.x()-1;
01780 result[*n_res].pos.y() = node.pos.y();
01781 result[*n_res].pos.z() = node.pos.z();
01782 result[*n_res].loc = Z_MAX_x_max_NODE;
01783 (*n_res)++;
01784 }
01785 }
01786 break;
01787
01788 case Z_MAX_x_max_NODE:
01789 case Z_MAX_y_min_NODE:
01790 case Z_MAX_y_max_NODE:
01791 break;
01792 }
01793
01794
01795
01796 switch(node.loc){
01797
01798 case CENTER_NODE:
01799 if ((node.pos.y()>0)
01800 &&(domain(node.pos.x(),node.pos.y()-1,node.pos.z())||
01801 ((node.pos.z()>0)&&
01802 domain(node.pos.x(),node.pos.y()-1,node.pos.z()-1)))){
01803 if(edge_residual(node,Y_MIN_DIR)>epsilon){
01804 result[*n_res].pos.x() = node.pos.x();
01805 result[*n_res].pos.y() = node.pos.y()-1;
01806 result[*n_res].pos.z() = node.pos.z();
01807 result[*n_res].loc = CENTER_NODE;
01808 (*n_res)++;
01809 }
01810 }
01811 break;
01812
01813 case Z_MAX_y_min_NODE:
01814 if ((node.pos.y()>0)
01815 &&(domain(node.pos.x(),node.pos.y()-1,node.pos.z()))){
01816 if(edge_residual(node,Y_MIN_DIR)>epsilon){
01817 result[*n_res].pos.x() = node.pos.x();
01818 result[*n_res].pos.y() = node.pos.y()-1;
01819 result[*n_res].pos.z() = node.pos.z();
01820 result[*n_res].loc = Z_MAX_y_max_NODE;
01821 (*n_res)++;
01822 }
01823 }
01824 break;
01825
01826 case Z_MAX_x_max_NODE:
01827 case Z_MAX_x_min_NODE:
01828 case Z_MAX_y_max_NODE:
01829 break;
01830 }
01831
01832
01833
01834
01835 switch(node.loc){
01836
01837 case CENTER_NODE:
01838
01839 if ((node.pos.z()>0)
01840 &&(domain(node.pos.x(),node.pos.y(),node.pos.z()-1))){
01841
01842 if (edge_residual(node,Z_MIN_x_min_DIR)>epsilon){
01843 result[*n_res].pos.x() = node.pos.x();
01844 result[*n_res].pos.y() = node.pos.y();
01845 result[*n_res].pos.z() = node.pos.z()-1;
01846 result[*n_res].loc = Z_MAX_x_min_NODE;
01847 (*n_res)++;
01848 }
01849
01850 if (edge_residual(node,Z_MIN_x_max_DIR)>epsilon){
01851 result[*n_res].pos.x() = node.pos.x();
01852 result[*n_res].pos.y() = node.pos.y();
01853 result[*n_res].pos.z() = node.pos.z()-1;
01854 result[*n_res].loc = Z_MAX_x_max_NODE;
01855 (*n_res)++;
01856 }
01857
01858 if (edge_residual(node,Z_MIN_y_min_DIR)>epsilon){
01859 result[*n_res].pos.x() = node.pos.x();
01860 result[*n_res].pos.y() = node.pos.y();
01861 result[*n_res].pos.z() = node.pos.z()-1;
01862 result[*n_res].loc = Z_MAX_y_min_NODE;
01863 (*n_res)++;
01864 }
01865
01866 if (edge_residual(node,Z_MIN_y_max_DIR)>epsilon){
01867 result[*n_res].pos.x() = node.pos.x();
01868 result[*n_res].pos.y() = node.pos.y();
01869 result[*n_res].pos.z() = node.pos.z()-1;
01870 result[*n_res].loc = Z_MAX_y_max_NODE;
01871 (*n_res)++;
01872 }
01873 }
01874 break;
01875
01876 case Z_MAX_x_min_NODE:
01877 if (edge_residual(node,Z_MIN_x_max_DIR)>epsilon){
01878 result[*n_res].pos.x() = node.pos.x();
01879 result[*n_res].pos.y() = node.pos.y();
01880 result[*n_res].pos.z() = node.pos.z();
01881 result[*n_res].loc = CENTER_NODE;
01882 (*n_res)++;
01883 }
01884 break;
01885
01886 case Z_MAX_x_max_NODE:
01887 if (edge_residual(node,Z_MIN_x_min_DIR)>epsilon){
01888 result[*n_res].pos.x() = node.pos.x();
01889 result[*n_res].pos.y() = node.pos.y();
01890 result[*n_res].pos.z() = node.pos.z();
01891 result[*n_res].loc = CENTER_NODE;
01892 (*n_res)++;
01893 }
01894 break;
01895
01896 case Z_MAX_y_min_NODE:
01897 if (edge_residual(node,Z_MIN_y_max_DIR)>epsilon){
01898 result[*n_res].pos.x() = node.pos.x();
01899 result[*n_res].pos.y() = node.pos.y();
01900 result[*n_res].pos.z() = node.pos.z();
01901 result[*n_res].loc = CENTER_NODE;
01902 (*n_res)++;
01903 }
01904 break;
01905
01906 case Z_MAX_y_max_NODE:
01907 if (edge_residual(node,Z_MIN_y_min_DIR)>epsilon){
01908 result[*n_res].pos.x() = node.pos.x();
01909 result[*n_res].pos.y() = node.pos.y();
01910 result[*n_res].pos.z() = node.pos.z();
01911 result[*n_res].loc = CENTER_NODE;
01912 (*n_res)++;
01913 }
01914 break;
01915 }
01916
01917
01918
01919 switch(node.loc){
01920
01921 case CENTER_NODE:
01922 if ((node.pos.x()<x_size-1)
01923 &&(domain(node.pos.x()+1,node.pos.y(),node.pos.z())||
01924 ((node.pos.z()>0)&&
01925 domain(node.pos.x()+1,node.pos.y(),node.pos.z()-1)))){
01926 if (edge_residual(node,X_MAX_DIR)>epsilon){
01927
01928 result[*n_res].pos.x() = node.pos.x()+1;
01929 result[*n_res].pos.y() = node.pos.y();
01930 result[*n_res].pos.z() = node.pos.z();
01931 result[*n_res].loc = CENTER_NODE;
01932 (*n_res)++;
01933 }
01934 }
01935 break;
01936
01937 case Z_MAX_x_max_NODE:
01938 if ((node.pos.x()<x_size-1)
01939 &&(domain(node.pos.x()+1,node.pos.y(),node.pos.z()))){
01940 if (edge_residual(node,X_MAX_DIR)>epsilon){
01941
01942 result[*n_res].pos.x() = node.pos.x()+1;
01943 result[*n_res].pos.y() = node.pos.y();
01944 result[*n_res].pos.z() = node.pos.z();
01945 result[*n_res].loc = Z_MAX_x_min_NODE;
01946 (*n_res)++;
01947 }
01948 }
01949 break;
01950
01951 case Z_MAX_x_min_NODE:
01952 case Z_MAX_y_min_NODE:
01953 case Z_MAX_y_max_NODE:
01954 break;
01955 }
01956
01957
01958
01959 switch(node.loc){
01960
01961 case CENTER_NODE:
01962 if ((node.pos.y()<y_size-1)
01963 &&(domain(node.pos.x(),node.pos.y()+1,node.pos.z())||
01964 ((node.pos.z()>0)&&
01965 domain(node.pos.x(),node.pos.y()+1,node.pos.z()-1)))){
01966
01967 if (edge_residual(node,Y_MAX_DIR)>epsilon){
01968
01969 result[*n_res].pos.x() = node.pos.x();
01970 result[*n_res].pos.y() = node.pos.y()+1;
01971 result[*n_res].pos.z() = node.pos.z();
01972 result[*n_res].loc = CENTER_NODE;
01973 (*n_res)++;
01974 }
01975 }
01976 break;
01977
01978 case Z_MAX_y_max_NODE:
01979 if ((node.pos.y()<y_size-1)
01980 &&(domain(node.pos.x(),node.pos.y()+1,node.pos.z()))){
01981 if (edge_residual(node,Y_MAX_DIR)>epsilon){
01982
01983 result[*n_res].pos.x() = node.pos.x();
01984 result[*n_res].pos.y() = node.pos.y()+1;
01985 result[*n_res].pos.z() = node.pos.z();
01986 result[*n_res].loc = Z_MAX_y_min_NODE;
01987 (*n_res)++;
01988 }
01989 }
01990 break;
01991
01992 case Z_MAX_x_min_NODE:
01993 case Z_MAX_x_max_NODE:
01994 case Z_MAX_y_min_NODE:
01995 break;
01996 }
01997
01998
01999
02000
02001 switch(node.loc){
02002 case CENTER_NODE:
02003 if ((node.pos.z()<z_size)
02004 &&(domain[node.pos])){
02005
02006 if (edge_residual(node,Z_MAX_x_min_DIR)>epsilon){
02007
02008 result[*n_res].pos.x() = node.pos.x();
02009 result[*n_res].pos.y() = node.pos.y();
02010 result[*n_res].pos.z() = node.pos.z();
02011 result[*n_res].loc = Z_MAX_x_min_NODE;
02012 (*n_res)++;
02013 }
02014
02015 if (edge_residual(node,Z_MAX_x_max_DIR)>epsilon){
02016
02017 result[*n_res].pos.x() = node.pos.x();
02018 result[*n_res].pos.y() = node.pos.y();
02019 result[*n_res].pos.z() = node.pos.z();
02020 result[*n_res].loc = Z_MAX_x_max_NODE;
02021 (*n_res)++;
02022 }
02023
02024 if (edge_residual(node,Z_MAX_y_min_DIR)>epsilon){
02025
02026 result[*n_res].pos.x() = node.pos.x();
02027 result[*n_res].pos.y() = node.pos.y();
02028 result[*n_res].pos.z() = node.pos.z();
02029 result[*n_res].loc = Z_MAX_y_min_NODE;
02030 (*n_res)++;
02031 }
02032
02033 if (edge_residual(node,Z_MAX_y_max_DIR)>epsilon){
02034
02035 result[*n_res].pos.x() = node.pos.x();
02036 result[*n_res].pos.y() = node.pos.y();
02037 result[*n_res].pos.z() = node.pos.z();
02038 result[*n_res].loc = Z_MAX_y_max_NODE;
02039 (*n_res)++;
02040 }
02041 }
02042 break;
02043
02044 case Z_MAX_x_min_NODE:
02045 if ((node.pos.z()<z_size)
02046 &&(domain[node.pos])
02047 &&(edge_residual(node,Z_MAX_x_max_DIR)>epsilon)){
02048
02049 result[*n_res].pos.x() = node.pos.x();
02050 result[*n_res].pos.y() = node.pos.y();
02051 result[*n_res].pos.z() = node.pos.z()+1;
02052 result[*n_res].loc = CENTER_NODE;
02053 (*n_res)++;
02054 }
02055 break;
02056
02057 case Z_MAX_x_max_NODE:
02058 if ((node.pos.z()<z_size)
02059 &&(domain[node.pos])
02060 &&(edge_residual(node,Z_MAX_x_min_DIR)>epsilon)){
02061
02062 result[*n_res].pos.x() = node.pos.x();
02063 result[*n_res].pos.y() = node.pos.y();
02064 result[*n_res].pos.z() = node.pos.z()+1;
02065 result[*n_res].loc = CENTER_NODE;
02066 (*n_res)++;
02067 }
02068 break;
02069
02070 case Z_MAX_y_min_NODE:
02071 if ((node.pos.z()<z_size)
02072 &&(domain[node.pos])
02073 &&(edge_residual(node,Z_MAX_y_max_DIR)>epsilon)){
02074
02075 result[*n_res].pos.x() = node.pos.x();
02076 result[*n_res].pos.y() = node.pos.y();
02077 result[*n_res].pos.z() = node.pos.z()+1;
02078 result[*n_res].loc = CENTER_NODE;
02079 (*n_res)++;
02080 }
02081 break;
02082
02083 case Z_MAX_y_max_NODE:
02084 if ((node.pos.z()<z_size)
02085 &&(domain[node.pos])
02086 &&(edge_residual(node,Z_MAX_y_min_DIR)>epsilon)){
02087
02088 result[*n_res].pos.x() = node.pos.x();
02089 result[*n_res].pos.y() = node.pos.y();
02090 result[*n_res].pos.z() = node.pos.z()+1;
02091 result[*n_res].loc = CENTER_NODE;
02092 (*n_res)++;
02093 }
02094 break;
02095 }
02096 }
02097
02098
02099 template<typename Real>
02100 typename PushRelabel_non_linear_grid_3D<Real>::real_type
02101 PushRelabel_non_linear_grid_3D<Real>::
02102 node_excess(const node_type& node) const{
02103
02104
02105 if ((node.loc == CENTER_NODE) && (node.pos.z() == min_bound(node.pos.x(),
02106 node.pos.y()))){
02107 #ifndef WITHOUT_LIMITS
02108 return (std::numeric_limits<real_type>::max());
02109 #else
02110 return 1e20;
02111 #endif
02112 }
02113
02114
02115 if ((node.loc == CENTER_NODE) && (node.pos.z() == max_bound(node.pos.x(),
02116 node.pos.y()))){
02117 return 0;
02118 }
02119
02120 node_type neighbor[12];
02121 size_type n_neighbor;
02122
02123 neighbors(node,neighbor,&n_neighbor);
02124
02125 real_type res = 0;
02126 for(size_type n=0;n<n_neighbor;n++){
02127 res -= edge_flow(node,direction(node,neighbor[n]));
02128 }
02129 return res;
02130 }
02131
02132
02136 template<typename Real>
02137 void
02138 PushRelabel_non_linear_grid_3D<Real>::
02139 add_flow(const node_type& from,
02140 const direction_type dir,
02141 const real_type flow){
02142
02143
02144 if (edge_residual(from,dir)+epsilon<flow){
02145 Message::error<<"not enough capacity ("<<edge_residual(from,dir)
02146 <<") for flow ("<<flow<<") ["
02147 <<(flow-edge_residual(from,dir))<<"]"
02148 <<Message::done;
02149 }
02150
02151 switch(from.loc){
02152
02153 case CENTER_NODE:
02154
02155 switch(dir){
02156 case X_MIN_DIR:
02157 x_major_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z())
02158 -= flow;
02159 return;
02160
02161 case X_MAX_DIR:
02162 x_major_derivative_edge_flow[from.pos]
02163 += flow;
02164 return;
02165
02166 case Y_MIN_DIR:
02167 y_major_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z())
02168 -= flow;
02169 return;
02170
02171 case Y_MAX_DIR:
02172 y_major_derivative_edge_flow[from.pos]
02173 += flow;
02174 return;
02175
02176 case Z_MIN_x_min_DIR:
02177 cost_edge_Z_MAX_MAX_x_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1)
02178 -= flow;
02179 return;
02180
02181 case Z_MIN_x_max_DIR:
02182 cost_edge_Z_MAX_MAX_x_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1)
02183 -= flow;
02184 return;
02185
02186 case Z_MIN_y_min_DIR:
02187 cost_edge_Z_MAX_MAX_y_min_flow(from.pos.x(),from.pos.y(),from.pos.z()-1)
02188 -= flow;
02189 return;
02190
02191 case Z_MIN_y_max_DIR:
02192 cost_edge_Z_MAX_MAX_y_max_flow(from.pos.x(),from.pos.y(),from.pos.z()-1)
02193 -= flow;
02194 return;
02195
02196 case Z_MAX_x_min_DIR:
02197 cost_edge_Z_MAX_x_min_flow[from.pos]
02198 += flow;
02199 return;
02200
02201 case Z_MAX_x_max_DIR:
02202 cost_edge_Z_MAX_x_max_flow[from.pos]
02203 += flow;
02204 return;
02205
02206 case Z_MAX_y_min_DIR:
02207 cost_edge_Z_MAX_y_min_flow[from.pos]
02208 += flow;
02209 return;
02210
02211 case Z_MAX_y_max_DIR:
02212 cost_edge_Z_MAX_y_max_flow[from.pos]
02213 += flow;
02214 return;
02215 }
02216
02217 Message::error<<"add_flow: bad end 1"<<Message::done;
02218
02219
02220 case Z_MAX_x_min_NODE:
02221
02222 switch(dir){
02223 case X_MIN_DIR:
02224 x_minor_derivative_edge_flow(from.pos.x()-1,from.pos.y(),from.pos.z())
02225 -= flow;
02226 return;
02227
02228 case X_MAX_DIR:
02229 Message::error<<"add_flow : bad end 2.1"<<Message::done;
02230
02231 case Y_MIN_DIR:
02232 Message::error<<"add_flow : bad end 2.2"<<Message::done;
02233
02234 case Y_MAX_DIR:
02235 Message::error<<"add_flow : bad end 2.3"<<Message::done;
02236
02237 case Z_MIN_x_min_DIR:
02238 Message::error<<"add_flow : bad end 2.4"<<Message::done;
02239
02240 case Z_MIN_x_max_DIR:
02241 cost_edge_Z_MAX_x_min_flow[from.pos]
02242 -= flow;
02243 return;
02244
02245 case Z_MIN_y_min_DIR:
02246 Message::error<<"add_flow : bad end 2.5"<<Message::done;
02247
02248 case Z_MIN_y_max_DIR:
02249 Message::error<<"add_flow : bad end 2.6"<<Message::done;
02250
02251 case Z_MAX_x_min_DIR:
02252 Message::error<<"add_flow : bad end 2.7"<<Message::done;
02253
02254 case Z_MAX_x_max_DIR:
02255 cost_edge_Z_MAX_MAX_x_min_flow[from.pos]
02256 += flow;
02257 return;
02258
02259 case Z_MAX_y_min_DIR:
02260 Message::error<<"add_flow : bad end 2.8"<<Message::done;
02261
02262 case Z_MAX_y_max_DIR:
02263 Message::error<<"add_flow : bad end 2.9"<<Message::done;
02264 }
02265
02266 Message::error<<"add_flow : bad end 2"<<Message::done;
02267
02268
02269 case Z_MAX_x_max_NODE:
02270
02271 switch(dir){
02272 case X_MIN_DIR:
02273 Message::error<<"add_flow : bad end 3.1"<<Message::done;
02274
02275 case X_MAX_DIR:
02276 x_minor_derivative_edge_flow[from.pos]
02277 += flow;
02278 return;
02279
02280 case Y_MIN_DIR:
02281 Message::error<<"add_flow : bad end 3.2"<<Message::done;
02282
02283 case Y_MAX_DIR:
02284 Message::error<<"add_flow : bad end 3.3"<<Message::done;
02285
02286 case Z_MIN_x_min_DIR:
02287 cost_edge_Z_MAX_x_max_flow[from.pos]
02288 -= flow;
02289 return;
02290
02291 case Z_MIN_x_max_DIR:
02292 Message::error<<"add_flow : bad end 3.4"<<Message::done;
02293
02294 case Z_MIN_y_min_DIR:
02295 Message::error<<"add_flow : bad end 3.5"<<Message::done;
02296
02297 case Z_MIN_y_max_DIR:
02298 Message::error<<"add_flow : bad end 3.6"<<Message::done;
02299
02300 case Z_MAX_x_min_DIR:
02301 cost_edge_Z_MAX_MAX_x_max_flow[from.pos]
02302 += flow;
02303 return;
02304
02305 case Z_MAX_x_max_DIR:
02306 Message::error<<"add_flow : bad end 3.7"<<Message::done;
02307
02308 case Z_MAX_y_min_DIR:
02309 Message::error<<"add_flow : bad end 3.8"<<Message::done;
02310
02311 case Z_MAX_y_max_DIR:
02312 Message::error<<"add_flow : bad end 3.9"<<Message::done;
02313 }
02314
02315 Message::error<<"add_flow : bad end 3"<<Message::done;
02316
02317
02318 case Z_MAX_y_min_NODE:
02319
02320 switch(dir){
02321 case X_MIN_DIR:
02322 Message::error<<"add_flow : bad end 4.1"<<Message::done;
02323
02324 case X_MAX_DIR:
02325 Message::error<<"add_flow : bad end 4.2"<<Message::done;
02326
02327 case Y_MIN_DIR:
02328 y_minor_derivative_edge_flow(from.pos.x(),from.pos.y()-1,from.pos.z())
02329 -= flow;
02330 return;
02331
02332 case Y_MAX_DIR:
02333 Message::error<<"add_flow : bad end 4.3"<<Message::done;
02334
02335 case Z_MIN_x_min_DIR:
02336 Message::error<<"add_flow : bad end 4.4"<<Message::done;
02337
02338 case Z_MIN_x_max_DIR:
02339 Message::error<<"add_flow : bad end 4.5"<<Message::done;
02340
02341 case Z_MIN_y_min_DIR:
02342 Message::error<<"add_flow : bad end 4.6"<<Message::done;
02343
02344 case Z_MIN_y_max_DIR:
02345 cost_edge_Z_MAX_y_min_flow[from.pos]
02346 -= flow;
02347 return;
02348
02349 case Z_MAX_x_min_DIR:
02350 Message::error<<"add_flow : bad end 4.7"<<Message::done;
02351
02352 case Z_MAX_x_max_DIR:
02353 Message::error<<"add_flow : bad end 4.8"<<Message::done;
02354
02355 case Z_MAX_y_min_DIR:
02356 Message::error<<"add_flow : bad end 4.9"<<Message::done;
02357
02358 case Z_MAX_y_max_DIR:
02359 cost_edge_Z_MAX_MAX_y_min_flow[from.pos]
02360 += flow;
02361 return;
02362 }
02363
02364 Message::error<<"add_flow : bad end 4"<<Message::done;
02365
02366
02367 case Z_MAX_y_max_NODE:
02368
02369 switch(dir){
02370 case X_MIN_DIR:
02371 Message::error<<"add_flow : bad end 5.1"<<Message::done;
02372
02373 case X_MAX_DIR:
02374 Message::error<<"add_flow : bad end 5.2"<<Message::done;
02375
02376 case Y_MIN_DIR:
02377 Message::error<<"add_flow : bad end 5.3"<<Message::done;
02378
02379 case Y_MAX_DIR:
02380 y_minor_derivative_edge_flow[from.pos]
02381 += flow;
02382 return;
02383
02384 case Z_MIN_x_min_DIR:
02385 Message::error<<"add_flow : bad end 5.4"<<Message::done;
02386
02387 case Z_MIN_x_max_DIR:
02388 Message::error<<"add_flow : bad end 5.5"<<Message::done;
02389
02390 case Z_MIN_y_min_DIR:
02391 cost_edge_Z_MAX_y_max_flow[from.pos]
02392 -= flow;
02393 return;
02394
02395 case Z_MIN_y_max_DIR:
02396 Message::error<<"add_flow : bad end 5.6"<<Message::done;
02397
02398 case Z_MAX_x_min_DIR:
02399 Message::error<<"add_flow : bad end 5.7"<<Message::done;
02400
02401 case Z_MAX_x_max_DIR:
02402 Message::error<<"add_flow : bad end 5.8"<<Message::done;
02403
02404 case Z_MAX_y_min_DIR:
02405 cost_edge_Z_MAX_MAX_y_max_flow[from.pos]
02406 += flow;
02407 return;
02408
02409 case Z_MAX_y_max_DIR:
02410 Message::error<<"add_flow : bad end 5.9"<<Message::done;
02411 }
02412
02413 Message::error<<"add_flow : bad end 5"<<Message::done;
02414 }
02415
02416 Message::error<<"add_flow : bad end 0"<<Message::done;
02417 }
02418
02419
02420 template<typename Real>
02421 typename PushRelabel_non_linear_grid_3D<Real>::direction_type
02422 PushRelabel_non_linear_grid_3D<Real>::
02423 direction(const node_type& from,
02424 const node_type& to) const{
02425
02426 switch(from.loc){
02427
02428
02429 case CENTER_NODE:
02430 switch(to.loc){
02431
02432 case CENTER_NODE:
02433 if (from.pos.z() == to.pos.z()){
02434
02435 if ((from.pos.x()+1 == to.pos.x()) && (from.pos.y() == to.pos.y())){
02436 return X_MAX_DIR;
02437 }
02438 else if ((from.pos.x() == to.pos.x()+1) && (from.pos.y() == to.pos.y())){
02439 return X_MIN_DIR;
02440 }
02441 else if ((from.pos.x() == to.pos.x()) && (from.pos.y()+1 == to.pos.y())){
02442 return Y_MAX_DIR;
02443 }
02444 else if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y()+1)){
02445 return Y_MIN_DIR;
02446 }
02447 else{
02448 Message::error<<"direction : bad end 1.1"<<Message::done;
02449 return 0;
02450 }
02451 }
02452 else{
02453 Message::error<<"direction : bad end 1.2"<<Message::done;
02454 return 0;
02455 }
02456
02457 case Z_MAX_x_min_NODE:
02458 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02459 if (to.pos.z()+1 == from.pos.z()){
02460 return Z_MIN_x_min_DIR;
02461 }
02462 else if (from.pos.z() == to.pos.z()){
02463 return Z_MAX_x_min_DIR;
02464 }
02465 else{
02466 Message::error<<"direction : bad end 1.3"<<Message::done;
02467 return 0;
02468 }
02469 }
02470 else{
02471 Message::error<<"direction : bad end 1.4"<<Message::done;
02472 return 0;
02473 }
02474
02475 case Z_MAX_x_max_NODE:
02476 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02477 if (to.pos.z()+1 == from.pos.z()){
02478 return Z_MIN_x_max_DIR;
02479 }
02480 else if (from.pos.z() == to.pos.z()){
02481 return Z_MAX_x_max_DIR;
02482 }
02483 else{
02484 Message::error<<"direction : bad end 1.5"<<Message::done;
02485 return 0;
02486 }
02487 }
02488 else{
02489 Message::error<<"direction : bad end 1.6"<<Message::done;
02490 return 0;
02491 }
02492
02493 case Z_MAX_y_min_NODE:
02494 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02495 if (to.pos.z()+1 == from.pos.z()){
02496 return Z_MIN_y_min_DIR;
02497 }
02498 else if (from.pos.z() == to.pos.z()){
02499 return Z_MAX_y_min_DIR;
02500 }
02501 else{
02502 Message::error<<"direction : bad end 1.7"<<Message::done;
02503 return 0;
02504 }
02505 }
02506 else{
02507 Message::error<<"direction : bad end 1.8\n"
02508 <<VALUE_OF(from.pos)<<"\t"<<VALUE_OF(int(from.loc))<<"\n"
02509 <<VALUE_OF(to.pos)<<"\t"<<VALUE_OF(int(to.loc))
02510 <<Message::done;
02511 return 0;
02512 }
02513
02514 case Z_MAX_y_max_NODE:
02515 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02516 if (to.pos.z()+1 == from.pos.z()){
02517 return Z_MIN_y_max_DIR;
02518 }
02519 else if (from.pos.z() == to.pos.z()){
02520 return Z_MAX_y_max_DIR;
02521 }
02522 else{
02523 Message::error<<"direction : bad end 1.9"<<Message::done;
02524 return 0;
02525 }
02526 }
02527 else{
02528 Message::error<<"direction : bad end 1.10"<<Message::done;
02529 return 0;
02530 }
02531 }
02532
02533
02534 case Z_MAX_x_min_NODE:
02535 switch(to.loc){
02536
02537 case CENTER_NODE:
02538 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02539 if (from.pos.z() == to.pos.z()){
02540 return Z_MIN_x_max_DIR;
02541 }
02542 else if (from.pos.z()+1 == to.pos.z()){
02543 return Z_MAX_x_max_DIR;
02544 }
02545 else{
02546 Message::error<<"direction : bad end 2.1"<<Message::done;
02547 return 0;
02548 }
02549 }
02550
02551 case Z_MAX_x_min_NODE:
02552 Message::error<<"direction : bad end 2.2"<<Message::done;
02553 return 0;
02554
02555
02556 case Z_MAX_x_max_NODE:
02557 if ((from.pos.y() == to.pos.y())&&
02558 (from.pos.z() == to.pos.z())&&
02559 (from.pos.x() == to.pos.x()+1)){
02560 return X_MIN_DIR;
02561 }
02562 else{
02563 Message::error<<"direction : bad end 2.3 "<<Message::done;
02564 return 0;
02565 }
02566
02567 case Z_MAX_y_min_NODE:
02568 Message::error<<"direction : bad end 2.4"<<Message::done;
02569 return 0;
02570
02571
02572 case Z_MAX_y_max_NODE:
02573 Message::error<<"direction : bad end 2.5 "<<Message::done;
02574 return 0;
02575 }
02576
02577
02578 case Z_MAX_x_max_NODE:
02579 switch(to.loc){
02580
02581 case CENTER_NODE:
02582 if ((from.pos.x() == to.pos.x())&&(from.pos.y() == to.pos.y())){
02583 if (from.pos.z() == to.pos.z()){
02584 return Z_MIN_x_min_DIR;
02585 }
02586 else if (from.pos.z()+1 == to.pos.z()){
02587 return Z_MAX_x_min_DIR;
02588 }
02589 else{
02590 Message::error<<"direction : bad end 3.1"<<Message::done;
02591 return 0;
02592 }
02593 }
02594
02595 case Z_MAX_x_min_NODE:
02596 if ((from.pos.y() == to.pos.y())&&
02597 (from.pos.z() == to.pos.z())&&
02598 (from.pos.x()+1 == to.pos.x())){
02599 return X_MAX_DIR;
02600 }
02601
02602 case Z_MAX_x_max_NODE:
02603 Message::error<<"direction : bad end 3.2"<<Message::done;
02604 return 0;
02605
02606 case Z_MAX_y_min_NODE:
02607 Message::error<<"direction : bad end 3.3"<<Message::done;
02608 return 0;
02609
02610 case Z_MAX_y_max_NODE:
02611 Message::error<<"direction : bad end 3.4"<<Message::done;
02612 return 0;
02613 }
02614
02615
02616 case Z_MAX_y_min_NODE:
02617 switch(to.loc){
02618
02619 case CENTER_NODE:
02620 if ((from.pos.x() == to.pos.x()) && (from.pos.y() == to.pos.y())){
02621 if (from.pos.z() == to.pos.z()){
02622 return Z_MIN_y_max_DIR;
02623 }
02624 else if (from.pos.z()+1 == to.pos.z()){
02625 return Z_MAX_y_max_DIR;
02626 }
02627 else{
02628 Message::error<<"direction : bad end 4.1"<<Message::done;
02629 return 0;
02630 }
02631 }
02632
02633 case Z_MAX_x_min_NODE:
02634 Message::error<<"direction : bad end 4.2"<<Message::done;
02635 return 0;
02636
02637
02638 case Z_MAX_x_max_NODE:
02639 Message::error<<"direction : bad end 4.3 "<<Message::done;
02640 return 0;
02641
02642 case Z_MAX_y_min_NODE:
02643 Message::error<<"direction : bad end 4.4"<<Message::done;
02644 return 0;
02645
02646
02647 case Z_MAX_y_max_NODE:
02648 if ((from.pos.x() == to.pos.x())&&
02649 (from.pos.z() == to.pos.z())&&
02650 (from.pos.y() == to.pos.y()+1)){
02651 return Y_MIN_DIR;
02652 }
02653 else{
02654 Message::error<<"direction : bad end 4.5 "<<Message::done;
02655 return 0;
02656 }
02657 }
02658
02659
02660 case Z_MAX_y_max_NODE:
02661 switch(to.loc){
02662
02663 case CENTER_NODE:
02664 if ((from.pos.x() == to.pos.x())&&(from.pos.y() == to.pos.y())){
02665 if (from.pos.z() == to.pos.z()){
02666 return Z_MIN_y_min_DIR;
02667 }
02668 else if (from.pos.z()+1 == to.pos.z()){
02669 return Z_MAX_y_min_DIR;
02670 }
02671 else{
02672 Message::error<<"direction : bad end 5.1"<<Message::done;
02673 return 0;
02674 }
02675 }
02676
02677 case Z_MAX_x_min_NODE:
02678 Message::error<<"direction : bad end 5.2"<<Message::done;
02679 return 0;
02680
02681 case Z_MAX_x_max_NODE:
02682 Message::error<<"direction : bad end 5.3"<<Message::done;
02683 return 0;
02684
02685 case Z_MAX_y_min_NODE:
02686 if ((from.pos.x() == to.pos.x())&&
02687 (from.pos.z() == to.pos.z())&&
02688 (from.pos.y()+1 == to.pos.y())){
02689 return Y_MAX_DIR;
02690 }
02691
02692 case Z_MAX_y_max_NODE:
02693 Message::error<<"direction : bad end 5.4"<<Message::done;
02694 return 0;
02695 }
02696
02697 }
02698
02699 Message::error<<"direction : bad end 0"<<Message::done;
02700 return 0;
02701 }
02702
02703
02704 template<typename Real>
02705 typename PushRelabel_non_linear_grid_3D<Real>::label_type
02706 PushRelabel_non_linear_grid_3D<Real>::
02707 relabel(const node_type& node){
02708
02709 node_type neighbor[12];
02710 size_type n_neighbor;
02711 #ifndef WITHOUT_LIMITS
02712 label_type min = std::numeric_limits<label_type>::max();
02713 #else
02714 label_type min = static_cast<size_type>(4294967295);
02715 #endif
02716
02717 static unsigned int NB_RELABEL = 0;
02718 static unsigned int NB_GLOBAL = 0;
02719
02720 if (NB_RELABEL >= nb_nodes/2 ){
02721 global_relabel();
02722 NB_GLOBAL++;
02723 NB_RELABEL = 0;
02724 return 0;
02725 }
02726
02727
02728 NB_RELABEL++;
02729
02730
02731 if ((node.pos.z() == min_bound(node.pos.x(),node.pos.y())) &&
02732 (node.loc == CENTER_NODE)){
02733
02734 label_CENTER[node.pos] = nb_nodes;
02735
02736 return nb_nodes;
02737 }
02738
02739 adjacent_nodes(node,neighbor,&n_neighbor);
02740
02741 for(size_type n=0;n<n_neighbor;n++){
02742 const label_type l = label(neighbor[n]);
02743
02744 if (l<min){
02745 min = l;
02746 }
02747 }
02748
02749 #ifndef WITHOUT_LIMITS
02750 if (min == std::numeric_limits<label_type>::max()){
02751 #else
02752 if (min == static_cast<size_type>(4294967295)){
02753 #endif
02754 Message::error<<"#neighbor = "<<n_neighbor<<std::endl;
02755 Message::error<<"labels :"<<std::endl;
02756 for(size_type n=0;n<n_neighbor;n++){
02757 const label_type l = label(neighbor[n]);
02758
02759 Message::error<<l<<std::endl;
02760 }
02761
02762 Message::error<<"relabel : no neighbor reachable !"<<Message::done;
02763 }
02764
02765 min++;
02766
02767 label(node) = min;
02768
02769 return min;
02770 }
02771
02772
02773 namespace{
02774 template<typename Real>
02775 struct Index_cmp{
02776 typedef typename PushRelabel_non_linear_grid_3D<Real>::index_type
02777 index_type;
02778
02779 typedef typename PushRelabel_non_linear_grid_3D<Real>::node_type
02780 node_type;
02781
02782 const PushRelabel_non_linear_grid_3D<Real>& graph;
02783
02784 Index_cmp(const PushRelabel_non_linear_grid_3D<Real>& g):graph(g){}
02785
02786 inline bool operator()(const index_type i1,
02787 const index_type i2){
02788 node_type n1,n2;
02789
02790 graph.from_index(&n1,i1);
02791 graph.from_index(&n2,i2);
02792
02793 return (graph.label(n1)>graph.label(n2));
02794 }
02795 };
02796 }
02797
02798
02799 template<typename Real>
02800 void
02801 PushRelabel_non_linear_grid_3D<Real>::
02802 global_relabel(){
02803
02804 using namespace std;
02805
02806 cout<<"GLOBAL RELABEL ";
02807
02808 deque<node_type>* to_explore = new deque<node_type>;
02809 deque<node_type>* exploring = new deque<node_type>;
02810 typedef typename deque<node_type>::iterator iterator;
02811
02812 for(size_type x=0;x<x_size;x++){
02813 for(size_type y=0;y<y_size;y++){
02814 for(size_type z=0;z<z_size+1;z++){
02815
02816 label_CENTER(x,y,z) = 0;
02817 label_Z_MAX_x_min(x,y,z) = 0;
02818 label_Z_MAX_x_max(x,y,z) = 0;
02819 label_Z_MAX_y_min(x,y,z) = 0;
02820 label_Z_MAX_y_max(x,y,z) = 0;
02821 }
02822 }
02823 }
02824
02825 for(size_type x=0;x<x_size;x++){
02826 for(size_type y=0;y<y_size;y++){
02827 const size_type z = max_bound(x,y);
02828 node_type node(position_type(x,y,z),CENTER_NODE);
02829
02830 if ((domain(x,y,z-1)) && (label_CENTER(x,y,z) < nb_nodes)){
02831
02832 exploring->push_back(node);
02833 label_CENTER(x,y,z) = 1;
02834 }
02835 }
02836 }
02837
02838 while(!exploring->empty()){
02839 for(iterator node=exploring->begin();node!=exploring->end();node++){
02840
02841 node_type neighbor[12];
02842 size_type n_neighbor;
02843
02844 neighbors(*node,neighbor,&n_neighbor);
02845
02846 const label_type l = label(*node);
02847
02848 for(size_type n=0;n<n_neighbor;n++){
02849 const node_type n_node = neighbor[n];
02850 const direction_type dir = direction(n_node,*node);
02851
02852 if (edge_residual(n_node,dir)>epsilon){
02853
02854 const label_type ln = label(n_node);
02855
02856 if (((ln > l+1)||(ln == 0))&&(ln < nb_nodes)){
02857
02858 label(n_node) = l+1;
02859
02860 to_explore->push_back(n_node);
02861
02862 }
02863 }
02864
02865 }
02866 }
02867
02868 delete exploring;
02869 exploring = to_explore;
02870 to_explore = new deque<node_type>;
02871 }
02872
02873 delete exploring;
02874 delete to_explore;
02875
02876 label_type max_l = 0;
02877 label_type max_L = 0;
02878 active_node_index.clear();
02879
02880 for(size_type x=0;x<x_size;x++){
02881 for(size_type y=0;y<y_size;y++){
02882 for(size_type z=0;z<z_size+1;z++){
02883
02884 if (label_CENTER(x,y,z) == 0) {
02885 label_CENTER(x,y,z) = nb_nodes;
02886 }
02887
02888 if (label_Z_MAX_x_min(x,y,z) == 0) {
02889 label_Z_MAX_x_min(x,y,z) = nb_nodes;
02890 }
02891
02892 if (label_Z_MAX_x_max(x,y,z) == 0) {
02893 label_Z_MAX_x_max(x,y,z) = nb_nodes;
02894 }
02895
02896 if (label_Z_MAX_y_min(x,y,z) == 0) {
02897 label_Z_MAX_y_min(x,y,z) = nb_nodes;
02898 }
02899
02900 if (label_Z_MAX_y_max(x,y,z) == 0) {
02901 label_Z_MAX_y_max(x,y,z) = nb_nodes;
02902 }
02903
02905 label_type l = label_CENTER(x,y,z);
02906 node_type n = node_type(position_type(x,y,z),CENTER_NODE);
02907 if ((l>max_L)&&(l<nb_nodes)){
02908 max_L = l;
02909 }
02910 if ((l<nb_nodes)&&(node_excess(n)>epsilon)){
02911 if (l>max_l){
02912 max_l = l;
02913 }
02914
02915 index_type i;
02916 to_index(n,&i);
02917 active_node_index.push_back(i);
02918 }
02919
02920 l = label_Z_MAX_x_min(x,y,z);
02921 n = node_type(position_type(x,y,z),Z_MAX_x_min_NODE);
02922 if ((l>max_L)&&(l<nb_nodes)){
02923 max_L = l;
02924 }
02925 if ((l<nb_nodes)&&(node_excess(n)>epsilon)){
02926 if (l>max_l){
02927 max_l = l;
02928 }
02929
02930 index_type i;
02931 to_index(n,&i);
02932 active_node_index.push_back(i);
02933 }
02934
02935 l = label_Z_MAX_x_max(x,y,z);
02936 n = node_type(position_type(x,y,z),Z_MAX_x_max_NODE);
02937 if ((l>max_L)&&(l<nb_nodes)){
02938 max_L = l;
02939 }
02940 if ((l<nb_nodes)&&(node_excess(n)>epsilon)){
02941 if (l>max_l){
02942 max_l = l;
02943 }
02944
02945 index_type i;
02946 to_index(n,&i);
02947 active_node_index.push_back(i);
02948 }
02949
02950 l = label_Z_MAX_y_min(x,y,z);
02951 n = node_type(position_type(x,y,z),Z_MAX_y_min_NODE);
02952 if ((l>max_L)&&(l<nb_nodes)){
02953 max_L = l;
02954 }
02955 if ((l<nb_nodes)&&(node_excess(n)>epsilon)){
02956 if (l>max_l){
02957 max_l = l;
02958 }
02959
02960 index_type i;
02961 to_index(n,&i);
02962 active_node_index.push_back(i);
02963 }
02964
02965 l = label_Z_MAX_y_max(x,y,z);
02966 n = node_type(position_type(x,y,z),Z_MAX_y_max_NODE);
02967 if ((l>max_L)&&(l<nb_nodes)){
02968 max_L = l;
02969 }
02970 if ((l<nb_nodes)&&(node_excess(n)>epsilon)){
02971 if (l>max_l){
02972 max_l = l;
02973 }
02974
02975 index_type i;
02976 to_index(n,&i);
02977 active_node_index.push_back(i);
02978 }
02980 }
02981 }
02982 }
02983
02984 Index_cmp<real_type> cmp(*this);
02985
02986 std::sort(active_node_index.begin(),active_node_index.end(),cmp);
02987
02988
02989 cout<<"label max : "<<max_l<<" "<<max_L<<endl;
02990
02991 }
02992
02993
02997 template<typename Real>
02998 bool
02999 PushRelabel_non_linear_grid_3D<Real>::
03000 discharge(const node_type& node){
03001
03002 using namespace std;
03003
03004 node_type neighbor[12];
03005 size_type n_neighbor;
03006
03007 adjacent_nodes(node,neighbor,&n_neighbor);
03008
03009 real_type excess = node_excess(node);
03010 size_type n = 0;
03011
03012 if (excess <= epsilon) {
03013 return false;
03014 }
03015
03016
03017 while ((n<n_neighbor)&&(excess>epsilon)){
03018
03019 const node_type n_node = neighbor[n];
03020
03021 if (is_admissible(node,n_node)){
03022
03023 const direction_type d = direction(node,n_node);
03024 const real_type to_push = std::min(excess,edge_residual(node,d));
03025
03026 add_flow(node,d,to_push);
03027
03028 excess -= to_push;
03029
03030 add_active_node(n_node);
03031 }
03032 n++;
03033 }
03034
03035 if ((n==n_neighbor)&&(excess>epsilon)){
03036 relabel(node);
03037 add_active_node(node);
03038 }
03039
03040 return true;
03041 }
03042
03043
03044
03045 template<typename Real>
03046 void
03047 PushRelabel_non_linear_grid_3D<Real>::
03048 add_active_node(const node_type& node){
03049
03050 index_type i;
03051
03052 to_index(node,&i);
03053
03054 active_node_index.push_back(i);
03055 }
03056
03057
03058
03059 template<typename Real>
03060 void
03061 PushRelabel_non_linear_grid_3D<Real>::
03062 get_next_node(node_type* const node,
03063 bool* const no_more_node){
03064
03065 label_type l;
03066
03067 do{
03068 *no_more_node = active_node_index.empty();
03069
03070 if(!(*no_more_node)){
03071
03072 const index_type i = active_node_index.front();
03073 active_node_index.pop_front();
03074
03075 from_index(node,i);
03076
03077 l = label(*node);
03078 }
03079
03080 } while(((l >= nb_nodes)||
03081 (node->pos.z() == max_bound(node->pos.x(),node->pos.y()))) &&
03082 (!(*no_more_node)));
03083 }
03084
03085
03086
03087 template<typename Real>
03088 bool
03089 PushRelabel_non_linear_grid_3D<Real>::
03090 is_active(const node_type& node) const{
03091
03092 return (node_excess(node) > epsilon);
03093 }
03094
03095
03096
03097 template<typename Real>
03098 bool
03099 PushRelabel_non_linear_grid_3D<Real>::
03100 is_admissible(const node_type& from,
03101 const node_type& to) const{
03102
03103 return (label(from) == (label(to)+1));
03104 }
03105
03106
03107 template<typename Real>
03108 typename PushRelabel_non_linear_grid_3D<Real>::real_type
03109 PushRelabel_non_linear_grid_3D<Real>::
03110 max_flow(){
03111
03112 using namespace std;
03113
03114 global_relabel();
03115
03116 node_type to_discharge;
03117 bool no_more_to_discharge;
03118
03119 get_next_node(&to_discharge,&no_more_to_discharge);
03120
03121 while(!no_more_to_discharge){
03122 discharge(to_discharge);
03123 get_next_node(&to_discharge,&no_more_to_discharge);
03124 }
03125
03126
03127 real_type flow = 0;
03128 for(size_type x=0;x<x_size;x++){
03129 for(size_type y=0;y<y_size;y++){
03130
03131 const size_type z = max_bound(x,y) - 1;
03132
03133 if (domain(x,y,z)){
03134 flow += cost_edge_Z_MAX_x_min_flow(x,y,z)
03135 + cost_edge_Z_MAX_x_max_flow(x,y,z)
03136 + cost_edge_Z_MAX_y_min_flow(x,y,z)
03137 + cost_edge_Z_MAX_y_max_flow(x,y,z);
03138 }
03139 }
03140 }
03141
03142 return flow;
03143 }
03144
03145
03146
03147 template<typename Real>
03148 void
03149 PushRelabel_non_linear_grid_3D<Real>::
03150 min_cut(Array_2D<size_type>* const cut) const{
03151
03152 using namespace std;
03153
03154 Array_3D<bool> visited_CENTER(x_size,y_size,z_size+1,false);
03155 Array_3D<bool> visited_Z_MAX_x_min(x_size,y_size,z_size+1,false);
03156 Array_3D<bool> visited_Z_MAX_x_max(x_size,y_size,z_size+1,false);
03157 Array_3D<bool> visited_Z_MAX_y_min(x_size,y_size,z_size+1,false);
03158 Array_3D<bool> visited_Z_MAX_y_max(x_size,y_size,z_size+1,false);
03159
03160 list<node_type> to_visit;
03161
03162
03163 for(size_type x=0;x<x_size;x++){
03164 for(size_type y=0;y<y_size;y++){
03165 node_type p(position_type(x,y,max_bound(x,y)),CENTER_NODE);
03166 if (domain(p.pos.x(),p.pos.y(),p.pos.z()-1)&&
03167 (label_CENTER[p.pos] < nb_nodes)){
03168 to_visit.push_back(p);
03169 visited_CENTER[p.pos] = true;
03170 }
03171 }
03172 }
03173
03174
03175 node_type neighbor[12];
03176 size_type n_neighbor;
03177
03178 while(!to_visit.empty()){
03179 node_type node = to_visit.front();
03180 to_visit.pop_front();
03181
03182 neighbors(node,neighbor,&n_neighbor);
03183
03184 for(size_type n=0;n<n_neighbor;n++){
03185
03186 const node_type n_node = neighbor[n];
03187 const direction_type dir = direction(n_node,node);
03188
03189 if ((edge_residual(n_node,dir)>epsilon)){
03190
03191 switch(n_node.loc){
03192 case CENTER_NODE:
03193
03194 if(!visited_CENTER[neighbor[n].pos]){
03195 visited_CENTER[neighbor[n].pos] = true;
03196 to_visit.push_back(neighbor[n]);
03197 }
03198 break;
03199
03200 case Z_MAX_x_min_NODE:
03201
03202 if(!visited_Z_MAX_x_min[neighbor[n].pos]){
03203 visited_Z_MAX_x_min[neighbor[n].pos] = true;
03204 to_visit.push_back(neighbor[n]);
03205 }
03206 break;
03207
03208 case Z_MAX_x_max_NODE:
03209
03210 if(!visited_Z_MAX_x_max[neighbor[n].pos]){
03211 visited_Z_MAX_x_max[neighbor[n].pos] = true;
03212 to_visit.push_back(neighbor[n]);
03213 }
03214 break;
03215
03216 case Z_MAX_y_min_NODE:
03217
03218 if(!visited_Z_MAX_y_min[neighbor[n].pos]){
03219 visited_Z_MAX_y_min[neighbor[n].pos] = true;
03220 to_visit.push_back(neighbor[n]);
03221 }
03222 break;
03223
03224 case Z_MAX_y_max_NODE:
03225
03226 if(!visited_Z_MAX_y_max[neighbor[n].pos]){
03227 visited_Z_MAX_y_max[neighbor[n].pos] = true;
03228 to_visit.push_back(neighbor[n]);
03229 }
03230 break;
03231 }
03232 }
03233 }
03234 }
03235
03236
03237 cut->assign(x_size,y_size,0);
03238
03239 for(size_type x=0;x<x_size;x++){
03240 for(size_type y=0;y<y_size;y++){
03241
03242 bool found = false;
03243
03244 for(size_type z=min_bound(x,y)+1;(z<z_size);z++){
03245 if ((visited_CENTER(x,y,z))&&(!found)){
03246 found = true;
03247 (*cut)(x,y) = z-1;
03248 }
03249 }
03250 }
03251 }
03252
03253 }
03254
03255
03256 }
03257
03258 #endif
03259
03260
