/*	$Id: chaleur.c,v 1.4 2008/03/20 08:16:10 minh Exp $	*/

/*-
 * Copyright (c) 2008 Nhat Minh Lê
 *
 * Permission is hereby granted, free of charge, to any person
 * obtaining a copy of this software and associated documentation
 * files (the "Software"), to deal in the Software without
 * restriction, including without limitation the rights to use, copy,
 * modify, merge, publish, distribute, sublicense, and/or sell copies
 * of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */


#include <float.h>
#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>


#ifndef MAXWIDTH
# define MAXWIDTH 400
#endif

#ifndef MAXHEIGHT
# define MAXHEIGHT MAXWIDTH
#endif

#ifndef PRECISION
# define PRECISION (DBL_EPSILON * 100.0)
#endif

struct cell {
	double value;
	unsigned disabled: 1;
};

struct state {
	struct cell *matrix;
	double maxval;
	unsigned int width, height;
	unsigned int niter;
};


static double precision = PRECISION;
static unsigned int maxiter = 100, refreshmod = 0;
static unsigned int oprec = 6;



/* Utilities */

static void
error(const char *fmt, ...)
{
	va_list ap;

	fprintf(stderr, "error: ");

	va_start(ap, fmt);
	vfprintf(stderr, fmt, ap);
	va_end(ap);

	fprintf(stderr, "\n");

	exit(EXIT_FAILURE);
}

static void *
xmalloc(size_t n)
{
	void *p;
	if ((p = malloc(n)) == NULL)
		error(0, "out of memory");
	return p;
}



/* Operations */

#define PT(m, w, x, y) (&(m)[(x)*(w) + (y)])
#define GETPT(m, w, x, y) ((m)[(x)*(w) + (y)].value)
#define SETPT(m, w, x, y, v) ((m)[(x)*(w) + (y)].value = (v))

static struct state *
clonestate(struct state *src)
{
	struct state *s;
	unsigned int msize;

	s = xmalloc(sizeof *s);
	s->width = src->width, s->height = src->height;
	msize = src->width * src->height;

	s->matrix = xmalloc(sizeof *s->matrix * msize);
	memcpy(s->matrix, src->matrix, sizeof *s->matrix * msize);

	s->niter = 0;

	return s;
}

static void
freestate(struct state *s)
{
	free(s->matrix);
	free(s);
}

static int
cmpmatrixes(struct cell *a, struct cell *b, unsigned int w, unsigned int h)
{
	unsigned int x, y;
	double r;

	for (x = 0; x < h; ++x)
		for (y = 0; y < w; ++y) {
			r = GETPT(a, w, x, y) - GETPT(b, w, x, y);
			if (r < -precision)
				return -1;
			else if (r > precision)
				return 1;
		}
	return 0;
}

static int
cmpstates(struct state *a, struct state *b)
{
	return cmpmatrixes(a->matrix, b->matrix, a->width, b->width);
}



/* Loader */

static struct cell *
loadmatrix(FILE *fp, unsigned int w, unsigned int h, double *maxptr)
{
	struct cell *m;
	double max, val;
	unsigned int x, y;

	max = 0.0;
	m = xmalloc(sizeof *m * w * h);
	for (x = 0; x < h; ++x)
		for (y = 0; y < w; ++y) {
			if (fscanf(fp, "%lf", &val) < 1)
				error("syntax error in matrix definition");

			SETPT(m, w, x, y, val);
			PT(m, w, x, y)->disabled = 0;

			if (val > max)
				max = val;
		}

	if (maxptr != NULL)
		*maxptr = max;
	return m;
}

static struct state *
loadstate(FILE *fp)
{
	struct state *s;
	double max;
	unsigned int w, h;

	s = xmalloc(sizeof *s);

	if (fscanf(fp, "%u %u", &w, &h) < 2)
		error("syntax error in matrix declaration");
	if (w < 2 || h < 2)
		error("matrix is too small (min: 2x2)");
	if (w > MAXWIDTH || h > MAXHEIGHT)
		error("matrix is too large (max: %ux%u)", MAXWIDTH, MAXHEIGHT);

	s->width = w, s->height = h;
	s->matrix = loadmatrix(fp, w, h, &max);

	s->maxval = max;
	s->niter = 0;

	return s;
}



/* Text output */

static void
printmatrix(struct cell *m, unsigned int w, unsigned int h)
{
	unsigned int x, y;

	for (x = 0; x < h; ++x) {
		for (y = 0; y < w; ++y)
			printf("%.*f\t", (int)oprec, GETPT(m, w, x, y));
		printf("\n");
	}
}

static void
printstate(struct state *s)
{
	printf("#%u:\n", s->niter);
	printmatrix(s->matrix, s->width, s->height);
}



/* Dithering */

static void
dithermatrix(struct cell *m, struct cell *src, unsigned int w, unsigned int h)
{
	unsigned int x, y;

	for (x = 1; x < h-1; ++x)
		for (y = 1; y < w-1; ++y)
			if (!PT(src, w, x, y)->disabled)
				SETPT(m, w, x, y,
				      (GETPT(src, w, x-1, y-1) +
				       GETPT(src, w, x-1, y+1) +
				       GETPT(src, w, x+1, y-1) +
				       GETPT(src, w, x+1, y+1)) / 4);
}

static void
dither(struct state *s)
{
	struct cell *m;

	m = xmalloc(sizeof *m * s->width * s->height);

	memcpy(m, s->matrix, sizeof *m * s->width * s->height);
	dithermatrix(m, s->matrix, s->width, s->height);

	free(s->matrix);
	s->matrix = m;

	++s->niter;
}



/* Main program */

static void
usage(void)
{
	fprintf(stderr, "Usage:\n%s%s%s%s%s",
		"\tchaleur [OPTIONS]\n",
		"-D prec	output precision\n",
		"-R modulo	output every modulo passes\n",
		"-n maxiter	limit on the number of iterations\n",
		"-p prec	requested precision\n");
}

static char *
optarg(char ***argvptr)
{
	if ((*argvptr)[1] != NULL)
		return *++*argvptr;
	else
		return NULL;
}

static void
parseopts(int *argcptr, char ***argvptr)
{
	char *arg;

	--*argcptr;
	++*argvptr;
	while (**argvptr != NULL && (**argvptr)[0] == '-') {
		switch ((**argvptr)[1]) {
		case 'D':
			if ((arg = optarg(argvptr)) != NULL)
				oprec = strtoul(arg, NULL, 0);
			break;
		case 'R':
			if ((arg = optarg(argvptr)) != NULL)
				refreshmod = strtoul(arg, NULL, 0);
			break;
		case 'n':
			if ((arg = optarg(argvptr)) != NULL)
				maxiter = strtoul(arg, NULL, 0);
			break;
		case 'p':
			if ((arg = optarg(argvptr)) != NULL)
				precision = strtod(arg, NULL);
			break;
		default:
			usage();
			exit(0);
		}
		--*argcptr;
		++*argvptr;
	}
}

int
main(int argc, char *argv[])
{
	struct state *s, *t;
	unsigned int i;

	parseopts(&argc, &argv);

	if (refreshmod == 0 || refreshmod > maxiter)
		refreshmod = maxiter / 10;

	s = loadstate(stdin);
	t = clonestate(s);

	printf("%u %u %f\n", s->width, s->height, s->maxval);
	for (i = 0; i < maxiter; ++i) {
		if (i % refreshmod == 0)
			printstate(s);

		dither(s);
		dither(t), dither(t);

		if (cmpstates(s, t) == 0) {
			printstate(s);
			break;
		}
	}

	freestate(s);
	freestate(t);

	return 0;
}