/*
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/



/* *
 * \page code_image.h \link image.h \endlink
 * \brief Code of  \ref code_image.h "image.h"
 * \includelineno image.h
 * \brief \b Code: \n
 * - \ref code_image.h "image.h"
*/

/**
 * \file image.h
 * \brief Class to store multichannel image
 * \author Luis Alvarez \n \n
*/
#ifndef image_H
#define image_H
#include <vector>
#include <iostream>
#include <istream>
#include <ostream>
#include <typeinfo>
#include <string>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef __AMI_OPENCV__
  #include <cv.h>
  #include <highgui.h>
  #include <cxcore.h>
#else
  #ifdef __AMI_IMAGE__
      //#include "ami_tif/ami_tif.h"
  #else
      #include <itk_tiff.h>
      #include <itk_png.h>
      extern "C"
      {
        // The regular jpeg lossy lib is the 8bits one:
        #include <itkjpeg/8/jpeglib.h>
        #include <setjmp.h>
        #include <itkjpeg/mangle_jpeg8bits.h>
      }
  #endif
  #include "ami_bmp/ami_bmp.h"
#endif

//FOR BUILDER C++
#ifdef _TOPENCV
  #include <opencv/cv.h>
  #include <opencv/highgui.h>
#endif

#ifdef AMI_OMP_H
  #include <omp.h>
#endif

#ifdef __AMIDEBUG__
  #include "wxAmiDebugLog.h"
#endif

#define IMAGE_DEBUG


namespace ami {

// ITERATOR ALIAS

typedef std::vector<unsigned char>::iterator usIt;
typedef std::vector<int>::iterator iIt;
typedef std::vector<short>::iterator sIt;
typedef std::vector<float>::iterator fIt;
typedef std::vector<double>::iterator dIt;

typedef std::vector<unsigned char>::const_iterator uscIt;
typedef std::vector<int>::const_iterator icIt;
typedef std::vector<short>::const_iterator scIt;
typedef std::vector<float>::const_iterator fcIt;
typedef std::vector<double>::const_iterator dcIt;

}


namespace ami {

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \class image
 * \brief Class to store multiChannel images and basic methods
 * \author Luis Alvarez
 */
template <class T>
class image{
  std::vector <T> image_ /** std vector to allocate image "channel 1 + channel 2 + ..." */;
  int width_ /** image width */;
  int height_ /** image height */;
  int nChannels_ /** number of image channels */;
  std::vector <int> roi_ /** ROI : Region Of Interest =(x0,x1,y0,y1,c0,c1) subwindow corners used (if defined) to operate only in a subwindow area */;

  std::vector <float> origin_ /** world coordinate origin */;
  std::vector <float> pixel_size_ /** pixel size */;

  public:

  // CONSTRUCTORS - DESTRUCTOR METHODS
  image();
  image(const int width,const int height);
  image(const T *vector,const int width,const int height,const int channels);
  image(const int width,const int height, const int nChannels, const T &a);
  image(const int width,const int height, const T &a, const T &b, const T &c);
  image(const int width,const int height, const int nChannels);
  image(std::string name /** INPUT FILE NAME */ );


  image(const image<bool> &img_in, bool isbool);
  image(const image<T> &img_in);      //-----------------------------------------------------------------------------------
  ~image();

  // CLASS ELEMENT ACCESS METHODS
  inline std::vector <T> *get_image(){return &image_;}
  inline int width() const {return width_;}
  inline int height() const {return height_;}
  inline int nChannels() const {return nChannels_;}
  inline int size() const {return image_.size();}
  inline T& operator[](const int &i);
  inline T& operator()(const int &x,const int &y,const int &channel);
  inline const T& operator[](const int &i) const;
//  inline T& operator()(const int &x,const int &y,const int &c);
//  inline const T& operator()(const int &x,const int &y,const int &c) const;
//  inline T& operator()(const int &x,const int &y);
//  inline const T& operator()(const int &x,const int &y) const;


  // BASIC OPERATION METHODS
  int write(std::string name);
  int write_bool(std::string name);
  int read(std::string name);
  image & operator=(const image &image2);
  void clear();
  void init(const int width,const int height, const int nChannels);
  void set_nchannels(int nchannel1);
  void set_size(int width,int height);
  void imageMirrored();
  #ifdef __AMI_OPENCV__
    image<T> resizeImage(int width, int height);
  #endif
  image<T> const resize(const int width, const int height) const;
  image<T> const resize_no_omp(const int width, const int height) const ;
  image<T> const linear_transform(const double H[3][3]) const ;
  image<T> const sampling() const ;

  // BASIC METHODS TO MANIPULATE SUBWINDOWS (ROI : REGION OF INTEREST )
  std::vector <int> get_roi() const{return roi_;}
  void roi_clear(){roi_.clear();}
  void set_roi(const int x0,const int x1,const int y0,const int y1);
  void set_roi(const int x0,const int x1,const int y0,const int y1,
               const int c0,const int c1);
  void set_roi(const std::vector<int> roi);
  template <class U> void  get_roi_image(image<U> &image2) const;
  image<T>  const get_roi_image(const std::vector <int> &roi2) const;
  template <class U> void  set_roi_image(const image<U> &image2);

  // BASIC METHODS TO MANIPULATE IMAGE WITH ITERATORS
  typedef typename std::vector <T>::iterator iterator;
  inline iterator begin() {return(image_.begin());}
  inline iterator end() {return(image_.end());}

  //CHANGE CHANNEL
  void rgb_to_hsv(image<unsigned char> &H,image<unsigned char> &S, image<unsigned char> &V);
  template <class U>
  void hsv_to_rgb(image<U> &image);

  // ACCESS VECTOR BOOL
  T get_value(const int x);
  void set_value(const int x,T value);

  // HISTOGRAM
  template<class U> void get_histogram(std::vector< U > &histograma,int channel);

  private:

    #ifdef __AMI_OPENCV__
        int read_with_opencv(std::string name);
        int write_with_opencv(std::string name);
    #else
        #ifdef __AMI_IMAGE__
            //int read_with_ami_tif (std::string name);
            //int write_with_ami_tif (std::string name);
        #else
            int read_with_jpeglib(std::string name);
            int write_with_jpeglib(std::string name);

            int read_with_pnglib(std::string name);
            int write_with_pnglib(std::string name);
            int write_with_pnglib_08(std::string name);
            int write_with_pnglib_16(std::string name);
            int write_with_pnglib_32(std::string name);

            int read_with_tifflib(std::string name);
            int write_with_tifflib(std::string name);
        #endif

        int read_with_ami_bmp (std::string name);
        int write_with_ami_bmp (std::string name);
    #endif

    int read_pgm(std::string name);
    int write_pgm(std::string name);


};

/**
* \fn template<class T> image<T> const image<T>::sampling() const
* \brief method to downsample the image by a factor of 2. We filter the image before downsampling
*/
template<class T>
image<T> const image<T>::sampling() const
{
  // WE DEFINE SAMPLED IMAGE
  image<T> sampled_img;
  int width_s=width_/2;
  int height_s=height_/2;
  sampled_img=image<T>(width_s,height_s,nChannels_);
  // WE FILTER THE IMAGE AND DOWNSAMPLING IT
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    int n_s=n*width_s*height_s;
    int i_end=height_s-1;
    #pragma omp parallel for shared(i_end,n0,n_s,sampled_img,width_s)
    for(int i=1;i<i_end;i++){
      int m0=n0+2*i*width_;
      int m_s=n_s+i*width_s;
      int j_end=width_s-1;
      for(int j=1;j<j_end;j++){
        int m=m0+2*j;
        sampled_img[m_s+j]=0.25*image_[m]+0.125*(image_[m+1]+image_[m-1]+image_[m+width_]+image_[m-width_])+
                                          0.0625*(image_[m+width_+1]+image_[m+width_-1]+image_[m-width_+1]+image_[m-width_-1]);
      }
    }
  }

  // WE SAMPLE IMAGE BORDER VERTICAL LINES
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    int n_s=n*width_s*height_s;
    for(int i=0;i<height_s;i++){
      int i0=2*i;
      int m0=n0+i0*width_;
      int m_s=n_s+i*width_s;
      for(int j=0;j<width_s;j+=width_s-1){
        int j0=2*j;
        int m=m0+j0;
        int mN=m,mS=m,mW=m,mE=m,mNE=m,mNW=m,mSE=m,mSW=m;
        if(i0>0){ mS-=width_; mSE-=width_; mSW=m-width_; }
        if(i0<(height_-1)){ mN+=width_; mNE+=width_; mNW+=width_; }
        if(j0>0) {mW-=1; mNW-=1; mSW-=1;}
        if(j0<(width_-1)){ mE+=1; mNE+=1; mSE+=1;}

        sampled_img[m_s+j]=0.25*image_[m]+0.125*(image_[mN]+image_[mS]+image_[mW]+image_[mE])+
                                          0.0625*(image_[mNW]+image_[mNE]+image_[mSW]+image_[mSE]);
      }
    }
  }

  // WE SAMPLE IMAGE BORDER HORIZONTAL LINES
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    int n_s=n*width_s*height_s;
    for(int i=0;i<height_s;i+=height_s-1){
      int i0=2*i;
      int m0=n0+i0*width_;
      int m_s=n_s+i*width_s;
      for(int j=1,j_end=width_s-1;j<j_end;j++){
        int j0=2*j;
        int m=m0+j0;
        int mN=m,mS=m,mW=m,mE=m,mNE=m,mNW=m,mSE=m,mSW=m;
        if(i0>0){ mS-=width_; mSE-=width_; mSW=m-width_; }
        if(i0<(height_-1)){ mN+=width_; mNE+=width_; mNW+=width_; }
        if(j0>0) {mW-=1; mNW-=1; mSW-=1;}
        if(j0<(width_-1)){ mE+=1; mNE+=1; mSE+=1;}

        sampled_img[m_s+j]=0.25*image_[m]+0.125*(image_[mN]+image_[mS]+image_[mW]+image_[mE])+
                                          0.0625*(image_[mNW]+image_[mNE]+image_[mSW]+image_[mSE]);
      }
    }
  }

//  static int cont =0;
//  cont++;
//  char name[300];
//  sprintf(name,"sampled_%d.jpg",100+cont);
//  sampled_img.write(name);
  return(sampled_img);
}

/**
 * \fn template <class T> template <class U> void image<T>::get_histogram(std::vector<U> &histograma, int channel)
 * \brief Get histogram of a channel of image
 * \author Carlos Falcon
*/
template <class T> template <class U>
void image<T>::get_histogram(std::vector<U> &histograma/** histograma de dimension 256*/,
                              int channel/**channel to obtain the histogram*/)
{
  if (nChannels_<channel)
    return;
  if (channel <= 0)
    return;
  int size=width_*height_;
  histograma.resize(256);
  #pragma omp parallel for
  for(int i=0;i<256;i++) histograma[i]=0;

  int i_end=size*channel;
  //#pragma omp parallel for default(none) shared(size, channel,i_end)
  for(int i=size*(channel-1);i<i_end;i++)
          histograma[image_[i]]++;
}




/**
 * \fn template <class T> void image<T>::init(const int width,const int height, const int nChannels)
 * \brief initializes the image taking memory
 * \author Carlos Falcon
*/
template <class T>
void image<T>::init(const int width,const int height, const int nChannels)
{
  width_=width;
  height_=height;
  nChannels_=nChannels;
  int size=width_*height_*nChannels_;
  image_.resize(size);
}


//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> void image<T>::set_value(const int x,T value)
 * \brief set value to x position of vector
 * \author Carlos Falcon
