/* -*- C++ -*- */
#ifndef __fft_Array2D_Header_File
#define __fft_Array2D_Header_File

#include "types.h"
#include "Image.h"
#include <cstring>
#include <cassert>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>

namespace fft
{

   template <typename T>
   class Array2D 
   {
      private:
	 Array2D<T>& self(){ return *this; }
	 const Array2D<T>& self() const { return *this; }
      public:
	 Array2D(u32 w, u32 h, const T* srcData = 0)
	    : data(0), w(w), h(h)
	 {
	    u32 size = w * h;
	    data = new T[size];
	    
	    if (srcData)
	    {
	       memcpy(data, srcData, size * sizeof(T));
	    }
	 }
	 
	 template <class ST>
	 explicit Array2D(const Array2D<ST>& src)
	    : data(0), w(src.W()), h(src.H())
	 {
	    u32 size = w * h;
	    data = new T[size];
	    for(u32 i = 0; i < size; i++)
	    {
	       data[i] = static_cast<T>(src[i]);
	    }
	 }

	 Array2D(const Array2D<T>& src)
	    : data(0), w(src.w), h(src.h)
	 {
	    u32 size = w * h;
	    data = new T[size];
	    memcpy(data, src.data, size * sizeof(T));
	    //std::cerr << "Copy Constr" << std::endl;
	 }

	 template <typename RhsT>
	 void operator=(const Array2D<RhsT>& rhs)
	 {
	    if (data)
	    {
	       delete [] data;
	    }

	    u32 size = w * h;
	    data = new T[size];
	    for(u32 i = 0; i < size; i++)
	    {
	       self()[i] = static_cast<T>(rhs.data[i]);
	    }
	 }


	 void operator=(const Array2D<T>& rhs)
	 {
	    if (data)
	    {
	       delete [] data;
	    }

	    u32 size = w * h;
	    data = new T[size];
	    for(u32 i = 0; i < size; i++)
	    {
	       data[i] = rhs.data[i];
	    }

	    //std::cerr << "Copy opr (" << min() << "," << max() <<")(" << rhs.min() << "," << rhs.max() << ")"  << std::endl;
	 }


	 ~Array2D()
	 {
	    delete [] data;
	    w = h = 0;
	    data = 0;
	 }

	 ////////////////////////////////////////////////////////////
	 //
	 // Accessors
	 //
	 ///////////////////////////////////////////////////////////
	 T& operator[](u32 i)
	 {
	    assert(i >= 0 && i < w * h);
	    return data[i];
	 }

	 const T& operator[](u32 i) const
	 {
	    assert(i >= 0 && i < w * h);
	    return data[i];
	 }


	 T& at(u32 i)
	 {
	    assert(i >= 0 && i < w * h);
	    return data[i];
	 }

	 const T& at(u32 i) const
	 {
	    assert(i >= 0 && i < w * h);
	    return data[i];
	 }
	

	 const T& at(u32 x, u32 y) const
	 {
	    assert(x >= 0 && y >= 0);
	    assert(x < w && y < h);
	    return data[y * w + x];
	 }

	 T& at(u32 x, u32 y)
	 {
	    assert(x >= 0 && y >= 0);
	    assert(x < w && y < h);
	    return data[y * w + x];
	 }

	 ////////////////////////////////////////////////////
	 // 
	 // Comparsion operators
	 //
	 ////////////////////////////////////////////////////
	 
	 template<typename ST>
	 Array2D<T> operator<(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = data[i] < rhs;
	    }
	    //std::cerr << "min: " << ret.min() << " max: " << ret.max() << std::endl;
	    return ret;
	 }

