#include "stdafx.h"

#include "ToneMappers.hxx"

// Tone Mapping, scales HDR colors with Ward's operator
void toneMapping(int resX, int resY, Image2D &img){
	cout << "starting tone mapping" << endl;

	// our normal display only shows color values \in [0,1]
	float maxDisplayLuminance=2.4;
	// decrease world adaption luminance
	float factor = 1;



	// check if we have HDR in the image (this is the case, if any color is >1)
	// because Radiance (as a phyiscally motivated measurement) can be extremely high
	bool toneMappingNeeded = false;
	for (int i=0; i<resY; i++){
		for (int j=0; j<resX; j++){
			if (img(j,i).x>1 || img(j,i).y>1 || img(j,i).z>1){
				toneMappingNeeded = true;
			}
		}
	}	
	if (!toneMappingNeeded) return;




	// calculate the world adaptation Luminance
	// because small values should not change much and due to Stevens' power law
	float LWA=0;
	for (int i=0; i<resY; i++){
		for (int j=0; j<resX; j++){
			LWA = LWA + log((img(j,i).x+img(j,i).y+img(j,i).z)/(3.0f));
		}
	}
	LWA = exp(LWA /(resX*resY)) * factor;



	// Ward's Operator (preserving perceived contrast), as explained in:
	// G. Ward. A contrast-based scalefactor for luminance display. 
	// In Graphics Gems IV, chapter VII.2, p. 415–421
	// (and http://www.cg.tuwien.ac.at/research/publications/2001/Prikryl-thesis/Prikryl-thesis-PDF.pdf)
	// m expresses the minimal discernible luminance changes
	float m = (1/maxDisplayLuminance) * powf(
		((1.219f + powf((maxDisplayLuminance/2),0.4f))
		/(1.219f + powf(LWA,0.4f)))				,2.5f);
	
	cout << m << " = m " << endl;


	// now map the tones
	for (int i=0; i<resY; i++){
		for (int j=0; j<resX; j++){
			img.set(j,i, img(j,i)*m);
		}
	}
}

void WardToneMapping(Image2D &img, Image2D& mask, bool foreground, float LwaFactor) {

	/*	
	When applying display equations such as Ward's Operator, the
	result depends on the choices of the adaptation states Lwa
	and Lda. In the absence of any obviously correct answer,
	we opt for the simplest choice. For the world adaptation we
	choose half the highest visible luminance. For the display observer
	we use half the maximum screen luminance (typically
	80=2 = 40cd=m2).
	*/
	
	// world adaptation level GIVEN
	//  will be computed as the average luminance value in the region of interest
	float Lwa; 

	// maximum display luminance GIVEN
	const float Ldmax = 40.0f; 

	// display adaptation level, deltaJND(Lda)=deltaJND(Ldw), JND: just noticeable difference
	const float Lda = Ldmax/2.0f; 	

	// luminance to be displayed for a pixel, Ld=m*Lw;
	Vec3f Ld; 

	// real world luminance of a pixel
	Vec3f Lw;  
	
	Vec3f imgMin(Infinity, Infinity, Infinity);
	Vec3f imgMax(-Infinity, -Infinity, -Infinity);
	
	// compute Lwa as half of the maximum luminance in the region of interest
	//float LwaAccum = 0, nSamples = 0;
	Vec3f iturY(0.2125, 0.7154, 0.0721);
	for(int c=0; c<img.resX; ++c) {
		for(int r=0; r<img.resY; r++) {
			if((bool)mask(c,r)[0] == foreground) {
				Ld = img(c,r);
				imgMin = Min(imgMin, Ld);
				imgMax = Max(imgMax, Ld);
				//LwaAccum += Dot(Ld, iturY);
				//++nSamples;			
			}
		}
	}
	cerr << "WardToneMapping: imgMin=" << imgMin << ", imgMax=" << imgMax << endl;
	Lwa = Dot(imgMax, iturY) / 2  * LwaFactor;
	//Lwa = LwaAccum / nSamples;	

	// Ward's Operator (preserving perceived contrast), as explained in:
	// G. Ward. A contrast-based scalefactor for luminance display. 
	// In Graphics Gems IV, chapter VII.2, p. 415–421
	// (and http://www.cg.tuwien.ac.at/research/publications/2001/Prikryl-thesis/Prikryl-thesis-PDF.pdf)
	// compute the factor m, which binds together the minimum discernible
	// luminance changes and the deltaJND at the display and world
	// adaptation levels
	float m = (1 / Ldmax) * powf(
								 (1.219f + powf(Lda, 0.4f))
								/(1.219f + powf(Lwa, 0.4f))
								,2.5f);
	
	cout << "WardToneMapping: LwaFactor=" << LwaFactor << ", Lwa=" << Lwa << ", m=" << m << endl;

	Vec3f imgAfterMin(Infinity, Infinity, Infinity);
	Vec3f imgAfterMax(-Infinity, -Infinity, -Infinity);

	// apply the ward's operator to the luminance of each pixel on the image
	for(int c=0; c<img.resX; ++c) {
		for(int r=0; r<img.resY; r++) {
			if((bool)mask(c,r)[0] == foreground) {			
				Lw = img(c,r);
				Ld = m * Lw;
				img.set(c, r, Ld);
				imgAfterMin = Min(imgMin, Ld);
				imgAfterMax = Max(imgAfterMax, Ld); 
			}			
		}
	}
	// transform dynamic range to the [0, 1] interval
	/*
	float sf = 1/imgAfterMax.Max();
	for(int c=0; c<img.resX; ++c) {
		for(int r=0; r<img.resY; r++) {
			if((bool)mask(c,r)[0] == foreground) {
				Ld = img(c,r) * sf;
				img.set(c, r, Ld);
			}
		}
	}
	*/		
	cerr << "WardToneMapping: imgAfterMin=" << imgAfterMin << ", imgAfterMax=" << imgAfterMax << endl; 
}