GridData3.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef _included_GridData_3_h
#define _included_GridData_3_h

#include "DAGHParams.h"

#ifdef DEBUG_PRINT_GD_MEMORY
#include "DAGHMemoryTrace.h"
#endif

#include "BBox.h"
#include "Coords.h"
#include "PackedGridDataBucket.h"

#ifdef DEBUG_PRINT
#include "CommServer.h"
#endif

#include "IndexGridData3.h"
#include "generic.h"

#ifndef gd_OperateRegion
#define gd_OperateRegion(op)    name2(gd_OperateRegion,op)
#endif

// #include "GDIterator.h"

#include <iosfwd>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cassert>

template <class Type,int dim> class GridData;
template <class Type>
std::ostream& operator<<(std::ostream&, const GridData<Type,3>&);
template <class Type>
std::ofstream& operator<<(std::ofstream&, const GridData<Type,3>&);
template <class Type>
std::ifstream& operator>>(std::ifstream&, GridData<Type,3>&);
template <class Type>
std::stringstream& operator<<(std::stringstream&, const GridData<Type,3>&);
template <class Type>
std::stringstream& operator>>(std::stringstream&, GridData<Type,3>&);

template <class Type> class GridData<Type,3> {
  BBox _bbox;
  Coords _extents;
  Coords _step;
  int _size;
  int _bottom; 
  Type *_data;

public:
  GridData(void) : 
    _bbox(), _extents(0,0), _step(3,1),
    _size(0), _bottom(0), _data((Type *) 0) {}

  GridData(const BBox &bb) :
    _bbox(bb), _extents(bb.extents()), _step(bb.stepsize()),
    _size(bb.size()), _bottom(bb.bottom()),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
  }

  GridData(const int i, const int j, const int k,
           const int ii, const int jj, const int kk) : 
    _bbox(3,i,j,k,ii,jj,kk,1), _extents(_bbox.extents()), _step(_bbox.stepsize()),
    _size(_bbox.size()), _bottom(_bbox.bottom()),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
  }

   GridData(const int i, const int j, const int k,
            const int ii, const int jj, const int kk, const int s) :
     _bbox(3,i,j,k,ii,jj,kk,s), _extents(_bbox.extents()), _step(_bbox.stepsize()),
     _size(_bbox.size()), _bottom(_bbox.bottom()),
     _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
   }

  GridData(const int i, const int j, const int k,
           const int ii, const int jj, const int kk, 
           const int s, const int ss, const int sss) :
    _bbox(3,i,j,k,ii,jj,kk,s,ss,sss),
    _extents(_bbox.extents()), _step(_bbox.stepsize()),
    _size(_bbox.size()), _bottom(_bbox.bottom()),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
   }

  GridData(const BBox &bb, Type *databuf) : 
    _bbox(bb), _extents(bb.extents()), _step(bb.stepsize()),
    _size(bb.size()), _bottom(bb.bottom()), _data(databuf) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
   }

  GridData(GridDataBucket<Type> &gdbkt) : 
    _bbox(gdbkt.bbox()), _extents((gdbkt.bbox()).extents()),
    _step((gdbkt.bbox()).stepsize()), _size((gdbkt.bbox()).size()),
    _bottom((gdbkt.bbox()).bottom()), _data(gdbkt.data()) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
   }

  GridData(GridDataBucket<Type> &gdbkt, const int n) : 
    _bbox(gdbkt.bbox(n)), _extents((gdbkt.bbox(n)).extents()),
    _step((gdbkt.bbox(n)).stepsize()), _size((gdbkt.bbox(n)).size()),
    _bottom((gdbkt.bbox(n)).bottom()), _data(gdbkt.data(n)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
   }

  /* A "pseudo" copy-constructor */
  GridData (const GridData<Type,3> &other) : 
    _bbox(other._bbox), _extents(other._extents),
    _step(other._step), _size(other._size),
    _bottom(other._bottom), _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
  }

  inline ~GridData(void) { 
    if (_data) {
      delete [] _data; 
#ifdef DEBUG_PRINT_GD_MEMORY
      DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
    }
  }

  void allocate(const BBox &bb) {
    _bbox = bb;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _size = _bbox.size();
    _bottom = _bbox.bottom();
    _data = _size ? new Type[_size] : (Type *) 0;
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size);
#endif
  }

  void allocate(const BBox &bb, Type *databuf) {
    _bbox = bb;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _size = _bbox.size();
    _bottom = _bbox.bottom();
#ifdef DEBUG_PRINT_GD
    assert(_data==0);
#endif
    _data = databuf;
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
  }

  inline void allocate(Type *databuf) {
    if (_data) {
      delete [] _data; _data = (Type*)0;
#ifdef DEBUG_PRINT_GD_MEMORY
      DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
    }
    _data = databuf;
  }

  inline void deallocate() {
    if (_data) {
      delete [] _data; _data = (Type*)0;
#ifdef DEBUG_PRINT_GD_MEMORY
      DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
    }
  }
  inline void deallocate(Type*& databuf) 
  { databuf = _data; _data = (Type*) 0; }
   
  void* databuffer() { return ((void*) _data); }

  /* Query funtions */
  
  inline const Coords& lower() const { return(_bbox.lower()); }
  inline const Coords& upper() const { return(_bbox.upper()); }
  inline const Coords& extents() const { return(_extents); }
  inline const Coords& stepsize() const { return(_step); }
  inline Coords lower() { return(_bbox.lower()); }
  inline Coords upper() { return(_bbox.upper()); }
  inline Coords extents() { return(_extents); }
  inline Coords stepsize() { return(_step); }
  inline int bottom() const { return(_bottom); }
  inline int lower(const int i) const { return(_bbox.lower(i)); }
  inline int upper(const int i) const { return(_bbox.upper(i)); }
  inline int extents(const int i) const { return(_extents(i)); }
  inline int stepsize(const int i) const { return(_step(i)); }
  inline const BBox& bbox() const { return(_bbox); }

  inline int ok_to_index() { return(_data != (Type *) 0); }
  inline int size() const { return(_size); }
  
  inline int idx(const int i, const int j, const int k) const
  { return( _bottom+(i/_step(0)+_extents(0)*(j/_step(1)+_extents(1)*k/_step(2))) ); }
  
  inline const Type& operator () (const int i, const int j, const int k) const
  { assert(idx(i,j,k) >= 0 && idx(i,j,k) < _size); return(_data[idx(i,j,k)]); }
  inline Type& operator () (const int i,  const int j,  const int k)
  { assert(idx(i,j,k) >= 0 && idx(i,j,k) < _size); return(_data[idx(i,j,k)]); }

  inline const Type& operator () (const Coords& c) const
  { assert(idx(c(0),c(1),c(2)) >= 0 && idx(c(0),c(1),c(2)) < _size);
    return(_data[idx(c(0),c(1),c(2))]); }
  inline Type& operator () (const Coords& c)
  { assert(idx(c(0),c(1),c(2)) >= 0 && idx(c(0),c(1),c(2)) < _size);
    return(_data[idx(c(0),c(1),c(2))]); }

  inline const Type *ptr(const int i,  const int j,  const int k) const
  { assert(idx(i,j,k) >= 0 && idx(i,j,k) < _size); return(_data + idx(i,j,k)); }
  inline Type *ptr(const int i,  const int j,  const int k)
  { assert(idx(i,j,k) >= 0 && idx(i,j,k) < _size); return(_data + idx(i,j,k)); }

  inline const Type *ptr(const Coords& c) const
  { assert(idx(c(0),c(1),c(2)) >= 0 && idx(c(0),c(1),c(2)) < _size);
    return(_data + idx(c(0),c(1),c(2))); }
  inline Type *ptr(const Coords& c)
  { assert(idx(c(0),c(1),c(2)) >= 0 && idx(c(0),c(1),c(2)) < _size);
    return(_data + idx(c(0),c(1),c(2))); }

  inline const Type *data() const { return(_data); }
  inline Type *data() { return(_data); }

  void fill(const Type &val)
  { register int i; for (i = 0; i < _size; i++) _data[i] = val; }

  /* Define the copy functions for communication between objects */
  void copy(const GridData<Type,3> &gd) {
    if (&gd != this) {
      BBox intersection = _bbox * gd._bbox;
      if (!intersection.empty()) { 
        Coords max_step = _step.getmax(gd._step);
        BBox to(intersection); to.setstepsize(_step); 
        BBox from(intersection); from.setstepsize(gd._step);
        gd_CopyRegion(gd, to, from, max_step); 
      }
    }
  }

  void copy(const GridData<Type,3> &gd, const BBox &where) {
    if (&gd != this) {
      BBox intersection = _bbox * gd._bbox * where;
      if (!intersection.empty()) { 
        Coords max_step = _step.getmax(gd._step);
        BBox to(intersection); to.setstepsize(_step); 
        BBox from(intersection); from.setstepsize(gd._step);
        gd_CopyRegion(gd, to, from, max_step); 
      }
    }
  }

  void copy(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_CopyRegion(gd, newto, newfrom, max_step);
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt) {
    BBox const &gdbktbb = gdbkt.bbox();
    BBox intersection = _bbox * gdbktbb;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const BBox &where) {
    BBox const &gdbktbb = gdbkt.bbox();
    BBox intersection = _bbox * gdbktbb * where;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const BBox &to, 
            const BBox &from) {
    BBox const &gdbktbb = gdbkt.bbox();
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      newto.setstepsize(_step);
      newfrom.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, newto, newfrom, max_step);
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const int n) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    BBox intersection = _bbox * gdbktbb;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &where) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    BBox intersection = _bbox * gdbktbb * where;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to, 
            const BBox &from) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      newto.setstepsize(_step);
      newfrom.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, newto, newfrom, max_step);
    }
  }

  /* Linear interpolation */
  void lin_interp(const GridData<Type,3> &gd1, double const frac1,
                  const GridData<Type,3> &gd2, double const frac2,
                  const BBox &where) {
    BBox intersection = _bbox * gd1._bbox * gd2._bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(gd1._step);
      max_step.max(gd2._step);

      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(src1, gd1._bbox, gd1._data, const Type);
      BeginFastIndex3(src2, gd2._bbox, gd2._data, const Type);
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, intersection, max_step)
        FastIndex3(dst,i,j,k) = (Type)  ((FastIndex3(src1,i,j,k) * frac1)
                                         + (FastIndex3(src2,i,j,k) * frac2));
      end_for

      EndFastIndex3(dst);
      EndFastIndex3(src1);
      EndFastIndex3(src2);
    }
  }

  inline void lin_interp(const GridData<Type,3> &gd1, double const frac1,
                         const GridData<Type,3> &gd2, double const frac2)
  { GridData<Type,3>::lin_interp(gd1,frac1,gd2,frac2,_bbox); }

  /* first moment (i think) sum(q(i,j,k)*i) */
  Type moment1(const int axis, const BBox &where) {
    Type m1 = (Type) 0;
    
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,3> &dst = *this;

      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      if (axis == DAGH_X) {
        for_3(i, j, k, intersection, max_step)
          m1 += FastIndex3(dst,i,j,k) * i;
        end_for
      }
      else if (axis == DAGH_Y) {
        for_3(i, j, k, intersection, max_step)
          m1 += FastIndex3(dst,i,j,k) * j;
        end_for
      }
      else if (axis == DAGH_Z) {
        for_3(i, j, k, intersection, max_step)
          m1 += FastIndex3(dst,i,j,k) * k;
        end_for
      }

      EndFastIndex3(dst);
    }
    return (m1);
  }

  inline Type moment1(const int axis)
  { return (GridData<Type,3>::moment1(axis,_bbox)); }

  /* sum of squares - for calculating L2-norm */
  double sumsqrd(const BBox &where) {
    double s = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);
      
      for_3(i, j, k, intersection, max_step)
        s += Abs(FastIndex3(dst,i,j,k)) * Abs(FastIndex3(dst,i,j,k));
      end_for

      EndFastIndex3(dst);
    }
    return s;
  }

  inline double sumsqrd()
  { return (GridData<Type,3>::sumsqrd(_bbox)); }

  /* sum of absolute values - for calculating L1-norm */
  double sumabs(const BBox &where) {
    double s = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);

      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, intersection, max_step)
        s += Abs(FastIndex3(dst,i,j,k));
      end_for

      EndFastIndex3(dst);

    }
    return s;
  }

  inline double sumabs()
  { return (GridData<Type,3>::sumabs(_bbox)); }

  /* Maximum of absolute values - for calculating Loo-norm */
  double maxabs(const BBox &where) {
    double m = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,3> &dst = *this;

      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, intersection, max_step)
        double t = Abs(FastIndex3(dst,i,j,k));
        m = Max(m,t);
      end_for

      EndFastIndex3(dst);

    }
    return m;
  }
  
  inline double maxabs()
  { return (GridData<Type,3>::maxabs(_bbox)); }

  /* operator == */
  BBox is_eq(const Type &val, const BBox &where) {
     BBox intersection = _bbox * where;
     if (!intersection.empty()) {
       Coords max_step = _step.getmax(where.stepsize());
       BBox to(intersection); to.setstepsize(_step);
       BBox bb(3,max_step);
       gd_OperateRegion(is_eq)(val, to, max_step, bb);
       return (bb);
     }
     else 
       { BBox bb(3,1); return (bb); }
  }

  inline BBox is_eq(const Type &val) 
  { return (GridData<Type,3>::is_eq(val, _bbox)); }

  BBox is_eq(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_eq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_eq(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_eq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_eq(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_eq(gd, where, where)); }

  /* operator != */
  BBox is_neq(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_neq)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_neq(const Type &val) 
  { return (GridData<Type,3>::is_neq(val, _bbox)); }

  BBox is_neq(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_neq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_neq(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_neq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_neq(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_neq(gd, where, where)); }

  /* operator > */
  BBox is_gt(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_gt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_gt(const Type &val) 
  { return (GridData<Type,3>::is_gt(val, _bbox)); }
  
  BBox is_gt(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_gt)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_gt(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_gt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_gt(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_gt(gd, where, where)); }

  /* operator >= */
  BBox is_ge(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_ge)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_ge(const Type &val) 
  { return (GridData<Type,3>::is_ge(val, _bbox)); }

  BBox is_ge(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_ge)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_ge(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_ge(gd, gd._bbox, gd._bbox)); }
  inline BBox is_ge(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_ge(gd, where, where)); }

  /* operator < */
  BBox is_lt(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_lt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_lt(const Type &val) 
  { return (GridData<Type,3>::is_lt(val, _bbox)); }

  BBox is_lt(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
   const Coords  toshift = from.lower() - to.lower();
   BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
   if (!newfrom.empty()) {
     BBox newto = shiftabs(newfrom, -toshift);
     Coords max_step = _step.getmax(gd._step);
     newto.setstepsize(_step);
     newfrom.setstepsize(gd._step);
     BBox bb(3,max_step);
     gd_OperateRegion(is_lt)(gd, newto, newfrom, max_step, bb);
     return (bb);
   }
   else 
     { BBox bb(3,1); return (bb); }
  }

  inline BBox is_lt(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_lt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_lt(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_lt(gd, where, where)); }

  /* operator <= */
  BBox is_le(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_le)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_le(const Type &val) 
  { return (GridData<Type,3>::is_le(val, _bbox)); }

  BBox is_le(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(3,max_step);
      gd_OperateRegion(is_le)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(3,1); return (bb); }
  }

  inline BBox is_le(const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_le(gd, gd._bbox, gd._bbox)); }
  inline BBox is_le(const GridData<Type,3> &gd, const BBox &where)
  { return (GridData<Type,3>::is_le(gd, where, where)); }

  /* operator = */
  void equals(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(equal)(val, to, max_step);
    }
  }

  inline void equals(const Type &val) 
  { GridData<Type,3>::equals(val, _bbox); }

  void equals(const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      // Coords max_step = max(_step, gd._step);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(equal)(gd, newto, newfrom, max_step);
    }
  }

  inline void equals(const GridData<Type,3> &gd)
  { GridData<Type,3>::equals(gd, gd._bbox, gd._bbox); }
  inline void equals(const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::equals(gd, where, where); }

  /* operator += */
  void plus (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(plus)(val, to, max_step);
    }
  }

  inline void plus (const Type &val) 
  { GridData<Type,3>::plus(val, _bbox); }

  void plus (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(plus)(gd, newto, newfrom, max_step);
    }
  }

  inline void plus (const GridData<Type,3> &gd)
  { GridData<Type,3>::plus(gd, gd._bbox, gd._bbox); }
  inline void plus (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::plus(gd, where, where); }

  /* operator -= */
  void minus (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(minus)(val, to, max_step);
    }
  }

  inline void minus (const Type &val) 
  { GridData<Type,3>::minus(val, _bbox); }

  void minus (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(minus)(gd, newto, newfrom, max_step);
    }
  }

  inline void minus (const GridData<Type,3> &gd)
  { GridData<Type,3>::minus(gd, gd._bbox, gd._bbox); }
  inline void minus (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::minus(gd, where, where); }

  /* operator *= */
  void multiply(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(mult)(val, to, max_step);
    }
  }

  inline void multiply(const Type &val) 
  { GridData<Type,3>::multiply(val, _bbox); }
  
  void multiply (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(mult)(gd, newto, newfrom, max_step);
    }
  }

  inline void multiply (const GridData<Type,3> &gd)
  { GridData<Type,3>::multiply(gd, gd._bbox, gd._bbox); }
  inline void multiply (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::multiply(gd, where, where); }

  /* operator /= */
  void divide (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(div)(val, to, max_step);
    }
  }

  inline void divide (const Type &val) 
  { GridData<Type,3>::divide(val, _bbox); }

  void divide (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(div)(gd, newto, newfrom, max_step);
    }
  }

  inline void divide (const GridData<Type,3> &gd)
  { GridData<Type,3>::divide(gd, gd._bbox, gd._bbox); }
  inline void divide (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::divide(gd, where, where); }
  
  /* point-wise minimum */
  void minimum (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(min)(val, to, max_step);
    }
  }

  inline void minimum (const Type &val) 
  { GridData<Type,3>::minimum(val, _bbox); }

  void minimum (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(min)(gd, newto, newfrom, max_step);
    }
  }

  inline void minimum (const GridData<Type,3> &gd)
  { GridData<Type,3>::minimum(gd, gd._bbox, gd._bbox); }
  inline void minimum (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::minimum(gd, where, where); }

  /* point-wise maximum */
  void maximum (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(max)(val, to, max_step);
    }
  }

  inline void maximum (const Type &val) 
  { GridData<Type,3>::maximum(val, _bbox); }

  void maximum (const GridData<Type,3> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(max)(gd, newto, newfrom, max_step);
    }
  }

  inline void maximum (const GridData<Type,3> &gd)
  { GridData<Type,3>::maximum(gd, gd._bbox, gd._bbox); }
  inline void maximum (const GridData<Type,3> &gd, const BBox &where)
  { GridData<Type,3>::maximum(gd, where, where); }

  /**********************************************************************/
  inline BBox operator == (const Type &val)
  { return (GridData<Type,3>::is_eq(val, _bbox)); }
  inline BBox operator == (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_eq(gd, gd._bbox, gd._bbox)); }

  inline BBox operator != (const Type &val)
  { return (GridData<Type,3>::is_neq(val, _bbox)); }
  inline BBox operator != (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_neq(gd, gd._bbox, gd._bbox)); }

  inline BBox operator > (const Type &val)
  { return (GridData<Type,3>::is_gt(val, _bbox)); }
  inline BBox operator > (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_gt(gd, gd._bbox, gd._bbox)); }

  inline BBox operator >= (const Type &val)
  { return (GridData<Type,3>::is_ge(val, _bbox)); }
  inline BBox operator >= (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_ge(gd, gd._bbox, gd._bbox)); }

  inline BBox operator < (const Type &val)
  { return (GridData<Type,3>::is_lt(val, _bbox)); }
  inline BBox operator < (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_lt(gd, gd._bbox, gd._bbox)); }
  
  inline BBox operator <= (const Type &val)
  { return (GridData<Type,3>::is_le(val, _bbox)); }
  inline BBox operator <= (const GridData<Type,3> &gd)
  { return (GridData<Type,3>::is_le(gd, gd._bbox, gd._bbox)); }

  inline void operator = (const Type &val)
  { GridData<Type,3>::equals(val, _bbox); }
  inline void operator = (const GridData<Type,3> &gd)
  { GridData<Type,3>::equals(gd, gd._bbox, gd._bbox); }

  inline void operator += (const Type &val)
  { GridData<Type,3>::plus(val, _bbox); }
  inline void operator += (const GridData<Type,3> &gd)
  { GridData<Type,3>::plus(gd, gd._bbox, gd._bbox); }

  inline void operator -= (const Type &val)
  { GridData<Type,3>::minus(val, _bbox); }
  inline void operator -= (const GridData<Type,3> &gd)
  { GridData<Type,3>::minus(gd, gd._bbox, gd._bbox); }
  
  inline void operator *= (const Type &val)
  { GridData<Type,3>::multiply(val, _bbox); }
  inline void operator *= (const GridData<Type,3> &gd)
  { GridData<Type,3>::multiply(gd, gd._bbox, gd._bbox); }

  inline void operator /= (const Type &val)
  { GridData<Type,3>::divide(val, _bbox); }
  inline void operator /= (const GridData<Type,3> &gd)
  { GridData<Type,3>::divide(gd, gd._bbox, gd._bbox); }
  /**********************************************************************/

  /**********************************************************************/
  // Reduction Operations 
  /**********************************************************************/

  /****************************** maxval ***************************************/
  Type maxval (const BBox &where) {
    Type max_val = (Type) 0;
    short flag = DAGHTrue;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);
      
      for_3(i, j, k, intersection, max_step)
        if (flag == DAGHTrue)
          { max_val = FastIndex3(dst,i,j,k); flag = DAGHFalse; }
        else 
          max_val = Max(max_val,FastIndex3(dst,i,j,k));
      end_for

      EndFastIndex3(dst);

    }
    return (max_val);
  }

  inline Type maxval (void)
  { return (GridData<Type,3>::maxval(_bbox)); }
  /************************************************************************/

  /****************************** minval ***************************************/
  Type minval (const BBox &where) {
    Type min_val = (Type) 0;
    short flag = DAGHTrue;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);
      
      for_3(i, j, k, intersection, max_step)
        if (flag == DAGHTrue)
          { min_val = FastIndex3(dst,i,j,k); flag = DAGHFalse; }
        else 
          min_val = Min(min_val,FastIndex3(dst,i,j,k));
      end_for

      EndFastIndex3(dst);

    }
    return (min_val);
  }

  inline Type minval (void)
  { return (GridData<Type,3>::minval(_bbox)); }
  /************************************************************************/

  /****************************** sum ***************************************/
  Type sum (const BBox &where) {
    Type sum_val = (Type) 0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);

      GridData<Type,3> &dst = *this;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, intersection, max_step)
        sum_val += FastIndex3(dst,i,j,k);
      end_for

      EndFastIndex3(dst);

    }
    return (sum_val);
  }

  inline Type sum (void)
  { return (GridData<Type,3>::sum(_bbox)); }
  /************************************************************************/

  /****************************** product ***************************************/
  Type product (const BBox &where) {
    Type prod_val = (Type) 1;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
       
      GridData<Type,3> &dst = *this;

      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, intersection, max_step)
        prod_val *= FastIndex3(dst,i,j,k);
      end_for

      EndFastIndex3(dst);

    }
    return (prod_val);
  }

  inline Type product (void)
  { return (GridData<Type,3>::product(_bbox)); }
  /**********************************************************************/

