/*
    Copyright (c) 2010, Luis Alvarez <lalvarez@dis.ulpgc.es>
	All rights reserved.
	
	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_filters.h \link filters.h \endlink
 * \brief Code of  \ref code_filters.h "filters.h"
 * \includelineno filters.h
 */

/**
 * \file filters.h
 * \brief filters class definition
 * \author Luis Alvarez \n \n
 * \brief \b Code: \n
 * - \ref code_filters.h "filters.h"
*/


#ifndef filters_H
#define filters_H

#ifndef NOMINMAX
#define NOMINMAX
#endif


#include <iostream>
#include <math.h>
#include <string.h>
#include <vector>
#include "utilities.h"
#include "image.h"

#define filters_DEBUG

using namespace std;

namespace ami
{


class filters
{
  public :

  template <class T,class U>
  void grad(const ami::image<T> &img,ami::image<U> &grad_x, ami::image<U> &grad_y, const bool NeigborhoodType=0);

  template <class T,class U>
  void gauss_conv(ami::image<T> &img,ami::image<U> &img_conv,double sigma_x,double sigma_y,double precision=1);

  template <class T,class U>
  void line_opening(ami::image<T> &img,ami::image<U> &img_out, const int radius_erosion, const int radius_dilation);

  template <class T,class U>
  void line_closing(ami::image<T> &img,ami::image<U> &img_out, const int radius_erosion, const int radius_dilation);
};

}


/**
 * \fn template <class T,class U> void ami::filters<T,U>::grad(const ami::image<T> &img,ami::image<U> &grad_x, ami::image<U> &grad_y, const bool NeigborhoodType)
 * \brief Gradient Computation
 * \author Luis Alvarez
*/
template <class T,class U>
void ami::filters::grad(
  const ami::image<T> &img /** INPUT IMAGE */,
  ami::image<U> &grad_x /** OUTPUT x-GRADIENT IMAGE */,
  ami::image<U> &grad_y /** OUTPUT y-GRADIENT IMAGE */,
  const bool NeigborhoodType /** =0 means 5 size neighborhood, =9  means 9 size neighborhood*/)
{
  if(img.get_roi().size()==6){
    ami::image<T> img2;
    img.get_roi_image(img2);
    grad(img2,grad_x,grad_y,NeigborhoodType);
    return;
  }

//  if(img.height()<100){
//    grad_no_omp(img,grad_x,grad_y,NeigborhoodType);
//    return;
//  }

  int width=img.width();
  int height=img.height();
  int nChannels=img.nChannels();
  int size_image=width*height;
  int size_image_width=size_image-width;
  int total_image_size=size_image*nChannels;
  //int height_1=height-1;
  int width_1=width-1;

  if(width!=grad_x.width() ||  height!=grad_x.height() || nChannels!=grad_x.nChannels()){
    grad_x=image<U>(width,height,nChannels);
  }

  if(width!=grad_y.width() ||  height!=grad_y.height() || nChannels!=grad_y.nChannels()){
    grad_y=image<U>(width,height,nChannels);
  }

  if(NeigborhoodType==0){

    ami::utilities u;
    vector < vector <int> > b=u.boundary_neighborhood_5n(width,height);
    int m,c,y,k,k_end,b_size=b.size();//l,

    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,b,b_size) \
    private(c,y,m,k,k_end)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          grad_y[k]=(U) img[k+width]- (U) img[k-width];
          grad_x[k]=(U) img[k-1]- (U) img[k+1];
        }
      }
      #pragma omp for nowait
      for(int k=0;k<b_size;k++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
        grad_y[b[k][0]+c]=(U) img[b[k][2]+c]- (U) img[b[k][1]+c];
        grad_x[b[k][0]+c]=(U) img[b[k][4]+c]- (U) img[b[k][3]+c];
      }
    }
  }
  else{
    ami::utilities u;
    vector < vector <int> > b=u.boundary_neighborhood_9n(width,height);
    double coef1,coef2,c1,d1;
    coef1=sqrt((double) 2.);
    coef2=0.25*(2.-coef1);
    coef1=0.5*(coef1-1);
    int m,c,y,k,l,k_end,b_size=b.size();
    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,b,b_size,coef1,coef2) \
    private(c,y,m,k,k_end,c1,d1)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          c1=img[k+width+1]-img[k-width-1];
          d1=img[k-width+1]-img[k+width-1];
          grad_y[k]=(U)(coef1*((U) img[k+width]- (U) img[k-width])+coef2*(c1-d1));
          grad_x[k]=(U)(-(coef1*((U) img[k+1]- (U) img[k-1])+coef2*(c1+d1)));
        }
      }
      #pragma omp for nowait
      for(l=0;l<b_size;l++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
        c1=img[c+b[l][6]]-img[c+b[l][7]];
        d1=img[c+b[l][5]]-img[c+b[l][8]];
        grad_y[c+b[l][0]]=(U)(coef1*((U) img[c+b[l][2]]- (U) img[c+b[l][1]])+coef2*(c1-d1));
        grad_x[c+b[l][0]]=(U)(-(coef1*((U) img[c+b[l][3]]- (U) img[c+b[l][4]])+coef2*(c1+d1)));
      }
    }
  }
}


