#include "stdafx.h"

#include "Scene.hxx"
#include "DirectionalLight.hxx"
#include <string.h>
#include <array2d.h>

void ComputeMedianCut(const Box& r, const Image2D& sat, int& d1, int& sp);

void ComputeSAT(const Image2D& img, Image2D& sat);

Vec3f GetAreaSum(const Image2D& sat, int u1, int v1, int u2, int v2);

/* Image based ligthing.
 *
 * Add directional lights to an scene based on a high dynamic image.
 * A fixed number of directional lights must be set.  The algorithm
 * computes the location and intensity on each RGB channel based
 * on the luminance given by a high dynamic range image.
 * In the input image, columns and rows correspond to the spherical
 * coordinates phi and theta of an spherical photography of the
 * environment.
 *
 * Reference:
 * Paul Debevec. Rendering Synthetic Objects into Real Scenes: Bridging Traditional and Image-based
 * Graphics with Global Illumination and High Dynamic Range Photography.
 * Computer Graphics (Proceedings of SIGGRAPH 98), 32(4):189-198, 1998.
 *
 * Miguel Granados, Richard Socher
 * Core Lecture in Computer Graphics by Prof. Dr. Philipp Slusallek
 * University of Saarland
 * WS2006/2007
 */

void Scene::AddImageBasedLighting(const Image2D& img0, int nLights) {
	std::list<Box> regions;

	/*
	Computing the total light energy is most naturally performed
	on a monochrome version of the lighting environment rather
	than the RGB pixel colors; such an image can be formed as a
	weighted average of the color channels of the light probe image,
	e.g. Y = 0.2125R+0.7154G+0.0721B following ITU-R Recommendation
	BT.709.
	*/
	Image2D imgMono(img0.resX, img0.resY);
	for(int r=0; r<imgMono.resY; ++r) {
		float cosPhi = cos(((float)r / imgMono.resY * M_PI) - (M_PI / 2));
		for(int c=0; c<imgMono.resX; ++c) {
			Vec3f rgb = img0(c, r);
			float avg = 0.2125*rgb[0] + 0.7154*rgb[1] + 0.0721*rgb[2];
			/*
			To compensate, the pixels of the probe image should first be
			scaled by cosf where f is the pixel’s angle of inclination.
			*/
			avg *= cosPhi;
			imgMono.set(c, r, Vec3f(avg));
		}
	}

	// 1 Add the entire light probe image to the region list as a single region.
	Box r0;
	r0.Extend(Vec3f(0,0,0));
	r0.Extend(Vec3f(imgMono.resX-1, imgMono.resY-1, 0));
	regions.push_back(r0);

	// Calculating the total energy within regions of the image can be
	// accelerated using a summed area table
	Image2D satMono(imgMono.resX, imgMono.resY);
	ComputeSAT(imgMono, satMono);

	// 2 For each region in the list, subdivide along the longest dimension
	// such that its light energy is divided evenly.
	int n = 1;
	while(!regions.empty() && n < nLights) {
		++n;
		Box r = regions.front();
		regions.pop_front();
		int d; // splitting dimension
		int sp; // splitting point
		ComputeMedianCut(r, satMono, d, sp);
		//cerr << "Scene::AddImageBasedLighting: region.min" << r.min << ", region.max=" << r.max << ", d=" << d << ", sp=" << sp << endl;
		Box rLeft, rRight;
		if( d == 0) {
			rLeft.Extend(Vec3f(r.min[0], r.min[1], 0));
			rLeft.Extend(Vec3f(sp-1, r.max[1], 0));
			rRight.Extend(Vec3f(sp, r.min[1], 0));
			rRight.Extend(Vec3f(r.max[0], r.max[1], 0));
		} else {
			rLeft.Extend(Vec3f(r.min[0], r.min[1], 0));
			rLeft.Extend(Vec3f(r.max[0], sp-1, 0));
			rRight.Extend(Vec3f(r.min[0], sp, 0));
			rRight.Extend(Vec3f(r.max[0], r.max[1], 0));
		}
		//cerr << "Scene::AddImageBasedLighting: rLeft.min" << rLeft.min << ", rLeft.max=" << rLeft.max << endl;
		//cerr << "Scene::AddImageBasedLighting: rRight.min" << rRight.min << ", rRight.max=" << rRight.max << endl;
		assert(r.contains(rLeft));
		regions.push_back(rLeft);
		assert(r.contains(rRight));
		regions.push_back(rRight);
	}

	/* Calculating the total energy within regions
	   of the image can be accelerated using a summed area table
	*/
	Image2D satColor(img0.resX, img0.resY);
	ComputeSAT(img0, satColor);

	// compute the centroid of the region and place a directional light
	Image2D tmpImg = img0; // DEBUG CODE

	for(std::list<Box>::const_iterator ri = regions.begin(); ri != regions.end(); ++ri) {
		Box b = *ri;
		Vec3f bOrg = Vec3f((int)b.min[0], (int)b.min[1], 0);
		int bmin[] = {(int)b.min[0], (int)b.min[1]};
		int bmax[] = {(int)b.max[0], (int)b.max[1]};

		/*
		Place a light source at the center or centroid of each region,
		and set the light source color to the sum of pixel values within
		the region.
		*/
		float regionSum = GetAreaSum(satMono, bmin[0], bmin[1], bmax[0], bmax[1])[0];
		if(regionSum > 0) {
			//cerr << "Scene::AddImageBasedLighting: regionSum=" << regionSum << endl;
			Vec3f pixelWeightedSum(0,0,0);
			for(int c = bmin[0]; c < bmax[0]; ++c) {
				for(int r = bmin[1]; r < bmax[1]; ++r) {
					pixelWeightedSum += (Vec3f(c, r, 0) - bOrg) * imgMono(c, r)[0];
				}
			}
			Vec3f centroid = bOrg + (pixelWeightedSum / regionSum);
			float theta = (float)centroid[0] / img0.resX * 2 * M_PI;
			float phi = ((float)centroid[1] / img0.resY * M_PI) - (M_PI / 2);
			Vec3f lightDir(cos(phi)*cos(theta), cos(phi)*sin(theta), sin(phi));
			lightDir = -lightDir;

			/*
			While the partitioning decisions are made on
			the monochrome image, the light source colors are computed using
			the corresponding regions in the original RGB image.
			*/
			Vec3f lightIntensity = GetAreaSum(satColor, bmin[0], bmin[1], bmax[0], bmax[1]) / (img0.resX*img0.resY);

			Add(new DirectionalLight(lightDir, lightIntensity));
			cerr << "Scene::AddImageBasedLighting: lightDir=" << lightDir << ", lightIntensity=" << lightIntensity << endl;
			/*
			cerr << "Scene::AddImageBasedLighting: b.min"
				<< b.min << ", b.max=" << b.max
				<< ", centroid=" << centroid
				<< ", (theta, phi)=(" << theta << ", " << phi << ")"
				<< ", lightDir=" << lightDir << ", lightIntensity=" << lightIntensity << endl;
			*/
			/* DEBUG CODE */
			for(int c = (int)b.min[0]; c < (int)b.max[0]; ++c) {
				for(int r = (int)b.min[1]; r < (int)b.max[1]; ++r) {
					Vec3f v = tmpImg(c, r);
					v = Max(Min(v, Vec3f(1,1,1)), Vec3f(0,0,0));
					tmpImg.set(c, r, v);
				}
			}
			Vec3f red(1,0,0), yellow(1, 1, 0);
			for(int c = (int)b.min[0]; c < (int)b.max[0]; ++c) {
				tmpImg.set(c, (int)b.min[1], red);
				tmpImg.set(c, (int)b.max[1], red);
			}
			for(int r = (int)b.min[1]; r < (int)b.max[1]; ++r) {
				tmpImg.set((int)b.min[0], r, red);
				tmpImg.set((int)b.max[0], r, red);
			}
			tmpImg.set((int)centroid[0], (int)centroid[1], yellow);

			tmpImg.set((int)centroid[0]-1, (int)centroid[1]-1, yellow);
			tmpImg.set((int)centroid[0]-1, (int)centroid[1]+1, yellow);
			tmpImg.set((int)centroid[0]+1, (int)centroid[1]-1, yellow);
			tmpImg.set((int)centroid[0]+1, (int)centroid[1]+1, yellow);

			tmpImg.set((int)centroid[0]-1, (int)centroid[1], yellow);
			tmpImg.set((int)centroid[0]+1, (int)centroid[1], yellow);
			tmpImg.set((int)centroid[0], (int)centroid[1]-1, yellow);
			tmpImg.set((int)centroid[0], (int)centroid[1]+1, yellow);
		}
	}
	/* DEBUG CODE */
	for(int r=0; r<tmpImg.resY; ++r) {
		for(int c=0; c<tmpImg.resX; ++c) {
			Vec3f v = tmpImg(c, r);
			v = Max(Min(v, Vec3f(1,1,1)), Vec3f(0,0,0));
			tmpImg.set(c, r, v);
		}
	}
	tmpImg.WritePPM("medianCutSampling.ppm");

}