private:
  void gd_CopyRegion(const GridData<Type,3> &src, const BBox &to,
                     const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    
    if (dst._bbox==src._bbox && to==from && to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) src._data, sizeof(Type)*to.size());
    else {
      const int di = from.lower(0)-to.lower(0);
      const int dj = from.lower(1)-to.lower(1);
      const int dk = from.lower(2)-to.lower(2);
      
      BeginFastIndex3(src, src._bbox, src._data, const Type);
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);

      for_3(i, j, k, to, step)
        FastIndex3(dst,i,j,k) = FastIndex3(src,i+di,j+dj,k+dk);
      end_for

      EndFastIndex3(dst);
      EndFastIndex3(src);
    }
  }

  void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const BBox &to,
                      const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;

    if (dst._bbox==gdbkt.bbox() && to==from && to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) gdbkt.data(), sizeof(Type)*to.size());
    else {
      const int di = from.lower(0)-to.lower(0);
      const int dj = from.lower(1)-to.lower(1);
      const int dk = from.lower(2)-to.lower(2);

      BeginFastIndex3(src, gdbkt.bbox(), gdbkt.data(), const Type);
      BeginFastIndex3(dst, dst.bbox(), dst.data(), Type);
      for_3(i, j, k, to, step)
        FastIndex3(dst,i,j,k) = FastIndex3(src,i+di,j+dj,k+dk);
      end_for
      EndFastIndex3(dst);
      EndFastIndex3(src);
    }
  }

  void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to,
                      const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    
    if (dst._bbox==gdbkt.bbox(n) && to==from && to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) gdbkt.data(n), sizeof(Type)*to.size());
    else {
      const int di = from.lower(0)-to.lower(0);
      const int dj = from.lower(1)-to.lower(1);
      const int dk = from.lower(2)-to.lower(2);
      
      BeginFastIndex3(src, gdbkt.bbox(n), gdbkt.data(n), const Type);
      BeginFastIndex3(dst, dst.bbox(), dst.data(), Type);
      for_3(i, j, k, to, step)
        FastIndex3(dst,i,j,k) = FastIndex3(src,i+di,j+dj,k+dk);
      end_for
      EndFastIndex3(dst);
      EndFastIndex3(src);
    }
  }

  /****************************** == ***************************************/
  void gd_OperateRegion(is_eq)(const Type &val, const BBox &to,
                               Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) == val) bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(is_eq)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);
    
    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) == FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** != ***************************************/
   void gd_OperateRegion(is_neq)(const Type &val, const BBox &to,
                                 Coords const &step, BBox &bb) {
     GridData<Type,3> &dst = *this;

     BeginFastIndex3(dst, dst._bbox, dst._data, Type);

     for_3(i, j, k, to, step)
       if (FastIndex3(dst,i,j,k) != val) bb += Coords(3,i,j,k);
     end_for

     EndFastIndex3(dst);
   }

  void gd_OperateRegion(is_neq)(const GridData<Type,3> &src, const BBox &to,
                                const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);

    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) != FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** > ***************************************/
  void gd_OperateRegion(is_gt)(const Type &val, const BBox &to,
                               Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) > val) bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(is_gt)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);
    
    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) > FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** >= ***************************************/
  void gd_OperateRegion(is_ge)(const Type &val, const BBox &to,
                               Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) >= val) bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(is_ge)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) >= FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** < ***************************************/
  void gd_OperateRegion(is_lt)(const Type &val, const BBox &to,
                               Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) < val) bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(is_lt)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);

    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) < FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** <= ***************************************/
  void gd_OperateRegion(is_le)(const Type &val, const BBox &to,
                               Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);
    
    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) <= val) bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(is_le)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step, BBox &bb) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      if (FastIndex3(dst,i,j,k) <= FastIndex3(src,i+di,j+dj,k+dk))
        bb += Coords(3,i,j,k);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  
  /****************************** = ***************************************/
  void gd_OperateRegion(equal)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);
    
    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = val;
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(equal)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = FastIndex3(src,i+di,j+dj,k+dk);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** + ***************************************/
  void gd_OperateRegion(plus)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) += val;
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(plus)(const GridData<Type,3> &src, const BBox &to,
                              const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);

    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) += FastIndex3(src,i+di,j+dj,k+dk);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** - ***************************************/
  void gd_OperateRegion(minus)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) -= val;
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(minus)(const GridData<Type,3> &src, const BBox &to,
                               const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) -= FastIndex3(src,i+di,j+dj,k+dk);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** * ***************************************/
  void gd_OperateRegion(mult)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) *= val;
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(mult)(const GridData<Type,3> &src, const BBox &to,
                              const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);
    
    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) *= FastIndex3(src,i+di,j+dj,k+dk);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** / ***************************************/
  void gd_OperateRegion(div)(const Type &val, const BBox &to, Coords const &step) {
    assert (val != (Type)0);
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) /= val;
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(div)(const GridData<Type,3> &src, const BBox &to,
                             const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) /= FastIndex3(src,i+di,j+dj,k+dk);
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** min ***************************************/
  void gd_OperateRegion(min)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = Min(val, FastIndex3(dst,i,j,k));
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(min)(const GridData<Type,3> &src, const BBox &to,
                             const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = Min(FastIndex3(dst,i,j,k), FastIndex3(src,i+di,j+dj,k+dk));
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

  /****************************** max ***************************************/
  void gd_OperateRegion(max)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,3> &dst = *this;

    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = Max(val, FastIndex3(dst,i,j,k));
    end_for

    EndFastIndex3(dst);
  }

  void gd_OperateRegion(max)(const GridData<Type,3> &src, const BBox &to,
                             const BBox &from, Coords const &step) {
    GridData<Type,3> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    const int dj = from.lower(1)-to.lower(1);
    const int dk = from.lower(2)-to.lower(2);
    
    BeginFastIndex3(src, src._bbox, src._data, const Type);
    BeginFastIndex3(dst, dst._bbox, dst._data, Type);

    for_3(i, j, k, to, step)
      FastIndex3(dst,i,j,k) = Max(FastIndex3(dst,i,j,k), FastIndex3(src,i+di,j+dj,k+dk));
    end_for

    EndFastIndex3(dst);
    EndFastIndex3(src);
  }
  /************************************************************************/

