Fix #140892: Cycles: equiangular sampling numerical issue

When `delta` is significantly larger than the difference between `tmin`
and `tmax`, the precision of `theta_a` and `theta_b` reduces, resulting
in banding artefacts.

For equiangular sampling, we use uniform sampling to fix this case.

For light tree, we use the equality
atan(a) + atan(b) = atan2(a + b, 1 - a*b).

Pull Request: https://projects.blender.org/blender/blender/pulls/142845
This commit is contained in:
Weizhen Huang
2025-07-23 12:56:11 +02:00
committed by Weizhen Huang
parent 5f230eec8e
commit bb689687a7
2 changed files with 34 additions and 14 deletions

View File

@@ -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 */

View File

@@ -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. */