/*
   RAW histogram creator v.0.93 (20070802) by Peter Ruevski
   Copyright 2006-2007 by Peter Ruevski (ruevs atsign hotmail point com)

   This is a command-line ANSI C program to create RAW histograms
   from raw photos converted to RAW PGM files by dcraw. Use it at your own risk.
*/

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

void UnknownOption(char *option)
{
	fprintf(stderr,"unknown option %s\n", option);
}

void Usage(char *programName)
{
	fprintf(stderr,"\n\
RAW histogram creator v.0.93 (20070802) by Peter Ruevski (ruevs\100hotmail\056com)\n\
\nUsage:\n\
%s [-?][-help] file\n\
\t-?, -help\tThis message\n\
\tfile\t\tPGM input file created with dcraw.\n\
\nThis program will build histograms of RAW files processed into 16 bit RAW PGM\n\
files by using dcraw (http://www.cybercom.net/~dcoffin/dcraw/). For example\n\
if you have a Canon RAW file called IMG_0005.CR2 first use dcraw:\n\
\n\tdcraw -D -4 -t 0 -o 0 -v IMG_0005.CR2\n\n\
this will create IMG_0005.pgm. Then run this program:\n\
\n\trawhistogram IMG_0005.pgm\n\n\
The result will be a (text) file called IMG_0005.csv that contains the RAW\n\
histogram data. This file can easily be opened in a spreadsheet program to plot\n\
the RGGB histograms.\n",programName);
}

/* returns the index of the first argument that is not an option; i.e.
   does not start with a dash or a slash
*/
int HandleOptions(int argc,char *argv[])
{
	int i,firstnonoption=0;

	for (i=1; i< argc;i++) {
		if (argv[i][0] == '/' || argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* An argument -? means help is requested */
				case '?':
					Usage(argv[0]);
					break;
				case 'h':
				case 'H':
					if (!stricmp(argv[i]+1,"help")) {
						Usage(argv[0]);
						break;
					}
					/* If the option -h means anything else
					 * in your application add code here
					 * Note: this falls through to the default
					 * to print an "unknow option" message
					*/
				/* add your option switches here */
				default:
					fprintf(stderr,"unknown option %s\n",argv[i]);
					break;
			}
		}
		else {
			firstnonoption = i;
			break;
		}
	}
	return firstnonoption;
}

typedef struct {
	unsigned long	width;
	unsigned long	height;
	unsigned long	bitspp;
	unsigned long	colors;
	void			*data;
} tRawData;

int FreeRaw(tRawData *r)
{
	if(!r)
		return 1;

	if(r->data) free(r->data);

	memset(r, 0, sizeof(*r));	// zero out the data pointer and other fields

	return 0;
}

int AllocRaw(tRawData *r, unsigned long width, unsigned long height, unsigned long bitspp, unsigned long colors)
{
	if(!r)
		return 1;

	r->data = malloc(width*height*colors*bitspp/8);

	if (!(r->data)) {
		FreeRaw(r);
		return 1;
	}

	r->width = width;
	r->height = height;
	r->bitspp = bitspp;
	r->colors = colors;

	return 0;
}

typedef struct {
	unsigned		size;							// Maximum value that the histogram can hold [0:size)
	unsigned long	*R, *G1, *G2, *B;
	unsigned		minR, minG1, minG2, minB;
	unsigned		maxR, maxG1, maxG2, maxB;
	unsigned long	sumR, sumG1, sumG2, sumB;
	double			meanR, meanG1, meanG2, meanB;	// Mean
	double			sdR, sdG1, sdG2, sdB;			// Standard deviation
} tHistogram;


int FreeHistogram(tHistogram *h)
{
	if(!h)
		return 1;

	if(h->R) free(h->R);
	if(h->G1) free(h->G1);
	if(h->G2) free(h->G2);
	if(h->B) free(h->B);

	memset(h, 0, sizeof(*h));	// zero out the pointers and size

	return 0;
}