public:
  /* Perform efficient buffer packing/unpacking and data movement */
  void PackRegion(Type *sendbuf, const BBox &from) const {
    GridData<Type,3> const &src = *this;

    if (from==src._bbox) 
      std::memcpy((void *) sendbuf, (void *) src._data, sizeof(Type)*from.size());
    else {
      register int p = 0;
      
      BeginFastIndex3(src, src._bbox, src._data, const Type);

      for_3(i, j, k, from, _step)
        sendbuf[p++] = FastIndex3(src,i,j,k);
      end_for
        
      EndFastIndex3(src);
    }
  }

  void UnPackRegion(const Type *recvbuf, const BBox &to) {
    GridData<Type,3> &dst = *this;
    
    if (to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) recvbuf, sizeof(Type)*to.size());
    else {
      register int p = 0;
      
      BeginFastIndex3(dst, dst._bbox, dst._data, Type);
      
      for_3(i, j, k, to, _step)
         FastIndex3(dst,i,j,k) = recvbuf[p++];
      end_for
        
        EndFastIndex3(dst);
    }
  }

  friend std::ostream& operator<< <>(std::ostream &os, 
                                     const GridData<Type,3> &gd);
  friend std::ofstream& operator<< <>(std::ofstream &ofs, 
                                      const GridData<Type,3> &gd);
  friend std::ifstream& operator>> <>(std::ifstream &ifs, 
                                      GridData<Type,3> &gd);  
  friend std::stringstream& operator<< <>(std::stringstream &ofs, 
                                          const GridData<Type,3> &gd);
  friend std::stringstream& operator>> <>(std::stringstream &ifs, 
                                          GridData<Type,3> &gd);
};

