/* Copyright 2007 TeX Users Group

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301 USA.  */

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

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

#define FAR
#define BYTE unsigned char
#define SHORT short
#define WORD unsigned short
#define LONG long
#define DWORD unsigned long

/* _CRTIMP size_t __cdecl fread(void *, size_t, size_t, FILE *); */

#include <pshpack2.h>
typedef struct tagBITMAPFILEHEADER {
	WORD    bfType;
	DWORD   bfSize;
	WORD    bfReserved1;
	WORD    bfReserved2;
	DWORD   bfOffBits;
} BITMAPFILEHEADER, FAR *LPBITMAPFILEHEADER, *PBITMAPFILEHEADER;
#include <poppack.h>

typedef struct tagBITMAPINFOHEADER{
	DWORD      biSize;
	LONG       biWidth;
	LONG       biHeight;
	WORD       biPlanes;
	WORD       biBitCount;
	DWORD      biCompression;
	DWORD      biSizeImage;
	LONG       biXPelsPerMeter;
	LONG       biYPelsPerMeter;
	DWORD      biClrUsed;
	DWORD      biClrImportant;
} BITMAPINFOHEADER, FAR *LPBITMAPINFOHEADER, *PBITMAPINFOHEADER;

typedef struct tagRGBQUAD {
	BYTE    rgbBlue;
	BYTE    rgbGreen;
	BYTE    rgbRed;
	BYTE    rgbReserved;
} RGBQUAD;
typedef RGBQUAD FAR* LPRGBQUAD;

typedef struct tagBITMAPINFO {
	BITMAPINFOHEADER    bmiHeader;
	RGBQUAD             bmiColors[1];
} BITMAPINFO, FAR *LPBITMAPINFO, *PBITMAPINFO;

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

int verboseflag=1;
int showcolorflag=0;
int showhistoflag=0;

double pi=3.141592653;

double tpi=2*3.141592653;

BITMAPFILEHEADER bmfh;

BITMAPINFOHEADER bmih;

RGBQUAD colortable[256];

/* int *image=NULL; */
/* unsigned char *image=NULL; */
double *image=NULL;

double *cc=NULL;
double *ss=NULL;

double *cosih=NULL;
double *sineh=NULL;

double *cosiv=NULL;
double *sinev=NULL;

unsigned int histogram[256];

void showcolortable (RGBQUAD colortable[], int n) {
	int k;
	for (k = 0; k < n; k++) {
		printf("%d\t%d\t%d\t%d\n",
			   k, colortable[k].rgbRed, colortable[k].rgbGreen, colortable[k].rgbBlue);
	}
}

int readhead (FILE *input) {
/*	BITMAPFILEHEADER bmfh; */
/*	BITMAPINFOHEADER bmih; */
	
/*	read file header */
	fread(&bmfh, sizeof(bmfh), 1, input);
	if (bmfh.bfType != (77 * 256 + 66)) {			/*  "BM" */
		printf("Not BMP file %X\n", bmfh.bfType);
/*		return -1; */
	}
	else if (verboseflag) printf("BMP file %X\n", bmfh.bfType);
	if (verboseflag) printf("File size in words %lu\n", bmfh.bfSize);
	if (verboseflag) printf("Offset in bytes to image %lu\n", bmfh.bfOffBits);
/*	read bitmap info header */	
	fread(&bmih, sizeof(bmih), 1, input);
	if (verboseflag) printf("Size of header %lu\n", bmih.biSize);	
	if (verboseflag) printf("Width in pixels %ld\n", bmih.biWidth);	
	if (verboseflag) printf("Height in pixels %ld\n", bmih.biHeight);	
	if (verboseflag) printf("Number of image planes %u\n", bmih.biPlanes);
	if (verboseflag) printf("Bits per pixels %u\n", bmih.biBitCount);
	if (verboseflag) printf("Compression %lu\n", bmih.biCompression);
	if (verboseflag) printf("Size of compressed image %lu\n", bmih.biSizeImage);
	if (verboseflag) printf("Horizontal pixel per meter %ld\n", bmih.biXPelsPerMeter);
	if (verboseflag) printf("Vertical pixel per meter %ld\n", bmih.biYPelsPerMeter);
	if (verboseflag) printf("Number of colors used %lu\n", bmih.biClrUsed);
	if (verboseflag) printf("Number of `important' colors %lu\n", bmih.biClrImportant);
/*	read color table */
	fread(&colortable, sizeof(RGBQUAD), bmih.biClrUsed, input);
	if (showcolorflag) showcolortable(colortable, bmih.biClrUsed);
	return 0;
}

void makehistogram (void) {
	int i, j, k;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;

	for (k = 0; k < 256; k++) histogram[k] = 0;
	for (i = 0; i < h; i++) {
		for (j = 0; j < w; j++) {
/*			k = image [i * w + j]; */
			k = (int) image [i * w + j];
			histogram[k]++;
		}
	}
}