int AllocHistogram(tHistogram *h, unsigned size, unsigned char clear)
{
	if(!h)
		return 1;

	if(clear) {
		memset(h, 0, sizeof(*h));	// zero out the pointers and other fields
		h->minR = h->minG1 = h->minG2 = h->minB = UINT_MAX;	// make min big so that the actual minimum is found
	}

	h->R = (unsigned long *) malloc(size*sizeof(h->R[0]));
	h->G1 = (unsigned long *) malloc(size*sizeof(h->G1[0]));
	h->G2 = (unsigned long *) malloc(size*sizeof(h->G2[0]));
	h->B = (unsigned long *) malloc(size*sizeof(h->B[0]));

	if (!( h->R && h->G1 && h->G2 && h->B)) {
		FreeHistogram(h);
		return 1;
	}

	if (clear) {
		memset(h->R, 0, size*sizeof(h->R[0]));
		memset(h->G1, 0, size*sizeof(h->G1[0]));
		memset(h->G2, 0, size*sizeof(h->G2[0]));
		memset(h->B, 0, size*sizeof(h->B[0]));
	}

	h->size = size;

	return 0;
}

int LoadRawPGM(FILE *fin, tRawData *r)
{
	fscanf(fin, "P%d\n%d %d\n%d\n", &(r->colors), &(r->width), &(r->height), &(r->bitspp));
	if( (5!=r->colors) || (65535!=r->bitspp) ) {
		fprintf(stderr,"This program needs a RAW (1 color) 16 bit PGM file\nUse \"dcraw -D -4 -t 0 -o 0 -v IMG_0005.CR2\"\n");
		return 1;
	}
	r->colors = 1;
	r->bitspp = 16;

	printf("Processing %d by %d image - %lu pixels.\n", r->width, r->height, (unsigned long)r->width*r->height);

	if(AllocRaw(r, r->width, r->height, r->bitspp, r->colors)) {
		fprintf(stderr,"Not enough memory!\n");
		return 1;
	}

	fread(r->data, r->width*r->height*r->colors*r->bitspp/8, 1, fin);

/*	printf("|");
	for(c=0; c<PROGRESSBAR_SIZE-2; ++c)
		printf("-");
	printf("|\n");
	progress=0;
	numvalues = width*height*colors;

	for (row=0; row < height; row++) {
		fread(ppm, 2*colors, width, fin);*/

	return 0;
}

int CreateRAWHistogram(const tRawData *r, tHistogram *h)
{
	const int HISTOGRAM_SIZE=16384;	//4096
	const int PROGRESSBAR_SIZE=70;

	int row, col, c;
	unsigned long numvalues, progress;

	h->size = HISTOGRAM_SIZE;
	if(AllocHistogram(h, h->size, 1)) {
		fprintf(stderr,"Not enough memory!\n");
		return 1;
	}

//	printf("Processing %d by %d image - %lu pixels.\n", r->width, r->height, (unsigned long)r->width*r->height);

	printf("|");
	for(c=0; c<PROGRESSBAR_SIZE-2; ++c)
		printf("-");
	printf("|\n");
	progress=0;
	numvalues = r->width * r->height * r->colors;

	for (row=0; row < r->height; row++) {
		unsigned long rowidx = row * r->width * r->colors;	/// * r->bitspp/8 gets "swallowed" by ((unsigned short*)(r->data)) below for 16bits only
		for (col=0; col < r->width; col++) {
			unsigned long colidx = col * r->colors;
			for (c=0; c < r->colors; c++) {
				unsigned short value = ntohs(((unsigned short*)(r->data))[rowidx+colidx+c]);
				if(value>=HISTOGRAM_SIZE) {
					printf("\nThere are values larger than %d in the file!", HISTOGRAM_SIZE-1);
					FreeHistogram(h);
					return 1;
				}

/*
   All RGB cameras use one of these Bayer grids:

	0x16161616:	0x61616161:	0x49494949:	0x94949494:

	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5
	0 B G B G B G	0 G R G R G R	0 G B G B G B	0 R G R G R G
	1 G R G R G R	1 B G B G B G	1 R G R G R G	1 G B G B G B
	2 B G B G B G	2 G R G R G R	2 G B G B G B	2 R G R G R G
	3 G R G R G R	3 B G B G B G	3 R G R G R G	3 G B G B G B

	CMOS Canons are 0x94949494
*/

				if(row&1)	//odd
					if(col&1) {	//odd
						h->B[value]++;
						if (h->minB>value) h->minB = value;
						if (h->maxB<value) h->maxB = value;
						++(h->sumB);
					} else {	//even
						h->G2[value]++;
						if (h->minG2>value) h->minG2 = value;
						if (h->maxG2<value) h->maxG2 = value;
						++(h->sumG2);
					}
				else		//even
					if(col&1) {	//odd
						h->G1[value]++;
						if (h->minG1>value) h->minG1 = value;
						if (h->maxG1<value) h->maxG1 = value;
						++(h->sumG1);
					} else {
						h->R[value]++;
						if (h->minR>value) h->minR = value;
						if (h->maxR<value) h->maxR = value;
						++(h->sumR);
					}

				// Progress "bar"
				if(PROGRESSBAR_SIZE*(row*r->width+col)*r->colors/numvalues >= progress) {
					progress++;
					printf(".");
				}
			}
		}
	}
	return 0;
}

