#ifndef IMAGE_HH
#define IMAGE_HH
/////////////////////////////////////////////////////////////////////
//! file="amma/Image/BasImage/Image.hh"
//! lib=Image
//! author="Radek Marik"
//! date="26.04.1994"
//! docentry="Image.Image Representations"
//! rcsid="$Id: Image.hh,v 1.52 2000/12/13 15:28:47 ees1cg Exp $"
//! example=exImage.cc

#include <values.h>
#include "amma/StdMath.hh"     //::Cos(), ::Sin()
#include "amma/Error.hh"
#include "amma/BasImage.hh"
#include "amma/Point2d.hh" // in Rotate()
#include "amma/RBfAcc.hh" 
#include "amma/TypeName.hh"
#include "amma/String.hh"
#include "amma/BufferR.hh"

class BinOStreamC;
class BinIStreamC;
class IStreamC;
class OStreamC;

typedef ByteT           ByteGreyValueT;  // byte value

template<class DataC> class ImageIterC;
template<class DataC> class ImagePIterC;
template<class Data1C,class Data2C> class Image2IterC;
template<class Data1C,class Data2C,class Data3C> class Image3IterC;

// ------------------------------------------------------------------------
// *** ImageC *************************************************************
// ------------------------------------------------------------------------

//! userlevel=Basic
//: Image template 
//
// The class ImageC is a container serving as a representation of rectangular
// images defined on uniform grid. Two dimensional indexing can be used
// for access of elements. The elements must have the default constructor,
// the copy constructor, the assigment.


template <class PValueT>
class ImageC: public BaseImageC
{
public:
  typedef PValueT PixelTypeT;
  //: Make pixel type known to routines using ImageC<x>
  
  //!section: Constructors, assigment, copy, and destructor.
  /*----------------------------------------------------------*/
  
  ImageC();
  //: Creates an empty image. 
  
  ImageC(ImageIStreamC & inS);
  //: Creates an image from the input stream.
  // Use of ImageIStreamC is OBSOLETE.
  
  ImageC(IStreamC & inS,BooleanT useHeader = false);
  //: Creates an image from the input stream.
  
  ImageC(ImageFileInfoC * info);
  //: Creates an image according to information in "info".
  
  ImageC(SizeT rows, SizeT columns);
  //: Creates a new image.
  
  ImageC(const ImageRectangleC & rect);
  //: Creates a new image.
  
  ImageC(const ImageC<PValueT> & image, const ImageRectangleC & rect);
  //: Creates a new access to the subimage.
  
  ImageC(const ImageC<PValueT> & image);
  //: Creates a new identifier for the image.
  
  ImageC(SizeT rows, SizeT columns, PValueT * value, BooleanT rm = false,BooleanT rowInv = false);
  //: Creates a new image which is connected to the memory pointed by "value".
  //
  // Use it to create an ImageC object from image data stored in a "C"-style
  // 1-D array.<br>
  //
  // If "rm" is TRUE, the memory pointed to by the pointer "value" is
  // deallocated using delete function during destruction.
  //
  // If row invert is true the image is created 'upside' down. i.e. bottom row first.
  // from the buffer.
  
  ImageC(SizeT rows, SizeT columns,const BufferRC<PValueT> &buff,BooleanT rowInv = false);
  //: Creates a new image which is connected to the memory in 'buff'".
  //
  // Can be used to access data stored in other amma arrays as an image.
  // (MatrixC and VectorC coming soon.)
  //
  // If row invert is true the image is created 'upside' down. i.e. bottom row first.
  // from the buffer.
  
  ImageC<PValueT> Copy() const;
  //: Creates another physical copy of this image.
  
  void Copy(const ImageC<PValueT> & image);
  //: Copies the contents of image "img" into this image.
  
  void Copy(const ImageRectangleC & rect, const PixelC & newOrigin);
  //: Copies the region of the image defined by rectangle "rect" into the 
  //: region whose origin is "newOrigin".
  
  ImageAC * ACopy() const;
  //: A physical copy of this image.
  
  ImageAC * AnotherImage() const;
  //: Creates a new empty image of the same type and size as this image is.
  
  const ImageC<PValueT> & operator=(const ImageC<PValueT> & image);
  //: Assigment of a big object.
  
  virtual ~ImageC();
   //: Destructor.
  
  //!section: Geometrical operations.
  /*---------------------------------*/
  
  void ReflectUp(const ImageRectangleC & rect);
  //: Mirrors the image rectangle along its top edge.
  
  void ReflectDown(const ImageRectangleC & rect);
  //: Mirrors the image rectangle along its bottom edge.
  
