Blender  V3.3
sky_model.cpp
Go to the documentation of this file.
1 /* SPDX-License-Identifier: BSD-3-Clause
2  * Copyright 2012-2013 Lukas Hosek and Alexander Wilkie. All rights reserved. */
3 
4 /* ============================================================================
5 
6 This file is part of a sample implementation of the analytical skylight and
7 solar radiance models presented in the SIGGRAPH 2012 paper
8 
9 
10  "An Analytic Model for Full Spectral Sky-Dome Radiance"
11 
12 and the 2013 IEEE CG&A paper
13 
14  "Adding a Solar Radiance Function to the Hosek Skylight Model"
15 
16  both by
17 
18  Lukas Hosek and Alexander Wilkie
19  Charles University in Prague, Czech Republic
20 
21 
22  Version: 1.4a, February 22nd, 2013
23 
24 Version history:
25 
26 1.4a February 22nd, 2013
27  Removed unnecessary and counter-intuitive solar radius parameters
28  from the interface of the color-space sky dome initialization functions.
29 
30 1.4 February 11th, 2013
31  Fixed a bug which caused the relative brightness of the solar disc
32  and the sky dome to be off by a factor of about 6. The sun was too
33  bright: this affected both normal and alien sun scenarios. The
34  coefficients of the solar radiance function were changed to fix this.
35 
36 1.3 January 21st, 2013 (not released to the public)
37  Added support for solar discs that are not exactly the same size as
38  the terrestrial sun. Also added support for suns with a different
39  emission spectrum ("Alien World" functionality).
40 
41 1.2a December 18th, 2012
42  Fixed a mistake and some inaccuracies in the solar radiance function
43  explanations found in ArHosekSkyModel.h. The actual source code is
44  unchanged compared to version 1.2.
45 
46 1.2 December 17th, 2012
47  Native RGB data and a solar radiance function that matches the turbidity
48  conditions were added.
49 
50 1.1 September 2012
51  The coefficients of the spectral model are now scaled so that the output
52  is given in physical units: W / (m^-2 * sr * nm). Also, the output of the
53  XYZ model is now no longer scaled to the range [0...1]. Instead, it is
54  the result of a simple conversion from spectral data via the CIE 2 degree
55  standard observer matching functions. Therefore, after multiplication
56  with 683 lm / W, the Y channel now corresponds to luminance in lm.
57 
58 1.0 May 11th, 2012
59  Initial release.
60 
61 
62 Please visit http://cgg.mff.cuni.cz/projects/SkylightModelling/ to check if
63 an updated version of this code has been published!
64 
65 ============================================================================ */
66 
67 /*
68 
69 All instructions on how to use this code are in the accompanying header file.
70 
71 */
72 
77 #include "sky_model.h"
78 #include "sky_model_data.h"
79 
80 #include <assert.h>
81 #include <math.h>
82 #include <stdio.h>
83 #include <stdlib.h>
84 
85 // Some macro definitions that occur elsewhere in ART, and that have to be
86 // replicated to make this a stand-alone module.
87 
88 #ifndef MATH_PI
89 # define MATH_PI 3.141592653589793
90 #endif
91 
92 #ifndef MATH_DEG_TO_RAD
93 # define MATH_DEG_TO_RAD (MATH_PI / 180.0)
94 #endif
95 
96 #ifndef DEGREES
97 # define DEGREES *MATH_DEG_TO_RAD
98 #endif
99 
100 #ifndef TERRESTRIAL_SOLAR_RADIUS
101 # define TERRESTRIAL_SOLAR_RADIUS ((0.51 DEGREES) / 2.0)
102 #endif
103 
104 #ifndef ALLOC
105 # define ALLOC(_struct) ((_struct *)malloc(sizeof(_struct)))
106 #endif
107 
108 // internal definitions
109 
110 typedef const double *ArHosekSkyModel_Dataset;
111 typedef const double *ArHosekSkyModel_Radiance_Dataset;
112 
113 // internal functions
114 
117  double turbidity,
118  double albedo,
119  double solar_elevation)
120 {
121  const double *elev_matrix;
122 
123  int int_turbidity = (int)turbidity;
124  double turbidity_rem = turbidity - (double)int_turbidity;
125 
126  solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
127 
128  // alb 0 low turb
129 
130  elev_matrix = dataset + (9 * 6 * (int_turbidity - 1));
131 
132  for (unsigned int i = 0; i < 9; ++i) {
133  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
134  config[i] =
135  (1.0 - albedo) * (1.0 - turbidity_rem) *
136  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[i] +
137  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[i + 9] +
138  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[i + 18] +
139  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[i + 27] +
140  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[i + 36] +
141  pow(solar_elevation, 5.0) * elev_matrix[i + 45]);
142  }
143 
144  // alb 1 low turb
145  elev_matrix = dataset + (9 * 6 * 10 + 9 * 6 * (int_turbidity - 1));
146  for (unsigned int i = 0; i < 9; ++i) {
147  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
148  config[i] +=
149  (albedo) * (1.0 - turbidity_rem) *
150  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[i] +
151  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[i + 9] +
152  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[i + 18] +
153  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[i + 27] +
154  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[i + 36] +
155  pow(solar_elevation, 5.0) * elev_matrix[i + 45]);
156  }
157 
158  if (int_turbidity == 10) {
159  return;
160  }
161 
162  // alb 0 high turb
163  elev_matrix = dataset + (9 * 6 * (int_turbidity));
164  for (unsigned int i = 0; i < 9; ++i) {
165  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
166  config[i] +=
167  (1.0 - albedo) * (turbidity_rem) *
168  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[i] +
169  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[i + 9] +
170  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[i + 18] +
171  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[i + 27] +
172  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[i + 36] +
173  pow(solar_elevation, 5.0) * elev_matrix[i + 45]);
174  }
175 
176  // alb 1 high turb
177  elev_matrix = dataset + (9 * 6 * 10 + 9 * 6 * (int_turbidity));
178  for (unsigned int i = 0; i < 9; ++i) {
179  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
180  config[i] +=
181  (albedo) * (turbidity_rem) *
182  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[i] +
183  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[i + 9] +
184  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[i + 18] +
185  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[i + 27] +
186  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[i + 36] +
187  pow(solar_elevation, 5.0) * elev_matrix[i + 45]);
188  }
189 }
190 
192  double turbidity,
193  double albedo,
194  double solar_elevation)
195 {
196  const double *elev_matrix;
197 
198  int int_turbidity = (int)turbidity;
199  double turbidity_rem = turbidity - (double)int_turbidity;
200  double res;
201  solar_elevation = pow(solar_elevation / (MATH_PI / 2.0), (1.0 / 3.0));
202 
203  // alb 0 low turb
204  elev_matrix = dataset + (6 * (int_turbidity - 1));
205  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
206  res = (1.0 - albedo) * (1.0 - turbidity_rem) *
207  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[0] +
208  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[1] +
209  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[2] +
210  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[3] +
211  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[4] +
212  pow(solar_elevation, 5.0) * elev_matrix[5]);
213 
214  // alb 1 low turb
215  elev_matrix = dataset + (6 * 10 + 6 * (int_turbidity - 1));
216  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
217  res += (albedo) * (1.0 - turbidity_rem) *
218  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[0] +
219  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[1] +
220  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[2] +
221  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[3] +
222  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[4] +
223  pow(solar_elevation, 5.0) * elev_matrix[5]);
224  if (int_turbidity == 10) {
225  return res;
226  }
227 
228  // alb 0 high turb
229  elev_matrix = dataset + (6 * (int_turbidity));
230  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
231  res += (1.0 - albedo) * (turbidity_rem) *
232  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[0] +
233  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[1] +
234  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[2] +
235  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[3] +
236  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[4] +
237  pow(solar_elevation, 5.0) * elev_matrix[5]);
238 
239  // alb 1 high turb
240  elev_matrix = dataset + (6 * 10 + 6 * (int_turbidity));
241  //(1-t).^3* A1 + 3*(1-t).^2.*t * A2 + 3*(1-t) .* t .^ 2 * A3 + t.^3 * A4;
242  res += (albedo) * (turbidity_rem) *
243  (pow(1.0 - solar_elevation, 5.0) * elev_matrix[0] +
244  5.0 * pow(1.0 - solar_elevation, 4.0) * solar_elevation * elev_matrix[1] +
245  10.0 * pow(1.0 - solar_elevation, 3.0) * pow(solar_elevation, 2.0) * elev_matrix[2] +
246  10.0 * pow(1.0 - solar_elevation, 2.0) * pow(solar_elevation, 3.0) * elev_matrix[3] +
247  5.0 * (1.0 - solar_elevation) * pow(solar_elevation, 4.0) * elev_matrix[4] +
248  pow(solar_elevation, 5.0) * elev_matrix[5]);
249  return res;
250 }
251 
253  const SKY_ArHosekSkyModelConfiguration configuration, const double theta, const double gamma)
254 {
255  const double expM = exp(configuration[4] * gamma);
256  const double rayM = cos(gamma) * cos(gamma);
257  const double mieM =
258  (1.0 + cos(gamma) * cos(gamma)) /
259  pow((1.0 + configuration[8] * configuration[8] - 2.0 * configuration[8] * cos(gamma)), 1.5);
260  const double zenith = sqrt(cos(theta));
261 
262  return (1.0 + configuration[0] * exp(configuration[1] / (cos(theta) + 0.01))) *
263  (configuration[2] + configuration[3] * expM + configuration[5] * rayM +
264  configuration[6] * mieM + configuration[7] * zenith);
265 }
266 
268 {
269  free(state);
270 }
271 
273  double theta,
274  double gamma,
275  double wavelength)
276 {
277  int low_wl = (int)((wavelength - 320.0) / 40.0);
278 
279  if (low_wl < 0 || low_wl >= 11) {
280  return 0.0;
281  }
282 
283  double interp = fmod((wavelength - 320.0) / 40.0, 1.0);
284 
285  double val_low = ArHosekSkyModel_GetRadianceInternal(state->configs[low_wl], theta, gamma) *
286  state->radiances[low_wl] * state->emission_correction_factor_sky[low_wl];
287 
288  if (interp < 1e-6) {
289  return val_low;
290  }
291 
292  double result = (1.0 - interp) * val_low;
293 
294  if (low_wl + 1 < 11) {
295  result += interp *
296  ArHosekSkyModel_GetRadianceInternal(state->configs[low_wl + 1], theta, gamma) *
297  state->radiances[low_wl + 1] * state->emission_correction_factor_sky[low_wl + 1];
298  }
299 
300  return result;
301 }
302 
303 // xyz and rgb versions
304 
306  const double albedo,
307  const double elevation)
308 {
310 
311  state->solar_radius = TERRESTRIAL_SOLAR_RADIUS;
312  state->turbidity = turbidity;
313  state->albedo = albedo;
314  state->elevation = elevation;
315 
316  for (unsigned int channel = 0; channel < 3; ++channel) {
318  datasetsXYZ[channel], state->configs[channel], turbidity, albedo, elevation);
319 
320  state->radiances[channel] = ArHosekSkyModel_CookRadianceConfiguration(
321  datasetsXYZRad[channel], turbidity, albedo, elevation);
322  }
323 
324  return state;
325 }
sqrt(x)+1/max(0
void BLI_kdtree_nd_() free(KDTree *tree)
Definition: kdtree_impl.h:102
typedef double(DMatrix)[4][4]
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
const int state
ccl_device_inline float2 interp(const float2 &a, const float2 &b, float t)
Definition: math_float2.h:232
ccl_device_inline float3 exp(float3 v)
Definition: math_float3.h:392
ccl_device_inline float3 pow(float3 v, float e)
Definition: math_float3.h:533
INLINE Rall1d< T, V, S > cos(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:319
const double * ArHosekSkyModel_Dataset
Definition: sky_model.cpp:110
static void ArHosekSkyModel_CookConfiguration(ArHosekSkyModel_Dataset dataset, SKY_ArHosekSkyModelConfiguration config, double turbidity, double albedo, double solar_elevation)
Definition: sky_model.cpp:115
void SKY_arhosekskymodelstate_free(SKY_ArHosekSkyModelState *state)
Definition: sky_model.cpp:267
#define TERRESTRIAL_SOLAR_RADIUS
Definition: sky_model.cpp:101
static double ArHosekSkyModel_CookRadianceConfiguration(ArHosekSkyModel_Radiance_Dataset dataset, double turbidity, double albedo, double solar_elevation)
Definition: sky_model.cpp:191
#define ALLOC(_struct)
Definition: sky_model.cpp:105
double SKY_arhosekskymodel_radiance(SKY_ArHosekSkyModelState *state, double theta, double gamma, double wavelength)
Definition: sky_model.cpp:272
#define MATH_PI
Definition: sky_model.cpp:89
const double * ArHosekSkyModel_Radiance_Dataset
Definition: sky_model.cpp:111
static double ArHosekSkyModel_GetRadianceInternal(const SKY_ArHosekSkyModelConfiguration configuration, const double theta, const double gamma)
Definition: sky_model.cpp:252
SKY_ArHosekSkyModelState * SKY_arhosek_xyz_skymodelstate_alloc_init(const double turbidity, const double albedo, const double elevation)
Definition: sky_model.cpp:305
double SKY_ArHosekSkyModelConfiguration[9]
Definition: sky_model.h:285
static const double * datasetsXYZ[]
static const double * datasetsXYZRad[]