/*
    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_snakes.h \link snakes.h \endlink
 * \brief \ref code_snakes.cpp "See snakes.cpp"
 * \brief Code of  \ref code_snakes.h "snakes.h"
 * \includelineno snakes.h
 * \brief \b Code: \n
 * - \ref code_snakes.cpp "snakes.cpp" \n
 * - \ref code_snakes.h "snakes.h"
 */

/**
 * \file snakes.h
 * \brief Class for basic snakes
 * \author Luis Alvarez \n \n
*/
#ifndef SNAKE_H
#define SNAKE_H

#ifndef NOMINMAX
#define NOMINMAX
#endif


#include "level_sets.h"

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


//#define SNAKE_DEBUG

//#define SNAKE_TIME_DEBUG

#ifdef SNAKE_TIME_DEBUG
#include <time.h>
#endif

// MACROS TO COMPUTE THE min/max OF 3 VALUES
#define min_3(a,b,c) (a<b?(a<c?a:c):(b<c?b:c))
#define max_3(a,b,c) (a>b?(a>c?a:c):(b>c?b:c))

namespace ami {

/**
 * \class snakes
 * \brief Class to implement snakes procedures
 * \author Luis Alvarez
 */
class snakes{

  public:

  template <class T,class U>
  void morphological_snakes_fast_marching(
    ami::image<U> &img,
    ami::contours_collection<T> &cc,
    const int snake_type,
    const int border_size,
    const float sigma,
    const int Niter,
    const int snake_ballon,
    const float snake_ballon_difference_radius,
    U snake_ballon_reference_value,
    const float g_I_threshold_speed_propagation=0.5
  );

  template <class T,class U>
  void morphological_snakes_fast_marching(
    ami::image<U> &img,
    ami::contours_collection<T> &cc,
    const int snake_type,
    const int border_size,
    const float sigma,
    const int Niter,
    const int snake_ballon,
    const float edge_threshold_ballon,
    const float edge_threshold_propagation
  );


  template <class T,class U>
  void morphological_snakes_gauss_regularization(
    ami::image<U> &img,
    ami::contours_collection<T> &cc,
    const int snake_type=0,
    const int border_size=10,
    const float sigma=2.,
    const int Niter=20,
    const int snake_ballon=0,
    const float snake_ballon_difference_radius=20.,
    U snake_ballon_reference_value=(U) 0,
    const float sigma2=2.
  );

  template <class T,class U>
  int fast_marching_dilation(
    ami::contours<T> &contour,
    ami::image<unsigned char> &snake_image,
    ami::image<int>  &snake_img_iteration_control,
    const ami::image<float>  &g_I_img,
    const float snake_ballon_difference_radius=20.,
    U snake_ballon_reference_value=(U) 0
  );

 template <class T,class U>
  int fast_marching_erosion(
    ami::contours<T> &contour,
    ami::image<unsigned char> &snake_image,
    ami::image<int>  &snake_img_iteration_control,
    const ami::image<float>  &g_I_img,
    const float snake_ballon_difference_radius=20.,
    U snake_ballon_reference_value=(U) 0
  );

  template <class T>
  int fast_marching_erosion(
    ami::contours<T> &contour,
    ami::image<unsigned char> &snake_image,
    ami::image<int>  &snake_img_iteration_control,
    const ami::image<float>  &g_I_img,
    const float g_I_threshold_speed_propagation
  );
  template <class T>
  int fast_marching_dilation(
    ami::contours<T> &c,
    ami::image<unsigned char> &snake_img,
    ami::image<int>  &snake_img_iteration_control,
    const ami::image<float>  &g_I_img,
    const float g_I_threshold_speed_propagation);

  template <class T>
  int fast_marching_line_erosion(
    ami::contours<T> &contour,
    ami::image<unsigned char> &snake_image,
    ami::image<int>  &snake_img_iteration_control
  );

  template <class T>
  int fast_marching_line_dilation(
    ami::contours<T> &contour,
    ami::image<unsigned char> &snake_image,
    ami::image<int>  &snake_img_iteration_control
  );

  template <class T>
  int fast_marching_speed_propagation(
    ami::contours<T> &c,
    ami::image<unsigned char> &snake_img,
    ami::image<int>  &snake_img_iteration_control,
    const ami::image<float>  &g_I_img,
    const float g_I_threshold_speed_propagation,
    const ami::image<float>  &g_I_img_x,
    const ami::image<float>  &g_I_img_y
  );


  template <class T>
  float contour_distance(
    const ami::contours<T> &c1,
    const ami::contours<T> &c2,
    ami::image<int>  &snake_img_iteration_control,
    ami::point2d<T> OffSet
  );

  template <class T,class U>
  U ballon_reference_value_estimation(
    const ami::contours<T> &contour,
    const ami::image<U> &img,
    const int snake_ballon
  );
};

}

/**
* \fn template <class T,class U> U ami::snakes::ballon_reference_value_estimation(const ami::contours<T> &contour,const ami::image<U> &img,const int snake_ballon)
* \brief function to compute the ballon reference value according with snake_ballon sign. If snake_ballon>0 (expansion) we compute
the median value in a inside contour neightborhood. If snake_ballon<0 (contraction) we compute the median value in
a outside contour neighborhood.
* \author Luis Alvarez
*/
template <class T,class U>
U ami::snakes::ballon_reference_value_estimation(
const ami::contours<T> &contour /** input contour */,
const ami::image<U> &img /** input image where the contour is defined */,
const int snake_ballon /** snake_ballon sign */
){
  int width=img.width();
  int height=img.height();
  if(width==0 || height==0 || snake_ballon==0) return((U)-1.);


  // WE ESTIMATE TRANSPARENCY IMAGE
  point2d<T> OffSet;
  level_sets ls;
  ami::image<unsigned char> img_t=ls.transparency_image(contour,OffSet);
  int width_t=img_t.width();
  int height_t=img_t.height();
  if(width_t==0 || height_t==0) return((U)-1.);


  // WE DEFINE AUXILIARY VECTOR TO ESTIMATE THE MEDIAN
  vector<U> v;
  // WE COMPUTE THE VECTOR WHEN snake_ballon>0
  if(snake_ballon>0){
    for(int m=0;m<contour.size();m++){
      int i=contour[m].y;
      int j=contour[m].x;
      for(int k=-2;k<=2;k++){
        for(int l=-2;l<2;l++){
          int i2=i+k;
          int j2=j+l;
          if(i2<0 || i2>=height || j2<0 || j2>=width) continue;
          int it=i2-OffSet.y;
          int jt=j2-OffSet.x;
          if(it<0 || it>=height_t || jt<0 || jt>=width_t) continue;
          if(img_t[it*width_t+jt]==255){
            v.push_back(img[i*width+j]);
            img_t[it*width_t+jt]=128;
          }
        }
      }
    }
  }

  // WE COMPUTE THE VECTOR WHEN snake_ballon<0
  if(snake_ballon<0){
    for(int m=0;m<contour.size();m++){
      int i=contour[m].y;
      int j=contour[m].x;
      for(int k=-2;k<=2;k++){
        for(int l=-2;l<2;l++){
          int i2=i+k;
          int j2=j+l;
          if(i2<0 || i2>=height || j2<0 || j2>=width) continue;
          int it=i2-OffSet.y;
          int jt=j2-OffSet.x;
          if(it<0 || it>=height_t || jt<0 || jt>=width_t) continue;
          if(img_t[it*width_t+jt]==0){
            v.push_back(img[i*width+j]);
            img_t[it*width_t+jt]=128;
          }
        }
      }
    }
  }

  // WE COMPUTE THE MEDIAN
  ami::utilities u;
  return(u.median(v,v.size()/2));

}
/**
 * \fn template <class T> float ami::snakes::contour_distance(
const ami::contours<T> &c1,
const ami::contours<T> &c2,
ami::image<int>  &snake_img_iteration_control,
ami::point2d<T> OffSet)
 * \brief Function to compute the relative distance between 2 contours. We perform the addition of distances divided by the contour size.
 * \return The contour distance
 * \author Luis Alvarez
 */