  void ReflectLeft(const ImageRectangleC & rect);
  //: Mirrors the image rectangle along its left edge.
  
  void ReflectRight(const ImageRectangleC & rect);
  //: Mirrors the image rectangle along its right edge.
  
  void ReflectFrom(const ImageRectangleC & rect);
  //: Mirrors the image rectangle "rect" along its sides to fill
  //: the whole image.
  
  void Rotate90(const ImageC<PValueT> & originalImage);
  //: The values of image "originalImage" rotated about PI/2 rad 
  //: and shifted to have the original upper-left corner and
  //: are saved into this image.
  
  void Rotate180(const ImageC<PValueT> & originalImage);
  //: The values of image "originalImage" rotated about PI rad 
  //: and shifted to have the original upper-left corner and
  //: are saved into this image.
  
  void Rotate270(const ImageC<PValueT> & originalImage);
  //: The values of image "originalImage" rotated about 3/2*PI rad 
  //: and shifted to have the original upper-left corner and
  //: are saved into this image.
  
  ImageC<PValueT> Rotate(const Point2dC &        center,
                         RealT                   angle, 
                         PValueT                 background) const;
  //: The values of this image are rotated about "angle" (in rad) 
  //: and saved to the result image.  
  //
  // Note the AMMA convention that the 1st or "x" coordinate of Point2dC is
  // associated with the image <b>row no.</b>, and the 2nd or "y" coordinate
  // with the image <b>column</b> no.
  
  void Rotate(const ImageC<PValueT> & originalImage,
              const Point2dC &        center,
              RealT                   angle,   // in radian
              PValueT                 background);
  //: The values of image "originalImage" rotated about "angle" rad 
  //: are saved into this image.
  //
  // Note the AMMA convention that the 1st or "x" coordinate of Point2dC is
  // associated with the image <b>row no.</b>, and the 2nd or "y" coordinate
  // with the image <b>column</b> no.
  
  inline 
    void ShiftRowIndexes(const IndexT offset);
  //: Row indices of the items will be shifted by "offset".
  //
  // There <i>must</i> be only one reference to this image.<br>
  //
  // For +ve "offset", the image is shifted <b>down</b>, or the coordinate
  // system is shifted <b>up</b>.
  
  inline
    void ShiftColIndexes(const IndexT offset);
  //: Column indices of the items will be shifted by "offset".
  //
  // There <i>must</i> be only one reference to this image.<br>
  //
  // For +ve "offset", the image is shifted <b>down</b>, or the coordinate
  // system is shifted <b>up</b>.

  ImageC<PValueT> Subsample (IntT factor=2) const;
  //: returns an image subsampled by a factor of "factor"

  inline
    void SetOrigin(const PixelC & newOrigin);
  //: The image coordinates are changed so that the top left pixel
  //: will have coordinates of "newOrigin".
  //
  // There <i>must</i> be only one reference to this image.
  
  //!section: Grey level operations.
  /*-------------------------------*/
  
  ImageC<PValueT> & Fill(PValueT value);
  //: Fills the image with "value"

  ImageC<PValueT> & Frame(PValueT value);
  //: Fills a 1-pixel-wide border around the image with "value"
  
  ImageC<PValueT> & Frame(PValueT v, const ImageRectangleC & rect);
  //: Fills a 1-pixel-wide rectangle with "v";
  // the rectangle is the intersection of "this" with "rect".

  ImageC<PValueT> Pad(const ImageRectangleC & rect) const;
  //: Creates a new image whose frame is from "rect", with pixels copied from "this".
  //
  // Where the new image lies outside "this", it is padded out with the nearest
  // value from "this".
  
  //!section:  Logical conditions.
  /*-----------------------------*/
  
  inline BooleanT IsValid() const;
  //: Returns FALSE if the image does not have a valid data area allocated
  
  //!section: Access to pixels values.
  /*----------------------------------------*/
  
  inline const RangeBufferAccessC<PValueT>  operator[](IndexT r) const;
  //: Access to the row of the row of pixels.
  
  inline RangeBufferAccessC<PValueT>  operator[](IndexT r);
  //: Access to the row of the row of pixels.
  
  inline PValueT & operator[](const Index2dC & pxl);
  //: Returns the content of the pixel.
  
  inline const PValueT & operator[](const Index2dC & pxl) const;
  //: Returns the content of the pixel.
  
  inline PValueT & operator[](const PixelC & pxl);
  //: Returns the content of the pixel.
  
  inline const PValueT & operator[](const PixelC & pxl) const;
  //: Returns the content of the pixel.
  
