/*
	This file is part of Mitsuba, a physically based rendering system.

	Copyright (c) 2007-2012 by Wenzel Jakob and others.

	Mitsuba is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License Version 3
	as published by the Free Software Foundation.

	Mitsuba 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, see <http://www.gnu.org/licenses/>.
*/

#include <mitsuba/render/scene.h>
#include <mitsuba/core/plugin.h>
#include <mitsuba/core/bitmap.h>
#include <mitsuba/core/qmc.h>
#include "sunsky/sunmodel.h"

#if SPECTRUM_SAMPLES == 3
# define SUNSKY_PIXELFORMAT Bitmap::ERGB
#else
# define SUNSKY_PIXELFORMAT Bitmap::ESpectrum
#endif

/* Apparent radius of the sun as seen from the earth (in degrees). 
   This is an approximation--the actual value is somewhere between 
   0.526 and 0.545 depending on the time of year */
#define SUN_APP_RADIUS 0.5358

MTS_NAMESPACE_BEGIN

/*!\plugin{sunsky}{Sun and sky emitter}
 * \icon{emitter_sunsky}
 * \order{8}
 * \parameters{
 *     \parameter{turbidity}{\Float}{
 *         This parameter determines the amount of aerosol present in the atmosphere. 
 *         Valid range: 1-10. \default{3, corresponding to a clear sky in a temperate climate}
 *     }
 *     \parameter{albedo}{\Spectrum}{Specifies the ground albedo \default{0.15}}
 *     \parameter{year, month, day}{\Integer}{Denote the date of the
 *      observation \default{2010, 07, 10}}
 *     \parameter{hour,minute,\showbreak second}{\Float}{Local time 
 *       at the location of the observer in 24-hour format\default{15, 00, 00,
 *       i.e. 3PM}}
 *     \parameter{latitude, longitude, timezone}{\Float}{
 *       These three parameters specify the oberver's latitude and longitude 
 *       in degrees, and the local timezone offset in hours, which are required 
 *       to compute the sun's position. \default{35.6894, 139.6917, 9 --- Tokyo, Japan}
 *     }
 *     \parameter{sunDirection}{\Vector}{Allows to manually 
 *       override the sun direction in world space. When this value
 *       is provided, parameters pertaining to the computation 
 *       of the sun direction (\code{year, hour, latitude,} etc. 
 *       are unnecessary. \default{none}
 *     }
 *     \parameter{stretch}{\Float}{
 *         Stretch factor to extend emitter below the horizon, must be
 *         in [1,2] \default{\code{1}, i.e. not used}
 *     }
 *     \parameter{resolution}{\Integer}{Specifies the horizontal resolution of the precomputed
 *         image that is used to represent the sun environment map \default{512, i.e. 512$\times$256}}
 *     \parameter{sunScale}{\Float}{
 *         This parameter can be used to separately scale the the amount of illumination
 *         emitted by the sun. \default{1}
 *     }
 *     \parameter{skyScale}{\Float}{
 *         This parameter can be used to separately scale the the amount of illumination
 *         emitted by the sky.\default{1}
 *     }
 * }
 *
 * \renderings{
 *   \medrendering{\pluginref{sky} emitter}{emitter_sunsky_sky}
 *   \medrendering{\pluginref{sun} emitter}{emitter_sunsky_sun}
 *   \medrendering{\pluginref{sunsky} emitter}{emitter_sunsky_sunsky}
 *   \caption{A coated rough copper test ball lit with the three
 *   provided daylight illumination models}
 * }
 * \vspace{5mm}
 * This convenience plugin has the sole purpose of instantiating 
 * \pluginref{sun} and \pluginref{sky} and merging them into a joint
 * environment map. Please refer to these plugins individually for more
 * details.
 */
