#ifndef GRID_H
#define GRID_H

// __________________________________________________________________________
// Include RCS Id.
// $Id: grid.h,v 1.3 2000/12/12 18:32:05 strahs Exp $

// _________________________________________________________________________
// Includes

#include <assert.h>
#include <fstream.h>
#include <stdio.h>

#include "vector.h"

// __________________________________________________________________________
// grid of elements of type T

template <class T>
class grid
{
  
 private:
  int nx, ny;
  T *sto;
  
  // Initialization Member Functions
  void initnew(int x, int y) { 
    assert(x >= 0 && y >= 0); 
    nx = x; ny = y; sto = (x == 0 || y == 0) ? 0 : new T[x * y];
  }
  void freegrid() { assert(invariant()); nx = ny = 0; delete [] sto; sto = 0; }
  
  // Consistency Checking Functions
  int inrange(int x, int y) const { 
    return (x >= 0 && y >= 0 && x < nx && y < ny && sto != 0); 
  }
  int invariant() const { 
    return (((nx == 0 || ny == 0) && sto == 0 && nx >= 0 && ny >= 0) || 
	    ((nx > 0 && ny > 0) && sto != 0));
  }
  
 public:
  // Constructors, Destructor, and Operator=
  grid() { initnew(0, 0); }
  grid(int x) { initnew(x, 1); }
  grid(int x, int y) { initnew(x, y); }
  grid(const grid &m) : nx(0), ny(0), sto(0) { *this = m; }
  grid<T>& operator=(const grid &m);
  void subgrid(grid& m, int x, int y, int xsize, int ysize) const;

  ~grid() { assert (invariant()); delete [] sto; }
  
  // Basic Member Access Functions
  int dimx() const { return nx; }
  int dimy() const { return ny; }
  T& operator()(int x) const { 
    assert(inrange(x, 0)); return sto[x]; 
  }
  T& operator()(int x, int y) const { 
    assert(inrange(x, y)); return sto[y * nx + x]; 
  }
  T& operator[](int x) const { return sto[x]; }
  void set(int x, int y, T val) { if (inrange(x, y)) (*this)(x, y) = val; }
  void set(double x, double y, T val) { set(int(x), int(y), val); }
  void set(int x, double y, T val) { set(x, int(y), val); }
  void set(double x, int y, T val) { set(int(x), y, val); }
  T get(int x, int y) const { return (inrange(x, y) ? (*this)(x, y) : 0); }
  T get(double x, double y) const { return get(int(x), int(y)); }
  T get(int x, double y) const { return get(x, int(y)); }
  T get(double x, int y) const { return get(int(x), y); }
  
  // Initialization Member Functions
  void init() { init(0, 0); }
  void init(int x) { init(x, 1); }
  void init(int x, int y) { 
    assert(invariant());
    if (x > 0 && y > 0 && x * y <= nx * ny && x * y * 2 >= nx * ny) 
    { nx = x; ny = y; }
    else 
    { freegrid() ; initnew(x, y); } 
    assert(invariant()); 
  }
  void clear(const T& val = 0)
  { assert(invariant()); for (int i = 0; i < nx * ny; ++i) sto[i] = val; }
  
  // I/O Functions
  int loadpgm(const char* pgmname);
  int size() const;
  void dump(int max = -1) const;
  void dump2(int max = -1) const;
  
  // Useful Utility Functions
  grid<T>& operator<<(grid &m);
  grid<T>& operator+=(const grid &m);
  grid<T>& operator-=(const grid &m);
  grid<T>& operator+=(int i);
  grid<T>& operator-=(int i) { return *this += -i;  }
  grid<T> operator*(const grid &m) const;
  void normalize(T val);
  void normalize(T minval, T maxval);
  const grid<T> transpose() const;
  grid<T> LU() const;
  grid<T> inverse() const;
  int offpixels() { 
    int c = 0; for (int i = 0; i < nx * ny; ++i) c += (sto[i] == 0); return c;
  }
  int onpixels() { return(nx * ny - offpixels()); }
  int compar(const void *a, const void *b) { 
    return (*((const int *)a) < *((const int *)b));
  }
  void sort(int left = -1, int right = -1);
  void sort_columns(int nel, int (*compar) (const void *, const void *));
  
}; // class grid

// __________________________________________________________________________
// Dump a grid; Will only work for classes with ostream<<(const T&).

template <class T>
void grid<T>::dump(int max) const
{
  if (max == -1)
    max = ny;
  for (int j = 0; j < max; ++j) {
    for (int i = 0; i < nx; ++i)
      cout << ((*this)(i, j)) << " ";
    cout << "\n";
  }

} // grid::dump