  inline const PValueT * CRow(PixelIndexT row) const;
  //: Access to the beginning of the row, eg. t
  //
  // The returned pointer points to the first pixel of the "row".
  
  inline PValueT * Row(PixelIndexT row);
  //: Access to the beginning of the constant "row"
  //
  // I.e. the returned pointer points to the first pixel of "row".
  
  void Save(ostream & outS,BooleanT useHeader = TRUE) const;
  //: Saves templated image to stream.
  //
  // For byte images it will save them in either PGM or PPM binary format.  For
  // other types it will save them in ascii if a suitable stream operator has
  // been written for the templated class.
  //
  // NB. These treat images of ByteGreyValueT specialy, the image values
  // are read as numbers not as charactors.
  
  void Load(istream & outS,BooleanT useHeader = TRUE);
  //: Loads templated image type from ascii stream.
  //
  // Should not be used for binary PGM/PPM files, but should work on most
  // templated types.
  //
  // NB. These treat images of ByteGreyValueT specialy, the image values
  // are written as numbers not as charactors.
  
  PValueT ** StartAddress() const;
  //: Returns the pointer to the table of pointers to image.
  // This being public is really only a hack!
  // EXPERTS ONLY!
  
  BufferRC<PValueT> &Buffer()
    { return buffer; }
  //: Access buffer where data is stored.
  
  const BufferRC<PValueT> &Buffer() const
    { return buffer; }
  //: Access buffer where data is stored.
  
protected:
  
  inline const PValueT * CMemRow(PixelIndexT row) const;
  //: Access to the beginning of the constant allocated row
  //
  // I.e. the returned pointer points to the first pixel of "row".
  
  inline PValueT * MemRow(PixelIndexT row);
  //: Access to the beginning of the allocated row
  //
  // I.e. the returned pointer points to the first pixel of "row".
    
  //!section: Error checking.
  /*------------------------*/
  
  inline void ErrNotValid(char * functionName) const;
  //: Triggers the error event if this image does not contain any data.
  
private:
  
  //!section: Memory management.
  /*---------------------------*/
  
  void AllocData();  
  //: Allocates a memory for image.
  
  void AllocData(PValueT * value,BooleanT rowInv = false);
  //: Allocates a access vector to the memory with data
  
  void DeallocAccess();
  //: Deallocates the access vector to the memory with data.
  
  void DeallocData();
  //: Deallocates the access vector and the memory containing data.
  
private:

  //!section: Object description.
  /*-----------------------------*/
  
  PValueT          **img;
  //: table of pointers to image rows
  
  BufferRC<PValueT>  buffer;
};

template <class PValueT> 
ostream &operator<< (ostream &out,const ImageC<PValueT> &img) {
  img.Save(out,FALSE);
  return out;
}

template <class PValueT> 
istream &operator>> (istream &in,ImageC<PValueT> &img) { 
  img.Load(in,FALSE);
  return in;
}

#include "amma/Error.hh"
#include "amma/ImageRec.hh"
#include "amma/ImagIter.hh"
#include "amma/Image2Iter.hh"
#include "amma/ImgPIter.hh"

// --------------------------------------------------------------------
// **************** ImageC ********************************************
// --------------------------------------------------------------------

template <class PValueT>
inline
void 
ImageC<PValueT>::ErrNotValid(char * functionName) const {
  if (img == 0)
  {
    errAMMA << "The image data are not valid.";
    errAMMA.Function(functionName).Exit();
  }
}

template <class PValueT>
inline
PValueT ** 
ImageC<PValueT>::StartAddress() const
{ return img; }


template <class PValueT>
inline
const PValueT *
ImageC<PValueT>::CRow(PixelIndexT row) const {
#ifdef AMMA_CHECK
  ErrNotInAllocRange(row, "ImageC<PValueT>::CRow(PixelIndexT row) const");
#endif
  return img[row.V()]+LCol().V();
}

template <class PValueT>
inline
PValueT *
ImageC<PValueT>::Row(PixelIndexT row) {
#ifdef AMMA_CHECK
  ErrNotInAllocRange(row, "ImageC<PValueT>::Row(PixelIndexT row)");
#endif
  return img[row.V()]+LCol().V();
}

template <class PValueT>
inline
const PValueT *
ImageC<PValueT>::CMemRow(PixelIndexT row) const {
#ifdef AMMA_CHECK
  ErrNotInAllocRange(row, "ImageC<PValueT>::CRow(PixelIndexT row) const");
#endif
  return img[row.V()]+allRect->LCol().V();
}

