#define MAX(X,Y) ((X)>(Y)?(X):(Y))#define MIN(X,Y) ((X)<(Y)?(X):(Y))#define ami_snake_tol_error 0.1 /* TOLERANCIA PARA EL ERROR DE CONVERGENCIA DE ami_snake() */#define ami_snake_dt 0.5        /* INCREMENTO TEMPORAL USADO PARA RESOLVER LA EDP EN ami_snake() */#define ami_snake_tol_grad2 10. /* TOLERANCIA PARA LA NORMA AL CUADRADRO DEL GRADIENTE DE LA IMAGEN                                   DE LEVEL SET PARA CAMBIAR EL VALOR DE LAIMAGEN. (ESTO ES SOLO PARA				   ACELERAR EL PROCESO Y EVITAR CALCULOS INNECESARIOS */void ami_semilla(long m,float *imagen,unsigned short width,long size, float nivel1,float nivel2);/* ALGORITMO DE SEMILLA PARA RELLENAR LA IMAGEN A PARTIR DE UN NIVEL CON OTRO NIVEL */int ami_rellenar_poligono(float **imagen,int *x,int *y,int Np,float nivel1,float nivel2,                          int borde,int *width,int *height,int *xd,int *yd);/* FUNCION QUE CONSTRUYE UNA IMAGEN RELLENANDO UN POLIGONO */int ami_snake(unsigned char *imagen,int width,int height,int *x,int *y,int Np,int borde,int *xd,int *yd,              float alfa,float lambda,float sigma,int Nescalas,float **imagen_level_set_r,int *width_g,int *height_g,int Niter);/* FUNCION QUE CALCULA LOS SNAKES GEODESICOS A PARTIR DE UN POLIGONO */void ami_snake_iteracion(float *imagen_level_set,float *imagen_level_set_x,float *imagen_level_set_y,                    float *imagen_g,float *imagen_g_x, float *imagen_g_y,int width2,int height2,                    float lambda,float *error,float max_grad_g);/* FUNCION QUE CALCULA UNA ITERACION EN EL PROCEDIMIENTO DE ami_snake() */void ami_dibujar_segmento(float *data,int width,int height,int x0,int y0,int x1,int y1,float color);/* FUNCION QUE DIBUJA UN SEGMENTO */void ami_gauss_conv_alma(float *rim,float sigma_x,float sigma_y,float precision,int width,int height);/* APPROXIMATION TO GAUSSIAN FILTER USING ALVAREZ-MAZORRA ALGORITHM */void ami_grad(float *image,int width,int height,float *image_x,float *image_y);/* FUNCION QUE CALCULA EL GRADIENTE DE UNA IMAGEN *//**************************************************//* FUNCION QUE CALCULA EL GRADIENTE DE UNA IMAGEN *//**************************************************/void ami_grad(float *image,/* ORIGINAL IMAGE */	      int width,int height,/* IMAGE SIZE */	      float *image_x,float *image_y /* OUTPUT GRADIENT */){  int i,j;  long m,size=width*height;  double coef1,coef2,c1,d1;  coef1=sqrt((double) 2.);  coef2=0.25*(2.-coef1);  coef1=0.5*(coef1-1);  /* BUCLE PRINCIPAL */  for(i=1;i<height-1;i++){    for(j=1;j<width-1;j++){      m=i*width+j;      c1=image[m+width+1]-image[m-width-1];      d1=image[m-width+1]-image[m+width-1];      image_y[m]=coef1*(image[m+width]-image[m-width])+coef2*(c1-d1);      image_x[m]=-(coef1*(image[m+1]-image[m-1])+coef2*(c1+d1));    }  }  /* EN LOS BORDES ASIGNAMOS 0 A LAS DERIVADAS */  for(m=0;m<width;m++){    image_y[m]=image_y[size-1-m]=0.;    image_x[m]=image_x[size-1-m]=0.;  }  for(m=0;m<(size-1);m+=width){    image_y[m]=image_y[m+width-1]=0.;    image_x[m]=image_x[m+width-1]=0.;  }}/********************************************************************//* APPROXIMATION TO GAUSSIAN FILTER USING ALVAREZ-MAZORRA ALGORITHM *//********************************************************************/void ami_gauss_conv_alma(float *rim,  /* Original picture and output picture*/float sigma_x,    /* standard deviation of the gaussian in the horizontal axis direction*/float sigma_y,    /* standard deviation of the gaussian in the vertical axis direction */float precision,  /* non-negative number to indicate the precisi髇 of the convolution */int width,int height) /* image dimensions */{   unsigned long m,M,cont;   unsigned short i,j,nx1,ny1;   double *z;   float l_x,v_x,l_y,v_y;   int nx=height;   int ny=width;   int Nc_x=(int) (precision*sigma_x>1?precision*sigma_x:1);   float t_x=sigma_x*sigma_x/(2*Nc_x);   int Nc_y=(int) (precision*sigma_y>1?precision*sigma_y:1);   float t_y=sigma_y*sigma_y/(2*Nc_y);   int Nc_a=(Nc_x>Nc_y)?Nc_x:Nc_y;   /* DEFINIMOS LOS PARAMETROS DEL ALGORITMO EN X E Y*/   if(t_x>0){     l_x=(1.+2.*t_x-sqrt((double) 4*t_x+1))/(2*t_x);     v_x=l_x/t_x;   }   if(t_y>0){     l_y=(1.+2.*t_y-sqrt((double) 4*t_y+1))/(2*t_y);     v_y=l_y/t_y;   }   nx1=nx-1; ny1=ny-1;   M=nx*ny;   /* COGEMOS MEMORIA PARA EL VECTOR DE PASO */   if(nx>ny)	z=(double *) malloc(nx*sizeof(double));   else 	z=(double *) malloc(ny*sizeof(double));   /* BUCLE PRINCIPAL DE LAS ITERACIONES */   for(cont=0;cont<Nc_a;cont++){     if(t_x>0 && cont<Nc_x){       for(j=0;j<ny;j++){         z[0]=rim[j]/(1-l_x);         m=j;         for(i=1;i<nx;i++){           z[i]=rim[m+=ny]+l_x*z[i-1];         }         rim[m]=z[nx1]/(1-l_x);         i=nx-1;         while(i>0){	       rim[m-ny]=(z[--i]+l_x*rim[m]);           m-=ny;         }       }	   m=0;       while(m<M) rim[m++]*=v_x;	 }	 if(t_y>0 && cont<Nc_y){       m=0;       for(i=0;i<nx;i++){         z[0]=rim[m]/(1-l_y);         for(j=1;j<ny;j++){           z[j]=rim[++m]+l_y*z[j-1];         }         rim[m]=z[ny1]/(1-l_y);         j=ny1;         while(j>0){           --m;	       rim[m]=z[--j]+l_y*rim[m+1];         }         m+=ny;       }       m=0;       while(m<M) rim[m++]*=v_y;     }  }  free(z);}/**********************************//* FUNCION QUE DIBUJA UN SEGMENTO *//*********************************/void ami_dibujar_segmento(   float *data, /* IMAGEN DONDE SE VA A DIBUJAR EL SEGMENTO */   int width, int height, /* DIMENSIONES DE LA IMAGEN */   int x0, int y0, /* COORDENADAS 1 PUNTO DEL SEGMENTO */   int x1, int y1,  /* COORDENADAS 2 PUNTO DEL SEGMENTO */   float color)     /* COLOR CON EL QUE SE RELLENA EL SEGMENTO */{ int incx = ami_abs(x1 - x0); int incy = ami_abs(y1 - y0); int p, dx, dy; int x, y; int temp; int nx=width; int ny=height; if (incx>incy) {   if (x0>x1){     temp=x0; x0=x1; x1=temp;     temp=y0; y0=y1; y1=temp;   }   x=x0; y=y0;   p = 2*incy - incx;   dy = (y1>y0)? 1 : -1;   data[y*nx+x] = color;   do {     if (p<0) {       p += 2*incy;     } else {       p += 2*incy - 2*incx;       y += dy;     }     x++;     data[y*nx+x] = color;   } while (x<x1); } else {     if (y0>y1) {       temp=x0; x0=x1; x1=temp;       temp=y0; y0=y1; y1=temp;     }     x=x0; y=y0;     p = 2*incx - incy;     dx = (x1>x0)? 1 : -1;     data[y*nx+x] = color;     do {       if (p<0) {         p += 2*incx;       } else {         p += 2*incx - 2*incy;	 x += dx;       }       y++;       data[y*nx+x] = color;    } while (y<y1);  }}/************************************************************************//* FUNCION QUE CALCULA UNA ITERACION EN EL PROCEDIMIENTO DE ami_snake() *//************************************************************************/void ami_snake_iteracion(  float *imagen_level_set, /* IMAGEN DEL CONJUNTO DE NIVEL */  float *imagen_level_set_x, /* DERIVADA EN X DE LA IMAGEN DEL CONJUNTO DE NIVEL */  float *imagen_level_set_y, /* DERIVADA EN y DE LA IMAGEN DE CONJUNTO DE NIVEL */  float *imagen_g,           /* IMAGEN g = 1/sqrt(1+alfa*grad(I)^2) */  float *imagen_g_x,         /* DERIVADA EN X DE LA IMAGEN g */  float *imagen_g_y,         /* DERIVADA EN Y DE LA IMAGEN g */  int width2,int height2,    /* DIMENSIONES IMAGEN QUE CONTIENE AL POLIGONO */  float lambda,              /* PARAMETRO DE BALANCE ENTRE TERMINOS DE REGULARIDAD Y ATRACCION */  float *error,              /* INCREMENTO MAXIMO DEL VALOR DE LA IMAGEN CONJUNTO DE NIVEL EN 1 ITER */  float max_grad_g)          /* MAXIMO GRADIENTE IMAGEN g PARA CALCULAR EL dt A UTILIZAR */{  int i,j;  long m,size2=width2*height2;  float error_max=0; /* VARIABLE ALMACENAR CAMBIO MAXIMO AL PASAR DE ITERACION */  double ux,uy,ux2,uy2,uxuy,l0,l1,l2,l3,l4,paso,p_escalar,norma_g,norma_u,uxg,uyg;  float *imagen_paso;  /* HACEMOS UNA COPIA DE LA IMAGEN imagen_level_set EN imagen_paso */  ami_malloc1d(imagen_paso,float,size2);  for(m=0;m<size2;m++) imagen_paso[m]=imagen_level_set[m];  /* AJUSTAMOS EL INCREMENTO TEMPORAL dt EN BASE A lambda Y max_grad_g */  double dt=ami_snake_dt/(1.+0.5*lambda*max_grad_g);  /* CALCULAMOS EL GRADIENTE DE LA SUPERFICIE DE NIVEL */  ami_grad(imagen_level_set,width2,height2,imagen_level_set_x,imagen_level_set_y);  /* CALCULAMOS EL OPERADOR DIFERENCIAL DE DIFUSION */  for(i=1;i<height2-1;i++){   for(j=1;j<width2-1;j++){     m=i*width2+j;     ux=imagen_level_set_x[m];     uy=imagen_level_set_y[m];     ux2=ux*ux;     uy2=uy*uy;     if((ux2+uy2)>ami_snake_tol_grad2){ /* SI EL GRADIENTE ES PEQUENO NO HACEMOS NADA */       uxuy=ux*uy;       if(ux2>uy2){	   if(uxuy>0){	     l0=uxuy-ux2-ux2-uy2;	     l1=ux2+ux2-uy2-uxuy;	     l1+=l1;	     l2=uy2-uxuy;	     l2+=l2;	     l3=uxuy+uxuy+uxuy+uy2;	     l4=uy2-uxuy;	   }	   else{	     l0=-ux2-ux2-uy2-uxuy;	     l1=ux2+ux2-uy2+uxuy;	     l1+=l1;	     l2=uy2+uxuy;	     l2+=l2;	     l3=uxuy+uy2;	     l4=uy2-uxuy-uxuy-uxuy;	   }       }       else{	   if(uxuy>0){	     l0=-uy2-uy2-ux2+uxuy;	     l1=ux2-uxuy;	     l1+=l1;	     l2=uy2+uy2-ux2-uxuy;	     l2+=l2;	     l3=uxuy+uxuy+uxuy+ux2;	     l4=ux2-uxuy;	   }	   else{	     l0=-uy2-uy2-ux2-uxuy;	     l1=ux2+uxuy;	     l1+=l1;