template <class Type>
std::ostream&  operator << (std::ostream& os, const GridData<Type,3> &gd) {
  if (&gd == (GridData<Type,3> *)0) return os;
  
  os << "BBox: " << gd.bbox() << " " << "Extents: " << gd.extents() << " ";
  os << "Step: " << gd.stepsize() << " ";
  os << "Size: " << gd.size() << " " << "Bottom: " << gd.bottom() << " ";
  
  os << "\n";
  
  if (gd.bbox().empty()) return os;

  const Coords &l = gd.lower();
  const Coords &u = gd.upper();
  const Coords &step = gd.stepsize();
  for (register int k=l(2);k<=u(2);k+=step(2)) {
    for (register int j=l(1);j<=u(1);j+=step(1)) {
      for (register int i=l(0);i<=u(0);i+=step(0)) 
        { os << "[" << i << "," << j << "," << k << "]=" << gd(i, j, k) << " "; }
      os << "\n"; } os << "\n";}

  return os;
}

template <class Type>
std::ofstream&  operator << (std::ofstream& ofs, const GridData<Type,3> &gd) {
  if (&gd == (GridData<Type,3> *)0) return ofs;
  
  ofs.write((char*)&gd.bbox(),sizeof(BBox));
  ofs.write((char*)gd.data(),gd.size()*sizeof(Type));
  
  return ofs;
}