*/
template <class T>
void image<T>::set_value(const int x,T value)
{
  if (x<=width_*height_*nChannels_)
    image_[x]=value;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> T image<T>::get_value(const int x)
 * \brief get value to x position of vector
 * \author Carlos Falcon
*/
template <class T>
T image<T>::get_value(const int x)
{
  T value=0;
  if (x<=width_*height_*nChannels_)
    return image_[x];
  return value;
}


//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T>::image(const T *vector,const int width,const int height,const int channels)
 * \brief image from T vector
 * \author Carlos Falcon
*/

template <class T>
image<T>::image(const T *vector,const int width,const int height,const int channels)
{
  width_=width ;
  height_=height;
  nChannels_=channels;
  int size=width_*height_*nChannels_;
  image_.resize(size);
  size=width_*height_;
  #pragma omp parallel \
    shared(size, vector)
    for(int n=0;n<nChannels();n++){
      int n2_=n*size;
      #pragma omp for nowait
      for(int i=0;i<size;i++)
      {
        image_[i+n2_]=vector[i+n2_];
      }
    }

//  for (int i=0;i< size;i++)
//    image_[i]=vector[i];
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T>::image(const image<T> &img_in)
 * \brief image copy constructor
 * \author Carlos Falcon
*/

template <class T>
image<T>::image(const image<T> &img_in)
{
  width_=img_in.width() ;
  height_=img_in.height();
  nChannels_=img_in.nChannels();
  int size=width_*height_*nChannels_;
  image_.resize(size);
  size=width_*height_;
  #pragma omp parallel \
    shared(size)
    for(int n=0;n<img_in.nChannels();n++){
      int n2_=n*size;
      #pragma omp for nowait
      for(int i=0;i<size;i++)
      {
        image_[i+n2_]=img_in[i+n2_];
      }
    }

    if(img_in.get_roi().size() == 6){
        roi_.resize(6);
        roi_.at(0)=img_in.get_roi().at(0);
        roi_.at(1)=img_in.get_roi().at(1);
        roi_.at(2)=img_in.get_roi().at(2);
        roi_.at(3)=img_in.get_roi().at(3);
        roi_.at(4)=img_in.get_roi().at(4);
        roi_.at(5)=img_in.get_roi().at(5);
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> template<class U> void image<T>::hsv_to_rgb(image<U> &image)
 * \brief function to get colour model RGB from HSV (unsigned char version)
 * \author Carlos Falcon
*/
template <class T> template<class U>
void image<T>::hsv_to_rgb(image<U> &image)
{
  if (nChannels_!=3) return;
  int size=width_*height_;
  for (int i=0; i<size;i++)
  {
    U hue=  image_[i];
    U value= image_[i+size];
    U saturation= image_[i+2*size];
    U red=0,blue=0,green=0;

    if(hue<=43){
    red=image[i]=value;
    image[i+size*2]=value-saturation*value/255;
    blue=image[i+size]=hue*(red-blue)/43+blue;
    }
    if(hue>43 && hue<=85){
    green=image[i+size]=value;
    blue=image[i+size*2]=value-saturation*value/255;
    image[i]=-((hue-85)*(green-blue)/43-blue);
    }
    if(hue>85 && hue<=128){
    red=image[i+size]=value;
    green=image[i]=value-saturation*value/255;
    image[i+size*2]=(hue-85)*(green-red)/43+red;
    }
    if(hue>128 && hue<=171){
    blue=image[i+2*size]=value;
    red=image[i]=value-saturation*value/255;
    image[i+size]=((hue-171)*(blue-red)/43+red);
    }
    if(hue>171 && hue<=214){
    blue=image[i+2*size]=value;
    green=image[i+size]=value-saturation*value/255;
    image[i]=-(hue-171)*(blue-green)/43+green;
    }
    if(hue>214){
    red=image[i]=value;
    green=image[i+size]=value-saturation*value/255;
    image[i+2*size]=((hue-214)*(red-green)/43+green);
    }
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> void image<T>::rgb_to_hsv(image<unsigned char> &H,image<unsigned char> &S, image<unsigned char> &V)
 * \brief function to get colour model HSV from RGB (unsigned char version)
 * \author Carlos Falcon
*/
template <class T>
void image<T>::rgb_to_hsv(image<unsigned char> &H,image<unsigned char> &S, image<unsigned char> &V)
{
  if (nChannels_!=3) return;
  if (H.width()!=width_ || S.width()!=width_ || V.width()!=width_)return;

  int size_=width_*height_;
  for (int i=0; i<size_;i++)
  {
    unsigned char red, green, blue;
    red= (unsigned char) image_[i];
    green=(unsigned char) image_[i+size_];
    blue=(unsigned char) image_[i+2*size_];
    unsigned char rgb_min,rgb_max;
    rgb_min=rgb_max=red;
    if(rgb_min>green) rgb_min=green;
    if(rgb_min>blue) rgb_min=blue;
    if(rgb_max<green) rgb_max=green;
    if(rgb_max<blue) rgb_max=blue;
    if(rgb_max == 0){
      H[i]=0;
      S[i]=0;
      V[i]=0;  //HSV
      continue;
    }
    V[i]=(unsigned char) rgb_max;   //V
    S[i]=(unsigned char) ((int) 255*(rgb_max - rgb_min)/rgb_max);  //S
    if(S[i]==0){
      H[i]==0; //H
      continue;
    }
      //H
    if(rgb_max==red){
      H[i]=(unsigned char) ((int) 0+43*((int) green - blue)/(rgb_max - rgb_min));
    }
    else if (rgb_max == green) {
      H[i]= (unsigned char) ((int) 85+ 43*((int) blue-red)/(rgb_max - rgb_min));
    }
    else{
      H[i]= (unsigned char) ((int) 171 + 43*((int) red - green)/(rgb_max - rgb_min));
    }
   }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> void image<T>::set_roi(const int x0,const int x1,const int y0,const int y1)
 * \brief function to set a subwindow border
 * \author Luis Alvarez
*/
template <class T>
void image<T>::set_roi(const int x0,const int x1,const int y0,const int y1)
{
  if(x0>x1 || y0>y1 || x0>width_ || y0>height_) return;

  roi_.resize(6);
  roi_.at(0)=x0;
  roi_.at(1)=x1>width_?width_:x1;
  roi_.at(2)=y0;
  roi_.at(3)=y1>height_?height_:y1;
  roi_.at(4)=0;
  roi_.at(5)=nChannels_;

}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn void image::set_roi(const int x0,const int x1,const int y0,
   const int y1,const int c0,const int c1)
 * \brief function to set a subwindow border
 * \author Luis Alvarez
*/
template <class T>
void image<T>::set_roi(const int x0,const int x1,const int y0,
                          const int y1,const int c0, const int c1)
{
  if(x0>x1 || y0>y1 || x0>width_ || y0>height_ || c0>nChannels_ || c0>c1 ) return;

  roi_.resize(6);
  roi_.at(0)=x0;
  roi_.at(1)=x1>width_?width_:x1;
  roi_.at(2)=y0;
  roi_.at(3)=y1>height_?height_:y1;
  roi_.at(4)=c0;
  roi_.at(5)=c1;

}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> void image<T>::set_roi(const std::vector<int> roi)
 * \brief function to set a subwindow border
 * \author Luis Alvarez
*/
template <class T>
void image<T>::set_roi(const std::vector<int> roi){  roi_=roi; }

/////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> template <class U> void  image<T>::set_roi_image(const image<U> &image2)
 * \brief function to fill the image with a subwindow image
 * \author Luis Alvarez
*/

template <class T> template <class U> void  image<T>::set_roi_image(const image<U> &image2)
{

  if(roi_.size()<6) return;
  int x0=roi_.at(0);
  int x1=roi_.at(1);
  int y0=roi_.at(2);
  int y1=roi_.at(3);
  int c0=roi_.at(4);
  int c1=roi_.at(5);

  if( x1<=width_ &&  y1<=height_ && x0<=x1 && y0<=y1 && c1<=nChannels_ && c0<=c1){
    int width=x1-x0;
    int height=y1-y0;
    int size=width*height;

    if(image2.width()!=width || image2.height()!=height || image2.nChannels()!=nChannels_) return;

    int i,j,n,i2,i2_,n2,n2_,size_=width_*height_;

    #pragma omp parallel \
    shared(x0,x1,y0,y1,c0,c1,width,size,size_) \
    private(n,i,j,n2,n2_,i2,i2_)
    for(n=c0;n<c1;n++){
      n2=(n-c0)*size;
      n2_=n*size_;
      #pragma omp for nowait
      for(i=y0;i<y1;++i){
        i2_=n2_+i*width_;
        i2=n2+(i-y0)*width-x0;
        for(j=x0;j<x1;++j){
          image_[i2_+j]=(T) image2[i2+j] ;
        }
      }
    }
    return;
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T>  const image<T>::get_roi_image(const std::vector <int> &roi2) const
 * \brief function to get a subwindow image
 * \author Luis Alvarez
*/
template <class T>
image<T>  const image<T>::get_roi_image(const std::vector <int> &roi2) const
{
  int x0,x1,y0,y1,c0,c1;

  // IF THERE IS NOT ROI DEFINED WE RETURN AN EMPTY IMAGE
  if(roi2.size()<6){
     return(image<T>());
  }
  else{
    x0=roi2.at(0);
    x1=roi2.at(1);
    y0=roi2.at(2);
    y1=roi2.at(3);
    c0=roi2.at(4);
    c1=roi2.at(5);
  }

  int width,height;
  if( x1<= width_ &&  y1<= height_ && x0<x1 && y0<y1 &&  c1<= nChannels_ && c0<c1){
    width=x1-x0;
    height=y1-y0;
  }
  else{
    return(image<T>());
  }

  image<T> image2(width,height,c1-c0);

  int size_=width_*height_;
  int size=width*height;
  int i,n,j,i2,i2_,n2,n2_;


  #pragma omp parallel \
  shared(image2,x0,x1,y0,y1,c0,c1,width,size,size_) \
  private(n,i,j,n2,n2_,i2,i2_)
  for(n=c0;n<c1;n++){
    n2=(n-c0)*size;
    n2_=n*size_;
    #pragma omp for nowait
    for(i=y0;i<y1;++i){
      i2_=n2_+i*width_;
      i2=n2+(i-y0)*width-x0;
      for(j=x0;j<x1;++j){
        image2[i2+j]= image_[i2_+j] ;
      }
    }
  }
  return(image2);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> template <class U> void  image<T>::get_roi_image(image<U> &image2) const
 * \brief function to fill a subwindow image
 * \author Luis Alvarez
*/

template <class T> template <class U> void  image<T>::get_roi_image(image<U> &image2) const
{
  int x0,x1,y0,y1,c0,c1;

  // IF THERE IS NOT ROI DEFINED WE TAKE THE WHOLE IMAGE
  if(roi_.size()<6){
    x0=0;
    x1=width_;
    y0=0;
    y1=height_;
    c0=0;
    c1=nChannels_;
  }
  else{
    x0=roi_.at(0);
    x1=roi_.at(1);
    y0=roi_.at(2);
    y1=roi_.at(3);
    c0=roi_.at(4);
    c1=roi_.at(5);
  }

  if( x1<= width_ &&  y1<= height_ && x0<=x1 && y0<=y1 &&  c1<= nChannels_ && c0<=c1){
    int width=x1-x0;
    int height=y1-y0;

    // WE ALLOCATE MEMORY IF NEEDED
    if(image2.width()!=width || image2.height()!=height || image2.nChannels()!=nChannels_){
      image2=image<U>(width,height,c1-c0);
    }

    int size_=width_*height_;
    int size=width*height;
    int i,n,j,i2,i2_,n2,n2_;


    #pragma omp parallel \
    shared(x0,x1,y0,y1,c0,c1,width,size,size_) \
    private(n,i,j,n2,n2_,i2,i2_)
    for(n=c0;n<c1;n++){
      n2=(n-c0)*size;
      n2_=n*size_;
      #pragma omp for nowait
      for(i=y0;i<y1;++i){
        i2_=n2_+i*width_;
        i2=n2+(i-y0)*width-x0;
        for(j=x0;j<x1;++j){
          image2[i2+j]=(U) image_[i2_+j] ;
        }
      }
    }
    return;
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T> & image<T>::operator=(const image &image2)
 * \brief operator = (equality of images of different size is not allowed )
 * \author Luis Alvarez
*/
template <class T>
image<T> & image<T>::operator=(const image &image2)
{
  if( width_!=image2.width() || height_!=image2.height() || nChannels_!=image2.nChannels()){
    width_=image2.width();
    height_=image2.height();
    nChannels_=image2.nChannels();
    image_.resize(image2.size());
  }

  int k,k_end=image2.size();
  #pragma omp parallel for shared(k_end) private(k)
  for(k=0;k<k_end;++k)
    image_[k]=  image2[k];

  return *this;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T> const image<T>::linear_transform(const double H[3][3]) const ;
 * \brief function to apply a linear transform given by an homography H to the image. H goes to the output image to the input image
 * \author Luis Alvarez
*/
template <class T>
image<T> const image<T>::linear_transform(const double H[3][3]) const
{
  image<T> image2(width_,height_,nChannels_,255);

  int width_1=width_-1;
  int height_1=height_-1;
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    #pragma omp parallel for shared(H,image2,n,n0,width_1,height_1)
    for(int i=0;i<height_;i++){
      for(int j=0;j<width_;j++){
        //WE COMPUTE THE HOMOGRAPHY TRANSFORMATION
        double z0=H[2][0]*j+H[2][1]*i+H[2][2];
        double x0=H[0][0]*j+H[0][1]*i+H[0][2];
        double y0=H[1][0]*j+H[1][1]*i+H[1][2];
        x0/=z0;
        y0/=z0;
        if(x0<0 || x0>=width_ || y0<0 || y0>=height_) continue;
        int y0i=y0;
        double dy0=y0-y0i;
        double dy0_1=1-dy0;
        int m0=n0+y0i*width_;
        int m2=n0+i*width_;
        int x0i=x0;
        double dx0=x0-x0i;
        double dx0_1=1-dx0;
        if( x0i<width_1){
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
        }
        else{
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
        }
      }
    }
  }
  return image2;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T> const image<T>::resize(const int width, const int height) const
 * \brief function to resize an image
 * \author Luis Alvarez
*/
template <class T>
image<T> const image<T>::resize(
const int width /** new image width */,
const int height /** new image height */) const
{
  printf("width=%d width_=%d height=%d height_=%d \n",width,width_,height,height_);

  if(width==width_ && height==height_) return (*this);

  if(width<=width_/2 && height<=height_/2){
     return(((*this).sampling()).resize(width,height));
  }

  if(height<100) return(resize_no_omp(width,height));

  if(width==0 || height==0) return(image<T> ());

  image<T> image2(width,height,nChannels_);


  // WE PROCESS IN THE CASE WE DO NOT NEED INTERPOLATION
  if( width_>width && width_%width==0  && height_>height && height_%height==0 ){
    int scale_x=width_/width;
    int scale_y=height_/height;
    for(int n=0;n<nChannels_;n++){
      int n0=n*width_*height_;
      int n2=n*width*height;
      #pragma omp parallel for shared(image2,n,n0,n2,scale_x,scale_y)
      for(int i=0;i<height;i++){
        int m0=n0+scale_y*i*width_;
        int m2=n2+i*width;
        for(int j=0;j<width;j++){
          image2[m2+j]=image_[m0+scale_x*j];
        }
      }
    }
    return image2;
  }

  // WE PROCESS THE CASE WHERE WE NEED A INTERPOLATION PROCEDURE
  double scale_x=(double) width_/width;
  double scale_y=(double) height_/height;

  int width_1=width_-1;
  int height_1=height_-1;
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    int n2=n*width*height;
    #pragma omp parallel for shared(image2,n,n0,n2,scale_x,scale_y,width_1,height_1)
    for(int i=0;i<height;i++){
      double y0=scale_y*i;
      int y0i=y0;
      double dy0=y0-y0i;
      double dy0_1=1-dy0;
      int m0=n0+y0i*width_;
      int m2=n2+i*width;
      for(int j=0;j<width;j++){
        double x0=scale_x*j;
        int x0i=x0;
        double dx0=x0-x0i;
        double dx0_1=1-dx0;
        if( x0i<width_1){
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
        }
        else{
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
        }
      }
    }
  }

  return image2;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T> const image<T>::resize_no_omp(const int width, const int height) const
 * \brief function to resize an image without omp
 * \author Luis Alvarez
*/
template <class T>
image<T> const image<T>::resize_no_omp(const int width, const int height) const
{

  if(width==0 || height==0) return(image<T> ());

  if(width==width_ && height==height_) return (*this);

  if(width<=width_/2 && height<=height_/2) return(((*this).sampling()).resize_no_omp(width,height));

  image<T> image2(width,height,nChannels_);


  // WE PROCESS IN THE CASE WE DO NOT NEED INTERPOLATION
  if( width_>width && width_%width==0  && height_>height && height_%height==0 ){
    int scale_x=width_/width;
    int scale_y=height_/height;
    for(int n=0;n<nChannels_;n++){
      int n0=n*width_*height_;
      int n2=n*width*height;
      for(int i=0;i<height;i++){
        int m0=n0+scale_y*i*width_;
        int m2=n2+i*width;
        for(int j=0;j<width;j++){
          image2[m2+j]=image_[m0+scale_x*j];
        }
      }
    }
    return image2;
  }

  // WE PROCESS THE CASE WHERE WE NEED A INTERPOLATION PROCEDURE
  double scale_x=(double) width_/width;
  double scale_y=(double) height_/height;

  int width_1=width_-1;
  int height_1=height_-1;
  for(int n=0;n<nChannels_;n++){
    int n0=n*width_*height_;
    int n2=n*width*height;
    for(int i=0;i<height;i++){
      double y0=scale_y*i;
      int y0i=y0;
      double dy0=y0-y0i;
      double dy0_1=1-dy0;
      if(y0i>=height_) y0i=height_-1;
      int m0=n0+y0i*width_;
      int m2=n2+i*width;
      for(int j=0;j<width;j++){
        double x0=scale_x*j;
        int x0i=x0;
        double dx0=x0-x0i;
        double dx0_1=1-dx0;
        if( x0i<width_1){
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i+1]+
              dy0_1*dx0*image_[m0+x0i+1]
            );
          }
        }
        else{
          if( y0i<height_1){
            image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+width_+x0i]+
              dy0*dx0*image_[m0+width_+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
          else{
              image2[m2+j]=(T) (
              dy0_1*dx0_1*image_[m0+x0i]+
              dy0*dx0_1*image_[m0+x0i]+
              dy0*dx0*image_[m0+x0i]+
              dy0_1*dx0*image_[m0+x0i]
            );
          }
        }
      }
    }
  }
  return image2;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_bool(std::string name )
 * \brief function to write to disk a 8/16 bit tif image from boolean image
 * \author Carlos Falcon
*/
template <class T>
int image<T>::write_bool(std::string name /** image file name */)
{
  ami::image<unsigned char> aux(width_,height_,nChannels_,0);
  int size=width_*height_;
  #pragma omp parallel for shared(size,aux)
  for (int i=0; i<size;i++)
  {
    if (get_value(i))
      aux[i]=255;
  }
  return aux.write(name);
}

////////////////////////////////////////////////////////////////////////////////////////////////
#ifdef __AMI_OPENCV__
/**
 * \fn template <class T> int image<T>::write_with_opencv( std::string name )
 * \brief Function to write an image to disk with opencv.
 * \param [in] name Filename of the image.
 * \author Pedro Henriquez
*/
template <class T>
int image<T>::write_with_opencv(std::string name /** image file name */)
{
  //DETERMINE IMAGE DATA TYPE
  int tipo;
  if (typeid( image_[0] ) == typeid( unsigned char )) tipo = IPL_DEPTH_8U;
  if (typeid( image_[0] ) == typeid( char ))         tipo = IPL_DEPTH_8S;
  if (typeid( image_[0] ) == typeid( unsigned short )) tipo= IPL_DEPTH_16U;
  if (typeid( image_[0] ) == typeid( short )) tipo= IPL_DEPTH_16U;
  if (typeid( image_[0] ) == typeid( int )) tipo = IPL_DEPTH_32S;
  if (typeid( image_[0] ) == typeid( float )) tipo = IPL_DEPTH_32F;
  if (typeid( image_[0] ) == typeid( double )) tipo = IPL_DEPTH_64F;

  //CREATE OPENCV IPLIMAGE
  IplImage* image2 = cvCreateImage(cvSize(width_,height_), tipo, nChannels_);

  if (image2 == NULL) return -1;

  long cont=0;
  long size_=width_*height_;

  if(nChannels_==3) //CONVERT RGB TO BGR
  {
  for(int i=height_-1;i>=0;i--)
  {
    for(int j=0;(int) j<width_;j++)
    {
      long m=i*width_+j;
      image2->imageData[cont++]= image_[m+2*size_]; //B
      image2->imageData[cont++]= image_[m+size_]; //G
      image2->imageData[cont++]= image_[m]; //R
    }
    int incr = 0;
    if((3*width_)%4 != 0) incr = (4-(3*width_)%4);
    for(int k=0; (int) k<incr; k++) image2->imageData[cont++]=0;
  }
  }
  else
  {
  for(int i=height_-1;i>=0;i--)
  {
    for(int j=0;(int)j<width_;j++)
    {
      long m=i*width_+j;
      for(int nc=0;(int) nc<nChannels_;nc++)
      {
        image2->imageData[cont++]= image_[m+nc*size_];
      }
    }
    int incr = 0;
    if((nChannels_*width_)%4 != 0) incr = (4-(nChannels_*width_)%4);
    for(int k=0; (int) k<incr; k++) image2->imageData[cont++]=0;
  }
  }

  cvSaveImage(name.c_str(),image2);
  cvReleaseImage (&image2);
  return 0;
}
#else

#ifdef __AMI_IMAGE__


#else
////////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_with_jpeglib( std::string name )
 * \brief Function to write an image to disk with jpeglib.
 * \param [in] name Filename of the image.
 * \author Pedro Henriquez
*/
template <class T>
int image<T>::write_with_jpeglib(std::string name /** image file name */)
{
  struct jpeg_compress_struct cinfo;
	struct jpeg_error_mgr jerr;
	cinfo.err = jpeg_std_error(&jerr);
	jpeg_create_compress(&cinfo);

	FILE * outfile;
	if ((outfile = fopen(name.c_str(), "wb")) == NULL) {
	  #ifdef __AMIDEBUG__
	  printf("can't open %s\n", name.c_str());
    #endif
	  return(-1);
	}
	jpeg_stdio_dest(&cinfo, outfile);

  cinfo.image_width = width_; 	/* image width and height, in pixels */
	cinfo.image_height = height_;
	cinfo.input_components = nChannels_;	/* # of color components per pixel */
	cinfo.in_color_space = JCS_RGB; /* colorspace of input image */
	jpeg_set_defaults(&cinfo);

  jpeg_start_compress(&cinfo, TRUE);

  JSAMPROW row_pointer[1];	/* pointer to a single row */
//	int row_stride;			/* physical row width in buffer */

  //CONSTRUCT ROW POINTERS
  JSAMPLE *row_pointers;
  row_pointers = (JSAMPLE*) malloc(sizeof(JSAMPLE) * width_ * nChannels_);

  long cont;
  long size_=width_*height_;
  if (typeid( image_[0] ) == typeid( bool ))
  {
    if(nChannels_==3)//CONVERT RGB BY PIXEL IN RGB CHANNELS
    {
      for(long i=0;i<height_;i++)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[cont++]=(JSAMPLE)(image_[m])?255:0;//R
          row_pointers[cont++]=(JSAMPLE)(image_[m+size_])?255:0;//G
          row_pointers[cont++]=(JSAMPLE)(image_[m+2*size_])?255:0; //B
        }
        row_pointer[0] =row_pointers;
        jpeg_write_scanlines(&cinfo, row_pointer, 1);
      }
    }
    else
    {
      for(long i=0;i<height_;i++)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          for(long n=0;n<(long) nChannels_;n++)
          {
            row_pointers[cont++]=(image_[m+n*size_])?255:0;
          }
        }
        jpeg_write_scanlines(&cinfo, row_pointer, 1);
      }
    }
  }
  else
  {
    if(nChannels_==3)//CONVERT RGB BY PIXEL IN RGB CHANNELS
    {
      for(long i=0;i<height_;i++)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[cont++]=(JSAMPLE)image_[m];//R
          row_pointers[cont++]=(JSAMPLE)image_[m+size_];//G
          row_pointers[cont++]=(JSAMPLE)image_[m+2*size_]; //B
        }
        row_pointer[0] =row_pointers;
        jpeg_write_scanlines(&cinfo, row_pointer, 1);
      }
    }
    else
    {
      for(long i=0;i<height_;i++)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          for(long n=0;n<(long) nChannels_;n++)
          {
            row_pointers[cont++]=image_[m+n*size_];
          }
        }
        jpeg_write_scanlines(&cinfo, row_pointer, 1);
       }
     }
  }

	jpeg_finish_compress(&cinfo);
	fclose(outfile);
	free (row_pointers);
	return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_with_pnglib( std::string name )
 * \brief Function to write an image to disk with pnglib.
 * \param [in] name Filename of the image.
 * \author Pedro Henriquez
*/
template <class T>
int image<T>::write_with_pnglib(std::string name /** image file name */)
{
  // SELECT BITs PER CHANNEL
  char bit_depth=16;
  if (typeid( image_[0] ) == typeid( unsigned char ))      bit_depth = 8;
  if (typeid( image_[0] ) == typeid( char ))               bit_depth = 8;
  if (typeid( image_[0] ) == typeid( bool ))               bit_depth = 8;
  if (typeid( image_[0] ) == typeid( short int ))          bit_depth = 16;
  if (typeid( image_[0] ) == typeid( unsigned short int )) bit_depth = 16;
  if (typeid( image_[0] ) == typeid( int ))                bit_depth = 16;
  if (typeid( image_[0] ) == typeid( float ))              bit_depth = 16;
  if (typeid( image_[0] ) == typeid( double ))             bit_depth = 16;
  if (typeid( image_[0] ) == typeid( int ))       bit_depth = 16;
  if (bit_depth== 8)return write_with_pnglib_08(name);
  if (bit_depth==16)return write_with_pnglib_16(name);
//  if (bit_depth==32)return write_with_pnglib_32(name);
  return -1;
}
//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_with_pnglib_08( std::string name )
 * \brief Function to write an 8 bits per channel image to disk with pnglib.
 * \param [in] name Filename of the image.
 * \author Jesus Rodriguez
*/
template <class T>
int image<T>::write_with_pnglib_08(std::string name /** image file name */)
{
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: SELECT COLOR \n");
  #endif
  // SELECT COLOR
  char color_type;
  switch (nChannels_)
  {
    case 0:
      return -1;
    case 1:
    case 2:
      color_type = PNG_COLOR_TYPE_GRAY;
      break;
    default:
      color_type = PNG_COLOR_TYPE_RGB;
      break;
  }
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: SELECT BITS PER CHANNEL \n");
  #endif
  // SELECT BITs PER CHANNEL
  char bit_depth=8;
//  if (typeid( image_[0] ) == typeid( unsigned char ))   bit_depth = 8;
//  if (typeid( image_[0] ) == typeid( char ))            bit_depth = 8;
  if (typeid( image_[0] ) == typeid( bool )){           bit_depth = 1; color_type = PNG_COLOR_TYPE_GRAY;}
  if (bit_depth==16){return write_with_pnglib_16(name);}
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: CREATE FILE \n");
  #endif
  // CREATE FILE
	FILE *imf = fopen(name.c_str(), "wb");
	if (!imf)
	{
    #ifdef __AMIDEBUG__
      printf("File %s could not be opened for writing", name.c_str());
    #endif
    return -1;
	}
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: INITIALIZE STUFF \n");
  #endif
	// INITIALIZE STUFF
	png_structp png_ptr;
  png_infop info_ptr;
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

	if (!png_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_write_struct failed");
    #endif
    return -2;
	}

	info_ptr = png_create_info_struct(png_ptr);
	if (!info_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_info_struct failed");
    #endif
    return -3;
	}

	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during init_io");
    #endif
    return -4;
	}
	png_init_io(png_ptr, imf);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: WRITE HEADER \n");
  #endif

	// WRITE HEADER
	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during writing header");
    #endif
    return -5;
	}
	png_set_IHDR(png_ptr, info_ptr, width_, height_,
		     bit_depth, color_type, PNG_INTERLACE_NONE,
		     PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: WRITE BYTES \n");
  #endif

	// WRITE BYTES
	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during writing bytes");
    #endif
    return -6;
	}
  //CONSTRUCT ROW POINTERS
  png_bytep *row_pointers;
  row_pointers = (png_bytep*) malloc(sizeof(png_bytep) * height_);
  if (row_pointers == NULL) return -7;
  int data_sizeb = sizeof(png_bytep) * width_ * ((nChannels_<3)?1:3);
	for (int y=0; y<height_; y++)
	{
		row_pointers[y] = (png_byte*) malloc(data_sizeb);
		if (row_pointers[y]==NULL) return -8;
	}
  long cont = 0;
  long size_=width_*height_;
	if(nChannels_>2)//CONVERT RGB BY PIXEL IN RGB CHANNELS
  {
    if (typeid( image_[0] ) == typeid( bool ))
    {
      for(long i=height_-1;i>=0;i--)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[i][cont++]=(png_byte)(image_[m])?255:0;//R
          row_pointers[i][cont++]=(png_byte)(image_[m+size_])?255:0;//G
          row_pointers[i][cont++]=(png_byte)(image_[m+2*size_])?255:0; //B
        }
      }
    }
    else
    {
      for(long i=height_-1;i>=0;i--)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[i][cont++]=(png_byte)image_[m];//R
          row_pointers[i][cont++]=(png_byte)image_[m+size_];//G
          row_pointers[i][cont++]=(png_byte)image_[m+2*size_]; //B
        }
      }
    }
  }
  else
  {
    //int numChannels = 1;
    if (typeid( image_[0] ) == typeid( bool ))
    {
      for(long i=height_-1;i>=0;i--)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[i][cont++]=(png_byte)(image_[m])?255:0;
        }
      }
    }
    else
    {
      for(long i=height_-1;i>=0;i--)
      {
        long m=(height_-1-i)*width_;
        cont = 0;
        for(long j=0;j<(long)width_;j++,m++)
        {
          row_pointers[i][cont++]=(png_byte)image_[m];
        }
      }
    }
  }
	png_write_image(png_ptr, row_pointers);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: END WRITE \n");
  #endif

	// END WRITE
	if (setjmp(png_jmpbuf(png_ptr)))
	{
	  png_destroy_write_struct(&png_ptr, &info_ptr);
    #ifdef __AMIDEBUG__
      printf("Error during end of write");
		#endif
		return -8;
	}
	png_write_end(png_ptr, NULL);

  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: CLEANUP HEAP ALLOCATION \n");
  #endif
  //CLEANUP HEAP ALLOCATION
	for (int y=0; y<height_; y++)free(row_pointers[y]);
	free(row_pointers);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_08: fclose \n");
  #endif
  png_destroy_write_struct(&png_ptr, &info_ptr);
  fclose(imf);
  return 0;
}

