amroc/hds/GridData2.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 #ifndef _included_GridData_2_h
00004 #define _included_GridData_2_h
00005 
00011 #include "DAGHParams.h"
00012 
00013 #ifdef DEBUG_PRINT_GD_MEMORY
00014 #include "DAGHMemoryTrace.h"
00015 #endif
00016 
00017 #include "BBox.h"
00018 #include "Coords.h"
00019 #include "PackedGridDataBucket.h"
00020 
00021 #ifdef DEBUG_PRINT
00022 #include "CommServer.h"
00023 #endif
00024 
00025 #include "IndexGridData2.h"
00026 #include "generic.h"
00027 
00028 #ifndef gd_OperateRegion
00029 #define gd_OperateRegion(op)    name2(gd_OperateRegion,op)
00030 #endif
00031 
00032 // #include "GDIterator.h" 
00033 
00034 #include <iosfwd>
00035 #include <iostream>
00036 #include <fstream>
00037 #include <sstream>
00038 #include <cstdlib>
00039 #include <cassert>
00040 
00041 template <class Type, int dim> class GridData;
00042 template <class Type>
00043 std::ostream& operator<<(std::ostream&, const GridData<Type,2>&);
00044 template <class Type>
00045 std::ofstream& operator<<(std::ofstream&, const GridData<Type,2>&);
00046 template <class Type>
00047 std::ifstream& operator>>(std::ifstream&, GridData<Type,2>&);
00048 template <class Type>
00049 std::stringstream& operator<<(std::stringstream&, const GridData<Type,2>&);
00050 template <class Type>
00051 std::stringstream& operator>>(std::stringstream&, GridData<Type,2>&);
00052 
00062 template <class Type> class GridData<Type,2> {
00063    BBox _bbox;
00064    Coords _extents;
00065    Coords _step;
00066    int _size;
00067    int _bottom;
00068    Type *_data;
00069 
00070 public:
00071   GridData(void) : _bbox(), _extents(0,0), _step(2,1),
00072                    _size(0), _bottom(0), _data((Type *) 0) {}
00073 
00074   GridData(const BBox &bb) : _bbox(bb) {
00075     _bbox.rank = 2; 
00076     _bbox.lower().rank = 2;
00077     _bbox.upper().rank = 2;
00078     _bbox.stepsize().rank = 2;
00079     _extents = _bbox.extents();
00080     _step = _bbox.stepsize();
00081     _bottom = _bbox.bottom();
00082     _size = _bbox.size();
00083     _data = _size ? new Type[_size] : ((Type *) 0); 
00084 #ifdef DEBUG_PRINT_GD
00085     assert (_data!=0);
00086     fill(9999);
00087     ( comm_service::log() << "Data Check"
00088       << _bbox 
00089       << ": " << _data[_size-1]
00090       << std::endl ).flush();
00091 #endif
00092 #ifdef DEBUG_PRINT_GD_MEMORY
00093     DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00094 #endif
00095   }
00096 
00097   GridData(const int i, const int j, const int ii, const int jj) : 
00098     _bbox(2,i,j,ii,jj,1), _extents(_bbox.extents()), _step(_bbox.stepsize()),
00099     _size(_bbox.size()), _bottom(_bbox.bottom()),
00100     _data(_size ? new Type[_size] : ((Type *) 0)) {
00101 #ifdef DEBUG_PRINT_GD
00102     assert (_data!=0);
00103     fill(9999);
00104     ( comm_service::log() << "Data Check"
00105       << _bbox 
00106       << ": " << _data[_size-1]
00107       << std::endl ).flush();
00108 #endif
00109 #ifdef DEBUG_PRINT_GD_MEMORY
00110     DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00111 #endif
00112   }
00113 
00114   GridData(const int i, const int j, const int ii, const int jj, const int s) : 
00115     _bbox(2,i,j,ii,jj,s), _extents(_bbox.extents()), _step(_bbox.stepsize()),
00116     _size(_bbox.size()), _bottom(_bbox.bottom()),
00117     _data(_size ? new Type[_size] : ((Type *) 0)) {
00118 #ifdef DEBUG_PRINT_GD
00119     assert (_data!=0);
00120     fill(9999);
00121     ( comm_service::log() << "Data Check"
00122       << _bbox 
00123       << ": " << _data[_size-1]
00124       << std::endl ).flush();
00125 #endif
00126 #ifdef DEBUG_PRINT_GD_MEMORY
00127     DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00128 #endif
00129    }
00130 
00131    GridData(const int i, const int j, const int ii, const int jj, 
00132             const int s, const int ss) : 
00133      _bbox(2,i,j,ii,jj,s,ss), _extents(_bbox.extents()), _step(_bbox.stepsize()),
00134      _size(_bbox.size()), _bottom(_bbox.bottom()),
00135      _data(_size ? new Type[_size] : ((Type *) 0)) {
00136 #ifdef DEBUG_PRINT_GD
00137      assert (_data!=0);
00138      fill(9999);
00139      ( comm_service::log() << "Data Check"
00140        << _bbox
00141        << ": " << _data[_size-1]
00142        << std::endl ).flush();
00143 #endif
00144 #ifdef DEBUG_PRINT_GD_MEMORY
00145      DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00146 #endif
00147    }
00148 
00149    GridData(const BBox &bb, Type *databuf) : 
00150      _bbox(bb), _data(databuf) {
00151      _bbox.rank = 2; 
00152      _bbox.lower().rank = 2;
00153      _bbox.upper().rank = 2;
00154      _bbox.stepsize().rank = 2;
00155      _extents = _bbox.extents();
00156      _step = _bbox.stepsize();
00157      _bottom = _bbox.bottom();
00158      _size = _bbox.size();
00159 #ifdef DEBUG_PRINT_GD
00160      assert (_data!=0);
00161      fill(9999);
00162      ( comm_service::log() << "Data Check"
00163        << _bbox 
00164        << ": " << _data[_size-1]
00165        << std::endl ).flush();
00166 #endif
00167    }
00168 
00169   GridData(GridDataBucket<Type> &gdbkt) :
00170     _bbox(gdbkt.bbox()), _data(gdbkt.data()) {
00171     _bbox.rank = 2; 
00172     _bbox.lower().rank = 2;
00173     _bbox.upper().rank = 2;
00174     _bbox.stepsize().rank = 2;
00175     _extents = _bbox.extents();
00176     _step = _bbox.stepsize();
00177     _bottom = _bbox.bottom();
00178     _size = _bbox.size();
00179 #ifdef DEBUG_PRINT_GD
00180     assert (_data!=0);
00181     fill(9999);
00182     ( comm_service::log() << "Data Check"
00183       << _bbox 
00184       << ": " << _data[_size-1]
00185       << std::endl ).flush();
00186 #endif
00187    }
00188 
00189    GridData(GridDataBucket<Type> &gdbkt, const int n) : 
00190      _bbox(gdbkt.bbox(n)), _data(gdbkt.data(n)) {
00191      _bbox.rank = 2; 
00192      _bbox.lower().rank = 2;
00193      _bbox.upper().rank = 2;
00194      _bbox.stepsize().rank = 2;
00195      _extents = _bbox.extents();
00196      _step = _bbox.stepsize();
00197      _bottom = _bbox.bottom();
00198      _size = _bbox.size();
00199 #ifdef DEBUG_PRINT_GD
00200      assert (_data!=0);
00201      fill(9999);
00202      ( comm_service::log() << "Data Check"
00203        << _bbox 
00204        << ": " << _data[_size-1]
00205        << std::endl ).flush();
00206 #endif
00207    }
00208 
00209    /* A "pseudo" copy-constructor */
00210    GridData (const GridData<Type,2> &other) : 
00211      _bbox(other._bbox), _extents(other._extents), _step(other._step),
00212      _size(other._size), _bottom(other._bottom),
00213      _data(_size ? new Type[_size] : ((Type *) 0)) {
00214 #ifdef DEBUG_PRINT_GD
00215      assert (_data!=0);
00216      fill(9999);
00217      ( comm_service::log() << "Data Check"
00218        << _bbox
00219        << ": " << _data[_size-1]
00220        << std::endl ).flush();
00221 #endif
00222 #ifdef DEBUG_PRINT_GD_MEMORY
00223      DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00224 #endif
00225    }
00226 
00227    inline ~GridData(void) { 
00228      if (_data) {
00229        delete [] _data; 
00230 #ifdef DEBUG_PRINT_GD_MEMORY
00231        DAGHMemoryTrace::free(sizeof(Type)*_size);
00232 #endif
00233      }
00234    }
00235 
00236    void allocate(const BBox &bb) {
00237      _bbox = bb;
00238      _bbox.rank = 2; 
00239      _bbox.lower().rank = 2;
00240      _bbox.upper().rank = 2;
00241      _bbox.stepsize().rank = 2;
00242      _extents = _bbox.extents();
00243      _step = _bbox.stepsize();
00244      _bottom = _bbox.bottom();
00245      _size = _bbox.size();
00246      _data = _size ? new Type[_size] : (Type *) 0;
00247 #ifdef DEBUG_PRINT_GD
00248      assert (_data!=0);
00249      fill(9999);
00250      ( comm_service::log() << "Data Check"
00251        << _bbox 
00252        << ": " << _data[_size-1]
00253        << std::endl ).flush();
00254 #endif
00255 #ifdef DEBUG_PRINT_GD_MEMORY
00256      DAGHMemoryTrace::alloc(sizeof(Type)*_size);
00257 #endif
00258    }
00259 
00260   void allocate(const BBox &bb, Type *databuf)  {
00261     _bbox = bb;
00262     _bbox.rank = 2; 
00263     _bbox.lower().rank = 2;
00264     _bbox.upper().rank = 2;
00265     _bbox.stepsize().rank = 2;
00266     _extents = _bbox.extents();
00267     _step = _bbox.stepsize();
00268     _bottom = _bbox.bottom();
00269     _size = _bbox.size();
00270 #ifdef DEBUG_PRINT_GD
00271     assert(_data==0);
00272 #endif
00273     _data = databuf;
00274 #ifdef DEBUG_PRINT_GD
00275     assert (_data!=0);
00276     fill(9999);
00277     ( comm_service::log() << "Data Check"
00278       << _bbox 
00279       << ": " << _data[_size-1]
00280       << std::endl ).flush();
00281 #endif
00282   }
00283   
00284   inline void allocate(Type *databuf) {
00285     if (_data) {
00286       delete [] _data; _data = (Type*)0;
00287 #ifdef DEBUG_PRINT_GD_MEMORY
00288       DAGHMemoryTrace::free(sizeof(Type)*_size);
00289 #endif
00290     }
00291     _data = databuf;
00292   }
00293 
00294   inline void deallocate() {
00295     if (_data) {
00296       delete [] _data; _data = (Type*)0;
00297 #ifdef DEBUG_PRINT_GD_MEMORY
00298       DAGHMemoryTrace::free(sizeof(Type)*_size);
00299 #endif
00300     }
00301   }
00302   inline void deallocate(Type*& databuf)
00303   { databuf = _data; _data = (Type*) 0; }
00304   
00305   void* databuffer() { return ((void*) _data); }
00306 
00307   /* Query funtions */
00308   /* Define some inline functions for getting at region extents */
00309   
00310   inline const Coords& lower() const { return(_bbox.lower()); }
00311   inline const Coords& upper() const { return(_bbox.upper()); }
00312   inline const Coords& extents() const { return(_extents); }
00313   inline const Coords& stepsize() const { return(_step); }
00314   inline Coords lower() { return(_bbox.lower()); }
00315   inline Coords upper() { return(_bbox.upper()); }
00316   inline Coords extents() { return(_extents); }
00317   inline Coords stepsize() { return(_step); }
00318   inline int bottom() const { return(_bottom); }
00319   inline int lower(const int i) const { return(_bbox.lower(i)); }
00320   inline int upper(const int i) const { return(_bbox.upper(i)); }
00321   inline int extents(const int i) const { return(_extents(i)); }
00322   inline int stepsize(const int i) const { return(_step(i)); }
00323   inline const BBox& bbox() const { return(_bbox); }
00324 
00325   inline int ok_to_index() { return(_data != (Type *) 0); }
00326   inline int size() const { return(_size); }
00327 
00328   inline int idx(const int i, const int j) const
00329   { return( _bottom+(i/_step(0)+_extents(0)*(j/_step(1))) ); }
00330 
00331   inline const Type& operator () (const int i, const int j) const
00332   { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data[idx(i,j)]); }
00333   inline Type& operator () (const int i,  const int j)
00334   { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data[idx(i,j)]); }
00335 
00336   inline const Type& operator () (const Coords& c) const
00337   { assert(idx(c(0),c(1)) >= 0 && idx(c(0),c(1)) < _size); return(_data[idx(c(0),c(1))]); }
00338   inline Type& operator () (const Coords& c)
00339   { assert(idx(c(0),c(1)) >= 0 && idx(c(0),c(1)) < _size); return(_data[idx(c(0),c(1))]); }
00340 
00341   inline const Type *ptr(const int i,  const int j) const
00342   { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data + idx(i,j)); }
00343   inline Type *ptr(const int i,  const int j)
00344   { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data + idx(i,j)); }
00345 
00346   inline const Type *ptr(const Coords& c) const
00347   { assert(idx(c(0),c(1)) >= 0 && idx(c(0),c(1)) < _size); return(_data + idx(c(0),c(1))); }
00348   inline Type *ptr(const Coords& c)
00349   { assert(idx(c(0),c(1)) >= 0 && idx(c(0),c(1)) < _size); return(_data + idx(c(0),c(1))); }
00350 
00351   inline const Type *data() const { return(_data); }
00352   inline Type *data() { return(_data); }
00353 
00354   void fill(const Type &val)  
00355   { register int i; for (i = 0; i < _size; i++) _data[i] = val; }
00356 
00357   /* Define the copy functions for communication between objects */
00358   void copy(const GridData<Type,2> &gd) {
00359     if (&gd != this) {
00360       BBox intersection = _bbox * gd._bbox;
00361       if (!intersection.empty()) { 
00362         Coords max_step = _step.getmax(gd._step);
00363         BBox to(intersection); to.setstepsize(_step); 
00364         BBox from(intersection); from.setstepsize(gd._step);
00365         gd_CopyRegion(gd, to, from, max_step); 
00366       }
00367     }
00368   }
00369 
00370   void copy(const GridData<Type,2> &gd, const BBox &where) {
00371     if (&gd != this) {
00372       BBox intersection = _bbox * gd._bbox * where;
00373       if (!intersection.empty()) { 
00374         Coords max_step = _step.getmax(gd._step);
00375         BBox to(intersection); to.setstepsize(_step); 
00376         BBox from(intersection); from.setstepsize(gd._step);
00377         gd_CopyRegion(gd, to, from, max_step); 
00378       }
00379     }
00380   }
00381 
00382   void copy(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00383     const Coords  toshift = from.lower() - to.lower();
00384     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00385     if (!newfrom.empty()) {
00386       BBox newto = shiftabs(newfrom, -toshift);
00387       Coords max_step = _step.getmax(gd._step);
00388       newto.setstepsize(_step);
00389       newfrom.setstepsize(gd._step);
00390       gd_CopyRegion(gd, newto, newfrom, max_step);
00391     }
00392   }
00393 
00394   void copy(const GridDataBucket<Type> &gdbkt) {
00395     const BBox &gdbktbb = gdbkt.bbox();
00396     BBox intersection = _bbox * gdbktbb;
00397     if (!intersection.empty()) { 
00398       Coords max_step = _step.getmax(gdbktbb.stepsize());
00399       BBox to(intersection); to.setstepsize(_step); 
00400       BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
00401       gdb_CopyRegion(gdbkt, to, from, max_step); 
00402     }
00403   }
00404 
00405   void copy(const GridDataBucket<Type> &gdbkt, const BBox &where) {
00406     const BBox &gdbktbb = gdbkt.bbox();
00407     BBox intersection = _bbox * gdbktbb * where;
00408     if (!intersection.empty()) { 
00409       Coords max_step = _step.getmax(gdbktbb.stepsize());
00410       BBox to(intersection); to.setstepsize(_step); 
00411       BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
00412       gdb_CopyRegion(gdbkt, to, from, max_step); 
00413     }
00414   }
00415 
00416   void copy(const GridDataBucket<Type> &gdbkt, const BBox &to, 
00417             const BBox &from) {
00418     const BBox &gdbktbb = gdbkt.bbox();
00419     const Coords  toshift = from.lower() - to.lower();
00420     BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
00421     if (!newfrom.empty()) {
00422       BBox newto = shiftabs(newfrom, -toshift);
00423       Coords max_step = _step.getmax(gdbktbb.stepsize());
00424       newto.setstepsize(_step);
00425       newfrom.setstepsize(gdbktbb.stepsize());
00426       gdb_CopyRegion(gdbkt, newto, newfrom, max_step);
00427     }
00428   }
00429 
00430   void copy(const GridDataBucket<Type> &gdbkt, const int n) {
00431     const BBox &gdbktbb = gdbkt.bbox(n);
00432     BBox intersection = _bbox * gdbktbb;
00433     if (!intersection.empty()) { 
00434       Coords max_step = _step.getmax(gdbktbb.stepsize());
00435       BBox to(intersection); to.setstepsize(_step); 
00436       BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
00437       gdb_CopyRegion(gdbkt, n, to, from, max_step); 
00438     }
00439   }
00440 
00441   void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &where)  {
00442     const BBox &gdbktbb = gdbkt.bbox(n);
00443     BBox intersection = _bbox * gdbktbb * where;
00444     if (!intersection.empty()) { 
00445       Coords max_step = _step.getmax(gdbktbb.stepsize());
00446       BBox to(intersection); to.setstepsize(_step); 
00447       BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
00448       gdb_CopyRegion(gdbkt, n, to, from, max_step); 
00449     }
00450   }
00451 
00452   void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to, 
00453             const BBox &from) {
00454     const BBox &gdbktbb = gdbkt.bbox(n);
00455     const Coords  toshift = from.lower() - to.lower();
00456     BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
00457     if (!newfrom.empty()) {
00458       BBox newto = shiftabs(newfrom, -toshift);
00459       Coords max_step = _step.getmax(gdbktbb.stepsize());
00460       newto.setstepsize(_step);
00461       newfrom.setstepsize(gdbktbb.stepsize());
00462       gdb_CopyRegion(gdbkt, n, newto, newfrom, max_step);
00463     }
00464   }
00465 
00466   /* Linear interpolation */
00467   void lin_interp(const GridData<Type,2> &gd1, double const frac1,
00468                   const GridData<Type,2> &gd2, double const frac2,
00469                   const BBox &where) {
00470     BBox intersection = _bbox * gd1._bbox * gd2._bbox * where;
00471     if (!intersection.empty()) {
00472       Coords max_step = _step.getmax(gd1._step);
00473       max_step.max(gd2._step);
00474 
00475       GridData<Type,2> &dst = *this;
00476 
00477       BeginFastIndex2(src1, gd1._bbox, gd1._data, const Type);
00478       BeginFastIndex2(src2, gd2._bbox, gd2._data, const Type);
00479       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
00480       
00481       for_2(i, j, intersection, max_step)
00482         FastIndex2(dst,i,j) = (Type) ((FastIndex2(src1,i,j) * frac1)
00483                                       + (FastIndex2(src2,i,j) * frac2));
00484       end_for
00485 
00486       EndFastIndex2(dst);
00487       EndFastIndex2(src1);
00488       EndFastIndex2(src2);
00489     }
00490   }
00491 
00492   inline void lin_interp(const GridData<Type,2> &gd1, double const frac1,
00493                          const GridData<Type,2> &gd2, double const frac2)
00494   { GridData<Type,2>::lin_interp(gd1,frac1,gd2,frac2,_bbox); }
00495 
00496   /* first moment (i think) sum(q(i,j)*i) */
00497   Type moment1(const int axis, const BBox &where) {
00498     Type m1 = (Type) 0;
00499     BBox intersection = _bbox * where;
00500     if (!intersection.empty()) {
00501       Coords max_step = _step.getmax(where.stepsize());
00502       intersection.setstepsize(_step);
00503       
00504       GridData<Type,2> &dst = *this;
00505       
00506       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
00507 
00508       if (axis == DAGH_X) {
00509         for_2(i, j, intersection, max_step)
00510           m1 += FastIndex2(dst,i,j) * i;
00511         end_for
00512       }
00513       else if (axis == DAGH_Y) {
00514         for_2(i, j, intersection, max_step)
00515           m1 += FastIndex2(dst,i,j) * j;
00516         end_for
00517       }
00518 
00519       EndFastIndex2(dst);
00520     }
00521     return (m1); 
00522   }
00523 
00524   inline Type moment1(const int axis)
00525   { return (GridData<Type,2>::moment1(axis,_bbox)); }
00526 
00527   /* sum of squares - for calculating L2-norm */
00528   double sumsqrd(const BBox &where) {
00529     double s = 0.0;
00530     BBox intersection = _bbox * where;
00531     if (!intersection.empty()) {
00532       Coords max_step = _step.getmax(where.stepsize());
00533       intersection.setstepsize(_step);
00534 
00535       GridData<Type,2> &dst = *this;
00536 
00537       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
00538 
00539       for_2(i, j, intersection, max_step)
00540         s += Abs(FastIndex2(dst,i,j)) * Abs(FastIndex2(dst,i,j));
00541       end_for
00542 
00543       EndFastIndex2(dst);
00544 
00545     }
00546     return s;
00547   }
00548 
00549   inline double sumsqrd()
00550   { return (GridData<Type,2>::sumsqrd(_bbox)); }
00551 
00552   /* sum of absolute values - for calculating L1-norm */
00553   double sumabs(const BBox &where) {
00554     double s = 0.0;
00555     BBox intersection = _bbox * where;
00556     if (!intersection.empty()) {
00557       Coords max_step = _step.getmax(where.stepsize());
00558       intersection.setstepsize(_step);
00559       
00560       GridData<Type,2> &dst = *this;
00561       
00562       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
00563       
00564       for_2(i, j, intersection, max_step)
00565         s += Abs(FastIndex2(dst,i,j));
00566       end_for
00567 
00568       EndFastIndex2(dst);
00569 
00570     }
00571     return s;
00572   }
00573 
00574   inline double sumabs()
00575   { return (GridData<Type,2>::sumabs(_bbox)); }
00576 
00577   /* Maximum of absolute values - for calculating Loo-norm */
00578   double maxabs(const BBox &where) {
00579     double m = 0.0;
00580     BBox intersection = _bbox * where;
00581     if (!intersection.empty()) {
00582       Coords max_step = _step.getmax(where.stepsize());
00583       intersection.setstepsize(_step);
00584       
00585       GridData<Type,2> &dst = *this;
00586 
00587       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
00588       
00589       for_2(i, j, intersection, max_step)
00590         double t = Abs(FastIndex2(dst,i,j));
00591         m = Max(m,t);
00592       end_for
00593 
00594       EndFastIndex2(dst);
00595 
00596     }
00597     return m;
00598   }
00599 
00600   inline double maxabs()
00601   { return (GridData<Type,2>::maxabs(_bbox)); }
00602 
00603   /* operator == */
00604   BBox is_eq(const Type &val, const BBox &where) {
00605     BBox intersection = _bbox * where;
00606     if (!intersection.empty()) {
00607       Coords max_step = _step.getmax(where.stepsize());
00608       BBox to(intersection); to.setstepsize(_step);
00609       BBox bb(2,max_step);
00610       gd_OperateRegion(is_eq)(val, to, max_step, bb);
00611       return (bb);
00612     }
00613     else 
00614       { BBox bb(2,1); return (bb); }
00615   }
00616 
00617   inline BBox is_eq(const Type &val) 
00618   { return (GridData<Type,2>::is_eq(val, _bbox)); }
00619 
00620   BBox is_eq(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00621     const Coords  toshift = from.lower() - to.lower();
00622     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00623     if (!newfrom.empty()) {
00624       BBox newto = shiftabs(newfrom, -toshift);
00625       Coords max_step = _step.getmax(gd._step);
00626       newto.setstepsize(_step);
00627       newfrom.setstepsize(gd._step);
00628       BBox bb(2,max_step);
00629       gd_OperateRegion(is_eq)(gd, newto, newfrom, max_step, bb);
00630       return (bb);
00631     }
00632     else 
00633       { BBox bb(2,1); return (bb); }
00634   }
00635 
00636   inline BBox is_eq(const GridData<Type,2> &gd)
00637   { return (GridData<Type,2>::is_eq(gd, gd._bbox, gd._bbox)); }
00638   inline BBox is_eq(const GridData<Type,2> &gd, const BBox &where)
00639   { return (GridData<Type,2>::is_eq(gd, where, where)); }
00640 
00641   /* operator != */
00642   BBox is_neq(const Type &val, const BBox &where) {
00643     BBox intersection = _bbox * where;
00644     if (!intersection.empty()) {
00645       Coords max_step = _step.getmax(where.stepsize());
00646       BBox to(intersection); to.setstepsize(_step);
00647       BBox bb(2,max_step);
00648       gd_OperateRegion(is_neq)(val, to, max_step, bb);
00649       return (bb);
00650     }
00651     else 
00652       { BBox bb(2,1); return (bb); }
00653   }
00654 
00655   inline BBox is_neq(const Type &val) 
00656   { return (GridData<Type,2>::is_neq(val, _bbox)); }
00657 
00658   BBox is_neq(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00659     const Coords  toshift = from.lower() - to.lower();
00660     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00661     if (!newfrom.empty()) {
00662       BBox newto = shiftabs(newfrom, -toshift);
00663       Coords max_step = _step.getmax(gd._step);
00664       newto.setstepsize(_step);
00665       newfrom.setstepsize(gd._step);
00666       BBox bb(2,max_step);
00667       gd_OperateRegion(is_neq)(gd, newto, newfrom, max_step, bb);
00668       return (bb);
00669     }
00670     else 
00671       { BBox bb(2,1); return (bb); }
00672   }
00673 
00674   inline BBox is_neq(const GridData<Type,2> &gd)
00675   { return (GridData<Type,2>::is_neq(gd, gd._bbox, gd._bbox)); }
00676   inline BBox is_neq(const GridData<Type,2> &gd, const BBox &where)
00677   { return (GridData<Type,2>::is_neq(gd, where, where)); }
00678 
00679   /* operator > */
00680   BBox is_gt(const Type &val, const BBox &where) {
00681     BBox intersection = _bbox * where;
00682     if (!intersection.empty()) {
00683       Coords max_step = _step.getmax(where.stepsize());
00684       BBox to(intersection); to.setstepsize(_step);
00685       BBox bb(2,max_step);
00686       gd_OperateRegion(is_gt)(val, to, max_step, bb);
00687       return (bb);
00688     }
00689     else 
00690       { BBox bb(2,1); return (bb); }
00691   }
00692 
00693   inline BBox is_gt(const Type &val) 
00694   { return (GridData<Type,2>::is_gt(val, _bbox)); }
00695 
00696   BBox is_gt(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00697     const Coords  toshift = from.lower() - to.lower();
00698     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00699     if (!newfrom.empty()) {
00700       BBox newto = shiftabs(newfrom, -toshift);
00701       Coords max_step = _step.getmax(gd._step);
00702       newto.setstepsize(_step);
00703       newfrom.setstepsize(gd._step);
00704       BBox bb(2,max_step);
00705       gd_OperateRegion(is_gt)(gd, newto, newfrom, max_step, bb);
00706       return (bb);
00707     }
00708     else 
00709       { BBox bb(2,1); return (bb); }
00710   }
00711 
00712   inline BBox is_gt(const GridData<Type,2> &gd)
00713   { return (GridData<Type,2>::is_gt(gd, gd._bbox, gd._bbox)); }
00714   inline BBox is_gt(const GridData<Type,2> &gd, const BBox &where)
00715   { return (GridData<Type,2>::is_gt(gd, where, where)); }
00716   
00717   /* operator >= */
00718   BBox is_ge(const Type &val, const BBox &where) {
00719     BBox intersection = _bbox * where;
00720     if (!intersection.empty()) {
00721       Coords max_step = _step.getmax(where.stepsize());
00722       BBox to(intersection); to.setstepsize(_step);
00723       BBox bb(2,max_step);
00724       gd_OperateRegion(is_ge)(val, to, max_step, bb);
00725       return (bb);
00726     }
00727     else 
00728       { BBox bb(2,1); return (bb); }
00729   }
00730 
00731   inline BBox is_ge(const Type &val) 
00732   { return (GridData<Type,2>::is_ge(val, _bbox)); }
00733 
00734   BBox is_ge(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00735     const Coords  toshift = from.lower() - to.lower();
00736     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00737     if (!newfrom.empty()) {
00738       BBox newto = shiftabs(newfrom, -toshift);
00739       Coords max_step = _step.getmax(gd._step);
00740       newto.setstepsize(_step);
00741       newfrom.setstepsize(gd._step);
00742       BBox bb(2,max_step);
00743       gd_OperateRegion(is_ge)(gd, newto, newfrom, max_step, bb);
00744       return (bb);
00745     }
00746     else 
00747       { BBox bb(2,1); return (bb); }
00748   }
00749 
00750   inline BBox is_ge(const GridData<Type,2> &gd)
00751   { return (GridData<Type,2>::is_ge(gd, gd._bbox, gd._bbox)); }
00752   inline BBox is_ge(const GridData<Type,2> &gd, const BBox &where)
00753   { return (GridData<Type,2>::is_ge(gd, where, where)); }
00754 
00755   /* operator < */
00756   BBox is_lt(const Type &val, const BBox &where) {
00757     BBox intersection = _bbox * where;
00758     if (!intersection.empty()) {
00759       Coords max_step = _step.getmax(where.stepsize());
00760       BBox to(intersection); to.setstepsize(_step);
00761       BBox bb(2,max_step);
00762       gd_OperateRegion(is_lt)(val, to, max_step, bb);
00763       return (bb);
00764     }
00765     else 
00766       { BBox bb(2,1); return (bb); }
00767   }
00768 
00769   inline BBox is_lt(const Type &val) 
00770   { return (GridData<Type,2>::is_lt(val, _bbox)); }
00771 
00772   BBox is_lt(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00773     const Coords  toshift = from.lower() - to.lower();
00774     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00775     if (!newfrom.empty()) {
00776       BBox newto = shiftabs(newfrom, -toshift);
00777       Coords max_step = _step.getmax(gd._step);
00778       newto.setstepsize(_step);
00779       newfrom.setstepsize(gd._step);
00780       BBox bb(2,max_step);
00781       gd_OperateRegion(is_lt)(gd, newto, newfrom, max_step, bb);
00782       return (bb);
00783     }
00784     else 
00785       { BBox bb(2,1); return (bb); }
00786   }
00787 
00788   inline BBox is_lt(const GridData<Type,2> &gd)
00789   { return (GridData<Type,2>::is_lt(gd, gd._bbox, gd._bbox)); }
00790   inline BBox is_lt(const GridData<Type,2> &gd, const BBox &where)
00791   { return (GridData<Type,2>::is_lt(gd, where, where)); }
00792 
00793   /* operator <= */
00794   BBox is_le(const Type &val, const BBox &where) {
00795     BBox intersection = _bbox * where;
00796     if (!intersection.empty()) {
00797       Coords max_step = _step.getmax(where.stepsize());
00798       BBox to(intersection); to.setstepsize(_step);
00799       BBox bb(2,max_step);
00800       gd_OperateRegion(is_le)(val, to, max_step, bb);
00801       return (bb);
00802     }
00803     else 
00804       { BBox bb(2,1); return (bb); }
00805   }
00806 
00807   inline BBox is_le(const Type &val) 
00808   { return (GridData<Type,2>::is_le(val, _bbox)); }
00809 
00810   BBox is_le(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00811     const Coords  toshift = from.lower() - to.lower();
00812     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00813     if (!newfrom.empty()) {
00814       BBox newto = shiftabs(newfrom, -toshift);
00815       Coords max_step = _step.getmax(gd._step);
00816       newto.setstepsize(_step);
00817       newfrom.setstepsize(gd._step);
00818       BBox bb(2,max_step);
00819       gd_OperateRegion(is_le)(gd, newto, newfrom, max_step, bb);
00820       return (bb);
00821     }
00822     else 
00823       { BBox bb(2,1); return (bb); }
00824   }
00825 
00826   inline BBox is_le(const GridData<Type,2> &gd)
00827   { return (GridData<Type,2>::is_le(gd, gd._bbox, gd._bbox)); }
00828   inline BBox is_le(const GridData<Type,2> &gd, const BBox &where)
00829   { return (GridData<Type,2>::is_le(gd, where, where)); }
00830 
00831   /* operator = */
00832   void equals(const Type &val, const BBox &where) {
00833     BBox intersection = _bbox * where;
00834     if (!intersection.empty()) {
00835       Coords max_step = _step.getmax(where.stepsize());
00836       BBox to(intersection); to.setstepsize(_step);
00837       gd_OperateRegion(equal)(val, to, max_step);
00838     }
00839   }
00840 
00841   inline void equals(const Type &val) 
00842   { GridData<Type,2>::equals(val, _bbox); }
00843 
00844   void equals(const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00845     const Coords  toshift = from.lower() - to.lower();
00846     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00847     if (!newfrom.empty()) {
00848       BBox newto = shiftabs(newfrom, -toshift);
00849       Coords max_step = _step.getmax(gd._step);
00850       newto.setstepsize(_step);
00851       newfrom.setstepsize(gd._step);
00852       gd_OperateRegion(equal)(gd, newto, newfrom, max_step);
00853     }
00854   }
00855   
00856   inline void equals(const GridData<Type,2> &gd)
00857   { GridData<Type,2>::equals(gd, gd._bbox, gd._bbox); }
00858   inline void equals(const GridData<Type,2> &gd, const BBox &where)
00859   { GridData<Type,2>::equals(gd, where, where); }
00860 
00861   /* operator += */
00862   void plus (const Type &val, const BBox &where) {
00863     BBox intersection = _bbox * where;
00864     if (!intersection.empty()) {
00865       Coords max_step = _step.getmax(where.stepsize());
00866       BBox to(intersection); to.setstepsize(_step);
00867       gd_OperateRegion(plus)(val, to, max_step);
00868     }
00869   }
00870 
00871   inline void plus (const Type &val) 
00872   { GridData<Type,2>::plus(val, _bbox); }
00873 
00874   void plus (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00875     const Coords  toshift = from.lower() - to.lower();
00876     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00877     if (!newfrom.empty()) {
00878       BBox newto = shiftabs(newfrom, -toshift);
00879       Coords max_step = _step.getmax(gd._step);
00880       newto.setstepsize(_step);
00881       newfrom.setstepsize(gd._step);
00882       gd_OperateRegion(plus)(gd, newto, newfrom, max_step);
00883     }
00884   }
00885 
00886   inline void plus (const GridData<Type,2> &gd)
00887   { GridData<Type,2>::plus(gd, gd._bbox, gd._bbox); }
00888   inline void plus (const GridData<Type,2> &gd, const BBox &where)
00889   { GridData<Type,2>::plus(gd, where, where); }
00890 
00891   /* operator -= */
00892   void minus (const Type &val, const BBox &where)  {
00893     BBox intersection = _bbox * where;
00894     if (!intersection.empty()) {
00895       Coords max_step = _step.getmax(where.stepsize());
00896       BBox to(intersection); to.setstepsize(_step);
00897       gd_OperateRegion(minus)(val, to, max_step);
00898     }
00899   }
00900 
00901   inline void minus (const Type &val) 
00902   { GridData<Type,2>::minus(val, _bbox); }
00903 
00904   void minus (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00905     const Coords  toshift = from.lower() - to.lower();
00906     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00907     if (!newfrom.empty()) {
00908       BBox newto = shiftabs(newfrom, -toshift);
00909       Coords max_step = _step.getmax(gd._step);
00910       newto.setstepsize(_step);
00911       newfrom.setstepsize(gd._step);
00912       gd_OperateRegion(minus)(gd, newto, newfrom, max_step);
00913     }
00914   }
00915 
00916   inline void minus (const GridData<Type,2> &gd)
00917   { GridData<Type,2>::minus(gd, gd._bbox, gd._bbox); }
00918   inline void minus (const GridData<Type,2> &gd, const BBox &where)
00919   { GridData<Type,2>::minus(gd, where, where); }
00920 
00921    /* operator *= */
00922    void multiply(const Type &val, const BBox &where) {
00923      BBox intersection = _bbox * where;
00924      if (!intersection.empty()) {
00925        Coords max_step = _step.getmax(where.stepsize());
00926        BBox to(intersection); to.setstepsize(_step);
00927        gd_OperateRegion(mult)(val, to, max_step);
00928      }
00929    }
00930 
00931   inline void multiply(const Type &val) 
00932   { GridData<Type,2>::multiply(val, _bbox); }
00933 
00934   void multiply (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00935     const Coords  toshift = from.lower() - to.lower();
00936     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00937     if (!newfrom.empty()) {
00938       BBox newto = shiftabs(newfrom, -toshift);
00939       Coords max_step = _step.getmax(gd._step);
00940       newto.setstepsize(_step);
00941       newfrom.setstepsize(gd._step);
00942       gd_OperateRegion(mult)(gd, newto, newfrom, max_step);
00943     }
00944   }
00945 
00946   inline void multiply (const GridData<Type,2> &gd)
00947   { GridData<Type,2>::multiply(gd, gd._bbox, gd._bbox); }
00948   inline void multiply (const GridData<Type,2> &gd, const BBox &where)
00949   { GridData<Type,2>::multiply(gd, where, where); }
00950 
00951   /* operator /= */
00952   void divide (const Type &val, const BBox &where) {
00953     BBox intersection = _bbox * where;
00954     if (!intersection.empty()) {
00955       Coords max_step = _step.getmax(where.stepsize());
00956       BBox to(intersection); to.setstepsize(_step);
00957       gd_OperateRegion(div)(val, to, max_step);
00958     }
00959   }
00960 
00961   inline void divide (const Type &val) 
00962   { GridData<Type,2>::divide(val, _bbox); }
00963 
00964   void divide (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00965     const Coords  toshift = from.lower() - to.lower();
00966     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00967     if (!newfrom.empty()) {
00968       BBox newto = shiftabs(newfrom, -toshift);
00969       Coords max_step = _step.getmax(gd._step);
00970       newto.setstepsize(_step);
00971       newfrom.setstepsize(gd._step);
00972       gd_OperateRegion(div)(gd, newto, newfrom, max_step);
00973     }
00974   }
00975 
00976   inline void divide (const GridData<Type,2> &gd)
00977   { GridData<Type,2>::divide(gd, gd._bbox, gd._bbox); }
00978   inline void divide (const GridData<Type,2> &gd, const BBox &where)
00979   { GridData<Type,2>::divide(gd, where, where); }
00980 
00981   /* point-wise minimum */
00982   void minimum (const Type &val, const BBox &where) {
00983     BBox intersection = _bbox * where;
00984     if (!intersection.empty()) {
00985       Coords max_step = _step.getmax(where.stepsize());
00986       BBox to(intersection); to.setstepsize(_step);
00987       gd_OperateRegion(min)(val, to, max_step);
00988     }
00989   }
00990 
00991   inline void minimum (const Type &val) 
00992   { GridData<Type,2>::minimum(val, _bbox); }
00993 
00994   void minimum (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
00995     const Coords  toshift = from.lower() - to.lower();
00996     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
00997     if (!newfrom.empty()) {
00998       BBox newto = shiftabs(newfrom, -toshift);
00999       Coords max_step = _step.getmax(gd._step);
01000       newto.setstepsize(_step);
01001       newfrom.setstepsize(gd._step);
01002       gd_OperateRegion(min)(gd, newto, newfrom, max_step);
01003     }
01004   }
01005 
01006   inline void minimum (const GridData<Type,2> &gd)
01007   { GridData<Type,2>::minimum(gd, gd._bbox, gd._bbox); }
01008   inline void minimum (const GridData<Type,2> &gd, const BBox &where)
01009   { GridData<Type,2>::minimum(gd, where, where); }
01010 
01011   /* point-wise maximum */
01012   void maximum (const Type &val, const BBox &where) {
01013     BBox intersection = _bbox * where;
01014     if (!intersection.empty()) {
01015       Coords max_step = _step.getmax(where.stepsize());
01016       BBox to(intersection); to.setstepsize(_step);
01017       gd_OperateRegion(max)(val, to, max_step);
01018     }
01019   }
01020 
01021   inline void maximum (const Type &val) 
01022   { GridData<Type,2>::maximum(val, _bbox); }
01023 
01024   void maximum (const GridData<Type,2> &gd, const BBox &to, const BBox &from) {
01025     const Coords  toshift = from.lower() - to.lower();
01026     BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
01027     if (!newfrom.empty()) {
01028       BBox newto = shiftabs(newfrom, -toshift);
01029       Coords max_step = _step.getmax(gd._step);
01030       newto.setstepsize(_step);
01031       newfrom.setstepsize(gd._step);
01032       gd_OperateRegion(max)(gd, newto, newfrom, max_step);
01033     }
01034   }
01035 
01036   inline void maximum (const GridData<Type,2> &gd)
01037   { GridData<Type,2>::maximum(gd, gd._bbox, gd._bbox); }
01038   inline void maximum (const GridData<Type,2> &gd, const BBox &where)
01039   { GridData<Type,2>::maximum(gd, where, where); }
01040 
01041   /**********************************************************************/
01042   inline BBox operator == (const Type &val)
01043   { return (GridData<Type,2>::is_eq(val, _bbox)); }
01044   inline BBox operator == (const GridData<Type,2> &gd)
01045   { return (GridData<Type,2>::is_eq(gd, gd._bbox, gd._bbox)); }
01046 
01047   inline BBox operator != (const Type &val)
01048   { return (GridData<Type,2>::is_neq(val, _bbox)); }
01049   inline BBox operator != (const GridData<Type,2> &gd)
01050   { return (GridData<Type,2>::is_neq(gd, gd._bbox, gd._bbox)); }
01051 
01052   inline BBox operator > (const Type &val)
01053   { return (GridData<Type,2>::is_gt(val, _bbox)); }
01054   inline BBox operator > (const GridData<Type,2> &gd)
01055   { return (GridData<Type,2>::is_gt(gd, gd._bbox, gd._bbox)); }
01056 
01057   inline BBox operator >= (const Type &val)
01058   { return (GridData<Type,2>::is_ge(val, _bbox)); }
01059   inline BBox operator >= (const GridData<Type,2> &gd)
01060   { return (GridData<Type,2>::is_ge(gd, gd._bbox, gd._bbox)); }
01061 
01062   inline BBox operator < (const Type &val)
01063   { return (GridData<Type,2>::is_lt(val, _bbox)); }
01064   inline BBox operator < (const GridData<Type,2> &gd)
01065   { return (GridData<Type,2>::is_lt(gd, gd._bbox, gd._bbox)); }
01066 
01067   inline BBox operator <= (const Type &val)
01068   { return (GridData<Type,2>::is_le(val, _bbox)); }
01069   inline BBox operator <= (const GridData<Type,2> &gd)
01070   { return (GridData<Type,2>::is_le(gd, gd._bbox, gd._bbox)); }
01071 
01072   inline void operator = (const Type &val)
01073   { GridData<Type,2>::equals(val, _bbox); }
01074   inline void operator = (const GridData<Type,2> &gd)
01075   { GridData<Type,2>::equals(gd, gd._bbox, gd._bbox); }
01076   
01077   inline void operator += (const Type &val)
01078   { GridData<Type,2>::plus(val, _bbox); }
01079   inline void operator += (const GridData<Type,2> &gd)
01080   { GridData<Type,2>::plus(gd, gd._bbox, gd._bbox); }
01081   
01082   inline void operator -= (const Type &val)
01083   { GridData<Type,2>::minus(val, _bbox); }
01084   inline void operator -= (const GridData<Type,2> &gd)
01085   { GridData<Type,2>::minus(gd, gd._bbox, gd._bbox); }
01086   
01087   inline void operator *= (const Type &val)
01088   { GridData<Type,2>::multiply(val, _bbox); }
01089   inline void operator *= (const GridData<Type,2> &gd)
01090   { GridData<Type,2>::multiply(gd, gd._bbox, gd._bbox); }
01091 
01092   inline void operator /= (const Type &val)
01093   { GridData<Type,2>::divide(val, _bbox); }
01094   inline void operator /= (const GridData<Type,2> &gd)
01095   { GridData<Type,2>::divide(gd, gd._bbox, gd._bbox); }
01096   /**********************************************************************/
01097   
01098   /**********************************************************************/
01099   // Reduction Operations
01100   /**********************************************************************/
01101   Type maxval (const BBox &where) {
01102     Type max_val = (Type) 0;
01103     short flag = DAGHTrue;
01104     BBox intersection = _bbox * where;
01105     if (!intersection.empty()) {
01106       Coords max_step = _step.getmax(where.stepsize());
01107       intersection.setstepsize(_step);
01108       
01109       GridData<Type,2> &dst = *this;
01110       
01111       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01112       
01113       for_2(i, j, intersection, max_step)
01114         if (flag == DAGHTrue) 
01115           { max_val = FastIndex2(dst,i,j); flag = DAGHFalse; }
01116         else 
01117            max_val = Max(max_val, FastIndex2(dst,i,j));
01118       end_for
01119 
01120       EndFastIndex2(dst);
01121 
01122     }
01123     return (max_val);
01124   }
01125 
01126   inline Type maxval (void)
01127   { return (GridData<Type,2>::maxval(_bbox)); }
01128   
01129   /****************************** minval ***************************************/
01130   Type minval (const BBox &where) {
01131     Type min_val = (Type) 0;
01132     short flag = DAGHTrue;
01133     BBox intersection = _bbox * where;
01134     if (!intersection.empty()) {
01135       Coords max_step = _step.getmax(where.stepsize());
01136       intersection.setstepsize(_step);
01137       
01138       GridData<Type,2> &dst = *this;
01139 
01140       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01141       
01142       for_2(i, j, intersection, max_step)
01143         if (flag == DAGHTrue) 
01144           { min_val = FastIndex2(dst,i,j); flag = DAGHFalse; }
01145         else 
01146           min_val = Min(min_val, FastIndex2(dst,i,j));
01147       end_for
01148 
01149       EndFastIndex2(dst);
01150 
01151     }
01152     return (min_val);
01153   }
01154 
01155   inline Type minval (void)
01156   { return (GridData<Type,2>::minval(_bbox)); }  
01157   /************************************************************************/
01158 
01159   /****************************** sum ***************************************/
01160   Type sum (const BBox &where) {
01161     Type sum_val = (Type) 0;
01162     BBox intersection = _bbox * where;
01163     if (!intersection.empty()) {
01164       Coords max_step = _step.getmax(where.stepsize());
01165       intersection.setstepsize(_step);
01166       
01167       GridData<Type,2> &dst = *this;
01168       
01169       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01170 
01171       for_2(i, j, intersection, max_step)
01172         sum_val += FastIndex2(dst,i,j);
01173       end_for
01174 
01175       EndFastIndex2(dst);
01176 
01177     }
01178     return (sum_val);
01179   }
01180 
01181   inline Type sum (void)
01182   { return (GridData<Type,2>::sum(_bbox)); }  
01183   /************************************************************************/
01184 
01185   /****************************** product ***************************************/
01186   Type product (const BBox &where) {
01187     Type prod_val = (Type) 1;
01188     BBox intersection = _bbox * where;
01189     if (!intersection.empty()) {
01190       Coords max_step = _step.getmax(where.stepsize());
01191       intersection.setstepsize(_step);
01192       
01193       GridData<Type,2> &dst = *this;
01194       
01195       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01196 
01197       for_2(i, j, intersection, max_step)
01198         prod_val *= FastIndex2(dst,i,j);
01199       end_for
01200 
01201       EndFastIndex2(dst);
01202 
01203     }
01204     return (prod_val);
01205   }
01206 
01207   inline Type product (void)
01208   { return (GridData<Type,2>::product(_bbox)); }
01209   /**********************************************************************/
01210 
01211 private:
01212   void gd_CopyRegion(const GridData<Type,2> &src, const BBox &to,
01213                      const BBox &from, Coords const &step) {
01214     GridData<Type,2> &dst = *this;
01215 
01216     if (dst._bbox==src._bbox && to==from && to==dst._bbox) 
01217       std::memcpy((void *) dst._data, (void *) src._data, sizeof(Type)*to.size());
01218     else {
01219       const int di = from.lower(0)-to.lower(0);
01220       const int dj = from.lower(1)-to.lower(1);
01221 
01222       BeginFastIndex2(src, src._bbox, src._data, const Type);
01223       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01224      
01225       for_2(i, j, to, step)
01226         FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
01227       end_for
01228 
01229       EndFastIndex2(dst);
01230       EndFastIndex2(src);
01231     }
01232   }
01233 
01234   void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const BBox &to,
01235                       const BBox &from, Coords const &step) {
01236     GridData<Type,2> &dst = *this;
01237 
01238     if (dst._bbox==gdbkt.bbox() && to==from && to==dst._bbox) 
01239       std::memcpy((void *) dst._data, (void *) gdbkt.data(), sizeof(Type)*to.size());
01240     else {
01241       const int di = from.lower(0)-to.lower(0);
01242       const int dj = from.lower(1)-to.lower(1);
01243 
01244       BeginFastIndex2(src, gdbkt.bbox(), gdbkt.data(), const Type);
01245       BeginFastIndex2(dst, dst.bbox(), dst.data(), Type);
01246       for_2(i, j, to, step)
01247         FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
01248       end_for
01249       EndFastIndex2(dst);
01250       EndFastIndex2(src);
01251     }
01252   }
01253 
01254   void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to,
01255                       const BBox &from, Coords const &step) {
01256     GridData<Type,2> &dst = *this;
01257 
01258     if (dst._bbox==gdbkt.bbox(n) && to==from && to==dst._bbox) 
01259       std::memcpy((void *) dst._data, (void *) gdbkt.data(n), sizeof(Type)*to.size());
01260     else {
01261       const int di = from.lower(0)-to.lower(0);
01262       const int dj = from.lower(1)-to.lower(1);
01263 
01264       BeginFastIndex2(src, gdbkt.bbox(n), gdbkt.data(n), const Type);
01265       BeginFastIndex2(dst, dst.bbox(), dst.data(), Type);
01266       for_2(i, j, to, step)
01267         FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
01268       end_for
01269       EndFastIndex2(dst);
01270       EndFastIndex2(src);
01271     }
01272   }
01273 
01274   /****************************** == ***************************************/
01275   void gd_OperateRegion(is_eq)(const Type &val, const BBox &to,
01276                                Coords const &step, BBox &bb) {
01277     GridData<Type,2> &dst = *this;
01278 
01279     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01280 
01281     for_2(i, j, to, step)
01282       if (FastIndex2(dst,i,j) == val) bb += Coords(2,i,j);
01283     end_for
01284 
01285     EndFastIndex2(dst);
01286   }
01287 
01288   void gd_OperateRegion(is_eq)(const GridData<Type,2> &src, const BBox &to,
01289                                const BBox &from, Coords const &step, BBox &bb) {
01290     GridData<Type,2> &dst = *this;
01291     const int di = from.lower(0)-to.lower(0);
01292     const int dj = from.lower(1)-to.lower(1);
01293     
01294     BeginFastIndex2(src, src._bbox, src._data, const Type);
01295     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01296 
01297     for_2(i, j, to, step)
01298       if (FastIndex2(dst,i,j) == FastIndex2(src,i+di,j+dj))
01299         bb += Coords(2,i,j);
01300     end_for
01301 
01302     EndFastIndex2(dst);
01303     EndFastIndex2(src);
01304   }
01305   /************************************************************************/
01306 
01307   /****************************** != ***************************************/
01308   void gd_OperateRegion(is_neq)(const Type &val, const BBox &to,
01309                                 Coords const &step, BBox &bb) {
01310     GridData<Type,2> &dst = *this;
01311     
01312     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01313 
01314     for_2(i, j, to, step)
01315       if (FastIndex2(dst,i,j) != val) bb += Coords(2,i,j);
01316     end_for
01317 
01318     EndFastIndex2(dst);
01319   }
01320 
01321   void gd_OperateRegion(is_neq)(const GridData<Type,2> &src, const BBox &to,
01322                                 const BBox &from, Coords const &step, BBox &bb) {
01323     GridData<Type,2> &dst = *this;
01324     const int di = from.lower(0)-to.lower(0);
01325     const int dj = from.lower(1)-to.lower(1);
01326     
01327     BeginFastIndex2(src, src._bbox, src._data, const Type);
01328     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01329 
01330     for_2(i, j, to, step)
01331       if (FastIndex2(dst,i,j) != FastIndex2(src,i+di,j+dj))
01332         bb += Coords(2,i,j);
01333     end_for
01334 
01335     EndFastIndex2(dst);
01336     EndFastIndex2(src);
01337   }
01338   /************************************************************************/
01339 
01340   /****************************** > ***************************************/
01341   void gd_OperateRegion(is_gt)(const Type &val, const BBox &to,
01342                                Coords const &step, BBox &bb) {
01343     GridData<Type,2> &dst = *this;
01344     
01345     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01346 
01347     for_2(i, j, to, step)
01348       if (FastIndex2(dst,i,j) > val) bb += Coords(2,i,j);
01349     end_for
01350 
01351     EndFastIndex2(dst);
01352   }
01353 
01354   void gd_OperateRegion(is_gt)(const GridData<Type,2> &src, const BBox &to,
01355                                const BBox &from, Coords const &step, BBox &bb) {
01356     GridData<Type,2> &dst = *this;
01357     const int di = from.lower(0)-to.lower(0);
01358     const int dj = from.lower(1)-to.lower(1);
01359     
01360     BeginFastIndex2(src, src._bbox, src._data, const Type);
01361     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01362     
01363     for_2(i, j, to, step)
01364       if (FastIndex2(dst,i,j) > FastIndex2(src,i+di,j+dj))
01365         bb += Coords(2,i,j);
01366     end_for
01367 
01368     EndFastIndex2(dst);
01369     EndFastIndex2(src);
01370   }
01371   /************************************************************************/
01372 
01373   /****************************** >= ***************************************/
01374   void gd_OperateRegion(is_ge)(const Type &val, const BBox &to,
01375                                Coords const &step, BBox &bb) {
01376     GridData<Type,2> &dst = *this;
01377     
01378     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01379 
01380     for_2(i, j, to, step)
01381       if (FastIndex2(dst,i,j) >= val) bb += Coords(2,i,j);
01382     end_for
01383 
01384     EndFastIndex2(dst);
01385   }
01386 
01387   void gd_OperateRegion(is_ge)(const GridData<Type,2> &src, const BBox &to,
01388                                const BBox &from, Coords const &step, BBox &bb) {
01389     GridData<Type,2> &dst = *this;
01390     const int di = from.lower(0)-to.lower(0);
01391     const int dj = from.lower(1)-to.lower(1);
01392 
01393     BeginFastIndex2(src, src._bbox, src._data, const Type);
01394     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01395 
01396     for_2(i, j, to, step)
01397       if (FastIndex2(dst,i,j) >= FastIndex2(src,i+di,j+dj))
01398         bb += Coords(2,i,j);
01399     end_for
01400 
01401     EndFastIndex2(dst);
01402     EndFastIndex2(src);
01403   }
01404   /************************************************************************/
01405 
01406   /****************************** < ***************************************/
01407   void gd_OperateRegion(is_lt)(const Type &val, const BBox &to,
01408                                Coords const &step, BBox &bb) {
01409     GridData<Type,2> &dst = *this;
01410     
01411     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01412 
01413     for_2(i, j, to, step)
01414       if (FastIndex2(dst,i,j) < val) bb += Coords(2,i,j);
01415     end_for
01416 
01417     EndFastIndex2(dst);
01418   }
01419 
01420   void gd_OperateRegion(is_lt)(const GridData<Type,2> &src, const BBox &to,
01421                                const BBox &from, Coords const &step, BBox &bb) {
01422     GridData<Type,2> &dst = *this;
01423     const int di = from.lower(0)-to.lower(0);
01424     const int dj = from.lower(1)-to.lower(1);
01425 
01426     BeginFastIndex2(src, src._bbox, src._data, const Type);
01427     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01428 
01429     for_2(i, j, to, step)
01430       if (FastIndex2(dst,i,j) < FastIndex2(src,i+di,j+dj))
01431         bb += Coords(2,i,j);
01432     end_for
01433 
01434     EndFastIndex2(dst);
01435     EndFastIndex2(src);
01436   }
01437   /************************************************************************/
01438 
01439   /****************************** <= ***************************************/
01440   void gd_OperateRegion(is_le)(const Type &val, const BBox &to,
01441                                Coords const &step, BBox &bb) {
01442     GridData<Type,2> &dst = *this;
01443 
01444     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01445 
01446     for_2(i, j, to, step)
01447       if (FastIndex2(dst,i,j) <= val) bb += Coords(2,i,j);
01448     end_for
01449 
01450     EndFastIndex2(dst);
01451   }
01452 
01453   void gd_OperateRegion(is_le)(const GridData<Type,2> &src, const BBox &to,
01454                                const BBox &from, Coords const &step, BBox &bb) {
01455     GridData<Type,2> &dst = *this;
01456     const int di = from.lower(0)-to.lower(0);
01457     const int dj = from.lower(1)-to.lower(1);
01458     
01459     BeginFastIndex2(src, src._bbox, src._data, const Type);
01460     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01461 
01462     for_2(i, j, to, step)
01463       if (FastIndex2(dst,i,j) <= FastIndex2(src,i+di,j+dj))
01464         bb += Coords(2,i,j);
01465     end_for
01466 
01467     EndFastIndex2(dst);
01468     EndFastIndex2(src);
01469   }
01470   
01471   /****************************** = ***************************************/
01472   void gd_OperateRegion(equal)(const Type &val, const BBox &to, Coords const &step) {
01473     GridData<Type,2> &dst = *this;
01474 
01475     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01476 
01477     for_2(i, j, to, step)
01478       FastIndex2(dst,i,j) = val;
01479     end_for
01480 
01481     EndFastIndex2(dst);
01482   }
01483 
01484   void gd_OperateRegion(equal)(const GridData<Type,2> &src, const BBox &to,
01485                                const BBox &from, Coords const &step) {
01486     GridData<Type,2> &dst = *this;
01487     const int di = from.lower(0)-to.lower(0);
01488     const int dj = from.lower(1)-to.lower(1);
01489     
01490     BeginFastIndex2(src, src._bbox, src._data, const Type);
01491     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01492 
01493     for_2(i, j, to, step)
01494       FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
01495     end_for
01496 
01497     EndFastIndex2(dst);
01498     EndFastIndex2(src);
01499   }
01500   /************************************************************************/
01501 
01502   /****************************** + ***************************************/
01503   void gd_OperateRegion(plus)(const Type &val, const BBox &to, Coords const &step) {
01504     GridData<Type,2> &dst = *this;
01505 
01506     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01507 
01508     for_2(i, j, to, step)
01509       FastIndex2(dst,i,j) += val;
01510     end_for
01511 
01512     EndFastIndex2(dst);
01513   }
01514 
01515   void gd_OperateRegion(plus)(const GridData<Type,2> &src, const BBox &to,
01516                               const BBox &from, Coords const &step) {
01517     GridData<Type,2> &dst = *this;
01518     const int di = from.lower(0)-to.lower(0);
01519     const int dj = from.lower(1)-to.lower(1);
01520     
01521     BeginFastIndex2(src, src._bbox, src._data, const Type);
01522     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01523 
01524     for_2(i, j, to, step)
01525       FastIndex2(dst,i,j) += FastIndex2(src,i+di,j+dj);
01526     end_for
01527 
01528     EndFastIndex2(dst);
01529     EndFastIndex2(src);
01530   }
01531   /************************************************************************/
01532 
01533   /****************************** - ***************************************/
01534   void gd_OperateRegion(minus)(const Type &val, const BBox &to, Coords const &step) {
01535     GridData<Type,2> &dst = *this;
01536 
01537     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01538 
01539     for_2(i, j, to, step)
01540       FastIndex2(dst,i,j) -= val;
01541     end_for
01542 
01543     EndFastIndex2(dst);
01544   }
01545 
01546   void gd_OperateRegion(minus)(const GridData<Type,2> &src, const BBox &to,
01547                                const BBox &from, Coords const &step) {
01548     GridData<Type,2> &dst = *this;
01549     const int di = from.lower(0)-to.lower(0);
01550     const int dj = from.lower(1)-to.lower(1);
01551     
01552     BeginFastIndex2(src, src._bbox, src._data, const Type);
01553     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01554 
01555     for_2(i, j, to, step)
01556       FastIndex2(dst,i,j) -= FastIndex2(src,i+di,j+dj);
01557     end_for
01558 
01559     EndFastIndex2(dst);
01560     EndFastIndex2(src);
01561   }
01562   /************************************************************************/
01563 
01564   /****************************** * ***************************************/
01565    void gd_OperateRegion(mult)(const Type &val, const BBox &to, Coords const &step) {
01566      GridData<Type,2> &dst = *this;
01567 
01568      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01569 
01570      for_2(i, j, to, step)
01571        FastIndex2(dst,i,j) *= val;
01572      end_for
01573 
01574      EndFastIndex2(dst);
01575    }
01576 
01577    void gd_OperateRegion(mult)(const GridData<Type,2> &src, const BBox &to,
01578                                const BBox &from, Coords const &step) {
01579      GridData<Type,2> &dst = *this;
01580      const int di = from.lower(0)-to.lower(0);
01581      const int dj = from.lower(1)-to.lower(1);
01582      
01583      BeginFastIndex2(src, src._bbox, src._data, const Type);
01584      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01585      
01586      for_2(i, j, to, step)
01587        FastIndex2(dst,i,j) *= FastIndex2(src,i+di,j+dj);
01588      end_for
01589 
01590      EndFastIndex2(dst);
01591      EndFastIndex2(src);
01592    }
01593   /************************************************************************/
01594 
01595   /****************************** / ***************************************/
01596   void gd_OperateRegion(div)(const Type &val, const BBox &to, Coords const &step) {
01597     assert (val != (Type)0);
01598     GridData<Type,2> &dst = *this;
01599 
01600     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01601 
01602     for_2(i, j, to, step)
01603       FastIndex2(dst,i,j) /= val;
01604     end_for
01605 
01606     EndFastIndex2(dst);
01607   }
01608 
01609   void gd_OperateRegion(div)(const GridData<Type,2> &src, const BBox &to,
01610                              const BBox &from, Coords const &step) {
01611     GridData<Type,2> &dst = *this;
01612     const int di = from.lower(0)-to.lower(0);
01613     const int dj = from.lower(1)-to.lower(1);
01614     
01615     BeginFastIndex2(src, src._bbox, src._data, const Type);
01616     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01617 
01618     for_2(i, j, to, step)
01619       FastIndex2(dst,i,j) /= FastIndex2(src,i+di,j+dj);
01620     end_for
01621 
01622     EndFastIndex2(dst);
01623     EndFastIndex2(src);
01624   }
01625   /************************************************************************/
01626 
01627   /****************************** min ***************************************/
01628   void gd_OperateRegion(min)(const Type &val, const BBox &to, Coords const &step) {
01629     GridData<Type,2> &dst = *this;
01630 
01631     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01632 
01633     for_2(i, j, to, step)
01634       FastIndex2(dst,i,j) = Min(val, FastIndex2(dst,i,j));
01635     end_for
01636 
01637     EndFastIndex2(dst);
01638   }
01639 
01640   void gd_OperateRegion(min)(const GridData<Type,2> &src, const BBox &to,
01641                              const BBox &from, Coords const &step) {
01642     GridData<Type,2> &dst = *this;
01643     const int di = from.lower(0)-to.lower(0);
01644     const int dj = from.lower(1)-to.lower(1);
01645     
01646     BeginFastIndex2(src, src._bbox, src._data, const Type);
01647     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01648 
01649     for_2(i, j, to, step)
01650       FastIndex2(dst,i,j) = Min(FastIndex2(dst,i,j), FastIndex2(src,i+di,j+dj));
01651     end_for
01652 
01653     EndFastIndex2(dst);
01654     EndFastIndex2(src);
01655   }
01656   /************************************************************************/
01657 
01658   /****************************** max ***************************************/
01659   void gd_OperateRegion(max)(const Type &val, const BBox &to, Coords const &step) {
01660     GridData<Type,2> &dst = *this;
01661 
01662     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01663 
01664     for_2(i, j, to, step)
01665       FastIndex2(dst,i,j) = Max(val, FastIndex2(dst,i,j));
01666     end_for
01667 
01668     EndFastIndex2(dst);
01669   }
01670 
01671   void gd_OperateRegion(max)(const GridData<Type,2> &src, const BBox &to,
01672                              const BBox &from, Coords const &step) {
01673     GridData<Type,2> &dst = *this;
01674     const int di = from.lower(0)-to.lower(0);
01675     const int dj = from.lower(1)-to.lower(1);
01676     
01677     BeginFastIndex2(src, src._bbox, src._data, const Type);
01678     BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01679 
01680     for_2(i, j, to, step)
01681       FastIndex2(dst,i,j) = Max(FastIndex2(dst,i,j), FastIndex2(src,i+di,j+dj));
01682     end_for
01683 
01684     EndFastIndex2(dst);
01685     EndFastIndex2(src);
01686   }
01687   /************************************************************************/
01688 
01689 public:
01690   /* Perform efficient buffer packing/unpacking and data movement */
01691   void PackRegion(Type *sendbuf, const BBox &from) const {
01692     const GridData<Type,2> &src = *this;
01693 
01694     if (from==src._bbox) 
01695       std::memcpy((void *) sendbuf, (void *) src._data, sizeof(Type)*from.size());
01696     else {
01697       register int p = 0;
01698 
01699       BeginFastIndex2(src, src._bbox, src._data, const Type);
01700 
01701       for_2(i, j, from, _step)
01702         sendbuf[p++] = FastIndex2(src,i,j);
01703       end_for
01704 
01705       EndFastIndex2(src);
01706     }
01707   }
01708 
01709   void UnPackRegion(const Type *recvbuf, const BBox &to) {
01710     GridData<Type,2> &dst = *this;
01711 
01712     if (to==dst._bbox) 
01713       std::memcpy((void *) dst._data, (void *) recvbuf, sizeof(Type)*to.size());
01714     else {
01715       register int p = 0;
01716 
01717       BeginFastIndex2(dst, dst._bbox, dst._data, Type);
01718      
01719       for_2(i, j, to, _step)
01720         FastIndex2(dst,i,j) = recvbuf[p++];
01721       end_for
01722 
01723       EndFastIndex2(dst);
01724     }
01725   }
01726 
01727   friend std::ostream& operator<< <>(std::ostream &os,
01728                                      const GridData<Type,2> &gd);
01729   friend std::ofstream& operator<< <>(std::ofstream &ofs,
01730                                       const GridData<Type,2> &gd);
01731   friend std::ifstream& operator>> <>(std::ifstream &ifs,
01732                                       GridData<Type,2> &gd);
01733   friend std::stringstream& operator<< <>(std::stringstream &ofs,
01734                                           const GridData<Type,2> &gd);
01735   friend std::stringstream& operator>> <>(std::stringstream &ifs,
01736                                           GridData<Type,2> &gd);
01737 };
01738 
01739 template <class Type>
01740 std::ostream&  operator << (std::ostream& os, const GridData<Type,2> &gd) {
01741   if (&gd == (GridData<Type,2> *)0) return os;
01742 
01743   os << "BBox: " << gd.bbox() << " " << "Extents: " << gd.extents() << " ";
01744   os << "Step: " << gd.stepsize() << " ";
01745   os << "Size: " << gd.size() << " " << "Bottom: " << gd.bottom() << " ";
01746   
01747   os << "\n";
01748   
01749   if (gd.bbox().empty()) return os;
01750   
01751   const Coords &l = gd.lower();
01752   const Coords &u = gd.upper();
01753   const Coords &step = gd.stepsize();
01754   for (register int j=l(1);j<=u(1);j+=step(1)) {
01755     for (register int i=l(0);i<=u(0);i+=step(0)) 
01756       { os << "[" << i << "," << j << "]=" << gd(i, j) << " "; }
01757     os << "\n"; }
01758 
01759   return os;
01760 }
01761 
01762 template <class Type>
01763 std::ofstream&  operator << (std::ofstream& ofs, const GridData<Type,2> &gd) {
01764   if (&gd == (GridData<Type,2> *)0) return ofs;
01765 
01766   ofs.write((char*)&gd.bbox(),sizeof(BBox));
01767   ofs.write((char*)gd.data(),gd.size()*sizeof(Type));
01768 
01769   return ofs;
01770 }
01771 
01772 template <class Type>
01773 std::ifstream&  operator >> (std::ifstream& ifs, GridData<Type,2> &gd) {
01774   if (&gd == (GridData<Type,2> *)0) return ifs;
01775 
01776   BBox bb;
01777   ifs.read((char*)&bb,sizeof(BBox));
01778   
01779   if (!gd.ok_to_index())
01780     { gd.allocate(bb); }
01781   
01782   if (bb == gd.bbox()) {
01783     ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
01784   }
01785   else {
01786     GridData<Type,2> tmpgd(bb);
01787     ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
01788     gd.copy(tmpgd);  
01789   } 
01790   
01791   return ifs;
01792 }
01793 
01794 template <class Type>
01795 std::stringstream& operator<<(std::stringstream& ofs, 
01796                               const GridData<Type,2> &gd) {
01797   if (&gd == (GridData<Type,2> *)0) return ofs;
01798 
01799   ofs.write((char*)&gd.bbox(),sizeof(BBox));
01800   ofs.write((char*)gd.data(),gd.size()*sizeof(Type));
01801   
01802   return ofs;
01803 }
01804 
01805 template <class Type>
01806 std::stringstream& operator>>(std::stringstream& ifs, 
01807                               GridData<Type,2> &gd) {
01808   if (&gd == (GridData<Type,2> *)0) return ifs;
01809 
01810   BBox bb;
01811   ifs.read((char*)&bb,sizeof(BBox));
01812   
01813   if (!gd.ok_to_index())
01814     { gd.allocate(bb); }
01815   
01816   if (bb == gd.bbox()) {
01817     ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
01818   }
01819   else {
01820     GridData<Type,2> tmpgd(bb);
01821     ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
01822     gd.copy(tmpgd);  
01823   } 
01824   
01825   return ifs;
01826 }
01827 
01828 /*********************************************************************/
01829 template <class Type>
01830 inline Coords lower(GridData<Type,2> &gd) { return(gd.lower()); }
01831 template <class Type>
01832 inline Coords upper(GridData<Type,2> &gd) { return(gd.upper()); }
01833 template <class Type>
01834 inline Coords extents(GridData<Type,2> &gd) { return(gd.extents()); }
01835 template <class Type>
01836 inline int lower(GridData<Type,2> &gd, const int i) { return(gd.lower(i)); }
01837 template <class Type>
01838 inline int upper(GridData<Type,2> &gd, const int i) { return(gd.upper(i)); }
01839 template <class Type>
01840 inline int extents(GridData<Type,2> &gd, const int i) { return(gd.extents(i)); }
01841 template <class Type>
01842 inline BBox bbox(GridData<Type,2> &gd) { return(gd.bbox()); }
01843 template <class Type>
01844 inline int size(GridData<Type,2> &gd) { return(gd.size()); }
01845 template <class Type>
01846 inline Coords stepsize(GridData<Type,2> &gd) { return(gd.stepsize()); }
01847   
01848 #endif
Generated on Thu Jun 30 02:18:59 2016 for AMROC's Hierachical Data Structures by 1.6.3