template <class PValueT>
inline
PValueT *
ImageC<PValueT>::MemRow(PixelIndexT row) {
#ifdef AMMA_CHECK
  ErrNotInAllocRange(row, "ImageC<PValueT>::Row(PixelIndexT row)");
#endif
  return img[row.V()]+allRect->LCol().V();
}

template <class PValueT>
inline
RangeBufferAccessC<PValueT>
ImageC<PValueT>::operator[](PixelIndexT row) {
#ifdef AMMA_CHECK
  ErrNotInRange(row, "ImageC<PValueT>::operator[](PixelIndexT row)");
#endif
  return RangeBufferAccessC<PValueT>(IndexRangeC(LCol(),RCol()),img[row.V()]);
}

template <class PValueT>
inline
const RangeBufferAccessC<PValueT>
ImageC<PValueT>::operator[](PixelIndexT row) const {
#ifdef AMMA_CHECK
  ErrNotInRange(row, "ImageC<PValueT>::operator[](PixelIndexT row) const");
#endif
  return RangeBufferAccessC<PValueT>(IndexRangeC(LCol(),RCol()),img[row.V()]);
}

template <class PValueT>
inline
PValueT &
ImageC<PValueT>::operator[](const Index2dC & pxl) {
#ifdef AMMA_CHECK
  ErrNotInRange(pxl, "ImageC<PValueT>::operator[](const Index2dC & pxl)");
#endif
  return img[pxl.Row().V()][pxl.Col().V()];
}

template <class PValueT>
inline
const PValueT &
ImageC<PValueT>::operator[](const Index2dC & pxl) const {
#ifdef AMMA_CHECK
  ErrNotInRange(pxl, "ImageC<PValueT>::operator[](const Index2dC &) const");
#endif
  return img[pxl.Row().V()][pxl.Col().V()];
}

template <class PValueT>
inline
PValueT &
ImageC<PValueT>::operator[](const PixelC & pxl) {
#ifdef AMMA_CHECK
  ErrNotInRange(pxl, "ImageC<PValueT>::operator[](const PixelC &) const");
#endif
  return img[pxl.Row().V()][pxl.Col().V()];
}

template <class PValueT>
inline
const PValueT &
ImageC<PValueT>::operator[](const PixelC & pxl) const {
#ifdef AMMA_CHECK
  ErrNotInRange(pxl, "ImageC<PValueT>::operator[](const PixelC &) const");
#endif
  return img[pxl.Row().V()][pxl.Col().V()];
}

template <class PValueT>
inline
void 
ImageC<PValueT>::ShiftRowIndexes(const IndexT offset) {
  if (!IsSingleReference()) { 
    // There is more references to this image.
    errMaMa << "there is more than one references to this image";
    errMaMa.Function("ImageC<PValueT>::ShiftRowIndexes(const IndexT)")
      .Exit();
  }
  img -= offset.V();
  ShiftRows(offset);
}

template <class PValueT>
inline
void 
ImageC<PValueT>::ShiftColIndexes(const IndexT offset) {
  if (!IsSingleReference()) { 
    // There is more references to this image.
    errMaMa << "there is more than one references to this image";
    errMaMa.Function("ImageC<PValueT>::ShiftColIndexes(const IndexT)")
           .Exit();
  }
  ShiftCols(offset);
  for (IndexT r = allRect->TRow(); r <= allRect->BRow(); ++r)
    img[r.V()] -= offset.V();
}

template <class PValueT>
ImageC<PValueT> 
ImageC<PValueT>::Subsample (IntT factor) const {
  ImageC<PValueT> subsampled (ImageRectangleC(
          (this->TRow()-1)/factor+1, this->BRow()/factor,
	  (this->LCol()-1)/factor+1, this->RCol()/factor));
  if (!subsampled.IsValid())  return ImageC<PValueT>();
  for (ImagePIterC<PValueT> it(subsampled);it.IsElm();it.Next())
    it.Data() = (*this)[factor*it.Pixel()];
  return subsampled;
}

template <class PValueT>
inline
void 
ImageC<PValueT>::SetOrigin(const PixelC & newOrigin) {
  PixelC offset(newOrigin - Rectangle().Origin());
  ShiftColIndexes(offset.Col());
  ShiftRowIndexes(offset.Row());
}

template <class PValueT>
inline 
BooleanT 
ImageC<PValueT>::IsValid() const 
{ return img != NULL; }

#include "amma/StdError.hh"