/**
 * \fn template <class T> int image<T>::write_with_pnglib_16( std::string name )
 * \brief Function to write an 16 bits per channel image to disk with pnglib.
 * \param [in] name Filename of the image.
 * \author Jesus Rodriguez
*/

template <class T>
int image<T>::write_with_pnglib_16(std::string name /** image file name */)
{
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: SELECT COLOR \n");
  #endif
  // SELECT COLOR
  char color_type;
  switch (nChannels_)
  {
    case 0:
      return -1;
    case 1:
    case 2:
      color_type = PNG_COLOR_TYPE_GRAY;
      break;
    default:
      color_type = PNG_COLOR_TYPE_RGB;
      break;
  }
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: SELECT BITS PER CHANNEL \n");
  #endif
  // SELECT BITs PER CHANNEL
  char bit_depth=16;
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib: CREATE FILE \n");
  #endif
  // CREATE FILE
	FILE *imf = fopen(name.c_str(), "wb");
	if (!imf)
	{
    #ifdef __AMIDEBUG__
      printf("File %s could not be opened for writing", name.c_str());
    #endif
    return -1;
	}
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: INITIALIZE STUFF \n");
  #endif
	// INITIALIZE STUFF
	png_structp png_ptr;
  png_infop info_ptr;
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

	if (!png_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_write_struct failed");
    #endif
    return -2;
	}

	info_ptr = png_create_info_struct(png_ptr);
	if (!info_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_info_struct failed");
    #endif
    return -3;
	}

	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during init_io");
    #endif
    return -4;
	}

	png_init_io(png_ptr, imf);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: WRITE HEADER \n");
  #endif

	// WRITE HEADER
	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during writing header");
    #endif
    return -5;
	}
	png_set_IHDR(png_ptr, info_ptr, width_, height_,
		     bit_depth, color_type, PNG_INTERLACE_NONE,
		     PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: WRITE BYTES \n");
  #endif

	// WRITE BYTES
	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during writing bytes");
    #endif
    return -6;
	}
  //CONSTRUCT ROW POINTERS
  png_uint_16 **row_pointers;
  row_pointers = (png_uint_16**) malloc(sizeof(png_uint_16*) * height_);
  if (row_pointers == NULL) return -7;
  int data_sizeb = sizeof(png_uint_16) * width_ * ((nChannels_<3)?1:3);
	for (int y=0; y<height_; y++){
		row_pointers[y] = (png_uint_16*) malloc(data_sizeb);
		if (row_pointers[y]==NULL) return -8;
	}
  long cont = 0;
  long size_=width_*height_;

  unsigned char *data_in, *data_out;
  T info_in, info_out;
  data_in  = (unsigned char*)&info_in;
  data_out = (unsigned char*)&info_out;
	if(nChannels_>2)//CONVERT RGB BY PIXEL IN RGB CHANNELS
  {
    for(long i=height_-1;i>=0;i--)
    {
      long m=(height_-1-i)*width_;
      cont = 0;
      for(long j=0;j<(long)width_;j++,m++)
      {
        info_in = image_[m]/* *256*/;
        data_out[0] = data_in[1];
        data_out[1] = data_in[0];
        row_pointers[i][cont++]=(png_uint_16)info_out;//R
        info_in = image_[m+size_]/* *256*/;
        data_out[0] = data_in[1];
        data_out[1] = data_in[0];
        row_pointers[i][cont++]=(png_uint_16)info_out;//G
        info_in = image_[m+2*size_]/* *256*/;
        data_out[0] = data_in[1];
        data_out[1] = data_in[0];
        row_pointers[i][cont++]=(png_uint_16)info_out; //B
      }
    }
  }
  else
  {
    for(long i=0;i< height_;i++)
    {
      long m=(height_-1-i)*width_;
      cont = 0;
      for(long j=0;j<(long)width_;j++,m++)
      {
        info_in = image_[m]/* *256*/;
        data_out[0] = data_in[1];
        data_out[1] = data_in[0];
        row_pointers[height_-1-i][width_-1-cont]=(png_uint_16)info_out;
        cont++;
      }
    }
  }
	png_write_image(png_ptr, (png_byte**)row_pointers);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: END WRITE \n");
  #endif

	// END WRITE
	if (setjmp(png_jmpbuf(png_ptr)))
	{
	  png_destroy_write_struct(&png_ptr, &info_ptr);
    #ifdef __AMIDEBUG__
      printf("Error during end of write");
		#endif
		return -8;
	}
	png_write_end(png_ptr, NULL);

  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: CLEANUP HEAP ALLOCATION \n");
  #endif
  //CLEANUP HEAP ALLOCATION
	for (int y=0; y<height_; y++)free(row_pointers[y]);
	free(row_pointers);
  #ifdef __AMIDEBUG__
    printf("write_with_pnglib_16: fclose \n");
  #endif
  png_destroy_write_struct(&png_ptr, &info_ptr);
  fclose(imf);
  return 0;
}

