diff --git a/intern/cycles/kernel/integrator/shade_volume.h b/intern/cycles/kernel/integrator/shade_volume.h index c840324c1eb..69962549030 100644 --- a/intern/cycles/kernel/integrator/shade_volume.h +++ b/intern/cycles/kernel/integrator/shade_volume.h @@ -278,21 +278,25 @@ ccl_device float volume_equiangular_sample(const ccl_private Ray *ccl_restrict r ccl_private float *pdf) { const float delta = dot((coeffs.P - ray->P), ray->D); - const float D = safe_sqrtf(len_squared(coeffs.P - ray->P) - delta * delta); + const float D = len(coeffs.P - ray->P - ray->D * delta); if (UNLIKELY(D == 0.0f)) { *pdf = 0.0f; return 0.0f; } const float tmin = coeffs.t_range.min; const float tmax = coeffs.t_range.max; + const float theta_a = atan2f(tmin - delta, D); const float theta_b = atan2f(tmax - delta, D); - const float t_ = D * tanf((xi * theta_b) + (1 - xi) * theta_a); - if (UNLIKELY(theta_b == theta_a)) { - *pdf = 0.0f; - return 0.0f; + const float theta_d = theta_b - theta_a; + if (UNLIKELY(theta_d < 1e-6f)) { + /* Use uniform sampling when `theta_d` is too small. */ + *pdf = safe_divide(1.0f, tmax - tmin); + return mix(tmin, tmax, xi); } - *pdf = D / ((theta_b - theta_a) * (D * D + t_ * t_)); + + const float t_ = D * tanf((xi * theta_b) + (1 - xi) * theta_a); + *pdf = D / (theta_d * (D * D + t_ * t_)); return clamp(delta + t_, tmin, tmax); /* clamp is only for float precision errors */ } @@ -302,22 +306,23 @@ ccl_device float volume_equiangular_pdf(const ccl_private Ray *ccl_restrict ray, const float sample_t) { const float delta = dot((coeffs.P - ray->P), ray->D); - const float D = safe_sqrtf(len_squared(coeffs.P - ray->P) - delta * delta); + const float D = len(coeffs.P - ray->P - ray->D * delta); if (UNLIKELY(D == 0.0f)) { return 0.0f; } const float tmin = coeffs.t_range.min; const float tmax = coeffs.t_range.max; - const float t_ = sample_t - delta; const float theta_a = atan2f(tmin - delta, D); const float theta_b = atan2f(tmax - delta, D); - if (UNLIKELY(theta_b == theta_a)) { - return 0.0f; + const float theta_d = theta_b - theta_a; + if (UNLIKELY(theta_d < 1e-6f)) { + return safe_divide(1.0f, tmax - tmin); } - const float pdf = D / ((theta_b - theta_a) * (D * D + t_ * t_)); + const float t_ = sample_t - delta; + const float pdf = D / (theta_d * (D * D + t_ * t_)); return pdf; } @@ -342,7 +347,7 @@ ccl_device_inline bool volume_equiangular_valid_ray_segment(KernelGlobals kg, /* Point light, the whole range of the ray is visible. */ kernel_assert(ls->type == LIGHT_POINT); - return true; + return !t_range->is_empty(); } /* Emission */ diff --git a/intern/cycles/kernel/light/tree.h b/intern/cycles/kernel/light/tree.h index a230ca5e506..f95623cab4f 100644 --- a/intern/cycles/kernel/light/tree.h +++ b/intern/cycles/kernel/light/tree.h @@ -341,9 +341,19 @@ ccl_device void light_tree_node_importance(const float3 P, const float3 closest_point = P + D * clamp(closest_t, 0.0f, t); /* Minimal distance of the ray to the cluster. */ distance = len(centroid - P - D * closest_t); + /* Estimate `theta_b - theta_a` using the centroid of the cluster and the complete ray * segment in volume. */ - theta_d = fast_atan2f(t - closest_t, distance) + fast_atan2f(closest_t, distance); + if (t == FLT_MAX) { + theta_d = fast_atan2f(closest_t, distance) + M_PI_2_F; + } + else { + /* Original equation is `theta_d = atan((t - closest) /d) + atan(closest / d)`, convert to + * the below equation using the equality `atan(a) + atan(b) = atan2(a + b, 1 - a*b)` for + * better precision at small angles. */ + theta_d = atan2f(t, distance - closest_t * safe_divide(t - closest_t, distance)); + } + /* Vector that forms a minimal angle with the emitter centroid. */ point_to_centroid = -compute_v(centroid, P, D, bcone.axis, t); cos_theta_u = light_tree_cos_bound_subtended_angle(bbox, centroid, closest_point); @@ -429,7 +439,12 @@ ccl_device void light_tree_emitter_importance(KernelGlobals kg, const float closest_t = dot(centroid - P, D); P_c += D * clamp(closest_t, 0.0f, t); const float d = len(centroid - P - D * closest_t); - theta_d = fast_atan2f(t - closest_t, d) + fast_atan2f(closest_t, d); + if (t == FLT_MAX) { + theta_d = fast_atan2f(closest_t, d) + M_PI_2_F; + } + else { + theta_d = atan2f(t, d - closest_t * safe_divide(t - closest_t, d)); + } } /* Early out if the emitter is guaranteed to be invisible. */