template <class PValueT>
void
ImageC<PValueT>::AllocData(PValueT * value,BooleanT rowInv) {
  if (value == NULL) {
    errAMMA << "no memory for image data";
    errAMMA.Function("ImageC<PValueT>::AllocData(PValueT * value)").Exit();
  }
  SizeT rows = Rnum();
  img = new PValueT *[rows];
  if(img == NULL) {
    errAMMA << "no memory for pointer table";
    errAMMA.Function("ImageC<PValueT>::AllocData(PValueT * value)").Exit();
  }
  
  // set pointers to the rows
  SizeT columns = Cnum();
  if(!rowInv) {
    for (MUIndexT r = 0; r < rows; r++) {
      img[r] = value - Origin().Col().V();
      value += columns;
    }
  } else {
    for (IntT r = (UIntT)(((int) rows)-1); r >= 0; r--) {
      img[r] = value - Origin().Col().V();
      value += columns;
    }
  }
  img -= TRow().V();
}

template <class PValueT>
void
ImageC<PValueT>::AllocData()
{
  if (IsBlock()) {
    LongIntT imageSize = Area();
    buffer = BufferRC<PValueT>(imageSize);
    AllocData(buffer.ReferenceElm());
  } else {
    SizeT rows = Rnum();
    img = new PValueT *[rows];
    if(img == NULL) {
      errAMMA << "no memory for pointer table";
      errAMMA.Function("ImageC<PValueT>::AllocData()").Exit();
    }
    
    // set pointers to the rows
    SizeT columns = Cnum();
    for (MUIndexT r = 0; r < rows; r++)
    {
      PValueT * px = new PValueT[columns];
      if (px == NULL)
      {
        errAMMA << "no memory for image data";
        errAMMA.Function("ImageC<PValueT>::AllocData()").Exit();
      }
      img[r] = px - Origin().Col().V();
    }
    img -= Origin().Row().V();
  }
}

template <class PValueT>
void 
ImageC<PValueT>::DeallocAccess() {
  delete [] (img + allRect->Origin().Row().V());
  #ifdef AMMA_CHECK
    img = 0;
  #endif
}

template <class PValueT>
void 
ImageC<PValueT>::DeallocData() {
  if (isInnerArray) {
    if (IsBlock()) {
      buffer.Invalidate();
      //delete [] MemRow(allRect->TRow());
#ifdef AMMA_CHECK
      //img[allRect->TRow().V()] = 0;
#endif
    } else {
      for (int r = allRect->TRow().V(); r <= allRect->BRow().V(); r++) {
	delete [] MemRow(r);
#ifdef AMMA_CHECK
	img[r] = 0;
#endif
      }
    }
  }
  delete [] (img + allRect->TRow().V());
#ifdef AMMA_CHECK
  img = 0;
#endif
}

template <class PValueT>
ImageC<PValueT>::ImageC()
  : BaseImageC(),
    img(NULL)
{}

template <class PValueT>
ImageC<PValueT>::ImageC(ImageIStreamC & inS)
  : BaseImageC(inS, TRUE)
{
  if(Rectangle().IsValid())
    AllocData();
}

template <class PValueT>
ImageC<PValueT>::ImageC(ImageFileInfoC * info)
  : BaseImageC(info, TRUE)
{
  if(Rectangle().IsValid())
    AllocData();
}

template <class PValueT>
ImageC<PValueT>::ImageC(SizeT rows, SizeT columns)
  : BaseImageC(rows, columns, rows>0 && columns>0),
    img(NULL)
{
  if (rows>0 && columns>0)
    AllocData();
}

template <class PValueT>
ImageC<PValueT>::ImageC(const ImageRectangleC & rect)
  : BaseImageC(rect, rect.IsValid()),
    img(NULL)
{
  if (rect.IsValid())
    AllocData();
}

template <class PValueT>
ImageC<PValueT>::ImageC(const ImageC<PValueT> & image)
  : BaseImageC(image),
    img(image.img),
    buffer(image.buffer)
{}

template <class PValueT>
ImageC<PValueT>::ImageC(const ImageC<PValueT> & image,
			const ImageRectangleC & rect)
  : BaseImageC(image, rect),
    img(image.img),
    buffer(image.buffer)
{}

template <class PValueT>
ImageC<PValueT>::ImageC(SizeT rows,SizeT columns, 
                        PValueT * value,BooleanT rm,BooleanT rowInv)
  : BaseImageC(rows, columns, rm)
{
  buffer = BufferRC<PValueT>(value,rm);
  if(rows > 0 && columns > 0) // Check image has some area...
    AllocData(value,rowInv);
}

