ami_snake.h

ami_snake算法源代码

	     l2=uy2+uy2-ux2+uxuy;
	     l2+=l2;
	     l3=uxuy+ux2;
	     l4=ux2-uxuy-uxuy-uxuy;
	   }
       }

       /* CALCULAMOS EL INCREMENTO PARA PASAR DE UNA ITERACION A LA SIGUIENTE */
       paso=l0*imagen_level_set[m] + 0.25*((imagen_level_set[m+width2]+imagen_level_set[m-width2])*l1 +
             (imagen_level_set[m+1]+imagen_level_set[m-1])*l2 + (imagen_level_set[m+width2+1]+imagen_level_set[m-width2-1])*l3 +
             (imagen_level_set[m+width2-1]+imagen_level_set[m-width2+1])*l4);
       paso/=(ux2+uy2);

      /* DISCRETIZAMOS EL TERMINO DE REACCION. TOMAMOS EL GRADIENTE DE u SIGUIENDO LA
         DIRECCION DEL FLUJO */
      if(imagen_g_x[m]<0) uxg=imagen_level_set[m+1]-imagen_level_set[m];
      else uxg=imagen_level_set[m]-imagen_level_set[m-1];
      if(imagen_g_y[m]>0) uyg=imagen_level_set[m+width2]-imagen_level_set[m];
      else uyg=imagen_level_set[m]-imagen_level_set[m-width2];

      /* ACTUALIZAMOS EL VALOR DEL PIXEL */
      paso=imagen_g[m]*paso-lambda*(imagen_g_x[m]*uxg-imagen_g_y[m]*uyg);
      imagen_paso[m]=imagen_level_set[m]+paso*dt;

      /* ACTUALIZAMOS EL ERROR */
      if(ami_abs(paso)>error_max) error_max=ami_abs(paso);
     }
   }
  }
  /* ASIGNAMOS EL ERROR ANTES DE SALIR */
  *error=error_max;

  /* ACTUALIZAMOS LA IMAGEN DEL CONJUNTO DE NIVEL */
  for(m=0;m<size2;m++) imagen_level_set[m]=imagen_paso[m];

  /* LIBERAMOS MEMORIA Y SALIMOS */
  free(imagen_paso);
}