//
//template <class T>
//int image<T>::write_with_pnglib_32(std::string name /** image file name */)
//{
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: SELECT COLOR \n");
//  #endif
//  // SELECT COLOR
//  char color_type;
//  switch (nChannels_)
//  {
//    case 0:
//      return -1;
//    case 1:
//    case 2:
//      color_type = PNG_COLOR_TYPE_GRAY;
//      break;
//    default:
//      color_type = PNG_COLOR_TYPE_RGB;
//      break;
//  }
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: SELECT BITS PER CHANNEL \n");
//  #endif
//  // SELECT BITs PER CHANNEL
//  char bit_depth=32;
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: CREATE FILE \n");
//  #endif
//  // CREATE FILE
//	FILE *imf = fopen(name.c_str(), "wb");
//	if (!imf)
//	{
//    #ifdef __AMIDEBUG__
//      printf("File %s could not be opened for writing", name.c_str());
//    #endif
//    return -1;
//	}
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: INITIALIZE STUFF \n");
//  #endif
//	// INITIALIZE STUFF
//	png_structp png_ptr;
//  png_infop info_ptr;
//	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
//
//	if (!png_ptr)
//	{
//    #ifdef __AMIDEBUG__
//      printf("png_create_write_struct failed");
//    #endif
//    return -2;
//	}
//
//	info_ptr = png_create_info_struct(png_ptr);
//	if (!info_ptr)
//	{
//    #ifdef __AMIDEBUG__
//      printf("png_create_info_struct failed");
//    #endif
//    return -3;
//	}
//
//	if (setjmp(png_jmpbuf(png_ptr)))
//	{
//    #ifdef __AMIDEBUG__
//      printf("Error during init_io");
//    #endif
//    return -4;
//	}
//
//	png_init_io(png_ptr, imf);
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: WRITE HEADER \n");
//  #endif
//
//	// WRITE HEADER
//	if (setjmp(png_jmpbuf(png_ptr)))
//	{
//    #ifdef __AMIDEBUG__
//      printf("Error during writing header");
//    #endif
//    return -5;
//	}
//	png_set_IHDR(png_ptr, info_ptr, width_, height_,
//		     bit_depth, color_type, PNG_INTERLACE_NONE,
//		     PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
//
//	png_write_info(png_ptr, info_ptr);
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: WRITE BYTES \n");
//  #endif
//
//	// WRITE BYTES
//	if (setjmp(png_jmpbuf(png_ptr)))
//	{
//    #ifdef __AMIDEBUG__
//      printf("Error during writing bytes");
//    #endif
//    return -6;
//	}
//  //CONSTRUCT ROW POINTERS
//  png_int_32 **row_pointers;
//  row_pointers = (png_int_32**) malloc(sizeof(png_int_32*) * height_);
//  if (row_pointers == NULL) return -7;
//  int data_sizeb = sizeof(png_int_32) * width_ * ((nChannels_<3)?1:3);
//	for (int y=0; y<height_; y++){
//		row_pointers[y] = (png_int_32*) malloc(data_sizeb);
//		if (row_pointers[y]==NULL) return -8;
//	}
//  long cont = 0;
//  long size_=width_*height_;
//
//  unsigned char *data_in, *data_out;
//  T info_in, info_out;
//  data_in  = (unsigned char*)&info_in;
//  data_out = (unsigned char*)&info_out;
//	if(nChannels_>2)//CONVERT RGB BY PIXEL IN RGB CHANNELS
//  {
//    for(long i=height_-1;i>=0;i--)
//    {
//      long m=(height_-1-i)*width_;
//      cont = 0;
//      for(long j=0;j<(long)width_;j++,m++)
//      {
//        info_in = image_[m]*10000;
//        data_out[0] = data_in[3];
//        data_out[1] = data_in[2];
//        data_out[2] = data_in[1];
//        data_out[3] = data_in[0];
//        row_pointers[i][cont++]=(png_int_32)info_out;//R
//        info_in = image_[m+size_]*10000;
//        data_out[0] = data_in[3];
//        data_out[1] = data_in[2];
//        data_out[2] = data_in[1];
//        data_out[3] = data_in[0];
//        row_pointers[i][cont++]=(png_int_32)info_out;//G
//        info_in = image_[m+2*size_]*10000;
//        data_out[0] = data_in[3];
//        data_out[1] = data_in[2];
//        data_out[2] = data_in[1];
//        data_out[3] = data_in[0];
//        row_pointers[i][cont++]=(png_int_32)info_out; //B
//      }
//    }
//  }
//  else
//  {
//    for(long i=0;i< height_;i++)
//    {
//      long m=(height_-1-i)*width_;
//      cont = 0;
//      for(long j=0;j<(long)width_;j++,m++)
//      {
//        info_in = image_[m]*10000;
//        data_out[0] = data_in[3];
//        data_out[1] = data_in[2];
//        data_out[2] = data_in[1];
//        data_out[3] = data_in[0];
//        row_pointers[height_-1-i][width_-1-cont]=(png_int_32)info_out;
//        cont++;
//      }
//    }
//  }
//	png_write_image(png_ptr, (png_byte**)row_pointers);
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: END WRITE \n");
//  #endif
//
//	// END WRITE
//	if (setjmp(png_jmpbuf(png_ptr)))
//	{
//	  png_destroy_write_struct(&png_ptr, &info_ptr);
//    #ifdef __AMIDEBUG__
//      printf("Error during end of write");
//		#endif
//		return -8;
//	}
//	png_write_end(png_ptr, NULL);
//
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: CLEANUP HEAP ALLOCATION \n");
//  #endif
//  //CLEANUP HEAP ALLOCATION
//	for (int y=0; y<height_; y++)free(row_pointers[y]);
//	free(row_pointers);
//  #ifdef __AMIDEBUG__
//    printf("write_with_pnglib_32: fclose \n");
//  #endif
//  png_destroy_write_struct(&png_ptr, &info_ptr);
//  fclose(imf);
//  return 0;
//}
//
//

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_with_tifflib( std::string name )
 * \brief Function to write an image to disk with tifflib.
 * \param [in] name Filename of the image.
 * \author Pedro Henriquez
