#include <stdio.h>
#include <math.h>

void traitement_image(float**& contour, float** image, int dimx, int dimy, float factor);
void lecture_image(float**& image, int& dimx, int& dimy);
void synthese_image(float**& image, int& dimx, int& dimy);
void ecriture_image(float** contour, int dimx, int dimy);
void ecriture_imageo(float** contour, int dimx, int dimy);
void affichage_image(float** image, int dimx, int dimy);
void lissage_image(float** image, int dimx, int dimy);

float a = -0.5;
float EPS = 0.08;
#define ABS(a) (((a) >= 0) ? (a) : (-(a)))

void main(int argc, char **argv)
{
  int dimx, dimy;
  float** image;
  float** contour;
  float factor = 1;
  
  printf("\nDetection de contours par splines\n");
  printf("Lecture de l'image\n");
  lecture_image(image,dimx,dimy);
  //synthese_image(image,dimx,dimy);

  lissage_image(image,dimx,dimy);
  ecriture_imageo(image,dimx,dimy);
  
  printf("Traitement de l'image\n");
  traitement_image(contour,image,dimx,dimy,factor);
  
  printf("Ecriture du resultat\n");
  ecriture_image(contour,int(dimx*factor),int(dimy*factor));  

}

void lecture_image(float**& image, int& dimx, int& dimy)
{
  FILE* dimFile;
  FILE* imaFile;
  
  dimFile = fopen("load.dim","r");
  fscanf(dimFile,"%d %d", &dimx, &dimy);
  fclose(dimFile);

  imaFile = fopen("load.ima","r");
  unsigned char** buffer = new unsigned char*[dimy];
  
  image = new float*[dimx];
  for (int i = 0; i<dimx; i++)
    image[i] = new float[dimy];
  
  

  float max = 0;
  float min = 256;
  for (int i = 0; i<dimy; i++)
    {
      buffer[i] = new unsigned char[dimx];
      fread(buffer[i], sizeof(char), dimx, imaFile);
    }
  
  for (int i = 0; i<dimy; i++)
    for (int j = 0; j<dimx; j++)
      if (buffer[i][j] > max)
	max = buffer[i][j];
      else if (buffer[i][j] < min)
	min = buffer[i][j];
  
  for (int k = 0; k<dimy; k++)
    {
      for (int l = 0; l<dimx; l++)
	{	  
	  image[l][k] = (buffer[k][l] - min)/(max-min);	  
	}
      printf("Ligne %d\r",k);
    }
  printf("\n");
  fclose(imaFile);
  
  for (int k = 0; k<dimy; k++)
    delete buffer[k];
  delete buffer;
}

void lissage_image(float** image, int dimx, int dimy)
{
  float** temp = new float*[dimx];

  for (int i = 1; i<dimx-1; i++)
    {
      temp[i] = new float[dimy];
      for (int j = 1; j<dimy-1; j++)
	{
	  temp[i][j] = (2*image[i][j] + image[i][j+1] + image[i][j-1]+
			image[i-1][j] + image[i+1][j] )/6;	  
	}
    }
  for (int i = 1; i<dimx-1; i++)
    for (int j = 1; j<dimy-1; j++)
      image[i][j] = temp[i][j];
}


void synthese_image(float**& image, int& dimx, int& dimy)
{
  dimx = 100;
  dimy = 100;
  image = new float*[dimx];
  for (int i = 0; i<dimx; i++)
    image[i] = new float[dimy];
  
  for (int x = 0; x<dimx; x++)
    for (int y = 0; y<dimy; y++)
      {
	image[x][y] = (exp(-0.01*((x-0.5*dimx)*(x-0.5*dimx) +
				    (y-0.5*dimy)*(y-0.5*dimy))));
	
      }
}

float c0(float t)
{
  return -a*t*t*t+a*t*t;
}

float c1(float t)
{
  return -(a+2)*t*t*t+(2*a+3)*t*t-a*t;
}

float c2(float t)
{
  return (a+2)*t*t*t-(a+3)*t*t+1;
}

float c3(float t)
{
  return a*t*t*t-2*a*t*t+a*t;
}

float c_0(float t)
{
  return -3*a*t*t+2*a*t;
}

float c_1(float t)
{
  return -3*(a+2)*t*t+2*(2*a+3)*t-a;
}