int CalculateHistogramStatistics(tHistogram *h)
{
	unsigned i;

	for (i=h->minR; i<=h->maxR; ++i)
		h->meanR += h->R[i]*i;
	h->meanR /= h->sumR;

	for (i=h->minR; i<=h->maxR; ++i)
		h->sdR += pow(h->meanR - i, 2) * h->R[i];
	h->sdR = sqrt(h->sdR/h->sumR);


	for (i=h->minG1; i<=h->maxG1; ++i)
		h->meanG1 += h->G1[i]*i;
	h->meanG1 /= h->sumG1;

	for (i=h->minG1; i<=h->maxG1; ++i)
		h->sdG1 += pow(h->meanG1 - i, 2) * h->G1[i];
	h->sdG1 = sqrt(h->sdG1/h->sumG1);


	for (i=h->minG2; i<=h->maxG2; ++i)
		h->meanG2 += h->G2[i]*i;
	h->meanG2 /= h->sumG2;

	for (i=h->minG2; i<=h->maxG2; ++i)
		h->sdG2 += pow(h->meanG2 - i, 2) * h->G2[i];
	h->sdG2 = sqrt(h->sdG2/h->sumG2);


	for (i=h->minB; i<=h->maxB; ++i)
		h->meanB += h->B[i]*i;
	h->meanB /= h->sumB;

	for (i=h->minB; i<=h->maxB; ++i)
		h->sdB += pow(h->meanB - i, 2) * h->B[i];
	h->sdB = sqrt(h->sdB/h->sumB);

	return 0;
}