template <class Type>
std::ifstream&  operator >> (std::ifstream& ifs, GridData<Type,3> &gd) {
  if (&gd == (GridData<Type,3> *)0) return ifs;
  
  BBox bb;
  ifs.read((char*)&bb,sizeof(BBox));
  
  if (!gd.ok_to_index())
    { gd.allocate(bb); }
  
  if (bb == gd.bbox()) {
    ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
  }
  else {
    GridData<Type,3> tmpgd(bb);
    ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
    gd.copy(tmpgd);  
  } 
  
  return ifs;
}

template <class Type>
std::stringstream& operator<<(std::stringstream& ofs, const GridData<Type,3> &gd) {
  if (&gd == (GridData<Type,3> *)0) return ofs;
  
  ofs.write((char*)&gd.bbox(),sizeof(BBox));
  ofs.write((char*)gd.data(),gd.size()*sizeof(Type));

  return ofs;
}

template <class Type>
std::stringstream& operator>>(std::stringstream& ifs, GridData<Type,3> &gd) {
  if (&gd == (GridData<Type,3> *)0) return ifs;
  
  BBox bb;
  ifs.read((char*)&bb,sizeof(BBox));
  
  if (!gd.ok_to_index())
    { gd.allocate(bb); }
  
  if (bb == gd.bbox()) {
    ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
  }
  else {
    GridData<Type,3> tmpgd(bb);
    ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
    gd.copy(tmpgd);  
  } 

  return ifs;
}

/*********************************************************************/
template <class Type>
inline Coords lower(GridData<Type,3> &gd) { return(gd.lower()); }
template <class Type>
inline Coords upper(GridData<Type,3> &gd) { return(gd.upper()); }
template <class Type>
inline Coords extents(GridData<Type,3> &gd) { return(gd.extents()); }
template <class Type>
inline int lower(GridData<Type,3> &gd, const int i) { return(gd.lower(i)); }
template <class Type>
inline int upper(GridData<Type,3> &gd, const int i) { return(gd.upper(i)); }
template <class Type>
inline int extents(GridData<Type,3> &gd, const int i) { return(gd.extents(i)); }
template <class Type>
inline BBox bbox(GridData<Type,3> &gd) { return(gd.bbox()); }
template <class Type>
inline int size(GridData<Type,3> &gd) { return(gd.size()); }
template <class Type>
inline Coords stepsize(GridData<Type,3> &gd) { return(gd.stepsize()); }
  
#endif

---