GridData2.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef _included_GridData_2_h
#define _included_GridData_2_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 "IndexGridData2.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,2>&);
template <class Type>
std::ofstream& operator<<(std::ofstream&, const GridData<Type,2>&);
template <class Type>
std::ifstream& operator>>(std::ifstream&, GridData<Type,2>&);
template <class Type>
std::stringstream& operator<<(std::stringstream&, const GridData<Type,2>&);
template <class Type>
std::stringstream& operator>>(std::stringstream&, GridData<Type,2>&);

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

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

  GridData(const BBox &bb) : _bbox(bb) {
    _bbox.rank = 2; 
    _bbox.lower().rank = 2;
    _bbox.upper().rank = 2;
    _bbox.stepsize().rank = 2;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
    _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 ii, const int jj) : 
    _bbox(2,i,j,ii,jj,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 ii, const int jj, const int s) : 
    _bbox(2,i,j,ii,jj,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 ii, const int jj, 
            const int s, const int ss) : 
     _bbox(2,i,j,ii,jj,s,ss), _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), _data(databuf) {
     _bbox.rank = 2; 
     _bbox.lower().rank = 2;
     _bbox.upper().rank = 2;
     _bbox.stepsize().rank = 2;
     _extents = _bbox.extents();
     _step = _bbox.stepsize();
     _bottom = _bbox.bottom();
     _size = _bbox.size();
#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()), _data(gdbkt.data()) {
    _bbox.rank = 2; 
    _bbox.lower().rank = 2;
    _bbox.upper().rank = 2;
    _bbox.stepsize().rank = 2;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
#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)), _data(gdbkt.data(n)) {
     _bbox.rank = 2; 
     _bbox.lower().rank = 2;
     _bbox.upper().rank = 2;
     _bbox.stepsize().rank = 2;
     _extents = _bbox.extents();
     _step = _bbox.stepsize();
     _bottom = _bbox.bottom();
     _size = _bbox.size();
#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,2> &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;
     _bbox.rank = 2; 
     _bbox.lower().rank = 2;
     _bbox.upper().rank = 2;
     _bbox.stepsize().rank = 2;
     _extents = _bbox.extents();
     _step = _bbox.stepsize();
     _bottom = _bbox.bottom();
     _size = _bbox.size();
     _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;
    _bbox.rank = 2; 
    _bbox.lower().rank = 2;
    _bbox.upper().rank = 2;
    _bbox.stepsize().rank = 2;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
#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 */
  /* Define some inline functions for getting at region extents */
  
  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
  { return( _bottom+(i/_step(0)+_extents(0)*(j/_step(1))) ); }

  inline const Type& operator () (const int i, const int j) const
  { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data[idx(i,j)]); }
  inline Type& operator () (const int i,  const int j)
  { assert(idx(i,j) >= 0 && idx(i,j) < _size); return(_data[idx(i,j)]); }

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

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

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

  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,2> &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,2> &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,2> &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) {
    const BBox &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) {
    const BBox &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) {
    const BBox &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) {
    const BBox &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)  {
    const BBox &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) {
    const BBox &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,2> &gd1, double const frac1,
                  const GridData<Type,2> &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,2> &dst = *this;

      BeginFastIndex2(src1, gd1._bbox, gd1._data, const Type);
      BeginFastIndex2(src2, gd2._bbox, gd2._data, const Type);
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
      
      for_2(i, j, intersection, max_step)
        FastIndex2(dst,i,j) = (Type) ((FastIndex2(src1,i,j) * frac1)
                                      + (FastIndex2(src2,i,j) * frac2));
      end_for

      EndFastIndex2(dst);
      EndFastIndex2(src1);
      EndFastIndex2(src2);
    }
  }

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

  /* first moment (i think) sum(q(i,j)*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,2> &dst = *this;
      
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);

      if (axis == DAGH_X) {
        for_2(i, j, intersection, max_step)
          m1 += FastIndex2(dst,i,j) * i;
        end_for
      }
      else if (axis == DAGH_Y) {
        for_2(i, j, intersection, max_step)
          m1 += FastIndex2(dst,i,j) * j;
        end_for
      }

      EndFastIndex2(dst);
    }
    return (m1); 
  }

  inline Type moment1(const int axis)
  { return (GridData<Type,2>::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,2> &dst = *this;

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

      for_2(i, j, intersection, max_step)
        s += Abs(FastIndex2(dst,i,j)) * Abs(FastIndex2(dst,i,j));
      end_for

      EndFastIndex2(dst);

    }
    return s;
  }

  inline double sumsqrd()
  { return (GridData<Type,2>::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,2> &dst = *this;
      
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
      
      for_2(i, j, intersection, max_step)
        s += Abs(FastIndex2(dst,i,j));
      end_for

      EndFastIndex2(dst);

    }
    return s;
  }

  inline double sumabs()
  { return (GridData<Type,2>::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,2> &dst = *this;

      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
      
      for_2(i, j, intersection, max_step)
        double t = Abs(FastIndex2(dst,i,j));
        m = Max(m,t);
      end_for

      EndFastIndex2(dst);

    }
    return m;
  }

  inline double maxabs()
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_eq)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_eq(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_eq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_eq(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_eq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_eq(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_neq)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_neq(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_neq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_neq(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_neq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_neq(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_gt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_gt(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_gt)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_gt(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_gt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_gt(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_ge)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_ge(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_ge)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_ge(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_ge(gd, gd._bbox, gd._bbox)); }
  inline BBox is_ge(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_lt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_lt(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_lt)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_lt(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_lt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_lt(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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(2,max_step);
      gd_OperateRegion(is_le)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

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

  BBox is_le(const GridData<Type,2> &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(2,max_step);
      gd_OperateRegion(is_le)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(2,1); return (bb); }
  }

  inline BBox is_le(const GridData<Type,2> &gd)
  { return (GridData<Type,2>::is_le(gd, gd._bbox, gd._bbox)); }
  inline BBox is_le(const GridData<Type,2> &gd, const BBox &where)
  { return (GridData<Type,2>::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,2>::equals(val, _bbox); }

  void equals(const GridData<Type,2> &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(equal)(gd, newto, newfrom, max_step);
    }
  }
  
  inline void equals(const GridData<Type,2> &gd)
  { GridData<Type,2>::equals(gd, gd._bbox, gd._bbox); }
  inline void equals(const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::plus(val, _bbox); }

  void plus (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::plus(gd, gd._bbox, gd._bbox); }
  inline void plus (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::minus(val, _bbox); }

  void minus (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::minus(gd, gd._bbox, gd._bbox); }
  inline void minus (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::multiply(val, _bbox); }

  void multiply (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::multiply(gd, gd._bbox, gd._bbox); }
  inline void multiply (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::divide(val, _bbox); }

  void divide (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::divide(gd, gd._bbox, gd._bbox); }
  inline void divide (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::minimum(val, _bbox); }

  void minimum (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::minimum(gd, gd._bbox, gd._bbox); }
  inline void minimum (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::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,2>::maximum(val, _bbox); }

  void maximum (const GridData<Type,2> &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,2> &gd)
  { GridData<Type,2>::maximum(gd, gd._bbox, gd._bbox); }
  inline void maximum (const GridData<Type,2> &gd, const BBox &where)
  { GridData<Type,2>::maximum(gd, where, where); }

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

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

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

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

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

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

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

  inline void operator /= (const Type &val)
  { GridData<Type,2>::divide(val, _bbox); }
  inline void operator /= (const GridData<Type,2> &gd)
  { GridData<Type,2>::divide(gd, gd._bbox, gd._bbox); }
  /**********************************************************************/
  
  /**********************************************************************/
  // Reduction Operations
  /**********************************************************************/
  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,2> &dst = *this;
      
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
      
      for_2(i, j, intersection, max_step)
        if (flag == DAGHTrue) 
          { max_val = FastIndex2(dst,i,j); flag = DAGHFalse; }
        else 
           max_val = Max(max_val, FastIndex2(dst,i,j));
      end_for

      EndFastIndex2(dst);

    }
    return (max_val);
  }

  inline Type maxval (void)
  { return (GridData<Type,2>::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,2> &dst = *this;

      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
      
      for_2(i, j, intersection, max_step)
        if (flag == DAGHTrue) 
          { min_val = FastIndex2(dst,i,j); flag = DAGHFalse; }
        else 
          min_val = Min(min_val, FastIndex2(dst,i,j));
      end_for

      EndFastIndex2(dst);

    }
    return (min_val);
  }

  inline Type minval (void)
  { return (GridData<Type,2>::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,2> &dst = *this;
      
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);

      for_2(i, j, intersection, max_step)
        sum_val += FastIndex2(dst,i,j);
      end_for

      EndFastIndex2(dst);

    }
    return (sum_val);
  }

  inline Type sum (void)
  { return (GridData<Type,2>::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,2> &dst = *this;
      
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);

      for_2(i, j, intersection, max_step)
        prod_val *= FastIndex2(dst,i,j);
      end_for

      EndFastIndex2(dst);

    }
    return (prod_val);
  }

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

