ami_snake.h

ami_snake算法源代码

#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;