/*********************************************************************/
/* FUNCION QUE CALCULA LOS SNAKES GEODESICOS A PARTIR DE UN POLIGONO */
/*********************************************************************/
int ami_snake(
  unsigned char *imagen, /* IMAGEN DONDE SE ENCUENTRA EL SNAKE */
  int width,int height,  /* DIMENSIONES DE LA IMAGEN */
  int *x,int *y,         /* COORDENADAS DEL POLIGONO */
  int Np,                /* NUMERO DE PUNTOS DEL POLIGONO */
  int borde,             /* BORDE DE SEPARACION DEL POLIGONO PARA CONSTRUIR LA IMAGEN DE CONJUNTO DE NIVEL */
  int *xd,int *yd,       /* DESPLAZAMIENTO DE LAS COORDENADAS DEL POLIGONO EN LA IMAGEN CONJUNTO DE NIVEL */
  float alfa,            /* PARAMETRO FUNCION g = 1/sqrt(1+alfa*grad(I)^2) */
  float lambda,          /* PARAMETRO DE BALANCE ENTRE TERMINOS DE REGULARIDAD Y ATRACCION */
  float sigma,           /* DESVIACION ESTANDARD DE LA GAUSSIANA PARA CALCULAR IMAGEN g */
  int Nescalas,          /* NUMERO DE ESCALAS DIFERENTES UTILIZADAS EN TODO EL PROCESO */
  float **imagen_level_set_r, /* IMAGEN DE CONJUNTO DE NIVEL RESULTADO (SU FRONTERA CORRESPONDE AL SNAKE) */
  int *width_g,int *height_g, /* DIMENSIONES IMAGEN CONJUNTO DE NIVEL */
  int Niter)                 /* NUMERO DE ITERACIONES MAXIMO EN CADA ESCALA */
{
  int i,j,k,width2,height2,xd2,yd2,iter;
  long m,size2;
  float *imagen_zoom,*imagen_level_set,*imagen_zoom_x,*imagen_zoom_y,*imagen_level_set_x,*imagen_level_set_y;
  float *imagen_g,*imagen_g_x,*imagen_g_y,**imagen_zoom2;
  float nivel1=200.; /* NIVEL EXTERIOR PARA LA IMAGEN DE LEVEL SET */
  float nivel2=50.; /* NIVEL INTERIOR PARA LA IMAGEN DE LEVEL SET */
  int cont1,cont2;
  float error,max_grad_g,sigma_0=sigma;

  if (Np<2){
    printf("ERROR funcion ami_snake: No puntos del poligono = %d muy bajo\n",Np);
    return(-1);
  }
  if (borde<2){
    printf("ERROR funcion ami_snake: Borde exterior de la imagen= %d muy bajo\n",borde);
    return(-1);
  }

  /* CREAMOS UNA NUEVA IMAGEN DONDE INCLUIR EL POLIGONO */
  if(ami_rellenar_poligono(&imagen_level_set,x,y,Np,nivel1,nivel2,borde,&width2,&height2,&xd2,&yd2)!=0)   return(-1);
  size2=width2*height2; *xd=xd2; *yd=yd2;
  /* COMPROBAMOS QUE LA IMAGEN DEL POLIGONO NO SE SALE DE LA IMAGEN ORIGINAL */
  if ((xd2+width2)>width || xd2<0 || (yd2+height2)>height || yd2<0 ){
    printf("ERROR funcion ami_snake: La subventana del poligono se sale de la imagen\n");
    return(-1);
  }

  /* COGEMOS MEMORIA PARA EL RESTO DE IMAGENES INVOLUCRADAS */
  ami_malloc1d(imagen_zoom,float,size2); ami_malloc1d(imagen_zoom_x,float,size2);
  ami_malloc1d(imagen_zoom_y,float,size2);
  ami_malloc1d(imagen_level_set_x,float,size2); ami_malloc1d(imagen_level_set_y,float,size2);
  ami_malloc1d(imagen_g,float,size2); ami_malloc1d(imagen_g_x,float,size2);
  ami_malloc1d(imagen_g_y,float,size2);
  ami_malloc2d(imagen_zoom2,float,width2,height2);

  /* RELLENAMOS LA IMAGEN imagen_zoom */
  for(i=0;i<height2;i++)
    for(j=0;j<width2;j++)
      imagen_zoom[i*width2+j]=imagen[(i+yd2)*width+xd2+j];

/* INICIO CODIGO XMW */
/*  level_set=CreateImage("level_set_image",RIM,width2,height2,"");
  zoom=CreateImage("snake_image",RIM,width2,height2,"");
  zoom_2=CreateImage("Gaussian_Image_convolution",RIM,width2,height2,"");
  g=CreateImage("image_g",RIM,width2,height2,"");


  for(m=0;m<size2;m++){
     level_set->RIMdata[m]=imagen_level_set[m];
     zoom->RIMdata[m]=imagen_zoom[m];
   }
   ShowImage(level_set);
   ShowImage(zoom); */
/* FIN CODIGO XMW */

  /* INICIAMOS EL PROCESO DE EVOLUCION DEL CONJUNTO DE NIVEL */
  for(k=0;k<Nescalas;k++){
    /* AJUSTAMOS DESVIACION TIPICA DE LA GAUSSIANA */
    sigma=sigma_0*(Nescalas-k);
    for(m=0;m<size2;m++) imagen_zoom2[0][m]=imagen_zoom[m];
/* INICIO CODIGO XMW */
 //  for(m=0;m<size2;m++){zoom_2->RIMdata[m]=imagen_zoom2[0][m];}
 //  ShowImage(zoom_2);
/* FIN CODIGO XMW */

    /* APLICAMOS LA CONVOLUCION CON UNA GAUSSIANA Y CALCULAMOS imagen_g*/
    ami_gauss_conv_alma(imagen_zoom2[0],sigma,sigma,10.,width2,height2);
    ami_grad(imagen_zoom2[0],width2,height2,imagen_zoom_x,imagen_zoom_y);
    for(m=0;m<size2;m++){
      double paso;
      paso=imagen_zoom_x[m]*imagen_zoom_x[m]+imagen_zoom_y[m]*imagen_zoom_y[m];
      imagen_g[m]=1./sqrt((double) 1. +alfa*paso);
    }
/* INICIO CODIGO XMW */
   /*for(m=0;m<size2;m++){g->RIMdata[m]=255*imagen_g[m];}
   ShowImage(g);
   for(m=0;m<size2;m++){zoom_2->RIMdata[m]=imagen_zoom2[0][m];}
   ShowImage(zoom_2); */
/* FIN CODIGO XMW */

    /* CALCULAMOS EL GRADIENTE DE g */
    ami_grad(imagen_g,width2,height2,imagen_g_x,imagen_g_y);

    /* CALCULAMOS EL MAXIMO DEL GRADIENTE DE g PARA LUEGO AJUSTAR el dt */
    for(m=0;m<size2;m++) imagen_level_set_x[m]=imagen_g_x[m]*imagen_g_x[m]+imagen_g_y[m]*imagen_g_y[m];
    ami_max_vector(imagen_level_set_x,m,max_grad_g,size2);
    max_grad_g=sqrt(max_grad_g);
    /* INICIALIZAMOS EL ERROR Y COMENZAMOS ITERACIONES */
    error=ami_snake_tol_error+1; iter=0;
    while(error>ami_snake_tol_error && Niter > iter){
      iter++;
      printf("iter=%d\n",iter);
      ami_snake_iteracion(imagen_level_set,imagen_level_set_x,imagen_level_set_y,
                          imagen_g,imagen_g_x,imagen_g_y,width2,height2,lambda,&error,max_grad_g);
/* INICIO CODIGO XMW */
 /*  for(m=0;m<size2;m++){level_set->RIMdata[m]=imagen_level_set[m];}
   ShowImage(level_set);
   for(m=0;m<size2;m++) zoom->RIMdata[m]=imagen_zoom[m];
   for(i=1;i<(height2-1);i++){
     for(j=1;j<(width2-1);j++){
       m=i*width2+j;
       if(imagen_level_set[m]<125.){
         cont1=cont2=0;
         if(imagen_level_set[m+1]>125.) cont1++;
	 else cont2++;
	 if(imagen_level_set[m-1]>125.) cont1++;
	 else cont2++;
	 if(imagen_level_set[m+width2]>125.) cont1++;
	 else cont2++;
	 if(imagen_level_set[m-width2]>125.) cont1++;
	 else cont2++;
	 if(cont1>0 && cont2>0) zoom->RIMdata[m]=255.;
      }
    }
  }
   ShowImage(zoom); */
/* FIN CODIGO XMW */

    }
  }
  /* AJUSTAMOS LA IMAGEN Y LAS DIMENSIONES DE LA IMAGEN DE CONJUNTO DE NIVEL */
  *width_g=width2; *height_g=height2;
  *imagen_level_set_r=imagen_level_set;
  /* LIBERAMOS MEMORIA */
  free(imagen_zoom); free(imagen_zoom_x); free(imagen_zoom_y); ami_free2d(imagen_zoom2);
  free(imagen_level_set_x); free(imagen_level_set_y);
  free(imagen_g); free(imagen_g_x); free(imagen_g_y);

  return(0);
}