*/
template <class T>
int image<T>::write_with_tifflib(std::string name /** image file name */)
{
  int width, height, page, pages = 1;

  width =  width_;
  height = height_;
  //int scomponents = nChannels_;

  //double resolution = -1;
  //uint32 rowsperstrip = (uint32) -1;
  int bps=16;

  if (typeid( image_[0] ) == typeid( unsigned char ))   bps = 8;
  if (typeid( image_[0] ) == typeid( char ))            bps = 8;
  if (typeid( image_[0] ) == typeid( bool )){           bps = 8;}

 // int predictor;

  TIFF *tif = TIFFOpen(name.c_str(), "w");
  if ( !tif )
  {
#ifdef __AMIDEBUG__
    printf("image<T>::write_with_tifflib cannot open file.\n");
#endif
    return -1;
  }
  if (typeid( image_[0] ) == typeid( char ))       TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( short ))      TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( long ))       TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( int ))        TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( float ))      TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( double ))     TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( long int ))   TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
  if (typeid( image_[0] ) == typeid( long double ))TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);

  uint32 w = width;
  uint32 h = height;

  for (page = 0; page < pages; page++)
  {
    TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, w);
    TIFFSetField(tif, TIFFTAG_IMAGELENGTH, h);
    TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
    TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, nChannels_);
    TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, bps); // Fix for stype
    TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);

    if (typeid( image_[0] ) == typeid( char ))       TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( short ))      TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( long ))       TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( int ))        TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( float ))      TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( double ))     TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( long int ))   TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);
    if (typeid( image_[0] ) == typeid( long double ))TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT,SAMPLEFORMAT_INT);

    uint16 photometric = (nChannels_ == 1) ? PHOTOMETRIC_MINISBLACK : PHOTOMETRIC_RGB;

    TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, photometric); // Fix for scomponents

    int rowLength=0; // in bytes

    if (typeid( image_[0] ) == typeid( unsigned char ))        rowLength = sizeof(unsigned char);
    if (typeid( image_[0] ) == typeid( char ))                 rowLength = sizeof(char);
    if (typeid( image_[0] ) == typeid( unsigned short ))       rowLength = sizeof(unsigned short);
    if (typeid( image_[0] ) == typeid( short ))                rowLength = sizeof(short);
    if (typeid( image_[0] ) == typeid( unsigned long ))        rowLength = sizeof(unsigned long);
    if (typeid( image_[0] ) == typeid( long ))                 rowLength = sizeof(long);
    if (typeid( image_[0] ) == typeid( int ))         rowLength = sizeof(int);
    if (typeid( image_[0] ) == typeid( int ))                  rowLength = sizeof(int);
    if (typeid( image_[0] ) == typeid( float ))                rowLength = sizeof(float);
    if (typeid( image_[0] ) == typeid( double ))               rowLength = sizeof(double);
    if (typeid( image_[0] ) == typeid( unsigned long int ))    rowLength = sizeof(unsigned long int);
    if (typeid( image_[0] ) == typeid( long int ))             rowLength = sizeof(long int);
    if (typeid( image_[0] ) == typeid( long double ))          rowLength = sizeof(long double);
    if (rowLength == 0)
    {
      printf("Unrecognized data type.\n");
      return -2;
    }

    rowLength = nChannels_;
    rowLength *= width;

    T *outPtr;
    outPtr = new T[rowLength];
    unsigned char *data ;
    data = (unsigned char *)outPtr;

    long size_  = width_ * height_;

    int row = 0;
    for (int idx2 = 0; idx2 < height; idx2++)
    {
      long cont = 0;
      long m=(height_-1-idx2)*width_;
      for(long j=0;j<(long)width_;j++,m++)
      {
        outPtr[cont++] = image_[m]/* *256*/;//r
        outPtr[cont++] = image_[m+size_]/* *256*/;//g
        outPtr[cont++] = image_[m+2*size_]/* *256*/;//b
      }
      if ( TIFFWriteScanline(tif, data, row, 0) < 0)
      {
        printf( "TIFFImageIO: error out of disk space" ); fflush(stdout);
        break;
      }
      row++;
    }
    delete[] outPtr;
  }
  TIFFClose(tif);
  return 0;
}
#endif //AMITIFF

