tmaq
/
tornavis
1
0
Fork 0
Code Issues 2 Pull Requests 29 Packages Projects 24 Releases Wiki Activity

Cycles: improve equiangular sampling in volume 

By restricting the sample range along the ray to the valid segment.

Supports

**Mesh Light**
- [x] restrict the ray segment to the side with MIS

**Area Light**
- [x] when the spread is zero, find the intersection of the ray and the bounding box/cylinder of the rectangle/ellipse area light beam
- [x] when the spread is non-zero, find the intersection of the ray and the minimal enclosing cone of the area light beam
*note the result is also unbiased when we just consider the cone from the sampled point in volume segment. Far away from the light source it's less noisy than the current solution, but near the light source it's much noisier. We have to restrict the sample region on the area light to the part that lits the ray then, I haven't tried yet to see if it would be less noisy.*

**Point Light**
- [x] the complete ray segment should be valid.

**Spot Light**
- [x] intersect the ray with the spot light cone
- [x] support non-zero radius

Pull Request: https://projects.blender.org/blender/blender/pulls/119438
This commit is contained in:
Weizhen Huang 2024-03-25 13:02:02 +01:00 committed by Weizhen Huang
parent ab93a426a0
commit 082b68fcb9
10 changed files with 343 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
intern/cycles/kernel/integrator/shade_volume.h
View File

