/*

  Programa para tirar tendencias de um conjunto amostral puntual

  Este programa é um filtro de média movel que calcula a média dos
  valores dos pontos que estejam dentro de uma região quadrada do
  espaco. O valor da aresta do quadrado é fornecido pelo usuário. É
  um valor em metros e deve ser fornecido quando o programa solicita
  "Entre com dx da janela:"

  A entrada é sempre um arquivo ASCII chamado in.dat onde temos os
  valores de x, y e z separados por espaço ou tabulação. O valor de
  z não é uma posicao espacial e sum um valor de atributo (temperatura,
  elevação, textura do solo, etc.... ).

  A saida é sempre um arquivo ASCII chamado out.dat. Esse arquivo
  contém várias informacoes (x, y, zf, z, zm, dp e npts) onde:
  x e y sao coordenadas espaciais
  zf é o valor final do atributo, sem a mediana local
  z é o valor original do atributo
  zm é o valor médio local do atributo dentro da janela de filtragem
  dp é o desvio padrão local da janela de filtragem
  npts é o número de pontos considerados dentro da janela de filtragem

*/




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

void
main()
{
char buf[300];
FILE *fpin, *fpout, *fpoutspr, *fpoutest, *fpoutdat;
float xmin, xmax, ymin, ymax, dx, dy;
float xp[1000],yp[1000],zp[1000];
double ztot, ztot2, media, desvp;
int i,j,npts,nptsjan;

printf("\nProgram to calculate the residual values");
printf("\nof a sample set based in a mean value ");
printf("\nof the samples in a window with size dx x dy.");
printf("\nThe input sample set have to be in a text file");
printf("\ncalled in.dat. The output residuals will be");
printf("\nstored in a text file called out.dat");
printf("\nThe program generates also two more files:");
printf("\nstatistic.dat with statistics of calculations, and");
printf("\nout.spr file ready to be imported in SPRING");

// Percorre arquivo, coloca dados no vetor e conta nro de pontos
if((fpin = fopen("in.dat", "rt"))==NULL){
	printf("\nFile in.dat not found");
	return;
}

printf("\n\nEnter the window size dx: ");
scanf("%f", &dx);
dy = dx;
printf("\nEnter the window size dy: ");
scanf("%f", &dy);

npts=0;

    while(fgets(buf,200,fpin)!= NULL){
    	sscanf(buf,"%f %f %f",&xp[npts], &yp[npts], &zp[npts]);
        npts++;
    }

if (fpin) fclose (fpin);

fpoutest=fopen("statistics.dat","wt");
fprintf(fpoutest,"Data in columns: xi yi zi-means zi mean std.dev npts_win");
fpoutdat=fopen("out.dat","wt");
fpoutspr=fopen("out.spr","wt");
fprintf(fpoutspr,"SAMPLE\nPOINT3D");

for (i=0; i<npts; i++){
	xmin = xp[i]-dx; xmax = xp[i]+dx;
	ymin = yp[i]-dy, ymax = yp[i]+dy;

	nptsjan=0;   ztot = ztot2 = (double)0.0;
	for (j=0; j<npts; j++)
		if(xp[j]<xmax && xp[j]>xmin && yp[j]<ymax && yp[j] > ymin){
			ztot += (double)zp[j];
			ztot2+= ((double)zp[j]*(double)zp[j]);
			nptsjan++;
		}
    if(nptsjan>0){
	   media= ztot/nptsjan;
    	desvp = ((double)nptsjan*ztot2-ztot*ztot)/(nptsjan*(nptsjan-1));
	   desvp = sqrt(desvp);
    	fprintf(fpoutest,"\n%10.4f %10.4f %10.4f %10.4f %f %f %d",xp[i],yp[i],zp[i]-media, zp[i],media,desvp,nptsjan);
    	fprintf(fpoutdat,"\n%10.4f %10.4f %10.4f",xp[i],yp[i],zp[i]-media);
    	fprintf(fpoutspr,"\n%10.4f %10.4f %10.4f",xp[i],yp[i],zp[i]-media);

     }

}

if (fpoutest) fclose(fpoutest);
if (fpoutdat) fclose(fpoutdat);
fprintf(fpoutspr,"\nEND");
fprintf(fpoutspr,"\nEND");
if (fpoutspr) fclose(fpoutspr);
printf("End of program");
}