Cycles: Make all #include statements relative to cycles source directory
[blender.git] / intern / cycles / util / util_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 "util/util_sky_model.h"
101 #include "util/util_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 MATH_PI
114 #define MATH_PI                     3.141592653589793
115 #endif
116
117 #ifndef MATH_DEG_TO_RAD
118 #define MATH_DEG_TO_RAD             ( MATH_PI / 180.0 )
119 #endif
120
121 #ifndef DEGREES
122 #define DEGREES                     * MATH_DEG_TO_RAD
123 #endif
124
125 #ifndef TERRESTRIAL_SOLAR_RADIUS
126 #define TERRESTRIAL_SOLAR_RADIUS    ( ( 0.51 DEGREES ) / 2.0 )
127 #endif
128
129 #ifndef ALLOC
130 #define ALLOC(_struct)              ((_struct *)malloc(sizeof(_struct)))
131 #endif
132
133 // internal definitions
134
135 typedef const double *ArHosekSkyModel_Dataset;
136 typedef const double *ArHosekSkyModel_Radiance_Dataset;
137
138 // internal functions
139
140 static void ArHosekSkyModel_CookConfiguration(
141         ArHosekSkyModel_Dataset dataset,
142         ArHosekSkyModelConfiguration config,
143         double turbidity,
144         double albedo,
145         double solar_elevation)
146 {
147         const double  * elev_matrix;
148
149         int int_turbidity = (int)turbidity;
150         double turbidity_rem = turbidity - (double)int_turbidity;
151
152         solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
153
154         // alb 0 low turb
155
156         elev_matrix = dataset + ( 9 * 6 * (int_turbidity-1));
157
158         for(unsigned int i = 0; i < 9; ++i) {
159                 //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
160                 config[i] =
161                         (1.0-albedo) * (1.0 - turbidity_rem)
162                         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  +
163                             5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
164                             10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
165                             10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
166                             5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
167                             pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
168         }
169
170         // alb 1 low turb
171         elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity-1));
172         for(unsigned int i = 0; i < 9; ++i) {
173                 //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
174                 config[i] +=
175                         (albedo) * (1.0 - turbidity_rem)
176                         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  +
177                             5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
178                             10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
179                             10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
180                             5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
181                             pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
182         }
183
184         if(int_turbidity == 10)
185                 return;
186
187         // alb 0 high turb
188         elev_matrix = dataset + (9*6*(int_turbidity));
189         for(unsigned int i = 0; i < 9; ++i) {
190                 //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
191                 config[i] +=
192                         (1.0-albedo) * (turbidity_rem)
193                         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  +
194                             5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
195                             10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
196                             10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
197                             5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
198                             pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
199         }
200
201         // alb 1 high turb
202         elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity));
203         for(unsigned int i = 0; i < 9; ++i) {
204                 //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
205                 config[i] +=
206                         (albedo) * (turbidity_rem)
207                         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  +
208                             5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
209                             10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
210                             10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
211                             5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
212                             pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
213         }
214 }
215
216 static double ArHosekSkyModel_CookRadianceConfiguration(
217         ArHosekSkyModel_Radiance_Dataset dataset,
218         double turbidity,
219         double albedo,
220         double solar_elevation)
221 {
222         const double* elev_matrix;
223
224         int int_turbidity = (int)turbidity;
225         double turbidity_rem = turbidity - (double)int_turbidity;
226         double res;
227         solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
228
229         // alb 0 low turb
230         elev_matrix = dataset + (6*(int_turbidity-1));
231         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
232         res = (1.0-albedo) * (1.0 - turbidity_rem) *
233                 ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
234                   5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
235                   10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
236                   10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
237                   5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
238                   pow(solar_elevation, 5.0) * elev_matrix[5]);
239
240         // alb 1 low turb
241         elev_matrix = dataset + (6*10 + 6*(int_turbidity-1));
242         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
243         res += (albedo) * (1.0 - turbidity_rem) *
244                 ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
245                   5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
246                   10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
247                   10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
248                   5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
249                   pow(solar_elevation, 5.0) * elev_matrix[5]);
250         if(int_turbidity == 10)
251                 return res;
252
253         // alb 0 high turb
254         elev_matrix = dataset + (6*(int_turbidity));
255         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
256         res += (1.0-albedo) * (turbidity_rem) *
257                 ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
258                   5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
259                   10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
260                   10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
261                   5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
262                   pow(solar_elevation, 5.0) * elev_matrix[5]);
263
264         // alb 1 high turb
265         elev_matrix = dataset + (6*10 + 6*(int_turbidity));
266         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
267         res += (albedo) * (turbidity_rem) *
268                 ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
269                   5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
270                   10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
271                   10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
272                   5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
273                   pow(solar_elevation, 5.0) * elev_matrix[5]);
274         return res;
275 }
276
277 static double ArHosekSkyModel_GetRadianceInternal(
278         ArHosekSkyModelConfiguration configuration,
279         double theta,
280         double gamma)
281 {
282         const double expM = exp(configuration[4] * gamma);
283         const double rayM = cos(gamma)*cos(gamma);
284         const double mieM = (1.0 + cos(gamma)*cos(gamma)) / pow((1.0 + configuration[8]*configuration[8] - 2.0*configuration[8]*cos(gamma)), 1.5);
285         const double zenith = sqrt(cos(theta));
286
287         return (1.0 + configuration[0] * exp(configuration[1] / (cos(theta) + 0.01))) *
288             (configuration[2] + configuration[3] * expM + configuration[5] * rayM + configuration[6] * mieM + configuration[7] * zenith);
289 }
290
291 void arhosekskymodelstate_free(ArHosekSkyModelState  * state)
292 {
293         free(state);
294 }
295
296 double arhosekskymodel_radiance(ArHosekSkyModelState  *state,
297                                 double theta,
298                                 double gamma,
299                                 double wavelength)
300 {
301         int low_wl = (int)((wavelength - 320.0) / 40.0);
302
303         if(low_wl < 0 || low_wl >= 11)
304                 return 0.0;
305
306         double interp = fmod((wavelength - 320.0 ) / 40.0, 1.0);
307
308         double val_low =
309                 ArHosekSkyModel_GetRadianceInternal(
310                         state->configs[low_wl],
311                         theta,
312                         gamma)
313                 * state->radiances[low_wl]
314                 * state->emission_correction_factor_sky[low_wl];
315
316         if(interp < 1e-6)
317                 return val_low;
318
319         double result = ( 1.0 - interp ) * val_low;
320
321         if(low_wl+1 < 11) {
322                 result +=
323                     interp
324                     * ArHosekSkyModel_GetRadianceInternal(
325                             state->configs[low_wl+1],
326                             theta,
327                             gamma)
328                     * state->radiances[low_wl+1]
329                     * state->emission_correction_factor_sky[low_wl+1];
330         }
331
332         return result;
333 }
334
335
336 // xyz and rgb versions
337
338 ArHosekSkyModelState * arhosek_xyz_skymodelstate_alloc_init(
339         const double turbidity,
340         const double albedo,
341         const double elevation)
342 {
343         ArHosekSkyModelState  * state = ALLOC(ArHosekSkyModelState);
344
345         state->solar_radius = TERRESTRIAL_SOLAR_RADIUS;
346         state->turbidity = turbidity;
347         state->albedo = albedo;
348         state->elevation = elevation;
349
350         for(unsigned int channel = 0; channel < 3; ++channel) {
351                 ArHosekSkyModel_CookConfiguration(
352                     datasetsXYZ[channel],
353                     state->configs[channel],
354                     turbidity,
355                     albedo,
356                     elevation);
357
358                 state->radiances[channel] =
359                 ArHosekSkyModel_CookRadianceConfiguration(
360                     datasetsXYZRad[channel],
361                     turbidity,
362                     albedo,
363                     elevation);
364     }
365
366     return state;
367 }
368
369 CCL_NAMESPACE_END
370