/*******************************************************************/
/* 		WARNING! WARNING! WARNING!			*/
/* Le programme ci-apres permet de traiter seulement les images */
/* presentant une dynamique inferieure ou egale a 8 bits	*/
/****************************************************************/


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

#define MAX(a,b) (((a) >= (b)) ? (a) : (b))
#define ABS(a) (((a) >= 0) ? (a) : (-(a)))


/******************************************************************/
/*common procedure*/

void NormeGradient(unsigned char *in, unsigned char *out,
			  int xsize, int ysize, float ALPHA);

void Laplacien(unsigned char *in, float *out, int xsize, int ysize, float ALPHA);

void PassageParZero(float *input, unsigned char *output, int dim);

void ImageZero(float *in, unsigned char *out, int xsize, int ysize);

void Lissage_XY(float *pointeur_1, int xsize, int ysize, float alpha);

void gradient_X(float *pointeur_1, int xsize,int ysize, float alpha);

void gradient_Y(float *pointeur_1, int xsize, int ysize, float alpha);

void Compute2dGradientNorme(float *gx, float *gy, float *gn, int size);

void CoefGradDeri(float alpha, float *afd1, float *bfd1, float *bfd2);

void GradientDeri(float *input, float *output, float *work,
float a1, float b1, float b2, int dim);


/**********************************************************/
/*static proc.*/


static void SmoothDeri(float *input, float *output, float *work, float a1, float a2, float a11, float a22, float b1, float b2, int dim);

static void LaplacienDeri(float *input, float *output, float *work, float a1, float a2, float b1, float b2, int dim);

static void CoefSmthDeri(float alpha, float *afl1, float *afl2, float *aflr1, float *aflr2, float *bfl1, float *bfl2);

static void CoefLaplDeri(float alpha, float *a1, float *a2, float *b1,
float *b2);

/*------------------------------------------------------------------*/

/************************************************************/

void gradient_X(
    float *pointeur_1,
    int xsize,
    int ysize,
    float alpha
					)
	/* calcul du gradient dans la direction X */

{
    int i,j;
	 float *p1, *p2, *out, *in;
    float *work,*input,*output,afl1,afl2,aflr1,aflr2,bfl1,bfl2;

    if (xsize>ysize)
    {
    input = (float*) calloc(xsize,sizeof(float));
    output = (float*) calloc(xsize,sizeof(float));
	 work = (float*) calloc(xsize,sizeof(float));
    } 
    else
    {
    input = (float*) calloc(ysize,sizeof(float));
    output = (float*) calloc(ysize,sizeof(float));
	 work = (float*) calloc(ysize,sizeof(float));
    }


  /* gradient suivant les lignes  */
 
	CoefGradDeri(alpha,&afl1,&bfl1,&bfl2);

   p1 = pointeur_1;
   for (j=0;j<ysize;j++)
   { 
		
		GradientDeri(p1,output,work,afl1,bfl1,bfl2,xsize);
		out = output;
		for (i=0;i<xsize;i++) *p1++ = *out++;
		
   }

   /* lissage suivant les colonnes    */

	CoefSmthDeri(alpha,&afl1,&afl2,&aflr1,&aflr2,&bfl1,&bfl2);

	p1 = pointeur_1;
	for (i=0;i<xsize;i++)
	{
   	p2 = p1;
		in = input;
		for(j=0;j<ysize;j++)
		{
			*in++ = *p2;
			p2 += xsize;
		}
		SmoothDeri(input,output,work,afl1,afl2,aflr1,aflr2,bfl1,bfl2,ysize);
		out = output;
		p2 = p1;
		for (j=0;j<ysize;j++) 
		{
			*p2 = *out++;
			p2 += xsize;
		}
		p1++;
	}

	free(input); 
   free(output);
	free(work);
}
/******************************************************/

void Lissage_XY(
    float *pointeur_1,
    int xsize,
    int ysize,
    float alpha
					)
	/* Lissage dans les directions X et Y */

