1 /********************************************************************************
2  *
3  * darksky.cc: SkyLight, "Real" Sunlight and Sky Background
4  *
5  *		Created on: 20/03/2009
6  *
7  *  	Based on the original implementation by Alejandro Conty (jandro), Mathias Wein (Lynx), anyone else???
8  *  	Actual implementation by Rodrigo Placencia (Darktide)
9  *
10  * 		Based on 'A Practical Analytic Model For DayLight" by Preetham, Shirley & Smits.
11  * 		http://www.cs.utah.edu/vissim/papers/sunsky/
12  * 		based on the actual code by Brian Smits
13  *
14  */
15 
16 #include <yafray_constants.h>
17 #include <core_api/logging.h>
18 #include <core_api/environment.h>
19 #include <core_api/background.h>
20 #include <core_api/params.h>
21 #include <core_api/scene.h>
22 #include <core_api/light.h>
23 
24 #include <utilities/ColorConv.h>
25 #include <utilities/spectralData.h>
26 #include <utilities/curveUtils.h>
27 
28 __BEGIN_YAFRAY
29 
30 class darkSkyBackground_t: public background_t
31 {
32 	public:
33 		darkSkyBackground_t(const point3d_t dir, float turb, float pwr, float skyBright, bool clamp, float av, float bv, float cv, float dv, float ev,
34 							float altitude, bool night, float exp, bool genc, ColorSpaces cs, bool ibl, bool with_caustic);
35 		virtual color_t operator() (const ray_t &ray, renderState_t &state, bool from_postprocessed=false) const;
36 		virtual color_t eval(const ray_t &ray, bool from_postprocessed=false) const;
37 		virtual ~darkSkyBackground_t();
38 		static background_t *factory(paraMap_t &,renderEnvironment_t &);
hasIBL()39 		bool hasIBL() { return withIBL; }
shootsCaustic()40 		bool shootsCaustic() { return shootCaustic; }
41 		color_t getAttenuatedSunColor();
42 
43 	protected:
44 		color_t getSkyCol(const ray_t &ray) const;
45 		double PerezFunction(const double *lam, double cosTheta, double gamma, double cosGamma, double lvz) const;
46 		double prePerez(const double *perez);
47 
48 		color_t getSunColorFromSunRad();
49 
50 		vector3d_t sunDir;
51 		double thetaS;
52 		double theta2, theta3;
53 		double sinThetaS, cosThetaS, cosTheta2;
54 		double T, T2;
55 		double zenith_Y, zenith_x, zenith_y;
56 		double perez_Y[6], perez_x[6], perez_y[6];
57 		float power;
58 		float skyBrightness;
59 		ColorConv convert;
60 		float alt;
61 		bool nightSky;
62 		bool withIBL;
63 		bool shootCaustic;
64 		bool shootDiffuse;
65 };
66 
darkSkyBackground_t(const point3d_t dir,float turb,float pwr,float skyBright,bool clamp,float av,float bv,float cv,float dv,float ev,float altitude,bool night,float exp,bool genc,ColorSpaces cs,bool ibl,bool with_caustic)67 darkSkyBackground_t::darkSkyBackground_t(const point3d_t dir, float turb, float pwr, float skyBright, bool clamp,float av, float bv, float cv, float dv, float ev,
68 										float altitude, bool night, float exp, bool genc, ColorSpaces cs, bool ibl, bool with_caustic):
69 									   power(pwr * skyBright), skyBrightness(skyBright), convert(clamp, genc, cs, exp), alt(altitude), nightSky(night), withIBL(ibl), shootCaustic(with_caustic)
70 {
71 
72 
73 	std::string act = "";
74 
75 	sunDir = vector3d_t(dir);
76 	sunDir.z += alt;
77 	sunDir.normalize();
78 
79 	thetaS = fAcos(sunDir.z);
80 
81 	act = (nightSky)?"ON":"OFF";
82 	Y_VERBOSE << "DarkSky: Night mode [ " << act << " ]" << yendl;
83 	Y_VERBOSE << "DarkSky: Solar Declination in Degrees (" << radToDeg(thetaS) << ")" << yendl;
84 	act = (clamp)?"active.":"inactive.";
85 	Y_VERBOSE << "DarkSky: RGB Clamping " << act << yendl;
86 	Y_VERBOSE << "DarkSky: Altitude " << alt << yendl;
87 
88 	cosThetaS = fCos(thetaS);
89 	cosTheta2 = cosThetaS * cosThetaS;
90 	sinThetaS = fSin(thetaS);
91 
92 	theta2 = thetaS*thetaS;
93 	theta3 = theta2*thetaS;
94 
95 	T = turb;
96 	T2 = turb*turb;
97 
98 	double chi = (0.44444444 - (T / 120.0)) * (M_PI - (2.0 * thetaS));
99 
100 	zenith_Y = (4.0453 * T - 4.9710) * tan(chi) - 0.2155 * T + 2.4192;
101 	zenith_Y *= 1000;  // conversion from kcd/m^2 to cd/m^2
102 
103 	zenith_x =
104 	( 0.00165*theta3 - 0.00374*theta2 + 0.00209*thetaS          ) * T2 +
105 	(-0.02902*theta3 + 0.06377*theta2 - 0.03202*thetaS + 0.00394) * T  +
106 	( 0.11693*theta3 - 0.21196*theta2 + 0.06052*thetaS + 0.25885);
107 
108 	zenith_y =
109 	( 0.00275*theta3 - 0.00610*theta2 + 0.00316*thetaS          ) * T2 +
110 	(-0.04214*theta3 + 0.08970*theta2 - 0.04153*thetaS + 0.00515) * T  +
111 	( 0.15346*theta3 - 0.26756*theta2 + 0.06669*thetaS + 0.26688);
112 
113 	perez_Y[0] = (( 0.17872 * T) - 1.46303) * av;
114 	perez_Y[1] = ((-0.35540 * T) + 0.42749) * bv;
115 	perez_Y[2] = ((-0.02266 * T) + 5.32505) * cv;
116 	perez_Y[3] = (( 0.12064 * T) - 2.57705) * dv;
117 	perez_Y[4] = ((-0.06696 * T) + 0.37027) * ev;
118 	perez_Y[5] = prePerez(perez_Y);
119 
120 	perez_x[0] = ((-0.01925 * T) - 0.25922);
121 	perez_x[1] = ((-0.06651 * T) + 0.00081);
122 	perez_x[2] = ((-0.00041 * T) + 0.21247);
123 	perez_x[3] = ((-0.06409 * T) - 0.89887);
124 	perez_x[4] = ((-0.00325 * T) + 0.04517);
125 	perez_x[5] = prePerez(perez_x);
126 
127 	perez_y[0] = ((-0.01669 * T) - 0.26078);
128 	perez_y[1] = ((-0.09495 * T) + 0.00921);
129 	perez_y[2] = ((-0.00792 * T) + 0.21023);
130 	perez_y[3] = ((-0.04405 * T) - 1.65369);
131 	perez_y[4] = ((-0.01092 * T) + 0.05291);
132 	perez_y[5] = prePerez(perez_y);
133 }
134 
getAttenuatedSunColor()135 color_t darkSkyBackground_t::getAttenuatedSunColor()
136 {
137 	color_t lightColor(1.0);
138 
139 	lightColor = getSunColorFromSunRad();
140 
141 	if(nightSky)
142 	{
143 		lightColor *= color_t(0.8,0.8,1.0);
144 	}
145 
146 	return lightColor;
147 }
148 
getSunColorFromSunRad()149 color_t darkSkyBackground_t::getSunColorFromSunRad()
150 {
151 	int L;
152 	double uL;
153 	double kgLm, kwaLmw, mw;
154     double Rayleigh, Angstrom, Ozone, Gas, Water, m, lm, m1, mB, am, m4;
155 	color_t sXYZ(0.0);
156 	color_t spdf(0.0);
157 
158 	double B = (0.04608365822050 * T) - 0.04586025928522;
159 	double a = 1.3;
160 	double l = 0.35;
161 	double w = 2.0;
162 
163 	IrregularCurve ko(koAmplitudes, koWavelengths, 64);
164     IrregularCurve kg(kgAmplitudes, kgWavelengths, 4);
165     IrregularCurve kwa(kwaAmplitudes, kwaWavelengths, 13);
166     RegularCurve sunRadianceCurve(sunRadiance, 380, 750, 38);
167 
168 	m = 1.0 / (cosThetaS + 0.15 * fPow(93.885f - radToDeg(thetaS), -1.253f));
169 	mw = m * w;
170 	lm = -m * l;
171 
172 	m1 = -0.008735;
173 	mB = -B;
174 	am = -a * m;
175 	m4 = -4.08 * m;
176 
177 	for(L = 380; L < 750; L+=5)
178 	{
179 		uL = L * 0.001;
180 		kgLm = kg(L) * m;
181 		kwaLmw = kwa(L) * mw;
182 
183 		Rayleigh = fExp(m1 * fPow(uL, m4));
184 		Angstrom = fExp(mB * fPow(uL, am));
185 		Ozone = fExp(ko(L) * lm);
186 		Gas = fExp((-1.41 * kgLm) / fPow(1 + 118.93 * kgLm, 0.45));
187 		Water = fExp((-0.2385 * kwaLmw) / fPow(1 + 20.07 * kwaLmw, 0.45));
188 		spdf = sunRadianceCurve(L) * Rayleigh * Angstrom * Ozone * Gas * Water;
189 		sXYZ += chromaMatch(L) * spdf * 0.013513514;
190 	}
191 
192 	return convert.fromXYZ(sXYZ, true);
193 }
194 
~darkSkyBackground_t()195 darkSkyBackground_t::~darkSkyBackground_t()
196 {
197 	// Empty
198 }
199 
prePerez(const double * perez)200 double darkSkyBackground_t::prePerez(const double *perez)
201 {
202 	double pNum = ((1 + perez[0] * fExp(perez[1])) * (1 +( perez[2] * fExp(perez[3] * thetaS) ) + (perez[4] * cosTheta2)));
203 	if(pNum == 0.0) return 0.0;
204 
205 	return 1.0 / pNum;
206 }
207 
PerezFunction(const double * lam,double cosTheta,double gamma,double cosGamma2,double lvz) const208 double darkSkyBackground_t::PerezFunction(const double *lam, double cosTheta, double gamma, double cosGamma2, double lvz) const
209 {
210 	double num = ( (1 + lam[0] * fExp(lam[1]/cosTheta) ) * (1 + lam[2] * fExp(lam[3]*gamma)  + lam[4] * cosGamma2));
211 	return lvz * num * lam[5];
212 }
213 
getSkyCol(const ray_t & ray) const214 inline color_t darkSkyBackground_t::getSkyCol(const ray_t &ray) const
215 {
216 	vector3d_t Iw = ray.dir;
217 	Iw.z += alt;
218 	Iw.normalize();
219 
220 	double cosTheta, gamma, cosGamma, cosGamma2;
221 	double x, y, Y;
222 	color_t skyCol(0.0);
223 
224 	cosTheta = Iw.z;
225 
226 	if(cosTheta <= 0.0) cosTheta = 1e-6;
227 
228 	cosGamma = Iw * sunDir;
229     cosGamma2 = cosGamma * cosGamma;
230 	gamma = fAcos(cosGamma);
231 
232 	x = PerezFunction(perez_x, cosTheta, gamma, cosGamma2, zenith_x);
233 	y = PerezFunction(perez_y, cosTheta, gamma, cosGamma2, zenith_y);
234 	Y = PerezFunction(perez_Y, cosTheta, gamma, cosGamma2, zenith_Y) * 6.66666667e-5;
235 
236 	skyCol = convert.fromxyY(x,y,Y);
237 
238 	if(nightSky) skyCol *= color_t(0.05,0.05,0.08);
239 
240 	return skyCol * skyBrightness;
241 }
242 
operator ()(const ray_t & ray,renderState_t & state,bool from_postprocessed) const243 color_t darkSkyBackground_t::operator() (const ray_t &ray, renderState_t &state, bool from_postprocessed) const
244 {
245 	color_t ret = getSkyCol(ray);
246 	return ret;
247 }
248 
eval(const ray_t & ray,bool from_postprocessed) const249 color_t darkSkyBackground_t::eval(const ray_t &ray, bool from_postprocessed) const
250 {
251 	color_t ret = getSkyCol(ray)  * power;
252 	return ret;
253 }
254 
factory(paraMap_t & params,renderEnvironment_t & render)255 background_t *darkSkyBackground_t::factory(paraMap_t &params,renderEnvironment_t &render)
256 {
257 	point3d_t dir(1,1,1);
258 	float turb = 4.0;
259 	float altitude = 0.0;
260 	int bgl_samples = 8;
261 	float power = 1.0; //bgLight Power
262 	float pw = 1.0;// sunLight power
263 	float bright = 1.0;
264 	bool add_sun = false;
265 	bool bgl = false;
266 	bool clamp = true;
267 	bool night = false;
268 	float av, bv, cv, dv, ev;
269 	bool caus = true;
270 	bool diff = true;
271 	av = bv = cv = dv = ev = 1.0;
272 	bool gammaEnc = true;
273 	std::string cs = "CIE (E)";
274 	float exp = 1.f;
275 	bool castShadows = true;
276 	bool castShadowsSun = true;
277 
278 	Y_VERBOSE << "DarkSky: Begin" << yendl;
279 
280 	params.getParam("from", dir);
281 	params.getParam("turbidity", turb);
282 	params.getParam("altitude", altitude);
283 	params.getParam("power", power);
284 	params.getParam("bright", bright);
285 
286 	//params.getParam("clamp_rgb", clamp);
287 	params.getParam("exposure", exp);
288 	//params.getParam("gamma_enc", gammaEnc);
289 	params.getParam("color_space", cs);
290 
291 	params.getParam("a_var", av); //Darkening or brightening towards horizon
292 	params.getParam("b_var", bv); //Luminance gradient near the horizon
293 	params.getParam("c_var", cv); //Relative intensity of circumsolar region
294 	params.getParam("d_var", dv); //Width of circumsolar region
295 	params.getParam("e_var", ev); //Relative backscattered light
296 
297 	params.getParam("add_sun", add_sun);
298 	params.getParam("sun_power", pw);
299 
300 	params.getParam("background_light", bgl);
301 	params.getParam("with_caustic", caus);
302 	params.getParam("with_diffuse", diff);
303 	params.getParam("light_samples", bgl_samples);
304 	params.getParam("cast_shadows", castShadows);
305 	params.getParam("cast_shadows_sun", castShadowsSun);
306 
307 	params.getParam("night", night);
308 
309 	ColorSpaces colorS = cieRGB_E_CS;
310 	if(cs == "CIE (E)") colorS = cieRGB_E_CS;
311 	else if(cs == "CIE (D50)") colorS = cieRGB_D50_CS;
312 	else if(cs == "sRGB (D65)") colorS = sRGB_D65_CS;
313 	else if(cs == "sRGB (D50)") colorS = sRGB_D50_CS;
314 
315 	if(night)
316 	{
317 		bright *= 0.5;
318 		pw *= 0.5;
319 	}
320 
321 	darkSkyBackground_t *darkSky = new darkSkyBackground_t(dir, turb, power, bright, clamp, av, bv, cv, dv, ev,
322 																altitude, night, exp, gammaEnc, colorS, bgl, caus);
323 
324 	if (add_sun && radToDeg(fAcos(dir.z)) < 100.0)
325 	{
326 		vector3d_t d(dir);
327 		d.normalize();
328 
329 		color_t suncol = darkSky->getAttenuatedSunColor();
330 		double angle = 0.5 * (2.0 - d.z);
331 
332 		Y_VERBOSE << "DarkSky: SunColor = " << suncol << yendl;
333 
334 		paraMap_t p;
335 		p["type"] = std::string("sunlight");
336 		p["direction"] = point3d_t(d);
337 		p["color"] = suncol;
338 		p["angle"] = parameter_t(angle);
339 		p["power"] = parameter_t(pw);
340 		p["samples"] = bgl_samples;
341 		p["cast_shadows"] = castShadowsSun;
342 		p["with_caustic"] = caus;
343 		p["with_diffuse"] = diff;
344 
345 		Y_VERBOSE << "DarkSky: Adding a \"Real Sun\"" << yendl;
346 
347 		light_t *light = render.createLight("DarkSky_RealSun", p);
348 
349 		if(light) render.getScene()->addLight(light);
350 	}
351 
352 	if(bgl)
353 	{
354 		paraMap_t bgp;
355 		bgp["type"] = std::string("bglight");
356 		bgp["samples"] = bgl_samples;
357 		bgp["with_caustic"] = caus;
358 		bgp["with_diffuse"] = diff;
359 		bgp["cast_shadows"] = castShadows;
360 
361 		Y_VERBOSE << "DarkSky: Adding background light" << yendl;
362 
363 		light_t *bglight = render.createLight("DarkSky_bgLight", bgp);
364 
365 		bglight->setBackground(darkSky);
366 
367 		if(bglight) render.getScene()->addLight(bglight);
368 	}
369 
370 	Y_VERBOSE << "DarkSky: End" << yendl;
371 
372 	return darkSky;
373 }
374 
375 extern "C"
376 {
377 
registerPlugin(renderEnvironment_t & render)378 	YAFRAYPLUGIN_EXPORT void registerPlugin(renderEnvironment_t &render)
379 	{
380 		render.registerFactory("darksky",darkSkyBackground_t::factory);
381 	}
382 
383 }
384 __END_YAFRAY
385