float c_2(float t)
{
  return 3*(a+2)*t*t-2*(a+3)*t;
}

float c_3(float t)
{
  return 3*a*t*t-4*a*t+a;
}

float c__0(float t)
{
  return -6*a*t+2*a;
}

float c__1(float t)
{
  return -6*(a+2)*t+2*(2*a+3);
}

float c__2(float t)
{
  return 6*(a+2)*t-2*(a+3);
}

float c__3(float t)
{
  return 6*a*t-4*a;
}


float H(float** image, int dimx, int dimy, int j, float x)
{
  int i = floor(x);
  float result = image[i][j] * c2(x-i);
  if (i>0) result += image[i-1][j] * c3(x-i);
  if (i<dimx-1)
    {
      result += image[i+1][j] * c1(x-i);
      if (i<dimx-2) result += image[i+2][j] * c0(x-i);
    }
  return result;
}

float H_(float** image, int dimx, int dimy, int j, float x)
{
  int i = floor(x);
  float result = image[i][j] * c_2(x-i);
  if (i>0) result += image[i-1][j] * c_3(x-i);
  if (i<dimx-1)
    {
      result += image[i+1][j] * c_1(x-i);
      if (i<dimx-2) result += image[i+2][j] * c_0(x-i);
    }
  return result;
}

float H__(float** image, int dimx, int dimy, int j, float x)
{
  int i = floor(x);
  float result = image[i][j] * c__2(x-i);
  if (i>0) result += image[i-1][j] * c__3(x-i);
  if (i<dimx-1)
    {
      result += image[i+1][j] * c__1(x-i);
      if (i<dimx-2) result += image[i+2][j] * c__0(x-i);
    }
  return result;
}


float f(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H(image,dimx,dimy,j,x) * c2(y-j);
  if (j>0) result += H(image,dimx,dimy,j-1,x) * c3(y-j);
  if (j<dimy-1)
    {
      result += H(image,dimx,dimy,j+1,x) * c1(y-j);
      if (j<dimy-2) result += H(image,dimx,dimy,j+2,x) * c0(y-j);
    }
  return result;
}

float f_x(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H_(image,dimx,dimy,j,x) * c2(y-j);
  if (j>0) result += H_(image,dimx,dimy,j-1,x) * c3(y-j);
  if (j<dimy-1)
    {
      result += H_(image,dimx,dimy,j+1,x) * c1(y-j);
      if (j<dimy-2) result += H_(image,dimx,dimy,j+2,x) * c0(y-j);
    }
  return result;
}

float f_y(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H(image,dimx,dimy,j,x) * c_2(y-j);
  if (j>0) result += H(image,dimx,dimy,j-1,x) * c_3(y-j);
  if (j<dimy-1)
    {
      result += H(image,dimx,dimy,j+1,x) * c_1(y-j);
      if (j<dimy-2) result += H(image,dimx,dimy,j+2,x) * c_0(y-j);
    }
  return result;
}

float f__x(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H__(image,dimx,dimy,j,x) * c2(y-j);
  if (j>0) result += H__(image,dimx,dimy,j-1,x) * c3(y-j);
  if (j<dimy-1)
    {
      result += H__(image,dimx,dimy,j+1,x) * c1(y-j);
      if (j<dimy-2) result += H__(image,dimx,dimy,j+2,x) * c0(y-j);
    }
  return result;
}

float f__y(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H(image,dimx,dimy,j,x) * c__2(y-j);
  if (j>0) result += H(image,dimx,dimy,j-1,x) * c__3(y-j);
  if (j<dimy-1)
    {
      result += H(image,dimx,dimy,j+1,x) * c__1(y-j);
      if (j<dimy-2) result += H(image,dimx,dimy,j+2,x) * c__0(y-j);
    }
  return result;
}

float f__xy(float x, float y, float** image, int dimx, int dimy)
{
  int j = floor(y);
  float result = H_(image,dimx,dimy,j,x) * c_2(y-j);
  if (j>0) result += H_(image,dimx,dimy,j-1,x) * c_3(y-j);
  if (j<dimy-1)
    {
      result += H_(image,dimx,dimy,j+1,x) * c_1(y-j);
      if (j<dimy-2) result += H_(image,dimx,dimy,j+2,x) * c_0(y-j);
    }
  return result;
}