/**
 * \fn template <class T,class U> void ami::filters<T,U>::gauss_conv(ami::image<T> &img,ami::image<U> &img_conv,double sigma_x,double sigma_y,double precision)
 * \brief Gaussian convolution filter (using Alvarez-Mazorra implementation)
 * \author Luis Alvarez
*/
template <class T,class U>
void ami::filters::gauss_conv(
ami::image<T> &img /** INPUT IMAGE (IT CAN BE = OUTPUT IMAGE */,
ami::image<U> &img_conv /** OUTPUT IMAGE */ ,
double sigma_x /** STANDARD DEVIATION IN THE x VARIABLE */,
double sigma_y /** STANDARD DEVIATION IN THE y VARIABLE */,
double precision /** PRECISION TO COMPUTE THE GAUSSIAN CONVOLUTION BIGGER IS THE PRECISION, MORE ITERATIONS
TO COMPUTE THE CONVOLUTION */)
{

//  if(img.height()<100){
//    gauss_conv_no_omp(img,img_conv,sigma_x,sigma_y,precision);
//    return;
//  }
  if(img.get_roi().size()==6){
    img.get_roi_image(img_conv);
  }
  else{
    if(img.width()!=img_conv.width() ||  img.height()!=img_conv.height() || img.nChannels()!=img_conv.nChannels()){
      img_conv=image<U>(img.width(),img.height(),img.nChannels());
    }
    int k,k_end=img.width()*img.height()*img.nChannels();
    #pragma omp parallel for shared(k_end) private(k)
    for(k=0;k<k_end;k++) img_conv[k]= (U) img[k];
  }

  int width=img_conv.width();
  int height=img_conv.height();
  int nChannels=img_conv.nChannels();
  int image_size=width*height;
  int image_size_total=image_size*nChannels;

  if(precision<0) precision=0;

  int Nc_x=(int) (precision*sigma_x>1?precision*sigma_x:1);
  int Nc_y=(int) (precision*sigma_y>1?precision*sigma_y:1);
  int Nc_a=(Nc_x>Nc_y)?Nc_x:Nc_y;
  Nc_x=Nc_y=Nc_a;
  double t_y=sigma_x*sigma_x/(2*Nc_x);
  double t_x=sigma_y*sigma_y/(2*Nc_y);
  double l_x,v_x,l_y,v_y,l_x_1,l_y_1;

  // WE CHECK IS STANDARD DEVIATIONS ARE 0
  if(t_x==0 && t_y==0) return;

  // PARAMETERS X AND Y
  if(t_x>0){
     l_x=(float)(1.+2.*t_x-sqrt((double) 4*t_x+1))/(2*t_x);
     v_x=l_x/t_x;
     l_x_1=(1-l_x);
  }
  if(t_y>0){
     l_y=(float)(1.+2.*t_y-sqrt((double) 4*t_y+1))/(2*t_y);
     v_y=l_y/t_y;
     l_y_1=(1-l_y);
  }

  //MAIN LOOP
  int cont;
  int c,y,x,x_end,m,y_end;
  for(cont=0;cont<Nc_a;cont++){
    if(t_y>0){
      #pragma omp parallel shared(width,image_size_total,image_size,l_y,v_y,l_y_1) \
      private(c,x,y,m,x_end)
      for(c=0;c<image_size_total;c+=image_size){
        #pragma omp for nowait
        for(y=0;y<image_size;y+=width){
          m=c+y;
          vector <U> paso(width);
          paso[0]= img_conv[m]/l_y_1;
          x_end=m+width;
          int n=1;
          for(x=m+1;x<x_end;x++,n++){
            paso[n]= img_conv[x]+l_y*paso[n-1];
          }
          img_conv[x_end-1]=paso[width-1]/l_y_1;
          n=width-2;
          for(x=x_end-2;x>m;x--,n--){
            img_conv[x]= paso[n]+l_y*img_conv[x+1];
          }
          img_conv[x]= paso[n]+l_y*img_conv[x+1];
          for(x=m;x<x_end;x++){
           img_conv[x]*=v_y;
          }
        }
      }
    }
    if(t_y>0){
     #pragma omp parallel shared(width,height,image_size_total,image_size,l_x,v_x,l_x_1) \
     private(c,x,y,m,y_end)
      for(c=0;c<image_size_total;c+=image_size){
        #pragma omp for nowait
        for(x=0;x<width;x++){
          m=c+x;
          vector <U> paso(height);
          paso[0]= img_conv[m]/l_x_1;
          y_end=m+image_size;
          int n=1;
          for(y=m+width;y<y_end;y+=width,n++){
            paso[n]= img_conv[y]+l_x*paso[n-1];
          }
          img_conv[y_end-width]=paso[height-1]/l_x_1;
          n=height-2;
          for(y=y_end-2*width;y>m;y-=width,n--){
            img_conv[y]= paso[n]+l_x*img_conv[y+width];
          }
          img_conv[y]= paso[n]+l_x*img_conv[y+width];
          for(y=m;y<y_end;y+=width){
           img_conv[y]*=v_x;
          }
        }
      }
    }
  }
}