template <class PValueT>
ImageC<PValueT>::ImageC(SizeT rows,SizeT columns, 
                        const BufferRC<PValueT> &buff,BooleanT rowInv)
  : BaseImageC(rows, columns,TRUE),
    buffer(buff)
{
  // Check image has some area...
  if(rows > 0 && columns > 0 && buffer.IsValid()) 
    AllocData(buffer.ReferenceElm(),rowInv);
}


template <class PValueT>
const ImageC<PValueT> & 
ImageC<PValueT>::operator=(const ImageC<PValueT> & image)
{
  if (this != &image) {
    if (ToBeDeleted())
      DeallocData();
    BaseImageC::operator=(image);
    img = image.img;
    buffer = image.buffer;
  }
  return *this;
}

template <class PValueT>
ImageC<PValueT>::~ImageC()
{
  if (ToBeDeleted() && allRect->IsValid())
    DeallocData();
}

template <class PValueT>
ImageAC * 
ImageC<PValueT>::ACopy() const
{ return new ImageC<PValueT>(Copy()); }

template <class PValueT>
ImageAC * 
ImageC<PValueT>::AnotherImage() const
{ return new ImageC<PValueT>(Rectangle()); }

template <class PValueT>
void     
ImageC<PValueT>::Copy(const ImageC<PValueT> & image)
{
  ImageRectangleC destination(Rectangle());
  destination.Clip(image.Rectangle());
  // FIXME:- Speed up !
  FOR_RECTANGLE(destination, px, m1,m2,m3)
    (*this)[px] = image[px];
}

template <class PValueT>
ImageC<PValueT> 
ImageC<PValueT>::Copy() const {
  ImageC<PValueT> copy(Rectangle());
  for(Image2IterC<PValueT,PValueT> it(copy,*this);it.IsElm();it.Next())
    it.Data1() = it.Data2();
  return copy;
}

template <class PValueT>
void 
ImageC<PValueT>::Copy(const ImageRectangleC & rect, const PixelC & newOrigin)
{
  PixelC offset(newOrigin - rect.Origin());
  if (offset == PixelC(0,0))
    return;

  // Construct the destination rectangle.
  ImageRectangleC destination(rect);
  destination.Clip(Rectangle());   // the source domain
  destination += offset;
  destination.Clip(Rectangle());   // destination.Area() <= source.Area()
  // Both 'destination' and 'source' rectangles are inside the image.

  const IndexT minRow = destination.Origin().Row();
  const IndexT maxRow = destination.End().Row();
  const IndexT minCol = destination.Origin().Col();
  const IndexT maxCol = destination.End().Col();
  PixelC pxl;
  IndexT & r = pxl.Row();
  IndexT & c = pxl.Col();
  if (offset.Row() <= 0)
    if (offset.Col() <= 0)
      for (r = minRow; r <= maxRow; r++)
        for (c = minCol; c <= maxCol; c++)
          (*this)[pxl] = (*this)[pxl-offset];
    else 
      // offset.Col() > 0
      for (r = minRow; r <= maxRow; r++)
        for (c = maxCol; c >= minCol; c--)
          (*this)[pxl] = (*this)[pxl-offset];
  else
    // offset.Row() > 0
    if (offset.Col() <= 0)
      for (r = maxRow; r >= minRow; r--)
        for (c = minCol; c <= maxCol; c++)
          (*this)[pxl] = (*this)[pxl-offset];
    else 
      // offset.Col() > 0
      for (r = maxRow; r >= minRow; r--)
        for (c = maxCol; c >= minCol; c--)
          (*this)[pxl] = (*this)[pxl-offset];
}

template <class PValueT>
void 
ImageC<PValueT>::ReflectUp(const ImageRectangleC & rect)
{
  // .... top - 1 - (row-top)
  const IndexT top21 = 2*rect.TRow() - 1;
  FOR_RECTANGLE(rect, pxl, m1,m2,m3)
    (*this)[PixelC(top21 - pxl.Row(), pxl.Col())] = (*this)[pxl];
}

template <class PValueT>
void 
ImageC<PValueT>::ReflectDown(const ImageRectangleC & rect)
{
  // .... bottom + 1 + (bottom-row)
  const IndexT bottom21 = 2*rect.BRow() + 1;
  FOR_RECTANGLE(rect, pxl, m1,m2,m3)
    (*this)[PixelC(bottom21 - pxl.Row(), pxl.Col())] = (*this)[pxl];
}

template <class PValueT>
void 
ImageC<PValueT>::ReflectLeft(const ImageRectangleC & rect)
{
  // .... left - 1 - (col-left)
  const IndexT left21 = 2*rect.LCol() - 1;
  FOR_RECTANGLE(rect, pxl, m1,m2,m3)
    (*this)[PixelC(pxl.Row(), left21 - pxl.Col())] = (*this)[pxl];
}

