diff options
24 files changed, 1847 insertions, 705 deletions
diff --git a/intern/cycles/blender/blender_shader.cpp b/intern/cycles/blender/blender_shader.cpp index f207d8ae07f..19d2730dc93 100644 --- a/intern/cycles/blender/blender_shader.cpp +++ b/intern/cycles/blender/blender_shader.cpp @@ -813,6 +813,14 @@ static ShaderNode *add_node(Scene *scene, sky->sun_direction = normalize(get_float3(b_sky_node.sun_direction())); sky->turbidity = b_sky_node.turbidity(); sky->ground_albedo = b_sky_node.ground_albedo(); + sky->sun_disc = b_sky_node.sun_disc(); + sky->sun_size = b_sky_node.sun_size(); + sky->sun_elevation = b_sky_node.sun_elevation(); + sky->sun_rotation = b_sky_node.sun_rotation(); + sky->altitude = b_sky_node.altitude(); + sky->air_density = b_sky_node.air_density(); + sky->dust_density = b_sky_node.dust_density(); + sky->ozone_density = b_sky_node.ozone_density(); BL::TexMapping b_texture_mapping(b_sky_node.texture_mapping()); get_tex_mapping(&sky->tex_mapping, b_texture_mapping); node = sky; diff --git a/intern/cycles/kernel/CMakeLists.txt b/intern/cycles/kernel/CMakeLists.txt index 2e839a616e9..35339abff45 100644 --- a/intern/cycles/kernel/CMakeLists.txt +++ b/intern/cycles/kernel/CMakeLists.txt @@ -113,6 +113,8 @@ set(SRC_HEADERS kernel_id_passes.h kernel_jitter.h kernel_light.h + kernel_light_background.h + kernel_light_common.h kernel_math.h kernel_montecarlo.h kernel_passes.h diff --git a/intern/cycles/kernel/kernel_emission.h b/intern/cycles/kernel/kernel_emission.h index 71b176a0a8f..4ac07d86dda 100644 --- a/intern/cycles/kernel/kernel_emission.h +++ b/intern/cycles/kernel/kernel_emission.h @@ -326,9 +326,7 @@ ccl_device_noinline_cpu float3 indirect_background(KernelGlobals *kg, /* Background MIS weights. */ # ifdef __BACKGROUND_MIS__ /* Check if background light exists or if we should skip pdf. */ - int res_x = kernel_data.integrator.pdf_background_res_x; - - if (!(state->flag & PATH_RAY_MIS_SKIP) && res_x) { + if (!(state->flag & PATH_RAY_MIS_SKIP) && kernel_data.background.use_mis) { /* multiple importance sampling, get background light pdf for ray * direction, and compute weight with respect to BSDF pdf */ float pdf = background_light_pdf(kg, ray->P, ray->D); diff --git a/intern/cycles/kernel/kernel_light.h b/intern/cycles/kernel/kernel_light.h index 04472212d0c..0448d0165b9 100644 --- a/intern/cycles/kernel/kernel_light.h +++ b/intern/cycles/kernel/kernel_light.h @@ -14,6 +14,8 @@ * limitations under the License. */ +#include "kernel_light_background.h" + CCL_NAMESPACE_BEGIN /* Light Sample result */ @@ -33,500 +35,6 @@ typedef struct LightSample { LightType type; /* type of light */ } LightSample; -/* Area light sampling */ - -/* Uses the following paper: - * - * Carlos Urena et al. - * An Area-Preserving Parametrization for Spherical Rectangles. - * - * https://www.solidangle.com/research/egsr2013_spherical_rectangle.pdf - * - * Note: light_p is modified when sample_coord is true. - */ -ccl_device_inline float rect_light_sample(float3 P, - float3 *light_p, - float3 axisu, - float3 axisv, - float randu, - float randv, - bool sample_coord) -{ - /* In our name system we're using P for the center, - * which is o in the paper. - */ - - float3 corner = *light_p - axisu * 0.5f - axisv * 0.5f; - float axisu_len, axisv_len; - /* Compute local reference system R. */ - float3 x = normalize_len(axisu, &axisu_len); - float3 y = normalize_len(axisv, &axisv_len); - float3 z = cross(x, y); - /* Compute rectangle coords in local reference system. */ - float3 dir = corner - P; - float z0 = dot(dir, z); - /* Flip 'z' to make it point against Q. */ - if (z0 > 0.0f) { - z *= -1.0f; - z0 *= -1.0f; - } - float x0 = dot(dir, x); - float y0 = dot(dir, y); - float x1 = x0 + axisu_len; - float y1 = y0 + axisv_len; - /* Compute internal angles (gamma_i). */ - float4 diff = make_float4(x0, y1, x1, y0) - make_float4(x1, y0, x0, y1); - float4 nz = make_float4(y0, x1, y1, x0) * diff; - nz = nz / sqrt(z0 * z0 * diff * diff + nz * nz); - float g0 = safe_acosf(-nz.x * nz.y); - float g1 = safe_acosf(-nz.y * nz.z); - float g2 = safe_acosf(-nz.z * nz.w); - float g3 = safe_acosf(-nz.w * nz.x); - /* Compute predefined constants. */ - float b0 = nz.x; - float b1 = nz.z; - float b0sq = b0 * b0; - float k = M_2PI_F - g2 - g3; - /* Compute solid angle from internal angles. */ - float S = g0 + g1 - k; - - if (sample_coord) { - /* Compute cu. */ - float au = randu * S + k; - float fu = (cosf(au) * b0 - b1) / sinf(au); - float cu = 1.0f / sqrtf(fu * fu + b0sq) * (fu > 0.0f ? 1.0f : -1.0f); - cu = clamp(cu, -1.0f, 1.0f); - /* Compute xu. */ - float xu = -(cu * z0) / max(sqrtf(1.0f - cu * cu), 1e-7f); - xu = clamp(xu, x0, x1); - /* Compute yv. */ - float z0sq = z0 * z0; - float y0sq = y0 * y0; - float y1sq = y1 * y1; - float d = sqrtf(xu * xu + z0sq); - float h0 = y0 / sqrtf(d * d + y0sq); - float h1 = y1 / sqrtf(d * d + y1sq); - float hv = h0 + randv * (h1 - h0), hv2 = hv * hv; - float yv = (hv2 < 1.0f - 1e-6f) ? (hv * d) / sqrtf(1.0f - hv2) : y1; - - /* Transform (xu, yv, z0) to world coords. */ - *light_p = P + xu * x + yv * y + z0 * z; - } - - /* return pdf */ - if (S != 0.0f) - return 1.0f / S; - else - return 0.0f; -} - -ccl_device_inline float3 ellipse_sample(float3 ru, float3 rv, float randu, float randv) -{ - to_unit_disk(&randu, &randv); - return ru * randu + rv * randv; -} - -ccl_device float3 disk_light_sample(float3 v, float randu, float randv) -{ - float3 ru, rv; - - make_orthonormals(v, &ru, &rv); - - return ellipse_sample(ru, rv, randu, randv); -} - -ccl_device float3 distant_light_sample(float3 D, float radius, float randu, float randv) -{ - return normalize(D + disk_light_sample(D, randu, randv) * radius); -} - -ccl_device float3 -sphere_light_sample(float3 P, float3 center, float radius, float randu, float randv) -{ - return disk_light_sample(normalize(P - center), randu, randv) * radius; -} - -ccl_device float spot_light_attenuation(float3 dir, - float spot_angle, - float spot_smooth, - LightSample *ls) -{ - float3 I = ls->Ng; - - float attenuation = dot(dir, I); - - if (attenuation <= spot_angle) { - attenuation = 0.0f; - } - else { - float t = attenuation - spot_angle; - - if (t < spot_smooth && spot_smooth != 0.0f) - attenuation *= smoothstepf(t / spot_smooth); - } - - return attenuation; -} - -ccl_device float lamp_light_pdf(KernelGlobals *kg, const float3 Ng, const float3 I, float t) -{ - float cos_pi = dot(Ng, I); - - if (cos_pi <= 0.0f) - return 0.0f; - - return t * t / cos_pi; -} - -/* Background Light */ - -#ifdef __BACKGROUND_MIS__ - -ccl_device float3 background_map_sample(KernelGlobals *kg, float randu, float randv, float *pdf) -{ - /* for the following, the CDF values are actually a pair of floats, with the - * function value as X and the actual CDF as Y. The last entry's function - * value is the CDF total. */ - int res_x = kernel_data.integrator.pdf_background_res_x; - int res_y = kernel_data.integrator.pdf_background_res_y; - int cdf_width = res_x + 1; - - /* this is basically std::lower_bound as used by pbrt */ - int first = 0; - int count = res_y; - - while (count > 0) { - int step = count >> 1; - int middle = first + step; - - if (kernel_tex_fetch(__light_background_marginal_cdf, middle).y < randv) { - first = middle + 1; - count -= step + 1; - } - else - count = step; - } - - int index_v = max(0, first - 1); - kernel_assert(index_v >= 0 && index_v < res_y); - - float2 cdf_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v); - float2 cdf_next_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v + 1); - float2 cdf_last_v = kernel_tex_fetch(__light_background_marginal_cdf, res_y); - - /* importance-sampled V direction */ - float dv = inverse_lerp(cdf_v.y, cdf_next_v.y, randv); - float v = (index_v + dv) / res_y; - - /* this is basically std::lower_bound as used by pbrt */ - first = 0; - count = res_x; - while (count > 0) { - int step = count >> 1; - int middle = first + step; - - if (kernel_tex_fetch(__light_background_conditional_cdf, index_v * cdf_width + middle).y < - randu) { - first = middle + 1; - count -= step + 1; - } - else - count = step; - } - - int index_u = max(0, first - 1); - kernel_assert(index_u >= 0 && index_u < res_x); - - float2 cdf_u = kernel_tex_fetch(__light_background_conditional_cdf, - index_v * cdf_width + index_u); - float2 cdf_next_u = kernel_tex_fetch(__light_background_conditional_cdf, - index_v * cdf_width + index_u + 1); - float2 cdf_last_u = kernel_tex_fetch(__light_background_conditional_cdf, - index_v * cdf_width + res_x); - - /* importance-sampled U direction */ - float du = inverse_lerp(cdf_u.y, cdf_next_u.y, randu); - float u = (index_u + du) / res_x; - - /* compute pdf */ - float sin_theta = sinf(M_PI_F * v); - float denom = (M_2PI_F * M_PI_F * sin_theta) * cdf_last_u.x * cdf_last_v.x; - - if (sin_theta == 0.0f || denom == 0.0f) - *pdf = 0.0f; - else - *pdf = (cdf_u.x * cdf_v.x) / denom; - - /* compute direction */ - return equirectangular_to_direction(u, v); -} - -/* TODO(sergey): Same as above, after the release we should consider using - * 'noinline' for all devices. - */ -ccl_device float background_map_pdf(KernelGlobals *kg, float3 direction) -{ - float2 uv = direction_to_equirectangular(direction); - int res_x = kernel_data.integrator.pdf_background_res_x; - int res_y = kernel_data.integrator.pdf_background_res_y; - int cdf_width = res_x + 1; - - float sin_theta = sinf(uv.y * M_PI_F); - - if (sin_theta == 0.0f) - return 0.0f; - - int index_u = clamp(float_to_int(uv.x * res_x), 0, res_x - 1); - int index_v = clamp(float_to_int(uv.y * res_y), 0, res_y - 1); - - /* pdfs in V direction */ - float2 cdf_last_u = kernel_tex_fetch(__light_background_conditional_cdf, - index_v * cdf_width + res_x); - float2 cdf_last_v = kernel_tex_fetch(__light_background_marginal_cdf, res_y); - - float denom = (M_2PI_F * M_PI_F * sin_theta) * cdf_last_u.x * cdf_last_v.x; - - if (denom == 0.0f) - return 0.0f; - - /* pdfs in U direction */ - float2 cdf_u = kernel_tex_fetch(__light_background_conditional_cdf, - index_v * cdf_width + index_u); - float2 cdf_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v); - - return (cdf_u.x * cdf_v.x) / denom; -} - -ccl_device_inline bool background_portal_data_fetch_and_check_side( - KernelGlobals *kg, float3 P, int index, float3 *lightpos, float3 *dir) -{ - int portal = kernel_data.integrator.portal_offset + index; - const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); - - *lightpos = make_float3(klight->co[0], klight->co[1], klight->co[2]); - *dir = make_float3(klight->area.dir[0], klight->area.dir[1], klight->area.dir[2]); - - /* Check whether portal is on the right side. */ - if (dot(*dir, P - *lightpos) > 1e-4f) - return true; - - return false; -} - -ccl_device_inline float background_portal_pdf( - KernelGlobals *kg, float3 P, float3 direction, int ignore_portal, bool *is_possible) -{ - float portal_pdf = 0.0f; - - int num_possible = 0; - for (int p = 0; p < kernel_data.integrator.num_portals; p++) { - if (p == ignore_portal) - continue; - - float3 lightpos, dir; - if (!background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) - continue; - - /* There's a portal that could be sampled from this position. */ - if (is_possible) { - *is_possible = true; - } - num_possible++; - - int portal = kernel_data.integrator.portal_offset + p; - const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); - float3 axisu = make_float3( - klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]); - float3 axisv = make_float3( - klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]); - bool is_round = (klight->area.invarea < 0.0f); - - if (!ray_quad_intersect(P, - direction, - 1e-4f, - FLT_MAX, - lightpos, - axisu, - axisv, - dir, - NULL, - NULL, - NULL, - NULL, - is_round)) - continue; - - if (is_round) { - float t; - float3 D = normalize_len(lightpos - P, &t); - portal_pdf += fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t); - } - else { - portal_pdf += rect_light_sample(P, &lightpos, axisu, axisv, 0.0f, 0.0f, false); - } - } - - if (ignore_portal >= 0) { - /* We have skipped a portal that could be sampled as well. */ - num_possible++; - } - - return (num_possible > 0) ? portal_pdf / num_possible : 0.0f; -} - -ccl_device int background_num_possible_portals(KernelGlobals *kg, float3 P) -{ - int num_possible_portals = 0; - for (int p = 0; p < kernel_data.integrator.num_portals; p++) { - float3 lightpos, dir; - if (background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) - num_possible_portals++; - } - return num_possible_portals; -} - -ccl_device float3 background_portal_sample(KernelGlobals *kg, - float3 P, - float randu, - float randv, - int num_possible, - int *sampled_portal, - float *pdf) -{ - /* Pick a portal, then re-normalize randv. */ - randv *= num_possible; - int portal = (int)randv; - randv -= portal; - - /* TODO(sergey): Some smarter way of finding portal to sample - * is welcome. - */ - for (int p = 0; p < kernel_data.integrator.num_portals; p++) { - /* Search for the sampled portal. */ - float3 lightpos, dir; - if (!background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) - continue; - - if (portal == 0) { - /* p is the portal to be sampled. */ - int portal = kernel_data.integrator.portal_offset + p; - const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); - float3 axisu = make_float3( - klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]); - float3 axisv = make_float3( - klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]); - bool is_round = (klight->area.invarea < 0.0f); - - float3 D; - if (is_round) { - lightpos += ellipse_sample(axisu * 0.5f, axisv * 0.5f, randu, randv); - float t; - D = normalize_len(lightpos - P, &t); - *pdf = fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t); - } - else { - *pdf = rect_light_sample(P, &lightpos, axisu, axisv, randu, randv, true); - D = normalize(lightpos - P); - } - - *pdf /= num_possible; - *sampled_portal = p; - return D; - } - - portal--; - } - - return make_float3(0.0f, 0.0f, 0.0f); -} - -ccl_device_inline float3 -background_light_sample(KernelGlobals *kg, float3 P, float randu, float randv, float *pdf) -{ - /* Probability of sampling portals instead of the map. */ - float portal_sampling_pdf = kernel_data.integrator.portal_pdf; - - /* Check if there are portals in the scene which we can sample. */ - if (portal_sampling_pdf > 0.0f) { - int num_portals = background_num_possible_portals(kg, P); - if (num_portals > 0) { - if (portal_sampling_pdf == 1.0f || randu < portal_sampling_pdf) { - if (portal_sampling_pdf < 1.0f) { - randu /= portal_sampling_pdf; - } - int portal; - float3 D = background_portal_sample(kg, P, randu, randv, num_portals, &portal, pdf); - if (num_portals > 1) { - /* Ignore the chosen portal, its pdf is already included. */ - *pdf += background_portal_pdf(kg, P, D, portal, NULL); - } - /* We could also have sampled the map, so combine with MIS. */ - if (portal_sampling_pdf < 1.0f) { - float cdf_pdf = background_map_pdf(kg, D); - *pdf = (portal_sampling_pdf * (*pdf) + (1.0f - portal_sampling_pdf) * cdf_pdf); - } - return D; - } - else { - /* Sample map, but with nonzero portal_sampling_pdf for MIS. */ - randu = (randu - portal_sampling_pdf) / (1.0f - portal_sampling_pdf); - } - } - else { - /* We can't sample a portal. - * Check if we can sample the map instead. - */ - if (portal_sampling_pdf == 1.0f) { - /* Use uniform as a fallback if we can't sample the map. */ - *pdf = 1.0f / M_4PI_F; - return sample_uniform_sphere(randu, randv); - } - else { - portal_sampling_pdf = 0.0f; - } - } - } - - float3 D = background_map_sample(kg, randu, randv, pdf); - /* Use MIS if portals could be sampled as well. */ - if (portal_sampling_pdf > 0.0f) { - float portal_pdf = background_portal_pdf(kg, P, D, -1, NULL); - *pdf = (portal_sampling_pdf * portal_pdf + (1.0f - portal_sampling_pdf) * (*pdf)); - } - return D; -} - -ccl_device float background_light_pdf(KernelGlobals *kg, float3 P, float3 direction) -{ - /* Probability of sampling portals instead of the map. */ - float portal_sampling_pdf = kernel_data.integrator.portal_pdf; - - float portal_pdf = 0.0f, map_pdf = 0.0f; - if (portal_sampling_pdf > 0.0f) { - /* Evaluate PDF of sampling this direction by portal sampling. */ - bool is_possible = false; - portal_pdf = background_portal_pdf(kg, P, direction, -1, &is_possible) * portal_sampling_pdf; - if (!is_possible) { - /* Portal sampling is not possible here because all portals point to the wrong side. - * If map sampling is possible, it would be used instead, - * otherwise fallback sampling is used. */ - if (portal_sampling_pdf == 1.0f) { - return kernel_data.integrator.pdf_lights / M_4PI_F; - } - else { - /* Force map sampling. */ - portal_sampling_pdf = 0.0f; - } - } - } - if (portal_sampling_pdf < 1.0f) { - /* Evaluate PDF of sampling this direction by map sampling. */ - map_pdf = background_map_pdf(kg, direction) * (1.0f - portal_sampling_pdf); - } - return (portal_pdf + map_pdf) * kernel_data.integrator.pdf_lights; -} -#endif - /* Regular Light */ ccl_device_inline bool lamp_light_sample( @@ -594,7 +102,7 @@ ccl_device_inline bool lamp_light_sample( /* spot light attenuation */ float3 dir = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]); ls->eval_fac *= spot_light_attenuation( - dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls); + dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls->Ng); if (ls->eval_fac == 0.0f) { return false; } @@ -732,7 +240,7 @@ ccl_device bool lamp_light_eval( /* spot light attenuation */ float3 dir = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]); ls->eval_fac *= spot_light_attenuation( - dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls); + dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls->Ng); if (ls->eval_fac == 0.0f) return false; diff --git a/intern/cycles/kernel/kernel_light_background.h b/intern/cycles/kernel/kernel_light_background.h new file mode 100644 index 00000000000..30e336f0f80 --- /dev/null +++ b/intern/cycles/kernel/kernel_light_background.h @@ -0,0 +1,448 @@ +/* + * Copyright 2011-2020 Blender Foundation + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "kernel_light_common.h" + +CCL_NAMESPACE_BEGIN + +/* Background Light */ + +#ifdef __BACKGROUND_MIS__ + +ccl_device float3 background_map_sample(KernelGlobals *kg, float randu, float randv, float *pdf) +{ + /* for the following, the CDF values are actually a pair of floats, with the + * function value as X and the actual CDF as Y. The last entry's function + * value is the CDF total. */ + int res_x = kernel_data.background.map_res_x; + int res_y = kernel_data.background.map_res_y; + int cdf_width = res_x + 1; + + /* this is basically std::lower_bound as used by pbrt */ + int first = 0; + int count = res_y; + + while (count > 0) { + int step = count >> 1; + int middle = first + step; + + if (kernel_tex_fetch(__light_background_marginal_cdf, middle).y < randv) { + first = middle + 1; + count -= step + 1; + } + else + count = step; + } + + int index_v = max(0, first - 1); + kernel_assert(index_v >= 0 && index_v < res_y); + + float2 cdf_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v); + float2 cdf_next_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v + 1); + float2 cdf_last_v = kernel_tex_fetch(__light_background_marginal_cdf, res_y); + + /* importance-sampled V direction */ + float dv = inverse_lerp(cdf_v.y, cdf_next_v.y, randv); + float v = (index_v + dv) / res_y; + + /* this is basically std::lower_bound as used by pbrt */ + first = 0; + count = res_x; + while (count > 0) { + int step = count >> 1; + int middle = first + step; + + if (kernel_tex_fetch(__light_background_conditional_cdf, index_v * cdf_width + middle).y < + randu) { + first = middle + 1; + count -= step + 1; + } + else + count = step; + } + + int index_u = max(0, first - 1); + kernel_assert(index_u >= 0 && index_u < res_x); + + float2 cdf_u = kernel_tex_fetch(__light_background_conditional_cdf, + index_v * cdf_width + index_u); + float2 cdf_next_u = kernel_tex_fetch(__light_background_conditional_cdf, + index_v * cdf_width + index_u + 1); + float2 cdf_last_u = kernel_tex_fetch(__light_background_conditional_cdf, + index_v * cdf_width + res_x); + + /* importance-sampled U direction */ + float du = inverse_lerp(cdf_u.y, cdf_next_u.y, randu); + float u = (index_u + du) / res_x; + + /* compute pdf */ + float sin_theta = sinf(M_PI_F * v); + float denom = (M_2PI_F * M_PI_F * sin_theta) * cdf_last_u.x * cdf_last_v.x; + + if (sin_theta == 0.0f || denom == 0.0f) + *pdf = 0.0f; + else + *pdf = (cdf_u.x * cdf_v.x) / denom; + + /* compute direction */ + return equirectangular_to_direction(u, v); +} + +/* TODO(sergey): Same as above, after the release we should consider using + * 'noinline' for all devices. + */ +ccl_device float background_map_pdf(KernelGlobals *kg, float3 direction) +{ + float2 uv = direction_to_equirectangular(direction); + int res_x = kernel_data.background.map_res_x; + int res_y = kernel_data.background.map_res_y; + int cdf_width = res_x + 1; + + float sin_theta = sinf(uv.y * M_PI_F); + + if (sin_theta == 0.0f) + return 0.0f; + + int index_u = clamp(float_to_int(uv.x * res_x), 0, res_x - 1); + int index_v = clamp(float_to_int(uv.y * res_y), 0, res_y - 1); + + /* pdfs in V direction */ + float2 cdf_last_u = kernel_tex_fetch(__light_background_conditional_cdf, + index_v * cdf_width + res_x); + float2 cdf_last_v = kernel_tex_fetch(__light_background_marginal_cdf, res_y); + + float denom = (M_2PI_F * M_PI_F * sin_theta) * cdf_last_u.x * cdf_last_v.x; + + if (denom == 0.0f) + return 0.0f; + + /* pdfs in U direction */ + float2 cdf_u = kernel_tex_fetch(__light_background_conditional_cdf, + index_v * cdf_width + index_u); + float2 cdf_v = kernel_tex_fetch(__light_background_marginal_cdf, index_v); + + return (cdf_u.x * cdf_v.x) / denom; +} + +ccl_device_inline bool background_portal_data_fetch_and_check_side( + KernelGlobals *kg, float3 P, int index, float3 *lightpos, float3 *dir) +{ + int portal = kernel_data.background.portal_offset + index; + const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); + + *lightpos = make_float3(klight->co[0], klight->co[1], klight->co[2]); + *dir = make_float3(klight->area.dir[0], klight->area.dir[1], klight->area.dir[2]); + + /* Check whether portal is on the right side. */ + if (dot(*dir, P - *lightpos) > 1e-4f) + return true; + + return false; +} + +ccl_device_inline float background_portal_pdf( + KernelGlobals *kg, float3 P, float3 direction, int ignore_portal, bool *is_possible) +{ + float portal_pdf = 0.0f; + + int num_possible = 0; + for (int p = 0; p < kernel_data.background.num_portals; p++) { + if (p == ignore_portal) + continue; + + float3 lightpos, dir; + if (!background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) + continue; + + /* There's a portal that could be sampled from this position. */ + if (is_possible) { + *is_possible = true; + } + num_possible++; + + int portal = kernel_data.background.portal_offset + p; + const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); + float3 axisu = make_float3( + klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]); + float3 axisv = make_float3( + klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]); + bool is_round = (klight->area.invarea < 0.0f); + + if (!ray_quad_intersect(P, + direction, + 1e-4f, + FLT_MAX, + lightpos, + axisu, + axisv, + dir, + NULL, + NULL, + NULL, + NULL, + is_round)) + continue; + + if (is_round) { + float t; + float3 D = normalize_len(lightpos - P, &t); + portal_pdf += fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t); + } + else { + portal_pdf += rect_light_sample(P, &lightpos, axisu, axisv, 0.0f, 0.0f, false); + } + } + + if (ignore_portal >= 0) { + /* We have skipped a portal that could be sampled as well. */ + num_possible++; + } + + return (num_possible > 0) ? portal_pdf / num_possible : 0.0f; +} + +ccl_device int background_num_possible_portals(KernelGlobals *kg, float3 P) +{ + int num_possible_portals = 0; + for (int p = 0; p < kernel_data.background.num_portals; p++) { + float3 lightpos, dir; + if (background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) + num_possible_portals++; + } + return num_possible_portals; +} + +ccl_device float3 background_portal_sample(KernelGlobals *kg, + float3 P, + float randu, + float randv, + int num_possible, + int *sampled_portal, + float *pdf) +{ + /* Pick a portal, then re-normalize randv. */ + randv *= num_possible; + int portal = (int)randv; + randv -= portal; + + /* TODO(sergey): Some smarter way of finding portal to sample + * is welcome. + */ + for (int p = 0; p < kernel_data.background.num_portals; p++) { + /* Search for the sampled portal. */ + float3 lightpos, dir; + if (!background_portal_data_fetch_and_check_side(kg, P, p, &lightpos, &dir)) + continue; + + if (portal == 0) { + /* p is the portal to be sampled. */ + int portal = kernel_data.background.portal_offset + p; + const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, portal); + float3 axisu = make_float3( + klight->area.axisu[0], klight->area.axisu[1], klight->area.axisu[2]); + float3 axisv = make_float3( + klight->area.axisv[0], klight->area.axisv[1], klight->area.axisv[2]); + bool is_round = (klight->area.invarea < 0.0f); + + float3 D; + if (is_round) { + lightpos += ellipse_sample(axisu * 0.5f, axisv * 0.5f, randu, randv); + float t; + D = normalize_len(lightpos - P, &t); + *pdf = fabsf(klight->area.invarea) * lamp_light_pdf(kg, dir, -D, t); + } + else { + *pdf = rect_light_sample(P, &lightpos, axisu, axisv, randu, randv, true); + D = normalize(lightpos - P); + } + + *pdf /= num_possible; + *sampled_portal = p; + return D; + } + + portal--; + } + + return make_float3(0.0f, 0.0f, 0.0f); +} + +ccl_device_inline float3 background_sun_sample(KernelGlobals *kg, + float randu, + float randv, + float *pdf) +{ + float3 D; + const float3 N = float4_to_float3(kernel_data.background.sun); + const float angle = kernel_data.background.sun.w; + sample_uniform_cone(N, angle, randu, randv, &D, pdf); + return D; +} + +ccl_device_inline float background_sun_pdf(KernelGlobals *kg, float3 D) +{ + const float3 N = float4_to_float3(kernel_data.background.sun); + const float angle = kernel_data.background.sun.w; + return pdf_uniform_cone(N, D, angle); +} + +ccl_device_inline float3 +background_light_sample(KernelGlobals *kg, float3 P, float randu, float randv, float *pdf) +{ + float portal_method_pdf = kernel_data.background.portal_weight; + float sun_method_pdf = kernel_data.background.sun_weight; + float map_method_pdf = kernel_data.background.map_weight; + + int num_portals = 0; + if (portal_method_pdf > 0.0f) { + /* Check if there are portals in the scene which we can sample. */ + num_portals = background_num_possible_portals(kg, P); + if (num_portals == 0) { + portal_method_pdf = 0.0f; + } + } + + float pdf_fac = (portal_method_pdf + sun_method_pdf + map_method_pdf); + if (pdf_fac == 0.0f) { + /* Use uniform as a fallback if we can't use any strategy. */ + *pdf = 1.0f / M_4PI_F; + return sample_uniform_sphere(randu, randv); + } + + pdf_fac = 1.0f / pdf_fac; + portal_method_pdf *= pdf_fac; + sun_method_pdf *= pdf_fac; + map_method_pdf *= pdf_fac; + + /* We have 100% in total and split it between the three categories. + * Therefore, we pick portals if randu is between 0 and portal_method_pdf, + * sun if randu is between portal_method_pdf and (portal_method_pdf + sun_method_pdf) + * and map if randu is between (portal_method_pdf + sun_method_pdf) and 1. */ + float sun_method_cdf = portal_method_pdf + sun_method_pdf; + + int method = 0; + float3 D; + if (randu < portal_method_pdf) { + method = 0; + /* Rescale randu. */ + if (portal_method_pdf != 1.0f) { + randu /= portal_method_pdf; + } + + /* Sample a portal. */ + int portal; + D = background_portal_sample(kg, P, randu, randv, num_portals, &portal, pdf); + if (num_portals > 1) { + /* Ignore the chosen portal, its pdf is already included. */ + *pdf += background_portal_pdf(kg, P, D, portal, NULL); + } + + /* Skip MIS if this is the only method. */ + if (portal_method_pdf == 1.0f) { + return D; + } + *pdf *= portal_method_pdf; + } + else if (randu < sun_method_cdf) { + method = 1; + /* Rescale randu. */ + if (sun_method_pdf != 1.0f) { + randu = (randu - portal_method_pdf) / sun_method_pdf; + } + + D = background_sun_sample(kg, randu, randv, pdf); + + /* Skip MIS if this is the only method. */ + if (sun_method_pdf == 1.0f) { + return D; + } + *pdf *= sun_method_pdf; + } + else { + method = 2; + /* Rescale randu. */ + if (map_method_pdf != 1.0f) { + randu = (randu - sun_method_cdf) / map_method_pdf; + } + + D = background_map_sample(kg, randu, randv, pdf); + + /* Skip MIS if this is the only method. */ + if (map_method_pdf == 1.0f) { + return D; + } + *pdf *= map_method_pdf; + } + + /* MIS weighting. */ + if (method != 0 && portal_method_pdf != 0.0f) { + *pdf += portal_method_pdf * background_portal_pdf(kg, P, D, -1, NULL); + } + if (method != 1 && sun_method_pdf != 0.0f) { + *pdf += sun_method_pdf * background_sun_pdf(kg, D); + } + if (method != 2 && map_method_pdf != 0.0f) { + *pdf += map_method_pdf * background_map_pdf(kg, D); + } + return D; +} + +ccl_device float background_light_pdf(KernelGlobals *kg, float3 P, float3 direction) +{ + float portal_method_pdf = kernel_data.background.portal_weight; + float sun_method_pdf = kernel_data.background.sun_weight; + float map_method_pdf = kernel_data.background.map_weight; + + float portal_pdf = 0.0f; + /* Portals are a special case here since we need to compute their pdf in order + * to find out if we can sample them. */ + if (portal_method_pdf > 0.0f) { + /* Evaluate PDF of sampling this direction by portal sampling. */ + bool is_possible = false; + portal_pdf = background_portal_pdf(kg, P, direction, -1, &is_possible); + if (!is_possible) { + /* Portal sampling is not possible here because all portals point to the wrong side. + * If other methods can be used instead, do so, otherwise uniform sampling is used as a + * fallback. */ + portal_method_pdf = 0.0f; + } + } + + float pdf_fac = (portal_method_pdf + sun_method_pdf + map_method_pdf); + if (pdf_fac == 0.0f) { + /* Use uniform as a fallback if we can't use any strategy. */ + return kernel_data.integrator.pdf_lights / M_4PI_F; + } + + pdf_fac = 1.0f / pdf_fac; + portal_method_pdf *= pdf_fac; + sun_method_pdf *= pdf_fac; + map_method_pdf *= pdf_fac; + + float pdf = portal_pdf * portal_method_pdf; + if (sun_method_pdf != 0.0f) { + pdf += background_sun_pdf(kg, direction) * sun_method_pdf; + } + if (map_method_pdf != 0.0f) { + pdf += background_map_pdf(kg, direction) * map_method_pdf; + } + + return pdf * kernel_data.integrator.pdf_lights; +} + +#endif + +CCL_NAMESPACE_END
\ No newline at end of file diff --git a/intern/cycles/kernel/kernel_light_common.h b/intern/cycles/kernel/kernel_light_common.h new file mode 100644 index 00000000000..39503a4b479 --- /dev/null +++ b/intern/cycles/kernel/kernel_light_common.h @@ -0,0 +1,159 @@ +/* + * Copyright 2011-2020 Blender Foundation + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +CCL_NAMESPACE_BEGIN + +/* Area light sampling */ + +/* Uses the following paper: + * + * Carlos Urena et al. + * An Area-Preserving Parametrization for Spherical Rectangles. + * + * https://www.solidangle.com/research/egsr2013_spherical_rectangle.pdf + * + * Note: light_p is modified when sample_coord is true. + */ +ccl_device_inline float rect_light_sample(float3 P, + float3 *light_p, + float3 axisu, + float3 axisv, + float randu, + float randv, + bool sample_coord) +{ + /* In our name system we're using P for the center, + * which is o in the paper. + */ + + float3 corner = *light_p - axisu * 0.5f - axisv * 0.5f; + float axisu_len, axisv_len; + /* Compute local reference system R. */ + float3 x = normalize_len(axisu, &axisu_len); + float3 y = normalize_len(axisv, &axisv_len); + float3 z = cross(x, y); + /* Compute rectangle coords in local reference system. */ + float3 dir = corner - P; + float z0 = dot(dir, z); + /* Flip 'z' to make it point against Q. */ + if (z0 > 0.0f) { + z *= -1.0f; + z0 *= -1.0f; + } + float x0 = dot(dir, x); + float y0 = dot(dir, y); + float x1 = x0 + axisu_len; + float y1 = y0 + axisv_len; + /* Compute internal angles (gamma_i). */ + float4 diff = make_float4(x0, y1, x1, y0) - make_float4(x1, y0, x0, y1); + float4 nz = make_float4(y0, x1, y1, x0) * diff; + nz = nz / sqrt(z0 * z0 * diff * diff + nz * nz); + float g0 = safe_acosf(-nz.x * nz.y); + float g1 = safe_acosf(-nz.y * nz.z); + float g2 = safe_acosf(-nz.z * nz.w); + float g3 = safe_acosf(-nz.w * nz.x); + /* Compute predefined constants. */ + float b0 = nz.x; + float b1 = nz.z; + float b0sq = b0 * b0; + float k = M_2PI_F - g2 - g3; + /* Compute solid angle from internal angles. */ + float S = g0 + g1 - k; + + if (sample_coord) { + /* Compute cu. */ + float au = randu * S + k; + float fu = (cosf(au) * b0 - b1) / sinf(au); + float cu = 1.0f / sqrtf(fu * fu + b0sq) * (fu > 0.0f ? 1.0f : -1.0f); + cu = clamp(cu, -1.0f, 1.0f); + /* Compute xu. */ + float xu = -(cu * z0) / max(sqrtf(1.0f - cu * cu), 1e-7f); + xu = clamp(xu, x0, x1); + /* Compute yv. */ + float z0sq = z0 * z0; + float y0sq = y0 * y0; + float y1sq = y1 * y1; + float d = sqrtf(xu * xu + z0sq); + float h0 = y0 / sqrtf(d * d + y0sq); + float h1 = y1 / sqrtf(d * d + y1sq); + float hv = h0 + randv * (h1 - h0), hv2 = hv * hv; + float yv = (hv2 < 1.0f - 1e-6f) ? (hv * d) / sqrtf(1.0f - hv2) : y1; + + /* Transform (xu, yv, z0) to world coords. */ + *light_p = P + xu * x + yv * y + z0 * z; + } + + /* return pdf */ + if (S != 0.0f) + return 1.0f / S; + else + return 0.0f; +} + +ccl_device_inline float3 ellipse_sample(float3 ru, float3 rv, float randu, float randv) +{ + to_unit_disk(&randu, &randv); + return ru * randu + rv * randv; +} + +ccl_device float3 disk_light_sample(float3 v, float randu, float randv) +{ + float3 ru, rv; + + make_orthonormals(v, &ru, &rv); + + return ellipse_sample(ru, rv, randu, randv); +} + +ccl_device float3 distant_light_sample(float3 D, float radius, float randu, float randv) +{ + return normalize(D + disk_light_sample(D, randu, randv) * radius); +} + +ccl_device float3 +sphere_light_sample(float3 P, float3 center, float radius, float randu, float randv) +{ + return disk_light_sample(normalize(P - center), randu, randv) * radius; +} + +ccl_device float spot_light_attenuation(float3 dir, float spot_angle, float spot_smooth, float3 N) +{ + float attenuation = dot(dir, N); + + if (attenuation <= spot_angle) { + attenuation = 0.0f; + } + else { + float t = attenuation - spot_angle; + + if (t < spot_smooth && spot_smooth != 0.0f) + attenuation *= smoothstepf(t / spot_smooth); + } + + return attenuation; +} + +ccl_device float lamp_light_pdf(KernelGlobals *kg, const float3 Ng, const float3 I, float t) +{ + float cos_pi = dot(Ng, I); + + if (cos_pi <= 0.0f) + return 0.0f; + + return t * t / cos_pi; +} + +CCL_NAMESPACE_END diff --git a/intern/cycles/kernel/kernel_montecarlo.h b/intern/cycles/kernel/kernel_montecarlo.h index 5c776e06547..0edcc1a5a14 100644 --- a/intern/cycles/kernel/kernel_montecarlo.h +++ b/intern/cycles/kernel/kernel_montecarlo.h @@ -98,6 +98,16 @@ ccl_device_inline void sample_uniform_cone( *pdf = M_1_2PI_F / (1.0f - zMin); } +ccl_device_inline float pdf_uniform_cone(const float3 N, float3 D, float angle) +{ + float zMin = cosf(angle); + float z = dot(N, D); + if (z > zMin) { + return M_1_2PI_F / (1.0f - zMin); + } + return 0.0f; +} + /* sample uniform point on the surface of a sphere */ ccl_device float3 sample_uniform_sphere(float u1, float u2) { diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h index 0a0cf1bd6c0..a692d7a844f 100644 --- a/intern/cycles/kernel/kernel_types.h +++ b/intern/cycles/kernel/kernel_types.h @@ -1291,6 +1291,24 @@ typedef struct KernelBackground { float ao_factor; float ao_distance; float ao_bounces_factor; + + /* portal sampling */ + float portal_weight; + int num_portals; + int portal_offset; + + /* sun sampling */ + float sun_weight; + /* xyz store direction, w the angle. float4 instead of float3 is used + * to ensure consistent padding/alignment across devices. */ + float4 sun; + + /* map sampling */ + float map_weight; + int map_res_x; + int map_res_y; + + int use_mis; } KernelBackground; static_assert_align(KernelBackground, 16); @@ -1302,15 +1320,8 @@ typedef struct KernelIntegrator { int num_all_lights; float pdf_triangles; float pdf_lights; - int pdf_background_res_x; - int pdf_background_res_y; float light_inv_rr_threshold; - /* light portals */ - float portal_pdf; - int num_portals; - int portal_offset; - /* bounces */ int min_bounce; int max_bounce; @@ -1372,7 +1383,7 @@ typedef struct KernelIntegrator { int max_closures; - int pad1; + int pad1, pad2; } KernelIntegrator; static_assert_align(KernelIntegrator, 16); diff --git a/intern/cycles/kernel/shaders/node_sky_texture.osl b/intern/cycles/kernel/shaders/node_sky_texture.osl index 4def237a2e0..08bc8f85120 100644 --- a/intern/cycles/kernel/shaders/node_sky_texture.osl +++ b/intern/cycles/kernel/shaders/node_sky_texture.osl @@ -44,13 +44,13 @@ float sky_perez_function(float lam[9], float theta, float gamma) (1.0 + lam[2] * exp(lam[3] * gamma) + lam[4] * cgamma * cgamma); } -color sky_radiance_old(normal dir, - float sunphi, - float suntheta, - color radiance, - float config_x[9], - float config_y[9], - float config_z[9]) +color sky_radiance_preetham(normal dir, + float sunphi, + float suntheta, + color radiance, + float config_x[9], + float config_y[9], + float config_z[9]) { /* convert vector to spherical coordinates */ vector spherical = sky_spherical_coordinates(dir); @@ -88,13 +88,13 @@ float sky_radiance_internal(float config[9], float theta, float gamma) (config[2] + config[3] * expM + config[5] * rayM + config[6] * mieM + config[7] * zenith); } -color sky_radiance_new(normal dir, - float sunphi, - float suntheta, - color radiance, - float config_x[9], - float config_y[9], - float config_z[9]) +color sky_radiance_hosek(normal dir, + float sunphi, + float suntheta, + color radiance, + float config_x[9], + float config_y[9], + float config_z[9]) { /* convert vector to spherical coordinates */ vector spherical = sky_spherical_coordinates(dir); @@ -116,16 +116,103 @@ color sky_radiance_new(normal dir, return xyz_to_rgb(x, y, z) * (M_2PI / 683); } +/* Nishita improved */ +vector geographical_to_direction(float lat, float lon) +{ + return vector(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat)); +} + +color sky_radiance_nishita(vector dir, float nishita_data[9], string filename) +{ + /* definitions */ + float sun_elevation = nishita_data[6]; + float sun_rotation = nishita_data[7]; + float angular_diameter = nishita_data[8]; + int sun_disc = angular_diameter > 0; + float alpha = 1.0; + color xyz; + /* convert dir to spherical coordinates */ + vector direction = sky_spherical_coordinates(dir); + + /* render above the horizon */ + if (dir[2] >= 0.0) { + /* definitions */ + vector sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2); + float sun_dir_angle = acos(dot(dir, sun_dir)); + float half_angular = angular_diameter / 2.0; + float dir_elevation = M_PI_2 - direction[0]; + + /* if ray inside sun disc render it, otherwise render sky */ + if (sun_dir_angle < half_angular && sun_disc == 1) { + /* get 3 pixels data */ + color pixel_bottom = color(nishita_data[0], nishita_data[1], nishita_data[2]); + color pixel_top = color(nishita_data[3], nishita_data[4], nishita_data[5]); + float y; + + /* sun interpolation */ + if (sun_elevation - half_angular > 0.0) { + if ((sun_elevation + half_angular) > 0.0) { + y = ((dir_elevation - sun_elevation) / angular_diameter) + 0.5; + xyz = mix(pixel_bottom, pixel_top, y); + } + } + else { + if (sun_elevation + half_angular > 0.0) { + y = dir_elevation / (sun_elevation + half_angular); + xyz = mix(pixel_bottom, pixel_top, y); + } + } + /* limb darkening, coefficient is 0.6f */ + float angle_fraction = sun_dir_angle / half_angular; + float limb_darkening = (1.0 - 0.6 * (1.0 - sqrt(1.0 - angle_fraction * angle_fraction))); + xyz *= limb_darkening; + } + /* sky */ + else { + /* sky interpolation */ + float x = (direction[1] + M_PI + sun_rotation) / M_2PI; + float y = 1.0 - (dir_elevation / M_PI_2); + if (x > 1.0) { + x = x - 1.0; + } + xyz = (color)texture(filename, x, y, "wrap", "clamp", "interp", "linear", "alpha", alpha); + } + } + /* ground */ + else { + if (dir[2] < -0.4) { + xyz = color(0, 0, 0); + } + else { + /* black ground fade */ + float mul = pow(1.0 + dir[2] * 2.5, 3.0); + /* interpolation */ + float x = (direction[1] + M_PI + sun_rotation) / M_2PI; + float y = 1.5; + if (x > 1.0) { + x = x - 1.0; + } + xyz = (color)texture( + filename, x, y, "wrap", "periodic", "interp", "linear", "alpha", alpha) * + mul; + } + } + /* convert to RGB and adjust strength */ + return xyz_to_rgb(xyz[0], xyz[1], xyz[2]) * 120000.0; +} + shader node_sky_texture(int use_mapping = 0, matrix mapping = matrix(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), vector Vector = P, string type = "hosek_wilkie", float theta = 0.0, float phi = 0.0, + string filename = "", color radiance = color(0.0, 0.0, 0.0), float config_x[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, float config_y[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, float config_z[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + float nishita_data[9] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, output color Color = color(0.0, 0.0, 0.0)) { vector p = Vector; @@ -133,8 +220,10 @@ shader node_sky_texture(int use_mapping = 0, if (use_mapping) p = transform(mapping, p); + if (type == "nishita_improved") + Color = sky_radiance_nishita(p, nishita_data, filename); if (type == "hosek_wilkie") - Color = sky_radiance_new(p, phi, theta, radiance, config_x, config_y, config_z); - else - Color = sky_radiance_old(p, phi, theta, radiance, config_x, config_y, config_z); + Color = sky_radiance_hosek(p, phi, theta, radiance, config_x, config_y, config_z); + if (type == "preetham") + Color = sky_radiance_preetham(p, phi, theta, radiance, config_x, config_y, config_z); } diff --git a/intern/cycles/kernel/svm/svm_sky.h b/intern/cycles/kernel/svm/svm_sky.h index 50fe0c8232f..e877bd9a5c8 100644 --- a/intern/cycles/kernel/svm/svm_sky.h +++ b/intern/cycles/kernel/svm/svm_sky.h @@ -37,16 +37,16 @@ ccl_device float sky_perez_function(float *lam, float theta, float gamma) (1.0f + lam[2] * expf(lam[3] * gamma) + lam[4] * cgamma * cgamma); } -ccl_device float3 sky_radiance_old(KernelGlobals *kg, - float3 dir, - float sunphi, - float suntheta, - float radiance_x, - float radiance_y, - float radiance_z, - float *config_x, - float *config_y, - float *config_z) +ccl_device float3 sky_radiance_preetham(KernelGlobals *kg, + float3 dir, + float sunphi, + float suntheta, + float radiance_x, + float radiance_y, + float radiance_z, + float *config_x, + float *config_y, + float *config_z) { /* convert vector to spherical coordinates */ float2 spherical = direction_to_spherical(dir); @@ -90,16 +90,16 @@ ccl_device float sky_radiance_internal(float *configuration, float theta, float configuration[6] * mieM + configuration[7] * zenith); } -ccl_device float3 sky_radiance_new(KernelGlobals *kg, - float3 dir, - float sunphi, - float suntheta, - float radiance_x, - float radiance_y, - float radiance_z, - float *config_x, - float *config_y, - float *config_z) +ccl_device float3 sky_radiance_hosek(KernelGlobals *kg, + float3 dir, + float sunphi, + float suntheta, + float radiance_x, + float radiance_y, + float radiance_z, + float *config_x, + float *config_y, + float *config_z) { /* convert vector to spherical coordinates */ float2 spherical = direction_to_spherical(dir); @@ -121,93 +121,206 @@ ccl_device float3 sky_radiance_new(KernelGlobals *kg, return xyz_to_rgb(kg, make_float3(x, y, z)) * (M_2PI_F / 683); } +/* Nishita improved sky model */ +ccl_device float3 geographical_to_direction(float lat, float lon) +{ + return make_float3(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat)); +} + +ccl_device float3 sky_radiance_nishita(KernelGlobals *kg, + float3 dir, + float *nishita_data, + uint texture_id) +{ + /* definitions */ + float sun_elevation = nishita_data[6]; + float sun_rotation = nishita_data[7]; + float angular_diameter = nishita_data[8]; + bool sun_disc = (angular_diameter > 0.0f); + float3 xyz; + /* convert dir to spherical coordinates */ + float2 direction = direction_to_spherical(dir); + + /* render above the horizon */ + if (dir.z >= 0.0f) { + /* definitions */ + float3 sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2_F); + float sun_dir_angle = acos(dot(dir, sun_dir)); + float half_angular = angular_diameter / 2.0f; + float dir_elevation = M_PI_2_F - direction.x; + + /* if ray inside sun disc render it, otherwise render sky */ + if (sun_disc && sun_dir_angle < half_angular) { + /* get 3 pixels data */ + float3 pixel_bottom = make_float3(nishita_data[0], nishita_data[1], nishita_data[2]); + float3 pixel_top = make_float3(nishita_data[3], nishita_data[4], nishita_data[5]); + float y; + + /* sun interpolation */ + if (sun_elevation - half_angular > 0.0f) { + if (sun_elevation + half_angular > 0.0f) { + y = ((dir_elevation - sun_elevation) / angular_diameter) + 0.5f; + xyz = interp(pixel_bottom, pixel_top, y); + } + } + else { + if (sun_elevation + half_angular > 0.0f) { + y = dir_elevation / (sun_elevation + half_angular); + xyz = interp(pixel_bottom, pixel_top, y); + } + } + /* limb darkening, coefficient is 0.6f */ + float limb_darkening = (1.0f - + 0.6f * (1.0f - sqrtf(1.0f - sqr(sun_dir_angle / half_angular)))); + xyz *= limb_darkening; + } + /* sky */ + else { + /* sky interpolation */ + float x = (direction.y + M_PI_F + sun_rotation) / M_2PI_F; + float y = dir_elevation / M_PI_2_F; + if (x > 1.0f) { + x -= 1.0f; + } + xyz = float4_to_float3(kernel_tex_image_interp(kg, texture_id, x, y)); + } + } + /* ground */ + else { + if (dir.z < -0.4f) { + xyz = make_float3(0.0f, 0.0f, 0.0f); + } + else { + /* black ground fade */ + float fade = 1.0f + dir.z * 2.5f; + fade = sqr(fade) * fade; + /* interpolation */ + float x = (direction.y + M_PI_F + sun_rotation) / M_2PI_F; + if (x > 1.0f) { + x -= 1.0f; + } + xyz = float4_to_float3(kernel_tex_image_interp(kg, texture_id, x, -0.5)) * fade; + } + } + + /* convert to rgb and adjust strength */ + return xyz_to_rgb(kg, xyz) * 120000.0f; +} + ccl_device void svm_node_tex_sky( KernelGlobals *kg, ShaderData *sd, float *stack, uint4 node, int *offset) { - /* Define variables */ - float sunphi, suntheta, radiance_x, radiance_y, radiance_z; - float config_x[9], config_y[9], config_z[9]; - /* Load data */ uint dir_offset = node.y; uint out_offset = node.z; int sky_model = node.w; - float4 data = read_node_float(kg, offset); - sunphi = data.x; - suntheta = data.y; - radiance_x = data.z; - radiance_y = data.w; - - data = read_node_float(kg, offset); - radiance_z = data.x; - config_x[0] = data.y; - config_x[1] = data.z; - config_x[2] = data.w; - - data = read_node_float(kg, offset); - config_x[3] = data.x; - config_x[4] = data.y; - config_x[5] = data.z; - config_x[6] = data.w; - - data = read_node_float(kg, offset); - config_x[7] = data.x; - config_x[8] = data.y; - config_y[0] = data.z; - config_y[1] = data.w; - - data = read_node_float(kg, offset); - config_y[2] = data.x; - config_y[3] = data.y; - config_y[4] = data.z; - config_y[5] = data.w; - - data = read_node_float(kg, offset); - config_y[6] = data.x; - config_y[7] = data.y; - config_y[8] = data.z; - config_z[0] = data.w; - - data = read_node_float(kg, offset); - config_z[1] = data.x; - config_z[2] = data.y; - config_z[3] = data.z; - config_z[4] = data.w; - - data = read_node_float(kg, offset); - config_z[5] = data.x; - config_z[6] = data.y; - config_z[7] = data.z; - config_z[8] = data.w; - float3 dir = stack_load_float3(stack, dir_offset); float3 f; - /* Compute Sky */ - if (sky_model == 0) { - f = sky_radiance_old(kg, - dir, - sunphi, - suntheta, - radiance_x, - radiance_y, - radiance_z, - config_x, - config_y, - config_z); + /* Preetham and Hosek share the same data */ + if (sky_model == 0 || sky_model == 1) { + /* Define variables */ + float sunphi, suntheta, radiance_x, radiance_y, radiance_z; + float config_x[9], config_y[9], config_z[9]; + + float4 data = read_node_float(kg, offset); + sunphi = data.x; + suntheta = data.y; + radiance_x = data.z; + radiance_y = data.w; + + data = read_node_float(kg, offset); + radiance_z = data.x; + config_x[0] = data.y; + config_x[1] = data.z; + config_x[2] = data.w; + + data = read_node_float(kg, offset); + config_x[3] = data.x; + config_x[4] = data.y; + config_x[5] = data.z; + config_x[6] = data.w; + + data = read_node_float(kg, offset); + config_x[7] = data.x; + config_x[8] = data.y; + config_y[0] = data.z; + config_y[1] = data.w; + + data = read_node_float(kg, offset); + config_y[2] = data.x; + config_y[3] = data.y; + config_y[4] = data.z; + config_y[5] = data.w; + + data = read_node_float(kg, offset); + config_y[6] = data.x; + config_y[7] = data.y; + config_y[8] = data.z; + config_z[0] = data.w; + + data = read_node_float(kg, offset); + config_z[1] = data.x; + config_z[2] = data.y; + config_z[3] = data.z; + config_z[4] = data.w; + + data = read_node_float(kg, offset); + config_z[5] = data.x; + config_z[6] = data.y; + config_z[7] = data.z; + config_z[8] = data.w; + + /* Compute Sky */ + if (sky_model == 0) { + f = sky_radiance_preetham(kg, + dir, + sunphi, + suntheta, + radiance_x, + radiance_y, + radiance_z, + config_x, + config_y, + config_z); + } + else { + f = sky_radiance_hosek(kg, + dir, + sunphi, + suntheta, + radiance_x, + radiance_y, + radiance_z, + config_x, + config_y, + config_z); + } } + /* Nishita */ else { - f = sky_radiance_new(kg, - dir, - sunphi, - suntheta, - radiance_x, - radiance_y, - radiance_z, - config_x, - config_y, - config_z); + /* Define variables */ + float nishita_data[9]; + + float4 data = read_node_float(kg, offset); + nishita_data[0] = data.x; + nishita_data[1] = data.y; + nishita_data[2] = data.z; + nishita_data[3] = data.w; + + data = read_node_float(kg, offset); + nishita_data[4] = data.x; + nishita_data[5] = data.y; + nishita_data[6] = data.z; + nishita_data[7] = data.w; + + data = read_node_float(kg, offset); + nishita_data[8] = data.x; + uint texture_id = __float_as_uint(data.y); + + /* Compute Sky */ + f = sky_radiance_nishita(kg, dir, nishita_data, texture_id); } stack_store_float3(stack, out_offset, f); diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h index e913d9e0489..f1ebb37e23e 100644 --- a/intern/cycles/kernel/svm/svm_types.h +++ b/intern/cycles/kernel/svm/svm_types.h @@ -414,7 +414,7 @@ typedef enum NodeWaveProfile { NODE_WAVE_PROFILE_TRI, } NodeWaveProfile; -typedef enum NodeSkyType { NODE_SKY_OLD, NODE_SKY_NEW } NodeSkyType; +typedef enum NodeSkyType { NODE_SKY_PREETHAM, NODE_SKY_HOSEK, NODE_SKY_NISHITA } NodeSkyType; typedef enum NodeGradientType { NODE_BLEND_LINEAR, diff --git a/intern/cycles/render/CMakeLists.txt b/intern/cycles/render/CMakeLists.txt index 472b5a0c101..e37a0407976 100644 --- a/intern/cycles/render/CMakeLists.txt +++ b/intern/cycles/render/CMakeLists.txt @@ -24,6 +24,7 @@ set(SRC hair.cpp image.cpp image_oiio.cpp + image_sky.cpp image_vdb.cpp integrator.cpp jitter.cpp @@ -64,6 +65,7 @@ set(SRC_HEADERS hair.h image.h image_oiio.h + image_sky.h image_vdb.h integrator.h light.h diff --git a/intern/cycles/render/image_sky.cpp b/intern/cycles/render/image_sky.cpp new file mode 100644 index 00000000000..3e7b491f609 --- /dev/null +++ b/intern/cycles/render/image_sky.cpp @@ -0,0 +1,95 @@ +/* + * Copyright 2011-2020 Blender Foundation + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "render/image_sky.h" + +#include "util/util_image.h" +#include "util/util_logging.h" +#include "util/util_path.h" +#include "util/util_sky_model.h" + +CCL_NAMESPACE_BEGIN + +SkyLoader::SkyLoader( + float sun_elevation, int altitude, float air_density, float dust_density, float ozone_density) + : sun_elevation(sun_elevation), + altitude(altitude), + air_density(air_density), + dust_density(dust_density), + ozone_density(ozone_density) +{ +} + +SkyLoader::~SkyLoader(){}; + +bool SkyLoader::load_metadata(ImageMetaData &metadata) +{ + metadata.width = 512; + metadata.height = 128; + metadata.channels = 3; + metadata.depth = 1; + metadata.type = IMAGE_DATA_TYPE_FLOAT4; + metadata.compress_as_srgb = false; + return true; +} + +bool SkyLoader::load_pixels(const ImageMetaData &metadata, + void *pixels, + const size_t /*pixels_size*/, + const bool /*associate_alpha*/) +{ + /* definitions */ + int width = metadata.width; + int height = metadata.height; + float *pixel_data = (float *)pixels; + float altitude_f = (float)altitude; + + /* precompute sky texture */ + const int num_chunks = TaskScheduler::num_threads(); + const int chunk_size = height / num_chunks; + TaskPool pool; + for (int chunk = 0; chunk < num_chunks; chunk++) { + const int chunk_start = chunk * chunk_size; + const int chunk_end = (chunk + 1 < num_chunks) ? (chunk + 1) * chunk_size : height; + pool.push(function_bind(&nishita_skymodel_precompute_texture, + pixel_data, + metadata.channels, + chunk_start, + chunk_end, + width, + height, + sun_elevation, + altitude_f, + air_density, + dust_density, + ozone_density)); + } + pool.wait_work(); + + return true; +} + +string SkyLoader::name() const +{ + return "sky_nishita"; +} + +bool SkyLoader::equals(const ImageLoader & /*other*/) const +{ + return false; +} + +CCL_NAMESPACE_END diff --git a/intern/cycles/render/image_sky.h b/intern/cycles/render/image_sky.h new file mode 100644 index 00000000000..cf4a3e8942c --- /dev/null +++ b/intern/cycles/render/image_sky.h @@ -0,0 +1,49 @@ +/* + * Copyright 2011-2020 Blender Foundation + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "render/image.h" + +CCL_NAMESPACE_BEGIN + +class SkyLoader : public ImageLoader { + private: + float sun_elevation; + int altitude; + float air_density; + float dust_density; + float ozone_density; + + public: + SkyLoader(float sun_elevation, + int altitude, + float air_density, + float dust_density, + float ozone_density); + ~SkyLoader(); + + bool load_metadata(ImageMetaData &metadata) override; + + bool load_pixels(const ImageMetaData &metadata, + void *pixels, + const size_t /*pixels_size*/, + const bool /*associate_alpha*/) override; + + string name() const override; + + bool equals(const ImageLoader & /*other*/) const override; +}; + +CCL_NAMESPACE_END diff --git a/intern/cycles/render/light.cpp b/intern/cycles/render/light.cpp index cb7474017fa..225cedfef55 100644 --- a/intern/cycles/render/light.cpp +++ b/intern/cycles/render/light.cpp @@ -450,6 +450,7 @@ void LightManager::device_update_distribution(Device *, /* update device */ KernelIntegrator *kintegrator = &dscene->data.integrator; + KernelBackground *kbackground = &dscene->data.background; KernelFilm *kfilm = &dscene->data.film; kintegrator->use_direct_light = (totarea > 0.0f); @@ -493,15 +494,18 @@ void LightManager::device_update_distribution(Device *, /* Portals */ if (num_portals > 0) { - kintegrator->portal_offset = light_index; - kintegrator->num_portals = num_portals; - kintegrator->portal_pdf = background_mis ? 0.5f : 1.0f; + kbackground->portal_offset = light_index; + kbackground->num_portals = num_portals; + kbackground->portal_weight = 1.0f; } else { - kintegrator->num_portals = 0; - kintegrator->portal_offset = 0; - kintegrator->portal_pdf = 0.0f; + kbackground->num_portals = 0; + kbackground->portal_offset = 0; + kbackground->portal_weight = 0.0f; } + + /* Map */ + kbackground->map_weight = background_mis ? 1.0f : 0.0f; } else { dscene->light_distribution.free(); @@ -511,9 +515,12 @@ void LightManager::device_update_distribution(Device *, kintegrator->pdf_triangles = 0.0f; kintegrator->pdf_lights = 0.0f; kintegrator->use_lamp_mis = false; - kintegrator->num_portals = 0; - kintegrator->portal_offset = 0; - kintegrator->portal_pdf = 0.0f; + + kbackground->num_portals = 0; + kbackground->portal_offset = 0; + kbackground->portal_weight = 0.0f; + kbackground->sun_weight = 0.0f; + kbackground->map_weight = 0.0f; kfilm->pass_shadow_scale = 1.0f; } @@ -562,7 +569,7 @@ void LightManager::device_update_background(Device *device, Scene *scene, Progress &progress) { - KernelIntegrator *kintegrator = &dscene->data.integrator; + KernelBackground *kbackground = &dscene->data.background; Light *background_light = NULL; /* find background light */ @@ -575,31 +582,79 @@ void LightManager::device_update_background(Device *device, /* no background light found, signal renderer to skip sampling */ if (!background_light || !background_light->is_enabled) { - kintegrator->pdf_background_res_x = 0; - kintegrator->pdf_background_res_y = 0; + kbackground->map_res_x = 0; + kbackground->map_res_y = 0; + kbackground->map_weight = 0.0f; + kbackground->sun_weight = 0.0f; + kbackground->use_mis = (kbackground->portal_weight > 0.0f); return; } progress.set_status("Updating Lights", "Importance map"); - assert(kintegrator->use_direct_light); + assert(dscene->data.integrator.use_direct_light); + + int2 environment_res = make_int2(0, 0); + Shader *shader = scene->background->get_shader(scene); + int num_suns = 0; + foreach (ShaderNode *node, shader->graph->nodes) { + if (node->type == EnvironmentTextureNode::node_type) { + EnvironmentTextureNode *env = (EnvironmentTextureNode *)node; + ImageMetaData metadata; + if (!env->handle.empty()) { + ImageMetaData metadata = env->handle.metadata(); + environment_res.x = max(environment_res.x, metadata.width); + environment_res.y = max(environment_res.y, metadata.height); + } + } + if (node->type == SkyTextureNode::node_type) { + SkyTextureNode *sky = (SkyTextureNode *)node; + if (sky->type == NODE_SKY_NISHITA && sky->sun_disc) { + /* Ensure that the input coordinates aren't transformed before they reach the node. + * If that is the case, the logic used for sampling the sun's location does not work + * and we have to fall back to map-based sampling. */ + const ShaderInput *vec_in = sky->input("Vector"); + if (vec_in && vec_in->link && vec_in->link->parent) { + ShaderNode *vec_src = vec_in->link->parent; + if ((vec_src->type != TextureCoordinateNode::node_type) || + (vec_in->link != vec_src->output("Generated"))) { + environment_res.x = max(environment_res.x, 4096); + environment_res.y = max(environment_res.y, 2048); + continue; + } + } + + float latitude = sky->sun_elevation; + float longitude = M_2PI_F - sky->sun_rotation + M_PI_2_F; + float half_angle = sky->sun_size * 0.5f; + kbackground->sun = make_float4(cosf(latitude) * cosf(longitude), + cosf(latitude) * sinf(longitude), + sinf(latitude), + half_angle); + kbackground->sun_weight = 4.0f; + environment_res.x = max(environment_res.x, 512); + environment_res.y = max(environment_res.y, 256); + num_suns++; + } + } + } + + /* If there's more than one sun, fall back to map sampling instead. */ + if (num_suns != 1) { + kbackground->sun_weight = 0.0f; + environment_res.x = max(environment_res.x, 4096); + environment_res.y = max(environment_res.y, 2048); + } + + /* Enable MIS for background sampling if any strategy is active. */ + kbackground->use_mis = (kbackground->portal_weight + kbackground->map_weight + + kbackground->sun_weight) > 0.0f; /* get the resolution from the light's size (we stuff it in there) */ int2 res = make_int2(background_light->map_resolution, background_light->map_resolution / 2); /* If the resolution isn't set manually, try to find an environment texture. */ if (res.x == 0) { - Shader *shader = scene->background->get_shader(scene); - foreach (ShaderNode *node, shader->graph->nodes) { - if (node->type == EnvironmentTextureNode::node_type) { - EnvironmentTextureNode *env = (EnvironmentTextureNode *)node; - ImageMetaData metadata; - if (!env->handle.empty()) { - ImageMetaData metadata = env->handle.metadata(); - res.x = max(res.x, metadata.width); - res.y = max(res.y, metadata.height); - } - } - } + res = environment_res; if (res.x > 0 && res.y > 0) { VLOG(2) << "Automatically set World MIS resolution to " << res.x << " by " << res.y << "\n"; } @@ -609,8 +664,8 @@ void LightManager::device_update_background(Device *device, res = make_int2(1024, 512); VLOG(2) << "Setting World MIS resolution to default\n"; } - kintegrator->pdf_background_res_x = res.x; - kintegrator->pdf_background_res_y = res.y; + kbackground->map_res_x = res.x; + kbackground->map_res_y = res.y; vector<float3> pixels; shade_background_pixels(device, dscene, res.x, res.y, pixels, progress); diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp index cdcaeb246dd..ab392839e52 100644 --- a/intern/cycles/render/nodes.cpp +++ b/intern/cycles/render/nodes.cpp @@ -19,6 +19,7 @@ #include "render/constant_fold.h" #include "render/film.h" #include "render/image.h" +#include "render/image_sky.h" #include "render/integrator.h" #include "render/light.h" #include "render/mesh.h" @@ -630,7 +631,7 @@ typedef struct SunSky { /* Parameter */ float radiance_x, radiance_y, radiance_z; - float config_x[9], config_y[9], config_z[9]; + float config_x[9], config_y[9], config_z[9], nishita_data[9]; } SunSky; /* Preetham model */ @@ -640,7 +641,7 @@ static float sky_perez_function(float lam[6], float theta, float gamma) (1.0f + lam[2] * expf(lam[3] * gamma) + lam[4] * cosf(gamma) * cosf(gamma)); } -static void sky_texture_precompute_old(SunSky *sunsky, float3 dir, float turbidity) +static void sky_texture_precompute_preetham(SunSky *sunsky, float3 dir, float turbidity) { /* * We re-use the SunSky struct of the new model, to avoid extra variables @@ -703,10 +704,10 @@ static void sky_texture_precompute_old(SunSky *sunsky, float3 dir, float turbidi } /* Hosek / Wilkie */ -static void sky_texture_precompute_new(SunSky *sunsky, - float3 dir, - float turbidity, - float ground_albedo) +static void sky_texture_precompute_hosek(SunSky *sunsky, + float3 dir, + float turbidity, + float ground_albedo) { /* Calculate Sun Direction and save coordinates */ float2 spherical = sky_spherical_coordinates(dir); @@ -743,6 +744,34 @@ static void sky_texture_precompute_new(SunSky *sunsky, arhosekskymodelstate_free(sky_state); } +/* Nishita improved */ +static void sky_texture_precompute_nishita(SunSky *sunsky, + bool sun_disc, + float sun_size, + float sun_elevation, + float sun_rotation, + int altitude, + float air_density, + float dust_density) +{ + /* sample 2 sun pixels */ + float pixel_bottom[3]; + float pixel_top[3]; + float altitude_f = (float)altitude; + nishita_skymodel_precompute_sun( + sun_elevation, sun_size, altitude_f, air_density, dust_density, pixel_bottom, pixel_top); + /* send data to svm_sky */ + sunsky->nishita_data[0] = pixel_bottom[0]; + sunsky->nishita_data[1] = pixel_bottom[1]; + sunsky->nishita_data[2] = pixel_bottom[2]; + sunsky->nishita_data[3] = pixel_top[0]; + sunsky->nishita_data[4] = pixel_top[1]; + sunsky->nishita_data[5] = pixel_top[2]; + sunsky->nishita_data[6] = sun_elevation; + sunsky->nishita_data[7] = M_2PI_F - sun_rotation; + sunsky->nishita_data[8] = sun_disc ? sun_size : 0.0f; +} + NODE_DEFINE(SkyTextureNode) { NodeType *type = NodeType::add("sky_texture", create, NodeType::SHADER); @@ -750,13 +779,22 @@ NODE_DEFINE(SkyTextureNode) TEXTURE_MAPPING_DEFINE(SkyTextureNode); static NodeEnum type_enum; - type_enum.insert("preetham", NODE_SKY_OLD); - type_enum.insert("hosek_wilkie", NODE_SKY_NEW); - SOCKET_ENUM(type, "Type", type_enum, NODE_SKY_NEW); + type_enum.insert("preetham", NODE_SKY_PREETHAM); + type_enum.insert("hosek_wilkie", NODE_SKY_HOSEK); + type_enum.insert("nishita_improved", NODE_SKY_NISHITA); + SOCKET_ENUM(type, "Type", type_enum, NODE_SKY_NISHITA); SOCKET_VECTOR(sun_direction, "Sun Direction", make_float3(0.0f, 0.0f, 1.0f)); SOCKET_FLOAT(turbidity, "Turbidity", 2.2f); SOCKET_FLOAT(ground_albedo, "Ground Albedo", 0.3f); + SOCKET_BOOLEAN(sun_disc, "Sun Disc", true); + SOCKET_FLOAT(sun_size, "Sun Size", 0.009512f); + SOCKET_FLOAT(sun_elevation, "Sun Elevation", M_PI_2_F); + SOCKET_FLOAT(sun_rotation, "Sun Rotation", 0.0f); + SOCKET_INT(altitude, "Altitude", 0); + SOCKET_FLOAT(air_density, "Air", 1.0f); + SOCKET_FLOAT(dust_density, "Dust", 1.0f); + SOCKET_FLOAT(ozone_density, "Ozone", 1.0f); SOCKET_IN_POINT( vector, "Vector", make_float3(0.0f, 0.0f, 0.0f), SocketType::LINK_TEXTURE_GENERATED); @@ -776,10 +814,32 @@ void SkyTextureNode::compile(SVMCompiler &compiler) ShaderOutput *color_out = output("Color"); SunSky sunsky; - if (type == NODE_SKY_OLD) - sky_texture_precompute_old(&sunsky, sun_direction, turbidity); - else if (type == NODE_SKY_NEW) - sky_texture_precompute_new(&sunsky, sun_direction, turbidity, ground_albedo); + if (type == NODE_SKY_PREETHAM) + sky_texture_precompute_preetham(&sunsky, sun_direction, turbidity); + else if (type == NODE_SKY_HOSEK) + sky_texture_precompute_hosek(&sunsky, sun_direction, turbidity, ground_albedo); + else if (type == NODE_SKY_NISHITA) { + sky_texture_precompute_nishita(&sunsky, + sun_disc, + sun_size, + sun_elevation, + sun_rotation, + altitude, + air_density, + dust_density); + /* precomputed texture image parameters */ + ImageManager *image_manager = compiler.scene->image_manager; + ImageParams impar; + impar.interpolation = INTERPOLATION_LINEAR; + impar.extension = EXTENSION_EXTEND; + + /* precompute sky texture */ + if (handle.empty()) { + SkyLoader *loader = new SkyLoader( + sun_elevation, altitude, air_density, dust_density, ozone_density); + handle = image_manager->add_image(loader, impar); + } + } else assert(false); @@ -787,38 +847,52 @@ void SkyTextureNode::compile(SVMCompiler &compiler) compiler.stack_assign(color_out); compiler.add_node(NODE_TEX_SKY, vector_offset, compiler.stack_assign(color_out), type); - compiler.add_node(__float_as_uint(sunsky.phi), - __float_as_uint(sunsky.theta), - __float_as_uint(sunsky.radiance_x), - __float_as_uint(sunsky.radiance_y)); - compiler.add_node(__float_as_uint(sunsky.radiance_z), - __float_as_uint(sunsky.config_x[0]), - __float_as_uint(sunsky.config_x[1]), - __float_as_uint(sunsky.config_x[2])); - compiler.add_node(__float_as_uint(sunsky.config_x[3]), - __float_as_uint(sunsky.config_x[4]), - __float_as_uint(sunsky.config_x[5]), - __float_as_uint(sunsky.config_x[6])); - compiler.add_node(__float_as_uint(sunsky.config_x[7]), - __float_as_uint(sunsky.config_x[8]), - __float_as_uint(sunsky.config_y[0]), - __float_as_uint(sunsky.config_y[1])); - compiler.add_node(__float_as_uint(sunsky.config_y[2]), - __float_as_uint(sunsky.config_y[3]), - __float_as_uint(sunsky.config_y[4]), - __float_as_uint(sunsky.config_y[5])); - compiler.add_node(__float_as_uint(sunsky.config_y[6]), - __float_as_uint(sunsky.config_y[7]), - __float_as_uint(sunsky.config_y[8]), - __float_as_uint(sunsky.config_z[0])); - compiler.add_node(__float_as_uint(sunsky.config_z[1]), - __float_as_uint(sunsky.config_z[2]), - __float_as_uint(sunsky.config_z[3]), - __float_as_uint(sunsky.config_z[4])); - compiler.add_node(__float_as_uint(sunsky.config_z[5]), - __float_as_uint(sunsky.config_z[6]), - __float_as_uint(sunsky.config_z[7]), - __float_as_uint(sunsky.config_z[8])); + /* nishita doesn't need this data */ + if (type != NODE_SKY_NISHITA) { + compiler.add_node(__float_as_uint(sunsky.phi), + __float_as_uint(sunsky.theta), + __float_as_uint(sunsky.radiance_x), + __float_as_uint(sunsky.radiance_y)); + compiler.add_node(__float_as_uint(sunsky.radiance_z), + __float_as_uint(sunsky.config_x[0]), + __float_as_uint(sunsky.config_x[1]), + __float_as_uint(sunsky.config_x[2])); + compiler.add_node(__float_as_uint(sunsky.config_x[3]), + __float_as_uint(sunsky.config_x[4]), + __float_as_uint(sunsky.config_x[5]), + __float_as_uint(sunsky.config_x[6])); + compiler.add_node(__float_as_uint(sunsky.config_x[7]), + __float_as_uint(sunsky.config_x[8]), + __float_as_uint(sunsky.config_y[0]), + __float_as_uint(sunsky.config_y[1])); + compiler.add_node(__float_as_uint(sunsky.config_y[2]), + __float_as_uint(sunsky.config_y[3]), + __float_as_uint(sunsky.config_y[4]), + __float_as_uint(sunsky.config_y[5])); + compiler.add_node(__float_as_uint(sunsky.config_y[6]), + __float_as_uint(sunsky.config_y[7]), + __float_as_uint(sunsky.config_y[8]), + __float_as_uint(sunsky.config_z[0])); + compiler.add_node(__float_as_uint(sunsky.config_z[1]), + __float_as_uint(sunsky.config_z[2]), + __float_as_uint(sunsky.config_z[3]), + __float_as_uint(sunsky.config_z[4])); + compiler.add_node(__float_as_uint(sunsky.config_z[5]), + __float_as_uint(sunsky.config_z[6]), + __float_as_uint(sunsky.config_z[7]), + __float_as_uint(sunsky.config_z[8])); + } + else { + compiler.add_node(__float_as_uint(sunsky.nishita_data[0]), + __float_as_uint(sunsky.nishita_data[1]), + __float_as_uint(sunsky.nishita_data[2]), + __float_as_uint(sunsky.nishita_data[3])); + compiler.add_node(__float_as_uint(sunsky.nishita_data[4]), + __float_as_uint(sunsky.nishita_data[5]), + __float_as_uint(sunsky.nishita_data[6]), + __float_as_uint(sunsky.nishita_data[7])); + compiler.add_node(__float_as_uint(sunsky.nishita_data[8]), handle.svm_slot(), 0, 0); + } tex_mapping.compile_end(compiler, vector_in, vector_offset); } @@ -828,10 +902,32 @@ void SkyTextureNode::compile(OSLCompiler &compiler) tex_mapping.compile(compiler); SunSky sunsky; - if (type == NODE_SKY_OLD) - sky_texture_precompute_old(&sunsky, sun_direction, turbidity); - else if (type == NODE_SKY_NEW) - sky_texture_precompute_new(&sunsky, sun_direction, turbidity, ground_albedo); + if (type == NODE_SKY_PREETHAM) + sky_texture_precompute_preetham(&sunsky, sun_direction, turbidity); + else if (type == NODE_SKY_HOSEK) + sky_texture_precompute_hosek(&sunsky, sun_direction, turbidity, ground_albedo); + else if (type == NODE_SKY_NISHITA) { + sky_texture_precompute_nishita(&sunsky, + sun_disc, + sun_size, + sun_elevation, + sun_rotation, + altitude, + air_density, + dust_density); + /* precomputed texture image parameters */ + ImageManager *image_manager = compiler.scene->image_manager; + ImageParams impar; + impar.interpolation = INTERPOLATION_LINEAR; + impar.extension = EXTENSION_EXTEND; + + /* precompute sky texture */ + if (handle.empty()) { + SkyLoader *loader = new SkyLoader( + sun_elevation, altitude, air_density, dust_density, ozone_density); + handle = image_manager->add_image(loader, impar); + } + } else assert(false); @@ -843,6 +939,11 @@ void SkyTextureNode::compile(OSLCompiler &compiler) compiler.parameter_array("config_x", sunsky.config_x, 9); compiler.parameter_array("config_y", sunsky.config_y, 9); compiler.parameter_array("config_z", sunsky.config_z, 9); + compiler.parameter_array("nishita_data", sunsky.nishita_data, 9); + /* nishita texture */ + if (type == NODE_SKY_NISHITA) { + compiler.parameter_texture("filename", handle.svm_slot()); + } compiler.add(this, "node_sky_texture"); } diff --git a/intern/cycles/render/nodes.h b/intern/cycles/render/nodes.h index 83c3ad071ae..846ba7423e5 100644 --- a/intern/cycles/render/nodes.h +++ b/intern/cycles/render/nodes.h @@ -168,7 +168,16 @@ class SkyTextureNode : public TextureNode { float3 sun_direction; float turbidity; float ground_albedo; + bool sun_disc; + float sun_size; + float sun_elevation; + float sun_rotation; + int altitude; + float air_density; + float dust_density; + float ozone_density; float3 vector; + ImageHandle handle; }; class OutputNode : public ShaderNode { diff --git a/intern/cycles/util/CMakeLists.txt b/intern/cycles/util/CMakeLists.txt index c1f71461dfd..16d47d57e69 100644 --- a/intern/cycles/util/CMakeLists.txt +++ b/intern/cycles/util/CMakeLists.txt @@ -100,6 +100,7 @@ set(SRC_HEADERS util_sky_model.cpp util_sky_model.h util_sky_model_data.h + util_sky_nishita.cpp util_avxf.h util_avxb.h util_semaphore.h diff --git a/intern/cycles/util/util_sky_model.h b/intern/cycles/util/util_sky_model.h index 84340614b2c..36f1079a16d 100644 --- a/intern/cycles/util/util_sky_model.h +++ b/intern/cycles/util/util_sky_model.h @@ -298,6 +298,8 @@ HINT #1: if you want to model the sky of an earth-like planet that orbits previous paragraph. */ +#include "util/util_types.h" + CCL_NAMESPACE_BEGIN #ifndef _SKY_MODEL_H_ @@ -426,4 +428,26 @@ double arhosekskymodel_solar_radiance(ArHosekSkyModelState *state, #endif // _SKY_MODEL_H_ +/* Nishita improved sky model */ + +void nishita_skymodel_precompute_texture(float *pixels, + int stride, + int start_y, + int end_y, + int width, + int height, + float sun_elevation, + float altitude, + float air_density, + float dust_density, + float ozone_density); + +void nishita_skymodel_precompute_sun(float sun_elevation, + float angular_diameter, + float altitude, + float air_density, + float dust_density, + float *pixel_bottom, + float *pixel_top); + CCL_NAMESPACE_END diff --git a/intern/cycles/util/util_sky_nishita.cpp b/intern/cycles/util/util_sky_nishita.cpp new file mode 100644 index 00000000000..92397804d43 --- /dev/null +++ b/intern/cycles/util/util_sky_nishita.cpp @@ -0,0 +1,371 @@ +/* + * Copyright 2011-2020 Blender Foundation + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "util/util_math.h" +#include "util/util_sky_model.h" + +CCL_NAMESPACE_BEGIN + +/* Constants */ +static const float rayleigh_scale = 8000.0f; // Rayleigh scale height (m) +static const float mie_scale = 1200.0f; // Mie scale height (m) +static const float mie_coeff = 2e-5f; // Mie scattering coefficient +static const float mie_G = 0.76f; // aerosols anisotropy +static const float earth_radius = 6360000.0f; // radius of Earth (m) +static const float atmosphere_radius = 6420000.0f; // radius of atmosphere (m) +static const int steps = 32; // segments per primary ray +static const int steps_light = 16; // segments per sun connection ray +static const int num_wavelengths = 21; // number of wavelengths +/* irradiance at top of atmosphere */ +static const float irradiance[] = { + 1.45756829855592995315f, 1.56596305559738380175f, 1.65148449067670455293f, + 1.71496242737209314555f, 1.75797983805020541226f, 1.78256407885924539336f, + 1.79095108475838560302f, 1.78541550133410664714f, 1.76815554864306845317f, + 1.74122069647250410362f, 1.70647127164943679389f, 1.66556087452739887134f, + 1.61993437242451854274f, 1.57083597368892080581f, 1.51932335059305478886f, + 1.46628494965214395407f, 1.41245852740172450623f, 1.35844961970384092709f, + 1.30474913844739281998f, 1.25174963272610817455f, 1.19975998755420620867f}; +/* Rayleigh scattering coefficient */ +static const float rayleigh_coeff[] = { + 0.00005424820087636473f, 0.00004418549866505454f, 0.00003635151910165377f, + 0.00003017929012024763f, 0.00002526320226989157f, 0.00002130859310621843f, + 0.00001809838025320633f, 0.00001547057129129042f, 0.00001330284977336850f, + 0.00001150184784075764f, 0.00000999557429990163f, 0.00000872799973630707f, + 0.00000765513700977967f, 0.00000674217203751443f, 0.00000596134125832052f, + 0.00000529034598065810f, 0.00000471115687557433f, 0.00000420910481110487f, + 0.00000377218381260133f, 0.00000339051255477280f, 0.00000305591531679811f}; +/* Ozone absorption coefficient */ +static const float ozone_coeff[] = { + 0.00000000325126849861f, 0.00000000585395365047f, 0.00000001977191155085f, + 0.00000007309568762914f, 0.00000020084561514287f, 0.00000040383958096161f, + 0.00000063551335912363f, 0.00000096707041180970f, 0.00000154797400424410f, + 0.00000209038647223331f, 0.00000246128056164565f, 0.00000273551299461512f, + 0.00000215125863128643f, 0.00000159051840791988f, 0.00000112356197979857f, + 0.00000073527551487574f, 0.00000046450130357806f, 0.00000033096079921048f, + 0.00000022512612292678f, 0.00000014879129266490f, 0.00000016828623364192f}; +/* CIE XYZ color matching functions */ +static const float cmf_xyz[][3] = {{0.00136800000f, 0.00003900000f, 0.00645000100f}, + {0.01431000000f, 0.00039600000f, 0.06785001000f}, + {0.13438000000f, 0.00400000000f, 0.64560000000f}, + {0.34828000000f, 0.02300000000f, 1.74706000000f}, + {0.29080000000f, 0.06000000000f, 1.66920000000f}, + {0.09564000000f, 0.13902000000f, 0.81295010000f}, + {0.00490000000f, 0.32300000000f, 0.27200000000f}, + {0.06327000000f, 0.71000000000f, 0.07824999000f}, + {0.29040000000f, 0.95400000000f, 0.02030000000f}, + {0.59450000000f, 0.99500000000f, 0.00390000000f}, + {0.91630000000f, 0.87000000000f, 0.00165000100f}, + {1.06220000000f, 0.63100000000f, 0.00080000000f}, + {0.85444990000f, 0.38100000000f, 0.00019000000f}, + {0.44790000000f, 0.17500000000f, 0.00002000000f}, + {0.16490000000f, 0.06100000000f, 0.00000000000f}, + {0.04677000000f, 0.01700000000f, 0.00000000000f}, + {0.01135916000f, 0.00410200000f, 0.00000000000f}, + {0.00289932700f, 0.00104700000f, 0.00000000000f}, + {0.00069007860f, 0.00024920000f, 0.00000000000f}, + {0.00016615050f, 0.00006000000f, 0.00000000000f}, + {0.00004150994f, 0.00001499000f, 0.00000000000f}}; + +static float3 geographical_to_direction(float lat, float lon) +{ + return make_float3(cosf(lat) * cosf(lon), cosf(lat) * sinf(lon), sinf(lat)); +} + +static float3 spec_to_xyz(float *spectrum) +{ + float3 xyz = make_float3(0.0f, 0.0f, 0.0f); + for (int i = 0; i < num_wavelengths; i++) { + xyz.x += cmf_xyz[i][0] * spectrum[i]; + xyz.y += cmf_xyz[i][1] * spectrum[i]; + xyz.z += cmf_xyz[i][2] * spectrum[i]; + } + return xyz * (20 * 683 * 1e-9f); +} + +/* Atmosphere volume models */ + +static float density_rayleigh(float height) +{ + return expf(-height / rayleigh_scale); +} + +static float density_mie(float height) +{ + return expf(-height / mie_scale); +} + +static float density_ozone(float height) +{ + float den = 0.0f; + if (height >= 10000.0f && height < 25000.0f) + den = 1.0f / 15000.0f * height - 2.0f / 3.0f; + else if (height >= 25000 && height < 40000) + den = -(1.0f / 15000.0f * height - 8.0f / 3.0f); + return den; +} + +static float phase_rayleigh(float mu) +{ + return 3.0f / (16.0f * M_PI_F) * (1.0f + sqr(mu)); +} + +static float phase_mie(float mu) +{ + static const float sqr_G = mie_G * mie_G; + + return (3.0f * (1.0f - sqr_G) * (1.0f + sqr(mu))) / + (8.0f * M_PI_F * (2.0f + sqr_G) * powf((1.0f + sqr_G - 2.0f * mie_G * mu), 1.5)); +} + +/* Intersection helpers */ +static bool surface_intersection(float3 pos, float3 dir) +{ + if (dir.z >= 0) + return false; + float t = dot(dir, -pos) / len_squared(dir); + float D = pos.x * pos.x - 2.0f * (-pos.x) * dir.x * t + dir.x * t * dir.x * t + pos.y * pos.y - + 2.0f * (-pos.y) * dir.y * t + (dir.y * t) * (dir.y * t) + pos.z * pos.z - + 2.0f * (-pos.z) * dir.z * t + dir.z * t * dir.z * t; + return (D <= sqr(earth_radius)); +} + +static float3 atmosphere_intersection(float3 pos, float3 dir) +{ + float b = -2.0f * dot(dir, -pos); + float c = len_squared(pos) - sqr(atmosphere_radius); + float t = (-b + sqrtf(b * b - 4.0f * c)) / 2.0f; + return make_float3(pos.x + dir.x * t, pos.y + dir.y * t, pos.z + dir.z * t); +} + +static float3 ray_optical_depth(float3 ray_origin, float3 ray_dir) +{ + /* This code computes the optical depth along a ray through the atmosphere. */ + float3 ray_end = atmosphere_intersection(ray_origin, ray_dir); + float ray_length = distance(ray_origin, ray_end); + + /* To compute the optical depth, we step along the ray in segments and + * accumulate the optical depth along each segment. */ + float segment_length = ray_length / steps_light; + float3 segment = segment_length * ray_dir; + + /* Instead of tracking the transmission spectrum across all wavelengths directly, + * we use the fact that the density always has the same spectrum for each type of + * scattering, so we split the density into a constant spectrum and a factor and + * only track the factors. */ + float3 optical_depth = make_float3(0.0f, 0.0f, 0.0f); + + /* The density of each segment is evaluated at its middle. */ + float3 P = ray_origin + 0.5f * segment; + for (int i = 0; i < steps_light; i++) { + /* Compute height above sea level. */ + float height = len(P) - earth_radius; + + /* Accumulate optical depth of this segment (density is assumed to be constant along it). */ + float3 density = make_float3( + density_rayleigh(height), density_mie(height), density_ozone(height)); + optical_depth += segment_length * density; + + /* Advance along ray. */ + P += segment; + } + + return optical_depth; +} + +/* Single Scattering implementation */ +static void single_scattering(float3 ray_dir, + float3 sun_dir, + float3 ray_origin, + float air_density, + float dust_density, + float ozone_density, + float *r_spectrum) +{ + /* This code computes single-inscattering along a ray through the atmosphere. */ + float3 ray_end = atmosphere_intersection(ray_origin, ray_dir); + float ray_length = distance(ray_origin, ray_end); + + /* To compute the inscattering, we step along the ray in segments and accumulate + * the inscattering as well as the optical depth along each segment. */ + float segment_length = ray_length / steps; + float3 segment = segment_length * ray_dir; + + /* Instead of tracking the transmission spectrum across all wavelengths directly, + * we use the fact that the density always has the same spectrum for each type of + * scattering, so we split the density into a constant spectrum and a factor and + * only track the factors. */ + float3 optical_depth = make_float3(0.0f, 0.0f, 0.0f); + + /* Zero out light accumulation. */ + for (int wl = 0; wl < num_wavelengths; wl++) { + r_spectrum[wl] = 0.0f; + } + + /* Compute phase function for scattering and the density scale factor. */ + float mu = dot(ray_dir, sun_dir); + float3 phase_function = make_float3(phase_rayleigh(mu), phase_mie(mu), 0.0f); + float3 density_scale = make_float3(air_density, dust_density, ozone_density); + + /* The density and in-scattering of each segment is evaluated at its middle. */ + float3 P = ray_origin + 0.5f * segment; + for (int i = 0; i < steps; i++) { + /* Compute height above sea level. */ + float height = len(P) - earth_radius; + + /* Evaluate and accumulate optical depth along the ray. */ + float3 density = density_scale * make_float3(density_rayleigh(height), + density_mie(height), + density_ozone(height)); + optical_depth += segment_length * density; + + /* If the earth isn't in the way, evaluate inscattering from the sun. */ + if (!surface_intersection(P, sun_dir)) { + float3 light_optical_depth = density_scale * ray_optical_depth(P, sun_dir); + float3 total_optical_depth = optical_depth + light_optical_depth; + + /* attenuation of light */ + for (int wl = 0; wl < num_wavelengths; wl++) { + float3 extinction_density = total_optical_depth * make_float3(rayleigh_coeff[wl], + 1.11f * mie_coeff, + ozone_coeff[wl]); + float attenuation = expf(-reduce_add(extinction_density)); + + float3 scattering_density = density * make_float3(rayleigh_coeff[wl], mie_coeff, 0.0f); + + /* The total inscattered radiance from one segment is: + * Tr(A<->B) * Tr(B<->C) * sigma_s * phase * L * segment_length + * + * These terms are: + * Tr(A<->B): Transmission from start to scattering position (tracked in optical_depth) + * Tr(B<->C): Transmission from scattering position to light (computed in + * ray_optical_depth) sigma_s: Scattering density phase: Phase function of the scattering + * type (Rayleigh or Mie) L: Radiance coming from the light source segment_length: The + * length of the segment + * + * The code here is just that, with a bit of additional optimization to not store full + * spectra for the optical depth. + */ + r_spectrum[wl] += attenuation * reduce_add(phase_function * scattering_density) * + irradiance[wl] * segment_length; + } + } + + /* Advance along ray. */ + P += segment; + } +} + +/* calculate texture array */ +void nishita_skymodel_precompute_texture(float *pixels, + int stride, + int start_y, + int end_y, + int width, + int height, + float sun_elevation, + float altitude, + float air_density, + float dust_density, + float ozone_density) +{ + /* calculate texture pixels */ + float spectrum[num_wavelengths]; + int half_width = width / 2; + float3 cam_pos = make_float3(0, 0, earth_radius + altitude); + float3 sun_dir = geographical_to_direction(sun_elevation, 0.0f); + + float latitude_step = M_PI_2_F / height; + float longitude_step = M_2PI_F / width; + + for (int y = start_y; y < end_y; y++) { + float latitude = latitude_step * y; + + float *pixel_row = pixels + (y * width) * stride; + for (int x = 0; x < half_width; x++) { + float longitude = longitude_step * x - M_PI_F; + + float3 dir = geographical_to_direction(latitude, longitude); + single_scattering(dir, sun_dir, cam_pos, air_density, dust_density, ozone_density, spectrum); + float3 xyz = spec_to_xyz(spectrum); + + pixel_row[x * stride + 0] = xyz.x; + pixel_row[x * stride + 1] = xyz.y; + pixel_row[x * stride + 2] = xyz.z; + int mirror_x = width - x - 1; + pixel_row[mirror_x * stride + 0] = xyz.x; + pixel_row[mirror_x * stride + 1] = xyz.y; + pixel_row[mirror_x * stride + 2] = xyz.z; + } + } +} + +/* Sun disc */ +static void sun_radiation(float3 cam_dir, + float altitude, + float air_density, + float dust_density, + float solid_angle, + float *r_spectrum) +{ + float3 cam_pos = make_float3(0, 0, earth_radius + altitude); + float3 optical_depth = ray_optical_depth(cam_pos, cam_dir); + + /* Compute final spectrum. */ + for (int i = 0; i < num_wavelengths; i++) { + /* Combine spectra and the optical depth into transmittance. */ + float transmittance = rayleigh_coeff[i] * optical_depth.x * air_density + + 1.11f * mie_coeff * optical_depth.y * dust_density; + r_spectrum[i] = (irradiance[i] / solid_angle) * expf(-transmittance); + } +} + +void nishita_skymodel_precompute_sun(float sun_elevation, + float angular_diameter, + float altitude, + float air_density, + float dust_density, + float *pixel_bottom, + float *pixel_top) +{ + /* definitions */ + float half_angular = angular_diameter / 2.0f; + float solid_angle = M_2PI_F * (1.0f - cosf(half_angular)); + float spectrum[num_wavelengths]; + float bottom = sun_elevation - half_angular; + float top = sun_elevation + half_angular; + float elevation_bottom, elevation_top; + float3 pix_bottom, pix_top, sun_dir; + + /* compute 2 pixels for sun disc */ + elevation_bottom = (bottom > 0.0f) ? bottom : 0.0f; + elevation_top = (top > 0.0f) ? top : 0.0f; + sun_dir = geographical_to_direction(elevation_bottom, 0.0f); + sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum); + pix_bottom = spec_to_xyz(spectrum); + sun_dir = geographical_to_direction(elevation_top, 0.0f); + sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum); + pix_top = spec_to_xyz(spectrum); + + /* store pixels */ + pixel_bottom[0] = pix_bottom.x; + pixel_bottom[1] = pix_bottom.y; + pixel_bottom[2] = pix_bottom.z; + pixel_top[0] = pix_top.x; + pixel_top[1] = pix_top.y; + pixel_top[2] = pix_top.z; +} + +CCL_NAMESPACE_END diff --git a/source/blender/editors/space_node/drawnode.c b/source/blender/editors/space_node/drawnode.c index 01ac3a80871..ab50158c8bc 100644 --- a/source/blender/editors/space_node/drawnode.c +++ b/source/blender/editors/space_node/drawnode.c @@ -849,12 +849,35 @@ static void node_shader_buts_tex_environment_ex(uiLayout *layout, bContext *C, P static void node_shader_buts_tex_sky(uiLayout *layout, bContext *UNUSED(C), PointerRNA *ptr) { uiItemR(layout, ptr, "sky_type", DEFAULT_FLAGS, "", ICON_NONE); - uiItemR(layout, ptr, "sun_direction", DEFAULT_FLAGS, "", ICON_NONE); - uiItemR(layout, ptr, "turbidity", DEFAULT_FLAGS, NULL, ICON_NONE); - if (RNA_enum_get(ptr, "sky_type") == SHD_SKY_NEW) { + if (RNA_enum_get(ptr, "sky_type") == SHD_SKY_PREETHAM) { + uiItemR(layout, ptr, "sun_direction", DEFAULT_FLAGS, "", ICON_NONE); + uiItemR(layout, ptr, "turbidity", DEFAULT_FLAGS, NULL, ICON_NONE); + } + if (RNA_enum_get(ptr, "sky_type") == SHD_SKY_HOSEK) { + uiItemR(layout, ptr, "sun_direction", DEFAULT_FLAGS, "", ICON_NONE); + uiItemR(layout, ptr, "turbidity", DEFAULT_FLAGS, NULL, ICON_NONE); uiItemR(layout, ptr, "ground_albedo", DEFAULT_FLAGS, NULL, ICON_NONE); } + if (RNA_enum_get(ptr, "sky_type") == SHD_SKY_NISHITA) { + uiItemR(layout, ptr, "sun_disc", DEFAULT_FLAGS, NULL, 0); + + if (RNA_boolean_get(ptr, "sun_disc")) { + uiItemR(layout, ptr, "sun_size", DEFAULT_FLAGS, NULL, ICON_NONE); + } + + uiLayout *col; + col = uiLayoutColumn(layout, true); + uiItemR(col, ptr, "sun_elevation", DEFAULT_FLAGS, NULL, ICON_NONE); + uiItemR(col, ptr, "sun_rotation", DEFAULT_FLAGS, NULL, ICON_NONE); + + uiItemR(layout, ptr, "altitude", DEFAULT_FLAGS, NULL, ICON_NONE); + + col = uiLayoutColumn(layout, true); + uiItemR(col, ptr, "air_density", DEFAULT_FLAGS, NULL, ICON_NONE); + uiItemR(col, ptr, "dust_density", DEFAULT_FLAGS, NULL, ICON_NONE); + uiItemR(col, ptr, "ozone_density", DEFAULT_FLAGS, NULL, ICON_NONE); + } } static void node_shader_buts_tex_gradient(uiLayout *layout, bContext *UNUSED(C), PointerRNA *ptr) diff --git a/source/blender/makesdna/DNA_node_types.h b/source/blender/makesdna/DNA_node_types.h index 4ff0e531168..993aad92564 100644 --- a/source/blender/makesdna/DNA_node_types.h +++ b/source/blender/makesdna/DNA_node_types.h @@ -844,6 +844,15 @@ typedef struct NodeTexSky { float sun_direction[3]; float turbidity; float ground_albedo; + float sun_size; + float sun_elevation; + float sun_rotation; + int altitude; + float air_density; + float dust_density; + float ozone_density; + char sun_disc; + char _pad[3]; } NodeTexSky; typedef struct NodeTexImage { @@ -1171,8 +1180,9 @@ enum { }; /* sky texture */ -#define SHD_SKY_OLD 0 -#define SHD_SKY_NEW 1 +#define SHD_SKY_PREETHAM 0 +#define SHD_SKY_HOSEK 1 +#define SHD_SKY_NISHITA 2 /* environment texture */ #define SHD_PROJ_EQUIRECTANGULAR 0 diff --git a/source/blender/makesrna/intern/rna_nodetree.c b/source/blender/makesrna/intern/rna_nodetree.c index 32999c91fad..13dee92c1cb 100644 --- a/source/blender/makesrna/intern/rna_nodetree.c +++ b/source/blender/makesrna/intern/rna_nodetree.c @@ -4405,8 +4405,9 @@ static void def_sh_tex(StructRNA *srna) static void def_sh_tex_sky(StructRNA *srna) { static const EnumPropertyItem prop_sky_type[] = { - {SHD_SKY_OLD, "PREETHAM", 0, "Preetham", ""}, - {SHD_SKY_NEW, "HOSEK_WILKIE", 0, "Hosek / Wilkie", ""}, + {SHD_SKY_PREETHAM, "PREETHAM", 0, "Preetham", "Preetham 1999"}, + {SHD_SKY_HOSEK, "HOSEK_WILKIE", 0, "Hosek / Wilkie", "Hosek / Wilkie 2012"}, + {SHD_SKY_NISHITA, "NISHITA", 0, "Nishita", "Nishita 1993 improved"}, {0, NULL, 0, NULL, NULL}, }; static float default_dir[3] = {0.0f, 0.0f, 1.0f}; @@ -4419,7 +4420,7 @@ static void def_sh_tex_sky(StructRNA *srna) prop = RNA_def_property(srna, "sky_type", PROP_ENUM, PROP_NONE); RNA_def_property_enum_sdna(prop, NULL, "sky_model"); RNA_def_property_enum_items(prop, prop_sky_type); - RNA_def_property_ui_text(prop, "Sky Type", ""); + RNA_def_property_ui_text(prop, "Sky Type", "Which sky model should be used"); RNA_def_property_update(prop, 0, "rna_Node_update"); prop = RNA_def_property(srna, "sun_direction", PROP_FLOAT, PROP_DIRECTION); @@ -4439,6 +4440,54 @@ static void def_sh_tex_sky(StructRNA *srna) RNA_def_property_ui_text( prop, "Ground Albedo", "Ground color that is subtly reflected in the sky"); RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "sun_disc", PROP_BOOLEAN, PROP_NONE); + RNA_def_property_ui_text(prop, "Sun Disc", "Include the sun itself in the output"); + RNA_def_property_boolean_sdna(prop, NULL, "sun_disc", 1); + RNA_def_property_boolean_default(prop, true); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "sun_size", PROP_FLOAT, PROP_ANGLE); + RNA_def_property_ui_text(prop, "Sun Size", "Size of sun disc (angular diameter)"); + RNA_def_property_range(prop, 0.0f, M_PI_2); + RNA_def_property_float_default(prop, DEG2RADF(0.545)); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "sun_elevation", PROP_FLOAT, PROP_ANGLE); + RNA_def_property_ui_text(prop, "Sun Elevation", "Angle between sun and horizon"); + RNA_def_property_range(prop, -M_PI_2, M_PI_2); + RNA_def_property_float_default(prop, M_PI_2); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "sun_rotation", PROP_FLOAT, PROP_ANGLE); + RNA_def_property_ui_text(prop, "Sun Rotation", "Rotation of sun around zenith"); + RNA_def_property_range(prop, 0.0f, 2.0f * M_PI); + RNA_def_property_float_default(prop, 0.0f); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "altitude", PROP_INT, PROP_NONE); + RNA_def_property_ui_text(prop, "Altitude", "Altitude height from sea level in meters"); + RNA_def_property_range(prop, 0, 60000); + RNA_def_property_int_default(prop, 0); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "air_density", PROP_FLOAT, PROP_FACTOR); + RNA_def_property_ui_text(prop, "Air", "Density of air molecules (Rayleigh scattering)"); + RNA_def_property_range(prop, 0.0f, 10.0f); + RNA_def_property_float_default(prop, 1.0f); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "dust_density", PROP_FLOAT, PROP_FACTOR); + RNA_def_property_ui_text(prop, "Dust", "Density of dust and water molecules (Mie scattering)"); + RNA_def_property_range(prop, 0.0f, 10.0f); + RNA_def_property_float_default(prop, 1.0f); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); + + prop = RNA_def_property(srna, "ozone_density", PROP_FLOAT, PROP_FACTOR); + RNA_def_property_ui_text(prop, "Ozone", "Density of Ozone layer (Ozone absorption)"); + RNA_def_property_range(prop, 0.0f, 10.0f); + RNA_def_property_float_default(prop, 1.0f); + RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update"); } static const EnumPropertyItem sh_tex_prop_interpolation_items[] = { diff --git a/source/blender/nodes/shader/nodes/node_shader_tex_sky.c b/source/blender/nodes/shader/nodes/node_shader_tex_sky.c index cb91e086cf4..0daa948c139 100644 --- a/source/blender/nodes/shader/nodes/node_shader_tex_sky.c +++ b/source/blender/nodes/shader/nodes/node_shader_tex_sky.c @@ -41,8 +41,15 @@ static void node_shader_init_tex_sky(bNodeTree *UNUSED(ntree), bNode *node) tex->sun_direction[2] = 1.0f; tex->turbidity = 2.2f; tex->ground_albedo = 0.3f; - tex->sky_model = SHD_SKY_NEW; - + tex->sun_disc = true; + tex->sun_size = DEG2RADF(0.545); + tex->sun_elevation = M_PI_2; + tex->sun_rotation = 0.0f; + tex->altitude = 0; + tex->air_density = 1.0f; + tex->dust_density = 1.0f; + tex->ozone_density = 1.0f; + tex->sky_model = SHD_SKY_NISHITA; node->storage = tex; } |