private:
  void gd_CopyRegion(const GridData<Type,2> &src, const BBox &to,
                     const BBox &from, Coords const &step) {
    GridData<Type,2> &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);

      BeginFastIndex2(src, src._bbox, src._data, const Type);
      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
     
      for_2(i, j, to, step)
        FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
      end_for

      EndFastIndex2(dst);
      EndFastIndex2(src);
    }
  }

  void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const BBox &to,
                      const BBox &from, Coords const &step) {
    GridData<Type,2> &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);

      BeginFastIndex2(src, gdbkt.bbox(), gdbkt.data(), const Type);
      BeginFastIndex2(dst, dst.bbox(), dst.data(), Type);
      for_2(i, j, to, step)
        FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
      end_for
      EndFastIndex2(dst);
      EndFastIndex2(src);
    }
  }

  void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to,
                      const BBox &from, Coords const &step) {
    GridData<Type,2> &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);

      BeginFastIndex2(src, gdbkt.bbox(n), gdbkt.data(n), const Type);
      BeginFastIndex2(dst, dst.bbox(), dst.data(), Type);
      for_2(i, j, to, step)
        FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
      end_for
      EndFastIndex2(dst);
      EndFastIndex2(src);
    }
  }

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) == val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) == FastIndex2(src,i+di,j+dj))
        bb += Coords(2,i,j);
    end_for

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) != val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) != FastIndex2(src,i+di,j+dj))
        bb += Coords(2,i,j);
    end_for

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) > val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) >= val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) >= FastIndex2(src,i+di,j+dj))
        bb += Coords(2,i,j);
    end_for

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) < val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) < FastIndex2(src,i+di,j+dj))
        bb += Coords(2,i,j);
    end_for

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

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

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) <= val) bb += Coords(2,i,j);
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      if (FastIndex2(dst,i,j) <= FastIndex2(src,i+di,j+dj))
        bb += Coords(2,i,j);
    end_for

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = val;
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = FastIndex2(src,i+di,j+dj);
    end_for

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

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) += val;
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) += FastIndex2(src,i+di,j+dj);
    end_for

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

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) -= val;
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) -= FastIndex2(src,i+di,j+dj);
    end_for

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

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

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

     for_2(i, j, to, step)
       FastIndex2(dst,i,j) *= val;
     end_for

     EndFastIndex2(dst);
   }

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

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

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) /= val;
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) /= FastIndex2(src,i+di,j+dj);
    end_for

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

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = Min(val, FastIndex2(dst,i,j));
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = Min(FastIndex2(dst,i,j), FastIndex2(src,i+di,j+dj));
    end_for

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

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

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = Max(val, FastIndex2(dst,i,j));
    end_for

    EndFastIndex2(dst);
  }

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

    for_2(i, j, to, step)
      FastIndex2(dst,i,j) = Max(FastIndex2(dst,i,j), FastIndex2(src,i+di,j+dj));
    end_for

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

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

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

      BeginFastIndex2(src, src._bbox, src._data, const Type);

      for_2(i, j, from, _step)
        sendbuf[p++] = FastIndex2(src,i,j);
      end_for

      EndFastIndex2(src);
    }
  }

  void UnPackRegion(const Type *recvbuf, const BBox &to) {
    GridData<Type,2> &dst = *this;

    if (to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) recvbuf, sizeof(Type)*to.size());
    else {
      register int p = 0;

      BeginFastIndex2(dst, dst._bbox, dst._data, Type);
     
      for_2(i, j, to, _step)
        FastIndex2(dst,i,j) = recvbuf[p++];
      end_for

      EndFastIndex2(dst);
    }
  }

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

template <class Type>
std::ostream&  operator << (std::ostream& os, const GridData<Type,2> &gd) {
  if (&gd == (GridData<Type,2> *)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 j=l(1);j<=u(1);j+=step(1)) {
    for (register int i=l(0);i<=u(0);i+=step(0)) 
      { os << "[" << i << "," << j << "]=" << gd(i, j) << " "; }
    os << "\n"; }

  return os;
}

template <class Type>
std::ofstream&  operator << (std::ofstream& ofs, const GridData<Type,2> &gd) {
  if (&gd == (GridData<Type,2> *)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,2> &gd) {
  if (&gd == (GridData<Type,2> *)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,2> 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,2> &gd) {
  if (&gd == (GridData<Type,2> *)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,2> &gd) {
  if (&gd == (GridData<Type,2> *)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,2> tmpgd(bb);
    ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
    gd.copy(tmpgd);  
  } 
  
  return ifs;
}

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

---