/*******************************************************/
/* FUNCION CONSTRUYE UNA IMAGEN RELLENANDO UN POLIGONO */
/*******************************************************/
/* DEVUELVE -1 SI ALGO VA MAL y 0 SI TERMINA BIEN */
int ami_rellenar_poligono(
   float **imagen,  /* IMAGEN QUE SE VA A RELLENAR */
   int *x,int *y,   /* COORDENADAS POLIGONO */
   int Np,          /* NUMERO DE PUNTOS DEL POLIGONO */
   float nivel1,     /* NIVEL EXTERIOR CON EL QUE SE RELLENA */
   float nivel2,    /* NIVEL INTERIOR CON EL QUE SE RELLENA */
   int borde,       /* MARGEN QUE SE DEJA DE LA IMAGEN AL POLIGONO */
   int *width,int *height, /* DIMENSIONES DE LA IMAGEN CREADA */
   int *xd,int *yd)  /* DESPLAZAMIENTO REALIZADO DE LAS COORDENADAS DE LOS POLIGONOS */
{
  int i,j,xmin,xmax,ymin,ymax,width2,height2,xd2,yd2;
  long m,size;
  float *imagen2;

  if (Np<2){
    printf("ERROR funcion ami_rellenar_poligono: No puntos del poligono = %d muy bajo\n",Np);
    return(-1);
  }
  if (borde<2){
    printf("ERROR funcion ami_rellenar_poligono: Borde exterior de la imagen= %d muy bajo\n",borde);
    return(-1);
  }
  /* CALCULAMOS LOS MAXIMOS Y MINIMOS DE LAS COORDENADAS DE LOS POLIGONOS */
  ami_max_vector(x,i,xmax,Np); ami_max_vector(y,i,ymax,Np);
  ami_min_vector(x,i,xmin,Np); ami_min_vector(y,i,ymin,Np);

  /* DEFINIMOS DIMENSIONES DE LA IMAGEN A CREAR */
  *width=width2=xmax-xmin+2*borde+1; *height=height2=ymax-ymin+2*borde+1; size=width2*height2;

  /* DEFINIMOS DESPLAZAMIENTOS A REALIZAR EN LAS COORDENADAS DE LOS POLIGONOS */
  *yd=yd2=ymin-borde; *xd=xd2=xmin-borde;
  if(yd2<0 || xd2<0){
   printf("ERROR funcion ami_rellenar_poligono: Borde exterior de la imagen= %d muy alto\n",borde);
   printf("La imagen del poligono + el borde se salen de la imagen \n",borde);
   return(-1);
  }

  /* RESERVAMOS MEMORIA PARA LA IMAGEN Y LA RELLENAMOS CON EL NIVEL INTERIOR */
  ami_malloc1d(imagen2,float,size); *imagen=imagen2;
  for(m=0;m<size;m++) imagen2[m]=nivel1;

  /* DIBUJAMOS EL POLIGONO EN LA IMAGEN */
  for(i=0;i<(Np-1);i++){
    /* DIBUJAR SEGMENTO COLOR nivel1 (x[i]-xd2,y[i]-yd2) - (x[i+1]-xd2,y[i+1]-yd2) */
    ami_dibujar_segmento(imagen2,width2,height2,x[i]-xd2,y[i]-yd2,x[i+1]-xd2,y[i+1]-yd2,nivel2);
  }
  /* DIBUJAR SEGMENTO COLOR nivel1 (x[Np-1]-xd2,y[Np-1]-yd2) - (x[0]-xd2,y[0]-yd2) */
  ami_dibujar_segmento(imagen2,width2,height2,x[Np-1]-xd2,y[Np-1]-yd2,x[0]-xd2,y[0]-yd2,nivel2);

  /* RELLENAMOS EL EXTERIOR DEL POLIGONO CON EL NIVEL EXTERIOR nivel1 */
  ami_semilla((long) 0,imagen2,(unsigned short) width2,size,nivel1,nivel2);

  return(0);
}