{
    int i,j;
	 float *p1, *p2, *out, *in;
    float *work,*input,*output,afl1,afl2,aflr1,aflr2,bfl1,bfl2;

    if (xsize>ysize)
    {
    input = (float*) calloc(xsize,sizeof(float));
    output = (float*) calloc(xsize,sizeof(float));
	 work = (float*) calloc(xsize,sizeof(float));
    } 
    else
    {
    input = (float*) calloc(ysize,sizeof(float));
    output = (float*) calloc(ysize,sizeof(float));
	 work = (float*) calloc(ysize,sizeof(float));
    }


  /* lissage suivant les lignes  */
 
	CoefSmthDeri(alpha,&afl1,&afl2,&aflr1,&aflr2,&bfl1,&bfl2);

   p1 = pointeur_1;
   for (j=0;j<ysize;j++)
   { 
		SmoothDeri(p1,output,work,afl1,afl2,aflr1,aflr2,bfl1,bfl2,xsize);
		out = output;
		for (i=0;i<xsize;i++) *p1++ = *out++;
		
   }

   /* lissage suivant les colonnes    */


	p1 = pointeur_1;
	for (i=0;i<xsize;i++)
	{
   	p2 = p1;
		in = input;
		for(j=0;j<ysize;j++)
		{
			*in++ = *p2;
			p2 += xsize;
		}
		SmoothDeri(input,output,work,afl1,afl2,aflr1,aflr2,bfl1,bfl2,ysize);
		out = output;
		p2 = p1;
		for (j=0;j<xsize;j++) 
		{
			*p2 = *out++;
			p2 += xsize;
		}
		p1++;
	}

	free(input); 
   free(output);
	free(work);
}


/******************************************************************/

void gradient_Y(
    float *pointeur_1,
    int xsize, int ysize,
	 float alpha
					)
/* calcul du gradient suivant les colonnes */

{
    int i,j;
    float *p1, *p2, *out, *in;
    float *work,*input,*output,afl1,afl2,aflr1,aflr2,bfl1,bfl2;

    if (xsize>ysize)
    {
    input = (float*) calloc(xsize,sizeof(float));
    output = (float*) calloc(xsize,sizeof(float));
    work = (float*) calloc(xsize,sizeof(float));
    } 
    else
    {
    input = (float*) calloc(ysize,sizeof(float));
    output = (float*) calloc(ysize,sizeof(float));
    work = (float*) calloc(ysize,sizeof(float));
    }



  /* gradient suivant les colonnes  */

    CoefGradDeri(alpha,&afl1,&bfl1,&bfl2);

    p1 = pointeur_1;
    for (i=0;i<xsize;i++)
    {
      p2 = p1;
      in = input;
      for(j=0;j<ysize;j++)
      {
            *in++ = *p2;
             p2 += xsize;
      }
            
      GradientDeri(input,output,work,afl1,bfl1,bfl2,ysize);
      out = output;
      p2 = p1;
      for (j=0;j<ysize;j++) 
      {
          *p2 = *out++;
           p2 += xsize;
      }
      p1++;
    }

    /* lissage suivant les lignes */

   CoefSmthDeri(alpha,&afl1,&afl2,&aflr1,&aflr2,&bfl1,&bfl2);

   p1 = pointeur_1;
   for (j=0;j<ysize;j++)
   { 
              
              SmoothDeri(p1,output,work,afl1,afl2,aflr1,aflr2,bfl1,bfl2,xsize); 
              out = output;
              for (i=0;i<xsize;i++) *p1++ = *out++;
              
   }

   free(input); 
   free(output);
   free(work);

}
/*****************************************************************************/
/*************************************************************/
/* filtre de deriche applique a une image 2D */
/* ALPHA: coef du filtre */

void NormeGradient(unsigned char *in, unsigned char *out,
			  int xsize, int ysize, float ALPHA)
{
   int i;
   float *gradx, *d1;
   float *grady, *d2;
   float *gradn, *gradptr;
/*	float maxgrad;*/

	if(!in || !out)
	{
		fprintf(stderr,"Probleme a l'entree de NormeGradient\n");
		exit(EXIT_FAILURE);
	}

/* lecture de l'image  */

   gradx = (float*) calloc(xsize*ysize,sizeof(float));
   grady = (float*) calloc(xsize*ysize,sizeof(float));
   gradn = (float*) calloc(xsize*ysize,sizeof(float));

	if(!gradx || !grady || !gradn)
	{
		fprintf(stderr,"Memory...\n");
		exit(EXIT_FAILURE);
	}

	d1 = gradx;
	d2 = grady;

	for(i=xsize*ysize; i>0; i--)
		*d1++ = *d2++  = (float)*in++;


/*   gradient en X   */

   gradient_X(gradx,xsize,ysize, ALPHA);

/*   gradient en Y   */

   gradient_Y(grady,xsize,ysize, ALPHA); 

/*   calcul de la norme   */

	Compute2dGradientNorme( gradx, grady, gradn, xsize*ysize);

	gradptr = gradn;
		
	for(i=0;i<xsize*ysize;i++) *out++ = (char)(*gradptr++);
	
	if(gradx) free(gradx);
	if(grady) free(grady);	
	if(gradn) free(gradn);
}


