Cycles / Sky Texture:
[blender.git] / intern / cycles / render / sky_model.cpp
1 /*
2 This source is published under the following 3-clause BSD license.
3
4 Copyright (c) 2012 - 2013, Lukas Hosek and Alexander Wilkie
5 All rights reserved.
6
7 Redistribution and use in source and binary forms, with or without 
8 modification, are permitted provided that the following conditions are met:
9
10     * Redistributions of source code must retain the above copyright
11       notice, this list of conditions and the following disclaimer.
12     * Redistributions in binary form must reproduce the above copyright
13       notice, this list of conditions and the following disclaimer in the
14       documentation and/or other materials provided with the distribution.
15     * None of the names of the contributors may be used to endorse or promote 
16       products derived from this software without specific prior written 
17       permission.
18
19 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
23 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31 /* ============================================================================
32
33 This file is part of a sample implementation of the analytical skylight and
34 solar radiance models presented in the SIGGRAPH 2012 paper
35
36
37            "An Analytic Model for Full Spectral Sky-Dome Radiance"
38
39 and the 2013 IEEE CG&A paper
40
41        "Adding a Solar Radiance Function to the Hosek Skylight Model"
42
43                                    both by 
44
45                        Lukas Hosek and Alexander Wilkie
46                 Charles University in Prague, Czech Republic
47
48
49                         Version: 1.4a, February 22nd, 2013
50                         
51 Version history:
52
53 1.4a  February 22nd, 2013
54       Removed unnecessary and counter-intuitive solar radius parameters 
55       from the interface of the colourspace sky dome initialisation functions.
56
57 1.4   February 11th, 2013
58       Fixed a bug which caused the relative brightness of the solar disc
59       and the sky dome to be off by a factor of about 6. The sun was too 
60       bright: this affected both normal and alien sun scenarios. The 
61       coefficients of the solar radiance function were changed to fix this.
62
63 1.3   January 21st, 2013 (not released to the public)
64       Added support for solar discs that are not exactly the same size as
65       the terrestrial sun. Also added support for suns with a different
66       emission spectrum ("Alien World" functionality).
67
68 1.2a  December 18th, 2012
69       Fixed a mistake and some inaccuracies in the solar radiance function
70       explanations found in ArHosekSkyModel.h. The actual source code is
71       unchanged compared to version 1.2.
72
73 1.2   December 17th, 2012
74       Native RGB data and a solar radiance function that matches the turbidity
75       conditions were added.
76
77 1.1   September 2012
78       The coefficients of the spectral model are now scaled so that the output
79       is given in physical units: W / (m^-2 * sr * nm). Also, the output of the
80       XYZ model is now no longer scaled to the range [0...1]. Instead, it is
81       the result of a simple conversion from spectral data via the CIE 2 degree
82       standard observer matching functions. Therefore, after multiplication
83       with 683 lm / W, the Y channel now corresponds to luminance in lm.
84      
85 1.0   May 11th, 2012
86       Initial release.
87
88
89 Please visit http://cgg.mff.cuni.cz/projects/SkylightModelling/ to check if
90 an updated version of this code has been published!
91
92 ============================================================================ */
93
94 /*
95
96 All instructions on how to use this code are in the accompanying header file.
97
98 */
99
100 #include "sky_model.h"
101 #include "sky_model_data.h"
102
103 #include <assert.h>
104 #include <stdio.h>
105 #include <stdlib.h>
106 #include <math.h>
107
108 CCL_NAMESPACE_BEGIN
109
110 //   Some macro definitions that occur elsewhere in ART, and that have to be
111 //   replicated to make this a stand-alone module.
112
113 #ifndef NIL
114 #define NIL                         0
115 #endif
116
117 #ifndef MATH_PI 
118 #define MATH_PI                     3.141592653589793
119 #endif
120
121 #ifndef MATH_DEG_TO_RAD
122 #define MATH_DEG_TO_RAD             ( MATH_PI / 180.0 )
123 #endif
124
125 #ifndef MATH_RAD_TO_DEG
126 #define MATH_RAD_TO_DEG             ( 180.0 / MATH_PI )
127 #endif
128
129 #ifndef DEGREES
130 #define DEGREES                     * MATH_DEG_TO_RAD
131 #endif
132
133 #ifndef TERRESTRIAL_SOLAR_RADIUS
134 #define TERRESTRIAL_SOLAR_RADIUS    ( ( 0.51 DEGREES ) / 2.0 )
135 #endif
136
137 #ifndef ALLOC
138 #define ALLOC(_struct)              ((_struct *)malloc(sizeof(_struct)))
139 #endif
140
141 // internal definitions
142
143 typedef double *ArHosekSkyModel_Dataset;
144 typedef double *ArHosekSkyModel_Radiance_Dataset;
145
146 // internal functions
147
148 void ArHosekSkyModel_CookConfiguration(
149         ArHosekSkyModel_Dataset       dataset, 
150         ArHosekSkyModelConfiguration  config, 
151         double                        turbidity, 
152         double                        albedo, 
153         double                        solar_elevation
154         )
155 {
156     double  * elev_matrix;
157
158     int     int_turbidity = (int)turbidity;
159     double  turbidity_rem = turbidity - (double)int_turbidity;
160
161     solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
162
163     // alb 0 low turb
164
165     elev_matrix = dataset + ( 9 * 6 * (int_turbidity-1) );
166     
167     
168     for( unsigned int i = 0; i < 9; ++i )
169     {
170         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
171         config[i] = 
172         (1.0-albedo) * (1.0 - turbidity_rem) 
173         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
174            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
175            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
176            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
177            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
178            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
179     }
180
181     // alb 1 low turb
182     elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity-1));
183     for(unsigned int i = 0; i < 9; ++i)
184     {
185         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
186         config[i] += 
187         (albedo) * (1.0 - turbidity_rem)
188         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
189            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
190            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
191            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
192            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
193            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
194     }
195
196     if(int_turbidity == 10)
197         return;
198
199     // alb 0 high turb
200     elev_matrix = dataset + (9*6*(int_turbidity));
201     for(unsigned int i = 0; i < 9; ++i)
202     {
203         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
204         config[i] += 
205         (1.0-albedo) * (turbidity_rem)
206         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
207            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
208            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
209            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
210            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
211            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
212     }
213
214     // alb 1 high turb
215     elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity));
216     for(unsigned int i = 0; i < 9; ++i)
217     {
218         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
219         config[i] += 
220         (albedo) * (turbidity_rem)
221         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
222            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
223            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
224            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
225            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
226            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
227     }
228 }
229
230 double ArHosekSkyModel_CookRadianceConfiguration(
231         ArHosekSkyModel_Radiance_Dataset  dataset, 
232         double                            turbidity, 
233         double                            albedo, 
234         double                            solar_elevation
235         )
236 {
237     double* elev_matrix;
238
239     int int_turbidity = (int)turbidity;
240     double turbidity_rem = turbidity - (double)int_turbidity;
241     double res;
242     solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
243
244     // alb 0 low turb
245     elev_matrix = dataset + (6*(int_turbidity-1));
246     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
247     res = (1.0-albedo) * (1.0 - turbidity_rem) *
248         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
249          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
250          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
251          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
252          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
253          pow(solar_elevation, 5.0) * elev_matrix[5]);
254
255     // alb 1 low turb
256     elev_matrix = dataset + (6*10 + 6*(int_turbidity-1));
257     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
258     res += (albedo) * (1.0 - turbidity_rem) *
259         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
260          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
261          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
262          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
263          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
264          pow(solar_elevation, 5.0) * elev_matrix[5]);
265     if(int_turbidity == 10)
266         return res;
267
268     // alb 0 high turb
269     elev_matrix = dataset + (6*(int_turbidity));
270     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
271     res += (1.0-albedo) * (turbidity_rem) *
272         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
273          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
274          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
275          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
276          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
277          pow(solar_elevation, 5.0) * elev_matrix[5]);
278
279     // alb 1 high turb
280     elev_matrix = dataset + (6*10 + 6*(int_turbidity));
281     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
282     res += (albedo) * (turbidity_rem) *
283         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
284          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
285          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
286          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
287          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
288          pow(solar_elevation, 5.0) * elev_matrix[5]);
289     return res;
290 }
291
292 double ArHosekSkyModel_GetRadianceInternal(
293         ArHosekSkyModelConfiguration  configuration, 
294         double                        theta, 
295         double                        gamma
296         )
297 {
298     const double expM = exp(configuration[4] * gamma);
299     const double rayM = cos(gamma)*cos(gamma);
300     const double mieM = (1.0 + cos(gamma)*cos(gamma)) / pow((1.0 + configuration[8]*configuration[8] - 2.0*configuration[8]*cos(gamma)), 1.5);
301     const double zenith = sqrt(cos(theta));
302
303     return (1.0 + configuration[0] * exp(configuration[1] / (cos(theta) + 0.01))) *
304             (configuration[2] + configuration[3] * expM + configuration[5] * rayM + configuration[6] * mieM + configuration[7] * zenith);
305 }
306
307 void arhosekskymodelstate_free(
308         ArHosekSkyModelState  * state
309         )
310 {
311     free(state);
312 }
313
314 double arhosekskymodel_radiance(
315         ArHosekSkyModelState  * state,
316         double                  theta, 
317         double                  gamma, 
318         double                  wavelength
319         )
320 {
321     int low_wl = (wavelength - 320.0 ) / 40.0;
322
323     if ( low_wl < 0 || low_wl >= 11 )
324         return 0.0f;
325
326     double interp = fmod((wavelength - 320.0 ) / 40.0, 1.0);
327
328     double val_low = 
329           ArHosekSkyModel_GetRadianceInternal(
330                 state->configs[low_wl],
331                 theta,
332                 gamma
333               )
334         * state->radiances[low_wl]
335         * state->emission_correction_factor_sky[low_wl];
336
337     if ( interp < 1e-6 )
338         return val_low;
339
340     double result = ( 1.0 - interp ) * val_low;
341
342     if ( low_wl+1 < 11 )
343     {
344         result +=
345               interp
346             * ArHosekSkyModel_GetRadianceInternal(
347                     state->configs[low_wl+1],
348                     theta,
349                     gamma
350                   )
351             * state->radiances[low_wl+1]
352             * state->emission_correction_factor_sky[low_wl+1];
353     }
354
355     return result;
356 }
357
358
359 // xyz and rgb versions
360
361 ArHosekSkyModelState  * arhosek_xyz_skymodelstate_alloc_init(
362         const double  turbidity, 
363         const double  albedo, 
364         const double  elevation
365         )
366 {
367     ArHosekSkyModelState  * state = ALLOC(ArHosekSkyModelState);
368
369     state->solar_radius = TERRESTRIAL_SOLAR_RADIUS;
370     state->turbidity    = turbidity;
371     state->albedo       = albedo;
372     state->elevation    = elevation;
373     
374     for( unsigned int channel = 0; channel < 3; ++channel )
375     {
376         ArHosekSkyModel_CookConfiguration(
377             datasetsXYZ[channel], 
378             state->configs[channel], 
379             turbidity, 
380             albedo, 
381             elevation
382             );
383         
384         state->radiances[channel] = 
385         ArHosekSkyModel_CookRadianceConfiguration(
386             datasetsXYZRad[channel],
387             turbidity, 
388             albedo,
389             elevation
390             );
391     }
392     
393     return state;
394 }
395
396 CCL_NAMESPACE_END
397