template <class T>
float ami::snakes::contour_distance(
const ami::contours<T> &c1 /** first contour */,
const ami::contours<T> &c2 /** second contour */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image where the contour c1 is drawing to speed up the procedure */,
ami::point2d<T> OffSet /** OffSet distance of contour coordinates with respect to snake_img_iteration_control image coordinates */
)
{
  // WE CHECK THE  DATA
  if(c1.size()==0 || c2.size()==0) return(0.);

  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img_iteration_control.width();
  int height2=snake_img_iteration_control.height();

  // WE DRAW THE CONTOUR c1 IN THE snake_img_iteration_control IMAGE
  for(int l=0,l_end=c1.size();l<l_end;l++){
    int m=(c1[l].y-OffSet.y)*width2+c1[l].x-OffSet.x;
    snake_img_iteration_control[m]=iter_fast_marching;
  }

  // WE COMPUTE THE DISTANCE BETWEEN c1 AND c2
  float distance=0;
  for(int l=0,width2_1=width2-1,height2_1=height2-1,l_end=c2.size();l<l_end;l++){
    int x=c2[l].x-OffSet.x;
    int y=c2[l].y-OffSet.y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]==iter_fast_marching) // point distance 0
      continue;
    if(snake_img_iteration_control[x<width2_1?m+1:m]==iter_fast_marching ||
       snake_img_iteration_control[x>0?m-1:m]==iter_fast_marching ||
       snake_img_iteration_control[y>0?m-width2:m]==iter_fast_marching ||
       snake_img_iteration_control[y<height2_1?m+width2:m]==iter_fast_marching){ //point distance 1
        distance+=1.;
        continue;
    }
    else{ // if the point is not to distance 0 or 1 we fit the distance to 2.
      distance+=2.;
    }
  }
  return(distance/c2.size());

}

/**
 * \fn template <class T> int ami::snakes::fast_marching_speed_propagation(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control,
const ami::image<float>  &g_I_img,
const float g_I_threshold_speed_propagation,
const ami::image<float>  &g_I_img_x,
const ami::image<float>  &g_I_img_y)
 * \brief Morphological fast marching procedure to move a contour using a speed function criterium
 * \author Luis Alvarez
 */

template <class T>
int ami::snakes::fast_marching_speed_propagation(
ami::contours<T> &c  /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image to control fast marching iteration step*/,
const ami::image<float>  &g_I_img /** g(I) : function to minimize. The contour moves to minimize g(I) */,
const float g_I_threshold_speed_propagation /** threshold of function g(I) to make evolve the contour g(I)<threshold */,
const ami::image<float>  &g_I_img_x /** x-derivate of g(I) */ ,
const ami::image<float>  &g_I_img_y /** y-derivative of g(I) */)
{
  // WEIGTH VARIABLES TO COMPUTE THE snake_img DEERIVATIVES
  double coef1,coef2,c1,d1;
  coef1=sqrt((double) 2.);
  coef2=0.25*(2.-coef1);
  coef1=0.5*(coef1-1);

  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img.width();
  int height2=snake_img.height();

  ami::contours<T> c3; // contour to store snake_image modified points

  int width2_2=width2-3,height2_2=height2-3;

/*  #pragma omp parallel for default(none) private(m,y,x) \
  shared(c_aux,size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)*/
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(g_I_img[m]<g_I_threshold_speed_propagation){ // WE CHECK IF g(I) IS UNDER THE THRESHOLD VALUE TO ATTRACT THE CONTOUR
          if(snake_img_iteration_control[m]!=iter_fast_marching){ // WE CHECK IF THE POINT HAS BEEN VISITED
            snake_img_iteration_control[m]=iter_fast_marching;
            // WE COMPUTE THE LEVEL SET IMAGE snake_img DERIVATIVES
            c1=snake_img[m+width2+1]-snake_img[m-width2-1];
            d1=snake_img[m-width2+1]-snake_img[m+width2-1];
            double snake_img_y= coef1*(snake_img[m+width2]- snake_img[m-width2])+coef2*(c1-d1);
            double snake_img_x= -(-coef1*(snake_img[m+1]- snake_img[m-1])+coef2*(c1+d1));
            // WE COMPUTE SPEED CRITERIUM TO DECIDE IF LEVEL SET SHOULD CHANGE
            double scalar_product=g_I_img_x[m]*snake_img_x+g_I_img_y[m]*snake_img_y;
            if(scalar_product<0 && snake_img[m]==255){
              c3.push_back(ami::point2d<T>(x,y));  // we keep coordinates of point where snake_img should change
            }
            if(scalar_product>0 && snake_img[m]==0 ){
              c3.push_back(ami::point2d<T>(x,y)); // we keep coordinates of point where snake_img should change
            }
          }
        }
      }
    }
  }
  //printf("c3.size()=%d\n",c3.size());

  // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  snake_img_iteration_control[0]++; // WE CHANGE VISITED POINT CONTROL VALUE
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    if(snake_img[m]==0) snake_img[m]=255;
    else if(snake_img[m]==255) snake_img[m]=0;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}

/**
 * \fn template <class T> int ami::snakes::fast_marching_line_dilation(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control)
 * \brief Morphological fast marching procedure to perform line dilation regularization
 * \author Luis Alvarez
 */

template <class T>
int ami::snakes::fast_marching_line_dilation(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control/** auxiliary image to control fast marching iteration step*/ )
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img.width();
  int height2=snake_img.height();

  ami::contours<T> c3; // contour to store snake_image modified points

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
  shared(c_aux,size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)*/
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if the point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0){ // erosion only deals with 0 level set value
		    // WE PERFORM LINE DILATION OPERATION IN THE PIXEL
            unsigned char paso=max_3(snake_img[m],snake_img[m+1],snake_img[m-1]);
            if(paso==0) break;
            paso=max_3(snake_img[m],snake_img[m+width2],snake_img[m-width2]);
            if(paso==0) break;
            paso=max_3(snake_img[m],snake_img[m+width2+1],snake_img[m-width2-1]);
            if(paso==0) break;
            paso=max_3(snake_img[m],snake_img[m+width2-1],snake_img[m-width2+1]);
            if(paso==255){
               c3.push_back(ami::point2d<T>(x,y)); // we keep coordinates of points where snake_img should change
            }
          }
        }
      }
    }
  }

  // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    if(snake_img[m]==0) snake_img[m]=255;
    else if(snake_img[m]==255) snake_img[m]=0;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}


/**
 * \fn template <class T> int ami::snakes::fast_marching_line_erosion(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control)
 * \brief Morphological fast marching procedure to perform line erosion regularization
 * \author Luis Alvarez
 */

template <class T>
int ami::snakes::fast_marching_line_erosion(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control/** auxiliary image to control fast marching iteration step*/ )
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img.width();
  int height2=snake_img.height();

  ami::contours<T> c3; // contour to store snake_image modified points

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
 shared(c_aux,size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value) */
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if the point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==255){ // erosion only deals with 255 level set value
            // WE PERFORM LINE DILATION OPERATION IN THE PIXEL
            unsigned char paso=min_3(snake_img[m],snake_img[m+1],snake_img[m-1]);
            if(paso==255) break;
            paso=min_3(snake_img[m],snake_img[m+width2],snake_img[m-width2]);
            if(paso==255) break;
            paso=min_3(snake_img[m],snake_img[m+width2+1],snake_img[m-width2-1]);
            if(paso==255) break;
            paso=min_3(snake_img[m],snake_img[m+width2-1],snake_img[m-width2+1]);
            if(paso==0){
               c3.push_back(ami::point2d<T>(x,y)); // we keep coordinates of points where snake_img should change
            }
          }
        }
      }
    }
  }

  // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    snake_img[m]=0;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}


/**
 * \fn template <class T,class U> int ami::snakes::fast_marching_dilation(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control,
const ami::image<float>  &g_I_img,
const float snake_ballon_difference_radius,
U snake_ballon_reference_value)
 * \brief Morphological fast marching procedure to perform dilation. To propagate the contour while the value gauss convolution of the image be in a given interval
 * \author Luis Alvarez
 */
template <class T,class U>
int ami::snakes::fast_marching_dilation(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image to control fast marching iteration step*/,
const ami::image<float>  &gauss_img /** convolution of the original image */,
const float snake_ballon_difference_radius /** radius of gauss_img[] value interval to propagate the contour */,
U snake_ballon_reference_value /** center of gauss_img[] value interval to propagate the contour */)
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];


  // WE COMPUTE BALLON PROPAGATION
  ami::contours<T> c3; // contour to store snake_image modified points
  int width2=snake_img.width();
  int height2=snake_img.height();
  float max_value=(float) snake_ballon_reference_value+snake_ballon_difference_radius;
  float min_value=(float) snake_ballon_reference_value-snake_ballon_difference_radius;

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
  shared(c_aux,size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value) */
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if the point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          // WE PERFORM DILATION OPERATION
          if(gauss_img[m]<max_value && gauss_img[m]>min_value && snake_img[m]==0 &&
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
            || snake_img[m+width2]==255 ) ){
              c3.push_back(ami::point2d<T>(c[l].x+j,c[l].y+i));// we keep coordinates of points where snake_img should change
          }
        }
      }
    }
  }

  // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    snake_img[m]=255;
  }

 // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}