/*****************************************************************************/
/*   computation of the smoothing for a line "input" of length "dim"   */
/*   the result is stored in "output", "work is used for computation   */
/*   the coefficients are already computed                             */
/***********************************************************************/
static void SmoothDeri(
float *input, float *output, float *work,
 float a1, float a2, float a11, float a22, float b1, float b2,
int dim)
{
	 int n;

	*(output)=a1*(*(input));
	*(output+1)=(a1*(*(input+1))) + (a2*(*(input))) + (b1*(*(output)));
	for (n=2 ; n<dim ; n++)
		*(output+n)=(a1*(*(input+n))) + (a2*(*(input+n-1))) + (b1*(*(output+n-1))) + (b2*(*(output+n-2)));

	*(work+dim-1)=0.0;
	*(work+dim-2)=(a11*(*(input+dim-1))) ;
	for (n=dim-3 ; n>=0 ; n--)
		*(work+n)=(a11*(*(input+n+1))) + (a22*(*(input+n+2))) + (b1*(*(work+n+1))) + (b2*(*(work+n+2)));

	for (n=0 ; n<dim ; n++)
		*(output+n)=(*(work+n))+(*(output+n));
}

/*****************************************************************************/

/* filtre de deriche applique a une image 2D */
/* ALPHA: coef du filtre */

void Laplacien(unsigned char *in, float *out, int xsize, int ysize, float ALPHA)
{
    int i, j;
    float *inFloat, *outInter, *temp, *tempi, *tempo;
    float a1, b1, a2, b2;

    if(!in || !out)
    {
        fprintf(stderr,"Probleme a l'entree de NormeGradient\n");
        exit(EXIT_FAILURE);
    }

/*  lecture de l'image  */

    inFloat = (float*) calloc(xsize*ysize,sizeof(float));
    outInter = (float*) calloc(xsize*ysize,sizeof(float));
    temp = (float*) calloc(MAX(xsize,ysize),sizeof(float));
    tempi = (float*) calloc(MAX(xsize,ysize),sizeof(float));
    tempo = (float*) calloc(MAX(xsize,ysize),sizeof(float));

    if(!inFloat || !outInter || !temp || !tempi || !tempo)
    {
        fprintf(stderr,"Memory...\n");
        exit(EXIT_FAILURE);
    }

    for (j = 0; j < xsize*ysize; j++)
	inFloat[j] = (float)in[j];

    CoefLaplDeri(ALPHA, &a1, &a2, &b1, &b2);

/*  laplacien en X   */

    for (i = 0; i < ysize; i++) {    
        LaplacienDeri(&inFloat[i*xsize], &outInter[i*xsize], temp,
                      a1, a2, b1, b2, xsize);
    };

/*  laplacien en Y   */

    for (i = 0; i < xsize; i++) {
        for (j = 0; j < ysize; j++) tempi[j] = outInter[i + j*xsize];
        LaplacienDeri(tempi, tempo, temp,
                      a1, a2, b1, b2, ysize);
        for (j = 0; j < ysize; j++) out[i + j*xsize] = tempo[j];
    };

    free(inFloat);
    free(outInter);	
    free(temp);
    free(tempi);
    free(tempo);
}


/******************************************************************************/
static void LaplacienDeri(
float *input, float *output, float *work,
 float a1, float a2, float b1, float b2,
int dim)
{
	 int n;
	
	*(output)=(*(input));
	*(output+1)=(*(input+1))+(a1*(*(input)))+(b1*(*(output)));
	for (n=2;n<dim;n++)
		*(output+n)=(*(input+n))+(a1*(*(input+n-1)))+(b1*(*(output+n-1)))+(b2*(*(output+n-2)));

	*(work+dim-1)=0.0;
	*(work+dim-2)=(a2*(*(input+dim-1)));
	for (n=dim-3;n>=0;n--)
		*(work+n)=(a2*(*(input+n+1)))+(b2*(*(input+n+2)))+(b1*(*(work+n+1)))+(b2*(*(work+n+2)));

	for (n=0;n<dim;n++)
		*(output+n)= -((*(output+n))+(*(work+n))); 
}

/*****************************************************************************/
/*   computation of the coefficients of the second derivative   */
/****************************************************************/
static void CoefLaplDeri(
float alpha, float *a1, float *a2, float *b1, float *b2
)
{
	float k;

	k=(1.0-(exp((double)(-2.0*alpha))))/(2.0*alpha*(exp((double)(-alpha))));
	*a1= -(1.0+(k*alpha))*(exp((double)(-alpha)));
	*a2= (1.0-(k*alpha))*(exp((double)(-alpha)));
	*b1= (2.0*exp((double)(-alpha)));
	*b2= -(exp((double)(-2.0*alpha)));
}