// __________________________________________________________________________
// Dump a grid; Will only work for classes with ostream<<(const T&).

template <class T>
void grid<T>::dump2(int max) const
{
  if (max == -1)
    max = nx;
  for (int i = 0; i < max; ++i) {
    for (int j = 0; j < ny; ++j)
      cout << ((*this)(i, j)) << " ";
    cout << "\n";
  }

} // grid::dump2

// __________________________________________________________________________
// Calculate the size of a grid which is to be written out

template<class T>
int grid<T>::size() const
{
  assert(invariant());
  return (sizeof(nx) + sizeof(ny) + nx * ny * sizeof(T));

} // grid::size

// __________________________________________________________________________
// Set a grid equal to a grid 'm'; leave 'm' unchanged.

template<class T>
grid<T>& grid<T>::operator=(const grid& m)
{
  assert(invariant() && m.invariant());
  if (this != &m) {
    init(m.nx, m.ny);
    if (nx > 0 && ny > 0) 
      memcpy(sto, m.sto, nx * ny * sizeof(T));
  }
  assert(invariant() && m.invariant());
  
  return *this;
  
} // grid::operator=

// __________________________________________________________________________
// Set 'm' equal to a subgrid of *this.

template<class T>
void grid<T>::subgrid(grid& m, int x, int y, int xsize, int ysize) const
{
  assert(invariant() && m.invariant());
  assert(x >= 0 && y >= 0 && xsize >= 0 && ysize >= 0);
  assert(x + xsize <= dimx() && y + ysize <= dimy());
  
  if (this != &m) {
    m.init(xsize, ysize);
    if (nx > 0 && ny > 0) 
      for (int i = 0; i < m.nx; ++i)
	for (int j = 0; j < m.ny; ++j)
	  m(i, j) = (*this)(i + x, j + y);
  }
  else {
    grid tmp;
    subgrid(tmp, x, y, xsize, ysize);
    m << tmp;
  }

  assert(invariant() && m.invariant());
  
} // grid::subgrid

// __________________________________________________________________________
// Set a grid equal to 'm'; obliterate 'm'.

template<class T>
grid<T>& grid<T>::operator<<(grid &m)
{
  assert(invariant() && m.invariant());
  if (this != &m) {
    freegrid(); 
    nx = m.nx; ny = m.ny; sto = m.sto;
    m.initnew(0, 0);
  }
  assert(invariant() && m.invariant());

  return *this;

} // grid::operator<<
  
// __________________________________________________________________________
// Add another grid to the current grid.

template<class T>
grid<T>& grid<T>::operator+=(const grid& m)
{
  assert(invariant() && m.invariant());
  assert(nx == m.nx && ny == m.ny);

  for (int i = 0; i < nx * ny; ++i)
    sto[i] += m.sto[i];

  assert(invariant() && m.invariant());

  return *this;

} // grid::operator+=

// __________________________________________________________________________
// Subtract a grid from the current grid.

template<class T>
grid<T>& grid<T>::operator-=(const grid& m)
{
  assert(invariant() && m.invariant());
  assert(nx == m.nx && ny == m.ny);

  for (int i = 0; i < nx * ny; ++i)
    sto[i] -= m.sto[i];

  assert(invariant() && m.invariant());

  return *this;

} // grid::operator-=

// __________________________________________________________________________
// Add a fixed value to each element of a grid

template<class T>
grid<T>& grid<T>::operator+=(int ii)
{
  for (int i = 0; i < nx * ny; ++i)
    sto[i] += ii;

  return *this;

} // grid::operator+=

// __________________________________________________________________________

template<class T>
void grid<T>::normalize(T val)
{
  T gridmax = 0;

  assert(invariant());
  int i, j;
  for (i = 0; i < nx; ++i)
    for (j = 0; j < ny; ++j)
      if ((*this)(i, j) > gridmax)
        gridmax = (*this)(i, j);
  if (gridmax > 0) {
    for (i = 0; i < nx; ++i)
      for (j = 0; j < ny; ++j)
        (*this)(i, j) = (*this)(i, j) * val / gridmax;
  }
  assert(invariant());

} // grid::normalize

// __________________________________________________________________________
// Perform a linear transformation on the values of a grid so that the
// values range from 'val1' to 'maxval'