/**
 * \fn template <class T,class U> int ami::snakes::fast_marching_erosion(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control,
const ami::image<float>  &g_I_img,
const float snake_ballon_difference_radius,
U snake_ballon_reference_value)
 * \brief Morphological fast marching procedure to perform erosion. To propagate the contour while the value gauss convolution of the image be in a given interval
 * \author Luis Alvarez
 */
template <class T,class U>
int ami::snakes::fast_marching_erosion(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image to control fast marching iteration step*/,
const ami::image<float>  &gauss_img /** convolution of the original image */,
const float snake_ballon_difference_radius /** radius of gauss_img[] value interval to propagate the contour */,
U snake_ballon_reference_value /** center of gauss_img[] value interval to propagate the contour */)
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  // WE COMPUTE BALLON PROPAGATION
  ami::contours<T> c3; // contour to store snake_image modified points
  int width2=snake_img.width();
  int height2=snake_img.height();
  float max_value=(float) snake_ballon_reference_value+snake_ballon_difference_radius;
  float min_value=(float) snake_ballon_reference_value-snake_ballon_difference_radius;

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
  shared(c_aux,size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value) */
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if the point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          // WE PERFORM DILATION OPERATION
          if(gauss_img[m]<max_value && gauss_img[m]>min_value && snake_img[m]==255 &&
            (snake_img[m+1]==0 || snake_img[m-1]==0 || snake_img[m-width2]==0
            || snake_img[m+width2]==0 ) ){
              c3.push_back(ami::point2d<T>(x,y));// we keep coordinates of points where snake_img should change
          }
        }
      }
    }
  }

  // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    snake_img[m]=0;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}




/**
 * \fn template <class T> int ami::snakes::fast_marching_erosion(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control,
const ami::image<float>  &g_I_img,
const float g_I_threshold_speed_propagation)
 * \brief Morphological fast marching procedure to perform erosion. We propagate the contour while g_I_img>g_I_threshold_speed_propagation
 * \author Luis Alvarez
 */
template <class T>
int ami::snakes::fast_marching_erosion(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image to control fast marching iteration step*/,
const ami::image<float>  &g_I_img /** g(I) : function to minimize. The contour moves to minimize g(I) */,
const float g_I_threshold_speed_propagation /** threshold of function g(I) to make evolve the contour g(I)>threshold */)
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img.width();
  int height2=snake_img.height();

  // WE COMPUTE BALLON PROPAGATION
  ami::contours<T> c3; // contour to store snake_image modified points

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
 shared(c_aux,size,float_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)*/
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){
          snake_img_iteration_control[m]=iter_fast_marching;
          if(g_I_img[m]>g_I_threshold_speed_propagation && snake_img[m]==255 &&
            (snake_img[m+1]==0 || snake_img[m-1]==0 || snake_img[m-width2]==0
            || snake_img[m+width2]==0 ) ){
              c3.push_back(ami::point2d<T>(x,y));
          }
        }
      }
    }
  }

 // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    snake_img[m]=0;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}


/**
 * \fn template <class T> int ami::snakes::fast_marching_dilation(
ami::contours<T> &c,
ami::image<unsigned char> &snake_img,
ami::image<int>  &snake_img_iteration_control,
const ami::image<float>  &g_I_img,
const float g_I_threshold_speed_propagation)
 * \brief Morphological fast marching procedure to perform dilation. We propagate the contour while g_I_img>g_I_threshold_speed_propagation
 * \author Luis Alvarez
 */
template <class T>
int ami::snakes::fast_marching_dilation(
ami::contours<T> &c /** contour to evolve */,
ami::image<unsigned char> &snake_img /** binary level set image */,
ami::image<int>  &snake_img_iteration_control /** auxiliary image to control fast marching iteration step*/,
const ami::image<float>  &g_I_img /** g(I) : function to minimize. The contour moves to minimize g(I) */,
const float g_I_threshold_speed_propagation /** threshold of function g(I) to make evolve the contour g(I)>threshold */)
{
  // TO CONTROL VISITED POINTS WE USE iter_fast_marching AND snake_img_iteration_control
  snake_img_iteration_control[0]++;
  int iter_fast_marching=snake_img_iteration_control[0];

  int width2=snake_img.width();
  int height2=snake_img.height();

  // WE COMPUTE BALLON PROPAGATION
  ami::contours<T> c3; // contour to store snake_image modified points

  int width2_2=width2-3,height2_2=height2-3;
/*  #pragma omp parallel for default(none) private(m,y,x) \
  shared(c_aux,size,float_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)*/
  for(int l=0;l<c.size();l++){
    for(int i=-1;i<2;i++){
      int y=c[l].y+i;
      if(y<2 || y>height2_2) continue;
      for(int j=-1;j<2;j++){
        int x=c[l].x+j;
        if(x<2 || x>width2_2) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){
          snake_img_iteration_control[m]=iter_fast_marching;
          if(g_I_img[m]>g_I_threshold_speed_propagation && snake_img[m]==0 &&
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
            || snake_img[m+width2]==255 ) ){
              c3.push_back(ami::point2d<T>(x,y));
          }
        }
      }
    }
  }

 // WE CHECK IF THERE IS NOT CHANGE
  if(c3.size()==0) return(0);

  // WE CHANGE VISITED POINT CONTROL VALUE
  snake_img_iteration_control[0]++;
  iter_fast_marching=snake_img_iteration_control[0];

  // WE MODIFY LEVEL SET IMAGE snake_img IN c3 CONTOUR POINTS
  for(int l=0;l<c3.size();l++){
    int m=c3[l].y*width2+c3[l].x;
    snake_img[m]=255;
  }

  // WE UPDATE THE CONTOUR JOINING NEW CONTOUR c3 AND PREVIOUS CONTOUR c
  ami::contours<T> c2; // auxiliary contour to update c
  for(int l=0,width2_1=width2-1,height2_1=height2-1;l<c.size();l++){
    int x=c[l].x;
    int y=c[l].y;
    int m=y*width2+x;
    if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
      snake_img_iteration_control[m]=iter_fast_marching ;
      if(snake_img[m]==0 &&  // we check if we deal with a level set edge point
            (snake_img[x<width2_1?m+1:m]==255 || snake_img[x>0?m-1:m]==255 || snake_img[y>0?m-width2:m]==255
              || snake_img[y<height2_1?m+width2:m]==255 ) ){
                c2.push_back(c[l]);
      }
    }
  }
  for(int l=0/*,width2_1=width2-1,height2_1=height2-1*/;l<c3.size();l++){
    for(int i=-1;i<2;i++){
      int y=c3[l].y+i;
      //if(y<0 || y>height2_1) continue;
      for(int j=-1;j<2;j++){
        int x=c3[l].x+j;
        //if(x<0 || x>width2_1) continue;
        int m=y*width2+x;
        if(snake_img_iteration_control[m]!=iter_fast_marching){ // we check if point has been visited
          snake_img_iteration_control[m]=iter_fast_marching;
          if(snake_img[m]==0 && // we check if we deal with a level set edge point
            (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
              || snake_img[m+width2]==255 ) ){
                c2.push_back(ami::point2d<T>(x,y));  // we add the contour point
          }
        }
      }
    }
  }
  // WE UPDATE THE CONTOUR
  c=c2;
  return(c3.size());
}



/**
 * \fn template <class T,class U> void ami::snakes::morphological_snakes_fast_marching(
ami::image<U> &img,
ami::contours_collection<T> &cc,
const int snake_type,
const int border_size,
const float sigma,
const int Niter,
const int snake_ballon,
const float snake_ballon_difference_radius,
U snake_ballon_reference_value,
const float g_I_threshold_speed_propagation)
 * \brief Morphological snake procedure using a fast marching strategy to speed up the procedure. Dilations and erosion are perfomed using an interval value of image
 * \author Luis Alvarez
 */