/*******************************************************/
/*   computation of the coefficients of the smoothing   */
/********************************************************/
static void CoefSmthDeri(
float alpha, float *afl1, float *afl2, float *bfl1, float *bfl2,
float *aflr1, float *aflr2)
{
	*afl1= (((1-exp((double)(-alpha)))*(1-exp((double)(-alpha))))/(1+((2*alpha)*(exp((double)(-alpha))))-(exp((double)(-2*alpha)))));
	*afl2= (exp((double)(-alpha)))*(alpha-1)*(*afl1);

	*bfl1= (2*exp((double)(-alpha)));
	*bfl2= (-exp((double)(-2*alpha)));

	*aflr1= (exp((double)(-alpha))*(alpha+1))*(*afl1);
	*aflr2= (-exp((double)(-2*alpha)))*(*afl1);
}


/******************************************************************************/
/*   computation of the first derivative for a line "input" of length "dim"   */
/*   the result is stored in "output", "work is used for computation          */
/*   the coefficients are already computed                                    */
/******************************************************************************/

void GradientDeri(
float *input, float *output, float *work,
 float a1, float b1, float b2,
int dim)
{
	 int n;

	*(output)=0.0;
	*(output+1)=(a1*(*(input))) ;
	for (n=2 ; n<dim ; n++)
		*(output+n)=(a1*(*(input+n-1))) + (b1*(*(output+n-1))) + (b2*(*(output+n-2)));

	*(work+dim-1)=0.0;
	*(work+dim-2)=(-a1*(*(input+dim-1)));
	for (n=dim-3 ; n>=0 ; n--)
		*(work+n)=(-a1*(*(input+n+1))) + (b1*(*(work+n+1))) + (b2*(*(work+n+2)));

	for (n=0 ; n<dim ; n++)
		*(output+n) = -((*(work+n))+(*(output+n)));
}

/**************************************************************/
/*   computation of the coefficients of the first derivative   */
/***************************************************************/
void CoefGradDeri(
 float alpha, float *afd1, float *bfd1, float *bfd2)
{
	*afd1 =  (1-exp((double)(-alpha)));
	*afd1=(*afd1)*(*afd1);

	*bfd1= (2*exp((double)(-alpha)));
	*bfd2= (-exp((double)(-2*alpha)));
}
/****************************************************************/

void Compute2dGradientNorme( 
	 float *gx, float *gy,
	float *gn,
	int size)
	/* Calcul la norme du gradient pour une liste de taille size */
{
	 int i;
	 float temp1, temp2;
	 float *ptr;

	if (!gx || !gy || !gn)
	{
		fprintf(stderr,"Pas de data pour calculer la norme du gradient\n");
		exit(EXIT_FAILURE);
	}

	ptr = gn;
	for(i=0;i<size;i++)
	{
		temp1 = *gx++;
		temp2 = *gy++;
		temp1 = (float)sqrt((double)(temp1*temp1+temp2*temp2));
		*ptr++ = temp1;
	}
}

/***********************************************************************/
/*   Calcul des passages par 0                                         */
/***********************************************************************/

#define EPS 1.0e-4

void PassageParZero(float *input, unsigned char *output, int dim)
{
    int i;
    float prod;

    for (i = 0; i < dim-1; i++) {
        if ((ABS(input[i]) < EPS) && (ABS(input[i+1]) < EPS)) continue;

        prod = input[i]*input[i+1];
        if (prod <= 0) {
            if (ABS(input[i] < ABS(input[i+1]))) {
                output[i] = 1;
	    } else {
                output[i+1] = 1;
            };
        };
    };
}

/*****************************************************************************/

/* filtre de deriche applique a une image 2D */
/* ALPHA: coef du filtre */

void ImageZero(float *in, unsigned char *out, int xsize, int ysize)
{
    int i, j;
    float *tempf;
    unsigned char *tempi;

    if(!in || !out)
    {
        fprintf(stderr,"Probleme a l'entree de NormeGradient\n");
        exit(EXIT_FAILURE);
    }

/*  lecture de l'image  */

    tempf = (float *) calloc(MAX(xsize,ysize),sizeof(float));
    tempi = (unsigned char *) calloc(MAX(xsize,ysize),sizeof(char));

    if(!tempf || !tempi)
    {
        fprintf(stderr,"Memory...\n");
        exit(EXIT_FAILURE);
    }

    for (i = 0; i < xsize*ysize; i++) out[i] = 0;

/*  zeros en X   */

    for (i = 0; i < ysize; i++) {
        PassageParZero(&in[i*xsize], &out[i*xsize], xsize);
    };

/*  zeros en Y   */

    for (i = 0; i < xsize; i++) {
        for (j = 0; j < ysize; j++) tempf[j] = in[i + j*xsize];
        for (j = 0; j < ysize; j++) tempi[j] = 0;
        PassageParZero(tempf, tempi, ysize);
        for (j = 0; j < ysize; j++) out[i + j*xsize] = tempi[j];
    };

    free(tempf);
    free(tempi);
}