// MACROS FOR MORPHOLOGICAL FILTERS
#define ami_min3(a,b,c) (a<b?(a<c?a:c):(b<c?b:c))
#define ami_max3(a,b,c) (a>b?(a>c?a:c):(b>c?b:c))
#define ami_max9(a0,b0,c0,a1,b1,c1,a2,b2,c2) ami_max3(ami_max3(a0,b0,c0),ami_max3(a1,b1,c1),ami_max3(a2,b2,c2))
#define ami_min9(a0,b0,c0,a1,b1,c1,a2,b2,c2) ami_min3(ami_min3(a0,b0,c0),ami_min3(a1,b1,c1),ami_min3(a2,b2,c2))
#define ami_max5(a0,b0,c0,a1,b1) ami_max3(ami_max3(a0,b0,c0),a1,b1)
#define ami_min5(a0,b0,c0,a1,b1) ami_min3(ami_min3(a0,b0,c0),a1,b1)



/**
 * \fn template <class T,class U> void ami::filters::line_opening(ami::image<T> &img,ami::image<U> &img_out, const int radius_erosion, const int radius_dilation)
 * \brief line opening morphological operator (first max of mins)
 * \author Luis Alvarez
*/
template <class T,class U>
void ami::filters::line_opening(
  ami::image<T> &img /** INPUT IMAGE (IT CAN BE ALSO = OUTPUT IMAGE)*/,
  ami::image<U> &img_out /** OUTPUT IMAGE */,
  const int radius_erosion /** RADIUS OF EROSION (max of mins) */,
  const int radius_dilation /** RADIUS OF DILATION (min of maxs)*/)
{
  // WE CHECK IF WE HAVE SELECTED A SUBWINDOW
  if(img.get_roi().size()==6){
    ami::image<T> img2;
    img.get_roi_image(img2);
    line_opening(img2,img_out,radius_erosion,radius_dilation);
    return;
  }
//  if(img.height()<100){
//    line_opening_no_omp(img,img_out,radius_erosion,radius_dilation);
//    return;
//  }

  // WE DEFINE AUXILARY VARIABLES
  int width=img.width();
  int height=img.height();
  int nChannels=img.nChannels();
  int size_image=width*height;
  int size_image_width=size_image-width;
  int total_image_size=size_image*nChannels;
  int height_1=height-1;
  int width_1=width-1;

  // WE MADE A COPY OF THE ORIGINAL IMAGE
  ami::image<U> img_paso=ami::image<U>(width,height,nChannels);
  {
    int i,k,k_end=img.width()*img.height()*img.nChannels();
    #pragma omp parallel for shared(k_end,img_paso) private(k)
    for(k=0;k<k_end;k++) img_paso[k]= (U) img[k];
  }

  if(width!=img_out.width() || height!=img_out.height() || nChannels!=img_out.nChannels())
    img_out=ami::image<U>(width,height,nChannels);

  // WE BUILD BOUNDARY COORDINATES TO MANAGE IMAGE BOUNDARY
  ami::utilities u;
  vector < vector <int> > b=u.boundary_neighborhood_9n(width,height);
  int m,c,y,k,l,k_end,b_size=b.size();
  U temp;

  // WE GO TRHOUGH ITERATIONS
  for(int l=0;l<radius_erosion;++l){
    if(l>0) img_paso=img_out;
    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,img_paso,b,b_size) \
    private(y,m,k,k_end,temp,c)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          U paso=ami_min3(img_paso[k],img_paso[k+1],img_paso[k-1]);
          temp=ami_min3(img_paso[k],img_paso[k+width],img_paso[k-width]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[k],img_paso[k+width+1],img_paso[k-width-1]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[k],img_paso[k+width-1],img_paso[k-width+1]);
          if(temp>paso) paso=temp;
          if(paso!=img_out[k]) img_out[k]=paso;
        }
      }

      #pragma omp for nowait
      for(int k=0;k<b_size;k++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
          m=c+b[k][0];
          U paso=ami_min3(img_paso[m],img_paso[c+b[k][2]],img_paso[c+b[k][1]]);
          temp=ami_min3(img_paso[m],img_paso[c+b[k][3]],img_paso[c+b[k][4]]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[m],img_paso[c+b[k][5]],img_paso[c+b[k][8]]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[m],img_paso[c+b[k][6]],img_paso[c+b[k][7]]);
          if(temp>paso) paso=temp;
          if(img_out[m]!=paso) img_out[m]=paso;
      }
    }
  }

  for(int l=0;l<radius_dilation;++l){
    img_paso=img_out;
    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,img_paso,b,b_size) \
    private(c,y,m,k,k_end,temp)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          U paso=ami_max3(img_paso[k],img_paso[k+1],img_paso[k-1]);
          temp=ami_max3(img_paso[k],img_paso[k+width],img_paso[k-width]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[k],img_paso[k+width+1],img_paso[k-width-1]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[k],img_paso[k+width-1],img_paso[k-width+1]);
          if(temp<paso) paso=temp;
          if(paso!=img_out[k]) img_out[k]=paso;
        }
      }

      #pragma omp for nowait
      for(int k=0;k<b_size;k++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
          m=c+b[k][0];
          U paso=ami_max3(img_paso[m],img_paso[c+b[k][2]],img_paso[c+b[k][1]]);
          temp=ami_max3(img_paso[m],img_paso[c+b[k][3]],img_paso[c+b[k][4]]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[m],img_paso[c+b[k][5]],img_paso[c+b[k][8]]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[m],img_paso[c+b[k][6]],img_paso[c+b[k][7]]);
          if(temp<paso) paso=temp;
          if(paso!=img_out[m]) img_out[m]=paso;
      }
    }
  }
}