// Given an image region, compute the splitting point and plane that
// makes the energy in both sides of the splitting most even
void ComputeMedianCut(const Box& r, const Image2D& sat, int& d1, int& sp) {
	/*
	Determining the longest dimension of a region should also take the
	overrepresentation into account; this can be accomplished by weighting
	a regions width by cosf for an inclination f at center of the region.
	*/
	Vec3f rCenter = r.min + Vec3f(r.dX()/2, r.dY()/2, 0);
	float cosPhi = cos((rCenter[1] / sat.resY * M_PI) - (M_PI / 2));
	//cerr << "ComputeMedianCut: r.dX=" << (r.dX()*cosPhi) << ", r.dY=" << r.dY() << endl;
	d1 = (r.dX()*cosPhi > r.dY()) ? 0: 1; // longest dimension
	int d2 = (d1+1)%2; // shortest dimension
	float minEDiff = Infinity; // minimum energy difference

	int dmin[2] = {(int)r.min[0], (int)r.min[1]};
	int dmax[2] = {(int)r.max[0], (int)r.max[1]};

	for(int di = dmin[d1]; di < dmax[d1]; ++di) {
		Vec3f eLeft; // energy in left subregion
		Vec3f eRight; // energy in right subregion
		if(d1 == 0) {
			eLeft = GetAreaSum(sat, dmin[d1], dmin[d2], di, dmax[d2]);
			eRight = GetAreaSum(sat, di, dmin[d2], dmax[d1], dmax[d2]);
		} else {
			eLeft = GetAreaSum(sat, dmin[d2], dmin[d1], dmax[d2], di);
			eRight = GetAreaSum(sat, dmin[d2], di, dmax[d2], dmax[d1]);
		}
		float EDiff = Length(eLeft - eRight); // energy difference
		if(EDiff < minEDiff) {
			minEDiff = EDiff;
			sp = di;
		}
	}
}

// Compute the sum of areas table for a three channel image
void ComputeSAT(const Image2D& img, Image2D& sat)  {
	for(int r=0; r<img.resY; ++r) {
		for(int c=0; c<img.resX; ++c) {
			Vec3f v = img(c, r);
			if(r > 0)
				v = v + sat(c, r-1);
			if(c > 0)
				v = v + sat(c-1, r);
			if(r > 0 && c > 0)
				v = v - sat(c-1, r-1);
			sat.set(c, r, v);
		}
	}
}

// Obtain the area sum using the SAT table
inline Vec3f GetAreaSum(const Image2D& sat, int u1, int v1, int u2, int v2) {
	return (sat(u2, v2) - sat(u1, v2) - sat(u2, v1) + sat(u1, v1));
}