template<class T>
void grid<T>::normalize(T val1, T val2)
{
  if (dimx() == 0 || dimy() == 0)
    return;
  
  T gridmin = this->sto[0];
  T gridmax = this->sto[0];

  assert(invariant());
  int i, j;
  for (i = 0; i < nx; ++i)
    for (j = 0; j < ny; ++j) {
      if ((*this)(i, j) > gridmax)
        gridmax = (*this)(i, j);
      if ((*this)(i, j) < gridmin)
        gridmin = (*this)(i, j);
    }
  
  if (gridmax > gridmin) {
    for (i = 0; i < nx; ++i)
      for (j = 0; j < ny; ++j)
        (*this)(i, j) = ((*this)(i, j) * (val2 - val1) +
			 (gridmax * val1 - gridmin * val2)) /
			   (gridmax - gridmin);
  }
  else if (gridmin < val1)
    clear(val1);
  else if (gridmax > val2)
    clear(val2);

  assert(invariant());

} // grid::normalize

// __________________________________________________________________________
// Multiply two grids together using matrix multiplication.

template<class T>
grid<T> grid<T>::operator*(const grid &m) const
{
  assert(invariant() && m.invariant());
  assert(ny == m.nx);

  grid<T> tmp(nx, m.ny);
  tmp.clear(0);
  for (int i = 0; i < nx; ++i)
    for (int j = 0; j < m.ny; ++j)
      for (int k = 0; k < ny; ++k)
	tmp(i, j) += (*this)(i, k) * m(k, j);

  assert(invariant() && m.invariant() && tmp.invariant());
  return tmp;

} // grid<T>::operator*

// __________________________________________________________________________
// Determine the LU decomposition of a square grid.

template<class T>
grid<T> grid<T>::LU() const
{
  assert(invariant());
  assert(nx == ny);

  grid<T> tmp = *this;
  int i, j;
  for (i = 0; i < nx - 1; ++i)
  {
    for (j = i + 1; j < nx; ++j)
      tmp(j, i) /= tmp(i, i);
    for (j = i + 1; j < nx; ++j)
      for (int k = i + 1; k < nx; ++k)
	tmp(j, k) -= tmp(j, i) * tmp(i, k);
  }

  assert(invariant() && tmp.invariant());
  return tmp;

} // grid::LU

// __________________________________________________________________________
// Calculate the matrix inverse of a grid.

template<class T>
grid<T> grid<T>::inverse() const
{
  assert(invariant());
  assert(nx == ny);
  
  grid<T> tmp = *this;

  grid<int> p(nx);
  int i, j, k;
  for (j = 0; j < nx; ++j)
    p(j) = j;
  grid<double> hv(nx);
  hv.clear(0);

  for (j = 0; j < nx; ++j)
  {
    T max = fabs(tmp(j, j));
    int r = j;
    for (i = j + 1; i < nx; ++i) {
      if (fabs(tmp(i, j)) > max) {
	max = fabs(tmp(i, j));
	r = i;
      }
    }
    if (max == 0.0) {
      cerr << "Unable to invert a matrix" << endl;
      exit(1);
    }
    if (r > j) {
      for (k = 0; k < nx; ++k)
	swap(tmp(j, k), tmp(r, k));
      swap(p(j), p(r));
    }
    T hr = 1 / tmp(j, j);
    for (i = 0; i < nx; ++i)
      tmp(i, j) *= hr;
    tmp(j, j) = hr;
    for (k = 0; k < nx; ++k) {
      if (k == j)
	continue;
      for (i = 0; i < nx; ++i) {
	if (i == j)
	  continue;
	tmp(i, k) -= tmp(i, j) * tmp(j, k);
      }
      tmp(j, k) *= (-hr);
    }
  }
  for (i = 0; i < nx; ++i) {
    for (k = 0; k < nx; ++k)
      hv(p(k)) = tmp(i, k);
    for (k = 0; k < nx; ++k)
      tmp(i, k) = hv(k);
  }
  
  assert(invariant() && tmp.invariant());
  return tmp;

} // grid::inverse

// __________________________________________________________________________

template<class T>
const grid<T> grid<T>::transpose() const
{
  assert(invariant());
  grid<T> tmp(ny, nx);
  for (int i = 0; i < ny; ++i)
    for (int j = 0; j < nx; ++j)
      tmp(i, j) = (*this)(j, i);
  assert(invariant() && tmp.invariant());

  return tmp;

} // grid::transpose

// __________________________________________________________________________

template<class T>
void grid<T>::sort_columns(int nel, int (*compar) (const void *, const void *))
{
  // sort nel columns of the grid according to compar
  qsort(sto, nel, nx * sizeof(T), compar);
}

// __________________________________________________________________________