	 template<typename ST>
	 Array2D<T> operator>(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] > rhs;
	    }
	    return ret;
	 }
	 
	 template<typename ST>
	 Array2D<T> operator<=(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] <= rhs;
	    }
	    return ret;
	 }

	 template<typename ST>
	 Array2D<T> operator>=(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] >= rhs;
	    }
	    return ret;
	 }

	 

	 //////////////////////////////////////////////////////////
	 //
	 // Scalar right hand side argument arithmetic operators
	 //
	 //////////////////////////////////////////////////////////

	 template<typename ST>
	 Array2D<T> operator+(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] + rhs;
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator-(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] - rhs;
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator*(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] * rhs;
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator/(ST rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] / rhs;
	    }
	    return ret;	    
	 }

	 //////////////////////////////////////////////////////////
	 //
	 // Array2D right hand side argument arithmetic operators
	 //
	 //////////////////////////////////////////////////////////

	 template<typename ST>
	 Array2D<T> operator+(Array2D<ST> rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] + rhs[i];
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator-(Array2D<ST> rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] - rhs[i];
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator*(Array2D<ST> rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] * rhs[i];
	    }
	    return ret;	    
	 }

	 template<typename ST>
	 Array2D<T> operator/(Array2D<ST> rhs) const
	 {
	    Array2D ret(w, h);
	    for (u32 i = 0; i < w * h; i ++)
	    {
	       ret[i] = self()[i] / rhs[i];
	    }
	    return ret;	    
	 }

	 ////////////////////////////////////////////////////
	 //
	 // Miscalious functions
	 //
	 ////////////////////////////////////////////////////
	 template <typename ST>
	 Array2D<T> convolution(const Array2D<T>& kernel) const
	 {
	    u32 kw2 = kernel.W()/2;
	    u32 kh2 = kernel.H()/2;
	    Array2D<T> ret(w - kw2 * 2, h - kh2 * 2);


	    return ret;
	 }

	 Array2D<T> medianFilter(u32 kernelSide)
	 {
	    u32 ks2 = kernelSide/2;
	    u32 rw = w - ks2 * 2;
	    u32 rh = h - ks2 * 2;
	    u32 rx = ks2;
	    u32 ry = ks2;
	    Array2D<T> ret(rw, rh);
	    std::vector<T> vals;

	    for (u32 i = 0; i < rw; i++)
	    {
	       for (u32 j = 0; j < rh; j++)
	       {
		  vals.clear();

		  for (u32 k = 0; k < kernelSide; k++)
		  {
		     for (u32 l = 0; l < kernelSide; l++)
		     {
			vals.push_back(at(i + k, j + l));
		     }
		  }

		  sort(vals.begin(), vals.end());
		  ret.at(i,j) = vals[vals.size()/2];
	       }
	    }

	    return ret;
	 }

	 T min() const
	 {
	    T mn = self()[0];
	    for (u32 i = 1; i < w * h; i++)
	    {
	       if (mn > self()[i])
	       {
		  mn = self()[i];
	       }
	    }
	    return mn;
	 }

	 T max() const
	 {
	    T mx = self()[0];
	    for (u32 i = 1; i < w * h; i++)
	    {
	       if (mx < self()[i])
	       {
		  mx = self()[i];
	       }
	    }
	    return mx;
	 }

	 void normalize()
	 {
	    T mn = min(); T mx = max();
	    self() = (self() - mn)/(mx - mn);
	 }

	 void logNormalize()
	 {
	    T* locData = new T[size()];
	    for (u32 i = 0; i < size(); i++)
	    {
	       locData[i] = log(1 + data[i]);
	    }

	    delete [] data;
	    data = locData;

	    normalize();
	 }



	void circle(u32 r, T inVal, T outVal)
        {
	   u32 size = w * h;
	   u32 rp2 = r * r;
	   u32 h2 = h / 2;
	   u32 w2 = w / 2;

	   for (u32 i = 0; i < w; i++)
	   {
	      for (u32 j = 0; j < h; j++)
	      {
		 s32 xd = i - w2;
		 s32 yd = j - h2;
		 u32 dist = xd * xd + yd * yd;
		 self()[j * w + i] = (dist <= rp2) ? inVal : outVal;
	      }
	   }

        }

	 void flipV()
	 {
	    //T* line = new T[w];
	    u32 h2 = h/2;
	    //u32 size = sizeof(T) * w;
	    for (u32 i = 0; i < h2; i++)
	    {
	       /*
	       u32 i1 = i * w;
	       u32 i2 = (h2 - i) * w;
	       memcpy(line, &at(i1), size);
	       memcpy(&at(i1), &at(i2), size);
	       memcpy(&at(i2), line, size);
	       */

	       for (u32 j = 0; j < w; j++)
	       {
		  T tmp = at(j,i);
		  at(j,i) = at(j,h - i - 1);
		  at(j,h - i - 1) = tmp;
	       }
	       
	    }
	    //delete line;
	 }

	 void storeGnuPlot(const String& fname, u32 xStep = 1, u32 yStep = 1)
	 {
	    std::ofstream ofs (fname.c_str(), std::ios::out);
	    for (u32 i = 0; i < h; i += yStep)
	    {
	       for (u32 j = 0; j < w; j += xStep)
	       {
		  //ofs << j << " " << i << " " << at(j,i) << std::endl;
		  ofs << at(j,i) << " ";
	       }
	       ofs << std::endl;
	    }
	    ofs.close();
	 }


	 u32 W() const { return w; }
	 u32 H() const { return h; }
	 u32 size() const { return w * h; }


      private:
	 T* data;
	 u32 w,h;
   };


}


#endif //__fft_Array_Header_File