template <class PValueT>
void 
ImageC<PValueT>::ReflectRight(const ImageRectangleC & rect)
{
  // .... right + 1 + (right-col)
  const IndexT right21 = 2*rect.RCol() + 1;
  FOR_RECTANGLE(rect, pxl, m1,m2,m3)
    (*this)[PixelC(pxl.Row(), right21 - pxl.Col())] = (*this)[pxl];
}

template <class PValueT>
void 
ImageC<PValueT>::ReflectFrom(const ImageRectangleC & rect)
{
  ImageRectangleC center(Rectangle());
  center.ClipCols(rect); 
  ImageRectangleC top(center.TopDiff(rect).ReflectDown());
  ImageRectangleC bot(center.BottomDiff(rect).ReflectUp()); 
  ImageRectangleC left(Rectangle().LeftDiff(center).ReflectRight()); 
  ImageRectangleC right(Rectangle().RightDiff(center).ReflectLeft());

  top.ErrIfNotInside(rect,"ImageC<PValueT>::ReflectFrom(...)"); 
  bot.ErrIfNotInside(rect,"ImageC<PValueT>::ReflectFrom(...)"); 
  left.ErrIfNotInside(center,"ImageC<PValueT>::ReflectFrom(...)"); 
  right.ErrIfNotInside(center,"ImageC<PValueT>::ReflectFrom(...)"); 

  ReflectUp(top);
  ReflectDown(bot);
  ReflectLeft(left);
  ReflectRight(right);
}

template <class PValueT>
ImageC<PValueT> &
ImageC<PValueT>::Fill(PValueT value)
{
  for(ImageIterC<PValueT> it(*this);it.IsElm();it.Next())
    it.Data() = value;
  return *this;
}

template <class PValueT>
ImageC<PValueT> &
ImageC<PValueT>::Frame(PValueT frame)
{
  // Do the top and bottom rows.
  PValueT *it1,*it2,*eor;
  
  it1 = &((*this)[Rectangle().TRow()][Rectangle().LCol().V()]);
  it2 = &((*this)[Rectangle().BRow()][Rectangle().LCol().V()]);
  IntT ColN = Cnum();
  eor = &(it1[ColN]);
  for(;it1 != eor;) {
    *(it1++) = frame;
    *(it2++) = frame;
  }
  ColN--;
  for(PixelIndexT r = Origin().Row()+1; r < End().Row(); r++) {
    it1 = &((*this)[r][Rectangle().LCol().V()]);
    it1[0] = frame;
    it1[ColN] = frame;
  }
  return *this;
}

template <class PValueT>
ImageC<PValueT> & 
ImageC<PValueT>::Frame(PValueT v, const ImageRectangleC & rect)
{
  ImageRectangleC cRect(rect);
  cRect.Clip(Rectangle());
  PixelC px1(cRect.Origin());
  PixelC px2(cRect.End());

  for(PixelIndexT c = cRect.Origin().Col(); c <= cRect.End().Col();c++) {
    operator[](px1) = v;
    operator[](px2) = v;
    px1.Right();
    px2.Left();
  }

  px1.Left();
  px2.Right();

  for(PixelIndexT r = cRect.Origin().Row(); r <= cRect.End().Row(); r++) {
    operator[](px1) = v;
    operator[](px2) = v;
    px1.Down();
    px2.Up();
  }
  return *this;
}

template <class PValueT>
ImageC<PValueT>  
ImageC<PValueT>::Pad(const ImageRectangleC & rect) const
{
  ImageC<PValueT> padded(rect);
#ifdef AMMA_CHECK
  if (IsEmpty()) {
    errAMMA << "empty original image";
    errAMMA.Function("ImageC<PValueT>ImageC<PValueT>::Pad(const ImageRectangleC & rect)").Exit(-1);
  }
#endif
  FOREACH_PIXEL(padded,px)
    padded[px] = operator[](PixelC(px).Cut(Rectangle()));
  return padded;
}


template <class PValueT>
void
ImageC<PValueT>::Rotate90(const ImageC<PValueT> & originalImage)
// Ref.: Radek Marik: Math, p.27
{
  ImageRectangleC rect(originalImage.Rectangle());
  rect.Rotate90().Clip(Rectangle()); // rect is area where the result is saved.
  const PixelIndexT offRow = originalImage.TRow()-originalImage.LCol();//xL-yL
  const PixelIndexT offCol = originalImage.TRow()+originalImage.RCol();//xL+yU
  FOR_RECTANGLE(rect, pxl, m1, m2, m3)
    (*this)[pxl] = originalImage[PixelC(offRow+pxl.Col(), offCol-pxl.Row())];
}

