amroc/hds/GridData3.h

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