template <class T,class U>
void ami::snakes::morphological_snakes_fast_marching(
ami::image<U> &img /** INPUT IMAGE */,
ami::contours_collection<T> &cc /** INPUT/OUTPUT CONTOUR COLLECTION. EACH SNAKE ITERATION GENERATES A NEW CONTOUR
                                    ADDED TO THE ORIGINAL cc CONTOUR COLLECTION. cc[0] IS THE INITIAL CONTOUR */,
const int snake_type /**=0 MEANS SNAKE LOOKS FOR CENTER LINES. =1 MEANS SNAKES LOOKS FOR REGION BORDERS */,
const int border_size /** BORDER SIZE ADDED TO THE CONTOUR IMAGE TO ALLOW SNAKE PROPAGATION */,
const float sigma /** STANDARD DEVIATION OF GAUSSIAN TO CONVOLVE THE ORIGINAL IMAGE */,
const int Niter /** NUMBER OF ITERATIONS */,
const int snake_ballon /**=1 MEANS DILATATION, =-1 MEANS EROSION*/,
const float snake_ballon_difference_radius /** DIFFERENC MAGNITUDE WITH RESPECT TO THE CENTRAL PIXEL VALUE
                                                TO PROPAGATE THE BALLON */,
U snake_ballon_reference_value /** SNAKE BALLON REFERENCE VALUE TO TAKE DIFFERENCE FOR PROPAGATION */,
const float g_I_threshold_speed_propagation /** THRESHOLD FOR THE g(I) VALUES IN SPEED PROPAGATION PROCEDURE */)
{

#ifdef SNAKE_TIME_DEBUG
clock_t  start,end;
start=clock();
#endif

  // WE CHECK THE NUMBER OF POINTS
  if(cc.size()==0 || Niter==0){  return; }

  // WE KEEP THE LAST CONTOUR
  ami::contours<T> c=cc[cc.size()-1];

  // WE KEEP THE ORIGINAL CONTOUR COLLECTION SIZE
  int cc_original_size=cc.size();

  // COORDINATES OF SUBWINDOW CORNERS TO PROCESS THE c
  int i0,j0,i1,j1;

  // WE COMPUTE THE MIN AND MAX OF THE c POINTS
  ami::point2d<T> p2d_min=c.min();
  ami::point2d<T> p2d_max=c.max();

  // WE COMPUTE THE img SIZE AND THE INITIAL POINT COORDINATES
  i0=p2d_min.y-border_size;
  j0=p2d_min.x-border_size;
  i1=p2d_max.y+border_size+1;
  j1=p2d_max.x+border_size+1;
  if(i0<0) i0=0;
  if(i1>=img.height()) i1=img.height();
  if(j0<0) j0=0;
  if(j1>=img.width()) j1=img.width();
  //cout << "i0: " << i0 << " j0: " << j0 << " i1: " << i1 << " j1: " << j1 << endl;
  img.set_roi(j0,j1,i0,i1);

  // WE DEFINE SUBWINDOW DIMENSIONS
  int width2=j1-j0;
  int height2=i1-i0;
  int size=width2*height2;

  //cout << "width2: " << width2 << " height2: " << height2 << " size: " << size <<endl;

  // WE DEFINE AND AUXILIARY LOCAL IMAGE
  ami::image<float> g_I_img;
  ami::image<float> g_I_img_x;
  ami::image<float> g_I_img_y;
  ami::image<float> snake_img_x;
  ami::image<float> snake_img_y;
//  ami::image<unsigned char>  snake_img(width2,height2,1,255);
//  ami::image<int>  snake_img_iteration_control(width2,height2,1,0);
//
//  // WE INITILIAZE VISITED POINT REFERENCE VALUE
//  snake_img_iteration_control[0]=0;


  // WE COMPUTE THE IMAGE CONVOLUTION AND DERIVATIVES
  ami::filters f;
  f.gauss_conv(img,g_I_img,sigma,sigma);
  f.grad(g_I_img,g_I_img_x,g_I_img_y,(bool) 1);

  ami::image<float> img_gauss=g_I_img;

  if(snake_type==1){
    //#pragma omp parallel for default(none) private(m) shared(size,g_I_img2,g_I_img_x,g_I_img_y)
    for(int m=0;m<size;++m){
      g_I_img[m]=1./sqrt(1e-20+g_I_img_x[m]*g_I_img_x[m]+g_I_img_y[m]*g_I_img_y[m]);
    }
    f.grad(g_I_img,g_I_img_x,g_I_img_y,(bool) 1);
  }

  // WE REMOVE TEMPORAL ROI
  img.roi_clear();

  // WE FILL THE AUXILIARY CONTOUR
  ami::contours<T> c_aux(c.size()); // translated contour to the subwindow
  ami::contours<T> c_aux2,c_aux3;
  c_aux[0]=c[0]-point2d<T>((T)j0,(T) i0);
  int m2=0;
  for(int m=0;m<c.size();++m){
    int i=(int) (c[m].y-i0);
    int j=(int) (c[m].x-j0);
    //cout << "(i,j) : " << i << " " << j << endl;
    //snake_img[i*width2+j]= 0;
    if(j!=c_aux[m2].x || i!=c_aux[m2].y)
    c_aux[++m2]=point2d<T>((T)j,(T)i);
  }
  c_aux.resize(m2+1);


  // WE BUILD SNAKE IMAGE AD WE FILL IT WITH A VALUE
   ami::image<unsigned char>  snake_img(width2,height2,1,2); // 255
  // WE DRAW THE CONTOUR IN THE IMAGE
  for(int m=0;m<c_aux.size();++m){
    int i=(int) (c_aux[m].y);
    int j=(int) (c_aux[m].x);
    //cout << "(i,j) : " << i << " " << j << endl;
    snake_img[i*width2+j]= 1; // 0
  }

  // WE FILL THE CONTOUR OUTSIDE
  ami::utilities u;
  int contador=u.scan_line_filling(snake_img,0,0, (int) 0,(unsigned char) 2, (unsigned char) 0); // 255 0

  // WE LOOK FOR HOLES
  for(int m=0;m<c_aux.size();++m){
    int i=(int) (c_aux[m].y);
    int j=(int) (c_aux[m].x);
    // WE ANALYZE NEIGBORHOOD
    int i3,j3;
    int cont0=0,cont255=0,cont2=0;
    for(int di=-1;di<2;di++){
      for(int dj=-1;dj<2;dj++){
        int i2=i+di;
        int j2=j+dj;
        if(i2>=0 && i2<height2 && j2>=0 && j2<width2){
          int m=i2*width2+j2;
          if(snake_img[m]==2) {cont2++; i3=i2; j3=j2;}
          else if(snake_img[m]==0) cont0++;
          else if(snake_img[m]==255) cont255++;
        }
      }
      // WE FILL HOLES
      if(cont2>0 && cont0>0)
        contador=u.scan_line_filling(snake_img,j3,i3, (int) 0,(unsigned char) 2, (unsigned char) 255);
      else if(cont2>0 && cont255>0)
        contador=u.scan_line_filling(snake_img,j3,i3, (int) 0,(unsigned char) 2, (unsigned char) 0);
    }
  }

  // WE FILL TO ZERO THE CONTOUR
  for(int m=0;m<c_aux.size();++m){
    int i=(int) (c_aux[m].y);
    int j=(int) (c_aux[m].x);
    snake_img[i*width2+j]=0;
  }


  // WE DEFINE AN AUXILIARY IMAGE TO CONTROL VISITED POINT REFERENCE VALUE
  ami::image<int>  snake_img_iteration_control(width2,height2,1,3);
  // WE INITIALIZE VISITED POINT REFERENCE VALUE
  snake_img_iteration_control[0]=0;



  // WE COMPUTE THRESHOLD VALUES FOR g(I)
  float g_I_threshold_speed_propagation_level;
  if(snake_type==0){
    g_I_threshold_speed_propagation_level=1e20;
  }
  else if(snake_type==1){
    if(g_I_threshold_speed_propagation<=0){
      g_I_threshold_speed_propagation_level=-1.;
    }
    else if (g_I_threshold_speed_propagation>=1){
      g_I_threshold_speed_propagation_level=1e20;
    }
    else{
      int level0=g_I_threshold_speed_propagation*width2*height2;
      g_I_threshold_speed_propagation_level=u.median(*g_I_img.get_image(),level0);
    }
  }

  #ifdef SNAKE_DEBUG
    printf("g_I_threshold_speed_propagation_level=%e\n",g_I_threshold_speed_propagation_level);
  #endif


  //WE GO TROUGH ITERATIONS TO COMPUTE THE EVOLUTION CONTOURS
  int ntries=0; //variable to control number of tries in convergence criterium
  for(int k=0;k<Niter;k++){
    cc.resize(cc_original_size+k+1); // we take memory for the new contour

    int NumberModifiedPoints; // we will compute the number of modified points after each iteration

    // WE COMPUTE BALLON PROPAGATION
    if(snake_ballon==1){
      NumberModifiedPoints=fast_marching_dilation(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        img_gauss,
        snake_ballon_difference_radius,
        snake_ballon_reference_value
      );
      #ifdef SNAKE_DEBUG
        printf("NumberModifiedPointsBallon+1=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
   }
   if(snake_ballon==-1){
      NumberModifiedPoints=fast_marching_erosion(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        img_gauss,
        snake_ballon_difference_radius,
        snake_ballon_reference_value
      );
      #ifdef SNAKE_DEBUG
        printf("NumberModifiedPointsBallon-1=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
   }

    // WE COMPUTE THE FRONT PROPAGATION ACCORDING TO SNAKE TYPE EVOLUTION
    if(snake_type<2){
      NumberModifiedPoints=fast_marching_speed_propagation(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        g_I_img,
        g_I_threshold_speed_propagation_level,
        g_I_img_x,
        g_I_img_y
      );
      #ifdef SNAKE_DEBUG
        printf("NumberModifiedPointsSpeedPropagation=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
    }

//    // WE REGULARIZE THE SNAKE
//    NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
    #ifdef SNAKE_DEBUG
      printf("NumberModifiedPointsRegularization=%d\n",NumberModifiedPoints);
    #endif
    // WE STORE THE CONTOUR AT ITERATION k
    cc[cc_original_size+k].resize(c_aux.size());
    point2d<T> OffSet((T) j0,(T) i0);
    for(int l=0,l_end=c_aux.size();l<l_end;l++){
      cc[cc_original_size+k][l]=c_aux[l]+OffSet;
    }

    // WE CHECK CONVERGENCE CRITERIUM
    if(k>1){
      float distance = contour_distance(cc[cc_original_size+k],cc[cc_original_size+k-2],snake_img_iteration_control,OffSet);
      #ifdef SNAKE_DEBUG
         printf("distance = %e\n",distance);
      #endif
      if(distance<0.01 || (distance<0.01 && snake_type==0) ){
         ntries++;
         if(ntries==3) break;
      }
      else if(ntries>0) ntries--;
    }

    #ifdef SNAKE_DEBUG
   {
      static int cont_debug5=10000;
      ami::image<unsigned char> img_debug(width2,height2,3,255);
      int size2=width2*height2;
      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]=(unsigned char) img_gauss[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }
      std::stringstream ss;//create a stringstream
      ss << "snake_gauss_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss.str());

      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]= snake_img[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }

      std::stringstream ss2;//create a stringstream
      ss2 << "snake_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss2.str());
    }
    #endif

//    ami::level_sets ls;
//    int error=ls.clean_contours(cc[cc_original_size+k]);
//    error=ls.clean_contours(cc[cc_original_size+k]);
//    error=ls.clean_contours(cc[cc_original_size+k]);
//    if(error!=0){ cout << "error clean_contours= " << error << endl; system("pause"); }
  }

  #ifdef SNAKE_TIME_DEBUG
  end=clock();
  double dif=(double) (end-start)/CLOCKS_PER_SEC;
  printf("Computational cost of morphological snakes : \n  Number of iterations =%d  time (in seconds) %e \n",cc.size()-cc_original_size,dif);
  #endif

}



/**
 * \fn template <class T,class U> void ami::snakes::morphological_snakes_fast_marching(
ami::image<U> &img,
ami::contours_collection<T> &cc,
const int snake_type,
const int border_size,
const float sigma,
const int Niter,
const int snake_ballon,
const float g_I_threshold_ballon_propagation,
const float g_I_threshold_speed_propagation)
 * \brief Morphological snake procedure using a fast marching strategy to speed up the procedure. Dilations and erosion are perfomed using g(I)>Threshold
 * \author Luis Alvarez
 */
template <class T,class U>
void ami::snakes::morphological_snakes_fast_marching(
ami::image<U> &img /** INPUT IMAGE */,
ami::contours_collection<T> &cc /** INPUT/OUTPUT CONTOUR COLLECTION. EACH SNAKE ITERATION GENERATES A NEW CONTOUR
                                    ADDED TO THE ORIGINAL cc CONTOUR COLLECTION. cc[0] IS THE INITIAL CONTOUR */,
const int snake_type /**=0 MEANS SNAKE LOOKS FOR CENTER LINES. =1 MEANS SNAKES LOOKS FOR REGION BORDERS */,
const int border_size /** BORDER SIZE ADDED TO THE CONTOUR IMAGE TO ALLOW SNAKE PROPAGATION */,
const float sigma /** STANDARD DEVIATION OF GAUSSIAN TO CONVOLVE THE ORIGINAL IMAGE */,
const int Niter /** NUMBER OF ITERATIONS */,
const int snake_ballon /** =TRUE MEANS WE USE SNAKE BALLON PROCEDURE */,
const float g_I_threshold_ballon_propagation /** THRESHOLD FOR THE g(I) VALUES IN BALLON PROPAGATION PROCEDURE */,
const float g_I_threshold_speed_propagation /** THRESHOLD FOR THE g(I) VALUES IN SPEED PROPAGATION PROCEDURE */)
{

#ifdef SNAKE_TIME_DEBUG
clock_t  start,end;
start=clock();
#endif

  // WE CHECK THE NUMBER OF POINTS
  if(cc.size()==0 || Niter==0){  return; }

  // WE KEEP THE LAST CONTOUR
  ami::contours<T> c=cc[cc.size()-1];

  // WE KEEP THE ORIGINAL CONTOUR COLLECTION SIZE
  int cc_original_size=cc.size();

  // COORDINATES OF SUBWINDOW CORNERS TO PROCESS THE c
  int i0,j0,i1,j1;

  // WE COMPUTE THE MIN AND MAX OF THE c POINTS
  ami::point2d<T> p2d_min=c.min();
  ami::point2d<T> p2d_max=c.max();

  // WE COMPUTE THE img SIZE AND THE INITIAL POINT COORDINATES
  i0=p2d_min.y-border_size;
  j0=p2d_min.x-border_size;
  i1=p2d_max.y+border_size+1;
  j1=p2d_max.x+border_size+1;
  if(i0<0) i0=0;
  if(i1>=img.height()) i1=img.height();
  if(j0<0) j0=0;
  if(j1>=img.width()) j1=img.width();
  //cout << "i0: " << i0 << " j0: " << j0 << " i1: " << i1 << " j1: " << j1 << endl;
  img.set_roi(j0,j1,i0,i1);

  // WE DEFINE SUBWINDOW DIMENSIONS
  int width2=j1-j0;
  int height2=i1-i0;
  int size=width2*height2;

  //cout << "width2: " << width2 << " height2: " << height2 << " size: " << size <<endl;

  // WE DEFINE AND AUXILIARY LOCAL IMAGE
  ami::image<float> g_I_img;
  ami::image<float> g_I_img_x;
  ami::image<float> g_I_img_y;
  ami::image<float> snake_img_x;
  ami::image<float> snake_img_y;
  ami::image<unsigned char>  snake_img(width2,height2,1,255);
  ami::image<int>  snake_img_iteration_control(width2,height2,1,0);

  // WE INITIALIZE VISITED POINT REFERENCE VALUE
  snake_img_iteration_control[0]=0;

  // WE COMPUTE THE IMAGE CONVOLUTION AND DERIVATIVES
  ami::filters f;
  f.gauss_conv(img,g_I_img,sigma,sigma);
  f.grad(g_I_img,g_I_img_x,g_I_img_y,(bool) 1);
  //ami::image<float> g_I_img2(width2,height2,1);
  #ifdef SNAKE_DEBUG
   ami::image<float> img_gauss=g_I_img;
  #endif

  if(snake_type==1){
    //#pragma omp parallel for default(none) private(m) shared(size,g_I_img2,g_I_img_x,g_I_img_y)
    for(int m=0;m<size;++m){
      g_I_img[m]=1./sqrt(1e-20+g_I_img_x[m]*g_I_img_x[m]+g_I_img_y[m]*g_I_img_y[m]);
    }
    f.grad(g_I_img,g_I_img_x,g_I_img_y,(bool) 1);
  }

  // WE REMOVE TEMPORAL ROI
  img.roi_clear();

  // WE FILL THE SNAKE IMAGE AND AUXILIARY CONTOUR
  ami::contours<T> c_aux(c.size()); // translated contour to the subwindow
  ami::contours<T> c_aux2,c_aux3;
  c_aux[0]=c[0]-point2d<T>((T)j0,(T) i0);
  int m2=0;
  for(int m=0;m<c.size();++m){
    int i=(int) (c[m].y-i0);
    int j=(int) (c[m].x-j0);
    //cout << "(i,j) : " << i << " " << j << endl;
    snake_img[i*width2+j]= 0;
    if(j!=c_aux[m2].x || i!=c_aux[m2].y)
    c_aux[++m2]=point2d<T>((T)j,(T)i);
  }
  c_aux.resize(m2+1);

  // WE FILL THE CONTOUR OUTSIDE
  ami::utilities u;
  //int contador=
  u.scan_line_filling(snake_img,0,0, (int) 0,(unsigned char) 255, (unsigned char) 0);

  // WE COMPUTE THRESHOLD VALUES FOR g(I)
  float g_I_threshold_ballon_propagation_level;
  if(g_I_threshold_ballon_propagation<=0){
    g_I_threshold_ballon_propagation_level=-1.;
  }
  else if (g_I_threshold_ballon_propagation>=1){
    g_I_threshold_ballon_propagation_level=1e20;
  }
  else{
    int level0=g_I_threshold_ballon_propagation*width2*height2;
    g_I_threshold_ballon_propagation_level=u.median(*g_I_img.get_image(),level0);
  }

  float g_I_threshold_speed_propagation_level;
  if(snake_type==0){
    g_I_threshold_speed_propagation_level=1e20;
  }
  else{
    if(g_I_threshold_speed_propagation<=0){
      g_I_threshold_speed_propagation_level=-1.;
    }
    else if (g_I_threshold_speed_propagation>=1){
      g_I_threshold_speed_propagation_level=1e20;
    }
    else{
      if(g_I_threshold_ballon_propagation==g_I_threshold_speed_propagation){
        g_I_threshold_speed_propagation_level=g_I_threshold_ballon_propagation_level;
      }
      else{
        int level1=g_I_threshold_speed_propagation*width2*height2;
        g_I_threshold_speed_propagation_level=u.median(*g_I_img.get_image(),level1);
      }
    }
  }

  #ifdef SNAKE_DEBUG
    printf("g_I_threshold_ballon_propagation_level=%e\n",g_I_threshold_ballon_propagation_level);
    printf("g_I_threshold_speed_propagation_level=%e\n",g_I_threshold_speed_propagation_level);
  #endif



  //WE GO TROUGH ITERATIONS TO COMPUTE THE EVOLUTION CONTOURS
  int ntries=0; //variable to control number of tries in convergence criterium
  for(int k=0;k<Niter;k++){
    cc.resize(cc_original_size+k+1); // we take memory for the new contour

    int NumberModifiedPoints; // we will compute the number of modified points after each iteration


    // WE COMPUTE BALLON PROPAGATION
    if(snake_ballon==1){
      NumberModifiedPoints=fast_marching_dilation(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        g_I_img,
        g_I_threshold_ballon_propagation_level
      );
      #ifdef SNAKE_DEBUG
      printf("NumberModifiedPointsBallon+1=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
   }

   if(snake_ballon==-1){
      NumberModifiedPoints=fast_marching_erosion(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        g_I_img,
        g_I_threshold_ballon_propagation_level
      );
      #ifdef SNAKE_DEBUG
        printf("NumberModifiedPointsBallon-1=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
   }

    // WE COMPUTE THE FRONT PROPAGATION ACCORDING TO SNAKE TYPE EVOLUTION
    if(snake_type<2){
      NumberModifiedPoints=fast_marching_speed_propagation(
        c_aux,
        snake_img,
        snake_img_iteration_control,
        g_I_img,
        g_I_threshold_speed_propagation_level,
        g_I_img_x,
        g_I_img_y
      );
      #ifdef SNAKE_DEBUG
        printf("NumberModifiedPointsSpeedPropagation=%d\n",NumberModifiedPoints);
      #endif
      // WE REGULARIZE THE SNAKE
      NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
      NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
    }

    // WE REGULARIZE THE SNAKE
//    NumberModifiedPoints=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_dilation(c_aux,snake_img,snake_img_iteration_control);
//    NumberModifiedPoints+=fast_marching_line_erosion(c_aux,snake_img,snake_img_iteration_control);
    #ifdef SNAKE_DEBUG
      printf("NumberModifiedPointsRegularization=%d\n",NumberModifiedPoints);
    #endif

    // WE STORE THE CONTOUR AT ITERATION k
    cc[cc_original_size+k].resize(c_aux.size());
    point2d<T> OffSet((T) j0,(T) i0);
    for(int l=0,l_end=c_aux.size();l<l_end;l++){
      cc[cc_original_size+k][l]=c_aux[l]+OffSet;
    }

    // WE CHECK CONVERGENCE CRITERIUM
    if(k>1){
      float distance = contour_distance(cc[cc_original_size+k],cc[cc_original_size+k-2],snake_img_iteration_control,OffSet);
     #ifdef SNAKE_DEBUG
         printf("distance = %e\n",distance);
      #endif
      if(distance<0.025 || (distance<0.05 && snake_type==0) ){
         ntries++;
         if(ntries==3) return;
      }
      else if(ntries>0) ntries--;
    }
    #ifdef SNAKE_DEBUG
   {
      static int cont_debug5=10000;
      ami::image<unsigned char> img_debug(width2,height2,3,255);
      int size2=width2*height2;
      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]=(unsigned char) img_gauss[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }
      std::stringstream ss;//create a stringstream
      ss << "snake_gauss_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss.str());

      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]= snake_img[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }

      std::stringstream ss2;//create a stringstream
      ss2 << "snake_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss2.str());
    }
    #endif

//    ami::level_sets ls;
//    int error=ls.clean_contours(cc[cc_original_size+k]);
//    error=ls.clean_contours(cc[cc_original_size+k]);
//    error=ls.clean_contours(cc[cc_original_size+k]);
//    if(error!=0){ cout << "error clean_contours= " << error << endl; system("pause"); }
  }

  #ifdef SNAKE_TIME_DEBUG
  end=clock();
  double dif=(double) (end-start)/CLOCKS_PER_SEC;
  printf("Computational cost of morphological snakes : \n  Number of iterations =%d  time (in seconds) %e \n",cc.size()-cc_original_size,dif);
  #endif

}


/**
 * \fn template <class T,class U> void ami::snakes::morphological_snakes_gauss_regularization(
ami::image<U> &img,
ami::contours_collection<T> &cc,
const int snake_type,
const int border_size,
const float sigma,
const int Niter,
const int snake_ballon,
const float snake_ballon_difference_radius,
U snake_ballon_reference_value,
const float sigma2)
 * \brief Morphological snake procedure using a gauss convolution regularization instead of morphological regularization
 * \author Luis Alvarez
 */
template <class T,class U>
void ami::snakes::morphological_snakes_gauss_regularization(
ami::image<U> &img /** INPUT IMAGE */,
ami::contours_collection<T> &cc /** INPUT/OUTPUT CONTOUR COLLECTION. EACH SNAKE ITERATION GENERATES A NEW CONTOUR
                                    ADDED TO THE ORIGINAL cc CONTOUR COLLECTION. cc[0] IS THE INITIAL CONTOUR */,
const int snake_type /**=0 MEANS SNAKE LOOKS FOR CENTER LINES. =1 MEANS SNAKES LOOKS FOR REGION BORDERS */,
const int border_size /** BORDER SIZE ADDED TO THE CONTOUR IMAGE TO ALLOW SNAKE PROPAGATION */,
const float sigma /** STANDARD DEVIATION OF GAUSSIAN TO CONVOLVE THE ORIGINAL IMAGE */,
const int Niter /** NUMBER OF ITERATIONS */,
const int snake_ballon /** =TRUE MEANS WE USE SNAKE BALLON PROCEDURE */,
const float snake_ballon_difference_radius /** DIFFERENC MAGNITUDE WITH RESPECT TO THE CENTRAL PIXEL VALUE
                                                TO PROPAGATE THE BALLON */,
U snake_ballon_reference_value /** SNAKE BALLON REFERENCE VALUE TO TAKE DIFFERENCE FOR PROPAGATION */,
const float sigma2/** . */)
{

  // WE CHECK THE NUMBER OF POINTS
  if(cc.size()==0 || Niter==0){  return; }

  ami::contours<T> c=cc[cc.size()-1];

  // WE KEEP THE ORIGINAL CONTOUR COLLECTION SIZE
  int cc_original_size=cc.size();

  // COORDINATES OF SUBWINDOW CORNERS TO PROCESS THE c
  int i0,j0,i1,j1;

  // VARIABLE TO CONTROL IF A ROI IS DEFINED
  bool roi_window=false;

  // AUXILIARY CONTOUR TO STORE THE CONTOUR IN THE CASE roi_window == true
  ami::contours<T> c2;

  // WE DEFINE SUBWINDOW IF THE ORIGINAL img HAVE NOT DEFINED A SUBWINDOW
  if(img.get_roi().size()<4){
    // WE COMPUTE THE MIN AND MAX OF THE c POINTS
    ami::point2d<T> p2d_min=c.min();
    ami::point2d<T> p2d_max=c.max();

    // WE COMPUTE THE img SIZE AND THE INITIAL POINT COORDINATES
    i0=p2d_min.y-border_size;
    j0=p2d_min.x-border_size;
    i1=p2d_max.y+border_size+1;
    j1=p2d_max.x+border_size+1;
    if(i0<0) i0=0;
    if(i1>=img.height()) i1=img.height();
    if(j0<0) j0=0;
    if(j1>=img.width()) j1=img.width();
    img.set_roi(j0,j1,i0,i1);
  }
  else{ // WE DEFINE SUBWINDOW USING THE ORIGINAL img SUBWINDOW
    roi_window=true;
    j0=img.get_roi().at(0);
    j1=img.get_roi().at(1);
    i0=img.get_roi().at(2);
    i1=img.get_roi().at(3);
  }

  // WE DEFINE SUBWINDOW DIMENSIONS
  int width2=j1-j0;
  int height2=i1-i0;
  int size=width2*height2;

  // WE DEFINE AND AUXILIARY LOCAL IMAGE
  ami::image<float> g_I_img;
  ami::image<float> g_I_img_x;
  ami::image<float> g_I_img_y;
  ami::image<float> snake_img_x;
  ami::image<float> snake_img_y;
  ami::image<unsigned char>  snake_img(width2,height2,1,255);

  // WE COMPUTE THE IMAGE CONVOLUTION AND DERIVATIVES
  ami::filters f;
  f.gauss_conv(img,g_I_img,sigma,sigma);
  f.grad(g_I_img,g_I_img_x,g_I_img_y,(bool) 1);
  if(snake_type==1){
    ami::image<float> g_I_img2(width2,height2,1);
    int m;
    #pragma omp parallel for private(m) shared(size,g_I_img2,g_I_img_x,g_I_img_y)
    for(m=0;m<size;++m){
      g_I_img2[m]=1./sqrt(1e-20+g_I_img_x[m]*g_I_img_x[m]+g_I_img_y[m]*g_I_img_y[m]);
    }

    f.grad(g_I_img2,g_I_img_x,g_I_img_y,(bool) 1);
  }
  // WE REMOVE TEMPORAL ROI
  if(roi_window==false) img.roi_clear();

  // WE FILL THE SNAKE IMAGE
  if(roi_window==false){
    for(int m=0;m<c.size();++m){
      int i=(int) (c[m].y-i0);
      int j=(int) (c[m].x-j0);
      snake_img[i*width2+j]= 0;
    }
  }
  else{
    // WE FILL THE SNAKE IMAGE INSIDE THE SUBWINDOW
    // WE COMPUTE THE CONTOUR POINT CLOSER TO THE WINDOW CENTER
    point2d<T> p2m=point2d<T>((T) (0.5*(j0+j1)),(T) (0.5*(i0+i1)));
    int pos_min,pos_max,pos=0;
    double min=(p2m-c[0]).norm2();
    for(int m=1;m<c.size();m++){
      double norm=(p2m-c[m]).norm2();
      if( norm<min){
        min=norm;
        pos=m;
      }
    }
    // WE COMPUTE BY PROPAGATION IN BOTH DIRECTION THE CONTOURS POINTS INSIDE THE WINDOW
    for(pos_max=pos;pos_max<c.size();++pos_max){
      pos_max=pos_max%(c.size());
      if(pos_max==(pos-1)) break; // ALL THE CONTOUR IS INSIDE THE SUBWINDOW
      int m=pos_max;
      int i=(int) (c[m].y-i0);
      int j=(int) (c[m].x-j0);
      if(i<0 || i>=height2 || j<0 || j>=width2) break;
      snake_img[i*width2+j]=0;
    }

    for(pos_min=pos;pos_min>=0;--pos_min){
      if(pos_min<0) pos_min=c.size()-1;
      if(pos_min==(pos+1)) break; // ALL THE CONTOUR IS INSIDE THE SUBWINDOW
      int m=pos_min;
      int i=(int) (c[m].y-i0);
      int j=(int) (c[m].x-j0);
      if(i<0 || i>=height2 || j<0 || j>=width2) break;
      snake_img[i*width2+j]=0;
    }

    // WE CONSTRUCT A NEW CONTOUR TAKING INTO ACCOUNT JUST THE POINT OUTSIDE THE SUBWINDOW
    if(pos_max!=(pos-1)){
      int m=pos_min;
      while(m!=pos_max){
        int i=(int) (c[m].y-i0);
        int j=(int) (c[m].x-j0);
        c2.push_back(point2d<T>((T) j, (T) i));
        m--;
        if(m<0) m=c.size()-1;
      }
    }
  }

  // WE CONSTRUCT THE 5 NEIGHBORHOOD IMAGE BOUNDARY INDEX VECTOR
  ami::utilities u;
  vector < vector <int> > b=u.boundary_neighborhood_5n(width2,height2);

  // WE FILL THE CONTOUR OUTSIDE
  {
    int k=0;
    for( ;k<b.size();k++){
      if(snake_img[b[k][0]]==255){
        int x2=b[k][0]%width2;
        int y2=b[k][0]/width2;
        int contador=u.scan_line_filling(snake_img,(int) x2, (int) y2, (int) 0,(unsigned char) 255, (unsigned char) 0);
        break;
      }
    }
    if(k==b.size()) {cc.resize(cc_original_size); return;}
  }


  //WE GO TROUGH ITERATIONS
  for(int k=0;k<Niter;k++){
    #ifdef SNAKE_DEBUG
      printf("Niter=%d iteration %d  cc.size()=%d cc[cc.size()-1].size()=%d c.size()=%d\n",
      Niter,k,cc.size(),cc[cc.size()-1].size(),c.size());
      #endif
      cc.resize(cc_original_size+k+1);

    // WE COMPUTE BALLON PROPAGATION
    if(snake_ballon==1){
      float max_value=(float) snake_ballon_reference_value+snake_ballon_difference_radius;
      float min_value=(float) snake_ballon_reference_value-snake_ballon_difference_radius;
      ami::image<unsigned char > snake_img2=snake_img;
      int m,x,y,width2_1=width2-1,height2_1=height2-1;
      #pragma omp parallel for private(m,y,x) \
      shared(size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)
      for(y=1;y<height2_1;++y){
        for(x=1;x<width2_1;++x){
          m=y*width2+x;
          if(g_I_img[m]<max_value && g_I_img[m]>min_value &&
             (snake_img2[m+1]==255 || snake_img2[m-1]==255 || snake_img2[m-width2]==255
             || snake_img2[m+width2]==255 ) ){
               if(snake_img[m]==0) snake_img[m]=255;
          }
        }
      }
      if(roi_window==true){ // WE COMPUTE BALLON PROPAGATION IN THE IMAGE BOUNDARY
        int n,b_size=b.size();
        #pragma omp parallel for private(n) shared(b,b_size,snake_img2,snake_img)
        for(n=0;n<b_size;++n){
          if( (snake_img2[b[n][1]]==255 || snake_img2[b[n][2]]==255 || snake_img2[b[n][3]]==255
             || snake_img2[b[n][4]]==255 ) ){
             if(snake_img[b[n][0]]==0) snake_img[b[n][0]]=255;
          }
        }
      }
    }

    // WE COMPUTE BALLON PROPAGATION
    if(snake_ballon==-1){
      float max_value=256.;
      float min_value=-1.;
      ami::image<unsigned char > snake_img2=snake_img;
      int m,x,y,width2_1=width2-1,height2_1=height2-1;
      #pragma omp parallel for private(m,y,x) \
      shared(size,g_I_img,snake_img,snake_img2,width2_1,height2_1,width2,max_value,min_value)
      for(y=1;y<height2_1;++y){
        for(x=1;x<width2_1;++x){
          m=y*width2+x;
          if(g_I_img[m]<max_value && g_I_img[m]>min_value &&
             (snake_img2[m+1]==0 || snake_img2[m-1]==0 || snake_img2[m-width2]==0
             || snake_img2[m+width2]==0 ) ){

             if(snake_img[m]==255) snake_img[m]=0;
          }
        }
      }
      if(roi_window==true){ // WE COMPUTE BALLON PROPAGATION IN THE IMAGE BOUNDARY
        int n,b_size=b.size();
        #pragma omp parallel for private(n) shared(b,b_size,snake_img2,snake_img)
        for(n=0;n<b_size;++n){
          if( (snake_img2[b[n][1]]==0 || snake_img2[b[n][2]]==0 || snake_img2[b[n][3]]==0
             || snake_img2[b[n][4]]==0 ) ){
             if(snake_img[b[n][0]]==255) snake_img[b[n][0]]=0;
          }
        }
      }
    }


    // WE COMPUTE THE GRADIENT OF THE LEVEL SET img
    if(snake_type<2){
      f.grad(snake_img,snake_img_x,snake_img_y,(bool) 1);
      int m,x,y,width2_1=width2-1,height2_1=height2-1;
      double scalar_product;
      #pragma omp parallel for private(m,y,x,scalar_product) \
      shared(size,g_I_img,snake_img,g_I_img_x,g_I_img_y,snake_img_x,snake_img_y,width2_1,height2_1,width2)
      for(y=1;y<height2_1;++y){
        for(x=1;x<width2_1;++x){
          m=y*width2+x;
          scalar_product=g_I_img_x[m]*snake_img_x[m]+g_I_img_y[m]*snake_img_y[m];
          if(scalar_product==0.) continue;
          if(scalar_product<0 ){ snake_img[m]=0; }
          if(scalar_product>0 ){ snake_img[m]=255; }
        }
      }
    }

    // WE REGULARIZE THE SNAKE
    ami::image<float> gauss_snake_img;
    f.gauss_conv(snake_img,gauss_snake_img,sigma2,sigma2);
    ami::image<unsigned char> snake_img4=snake_img;
    for(int y=1;y<(snake_img.height()-1);++y){
        int n=y*width2;
        for(int x=1;x<(snake_img4.width()-1);++x){
          int m=n+x;
          if(snake_img4[m]==0 && (snake_img4[m+1]==255 || snake_img4[m-1]==255 || snake_img4[m-width2]==255
             || snake_img4[m+width2]==255 ) ){
               if(gauss_snake_img[m]>127.5) snake_img[m]=255;
          }
          if(snake_img4[m]==255 && (snake_img4[m+1]==0 || snake_img4[m-1]==0 || snake_img4[m-width2]==0
             || snake_img4[m+width2]==0 ) ){
               if(gauss_snake_img[m]<127.5) snake_img[m]=0;
          }
        }
    }


    { // WE BUILD THE CONTOUR AT ITERATION k
      cc[cc_original_size+k].clear();
      int n,m,x,y,width2_1=width2-1,height2_1=height2-1;
      double scalar_product;
      #pragma omp parallel for schedule(dynamic) private(m,n,y,x) \
      shared(roi_window,i0,j0,size,cc_original_size,snake_img,width2_1,height2_1,width2,k)
      for(y=1;y<height2_1;++y){
        n=y*width2;
        for(x=1;x<width2_1;++x){
          m=n+x;
          if(snake_img[m]==0 && (snake_img[m+1]==255 || snake_img[m-1]==255 || snake_img[m-width2]==255
             || snake_img[m+width2]==255 ) ){
            if((x==1 || y==1  || y==(height2_1-1) || x==(width2_1-1)) ){
              if(roi_window==false){
                #pragma omp critical (setting_dimensions_0)
                { height2_1=width2_1=0; }
              }
              else continue;
            }
            #pragma omp critical (adding_contour_point)
            {
              cc[cc_original_size+k].push_back(point2d<T>((T) j0+x, (T) i0+y));
            }
          }
        }
      }
      if( height2_1==0){
        cc.resize(cc_original_size+k);
        return;
      }
    }

    if(roi_window==true && c2.size()>3){
      // WE COMPLETE THE CONTOUR
      // WE COMPUTE THE NEAREST NEW CONTOUR POINT TO OLD CONTOUR FIRST POINT
      int pos_1=0;
      double min=1e30;

      for(int m=0;m<cc[cc_original_size+k].size();m++){

        double n_j0=(cc[cc_original_size+k][m].x-j0)*(cc[cc_original_size+k][m].x-j0);
        double n_j1=(cc[cc_original_size+k][m].x-j1)*(cc[cc_original_size+k][m].x-j1);
        double n_i0=(cc[cc_original_size+k][m].y-i0)*(cc[cc_original_size+k][m].y-i0);
        double n_i1=(cc[cc_original_size+k][m].y-i1)*(cc[cc_original_size+k][m].y-i1);
        double norm2=n_j0<n_j1?n_j0:n_j1;
        norm2=norm2<n_i0?norm2:n_i0;
        norm2=norm2<n_i1?norm2:n_i1; // WE TAKE INTO ACCOUNT THE DISTANCE TO THE BOUNDARY
        double norm=(c2[0]-cc[cc_original_size+k][m]).norm2()+10.*norm2;
        if( norm<min){
          min=norm;
          pos_1=m;
        }
      }
      // WE COMPUTE THE NEAREST NEW CONTOUR POINT TO OLD CONTOUR SECOND POINT
      int pos_2=0;
      int pos_3=c2.size()-1;
      min=1e30;
      if(cc[cc_original_size+k].size()>0){
        for(int m=0;m<cc[cc_original_size+k].size();m++){
          double n_j0=(cc[cc_original_size+k][m].x-j0)*(cc[cc_original_size+k][m].x-j0);
          double n_j1=(cc[cc_original_size+k][m].x-j1)*(cc[cc_original_size+k][m].x-j1);
          double n_i0=(cc[cc_original_size+k][m].y-i0)*(cc[cc_original_size+k][m].y-i0);
          double n_i1=(cc[cc_original_size+k][m].y-i1)*(cc[cc_original_size+k][m].y-i1);
          double norm2=n_j0<n_j1?n_j0:n_j1;
          norm2=norm2<n_i0?norm2:n_i0;
          norm2=norm2<n_i1?norm2:n_i1;  // WE TAKE INTO ACCOUNT THE DISTANCE TO THE BOUNDARY
          double norm=(c2[pos_3]-cc[cc_original_size+k][m]).norm2()+10.*norm2;
          if( norm<min){
            min=norm;
            pos_2=m;
          }
        }

        // WE ADD SEGMENTS TO CLOSE THE CONTOUR
        cc[cc_original_size+k].add_segment(c2[0],cc[cc_original_size+k][pos_1]);
        cc[cc_original_size+k].add_segment(c2[pos_3],cc[cc_original_size+k][pos_2]);

        // WE ADD THE REMAINING OLD CONTOUR POINTS
        for(int m=0;m<c2.size();m++) cc[cc_original_size+k].push_back(c2[m]);
      }
      else{
        cc.resize(cc_original_size+k);
        return;
      }

    }


    #ifdef SNAKE_DEBUG
   {
      static int cont_debug5=1000;
      ami::image<unsigned char> img_debug(width2,height2,3,255);
      int size2=width2*height2;
      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]=(unsigned char) g_I_img[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }
      std::stringstream ss;//create a stringstream
      ss << "snake_gauss_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss.str());

      for(int m=0;m!=size2;++m){
        img_debug[m]=img_debug[m+size2]=img_debug[m+2*size2]= snake_img[m];
      }
      for(int m=0;m<cc[cc_original_size+k].size();m++){
        int x=cc[cc_original_size+k][m].x-j0;
        int y=cc[cc_original_size+k][m].y-i0;
        if( x<width2 && y<height2){
          int n=y*width2+x;
          img_debug[n+size2]=255;
          img_debug[n]=0;
          img_debug[n+2*size2]=0;
        }
      }

      std::stringstream ss2;//create a stringstream
      ss2 << "snake_contours" << cont_debug5++ <<".jpg";
      img_debug.write(ss2.str());
    }
    #endif

  }

}







#endif
