Cycles: Code cleanup, spaces around keywords
[blender-staging.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 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 {
148     const double  * elev_matrix;
149
150     int     int_turbidity = (int)turbidity;
151     double  turbidity_rem = turbidity - (double)int_turbidity;
152
153     solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
154
155     // alb 0 low turb
156
157     elev_matrix = dataset + ( 9 * 6 * (int_turbidity-1) );
158     
159     
160     for( unsigned int i = 0; i < 9; ++i )
161     {
162         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
163         config[i] = 
164         (1.0-albedo) * (1.0 - turbidity_rem) 
165         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
166            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
167            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
168            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
169            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
170            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
171     }
172
173     // alb 1 low turb
174     elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity-1));
175     for(unsigned int i = 0; i < 9; ++i)
176     {
177         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
178         config[i] += 
179         (albedo) * (1.0 - turbidity_rem)
180         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
181            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
182            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
183            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
184            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
185            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
186     }
187
188     if(int_turbidity == 10)
189         return;
190
191     // alb 0 high turb
192     elev_matrix = dataset + (9*6*(int_turbidity));
193     for(unsigned int i = 0; i < 9; ++i)
194     {
195         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
196         config[i] += 
197         (1.0-albedo) * (turbidity_rem)
198         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
199            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
200            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
201            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
202            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
203            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
204     }
205
206     // alb 1 high turb
207     elev_matrix = dataset + (9*6*10 + 9*6*(int_turbidity));
208     for(unsigned int i = 0; i < 9; ++i)
209     {
210         //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
211         config[i] += 
212         (albedo) * (turbidity_rem)
213         * ( pow(1.0-solar_elevation, 5.0) * elev_matrix[i]  + 
214            5.0  * pow(1.0-solar_elevation, 4.0) * solar_elevation * elev_matrix[i+9] +
215            10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[i+18] +
216            10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[i+27] +
217            5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[i+36] +
218            pow(solar_elevation, 5.0)  * elev_matrix[i+45]);
219     }
220 }
221
222 static double ArHosekSkyModel_CookRadianceConfiguration(
223         ArHosekSkyModel_Radiance_Dataset  dataset, 
224         double                            turbidity, 
225         double                            albedo, 
226         double                            solar_elevation
227         )
228 {
229     const double* elev_matrix;
230
231     int int_turbidity = (int)turbidity;
232     double turbidity_rem = turbidity - (double)int_turbidity;
233     double res;
234     solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
235
236     // alb 0 low turb
237     elev_matrix = dataset + (6*(int_turbidity-1));
238     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
239     res = (1.0-albedo) * (1.0 - turbidity_rem) *
240         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
241          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
242          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
243          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
244          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
245          pow(solar_elevation, 5.0) * elev_matrix[5]);
246
247     // alb 1 low turb
248     elev_matrix = dataset + (6*10 + 6*(int_turbidity-1));
249     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
250     res += (albedo) * (1.0 - turbidity_rem) *
251         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
252          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
253          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
254          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
255          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
256          pow(solar_elevation, 5.0) * elev_matrix[5]);
257     if(int_turbidity == 10)
258         return res;
259
260     // alb 0 high turb
261     elev_matrix = dataset + (6*(int_turbidity));
262     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
263     res += (1.0-albedo) * (turbidity_rem) *
264         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
265          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
266          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
267          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
268          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
269          pow(solar_elevation, 5.0) * elev_matrix[5]);
270
271     // alb 1 high turb
272     elev_matrix = dataset + (6*10 + 6*(int_turbidity));
273     //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
274     res += (albedo) * (turbidity_rem) *
275         ( pow(1.0-solar_elevation, 5.0) * elev_matrix[0] +
276          5.0*pow(1.0-solar_elevation, 4.0)*solar_elevation * elev_matrix[1] +
277          10.0*pow(1.0-solar_elevation, 3.0)*pow(solar_elevation, 2.0) * elev_matrix[2] +
278          10.0*pow(1.0-solar_elevation, 2.0)*pow(solar_elevation, 3.0) * elev_matrix[3] +
279          5.0*(1.0-solar_elevation)*pow(solar_elevation, 4.0) * elev_matrix[4] +
280          pow(solar_elevation, 5.0) * elev_matrix[5]);
281     return res;
282 }
283
284 static double ArHosekSkyModel_GetRadianceInternal(
285         ArHosekSkyModelConfiguration  configuration, 
286         double                        theta, 
287         double                        gamma
288         )
289 {
290     const double expM = exp(configuration[4] * gamma);
291     const double rayM = cos(gamma)*cos(gamma);
292     const double mieM = (1.0 + cos(gamma)*cos(gamma)) / pow((1.0 + configuration[8]*configuration[8] - 2.0*configuration[8]*cos(gamma)), 1.5);
293     const double zenith = sqrt(cos(theta));
294
295     return (1.0 + configuration[0] * exp(configuration[1] / (cos(theta) + 0.01))) *
296             (configuration[2] + configuration[3] * expM + configuration[5] * rayM + configuration[6] * mieM + configuration[7] * zenith);
297 }
298
299 void arhosekskymodelstate_free(
300         ArHosekSkyModelState  * state
301         )
302 {
303     free(state);
304 }
305
306 double arhosekskymodel_radiance(
307         ArHosekSkyModelState  * state,
308         double                  theta, 
309         double                  gamma, 
310         double                  wavelength
311         )
312 {
313     int low_wl = (int)((wavelength - 320.0) / 40.0);
314
315     if( low_wl < 0 || low_wl >= 11 )
316         return 0.0f;
317
318     double interp = fmod((wavelength - 320.0 ) / 40.0, 1.0);
319
320     double val_low = 
321           ArHosekSkyModel_GetRadianceInternal(
322                 state->configs[low_wl],
323                 theta,
324                 gamma
325               )
326         * state->radiances[low_wl]
327         * state->emission_correction_factor_sky[low_wl];
328
329     if( interp < 1e-6 )
330         return val_low;
331
332     double result = ( 1.0 - interp ) * val_low;
333
334     if( low_wl+1 < 11 )
335     {
336         result +=
337               interp
338             * ArHosekSkyModel_GetRadianceInternal(
339                     state->configs[low_wl+1],
340                     theta,
341                     gamma
342                   )
343             * state->radiances[low_wl+1]
344             * state->emission_correction_factor_sky[low_wl+1];
345     }
346
347     return result;
348 }
349
350
351 // xyz and rgb versions
352
353 ArHosekSkyModelState  * arhosek_xyz_skymodelstate_alloc_init(
354         const double  turbidity, 
355         const double  albedo, 
356         const double  elevation
357         )
358 {
359     ArHosekSkyModelState  * state = ALLOC(ArHosekSkyModelState);
360
361     state->solar_radius = TERRESTRIAL_SOLAR_RADIUS;
362     state->turbidity    = turbidity;
363     state->albedo       = albedo;
364     state->elevation    = elevation;
365     
366     for( unsigned int channel = 0; channel < 3; ++channel )
367     {
368         ArHosekSkyModel_CookConfiguration(
369             datasetsXYZ[channel], 
370             state->configs[channel], 
371             turbidity, 
372             albedo, 
373             elevation
374             );
375         
376         state->radiances[channel] = 
377         ArHosekSkyModel_CookRadianceConfiguration(
378             datasetsXYZRad[channel],
379             turbidity, 
380             albedo,
381             elevation
382             );
383     }
384     
385     return state;
386 }
387
388 CCL_NAMESPACE_END
389