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);   }}