/************************************************************************************/
/* ALGORITMO DE SEMILLA PARA RELLENAR LA IMAGEN A PARTIR DE UN NIVEL CON OTRO NIVEL */
/************************************************************************************/
void ami_semilla(
   long m /* POSICION DEL PUNTO ACTUAL */,
   float *imagen /* IMAGEN DE ENTRADA */ ,
   unsigned short width /* ANCHO DE LA IMAGEN*/ ,
   long size/* TAMANO DE LA IMAGEN */ ,
   float nivel1 /* NIVEL QUE SE VA MODIFICAR */,
   float nivel2 /* NIVEL CON EL QUE SE RELLENA*/)
{
   imagen[m]=nivel2;
   if((m-1)>=0){
   	if(imagen[m-1]==nivel1 )
   	    ami_semilla(m-1,imagen,width,size,nivel1,nivel2);
   }
   if((m+1)<size){
   	if(imagen[m+1]==nivel1)
   	    ami_semilla(m+1,imagen,width,size,nivel1,nivel2);
   }
   if((m-width)>=0){
   	if(imagen[m-width]==nivel1)
   	    ami_semilla(m-width,imagen,width,size,nivel1,nivel2);
   }
   if((m+width)<size){
   	if(imagen[m+width]==nivel1)
   	    ami_semilla(m+width,imagen,width,size,nivel1,nivel2);
   }
}