template<class T>
void grid<T>::sort(int left, int right)
{
  if (left == -1 && right == -1) {
    left = 0;
    right = nx * ny - 1;
  }

  int i = left;
  int j = right;

  T midval = (*this)((left + right) / 2);

  do {
    while ((*this)(i) < midval && i < right)
      ++i;
    while (midval < (*this)(j) && j > left)
      --j;

    if (i <= j) 
    {
      T tmpval = (*this)(i);
      (*this)(i) = (*this)(j);
      (*this)(j) = tmpval;
      ++i;
      --j;
    }

  } 
  while (i <= j);

  if (left < j)
    sort(left, j);
  if (i < right)
    sort(i, right);

} // grid::sort

// __________________________________________________________________________

template <class T>
int grid<T>::loadpgm(const char* pgmname)
{
  // Open the pgm file
  Ifstream(ifs, pgmname);
  if (!ifs)
    return 0;

  // Find the first line that doesn't begin with white space.
  char buf[255];
  char pchar;
  int mode, dimx, dimy, maxval, matches;
  do {
    ifs.Getline(buf, 255);
  } while (!ifs.eof() && (buf[0] == 0 || buf[0] == ' '));
  if (ifs.eof())
    return 0;

  // Make sure the file is a pgm file.  Determine the pgm mode, the
  // dimensions, and the range of pixel values
  matches = sscanf(buf, "%c%1d%d%d%d", &pchar, &mode, &dimx, &dimy, &matches);
  if (matches < 2 || mode < 2 || mode > 6)
    return 0;

  if (matches < 5)
  {
    ifs.Getline(buf, 255);
    while (!ifs.eof() && (buf[0] == '#' || buf[0] == 0))
    ifs.Getline(buf, 255);
    sscanf(buf, "%d %d", &dimx, &dimy);
    ifs.Getline(buf, 255);
    sscanf(buf, "%d", &maxval);
    if (ifs.eof() || dimx <= 0 || dimy <= 0 || maxval <= 0)
      return 0;
  }
  
  init(dimx, dimy);
  int y;
  switch(mode)
  {
   case 2:
    for(y = 0; y < dimy; ++y) {
      int tmpi;
      for(int x = 0; x < dimx; ++x) {
	ifs >> tmpi;
	(*this)(x, y) = tmpi;
      }
    };
    break;
    
   case 3:
    for(y = 0; y < dimy; y++) {
      int tmpi[3];
      for(int x = 0; x < dimx; ++x) {
	ifs >> tmpi[0] >> tmpi[1] >> tmpi[2];
	(*this)(x, y) = uchar(0.212671 * tmpi[0] + 0.715160 * tmpi[1] + 
                              0.072169 * tmpi[2]);
      }
    };
    break;
    
   case 5:
    for(y = 0; y < dimy; ++y) {
      for(int x = 0; x < dimx; ++x)
	ifs.read((char *)&((*this)(x, y)), 1);
    };
    break;
        
   case 6:
    for(y=0; y < dimy; y++) {
      unsigned char tmpc[3];
      for(int x = 0; x < dimx; ++x) {
        ifs.read((char *)tmpc, 3);
	(*this)(x, y) = uchar(0.212671 * tmpc[0] + 0.715160 * tmpc[1] + 
                              0.072169 * tmpc[2]);
      }
    };
    break;
   
   default:
    return 0;
    break;

  };
  
  return 1;

} // grid::loadpgm

template<class T>
const grid<T> operator+(const grid<T>& m, const grid<T>& n) 
{ grid<T> p = m; return (p += n); }

template<class T>
const grid<T> operator-(const grid<T>& m, const grid<T>& n) 
{ grid<T> p = m; return (p -= n); }

//template<class T1, class T2>
//const grid<T1>& operator=(const grid<T1>& m, const grid<T2>& n) 
//{ 
  //n.init(m.nx, m.ny);
  //for (int i = 0; i < m.nx * m.ny; ++i)
    //m[i] = T1(n[i]);
//}

// __________________________________________________________________________
// Local defines

typedef grid<char> cgrid;
typedef grid<unsigned char> ucgrid;
typedef grid<int> igrid;
typedef grid<unsigned int> uigrid;
typedef grid<long> lgrid;
typedef grid<float> fgrid;
typedef grid<double> dgrid;
typedef vector<cgrid> cgrids;
typedef vector<ucgrid> ucgrids;
typedef vector<ucgrids> ucgridss;
typedef vector<igrid> igrids;
typedef vector<uigrid> uigrids;
typedef vector<lgrid> lgrids;
typedef vector<fgrid> fgrids;
typedef vector<dgrid> dgrids;

// __________________________________________________________________________
// grid.h

#endif // GRID_H