#include "main.h"

void filtre_AR(double *rcoeffs, double *icoeffs, double *rresult, double *iresult)
{
	double *rdptr, *idptr, *rrptr, *irptr;
	double rfiltre, ifiltre, norm;
	int l, m;
	
	rdptr = rcoeffs;
	idptr = icoeffs;
	rrptr = rresult;
	irptr = iresult;
	
	for (m=0; m<B; m++)
	{
		for (l=m; l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			rfiltre = 1/norm;
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-m], idptr[l-m], rrptr[l-m], irptr[l-m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += B-m;
		idptr += B-m;
		rrptr += B-m;
		irptr += B-m;
	}
	
	for (m=(B+1); m<(2*B); m++)
	{
		for (l=(2*B-m); l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			rfiltre = 1/norm;
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-2*B+m], idptr[l-2*B+m], rrptr[l-2*B+m], irptr[l-2*B+m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += m-B;
		idptr += m-B;
		rrptr += m-B;
		irptr += m-B;
	}
}


void filtre_PB(double *rcoeffs, double *icoeffs, double *rresult, double *iresult, int Ka, int Kb)
{
	double *rdptr, *idptr, *rrptr, *irptr;
	double rfiltre, ifiltre, norm;
	int l, m;
	
	rdptr = rcoeffs;
	idptr = icoeffs;
	rrptr = rresult;
	irptr = iresult;
	
	for (m=0; m<B; m++)
	{
		for (l=m; l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			if (l>Kb)
				rfiltre = 0;
			else if (l<Ka)
				rfiltre = 1/norm;
			else
				rfiltre = (1/norm)*((Kb-l)/(Kb-Ka));
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-m], idptr[l-m], rrptr[l-m], irptr[l-m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += B-m;
		idptr += B-m;
		rrptr += B-m;
		irptr += B-m;
	}
	
	for (m=(B+1); m<(2*B); m++)
	{
		for (l=(2*B-m); l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			if (sqrt(CARRE(l>Kb)))
				rfiltre = 0;
			else if (l<Ka)
				rfiltre = 1/norm;
			else
				rfiltre = (1/norm)*((Kb-l)/(Kb-Ka));
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-2*B+m], idptr[l-2*B+m], rrptr[l-2*B+m], irptr[l-2*B+m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += m-B;
		idptr += m-B;
		rrptr += m-B;
		irptr += m-B;
	}
}

void filtre_KH(double *rcoeffs, double *icoeffs, double *rresult, double *iresult, double H)
{
	double *rdptr, *idptr, *rrptr, *irptr;
	double rfiltre, ifiltre, norm;
	int l, m;
	
	rdptr = rcoeffs;
	idptr = icoeffs;
	rrptr = rresult;
	irptr = iresult;
	
	for (m=0; m<B; m++)
	{
		for (l=m; l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			if (l == 0)
				rfiltre = 1/norm;
			else 
				rfiltre = (1/norm)*(1/pow(l,H));
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-m], idptr[l-m], rrptr[l-m], irptr[l-m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += B-m;
		idptr += B-m;
		rrptr += B-m;
		irptr += B-m;
	}
	
	for (m=(B+1); m<(2*B); m++)
	{
		for (l=(2*B-m); l<B; l++)
		{
			norm = D_PI*sqrt(4*T_PI/(2*l+1));
			
			if (l == 0)
				rfiltre = 1/norm;
			else 
				rfiltre = (1/norm)*(1/pow(l,H));
			ifiltre = 0.0;
			
			compmult(rfiltre, ifiltre, rdptr[l-2*B+m], idptr[l-2*B+m], rrptr[l-2*B+m], irptr[l-2*B+m]);
			
			rrptr[l-m] *= norm;
			irptr[l-m] *= norm;
		}
		rdptr += m-B;
		idptr += m-B;
		rrptr += m-B;
		irptr += m-B;
	}
}

void convolution(int type)
{
	int cutoff, i, B2, Ka, Kb;
	double *tex, *itex, *rcoeffs, *icoeffs, *rresult, *iresult, *rresultat, *iresultat;
	double *workspace, *scratchpad, H;
	
	int rank, howmany_rank;
	fftw_iodim dims[1], howmany_dims[1];
	
	fftw_plan dctPlan, fftPlan;
	fftw_plan idctPlan, ifftPlan;
	double *weights;


	cutoff = (int)(B);
	B2 = CARRE(B);
	
	printf("Allocation memoire des differentes donnees necessaires\n");
	
	tex = (double *)malloc(4*B2 * sizeof(double));
	if (tex == NULL)
	{
		perror("Erreur allocation memoire tex\n");
		exit(1);
	}
	itex = (double *)malloc(4*B2 * sizeof(double));
	if (itex == NULL)
	{
		perror("Erreur allocation memoire itex\n");
		exit(1);
	}
	for (i=0; i<4*B2; i++)
	{
		tex[i] = liste_sommets[i].t;
		itex[i] = 0;
	}
	
	workspace = (double *)malloc((14*B2+24*B)*sizeof(double));
	if(workspace == NULL)
	{
		perror("Erreur allocation memoire workspace\n");
		exit(1);
	}
	
	rcoeffs = workspace;
	icoeffs = rcoeffs + B2;
	rresult = icoeffs + B2;
	iresult = rresult + B2;
	scratchpad = iresult + B2;
	
	weights = (double *)malloc((4*B)*sizeof(double));
	if (weights == NULL)
	{
		perror("Erreur allocation memoire weights\n");
		exit(1);
	}
	
	makeweights(B, weights);
	
	dctPlan = fftw_plan_r2r_1d(2*B, weights, tex, FFTW_REDFT10, FFTW_ESTIMATE);
	
	rank = 1;
	dims[0].n = 2*B;
	dims[0].is = 1;
	dims[0].os = 2*B;
	howmany_rank = 1;
	howmany_dims[0].n = 2*B;
	howmany_dims[0].is = 2*B;
	howmany_dims[0].os = 1;
	
	fftPlan = fftw_plan_guru_split_dft(rank, dims, howmany_rank, howmany_dims, tex, itex, workspace, workspace+(4*B2), FFTW_ESTIMATE);
	
	printf("Transformee en harmoniques spheriques directe\n");
	
	FST_semi_fly(tex, itex, rcoeffs, icoeffs, B, scratchpad, 1, cutoff, &dctPlan, &fftPlan, weights);
	
	fftw_destroy_plan(fftPlan);
	fftw_destroy_plan(dctPlan);
	free(tex);
	free(itex);
	
	
	if (type == AR)
	{
		printf("Aller-retour dans l'espace des harmoniques spheriques\n");
		filtre_AR(rcoeffs, icoeffs, rresult, iresult);
	}
	else if (type == PB)
	{	
		printf("Filtrage passe-bas\n");
		printf("Entrer les frequences min et max :	\n");
		scanf("%d %d", &Ka, &Kb);
		filtre_PB(rcoeffs, icoeffs, rresult, iresult, Ka, Kb);
	}
	else if (type == KH)
	{	
		printf("Filtrage en 1/k^H\n");
		printf("Entrer H :	\n");
		scanf("%lf", &H);
		filtre_KH(rcoeffs, icoeffs, rresult, iresult, H);
	}
			
	rresultat = (double *)malloc(4*B2 * sizeof(double));
	if (rresultat == NULL)
	{
		perror("Erreur allocation memoire rresultat\n");
		exit(1);
	}
	iresultat = (double *)malloc(4*B2 * sizeof(double));
	if (iresultat == NULL)
	{
		perror("Erreur allocation memoire iresultat\n");
		exit(1);
	}

	idctPlan = fftw_plan_r2r_1d(2*B, weights, rresultat, FFTW_REDFT01, FFTW_ESTIMATE);
	
	rank = 1;
	dims[0].n = 2*B;
	dims[0].is = 2*B;
	dims[0].os = 1;
	howmany_rank = 1;
	howmany_dims[0].n = 2*B;
	howmany_dims[0].is = 1;
	howmany_dims[0].os = 2*B;

	
	ifftPlan = fftw_plan_guru_split_dft(rank, dims, howmany_rank, howmany_dims, rresultat, iresultat, workspace, workspace+(4*B2), FFTW_ESTIMATE);
	
	printf("Transformee en harmoniques spheriques inverse\n");
	
	InvFST_semi_fly(rresult, iresult, rresultat, iresultat, B, scratchpad, 1, cutoff, &idctPlan, &ifftPlan);
	
	fftw_destroy_plan(ifftPlan);
	fftw_destroy_plan(idctPlan);
	free(weights);
	free(workspace);
	
	for (i=0; i<nb_sommets; i++)
		liste_sommets[i].t = rresultat[i];
	
	free(rresultat);
	free(iresultat);
	
}


void convolution_S2Kit(void)
{
	FILE *fp;
	char name[100];
	int i, maillage, taille;
	double *rsignal, *isignal, *rfiltre, *ifiltre, *rresultat, *iresultat;
	double *workspace;
	double x, y, z, t;
	
	//printf("Fichier filtre :\n");
	//printf("fichier de filtre : %s", name);
	
	rsignal = (double *) malloc(sizeof(double) * CARRE(2*B));
	isignal = (double *) malloc(sizeof(double) * CARRE(2*B));
	rfiltre = (double *) malloc(sizeof(double) * CARRE(2*B));
	ifiltre = (double *) malloc(sizeof(double) * CARRE(2*B));
	rresultat = (double *) malloc(sizeof(double) * CARRE(2*B));
	iresultat = (double *) malloc(sizeof(double) * CARRE(2*B));
	workspace = (double *) malloc(sizeof(double) * (14*CARRE(B)+26*B));
	
	printf("convolution S2Kit\n");
	
	fp = fopen("texture_28_09h54", "r");
	if (!fp)
	{
		fprintf(stderr, "Erreur d'ouverture du fichier\n");
		exit(1);
	}
	fscanf(fp, "%d\n", &maillage);
	fscanf(fp, "%d\n", &maillage);
	fscanf(fp, "%d\n", &taille);
	printf("test_0\n");
	for (i=0; i<CARRE(2*B); i++)
	{
		fscanf(fp, "%lf %lf %lf %lf\n", &x, &y, &z, &t);
		rfiltre[i] = t;
		ifiltre[i] = 0;
		rsignal[i] = liste_sommets[i].t;
		isignal[i] = 0;
		rresultat[i] = 0;
		iresultat[i] = 0;
	}
	fclose(fp);
	printf("test 1\n");
	
	Conv2Sphere_semi_fly(rsignal, isignal, rfiltre, ifiltre, rresultat, iresultat, B, workspace);
	printf("test 2\n");
	for (i=0; i<nb_sommets; i++)
		liste_sommets[i].t = rresultat[i];
	printf("test 3\n");
	free(rsignal);
	free(isignal);
	free(rfiltre);
	free(ifiltre);
	free(rresultat);
	free(iresultat);
	free(workspace);
	
}

void inv_filtre(void)
{
	int i, l, m;
	double *workspace, *rdata, *idata, *rresultat, *iresultat, *weights;
	double *rdptr, *idptr;
	
	int rank, howmany_rank;
	fftw_iodim dims[1], howmany_dims[1];
	
	fftw_plan idctPlan, ifftPlan;
	
	workspace = (double *)malloc((CARRE(B)*10+24*B)*sizeof(double));
	rdata = (double *)malloc(CARRE(B)*sizeof(double));
	idata = (double *)malloc(CARRE(B)*sizeof(double));
	rresultat = (double *)malloc(CARRE(2*B)*sizeof(double));
	iresultat = (double *)malloc(CARRE(2*B)*sizeof(double));
	weights = (double *)malloc((4*B)*sizeof(double));
	
	rdptr = rdata;
	idptr = idata;
	
	for (m=0; m<B; m++)
	{
		for (l=m; l<B; l++)
		{
			//if (m==0)
				rdptr[l-m] = 1;
			//else 
			//	rdptr[l-m] = 0;
			idptr[l-m] = 0;
		}
		rdptr += B-m;
		idptr += B-m;
	}
	
	for (m=(B+1); m<(2*B); m++)
	{
		for (l=(2*B-m); l<B; l++)
		{			
			//if (m==0)
				rdptr[l-m] = 1;
			//else 
			//	rdptr[l-m] = 0;
			idptr[l-m] = 0;
		}
		rdptr += m-B;
		idptr += m-B;
	}
	
	makeweights(B, weights);
	
	idctPlan = fftw_plan_r2r_1d(2*B, weights, rresultat, FFTW_REDFT01, FFTW_ESTIMATE);
	
	rank = 1;
	dims[0].n = 2*B;
	dims[0].is = 2*B;
	dims[0].os = 1;
	howmany_rank = 1;
	howmany_dims[0].n = 2*B;
	howmany_dims[0].is = 1;
	howmany_dims[0].os = 2*B;

	
	ifftPlan = fftw_plan_guru_split_dft(rank, dims, howmany_rank, howmany_dims, rresultat, iresultat, workspace, workspace+(4*CARRE(B)), FFTW_ESTIMATE);
	
	InvFST_semi_fly(rdata, idata, rresultat, iresultat, B, workspace, 0, B, &idctPlan, &ifftPlan);
	
	fftw_destroy_plan(ifftPlan);
	fftw_destroy_plan(idctPlan);
	free(weights);
	free(workspace);
	
	for (i=0; i<nb_sommets; i++)
		liste_sommets[i].t = rresultat[i];
	
	free(rresultat);
	free(iresultat);
	free(rdata);
	free(idata);

	
}