template <class PValueT>
void
ImageC<PValueT>::Rotate180(const ImageC<PValueT> & originalImage)
// Ref.: Radek Marik: Math, p.27
{
  ImageRectangleC rect(originalImage.Rectangle());
  rect.Rotate180().Clip(Rectangle());
  //xL+xU, yL+yU
  const PixelC offset(originalImage.TRow() + originalImage.BRow(), 
                      originalImage.LCol() + originalImage.RCol());
  FOR_RECTANGLE(rect, pxl, m1, m2, m3)
    (*this)[pxl] = originalImage[offset - pxl];
}

template <class PValueT>
void
ImageC<PValueT>::Rotate270(const ImageC<PValueT> & originalImage)
// Ref.: Radek Marik: Math, p.27
{
  ImageRectangleC rect(originalImage.Rectangle());
  rect.Rotate270().Clip(Rectangle());
  const PixelIndexT offRow = originalImage.LCol()+originalImage.BRow();//yL+xU
  const PixelIndexT offCol = originalImage.LCol()-originalImage.TRow();//yL-xL
  FOR_RECTANGLE(rect, pxl, m1, m2, m3)
    (*this)[pxl] = originalImage[PixelC(offRow-pxl.Col(), offCol+pxl.Row())];
}

template <class PValueT>
ImageC<PValueT>
ImageC<PValueT>::Rotate(const Point2dC &        center,
                        RealT                   angle, 
                        PValueT                 background) const
{
  ImageC<PValueT> rotImage(Rectangle());
  rotImage.Rotate(*this,center, angle, background);
  return rotImage;
}
  
template <class PValueT>
void
ImageC<PValueT>::Rotate(const ImageC<PValueT> & originalImage,
                        const Point2dC &        center,
                        RealT                   angle,
                        PValueT                 background)
// modification of a function written by Nicos Georgis
{
  RealT co = ::Cos(angle);
  RealT si = ::Sin(angle);

  FOR_IMAGE_PIXEL(*this, pxl, m1, m2, m3)
  {
    Point2dC newPos(pxl);
    newPos -=center;
    Point2dC oldPos( newPos.X() * co + newPos.Y() * si,
                    -newPos.X() * si + newPos.Y() * co );
    oldPos +=center;
    PixelC oldPixel(oldPos);
    (*this)[pxl] = (oldPixel.IsInside(originalImage.Rectangle()))
                    ? originalImage[oldPixel]
                    : background;
  }
}

template <class PValueT>
void
ImageC<PValueT>::Save(ostream & outS,BooleanT useHeader) const
{

  // Check output stream is good
  // ---------------------------
  if (!outS)
    throw ExceptionC("ImageC<PValueT>::Save(), Passed a bad stream. ");
  if(useHeader)
    outS << TypeName(typeid(ImageC<PValueT>)) << '\n';
  outS << Rectangle() << '\n';
  ImageIterC<PValueT> it(*this);
  if(!it.IsElm()) // Empty image ?
    return ;
  while(1) {
    outS << it.Data();
    if(!it.RNext()) { // End of row/image ?
      if(!it.IsElm()) // Check for end of image.
	break;
      outS << '\n'; // Put carage returns at end of rows.
    } else
      outS << ' '; // Put in pixel seperator.
  }
  // End routine
}

template <class PValueT>
void
ImageC<PValueT>::Load(istream & inS,BooleanT useHeader) 
{
  if (!inS)
    throw ExceptionC("ImageC<PValueT>::Load(), Passed a bad stream. ");
  if(useHeader) {
    StringC header(inS);
    if(header != TypeName(typeid(ImageC<PValueT>)))
      throw ExceptionC("ImageC<PValueT>::Load(), file format not recognised. ");
    //cerr << "Image in file not of format: " << TypeName(typeid(ImageC<PValueT>)) << "\n";
  }
  (*this) = ImageC<PValueT>(LoadImageRect(inS));;
  for(ImageIterC<PValueT> it(*this);it.IsElm();it.Next())
    inS >> it.Data();
}

//: Creates an image from the input stream.

template <class PValueT>
ImageC<PValueT>::ImageC(IStreamC & inS,BooleanT useHeader) 
{ Load(Ins,useHeader); }

#endif

// IAPS - Image analysis program system.
// End of include file Image.hh