/**
* \fn int image<T>::write_with_ami_bmp (std::string name)
 * \brief Function to write an image to disk with ami_bmp.
 * \param [in] name Filename of the image.
* \author Jesus Rodriguez
*/
template <class T>
int image<T>::write_with_ami_bmp (std::string name/** image file name */)
{
  // BUILD NAME
  char name_c[200];
  snprintf(name_c, 199, "%s", name.c_str());
  // ALLOC MEMORY
  int size_ = width_ * height_;
  unsigned char *red, *green, *blue;
  red   = new unsigned char[size_];
  green = new unsigned char[size_];
  blue  = new unsigned char[size_];
  // INITIALIZE INDIXES
  int mr = 0;
  int mg = size_;
  int mb = 2*size_;
  // REFILL DATA VECTORS
  if(nChannels_ < 1) return 0;
  if(nChannels_ > 2)
  {
    for(long i=height_-1;i>-1;i--)
    {
      int m = (height_-1-i)*width_;
      for(long j=0;j<width_;j++,mr++,mg++,mb++,m++)
      {
        red[m]  = (unsigned char)image_[mr];//R
        green[m]= (unsigned char)image_[mg];//G
        blue[m] = (unsigned char)image_[mb];//B
      }
    }
  }
  else
  {
    for(long i=height_-1;i>-1;i--)
    {
      int m = (height_-1-i)*width_;
      for(long j=0;j<width_;j++,mr++,m++)
      {
        red[m]  = green[m]= blue[m] = (unsigned char)image_[mr];
      }
    }
  }
  // WRITE IMAGE
  ami_write_bmp(name_c, red, green, blue, width_, height_);
  // FREE MEMORY
  delete []red;
  delete []green;
  delete []blue;
  return 0;
}
#endif //OPENCV
////////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write_pgm( std::string name )
 * \brief Function to write a pgm image to disk.
 * \param [in] name Filename of the image.
 * \author Javier Martin
*/
template <class T>
int image<T>::write_pgm(std::string name /** image file name */)
{
  // IS PGM
  imageMirrored();  //IN PGM THE IMAGE IS MIRRORED
  FILE *fp;
  fp=fopen(name.c_str(),"wb");
  int size = width_*height_;
  if ((typeid( image_[0] ) == typeid( unsigned char )) || (typeid( image_[0] ) == typeid( char ))){
  fprintf(fp,"P5\n%d %d\n%d\n",width_,height_,255); //ONE BYTE TYPE
  char * im = (char*) &image_[0];
  for (int i=size-1;i>=0;i--){
    fprintf(fp,"%c", im[i]);
  }
  fclose(fp);
  imageMirrored();  //IN PGM THE IMAGE IS MIRRORED
  return 0;
  }
  if ((typeid( image_[0] ) == typeid( unsigned short )) || (typeid( image_[0] ) == typeid( unsigned short ))){
  fprintf(fp,"P5\n%d %d\n%d\n",width_,height_,1023);//TWO BYTE TYPE, 10 BITS MSG IMAGE
  char *im = (char*) &image_[0];
  for (int i=2*size-1;i>=0;i--){
    fprintf(fp,"%c", im[i]);
  }
  fclose(fp);
  imageMirrored(); //IN PGM THE IMAGE IS MIRRORED
  return 0;
  }
  imageMirrored(); //IN PGM THE IMAGE IS MIRRORED
  return -1;
}
////////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> int image<T>::write( std::string name )
 * \brief Function to write an image to disk.
 * \param [in] name Filename of the image.
 * \author Pedro Henriquez
