/* to compile this unit, use the following as a Makefile:
----------------------------------------------
prefix     = /usr/local
includedir = $(prefix)/include
libdir     = $(prefix)/lib

CFLAGS=-O3

all: extractstat

extractstat: extractstat.c
	gcc $(CFLAGS) -I$(includedir) -L$(libdir) -lmba -lm -Wall -o extractstat `Wand-config --cflags --cppflags` extractstat.c `Wand-config --ldflags --libs`
----------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/stat.h>
#include <errno.h>
#include <math.h>
#include <mba/csv.h>
#include <wand/magick-wand.h>

#define BUF_MAX 0x7FFF
#define ROW_MAX 10

#define ROW_FRAME 0
#define ROW_QUANT 1
#define ROW_PSNR 6

#define ThrowWandException(wand) \
{ \
	char *description; \
	ExceptionType severity; \
\
	description = MagickGetException(wand,&severity); \
	(void) fprintf(stderr,"%s %s %ld %s\n",GetMagickModule(),description); \
	description = (char *) MagickRelinquishMemory(description); \
	exit(-1); \
}

int main(int argc, char **argv)
{
	struct stat statbuf;
	char *filename = NULL;
	FILE *in;
	unsigned char buf[BUF_MAX];
	unsigned char *row[ROW_MAX];
	char pictfilename[512];
	int n;
	unsigned long i, q;

	MagickBooleanType status;
	MagickWand *magick_wand, *reference_wand;
	double distortion;
	long double *allmae;

	unsigned long count = 0, pictcount = 0;
	long double min_psnr = 1000.0, mean_psnr = 0.0, max_psnr = 0.0, sd_psnr = 0.0, psnr;
	long double min_mae = 100000.0, mean_mae = 0.0, max_mae = 0.0, sd_mae = 0.0, mae;
	long double min_quant = 1000.0, mean_quant = 0.0, max_quant = 0.0, sd_quant = 0.0, quant;
	unsigned long q_freq[25];
	for ( q=0; q<25; q++ ) q_freq[q] = 0;

	if ((argc != 4) && (argc != 2))
	{
		fprintf(stderr, "usage: %s <filename> [<refpictpath> <comppictpath>]\n", argv[0]);
		return EXIT_FAILURE;
	}

	filename = argv[1];

	in = fopen(filename, "r");
	if (in == NULL)
	{
		perror(filename);
		return EXIT_FAILURE;
	}

	while ((n = csv_row_fread(in, buf, BUF_MAX, row, ROW_MAX, ',', CSV_TRIM | CSV_QUOTES)) > 0)
	{
		psnr = (long double)atof(row[ROW_PSNR]);
		if (psnr<1.0 || psnr>100.0) continue;
		quant = (long double)atof(row[ROW_QUANT]);
		if (quant<1.0 || quant>24.0) continue;
		if (psnr > max_psnr) max_psnr = psnr;
		if (psnr < min_psnr) min_psnr = psnr;
		mean_psnr += psnr;
		q_freq[q-1]++;
		if (quant > max_quant) max_quant = quant;
		if (quant < min_quant) min_quant = quant;
		mean_quant += quant;
		count++;
	}
	fclose(in);

	mean_psnr = mean_psnr / (long double)count;
	mean_quant = mean_quant / (long double)count;

	in = fopen(filename, "r");
	if (in == NULL)
	{
		perror(filename);
		return EXIT_FAILURE;
	}

	while ((n = csv_row_fread(in, buf, BUF_MAX, row, ROW_MAX, ',', CSV_TRIM | CSV_QUOTES)) > 0)
	{
		psnr = (long double)atof(row[ROW_PSNR]);
		if (psnr<1.0 || psnr>100.0) continue;
		quant = (long double)atof(row[ROW_QUANT]);
		if (quant<1.0 || quant>24.0) continue;
		q = (int)floor(quant);
		if ((long double)q != quant) continue;
		psnr -= mean_psnr;
		sd_psnr += psnr*psnr;
		quant -= mean_quant;
		sd_quant += quant*quant;
	}
	fclose(in);

	sd_psnr = sqrt(sd_psnr / (long double)(count-1));
	sd_quant = sqrt(sd_quant / (long double)(count-1));

	if (argc>2)
	{
		allmae = calloc( count, sizeof(long double) );
		if (allmae)
		{
			for (i=1;i<=count;i++)
			{
				if ((i%10)==0) fprintf(stderr, "%lu:%lu\r", i, count);
	
				sprintf( pictfilename, "%s/%8.8lu.png", argv[2], i );
				if (stat(pictfilename, &statbuf)) break;
				reference_wand = NewMagickWand();
				status = MagickReadImage(reference_wand, pictfilename);
				if (status == MagickFalse) ThrowWandException(reference_wand);
	
				sprintf( pictfilename, "%s/%8.8lu.png", argv[3], i );
				if (stat(pictfilename, &statbuf))
				{
					reference_wand = DestroyMagickWand(reference_wand);
					break;
				}
				magick_wand = NewMagickWand();
				status = MagickReadImage(magick_wand, pictfilename);
				if (status == MagickFalse) ThrowWandException(magick_wand);
	
				status = MagickGetImageDistortion(magick_wand, reference_wand, MeanAbsoluteErrorMetric, &distortion);
				if (status == MagickFalse) ThrowWandException(magick_wand);
				reference_wand = DestroyMagickWand(reference_wand);
				magick_wand = DestroyMagickWand(magick_wand);
	
				mae = distortion;
				if (mae > max_mae) max_mae = mae;
				if (mae < min_mae) min_mae = mae;
				mean_mae += mae;
				allmae[i-1] = mae;
				pictcount++;
			}
			mean_mae = mean_mae / (long double)pictcount;
			for (i=0;i<pictcount;i++)
			{
				mae = allmae[i] - mean_mae;
				sd_mae += mae*mae;
			}
			sd_mae = sqrt(sd_mae / (long double)(pictcount-1));
	
			free(allmae);
		}
		else
		{
			fputs("Not enough memory.\n", stderr);
			return EXIT_FAILURE;
		}
	}

	fprintf(stderr, "COUNT: %lu/%lu\nFREQ: ", pictcount, count);
	for ( q=0; q<25; q++ ) fprintf(stderr, " %lu", q_freq[q]);
	fputs("\n#\t\t\tPSNR\t\t\t\tMAE\t\t\t\tQUANT\n\t\t\tmin\tmean\tmax\tsd\tmin\tmean\tmax\tsd\tmin\tmean\tmax\tsd\n", stderr);
	printf("%s\t", argv[1]);
	printf("%.4Lg\t%.4Lg\t%.4Lg\t%.4Lg\t", min_psnr, mean_psnr, max_psnr, sd_psnr);
	printf("%.4Lg\t%.4Lg\t%.4Lg\t%.4Lg\t", min_mae, mean_mae, max_mae, sd_mae);
	printf("%.4Lg\t%.4Lg\t%.4Lg\t%.4Lg\n", min_quant, mean_quant, max_quant, sd_quant);
	return EXIT_SUCCESS;
}