class SunSkyEmitter : public Emitter {
public:
	SunSkyEmitter(const Properties &props)
		: Emitter(props) {
		Float scale = props.getFloat("scale", 1.0f),
		      sunScale = props.getFloat("sunScale", scale),
		      skyScale = props.getFloat("skyScale", scale);

		const Transform &trafo = m_worldTransform->eval(0);

		Properties skyProps(props);
		skyProps.removeProperty("toWorld");
		if (props.hasProperty("sunDirection"))
			skyProps.setVector("sunDirection", trafo.inverse()(props.getVector("sunDirection")));
		skyProps.setPluginName("sky");
		skyProps.setFloat("scale", skyScale, false);

		ref<Emitter> sky = static_cast<Emitter *>(
			PluginManager::getInstance()->createObject(
			MTS_CLASS(Emitter), skyProps));
		sky->configure();
		props.markQueried("albedo");

		int resolution = props.getInteger("resolution", 512);
		ref<Bitmap> bitmap = new Bitmap(SUNSKY_PIXELFORMAT, Bitmap::EFloat, 
			Vector2i(resolution, resolution/2));

		Point2 factor((2*M_PI) / bitmap->getWidth(),
			M_PI / bitmap->getHeight());

		ref<Timer> timer = new Timer();
		Log(EDebug, "Rasterizing sun & skylight emitter to an %ix%i environment map ..",
				resolution, resolution/2);
			
		Spectrum *data = (Spectrum *) bitmap->getFloatData();

		/* First, rasterize the sky */
		#if defined(MTS_OPENMP)
			#pragma omp parallel for 
		#endif
		for (int y=0; y<bitmap->getHeight(); ++y) {
			Float theta = (y+.5f) * factor.y;
			Spectrum *target = data + y * bitmap->getWidth();

			for (int x=0; x<bitmap->getWidth(); ++x) {
				Float phi = (x+.5f) * factor.x;

				RayDifferential ray(Point(0.0f),
					toSphere(SphericalCoordinates(theta, phi)), 0.0f);
				
				*target++ = sky->evalEnvironment(ray);
			}
		}

		/* Rasterizing the sphere to an environment map and checking the
		   individual pixels for coverage (which is what Mitsuba 0.3.0 did) 
		   was slow and not very effective; for instance the power varied 
		   dramatically with resolution changes. Since the sphere generally
		   just covers a few pixels, the code below rasterizes it much more
		   efficiently by generating a few thousand QMC samples.

		   Step 1: compute a *very* rough estimate of how many
		   pixel in the output environment map will be covered
		   by the sun */

		SphericalCoordinates sun = computeSunCoordinates(props);
		Spectrum sunRadiance = computeSunRadiance(sun.elevation, 
			props.getFloat("turbidity", 3.0f)) * sunScale;
		sun.elevation *= props.getFloat("stretch", 1.0f);
		Frame sunFrame = Frame(toSphere(sun));
		Float theta = degToRad(SUN_APP_RADIUS * 0.5f);

		size_t pixelCount = resolution*resolution/2;
		Float cosTheta = std::cos(theta * props.getFloat("sunRadiusScale", 1.0f));

		/* Ratio of the sphere that is covered by the sun */
		Float coveredPortion = 0.5f * (1 - cosTheta);

		/* Approx. number of samples that need to be generated, 
		   be very conservative */
		size_t nSamples = (size_t) std::max((Float) 100,
			(pixelCount * coveredPortion * 1000));

		factor = Point2(bitmap->getWidth() / (2*M_PI),
			bitmap->getHeight() / M_PI);

		Spectrum value = 
			sunRadiance * (2 * M_PI * (1-std::cos(theta))) *
			(bitmap->getWidth() * bitmap->getHeight()) 
			/ (2.0f * M_PI * M_PI * (Float) nSamples);

		for (size_t i=0; i<nSamples; ++i) {
			Vector dir = sunFrame.toWorld(
				Warp::squareToUniformCone(cosTheta, sample02(i)));

			Float sinTheta = math::safe_sqrt(1-dir.y*dir.y);
			SphericalCoordinates sphCoords = fromSphere(dir);

			Point2i pos(
				std::min(std::max(0, (int) (sphCoords.azimuth * factor.x)), bitmap->getWidth()-1),
				std::min(std::max(0, (int) (sphCoords.elevation * factor.y)), bitmap->getHeight()-1));

			data[pos.x + pos.y * bitmap->getWidth()] += value / std::max((Float) 1e-3f, sinTheta);
		}

		Log(EDebug, "Done (took %i ms)", timer->getMilliseconds());

		/* Instantiate a nested envmap plugin */
		Properties envProps("envmap");
		Properties::Data bitmapData;
		bitmapData.ptr = (uint8_t *) bitmap.get();
		bitmapData.size = sizeof(Bitmap);
		envProps.setData("bitmap", bitmapData);
		envProps.setTransform("toWorld", trafo);
		envProps.setFloat("samplingWeight", m_samplingWeight);
		m_emitter = static_cast<Emitter *>(
			PluginManager::getInstance()->createObject(
			MTS_CLASS(Emitter), envProps));

		#if 0
			/* For debugging purposes */
			ref<FileStream> fs = new FileStream("debug.exr", FileStream::ETruncReadWrite);
			bitmap->write(Bitmap::EOpenEXR, fs);
		#endif
	}

	SunSkyEmitter(Stream *stream, InstanceManager *manager)
		: Emitter(stream, manager) {
		m_emitter = static_cast<Emitter *>(manager->getInstance(stream));
	}

	void serialize(Stream *stream, InstanceManager *manager) const {
		Emitter::serialize(stream, manager);
		manager->serialize(stream, m_emitter.get());
	}

	void configure() {
		Emitter::configure();
		m_emitter->configure();
	}

	bool isCompound() const {
		return true;
	}

	AABB getAABB() const {
		NotImplementedError("getAABB");
	}

	Emitter *getElement(size_t i) {
		if (i == 0)
			return m_emitter;
		else
			return NULL;
	}

	MTS_DECLARE_CLASS()
private:
	ref<Emitter> m_emitter;
};

MTS_IMPLEMENT_CLASS_S(SunSkyEmitter, false, Emitter)
MTS_EXPORT_PLUGIN(SunSkyEmitter, "Sun & sky emitter");
MTS_NAMESPACE_END