*/
template <class T>
int image<T>::write(std::string name /** image file name */)
{
  if(image_.size()==0) return -1;
  if((strcmp(&name.c_str()[strlen(name.c_str())-3],"pgm")==0) || (strcmp(&name.c_str()[strlen(name.c_str())-3],"PGM")==0))
  {
    return write_pgm(name);
  }
  else
  {
    //OTHER FORMATS
    #ifdef __AMI_OPENCV__
      return write_with_opencv(name);
    #else
      #ifdef __AMI_IMAGE__
          int pos=name.find_last_of('.');
          if(pos == (int)std::string::npos) return -1;

          std::string extension=name.substr(pos+1);
          if( (extension == std::string("bmp")) || (extension == std::string("BMP")))
          {
            return write_with_ami_bmp(name);
          }
          printf("WRITE::Unrecognized image format\n");
          return -1;
      #else
          int pos=name.find_last_of('.');
          if(pos == (int)std::string::npos) return -1;

          std::string extension=name.substr(pos+1);
          if( (extension == std::string("PNG")) || (extension == std::string("png")))
          {
            return write_with_pnglib(name);
          }
          if( (extension == std::string("JPG")) || (extension == std::string("jpg")) ||
              (extension == std::string("JPEG")) || (extension == std::string("jpeg")))
          {
            return write_with_jpeglib(name);
          }
          if( (extension == std::string("TIF")) || (extension == std::string("tif")) ||
              (extension == std::string("TIFF")) || (extension == std::string("tiff")))
          {
            return write_with_tifflib(name);
          }
          if( (extension == std::string("bmp")) || (extension == std::string("BMP")))
          {
            return write_with_ami_bmp(name);
          }
          printf("WRITE::Unrecognized image format\n");
          return -1;
      #endif
    #endif
  }
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> T& image<T>::operator()(const int &x,const int &y,const int &channel)
 * \brief operator () to acces image value
 * \author Javier Martin
*/
template <class T>
T& image<T>::operator()(const int &x,const int &y,const int &channel){
  #ifdef IMAGE_DEBUG
    if((x+y*width_+channel*width_*height_)>=(image_.size())){
      printf("image<T>: bounds error vector access\n");
      printf("image size()=%d index to accces the vector =%d\n",
                                        image_.size(),x+y*width_+channel*width_*height_);
      int j; scanf("%d",&j);
      exit(0);
    }
  #endif
  return image_[x+y*width_+channel*width_*height_];
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> const T& image<T>::operator[](const int &i) const
 * \brief operator [] to acces image value
 * \author Luis Alvarez
*/
template <class T>
const T& image<T>::operator[](const int &i) const
{
  #ifdef IMAGE_DEBUG
    if(i>=(int)(image_.size())){
      printf("image<T>: bounds error vector access\n");
      printf("image size()=%d index to accces the vector =%d\n",(int)image_.size(),i);
      int j; scanf("%d",&j);
      exit(0);
    }
  #endif
  return image_[i];
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> T& image<T>::operator[](const int &i)
 * \brief operator [] to acces image value
 * \author Luis Alvarez
*/
template <class T>
T& image<T>::operator[](const int &i)
{
  #ifdef IMAGE_DEBUG
    if(i>=(int)(image_.size())){
      printf("image<T>: bounds error vector access\n");
      printf("image size()=%d index to accces the vector =%d\n",(int)image_.size(),i);
      int j; scanf("%d",&j);
      exit(0);
    }
  #endif
  return  (image_.at(i));
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T>::image(int width,int height)
 * \brief image constructor taking memory
 * \author Luis Alvarez
*/
template <class T>
image<T>::image(int width,int height)
{
  width_=width ;
  height_=height;
  nChannels_=1;
  image_.resize(width*height);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T>::image(const int width,const int height, const T &a, const T &b, const T &c)
 * \brief image constructor taking memory and initialiting the channels: R with a value, G with b value and B with c value
 * \author Carlos Falcon
*/
template <class T>
image<T>::image(const int width,const int height, const T &a, const T &b, const T &c)
{
  width_=width ;
  height_=height;
  nChannels_=3;
  int size=width_*height_;
  image_.resize(size*3);
  int i,j;
  #pragma omp parallel for shared(size) private(i)
  for(i=0;i<size;i++)
    image_.at(i)=a;
  #pragma omp parallel for shared(size) private(i)
  for(i=size;i<2*size;i++)
    image_.at(i)=b;
  #pragma omp parallel for shared(size) private(i)
  for(i=2*size;i<3*size;i++)
    image_.at(i)=c;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T>::image(const int width,const int height, const int nChannels, const T &a)
 * \brief image constructor taking memory and initialiting with a value
 * \author Luis Alvarez
*/
template <class T>
image<T>::image(const int width,const int height, const int nChannels, const T &a)
{
  width_=width ;
  height_=height;
  nChannels_=nChannels;
  int size=width_*height_*nChannels_;
  image_.resize(size);
  int i;
  #pragma omp parallel for shared(size) private(i)
  for(i=0;i<size;i++) image_.at(i)=a;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn template <class T> image<T>::image(const int width,const int height, const int nChannels)
 * \brief image constructor taking memory and initialiting with a value
 * \author Luis Alvarez
*/
template <class T>
image<T>::image(const int width,const int height, const int nChannels)
{
  width_=width;
  height_=height;
  nChannels_=nChannels;
  int size=width_*height_*nChannels_;
  image_.resize(size);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T>::~image()
 * \brief basic constructor
 * \author Luis Alvarez
*/
template <class T>
image<T>::image() : width_(0), height_(0), nChannels_(0) {}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn image<T>::image()
 * \brief basic destructor
 * \author Luis Alvarez
*/
template <class T>
image<T>::~image(){}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
 * \fn  template <class T> image<T>::image(std::string name)
 * \brief Constructor from an image file
 * \param [in] name Input file name.
 * \author Pedro Henriquez
*/
template <class T>
image<T>::image(std::string name /** INPUT FILE NAME */ )
{
  read(name);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn void image<T>::clear()
* \brief Function to put the image to 0 and clear its vectors.
* \author Pedro Henríquez
*/
template <class T>
void image<T>::clear()
{
  width_=0;
  height_=0;
  nChannels_=0;
  roi_clear();
  image_.clear();
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn void image<T>::set_nchannels(int nchannel1)
* \brief Change number of channels
* \param [in] nchannel1 New number of channels.
* \author Pedro Henríquez
*/
template <class T>
void image<T>::set_nchannels(int nchannel1)
{
  nChannels_=nchannel1;
  image_.resize(nChannels_*width_*height_);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn void image<T>::set_size(int width,int height)
* \brief Change the image dimensions
* \param [in] width,height New image dimensions.
* \author Pedro Henríquez
*/
template <class T>
void image<T>::set_size(int width,int height)
{
  width_=width;
  height_=height;
  image_.resize(nChannels_*width_*height_);
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn image<T>::imageMirrored()
* \brief mirrored the image
* \author Javier Martín Abasolo
*/
template <class T>
void image<T>::imageMirrored()
{
  //#pragma omp parallel for
  for(int k=0; k<nChannels_;k++){                                                                     // all of channels
    for(int i=0; i< height_;i++){                                                                   // reading for columns
      for(int j=0; j<(width_/2);j++){                                                             // reading for rows
        T aux = image_[(k*height_*width_)+i*width_+j];                                          //
        image_[(k*height_*width_)+i*width_+j] = image_[(k*height_*width_)+i*width_+width_-j-1]; // swapping mirrored values in the row
        image_[(k*height_*width_)+i*width_+width_-j-1] = aux;                                   //
      }
    }
  }
  return;
}

#ifdef __AMI_OPENCV__
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn template <class T> image<T> image<T>::resizeImage(int width, int height)
* \brief resize the image
* \author Javier Martín Abasolo
*/
template <class T>
image<T> image<T>::resizeImage(int width, int height){
  int tipo;
  if (typeid( image_[0] ) == typeid( unsigned char )) tipo = IPL_DEPTH_8U;
  if (typeid( image_[0] ) == typeid( char ))         tipo = IPL_DEPTH_8S;
  if (typeid( image_[0] ) == typeid( unsigned short )) tipo= IPL_DEPTH_16U;
  if (typeid( image_[0] ) == typeid( short )) tipo= IPL_DEPTH_16U;
  if (typeid( image_[0] ) == typeid( int )) tipo = IPL_DEPTH_32S;
  if (typeid( image_[0] ) == typeid( float )) tipo = IPL_DEPTH_32F;
  if (typeid( image_[0] ) == typeid( double )) tipo = IPL_DEPTH_64F;

  //CREATE OPENCV IPLIMAGE
  IplImage* image2 = cvCreateImage(cvSize(width_,height_), tipo, nChannels_);

  if (image2 == NULL) return image<T>();

  long cont=0;
  long size_=width_*height_;

  if(nChannels_==3) //CONVERT RGB TO BGR
  {
    for(int i=height_-1;i>=0;i--)
    {
      for(int j=0;(int) j<width_;j++)
      {
        long m=i*width_+j;
        image2->imageData[cont++]= image_[m+2*size_]; //B
        image2->imageData[cont++]= image_[m+size_]; //G
        image2->imageData[cont++]= image_[m]; //R
      }
      int incr = 0;
      if((3*width_)%4 != 0) incr = (4-(3*width_)%4);
      for(int k=0; (int) k<incr; k++) image2->imageData[cont++]=0;
    }
  }
  else
  {
    for(int i=height_-1;i>=0;i--)
    {
      for(int j=0;(int)j<width_;j++)
      {
        long m=i*width_+j;
        for(int nc=0;(int) nc<nChannels_;nc++)
        {
          image2->imageData[cont++]= image_[m+nc*size_];
        }
      }
      int incr = 0;
      if((nChannels_*width_)%4 != 0) incr = (4-(nChannels_*width_)%4);
      for(int k=0; (int) k<incr; k++) image2->imageData[cont++]=0;
    }
  }
  int px = width;
  int py = height;
  cont = 0;
  long size = ((width))*(height);
  IplImage *resized = cvCreateImage(cvSize(px,py),tipo, nChannels_);
  cvResize(image2,resized,CV_INTER_CUBIC); // CV_INTER_LINEAR - CV_INTER_CUBIC///////////////////////////////////

  image<T> imResized(resized->width, resized->height, resized->nChannels);


  if(nChannels_==3)//CONVERT BGR TO RGB
  {
    for(long i=height-1;i>=0;i--)
    {
      long m=i*width;
      for(long j=0;j<(long)width;j++,m++)
      {
        imResized[m+2*size]=resized->imageData[cont++]; //B
        imResized[m+size]=resized->imageData[cont++];//G
        imResized[m]=resized->imageData[cont++];//R
      }
      cont+=width_%4;
    }
  }
  else
  {
    for(long i=height-1;i>=0;i--)
    {
      long m=i*width;
      for(long j=0;j<(long)width;j++,m++)
      {
        for(long n=0;n<(long) nChannels_;n++)
        {
          imResized[m+n*size]=resized->imageData[cont++];
        }
      }
      cont+=width_%4;
    }
  }
  cvReleaseImage (&resized);
  cvReleaseImage (&image2);
  return imResized;
}
#endif

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn int image<T>::read_pgm (std::string name)
* \brief Read a pgm image
* \author Javier Martín
*/
template <class T>
int image<T>::read_pgm (std::string name)
{
    // PGM
    char tmp[100],*p,c;
    int wpp, Width, Height, Depth=0;
    FILE *fp;
    fp=fopen(name.c_str(),"rb");
    if(fp == NULL){
        printf("Image can't be loaded %s\n", name.c_str());
        return -1;
    }
    fread(tmp,2,1,fp);
    /* ---- Check sign ---- */
    if (!memcmp(tmp,"P6",strlen("P6")))      wpp=3; //PPM no used
    else if (!memcmp(tmp,"P5",strlen("P5"))) wpp=1;
    else{
        printf("Image can't be loaded %s\n", name.c_str());
        fclose(fp);
        return -1;
    }
    /* ---- Get comment ---- */
    while ((c=fgetc(fp)) && (c!=0x0a))       /* goto new line */
        if (c==EOF){
            printf("Image can't be loaded %s\n", name.c_str());
            fclose(fp);
            return -1;
        }
    while ((c=fgetc(fp)) && (c==' ')){};       /* skip spaces */
    fseek(fp,-1,SEEK_CUR);
    while ((c=fgetc(fp)) && (c=='#'))        /* detect comment */
    {
        while ((c=fgetc(fp)) && (c!=0x0a))
        {
            if (c==EOF){
                printf("Image can't be loaded %s\n", name.c_str());
                fclose(fp);
                return -1;
            }
        }
    }
    fseek(fp,-1,SEEK_CUR);
    /* ---- get size ---- */
    p=tmp;
    while ((c=fgetc(fp)) && (c!=0x0a))
    {
        if (c==EOF){
            printf("Image can't be loaded %s\n", name.c_str());
            fclose(fp);
            return -1;
        }
        *p=c; p++;
        if (p-tmp >= 99){
            printf("Image can't be loaded %s\n", name.c_str());
            fclose(fp);
            return -1;
        }
    }
    *p=0;
    if (!(p=strtok(tmp," "))){
        printf("Image can't be loaded %s\n", name.c_str());
        fclose(fp);
        return -1;
    }
    Width=atoi(p);
    if (!(p=strtok(NULL," "))){
        printf("Image can't be loaded %s\n", name.c_str());
        fclose(fp);
        return -1;
    }
    Height=atoi(p);
    /* get max. value; goto end header */
    p=tmp;
    while ((c=fgetc(fp)) && (c!=0x0a))
    {
        if (c==EOF){
            printf("Image can't be loaded %s\n", name.c_str());
            fclose(fp);
            return -1;
        }
        *p=c; p++;
        if (p-tmp >= 99){
            printf("Image can't be loaded %s\n", name.c_str());
            fclose(fp);
            return -1;
        }
    }
    *p=0;
    if (!(p=strtok(tmp," "))){
        printf("Image can't be loaded %s\n", name.c_str());
        fclose(fp);
        return -1;
    }
    Depth=atoi(p);

    width_ = Width;
    height_= Height;
    nChannels_=1; // default pgm
    int size = width_*height_;

    if(Depth > 255){ // 2 bytes pp
        char *data = (char*) malloc(2*size);
        for(int i = 2*size-1;i>=0;i--){
            data[i] = fgetc(fp);
        }
        unsigned short int *pdata = (unsigned short int*) data;
        image_.resize(size);
        for (int i=0;i<size;i++){
            image_[i] = (T) pdata[i];
        }
    }else{  // 1 bytes pp
        char *data = (char*) malloc(size);
        for(int i = size-1;i>=0;i--){
            data[i] = fgetc(fp);
        }
        unsigned char *pdata = (unsigned char*) data;
        image_.resize(size);
        for (int i=0;i<size;i++){
            image_[i] = (T) pdata[i];
        }
        fclose(fp);
    }

    return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////
#ifdef __AMI_OPENCV__
/**
* \fn int image<T>::read_with_opencv (std::string name)
* \brief Read an image
* \author Pedro Henríquez
*/
template <class T>
int image<T>::read_with_opencv (std::string name)
{
    #ifdef __AMIDEBUG__
        printf("read_with_opencv");
    #endif
    //OTHER FORMATS
    IplImage *view = 0;
    //WHEN WE USE THIS FUNCTION img=cvLoadImage(fileName,flag);
    //THE FLAG VALUE DETERMINES THE NUMBER OF CHANNELS WE READ.
    //flag <0 THE LOADED IMAGE IS LOADED AS IS (WITH NUMBER OF CHANNELS IN THE FILE)
    view = cvLoadImage(name.c_str(),-1);
    if (!view)
    {
        //IF WE CAN'T LOAD THE IMAGE, RETURN DIMENSIONS 0
        printf("Image can't be loaded %s\n", name.c_str());
        int x;
        printf("\nPres any key to continue...");
        scanf("%d",&x);
        width_ = 0;
        height_= 0;
        nChannels_=0;
        return 0;
    }
    width_ = view->width;
    height_= view->height;
    nChannels_=view->nChannels;
    long size_ = ((width_))*(height_);
    long cont=0;
    image_.resize(nChannels_*size_);
    if(nChannels_==3)//CONVERT BGR TO RGB
    {
     for(long i=height_-1;i>=0;i--)
     {
        long m=i*width_;
        for(long j=0;j<(long)width_;j++,m++)
        {
            image_[m+2*size_]=view->imageData[cont++]; //B
            image_[m+size_]=view->imageData[cont++];//G
            image_[m]=view->imageData[cont++];//R
        }
        cont+=width_%4;
     }
    }
    else
    {
        for(long i=height_-1;i>=0;i--)
        {
          long m=i*width_;
          for(long j=0;j<(long)width_;j++,m++)
          {
            for(long n=0;n<(long) nChannels_;n++)
            {
                image_[m+n*size_]=view->imageData[cont++];
            }
          }
          cont+=width_%4;
        }
    }
    cvReleaseImage (&view);

    return 1;
}
#else

#ifdef __AMI_IMAGE__

#else
//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn int image<T>::read_with_jpeglib (std::string name)
* \brief Read an image with jpeglib
* \author Pedro Henríquez
*/
template <class T>
int image<T>::read_with_jpeglib (std::string name)
{
  struct jpeg_decompress_struct cinfo;
	struct jpeg_error_mgr jerr;
	cinfo.err = jpeg_std_error(&jerr);
	jpeg_create_decompress(&cinfo);

	FILE * infile;
	if ((infile = fopen(name.c_str(), "rb")) == NULL)
	{
	  fprintf(stderr, "can't open %s\n", name.c_str());
	  return(-1);
	}
	jpeg_stdio_src(&cinfo, infile);
  jpeg_read_header(&cinfo, TRUE);
  jpeg_start_decompress(&cinfo);

  // Make a one-row-high sample array that will go away when done with image
	int rowstride = cinfo.output_width * cinfo.output_components;
	JSAMPARRAY buffer = (*cinfo.mem->alloc_sarray)((j_common_ptr) &cinfo, JPOOL_IMAGE, rowstride, 1);

	// Get memory for local image
	char *data = new char[cinfo.output_height * rowstride];

	// Read all scanlines
	int y = 0;
 	while( cinfo.output_scanline < cinfo.output_height )
 	{
	  jpeg_read_scanlines(&cinfo, buffer, 1);
		// insert the row into the global image
		memcpy(&data[y*rowstride], buffer[0], rowstride);
		y++;
	}

	// Fill in ImageInfo struct
	width_ = cinfo.output_width;
	height_ = cinfo.output_height;
	nChannels_ = cinfo.output_components;
	long size_ = ((width_))*(height_);
	image_.resize(nChannels_*size_);

  long cont=0;
  if(nChannels_==3)//CONVERT RGB BY PIXEL IN RGB CHANNELS
  {
    for(long i=0;i<height_;i++)
    {
      long m=(height_-1-i)*width_;
      for(long j=0;j<(long)width_;j++,m++)
      {
        image_[m]=        (T)data[cont++];//R
        image_[m+size_]=  (T)data[cont++];//G
        image_[m+2*size_]=(T)data[cont++]; //B
      }
    }
  }
  else
  {
    for(long i=0;i<height_;i++)
    {
      long m=(height_-1-i)*width_;
      for(long j=0;j<(long)width_;j++,m++)
      {
        for(long n=0;n<(long) nChannels_;n++)
        {
          image_[m+n*size_]=(T)data[cont++];
        }
      }
    }
  }

  // Release JPEG decompression object
	jpeg_destroy_decompress(&cinfo);

	// Close the input file
	fclose(infile);

  delete []data;
	return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn int image<T>::read_with_pnglib (std::string name)
* \brief Read an image with pnglib
* \author Pedro Henríquez
*/
template <class T>
int image<T>::read_with_pnglib (std::string name)
{
  char header[8];	// 8 IS THE MAXIMUM SIZE THAT CAN BE CHECKED

	// OPEN FILE AND TEST FOR IT BEING A PNG
	FILE *imf = fopen(name.c_str(), "rb");
	if (!imf)
	{
    #ifdef __AMIDEBUG__
      printf("File %s could not be opened for reading\n",name.c_str()); fflush(stdout);
    #endif
    return -1;
	}
  fread(header, 1, 8, imf);
  //INITIALIZE STUFF
  png_structp png_ptr;
  png_infop info_ptr;
	png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

	if (!png_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_read_struct failed\n");
		#endif
		return -2;
	}

	info_ptr = png_create_info_struct(png_ptr);
	if (!info_ptr)
	{
    #ifdef __AMIDEBUG__
      printf("png_create_info_struct failed\n");
    #endif
    return -3;
	}

	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during init_io\n");
    #endif
    return -4;
	}
	png_init_io(png_ptr, imf);
	png_set_sig_bytes(png_ptr, 8);

	png_read_info(png_ptr, info_ptr);

	width_ = info_ptr->width;
	height_ = info_ptr->height;
	nChannels_ = png_get_channels(png_ptr, info_ptr);
	image_.resize(width_*height_*nChannels_);

	//int number_of_passes =
	png_set_interlace_handling(png_ptr);
	png_read_update_info(png_ptr, info_ptr);

	// READ FILE
	if (setjmp(png_jmpbuf(png_ptr)))
	{
    #ifdef __AMIDEBUG__
      printf("Error during read_image\n");
    #endif
    return -5;
	}

	switch (info_ptr->bit_depth)
	{
	  case 8:
	  {
      #ifdef __AMIDEBUG__
        printf("\nread_with_pnglib 8: Entro\n");
      #endif
      png_bytep * row_pointers;
      row_pointers = (png_bytep*) malloc(sizeof(png_bytep) * height_);
      for (int y=0; y<height_; y++)
        row_pointers[y] = (png_byte*) malloc(info_ptr->rowbytes);

      png_read_image(png_ptr, row_pointers);

      fclose(imf);
      long cont=0;
      long size_=width_*height_;
      if(nChannels_==3)//CONVERT RGB BY PIXEL IN RGB CHANNELS
      {
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 8: 3 canales\n");
        #endif
        for(long i=height_-1;i>=0;i--)
        {
          long m=(height_-1-i)*width_;
          cont = 0;
          for(long j=0;j<(long)width_;j++,m++)
          {
            T value;
            value = (T)row_pointers[i][cont++];//R
            image_[m]=value;
            value = (T)row_pointers[i][cont++];//G
            image_[m+size_]=value;//G
            value = (T)row_pointers[i][cont++];//B
            image_[m+2*size_]=value; //B
          }
        }
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 8: FIN 3 canales\n");
        #endif
      }
      else
      {
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 8: %d canales\n", nChannels_);
        #endif
        for(long i=height_-1;i>=0;i--)
        {
          long m=(height_-1-i)*width_;
          cont = 0;
          for(long j=0;j<(long)width_;j++,m++)
          {
            for(long n=0;n<(long) nChannels_;n++)
            {
              image_[m+n*size_]=row_pointers[i][cont++];
            }
          }
        }
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 8: FIN %d canales\n", nChannels_);
        #endif
      }
      //CLEANUP HEAP ALLOCATION
      for (int y=0; y<height_; y++)free(row_pointers[y]);
      free(row_pointers);
      break;
	  }
	  case 16:
	  {
      png_uint_16 **row_pointers;
      row_pointers = (png_uint_16**) malloc(sizeof(png_uint_16*) * height_);
      for (int y=0; y<height_; y++)
        row_pointers[y] = (png_uint_16*) malloc(sizeof(png_uint_16) * info_ptr->rowbytes);

      png_read_image(png_ptr, (png_bytep *)row_pointers);

      fclose(imf);
      long cont=0;
      long size_=width_*height_;
      T value_d, value_m;
      unsigned char *data_d, *data_m;
      data_m = (unsigned char*)&value_m;
      data_d = (unsigned char*)&value_d;
      if(nChannels_==3)//CONVERT RGB BY PIXEL IN RGB CHANNELS
      {
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 16: 3 canales\n");
        #endif
        for(long i=height_-1;i>=0;i--)
        {
          long m=(height_-1-i)*width_;
          cont = 0;
          for(long j=0;j<(long)width_;j++,m++)
          {
//            printf("%x ", row_pointers[i][cont]);
            value_d = (T)row_pointers[i][cont++];//R
            data_m[0] = data_d[1];
            data_m[1] = data_d[0];
            image_[m]=value_m;
//            printf("%x ", row_pointers[i][cont]);
            value_d = (T)row_pointers[i][cont++];//G
            data_m[0] = data_d[1];
            data_m[1] = data_d[0];
            image_[m+size_]=value_m;//G
//            printf("%x\n", row_pointers[i][cont]);
            value_d = (T)row_pointers[i][cont++];//B
            data_m[0] = data_d[1];
            data_m[1] = data_d[0];
            image_[m+2*size_]=value_m; //B
          }
        }
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 16: end 3 canales\n");
        #endif
      }
      else
      {
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 16: %d canales\n", nChannels_);
        #endif
        for(long i=height_-1;i>=0;i--)
        {
          long m=(height_-1-i)*width_;
          cont = 0;
          for(long j=0;j<(long)width_;j++,m++)
          {
            for(long n=0;n<(long) nChannels_;n++)
            {
              image_[m+n*size_]=row_pointers[i][cont++];
            }
          }
        }
        #ifdef __AMIDEBUG__
            printf("\nread_with_pnglib 16: %d canales\n", nChannels_);
        #endif
      }
      //CLEANUP HEAP ALLOCATION
      for (int y=0; y<height_; y++)free(row_pointers[y]);
      free(row_pointers);
      break;
	  }
	  case 1:
	  case 2:
	  case 4:
	  default:
	    return -6;
	}
  return 0;
}

//////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn int image<T>::read_with_tifflib (std::string name)
* \brief Read an image with tifflib
* \author Pedro Henríquez
*/
template <class T>
int image<T>::read_with_tifflib (std::string name)
{
  char *name_c;
  name_c = (char *)malloc(sizeof(char)*(name.length()+1));
  snprintf(name_c, name.length()+1, "%s", name.c_str());
  TIFF* tif = TIFFOpen(name_c, "r");
  if (tif)
  {
    uint32 w, h;
    size_t npixels;
    uint32* raster;

    TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &w);
    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &h);

    height_ = h;
    width_  = w;
    nChannels_ = 3;
    long size_ = width_ * height_;

    image_.resize(w*h*nChannels_);

    npixels = w * h;
    raster = (uint32*) _TIFFmalloc(npixels * sizeof (uint32));
    uint32 value;
    unsigned char *data;
    if (raster != NULL)
    {
      long cont = 0;
      if (TIFFReadRGBAImage(tif, w, h, raster, 0))
      {
         for(long i=height_-1;i>=0;i--)
         {
            long m=(height_-1-i)*width_;
            for(long j=0;j<(long)width_;j++,m++)
            {
                value = raster[cont++];
                data = (unsigned char *)&value;
                image_[m]         =(T)data[0];//R
                image_[m+size_]   =(T)data[1];//G
                image_[m+2*size_] =(T)data[2];//B
            }
         }
      }
      _TIFFfree(raster);
    }
    TIFFClose(tif);
    return 0;
  }
  else
    return -1;
}

#endif //AMITIFF

/**
* \fn int image<T>::read_with_ami_bmp (std::string name)
* \brief Read an image with ami_bmp
* \author Jesus Rodriguez
*/
template <class T>
int image<T>::read_with_ami_bmp (std::string name)
{
  // BUILD NAME
  char name_c[200];
  snprintf(name_c, 199, "%s", name.c_str());
  // READ IMAGE
  int width, height;
  unsigned char *red, *green, *blue;
  int output_value = ami_read_bmp(name_c, &red, &green, &blue, &width, &height);
  if (output_value != 0) return output_value;
  // ACTUALIZE INTERNAL INFORMATION OF IMAGE
  width_ = width;
  height_= height;
  nChannels_ = 3;
  // INITIALIZE INDIXES
  int size_ = width_ * height_;
  int mr = 0;
  int mg = size_;
  int mb = 2*size_;
  // RESIZE IMAGE
  image_.resize(size_*nChannels_);
  // REFILL IMAGE WITH DATA
  for(long i=height_-1;i>-1;i--)
  {
    int m = (height_-1-i)*width_;
    for(long j=0;j<width_;j++,mr++,mg++,mb++,m++)
    {
        image_[mr] =(T)red[m];//R
        image_[mg] =(T)green[m];//G
        image_[mb] =(T)blue[m];//B
    }
  }
  // FREE MEMORY
  free(red);  free(green);  free(blue);
  return 0;
}

#endif //OPENCV
////////////////////////////////////////////////////////////////////////////////////////////////
/**
* \fn int image<T>::read (std::string name)
* \brief Read an image selecting the library depend on the image format, returns 0 when it can't loads the image.
* \author Pedro Henríquez
*/
template <class T>
int image<T>::read (std::string name)
{

    #ifdef __AMIDEBUG__
        printf("read generico\n");
    #endif
    char path[512];
    strcpy(path,name.c_str());
    if ((strcmp(&name.c_str()[strlen(name.c_str())-3],"pgm")==0) || (strcmp(&name.c_str()[strlen(name.c_str())-3],"PGM")==0))
    {
        return read_pgm(name);
    }
    else
    {
        #ifdef __AMI_OPENCV__
            return read_with_opencv(name);
        #else
            #ifdef __AMI_IMAGE__
                FILE *mf;
                mf=fopen(name.c_str(),"rb");
                if (!mf) return -1;
                int mn;
                fread (&mn, 2, 2,mf);
                #ifdef __AMIDEBUG__
                    printf("%x: %s: ",mn, name.c_str());
                #endif
                fclose(mf);
                switch(mn)
                {
                    case 0x2a004d4d://big endian
                    case 0x002a4949://little endian
                        #ifdef __AMIDEBUG__
                            printf("READ::Unreadable tif file\n");
                        #endif
                        //return read_with_ami_tif(name);
                        width_=0;
                        height_=0;
                        return -1;
                        break;

                    case 0xe0ffd8ff:
                    case 0x474e5089:
                         #ifdef __AMIDEBUG__
                              printf("READ::Unreadable image format\n");
                          #endif
                          width_=0;
                          height_=0;
                          return -1;
                          break;

                    default:
                        if ((mn&0x0000ffff) != 0x00004d42)
                        {
                          #ifdef __AMIDEBUG__
                              printf("READ::Unrecognized image format\n");
                          #endif
                          return -1;
                        }
                        #ifdef __AMIDEBUG__
                            printf("bmp file\n");
                        #endif
                        return read_with_ami_bmp(name);
                        break;
                }
            #else
                FILE *mf;
                mf=fopen(name.c_str(),"rb");
                if (!mf) return -1;
                int mn;
                fread (&mn, 2, 2,mf);
                #ifdef __AMIDEBUG__
                    printf("%x: %s: ",mn, name.c_str());
                #endif
                fclose(mf);
                switch(mn)
                {
                    case 0xe0ffd8ff:
                        #ifdef __AMIDEBUG__
                            printf("jpg file\n");
                        #endif
                        return read_with_jpeglib(name);
                    break;

                    case 0x474e5089:
                        #ifdef __AMIDEBUG__
                            printf("png file\n");
                        #endif
                        return read_with_pnglib(name);
                    break;

                    case 0x2a004d4d://big endian
                    case 0x002a4949://little endian
                        #ifdef __AMIDEBUG__
                            printf("tif file\n");
                        #endif
                        return read_with_tifflib(name);
                        break;
                    default:
                        if ((mn&0x0000ffff) != 0x00004d42)
                        {
                          #ifdef __AMIDEBUG__
                              printf("READ::Unrecognized image format\n");
                          #endif
                          return -1;
                        }
                        #ifdef __AMIDEBUG__
                            printf("bmp file\n");
                        #endif
                        return read_with_ami_bmp(name);
                        break;
                }
            #endif
        #endif
    }
}

}

#endif
