/*
 * scene.cc: screen sample generation and image reconstruction
 *
 * This file is part of Pyrit Ray Tracer.
 *
 * Copyright 2008  Radek Brich
 *
 * 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 <math.h>

#include "common.h"
#include "scene.h"

void DefaultSampler::init()
{
	phase = 0;
	packetable = (subsample <= 1);
}

int DefaultSampler::initSampleSet()
{
	static const int gridsamples[] = {1,5,9,16};
	const int samples = gridsamples[oversample];
	if ( phase == 0 )
	{
		if (subsample > 1)
		{
			phase = 1;
			sx = -1;
			return (w/subsample+1)*(h/subsample+1);
		}
		else
		{
			phase = 2;
			sx = -1;
			return w*h*samples;
		}
	}
	if ( phase == 1 )
	{
		// finalize subsampling
		const Float subsample2 = 1.0/(subsample*subsample);
		int num_samples = 0;
		Colour ic;
		phase = 2;
		sx = -1;
		for (int y = 0; y < h/subsample; y++)
			for (int x = 0; x < w/subsample; x++)
			{
				int x1 = x*subsample;
				int y1 = y*subsample;
				int x2 = (x+1)*subsample;
				int y2 = (y+1)*subsample;
				if (x2 > w-1)	x2 = w-1;
				if (y2 > h-1)	y2 = h-1;
				if (x1 == x2 || y1 == y2)
					continue;
				Float *p;
				p = buffer + 3*(y1*w + x1);
				Colour c1(*p, *(p+1), *(p+2));
				p = buffer + 3*(y1*w + x2);
				Colour c2(*p, *(p+1), *(p+2));
				p = buffer + 3*(y2*w + x1);
				Colour c3(*p, *(p+1), *(p+2));
				p = buffer + 3*(y2*w + x2);
				Colour c4(*p, *(p+1), *(p+2));
				Float m = (c1-c2).mag2();
				m = max(m, (c2-c3).mag2());
				m = max(m, (c3-c4).mag2());
				m = max(m, (c4-c1).mag2());
				if (m < 0.002)
				{
					// interpolate
					for (int i = 0; i < subsample; i++)
						for (int j = 0; j < subsample; j++)
						{
							ic = c1*(subsample-i)*(subsample-j)*subsample2
								+ c2*(i)*(subsample-j)*subsample2
								+ c3*(subsample-i)*(j)*subsample2
								+ c4*(i)*(j)*subsample2;
							p = buffer + 3*((y1+j)*w + x1+i);
							*(p + 0) = ic.r;
							*(p + 1) = oversample ? -ic.g : ic.g;
							*(p + 2) = ic.b;
						}
				}
				else
				{
					// mark as to be computed
					num_samples += subsample * subsample;
					for (int i = 0; i < subsample; i++)
						for (int j = 0; j < subsample; j++)
							if (oversample || i != 0 || j != 0)
								*(buffer + 3*((y1+j)*w + x1+i)) = -1.;
				}
			}
		return num_samples;
	}
	if ( phase == 2 && oversample )
	{
		// finalize oversampling
		Float *buf;
		if (subsample > 1)
			for (buf = buffer; buf != buffer + w*h*3; buf += 3)
				if (*(buf+1) < 0)
				{
					// interpolated
					*(buf+1) = -*(buf+1);
				}
				else
				{
					*buf = *buf * (1.0/samples);
					*(buf+1) = *(buf+1) * (1.0/samples);
					*(buf+2) = *(buf+2) * (1.0/samples);
				}
		else
			for (buf = buffer; buf != buffer + w*h*3; buf++)
				*buf = *buf * (1.0/samples);
	}
	phase = -1;
	return 0;
}

bool DefaultSampler::nextSample(Sample* s)
{
	if (phase == 1)
	{
		// subsampling
		if (sx < 0)
		{
			// first sample
			s->x = -(Float)w/h/2.0;
			s->y = -0.5;
			sx = 0;
			sy = 0;
			osa_samp = 0;
		}
		else
		{
			if (sx == w-1)
			{
				if (sy == h-1)
					return false;
				sy += subsample;
				if (sy > h-1)
					sy = h-1;
				sx = 0;
			}
			else
			{
				sx += subsample;
				if (sx > w-1)
					sx = w-1;
			}

			s->x = (Float)sx/h - (Float)w/h/2.0;
			s->y = (Float)sy/h - 0.5;
		}
	}
	else if (phase == 2)
	{
		/* grid oversampling */
		static const int gridsamples[] = {1,5,9,16};
		static const Float osa5x[] = {0.0, -0.4, +0.4, +0.4, -0.4};
		static const Float osa5y[] = {0.0, -0.4, -0.4, +0.4, +0.4};
		static const Float osa9x[] = {-0.34,  0.00, +0.34,
			-0.34,  0.00, +0.34, -0.34,  0.00, +0.34};
		static const Float osa9y[] = {-0.34, -0.34, -0.34,
				0.00,  0.00,  0.00, +0.34, +0.34, +0.34};
		static const Float osa16x[] = {-0.375, -0.125, +0.125, +0.375,
			-0.375, -0.125, +0.125, +0.375, -0.375, -0.125, +0.125, +0.375,
			-0.375, -0.125, +0.125, +0.375};
		static const Float osa16y[] = {-0.375, -0.375, -0.375, -0.375,
			-0.125, -0.125, -0.125, -0.125, +0.125, +0.125, +0.125, +0.125,
			+0.375, +0.375, +0.375, +0.375};
		static const Float *osaSx[] = {NULL, osa5x, osa9x, osa16x};
		static const Float *osaSy[] = {NULL, osa5y, osa9y, osa16y};
		const int samples = gridsamples[oversample];
		const Float *osax = osaSx[oversample];
		const Float *osay = osaSy[oversample];

		if (sx < 0)
		{
			// first sample
			s->x = -(Float)w/h/2.0;
			s->y = -0.5;
			sx = 0;
			sy = 0;
			osa_samp = 0;
		}
		else
		{
			osa_samp++;

			if (oversample && oversample <= 3 && osa_samp < samples)
			{
				s->x = osax[osa_samp]/h + (Float)sx/h - (Float)w/h/2.0;
				s->y = osay[osa_samp]/h + (Float)sy/h - 0.5;
			}
			else
			{
				sx++;
				if (sx >= w)
				{
					sx = 0;
					sy++;
				}
				if (sy >= h)
					return false;

				if (subsample > 1)
				{
					// find next not interpolated pixel
					while ( *(buffer + 3*(sy*w + sx)) >= 0. )
					{
						sx++;
						if (sx >= w)
						{
							sx = 0;
							sy++;
						}
						if (sy >= h)
							return false;
					}
				}

				s->x = (Float)sx/h - (Float)w/h/2.0;
				s->y = (Float)sy/h - 0.5;
				osa_samp = 0;
			}
		}

		if (osa_samp == 0 && oversample && oversample <= 3)
		{
			s->x += osax[0]/h;
			s->y += osay[0]/h;
			Float *buf = buffer + 3*(sy*w + sx);
			*(buf++) = 0;
			*(buf++) = 0;
			*(buf++) = 0;
		}
	}

	s->sx = sx;
	s->sy = sy;
	s->osa_samp = osa_samp;

	return true;
}

void DefaultSampler::saveSample(Sample &samp, Colour &col)
{
	Float *buf = buffer + 3*(samp.sy * w + samp.sx);
	if (phase == 2 && oversample)
	{
		*(buf+0) += col.r;
		*(buf+1) += col.g;
		*(buf+2) += col.b;
	}
	else
	{
		*(buf++) = col.r;
		*(buf++) = col.g;
		*(buf++) = col.b;
	}
}