/**
 * \fn template <class T,class U> void ami::filters<T,U>::line_closing(ami::image<T> &img,ami::image<U> &img_out, const int radius_erosion, const int radius_dilation)
 * \brief line closing morphological operator (first min of maxs)
 * \author Luis Alvarez
*/
template <class T,class U>
void ami::filters::line_closing(
  ami::image<T> &img /** INPUT IMAGE (IT CAN BE ALSO = OUTPUT IMAGE)*/,
  ami::image<U> &img_out /** OUTPUT IMAGE */,
  const int radius_erosion /** RADIUS OF EROSION (max of mins) */,
  const int radius_dilation /** RADIUS OF DILATION (min of maxs)*/)
{
  // WE CHECK IF WE HAVE SELECTED A SUBWINDOW
  if(img.get_roi().size()==6){
    ami::image<T> img2;
    img.get_roi_image(img2);
    line_closing(img2,img_out,radius_erosion,radius_dilation);
    return;
  }
//  if(img.height()<100){
//    line_closing_no_omp(img,img_out,radius_erosion,radius_dilation);
//    return;
//  }

  // WE DEFINE AUXILARY VARIABLES
  int width=img.width();
  int height=img.height();
  int nChannels=img.nChannels();
  int size_image=width*height;
  int size_image_width=size_image-width;
  int total_image_size=size_image*nChannels;
  int height_1=height-1;
  int width_1=width-1;

  // WE MADE A COPY OF THE ORIGINAL IMAGE
  ami::image<U> img_paso=ami::image<U>(width,height,nChannels);
  {
    int i,k,k_end=img.width()*img.height()*img.nChannels();
    #pragma omp parallel for shared(k_end,img_paso) private(k)
    for(k=0;k<k_end;k++) img_paso[k]= (U) img[k];
  }

  if(width!=img_out.width() || height!=img_out.height() || nChannels!=img_out.nChannels())
    img_out=ami::image<U>(width,height,nChannels);

  // WE BUILD BOUNDARY COORDINATES TO MANAGE IMAGE BOUNDARY
  ami::utilities u;
  vector < vector <int> > b=u.boundary_neighborhood_9n(width,height);
  int m,c,y,k,l,k_end,b_size=b.size();
  U temp;

  // WE GO TRHOUGH ITERATIONS
  for(int l=0;l<radius_dilation;++l){
    if(l>0) img_paso=img_out;
    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,img_paso,b,b_size) \
    private(y,m,k,k_end,temp,c)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          U paso=ami_max3(img_paso[k],img_paso[k+1],img_paso[k-1]);
          temp=ami_max3(img_paso[k],img_paso[k+width],img_paso[k-width]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[k],img_paso[k+width+1],img_paso[k-width-1]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[k],img_paso[k+width-1],img_paso[k-width+1]);
          if(temp<paso) paso=temp;
          if(paso!=img_out[k]) img_out[k]=paso;
        }
      }

      #pragma omp for nowait
      for(int k=0;k<b_size;k++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
          m=c+b[k][0];
          U paso=ami_max3(img_paso[m],img_paso[c+b[k][2]],img_paso[c+b[k][1]]);
          temp=ami_max3(img_paso[m],img_paso[c+b[k][3]],img_paso[c+b[k][4]]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[m],img_paso[c+b[k][5]],img_paso[c+b[k][8]]);
          if(temp<paso) paso=temp;
          temp=ami_max3(img_paso[m],img_paso[c+b[k][6]],img_paso[c+b[k][7]]);
          if(temp<paso) paso=temp;
          if(paso!=img_out[m]) img_out[m]=paso;
      }
    }
  }

  for(int l=0;l<radius_erosion;++l){
    img_paso=img_out;
    #pragma omp parallel shared(width,width_1,total_image_size,size_image,size_image_width,img_paso,b,b_size) \
    private(c,y,m,k,k_end,temp)
    for(c=0;c<total_image_size;c+=size_image){
      #pragma omp for nowait
      for(y=width;y<size_image_width;y+=width){
        m=c+y;
        k_end=m+width_1;
        for(k=m+1; k<k_end; k++){
          U paso=ami_min3(img_paso[k],img_paso[k+1],img_paso[k-1]);
          temp=ami_min3(img_paso[k],img_paso[k+width],img_paso[k-width]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[k],img_paso[k+width+1],img_paso[k-width-1]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[k],img_paso[k+width-1],img_paso[k-width+1]);
          if(temp>paso) paso=temp;
          if(paso!=img_out[k]) img_out[k]=paso;
        }
      }

      #pragma omp for nowait
      for(int k=0;k<b_size;k++){ // IMAGE BOUNDARY GRADIENT ESTIMATION
          m=c+b[k][0];
          U paso=ami_min3(img_paso[m],img_paso[c+b[k][2]],img_paso[c+b[k][1]]);
          temp=ami_min3(img_paso[m],img_paso[c+b[k][3]],img_paso[c+b[k][4]]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[m],img_paso[c+b[k][5]],img_paso[c+b[k][8]]);
          if(temp>paso) paso=temp;
          temp=ami_min3(img_paso[m],img_paso[c+b[k][6]],img_paso[c+b[k][7]]);
          if(temp>paso) paso=temp;
          if(paso!=img_out[m]) img_out[m]=paso;
      }
    }
  }
}

#endif