void showhistogram (void) {
	int k;
	for (k = 0; k < 256; k++) {
		if (histogram[k] > 0) printf("%d\t%d\n", k, histogram[k]);
	}
}

void trigtables (void) {
	int i, j;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;
	cosih = (double *) malloc (sizeof(double) * w);
	sineh = (double *) malloc (sizeof(double) * w);
	cosiv = (double *) malloc (sizeof(double) * h);
	sinev = (double *) malloc (sizeof(double) * h);
	for (j = 0; j < w; j++) {
		cosih[j] = cos (tpi * j / w);
		sineh[j] = sin (tpi * j / w);
	}
	for (i = 0; i < h; i++) {
		cosiv[i] = cos (tpi * i / h);
		sinev[i] = sin (tpi * i / h);
	}
}

/* double coscomp (double u, double v) { */
double coscomp (int id, int jd) {
	int i, j, k;
	int ii, jj;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;
/*	double theta; */
	double c;
	double sum = 0.0;

/*	u = tpi * jd / w; */
/*	v = tpi * id / w; */
	for (i = 0; i < h; i++) {
		k = i * w;
		ii = (i * id) % h;
/*		theta = u * i; */
		for (j = 0; j < w; j++) {
/*			k = i * w + j; */
/*			theta = u * i + v * j; */
/*			c = cos (theta); */
			jj = (j * jd) % w;
			c = cosih [ii] * cosiv [jj] - sineh [ii] * sinev [jj];
			sum += c * image [k];
			k++;
/*			theta += v; */
		}
/*		putc('.', stdout); */
	}
	return sum / (w * h);
}

/* double sincomp (double u, double v) { */
double sincomp (int id, int jd) {
	int i, j, k;
	int ii, jj;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;
/*	double theta; */
	double s;
	double sum = 0.0;

/*	u = tpi * jd / w; */
/*	v = tpi * id / w; */
	for (i = 0; i < h; i++) {
		k = i * w;
		ii = (i * id) % h;
/*		theta = u * i; */
		for (j = 0; j < w; j++) {
/*			k = i * w + j; */
/*			theta = u * i + v * j; */
/*			s = sin (theta); */
			jj = (j * jd) % w;
			s = sineh [ii] * cosiv [jj] + cosih [ii] * sinev [jj];
			sum += s * image [k];
			k++;
/*			theta += v; */
		}
/*		putc('!', stdout); */
	}
	return sum / (w * h);
}

void dotransform (void) {
	int i, j, k;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;
	double u, v;
	double c, s;

	cc = (double *) malloc (sizeof(double) * w * h);
	ss = (double *) malloc (sizeof(double) * w * h);

	for (i = 0; i < h; i++) {
		k = i * w;
		v = tpi * i / h;
		for (j = 0; j < w; j++) {
/*			k = i * w + j; */
			u = tpi * j / w;
/*			c = coscomp (u, v);*/
			c = coscomp (i, j);
/*			s = sincomp (u, v); */
			s = sincomp (i, j);
			cc [k] = c;
			ss [k] = s;
			printf("%d\t%d\t%lg\t%lg\t%lg\t%lg\n",  i, j, u, v, c, s);
			k++;
		}
	}
}

int readimage(FILE *input) {
	long i, j, k;
	unsigned char color;
	int w = (int) bmih.biWidth;
	int h = (int) bmih.biHeight;
	long n = w * h;
	int nbytes=1;

	if (n == 0) return -1;
/*	image = (unsigned char *) malloc(sizeof(unsigned char) * n); */
	image = (double *) malloc(sizeof(double) * n);
	if (image == NULL) {
		fprintf(stderr, "Unable to allocate %d pixel image\n", n);
		return -1;
	}
	fseek(input, sizeof(BITMAPFILEHEADER) + bmfh.bfOffBits, SEEK_SET);
	for (i = 0; i < h; i++) {
		k = i * w;
		for (j = 0; j < w; j++) {
			fread(&color, nbytes, 1, input);
/*			k = i * w + j; */
/*			image [i * w + j] = color; */
			image [k] = (double) color;			/* crude for now */
			k++;
		}
	}
	return 0;
}

int main (int argc, char *argv[]) {
	FILE *input;
	int firstarg=1;

	if (argc < firstarg + 1) exit(1);
	if ((input = fopen(argv[firstarg], "rb")) == NULL) {
		perror(argv[firstarg]);
		exit(1);
	}
	readhead(input);
	readimage(input);
	if (showhistoflag) {
		makehistogram();
		showhistogram();
	}
/*	return 0; */
	trigtables();
	dotransform();
	fclose (input);
	if (cc != NULL) free(cc);
	if (ss != NULL) free(ss);
	if (cosih != NULL) free (cosih);
	if (sineh != NULL) free (sineh);
	if (cosiv != NULL) free (cosiv);
	if (sinev != NULL) free (sinev);
	if (image != NULL) free (image);
	return 0;
}