void traitement_image(float**& contour, float** image, int dimx, int dimy, float factor)
{
  int sdimx = int(factor * dimx);
  int sdimy = int(factor * dimy);

  contour = new float*[sdimx];
  for (int j = 0; j<sdimx; j++)
    contour[j] = new float[sdimy];
  
  for (int i = 0; i<sdimx; i++)
    {
      float I = i/factor;
      for (int j = 0; j<sdimy; j++)
	{
	  float J = j/factor;
	  float fx = f_x(I,J,image,dimx,dimy);
	  float fy = f_y(I,J,image,dimx,dimy);
	  if (fx*fx+fy*fy != 0)
	    contour[i][j] = (f__x(I,J,image,dimx,dimy)*fx*fx+
			     f__y(I,J,image,dimx,dimy)*fy*fy+
			     2*f__xy(I,J,image,dimx,dimy)*fx*fy)/(fx*fx+fy*fy);
	  else
	    contour[i][j] = 0;
	 
	}
      printf("Ligne %d\r",i);
      fflush(stdout);
    }
  printf("\n");
}

void ecriture_imageo(float** contour, int dimx, int dimy)
{
  FILE* dimFile;
  FILE* imaFile;

  dimFile = fopen("saveo.dim","w");
  fprintf(dimFile,"%d %d", dimx, dimy);
  fclose(dimFile);

  
  imaFile = fopen("saveo.ima","w");
  unsigned char* buffer = new unsigned char[dimx];
  
  for (int i = 0; i<dimy; i++)
    {
      for (int j = 0; j<dimx; j++)
	{
	  buffer[j] = (char)(255*contour[j][i]);
	  
	}
      fwrite(buffer, sizeof(char), dimx, imaFile);
      
    }
  
  fclose(imaFile);
  
  delete buffer;
}

void ecriture_image(float** contour, int dimx, int dimy)
{
  FILE* dimFile;
  FILE* imaFile;

  dimFile = fopen("save.dim","w");
  fprintf(dimFile,"%d %d", dimx, dimy);
  fclose(dimFile);

  imaFile = fopen("save.ima","w");
  unsigned char* buffer = new unsigned char[dimx];
  float max = 0;
  float min = 0;
  for (int i = 0; i<dimy; i++)
    for (int j = 0; j<dimx; j++)
      if (contour[j][i] > max)
	max = contour[j][i];
      else if (contour[j][i] < min)
	min = contour[j][i];
  for (int i = 0; i<dimy; i++)
    {
      for (int j = 0; j<dimx; j++)
	{
	  if (j>0 && i>0 && j<dimx-1 && i<dimy-1)
	    if ( ( (ABS(contour[j-1][i]) > EPS || ABS(contour[j][i]) > EPS)
		   && contour[j-1][i] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j-1][i]) )
		 || ( (ABS(contour[j+1][i]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j+1][i] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j+1][i]) )
		 || ( (ABS(contour[j][i-1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j][i-1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j][i-1]) )
		 || ( (ABS(contour[j][i+1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j][i+1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j][i+1]) ) 
		 || ( (ABS(contour[j+1][i+1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j+1][i+1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j+1][i+1]) )
		 || ( (ABS(contour[j+1][i-1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j+1][i-1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j+1][i-1]) )
		 || ( (ABS(contour[j-1][i-1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j-1][i-1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j-1][i-1]) )
		 || ( (ABS(contour[j-1][i+1]) > EPS || ABS(contour[j][i]) > EPS)
		      && contour[j-1][i+1] * contour[j][i] < 0 && 
		   ABS(contour[j][i]) > ABS(contour[j-1][i+1]) ) ) 
	      buffer[j] = 1;
	    else 
	      buffer[j] = 0;
	  else
	    buffer[j] = 0;

	  
	  
	  
	}
      fwrite(buffer, sizeof(char), dimx, imaFile);
      
    }
  
  fclose(imaFile);
  
  delete buffer;
}


void affichage_image(float** image, int dimx, int dimy)
{
  for (int i = 0; i<dimx; i++)
    {
      for (int j = 0; j<dimy; j++)
	printf("%.3f\t",image[i][j]);
      printf("\n");
    }
  printf("\n");
}