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

/********************************************************
 *	Takes co-ordinates of points in from stdin	*
 *	out writes to stdout a pgm file of the		*
 *	fourier transform				*
 *	Based on fourier.c by Stuart Levy from 		*
 *	"Quasicrystal and Geometry" by Marjorie Senchal	*
 *	David Longbottom				*
 ********************************************************/



struct point {
	double x,y;
};

int getnumber (double *v) {
	int c;
	for(;;) {
		switch (c = getchar()) {
		case' ':case'\t':case'\n':case'\r':
		case',':case'(':case')':case'{':case'}':case'[':case']':
			continue;	//ignore whitespace etc
			
		case'#':		//commented-out line
			while((c = getchar()) != '\n' && c != EOF);
			continue;

		default:
			ungetc(c,stdin);
			return (scanf("%lf",v) == 1);
		}
	}
}

int main (int argc, char *argv[]) {
	double 	srange = 100.0;
	int 	n = 100;
	int	m = 100;
	double 	brightness = 1;
	double 	scale;
	int	i, j, k;
	int	i1, j1;
	struct point	s, *source;
	double	re, im;
	int	nsource, allocated;
	double	dot, mag, maxmag;

	if (argc <= 1 || (n = atoi(argv[1])) <= 0 ) {
		fprintf(stderr,"\
		Useage is:\n%s n [m] [srange] [brightness] < data points\n",argv[0]);
		return(EXIT_FAILURE);
	}

	if (n==1) n=2;
	m = n;
	if (argc > 2) m = atoi(argv[2]);
	if (argc > 3) srange = atof(argv[3]);
	if (argc > 4) brightness = atof(argv[4]);



	allocated = 15;
	source = (struct point *)malloc(allocated * sizeof(struct point));
	nsource = 0;
	while (getnumber(&source[nsource].x) && getnumber(&source[nsource].y)) {
		nsource++;
		if(nsource >= allocated) {
			allocated *=2;
			source = (struct point *)realloc(source,allocated*sizeof(struct point));
		}
	}
	fprintf(stderr,"%d light sources.\n",nsource);

	scale = nsource > 0 ? brightness * 256 / nsource : 1;

	printf("P5\n%d %d\n255\n",n,m);

	
	for(k=0; k<nsource;k++)	fprintf(stderr,"(%f,%f)\n",source[k].x,source[k].y);

	maxmag = 0;
	for(i1 = 0; i1<m; i1++) {
	i = (i1 + m/2) % m;		//hack to center the fourier transform	
	s.y = srange*((float)i/(n-1) - 0.5 );
	
	for(j1=0; j1<n;j1++) {
			j = (j1 + n/2) % n;	//hack to center the fourier transform	
			s.x = srange*((float)j/(n-1) - 0.5);

			re = im = 0;
			for(k=0;k < nsource;k++) {
				dot= 2*M_PI*(s.x*source[k].x + s.y*source[k].y);
				re += cos(dot);
				im += sin(dot);
			}

			mag = scale * sqrt(re*re + im*im);
			//mag = 255*log10(9*(re*re + im*im)/(255*255) + 1.0);
			putchar(mag >= 255 ? 255:(int)mag);
			if(maxmag < mag) maxmag = mag;
		}
	}

	if (maxmag > 0 )
		fprintf(stderr,"Max mag %.1f; max brightness which won't saturate:%g\n",maxmag,brightness*256/maxmag);

	free(source);

	return(EXIT_SUCCESS);
}