@ -64,6 +64,11 @@ typedef struct VolumeShaderCoefficients {
Spectrum emission;
} VolumeShaderCoefficients;
typedef struct EquiangularCoefficients {
float3 P;
float2 t_range;
} EquiangularCoefficients;
/* Evaluate shader to get extinction coefficient at P. */
ccl_device_inline bool shadow_volume_shader_sample(KernelGlobals kg,
IntegratorShadowState state,
@ -264,18 +269,18 @@ ccl_device void volume_shadow_heterogeneous(KernelGlobals kg,
# define VOLUME_SAMPLE_PDF_CUTOFF 1e-8f
ccl_device float volume_equiangular_sample(ccl_private const Ray *ccl_restrict ray,
const float3 light_P,
ccl_private const EquiangularCoefficients &coeffs,
const float xi,
ccl_private float *pdf)
{
const float tmin = ray->tmin;
const float tmax = ray->tmax;
const float delta = dot((light_P - ray->P), ray->D);
const float D = safe_sqrtf(len_squared(light_P - ray->P) - delta * delta);
const float delta = dot((coeffs.P - ray->P), ray->D);
const float D = safe_sqrtf(len_squared(coeffs.P - ray->P) - delta * delta);
if (UNLIKELY(D == 0.0f)) {
*pdf = 0.0f;
return 0.0f;
}
const float tmin = coeffs.t_range.x;
const float tmax = coeffs.t_range.y;
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);
@ -289,17 +294,17 @@ ccl_device float volume_equiangular_sample(ccl_private const Ray *ccl_restrict r
}
ccl_device float volume_equiangular_pdf(ccl_private const Ray *ccl_restrict ray,
const float3 light_P,
ccl_private const EquiangularCoefficients &coeffs,
const float sample_t)
{
const float delta = dot((light_P - ray->P), ray->D);
const float D = safe_sqrtf(len_squared(light_P - ray->P) - delta * delta);
const float delta = dot((coeffs.P - ray->P), ray->D);
const float D = safe_sqrtf(len_squared(coeffs.P - ray->P) - delta * delta);
if (UNLIKELY(D == 0.0f)) {
return 0.0f;
}
const float tmin = ray->tmin;
const float tmax = ray->tmax;
const float tmin = coeffs.t_range.x;
const float tmax = coeffs.t_range.y;
const float t_ = sample_t - delta;
const float theta_a = atan2f(tmin - delta, D);
@ -313,6 +318,36 @@ ccl_device float volume_equiangular_pdf(ccl_private const Ray *ccl_restrict ray,
return pdf;
}
ccl_device_inline bool volume_equiangular_valid_ray_segment(KernelGlobals kg,
const float3 ray_P,
const float3 ray_D,
ccl_private float2 *t_range,
const ccl_private LightSample *ls)
{
# ifdef __LIGHT_TREE__
/* Do not compute ray segment until #119389 is landed. */
if (kernel_data.integrator.use_light_tree) {
return true;
}
# endif
if (ls->type == LIGHT_SPOT) {
ccl_global const KernelLight *klight = &kernel_data_fetch(lights, ls->lamp);
return spot_light_valid_ray_segment(klight, ray_P, ray_D, t_range);
}
if (ls->type == LIGHT_AREA) {
ccl_global const KernelLight *klight = &kernel_data_fetch(lights, ls->lamp);
return area_light_valid_ray_segment(&klight->area, ray_P - klight->co, ray_D, t_range);
}
if (ls->type == LIGHT_TRIANGLE) {
return triangle_light_valid_ray_segment(kg, ray_P - ls->P, ray_D, t_range, ls);
}
/* Point light, the whole range of the ray is visible. */
kernel_assert(ls->type == LIGHT_POINT);
return true;
}
/* Distance sampling */
ccl_device float volume_distance_sample(float max_t,
@ -403,7 +438,7 @@ typedef struct VolumeIntegrateState {
ccl_device_forceinline void volume_integrate_step_scattering(
ccl_private const ShaderData *sd,
ccl_private const Ray *ray,
const float3 equiangular_light_P,
ccl_private const EquiangularCoefficients &equiangular_coeffs,
ccl_private const VolumeShaderCoefficients &ccl_restrict coeff,
const Spectrum transmittance,
ccl_private VolumeIntegrateState &ccl_restrict vstate,
@ -474,7 +509,7 @@ ccl_device_forceinline void volume_integrate_step_scattering(
/* Multiple importance sampling. */
if (vstate.use_mis) {
const float equiangular_pdf = volume_equiangular_pdf(ray, equiangular_light_P, new_t);
const float equiangular_pdf = volume_equiangular_pdf(ray, equiangular_coeffs, new_t);
const float mis_weight = power_heuristic(vstate.distance_pdf * distance_pdf,
equiangular_pdf);
result.direct_throughput *= 2.0f * mis_weight;
@ -509,7 +544,7 @@ ccl_device_forceinline void volume_integrate_heterogeneous(
ccl_global float *ccl_restrict render_buffer,
const float object_step_size,
const VolumeSampleMethod direct_sample_method,
const float3 equiangular_light_P,
ccl_private const EquiangularCoefficients &equiangular_coeffs,
ccl_private VolumeIntegrateResult &result)
{
PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_INTEGRATE);
@ -560,7 +595,7 @@ ccl_device_forceinline void volume_integrate_heterogeneous(
/* Equiangular sampling: compute distance and PDF in advance. */
if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
result.direct_t = volume_equiangular_sample(
ray, equiangular_light_P, vstate.rscatter, &vstate.equiangular_pdf);
ray, equiangular_coeffs, vstate.rscatter, &vstate.equiangular_pdf);
}
# ifdef __PATH_GUIDING__
result.direct_sample_method = vstate.direct_sample_method;
@ -614,7 +649,7 @@ ccl_device_forceinline void volume_integrate_heterogeneous(
/* Scattering and absorption. */
volume_integrate_step_scattering(
sd, ray, equiangular_light_P, coeff, transmittance, vstate, result);
sd, ray, equiangular_coeffs, coeff, transmittance, vstate, result);
}
else {
/* Absorption only. */
@ -673,7 +708,7 @@ ccl_device_forceinline bool integrate_volume_equiangular_sample_light(
ccl_private const Ray *ccl_restrict ray,
ccl_private const ShaderData *ccl_restrict sd,
ccl_private const RNGState *ccl_restrict rng_state,
ccl_private float3 *ccl_restrict P)
ccl_private EquiangularCoefficients *ccl_restrict equiangular_coeffs)
{
/* Test if there is a light or BSDF that needs direct light. */
if (!kernel_data.integrator.use_direct_light) {
@ -708,9 +743,10 @@ ccl_device_forceinline bool integrate_volume_equiangular_sample_light(
return false;
}
*P = ls.P;
equiangular_coeffs->P = ls.P;
return true;
return volume_equiangular_valid_ray_segment(
kg, ray->P, ray->D, &equiangular_coeffs->t_range, &ls);
}
/* Path tracing: sample point on light and evaluate light shader, then
@ -990,10 +1026,12 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
/* Sample light ahead of volume stepping, for equiangular sampling. */
/* TODO: distant lights are ignored now, but could instead use even distribution. */
const bool need_light_sample = !(INTEGRATOR_STATE(state, path, flag) & PATH_RAY_TERMINATE);
float3 equiangular_P = zero_float3();
EquiangularCoefficients equiangular_coeffs = {zero_float3(), make_float2(ray->tmin, ray->tmax)};
const bool have_equiangular_sample = need_light_sample &&
integrate_volume_equiangular_sample_light(
kg, state, ray, &sd, &rng_state, &equiangular_P);
kg, state, ray, &sd, &rng_state, &equiangular_coeffs);
VolumeSampleMethod direct_sample_method = (have_equiangular_sample) ?
volume_stack_sample_method(kg, state) :
@ -1023,7 +1061,7 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
render_buffer,
step_size,
direct_sample_method,
equiangular_P,
equiangular_coeffs,
result);
/* Perform path termination. The intersect_closest will have already marked this path

58
intern/cycles/kernel/light/area.h
View File

@ -233,6 +233,11 @@ ccl_device bool area_light_spread_clamp_light(const float3 P,
return true;
}
ccl_device_forceinline bool area_light_is_ellipse(const ccl_global KernelAreaLight *light)
{
return light->invarea < 0.0f;
}
/* Common API. */
/* Compute `eval_fac` and `pdf`. Also sample a new position on the light if `sample_coord`. */
template<bool in_volume_segment>
@ -338,7 +343,7 @@ ccl_device_inline bool area_light_sample(const ccl_global KernelLight *klight,
const float light_v = dot(inplane, klight->area.axis_v) / klight->area.len_v;
if (!in_volume_segment) {
const bool is_ellipse = (klight->area.invarea < 0.0f);
const bool is_ellipse = area_light_is_ellipse(&klight->area);
/* Sampled point lies outside of the area light. */
if (is_ellipse && (sqr(light_u) + sqr(light_v) > 0.25f)) {
@ -380,7 +385,7 @@ ccl_device_inline bool area_light_intersect(const ccl_global KernelLight *klight
{
/* Area light. */
const float invarea = fabsf(klight->area.invarea);
const bool is_ellipse = (klight->area.invarea < 0.0f);
const bool is_ellipse = area_light_is_ellipse(&klight->area);
if (invarea == 0.0f) {
return false;
}
@ -428,6 +433,55 @@ ccl_device_inline bool area_light_sample_from_intersection(
return area_light_eval<false>(klight, ray_P, &light_P, ls, zero_float2(), false);
}
/* Returns the maximal distance between the light center and the boundary. */
ccl_device_forceinline float area_light_max_extent(const ccl_global KernelAreaLight *light)
{
return 0.5f * (area_light_is_ellipse(light) ? fmaxf(light->len_u, light->len_v) :
len(make_float2(light->len_u, light->len_v)));
}
/* Find the ray segment lit by the area light. */
ccl_device_inline bool area_light_valid_ray_segment(const ccl_global KernelAreaLight *light,
float3 P,
float3 D,
ccl_private float2 *t_range)
{
bool valid;
const float tan_half_spread = light->tan_half_spread;
float3 axis = light->dir;
const bool angle_almost_zero = (tan_half_spread < 1e-5f);
if (angle_almost_zero) {
/* Map to local coordinate of the light. Do not use `itfm` in `KernelLight` as there might be
* additional scaling in the light size. */
const Transform tfm = make_transform(light->axis_u, light->axis_v, axis);
P = transform_point(&tfm, P);
D = transform_direction(&tfm, D);
axis = make_float3(0.0f, 0.0f, 1.0f);
const float half_len_u = 0.5f * light->len_u;
const float half_len_v = 0.5f * light->len_v;
if (area_light_is_ellipse(light)) {
valid = ray_infinite_cylinder_intersect(P, D, half_len_u, half_len_v, t_range);
}
else {
const float3 bbox_min = make_float3(-half_len_u, -half_len_v, 0.0f);
const float3 bbox_max = make_float3(half_len_u, half_len_v, FLT_MAX);
valid = ray_aabb_intersect(bbox_min, bbox_max, P, D, t_range);
}
}
else {
/* Conservative estimation with the smallest possible cone covering the whole spread. */
const float3 apex_to_point = P + area_light_max_extent(light) / tan_half_spread * axis;
const float cos_angle_sq = 1.0f / (1.0f + sqr(tan_half_spread));
valid = ray_cone_intersect(axis, apex_to_point, D, cos_angle_sq, t_range);
}
/* Limit the range to the positive side of the area light. */
return valid && ray_plane_intersect(axis, P, D, t_range);
}
template<bool in_volume_segment>
ccl_device_forceinline bool area_light_tree_parameters(const ccl_global KernelLight *klight,
const float3 centroid,

18
intern/cycles/kernel/light/spot.h
View File

@ -265,6 +265,24 @@ ccl_device_inline bool spot_light_sample_from_intersection(
return true;
}
/* Find the ray segment lit by the spot light. */
ccl_device_inline bool spot_light_valid_ray_segment(const ccl_global KernelLight *klight,
const float3 P,
const float3 D,
ccl_private float2 *t_range)
{
/* Convert to local space of the spot light. */
const Transform itfm = klight->itfm;
float3 local_P = P + klight->spot.dir * klight->spot.ray_segment_dp;
local_P = transform_point(&itfm, local_P);
const float3 local_D = transform_direction(&itfm, D);
const float3 axis = make_float3(0.0f, 0.0f, -1.0f);
/* Intersect the ray with the smallest enclosing cone of the light spread. */
return ray_cone_intersect(
axis, local_P, local_D, sqr(klight->spot.cos_half_spot_angle), t_range);
}
template<bool in_volume_segment>
ccl_device_forceinline bool spot_light_tree_parameters(const ccl_global KernelLight *klight,
const float3 centroid,

19
intern/cycles/kernel/light/triangle.h
View File

@ -269,6 +269,25 @@ ccl_device_forceinline bool triangle_light_sample(KernelGlobals kg,
return (ls->pdf > 0.0f);
}
/* Find the ray segment lit by the triangle light. */
ccl_device_inline bool triangle_light_valid_ray_segment(KernelGlobals kg,
const float3 P,
const float3 D,
ccl_private float2 *t_range,
const ccl_private LightSample *ls)
{
const int shader_flag = kernel_data_fetch(shaders, ls->shader & SHADER_MASK).flags;
const int SD_MIS_BOTH = SD_MIS_BACK | SD_MIS_FRONT;
if ((shader_flag & SD_MIS_BOTH) == SD_MIS_BOTH) {
/* Both sides are sampled, the complete ray segment is visible. */
return true;
}
/* Only one side is sampled, intersect the ray and the triangle light plane to find the visible
* ray segment. Flip normal if Emission Sampling is set to back. */
return ray_plane_intersect((shader_flag & SD_MIS_BACK) ? -ls->Ng : ls->Ng, P, D, t_range);
}
template<bool in_volume_segment>
ccl_device_forceinline bool triangle_light_tree_parameters(
KernelGlobals kg,

2
intern/cycles/kernel/types.h
View File

@ -1376,6 +1376,8 @@ typedef struct KernelSpotLight {
int is_sphere;
/* For non-uniform object scaling, the actual spread might be different. */
float cos_half_larger_spread;
/* Distance from the apex of the smallest enclosing cone of the light spread to light center. */
float ray_segment_dp;
} KernelSpotLight;
/* PointLight is SpotLight with only radius and invarea being used. */

3
intern/cycles/scene/light.cpp
View File

@ -1362,6 +1362,9 @@ void LightManager::device_update_lights(Device *device, DeviceScene *dscene, Sce
/* Choose the angle which spans a larger cone. */
klights[light_index].spot.cos_half_larger_spread = inversesqrtf(
1.0f + tan_sq * fmaxf(len_u_sq, len_v_sq) / len_w_sq);
/* radius / sin(half_angle_small) */
klights[light_index].spot.ray_segment_dp =
light->size * sqrtf(1.0f + len_w_sq / (tan_sq * fminf(len_u_sq, len_v_sq)));
}
klights[light_index].shader_id = shader_id;

40
intern/cycles/util/math.h
View File

@ -1030,6 +1030,46 @@ ccl_device_inline uint32_t reverse_integer_bits(uint32_t x)
#endif
}
/* Check if intervals (first->x, first->y) and (second.x, second.y) intersect, and replace the
* first interval with their intersection. */
ccl_device_inline bool intervals_intersect(ccl_private float2 *first, const float2 second)
{
first->x = fmaxf(first->x, second.x);
first->y = fminf(first->y, second.y);
return first->x < first->y;
}
/* Solve quadratic equation a*x^2 + b*x + c = 0, adapted from Mitsuba 3
* The solution is ordered so that x1 <= x2.
* Returns true if at least one solution is found. */
ccl_device_inline bool solve_quadratic(
const float a, const float b, const float c, ccl_private float &x1, ccl_private float &x2)
{
/* If the equation is linear, the solution is -c/b, but b has to be non-zero. */
const bool valid_linear = (a == 0.0f) && (b != 0.0f);
x1 = x2 = -c / b;
const float discriminant = sqr(b) - 4.0f * a * c;
/* Allow slightly negative discriminant in case of numerical precision issues. */
const bool valid_quadratic = (a != 0.0f) && (discriminant > -1e-5f);
if (valid_quadratic) {
/* Numerically stable version of (-b ± sqrt(discriminant)) / (2 * a), avoiding catastrophic
* cancellation when `b` is very close to `sqrt(discriminant)`, by finding the solution of
* greater magnitude which does not suffer from loss of precision, then using the identity
* x1 * x2 = c / a. */
const float temp = -0.5f * (b + copysignf(safe_sqrtf(discriminant), b));
const float r1 = temp / a;
const float r2 = c / temp;
x1 = fminf(r1, r2);
x2 = fmaxf(r1, r2);
}
return (valid_linear || valid_quadratic);
}
CCL_NAMESPACE_END
#endif /* __UTIL_MATH_H__ */

134
intern/cycles/util/math_intersect.h
View File

@ -302,6 +302,140 @@ ccl_device bool ray_quad_intersect(float3 ray_P,
return true;
}
/* Find the ray segment that lies in the same side as the normal `N` of the plane.
* `P` is the vector pointing from any point on the plane to the ray origin. */
ccl_device bool ray_plane_intersect(const float3 N,
const float3 P,
const float3 ray_D,
ccl_private float2 *t_range)
{
const float DN = dot(ray_D, N);
/* Distance from P to the plane. */
const float t = -dot(P, N) / DN;
/* Limit the range to the positive side. */
if (DN > 0.0f) {
t_range->x = fmaxf(t_range->x, t);
}
else {
t_range->y = fminf(t_range->y, t);
}
return t_range->x < t_range->y;
}
/* Find the ray segment inside an axis-aligned bounding box. */
ccl_device bool ray_aabb_intersect(const float3 bbox_min,
const float3 bbox_max,
const float3 ray_P,
const float3 ray_D,
ccl_private float2 *t_range)
{
const float3 inv_ray_D = rcp(ray_D);
/* Absolute distances to lower and upper box coordinates; */
const float3 t_lower = (bbox_min - ray_P) * inv_ray_D;
const float3 t_upper = (bbox_max - ray_P) * inv_ray_D;
/* The four t-intervals (for x-/y-/z-slabs, and ray p(t)). */
const float4 tmins = float3_to_float4(min(t_lower, t_upper), t_range->x);
const float4 tmaxes = float3_to_float4(max(t_lower, t_upper), t_range->y);
/* Max of mins and min of maxes. */
const float tmin = reduce_max(tmins);
const float tmax = reduce_min(tmaxes);
*t_range = make_float2(tmin, tmax);
return tmin < tmax;
}
/* Find the segment of a ray defined by P + D * t that lies inside a cylinder defined by
* (x / len_u)^2 + (y / len_v)^2 = 1. */
ccl_device_inline bool ray_infinite_cylinder_intersect(const float3 P,
const float3 D,
const float len_u,
const float len_v,
ccl_private float2 *t_range)
{
/* Convert to a 2D problem. */
const float2 inv_len = 1.0f / make_float2(len_u, len_v);
float2 P_proj = float3_to_float2(P) * inv_len;
const float2 D_proj = float3_to_float2(D) * inv_len;
/* Solve quadratic equation a*t^2 + 2b*t + c = 0. */
const float a = dot(D_proj, D_proj);
float b = dot(P_proj, D_proj);
/* Move ray origin closer to the cylinder to prevent precision issue when the ray is far away. */
const float t_mid = -b / a;
P_proj += D_proj * t_mid;
/* Recompute b from the shifted origin. */
b = dot(P_proj, D_proj);
const float c = dot(P_proj, P_proj) - 1.0f;
float tmin, tmax;
const bool valid = solve_quadratic(a, 2.0f * b, c, tmin, tmax);
return valid && intervals_intersect(t_range, make_float2(tmin, tmax) + t_mid);
}
/* *
* Find the ray segment inside a single-sided cone.
*
* \param axis: a unit-length direction around which the cone has a circular symmetry
* \param P: the vector pointing from the cone apex to the ray origin
* \param D: the direction of the ray, does not need to have unit-length
* \param cos_angle_sq: `sqr(cos(half_aperture_of_the_cone))`
* \param t_range: the lower and upper bounds between which the ray lies inside the cone
* \return whether the intersection exists and is in the provided range
*
* See https://www.geometrictools.com/Documentation/IntersectionLineCone.pdf for illustration
*/
ccl_device_inline bool ray_cone_intersect(const float3 axis,
const float3 P,
float3 D,
const float cos_angle_sq,
ccl_private float2 *t_range)
{
if (cos_angle_sq < 1e-4f) {
/* The cone is nearly a plane. */
return ray_plane_intersect(axis, P, D, t_range);
}
const float inv_len = inversesqrtf(len_squared(D));
D *= inv_len;
const float AD = dot(axis, D);
const float AP = dot(axis, P);
const float a = sqr(AD) - cos_angle_sq;
const float b = 2.0f * (AD * AP - cos_angle_sq * dot(D, P));
const float c = sqr(AP) - cos_angle_sq * dot(P, P);
float tmin = 0.0f, tmax = FLT_MAX;
bool valid = solve_quadratic(a, b, c, tmin, tmax);
/* Check if the intersections are in the same hemisphere as the cone. */
const bool tmin_valid = AP + tmin * AD > 0.0f;
const bool tmax_valid = AP + tmax * AD > 0.0f;
valid &= (tmin_valid || tmax_valid);
if (!tmax_valid) {
tmax = tmin;
tmin = 0.0f;
}
else if (!tmin_valid) {
tmin = tmax;
tmax = FLT_MAX;
}
return valid && intervals_intersect(t_range, make_float2(tmin, tmax) * inv_len);
}
CCL_NAMESPACE_END
#endif /* __UTIL_MATH_INTERSECT_H__ */

11
intern/cycles/util/transform.h
View File

@ -161,6 +161,17 @@ ccl_device_inline Transform make_transform(float a,
return t;
}
ccl_device_inline Transform make_transform(const float3 x, const float3 y, const float3 z)
{
Transform t;
t.x = float3_to_float4(x, 0.0f);
t.y = float3_to_float4(y, 0.0f);
t.z = float3_to_float4(z, 0.0f);
return t;
}
ccl_device_inline Transform euler_to_transform(const float3 euler)
{
float cx = cosf(euler.x);

2
tests/data

@ -1 +1 @@
Subproject commit 5038ad7165fd1a77e61e0d2d6efdadd6ea7c0dfb
Subproject commit 00af9c65712b6aa78ce6eb0c62c5aafb7a867f18