int SaveHistogramToFile(const tHistogram *h, FILE *fout)
{
// Tolerance for number of combing period matches that will make the histogram considered "bad"
#define NPER0TOLR	50
#define NPER0TOLG1	50
#define NPER0TOLG2	50
#define NPER0TOLB	50

// Threshholds for combing period length that will be counted
#define PER0TR_R	5
#define PER0TR_G1	5
#define PER0TR_G2	5
#define PER0TR_B	5

// Combing period length fluctuations that will be accepted
#define PER0FL_R	2
#define PER0FL_G1	2
#define PER0FL_G2	2
#define PER0FL_B	2

	unsigned i;

	long ind0R, ind0G1, ind0G2, ind0B;
	long per0R, per0G1, per0G2, per0B;
	long nper0R, nper0G1, nper0G2, nper0B;

	// Write the histograms to the output file
	ind0R = ind0G1 = ind0G2 = ind0B = 0;
	nper0R = nper0G1 = nper0G2 = nper0B = 0;
	per0R = per0G1 = per0G2 = per0B = 0;

	fprintf(fout, "%5s, %10s, %10s, %10s, %10s\r\n", "val", "R", "G1", "G2", "B");
	for(i=0; i<h->size; ++i) {
		if (!h->R[i]) {
			if ((i-ind0R)>1)
				if ( (abs(i-ind0R - per0R) <= PER0FL_R) && (PER0TR_R <= per0R) )
					nper0R++;
				per0R = i-ind0R;
			ind0R=i;
		}

		if (!h->G1[i]) {
			if ((i-ind0G1)>1)
				if ( (abs(i-ind0G1 - per0G1) <= PER0FL_G1) && (PER0TR_G1 <= per0G1) )
					nper0G1++;
				per0G1 = i-ind0G1;
			ind0G1=i;
		}

		if (!h->G2[i]) {
			if ((i-ind0G2)>1)
				if ( (abs(i-ind0G2 - per0G2) <= PER0FL_G1) && (PER0TR_G2 <= per0G2) )
					nper0G2++;
				per0G2 = i-ind0G2;
			ind0G2=i;
		}

		if (!h->B[i]) {
			if ((i-ind0B)>1)
				if ( (abs(i-ind0B - per0B) <= PER0FL_B) && (PER0TR_B <= per0B) )
					nper0B++;
				per0B = i-ind0B;
			ind0B=i;
		}

		fprintf(fout, "%5d, %10lu, %10lu, %10lu, %10lu\r\n", i, h->R[i], h->G1[i],	h->G2[i], h->B[i]);
	}

	// Print statistics at the end
	fprintf(fout, "%5s, %10lu, %10lu, %10lu, %10lu\r\n", "Sum", h->sumR, h->sumG1, h->sumG2, h->sumB);
	fprintf(fout, "%5s, %10u, %10u, %10u, %10u\r\n", "Min", h->minR, h->minG1, h->minG2, h->minB);
	fprintf(fout, "%5s, %10u, %10u, %10u, %10u\r\n", "Max", h->maxR, h->maxG1, h->maxG2, h->maxB);
	fprintf(fout, "%5s, %10f, %10f, %10f, %10f\r\n", "Mean", h->meanR, h->meanG1, h->meanG2, h->meanB);
	fprintf(fout, "%5s, %10f, %10f, %10f, %10f\r\n", "StDev", h->sdR, h->sdG1, h->sdG2, h->sdB);

	return ( (nper0R>NPER0TOLR) || (nper0G1>NPER0TOLG1) || (nper0G2>NPER0TOLG2) || (nper0B>NPER0TOLB) );
}


int main(int argc,char *argv[])
{
	int error;
	int i;
	FILE *fin, *fout;
	char outfname[1024];
	tRawData r;
	tHistogram h;

/*	if (argc == 1) {
		Usage(argv[0]);
		return 1;
	}*/

	/* handle the program options */
	if (! (i=HandleOptions(argc, argv)) ) {
		Usage(argv[0]);
		return 1;
	}

	if( !(fin = fopen( argv[i], "r+b" )))
	{
		printf( "Can not open input file: %s\r\n", argv[i] );
		return 1;
	};

	strncpy(outfname, argv[i], 1024);
	strncpy(strstr(outfname, ".")+1, "csv", 3);
	if( !(fout = fopen( outfname, "w+b" )))
	{
		if(fin)
			fclose(fin);
		printf( "Can not create output file: %s\r\n", outfname );
		return 1;
	};

	error = LoadRawPGM(fin, &r);
	error = CreateRAWHistogram(&r, &h);
	error = CalculateHistogramStatistics(&h);
	error = SaveHistogramToFile(&h, fout);

	if (error)
		printf("Combed\n");
	else
		printf("OK\n");

	FreeRaw(&r);
	FreeHistogram(&h);

	if(fin)
		fclose(fin);
	if(fout)
		fclose(fout);

	return error;
}