884 lines
41 KiB
GLSL
884 lines
41 KiB
GLSL
// Copyright (C) 2021, Christoph Peters, Karlsruhe Institute of Technology
|
|
//
|
|
// This source code file is licensed under both the three-clause BSD license
|
|
// and the GPLv3. You may select, at your option, one of the two licenses. The
|
|
// corresponding license headers follow:
|
|
//
|
|
//
|
|
// Three-clause BSD license:
|
|
//
|
|
// Redistribution and use in source and binary forms, with or without
|
|
// modification, are permitted provided that the following conditions are met:
|
|
//
|
|
// 1. Redistributions of source code must retain the above copyright notice,
|
|
// this list of conditions and the following disclaimer.
|
|
//
|
|
// 2. Redistributions in binary form must reproduce the above copyright
|
|
// notice, this list of conditions and the following disclaimer in the
|
|
// documentation and/or other materials provided with the distribution.
|
|
//
|
|
// 3. Neither the name of the copyright holder nor the names of its
|
|
// contributors may be used to endorse or promote products derived from
|
|
// this software without specific prior written permission.
|
|
//
|
|
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
|
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
// POSSIBILITY OF SUCH DAMAGE.
|
|
//
|
|
//
|
|
// GPLv3:
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License, or
|
|
// (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
|
|
|
|
#include "math_constants.glsl"
|
|
|
|
|
|
/*! This structure carries intermediate results that only need to be computed
|
|
once per polygon and shading point to take samples proportional to solid
|
|
angle. Sampling is performed by subdividing the convex polygon into
|
|
triangles as triangle fan around vertex 0 and then using our variant of
|
|
Arvo's method.*/
|
|
struct solid_angle_polygon_t {
|
|
//! The number of vertices that form the polygon
|
|
uint vertex_count;
|
|
//! Normalized direction vectors from the shading point to each vertex
|
|
vec3 vertex_dirs[MAX_POLYGON_VERTEX_COUNT];
|
|
/*! A few intermediate quantities about the triangle consisting of vertices
|
|
i + 1, 0, and i + 2. If the three vertices are v0, v1, v2, the entries
|
|
are determinant(mat3(v0, v1, v2)), dot(v0 + v1, v2) and
|
|
1.0f + dot(v0, v1).*/
|
|
vec3 triangle_parameters[MAX_POLYGON_VERTEX_COUNT - 2];
|
|
//! At index i, this array holds the solid angle of the triangle fan formed
|
|
//! by vertices 0 to i + 2.
|
|
float fan_solid_angles[MAX_POLYGON_VERTEX_COUNT - 2];
|
|
//! The total solid angle of the polygon
|
|
float solid_angle;
|
|
};
|
|
|
|
|
|
/*! A piecewise polynomial approximation to positive_atan(y). The maximal
|
|
absolute error is 1.16e-05f. At least on Turing GPUs, it is faster but also
|
|
significantly less accurate. The proper atan has at most 2 ulps of error
|
|
there.*/
|
|
float fast_positive_atan(float y) {
|
|
float rx;
|
|
float ry;
|
|
float rz;
|
|
rx = (abs(y) > 1.0f) ? (1.0f / abs(y)) : abs(y);
|
|
ry = rx * rx;
|
|
rz = fma(ry, 0.02083509974181652f, -0.08513300120830536);
|
|
rz = fma(ry, rz, 0.18014100193977356f);
|
|
rz = fma(ry, rz, -0.3302994966506958f);
|
|
ry = fma(ry, rz, 0.9998660087585449f);
|
|
rz = fma(-2.0f * ry, rx, M_HALF_PI);
|
|
rz = (abs(y) > 1.0f) ? rz : 0.0f;
|
|
rx = fma(rx, ry, rz);
|
|
return (y < 0.0f) ? (M_PI - rx) : rx;
|
|
}
|
|
|
|
|
|
/*! Returns an angle between 0 and M_PI such that tan(angle) == tangent. In
|
|
other words, it is a version of atan() that is offset to be non-negative.
|
|
Note that it may be switched to an approximate mode by the
|
|
USE_BIASED_PROJECTED_SOLID_ANGLE_SAMPLING flag.*/
|
|
float positive_atan(float tangent) {
|
|
#ifdef USE_BIASED_PROJECTED_SOLID_ANGLE_SAMPLING
|
|
return fast_positive_atan(tangent);
|
|
#else
|
|
float offset = (tangent < 0.0f) ? M_PI : 0.0f;
|
|
return atan(tangent) + offset;
|
|
#endif
|
|
}
|
|
|
|
|
|
/*! Prepares all intermediate values to sample a triangle fan around vertex 0
|
|
(e.g. a convex polygon) proportional to solid angle using our method.
|
|
\param vertex_count Number of vertices forming the polygon.
|
|
\param vertices List of vertex locations.
|
|
\param shading_position The location of the shading point.
|
|
\return Input for sample_solid_angle_polygon().*/
|
|
solid_angle_polygon_t prepare_solid_angle_polygon_sampling(uint vertex_count, vec3 vertices[MAX_POLYGON_VERTEX_COUNT], vec3 shading_position) {
|
|
solid_angle_polygon_t polygon;
|
|
polygon.vertex_count = vertex_count;
|
|
// Normalize vertex directions
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT; ++i) {
|
|
polygon.vertex_dirs[i] = normalize(vertices[i] - shading_position);
|
|
}
|
|
// Prepare a Householder transform that maps vertex 0 onto (+/-1, 0, 0). We
|
|
// only store the yz-components of that Householder vector and a factor of
|
|
// 2.0f / sqrt(abs(polygon.vertex_dirs[0].x) + 1.0f) is pulled in there to
|
|
// save on multiplications later. This approach is necessary to avoid
|
|
// numerical instabilities in determinant computation below.
|
|
float householder_sign = (polygon.vertex_dirs[0].x > 0.0f) ? -1.0f : 1.0f;
|
|
vec2 householder_yz = polygon.vertex_dirs[0].yz * (1.0f / (abs(polygon.vertex_dirs[0].x) + 1.0f));
|
|
// Compute solid angles and prepare sampling
|
|
polygon.solid_angle = 0.0f;
|
|
float previous_dot_1_2 = dot(polygon.vertex_dirs[0], polygon.vertex_dirs[1]);
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT - 2; ++i) {
|
|
if (i >= 1 && i + 2 >= vertex_count) break;
|
|
// We look at one triangle of the triangle fan at a time
|
|
vec3 vertices[3] = {
|
|
polygon.vertex_dirs[i + 1],
|
|
polygon.vertex_dirs[0],
|
|
polygon.vertex_dirs[i + 2]};
|
|
float dot_0_1 = previous_dot_1_2;
|
|
float dot_0_2 = dot(vertices[0], vertices[2]);
|
|
float dot_1_2 = dot(vertices[1], vertices[2]);
|
|
previous_dot_1_2 = dot_1_2;
|
|
// Compute the bottom right minor of vertices after application of the
|
|
// Householder transform
|
|
float dot_householder_0 = fma(-householder_sign, vertices[0].x, dot_0_1);
|
|
float dot_householder_2 = fma(-householder_sign, vertices[2].x, dot_1_2);
|
|
mat2 bottom_right_minor = mat2(
|
|
fma(vec2(-dot_householder_0), householder_yz, vertices[0].yz),
|
|
fma(vec2(-dot_householder_2), householder_yz, vertices[2].yz));
|
|
// The absolute value of the determinant of vertices equals the 2x2
|
|
// determinant because the Householder transform turns the first column
|
|
// into (+/-1, 0, 0)
|
|
float simplex_volume = abs(determinant(bottom_right_minor));
|
|
// Compute the solid angle of the triangle using a formula proposed by:
|
|
// A. Van Oosterom and J. Strackee, 1983, The Solid Angle of a
|
|
// Plane Triangle, IEEE Transactions on Biomedical Engineering 30:2
|
|
// https://doi.org/10.1109/TBME.1983.325207
|
|
float dot_0_2_plus_1_2 = dot_0_2 + dot_1_2;
|
|
float one_plus_dot_0_1 = 1.0f + dot_0_1;
|
|
float tangent = simplex_volume / (one_plus_dot_0_1 + dot_0_2_plus_1_2);
|
|
float triangle_solid_angle = 2.0f * positive_atan(tangent);
|
|
polygon.solid_angle += triangle_solid_angle;
|
|
polygon.fan_solid_angles[i] = polygon.solid_angle;
|
|
// Some intermediate results from above help us with sampling
|
|
polygon.triangle_parameters[i] = vec3(simplex_volume, dot_0_2_plus_1_2, one_plus_dot_0_1);
|
|
}
|
|
return polygon;
|
|
}
|
|
|
|
|
|
/*! An implementation of mix() using two fused-multiply add instructions. Used
|
|
because the native mix() implementation had stability issues in a few
|
|
spots. Credit to Fabian Giessen's blog, see:
|
|
https://fgiesen.wordpress.com/2012/08/15/linear-interpolation-past-present-and-future/
|
|
*/
|
|
float mix_fma(float x, float y, float a) {
|
|
return fma(a, y, fma(-a, x, x));
|
|
}
|
|
|
|
|
|
/*! Given the output of prepare_solid_angle_polygon_sampling(), this function
|
|
maps the given random numbers in the range from 0 to 1 to a normalized
|
|
direction vector providing a sample of the solid angle of the polygon in
|
|
the original space (used for arguments of
|
|
prepare_solid_angle_polygon_sampling()). Samples are distributed in
|
|
proportion to solid angle assuming uniform inputs.*/
|
|
vec3 sample_solid_angle_polygon(solid_angle_polygon_t polygon, vec2 random_numbers) {
|
|
// Decide which triangle needs to be sampled
|
|
float target_solid_angle = polygon.solid_angle * random_numbers[0];
|
|
float subtriangle_solid_angle = target_solid_angle;
|
|
vec3 parameters = polygon.triangle_parameters[0];
|
|
vec3 vertices[3] = {
|
|
polygon.vertex_dirs[1], polygon.vertex_dirs[0], polygon.vertex_dirs[2]
|
|
};
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT - 3; ++i) {
|
|
if (i + 3 >= polygon.vertex_count || polygon.fan_solid_angles[i] >= target_solid_angle) break;
|
|
subtriangle_solid_angle = target_solid_angle - polygon.fan_solid_angles[i];
|
|
vertices[0] = polygon.vertex_dirs[i + 2];
|
|
vertices[2] = polygon.vertex_dirs[i + 3];
|
|
parameters = polygon.triangle_parameters[i + 1];
|
|
}
|
|
// Construct a new vertex 2 on the arc between vertices 0 and 2 such that
|
|
// the resulting triangle has solid angle subtriangle_solid_angle
|
|
vec2 cos_sin = vec2(cos(0.5f * subtriangle_solid_angle), sin(0.5f * subtriangle_solid_angle));
|
|
vec3 offset = vertices[0] * (parameters[0] * cos_sin.x - parameters[1] * cos_sin.y) + vertices[2] * (parameters[2] * cos_sin.y);
|
|
vec3 new_vertex_2 = fma(2.0f * vec3(dot(vertices[0], offset) / dot(offset, offset)), offset, -vertices[0]);
|
|
// Now sample the line between vertex 1 and the newly created vertex 2
|
|
float s2 = dot(vertices[1], new_vertex_2);
|
|
float s = mix_fma(1.0f, s2, random_numbers[1]);
|
|
float denominator = fma(-s2, s2, 1.0f);
|
|
float t_normed = sqrt(fma(-s, s, 1.0f) / denominator);
|
|
// s2 may exceed one due to rounding error. random_numbers[1] is the
|
|
// limit of t_normed for s2 -> 1.
|
|
t_normed = (denominator > 0.0f) ? t_normed : random_numbers[1];
|
|
return fma(-t_normed, s2, s) * vertices[1] + t_normed * new_vertex_2;
|
|
}
|
|
|
|
|
|
/*! This structure carries intermediate results that only need to be computed
|
|
once per polygon and shading point to take samples proportional to
|
|
projected solid angle.*/
|
|
struct projected_solid_angle_polygon_t {
|
|
//! The number of vertices that form the polygon
|
|
uint vertex_count;
|
|
/*! The x- and y-coordinates of each polygon vertex in a coordinate system
|
|
where the normal is the z-axis. The vertices are sorted
|
|
counterclockwise.*/
|
|
vec2 vertices[MAX_POLYGON_VERTEX_COUNT];
|
|
/*! For each vertex in vertices, this vector describes the ellipse for the
|
|
next edge in counterclockwise direction. The last entry is meaningless,
|
|
except in the central case. For vertex 0, it holds the outer ellipse.
|
|
\see ellipse_from_edge() */
|
|
vec2 ellipses[MAX_POLYGON_VERTEX_COUNT];
|
|
//! The inner ellipse adjacent to vertex 0. If the x-component is positive,
|
|
//! the central case is present.
|
|
vec2 inner_ellipse_0;
|
|
/*! At index i, this array holds the projected solid angle of the polygon
|
|
in the sector between (sorted) vertices i and (i + 1) % vertex_count
|
|
In the central case, entry vertex_count - 1 is meaningful, otherwise
|
|
not.*/
|
|
float sector_projected_solid_angles[MAX_POLYGON_VERTEX_COUNT];
|
|
//! The total projected solid angle of the polygon
|
|
float projected_solid_angle;
|
|
};
|
|
|
|
|
|
//! Computes a * b - c * d with at most 1.5 ulps of error in the result. See
|
|
//! https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html or
|
|
//! Claude-Pierre Jeannerod, Nicolas Louvet and Jean-Michel Muller, 2013,
|
|
//! Further analysis of Kahan's algorithm for the accurate computation of 2x2
|
|
//! determinants, AMS Mathematics of Computation 82:284,
|
|
//! https://doi.org/10.1090/S0025-5718-2013-02679-8
|
|
float kahan(float a, float b, float c, float d) {
|
|
// Uncomment the line below to improve efficiency but reduce accuracy
|
|
// return a * b - c * d;
|
|
float cd = c * d;
|
|
float error = fma(c, d, -cd);
|
|
float result = fma(a, b, -cd);
|
|
return result - error;
|
|
}
|
|
|
|
|
|
//! Implements a cross product using Kahan's algorithm for every single entry,
|
|
//! i.e. the error in each output entry is at most 1.5 ulps
|
|
vec3 cross_stable(vec3 lhs, vec3 rhs) {
|
|
return vec3(
|
|
kahan(lhs.y, rhs.z, lhs.z, rhs.y),
|
|
kahan(lhs.z, rhs.x, lhs.x, rhs.z),
|
|
kahan(lhs.x, rhs.y, lhs.y, rhs.x)
|
|
);
|
|
}
|
|
|
|
|
|
//! \return The given vector, rotated 90 degrees counterclockwise around the
|
|
//! origin
|
|
vec2 rotate_90(vec2 input_vector) {
|
|
return vec2(-input_vector.y, input_vector.x);
|
|
}
|
|
|
|
|
|
//! \return true iff the given ellipse is marked as inner ellipse, i.e. iff the
|
|
//! spherical polygon that is bounded by it is further away from the zenith
|
|
//! than this ellipse.
|
|
bool is_inner_ellipse(vec2 ellipse) {
|
|
// If the implementation below causes you trouble, e.g. because you want to
|
|
// port to a different language, you may replace it by the commented line
|
|
// below but it will lead to seldom artifacts when ellipse.x == 0.0f.
|
|
// return ellipse.x < 0.0f;
|
|
// Extract the sign bit from ellipse.x (to be able to tell apart +0 and -0)
|
|
return (floatBitsToUint(ellipse.x) & 0x80000000) != 0;
|
|
}
|
|
|
|
|
|
//! \return true iff the given polygon contains the zenith (also known as
|
|
//! normal vector).
|
|
bool is_central_case(projected_solid_angle_polygon_t polygon) {
|
|
return polygon.inner_ellipse_0.x > 0.0f;
|
|
}
|
|
|
|
|
|
/*! Takes the great circle for the plane through the origin and the given two
|
|
points and constructs an ellipse for its projection to the xy-plane.
|
|
\return A vector ellipse such that a point is on the ellipse if and only if
|
|
dot(ellipse, point) * dot(ellipse, point) + dot(point, point) == 1.0f.
|
|
In other words, it is a normal vector of the great circle in half-
|
|
vector space. The sign bit of x encodes whether the edge runs clockwise
|
|
from vertex_0 to vertex_1 (inner ellipse) or not.
|
|
\see is_inner_ellipse() */
|
|
vec2 ellipse_from_edge(vec3 vertex_0, vec3 vertex_1) {
|
|
vec3 normal = cross_stable(vertex_0, vertex_1);
|
|
float scaling = 1.0f / normal.z;
|
|
scaling = is_inner_ellipse(normal.xy) ? -scaling : scaling;
|
|
vec2 ellipse = normal.xy * scaling;
|
|
// By convention, degenerate ellipses are outer ellipses, i.e. the first
|
|
// component is infinite
|
|
ellipse.x = (normal.z != 0.0f) ? ellipse.x : M_INFINITY;
|
|
return ellipse;
|
|
}
|
|
|
|
|
|
//! Transforms the given point using the matrix that characterizes the given
|
|
//! ellipse (as produced by ellipse_from_edge()). To be precise, this matrix is
|
|
//! identity + outerProduct(ellipse, ellipse).
|
|
vec2 ellipse_transform(vec2 ellipse, vec2 point) {
|
|
return fma(vec2(dot(ellipse, point)), ellipse, point);
|
|
}
|
|
|
|
|
|
//! Given an ellipse in the format produced by ellipse_from_edge(), this
|
|
//! function returns the determinant of the matrix characterizing this
|
|
//! ellipse.
|
|
float get_ellipse_det(vec2 ellipse) {
|
|
return fma(ellipse.x, ellipse.x, fma(ellipse.y, ellipse.y, 1.0f));
|
|
}
|
|
|
|
//! Returns the reciprocal square root of the ellipse determinant produced by
|
|
//! get_ellipse_det().
|
|
float get_ellipse_rsqrt_det(vec2 ellipse) {
|
|
return inversesqrt(get_ellipse_det(ellipse));
|
|
}
|
|
|
|
//! \return Reciprocal square of get_ellipse_direction_factor(ellipse, dir)
|
|
float get_ellipse_direction_factor_rsq(vec2 ellipse, vec2 dir) {
|
|
float ellipse_dot_dir = dot(ellipse, dir);
|
|
float dir_dot_dir = dot(dir, dir);
|
|
return fma(ellipse_dot_dir, ellipse_dot_dir, dir_dot_dir);
|
|
}
|
|
|
|
/*! Computes a factor by which a direction vector has to be multiplied to
|
|
obtain a point on the given ellipse.
|
|
\param ellipse An ellipse as produced by ellipse_from_edge().
|
|
\param dir The direction vector to be scaled onto the ellipse.
|
|
\return get_ellipse_direction_factor(ellipse, dir) * dir is a point on
|
|
the ellipse.*/
|
|
float get_ellipse_direction_factor(vec2 ellipse, vec2 dir) {
|
|
return inversesqrt(get_ellipse_direction_factor_rsq(ellipse, dir));
|
|
}
|
|
|
|
//! Like get_ellipse_direction_factor() but assumes that the given direction is
|
|
//! normalized. Faster.
|
|
float get_ellipse_normalized_direction_factor(vec2 ellipse, vec2 normalized_dir) {
|
|
float ellipse_dot_dir = dot(ellipse, normalized_dir);
|
|
return inversesqrt(fma(ellipse_dot_dir, ellipse_dot_dir, 1.0f));
|
|
}
|
|
|
|
|
|
//! Helper for get_area_between_ellipses_in_sector() and
|
|
//! sample_sector_between_ellipses()
|
|
float get_area_between_ellipses_in_sector_from_tangents(float inner_rsqrt_det, float inner_tangent, float outer_rsqrt_det, float outer_tangent) {
|
|
float inner_area = inner_rsqrt_det * positive_atan(inner_tangent);
|
|
float result = fma(outer_rsqrt_det, positive_atan(outer_tangent), -inner_area);
|
|
// Sort out NaNs and negative results
|
|
return (result > 0.0f) ? (0.5f * result) : 0.0f;
|
|
}
|
|
|
|
|
|
/*! Returns the signed area between the given outer and inner ellipses within
|
|
the sector enclosed by dir_0 and dir_1. Besides ellipses as produced by
|
|
ellipse_from_edge(), you also have to pass output of
|
|
get_ellipse_rsqrt_det(). Faster than calling get_ellipse_area_in_sector()
|
|
twice.*/
|
|
float get_area_between_ellipses_in_sector(vec2 inner_ellipse, float inner_rsqrt_det, vec2 outer_ellipse, float outer_rsqrt_det, vec2 dir_0, vec2 dir_1) {
|
|
float det_dirs = max(+0.0f, dot(dir_1, rotate_90(dir_0)));
|
|
float inner_dot = inner_rsqrt_det * dot(dir_0, ellipse_transform(inner_ellipse, dir_1));
|
|
float outer_dot = outer_rsqrt_det * dot(dir_0, ellipse_transform(outer_ellipse, dir_1));
|
|
return get_area_between_ellipses_in_sector_from_tangents(
|
|
inner_rsqrt_det, det_dirs / inner_dot,
|
|
outer_rsqrt_det, det_dirs / outer_dot);
|
|
}
|
|
|
|
|
|
/*! Computes the area for the intersection of the given ellipse and the sector
|
|
between the given two directions (going counterclockwise from dir_0 to
|
|
dir_1 for at most 180 degrees). The scaling of the directions is
|
|
irrelevant.
|
|
\see ellipse_from_edge() */
|
|
float get_ellipse_area_in_sector(vec2 ellipse, vec2 dir_0, vec2 dir_1) {
|
|
float ellipse_rsqrt_det = get_ellipse_rsqrt_det(ellipse);
|
|
float det_dirs = max(+0.0f, dot(dir_1, rotate_90(dir_0)));
|
|
float ellipse_dot = ellipse_rsqrt_det * dot(dir_0, ellipse_transform(ellipse, dir_1));
|
|
float area = 0.5f * ellipse_rsqrt_det * positive_atan(det_dirs / ellipse_dot);
|
|
// For degenerate ellipses, the result may be NaN but must be 0.0f
|
|
return (ellipse_rsqrt_det > 0.0f) ? area : 0.0f;
|
|
}
|
|
|
|
|
|
/*! Swaps vertices lhs and rhs (along with corresponding ellipses) of the given
|
|
polygon if the shorter path from lhs to rhs is clockwise. If the vertices
|
|
have identical directions in the xy-plane, vertices with degenerate
|
|
ellipses come first.
|
|
\note To avoid costly register spilling, lhs and rhs must be compile time
|
|
constants.*/
|
|
void compare_and_swap(inout projected_solid_angle_polygon_t polygon, uint lhs, uint rhs) {
|
|
vec2 lhs_copy = polygon.vertices[lhs];
|
|
// This line is designed to agree with the implementation of cross_stable
|
|
// for the z-coordinate, which determines if ellipses are inner or outer
|
|
float normal_z = kahan(lhs_copy.x, -polygon.vertices[rhs].y, lhs_copy.y, -polygon.vertices[rhs].x);
|
|
// Tie breaker: If both vertices are at the same angle (i.e. on a common
|
|
// great circle through the zenith), the one with the degenerate ellipse
|
|
// comes first
|
|
bool swap = (normal_z == 0.0f) ? isinf(polygon.ellipses[rhs].x) : (normal_z > 0.0f);
|
|
polygon.vertices[lhs] = swap ? polygon.vertices[rhs] : lhs_copy;
|
|
polygon.vertices[rhs] = swap ? lhs_copy : polygon.vertices[rhs];
|
|
lhs_copy = polygon.ellipses[lhs];
|
|
polygon.ellipses[lhs] = swap ? polygon.ellipses[rhs] : lhs_copy;
|
|
polygon.ellipses[rhs] = swap ? lhs_copy : polygon.ellipses[rhs];
|
|
}
|
|
|
|
|
|
//! Sorts the vertices of the given convex polygon counterclockwise using a
|
|
//! special sorting network. For non-convex polygons, the method may fail.
|
|
void sort_convex_polygon_vertices(inout projected_solid_angle_polygon_t polygon) {
|
|
if (polygon.vertex_count == 3) {
|
|
compare_and_swap(polygon, 1, 2);
|
|
}
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 4
|
|
else if (polygon.vertex_count == 4) {
|
|
compare_and_swap(polygon, 1, 3);
|
|
}
|
|
#endif
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 5
|
|
else if (polygon.vertex_count == 5) {
|
|
compare_and_swap(polygon, 2, 4);
|
|
compare_and_swap(polygon, 1, 3);
|
|
compare_and_swap(polygon, 1, 2);
|
|
compare_and_swap(polygon, 0, 3);
|
|
compare_and_swap(polygon, 3, 4);
|
|
}
|
|
#endif
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 6
|
|
else if (polygon.vertex_count == 6) {
|
|
compare_and_swap(polygon, 3, 5);
|
|
compare_and_swap(polygon, 2, 4);
|
|
compare_and_swap(polygon, 1, 5);
|
|
compare_and_swap(polygon, 0, 4);
|
|
compare_and_swap(polygon, 4, 5);
|
|
compare_and_swap(polygon, 1, 3);
|
|
}
|
|
#endif
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 7
|
|
else if (polygon.vertex_count == 7) {
|
|
compare_and_swap(polygon, 2, 5);
|
|
compare_and_swap(polygon, 1, 6);
|
|
compare_and_swap(polygon, 5, 6);
|
|
compare_and_swap(polygon, 3, 4);
|
|
compare_and_swap(polygon, 0, 4);
|
|
compare_and_swap(polygon, 4, 6);
|
|
compare_and_swap(polygon, 1, 3);
|
|
compare_and_swap(polygon, 3, 5);
|
|
compare_and_swap(polygon, 4, 5);
|
|
}
|
|
#endif
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 8
|
|
else if (polygon.vertex_count == 8) {
|
|
compare_and_swap(polygon, 2, 6);
|
|
compare_and_swap(polygon, 3, 7);
|
|
compare_and_swap(polygon, 1, 5);
|
|
compare_and_swap(polygon, 0, 4);
|
|
compare_and_swap(polygon, 4, 6);
|
|
compare_and_swap(polygon, 5, 7);
|
|
compare_and_swap(polygon, 6, 7);
|
|
compare_and_swap(polygon, 4, 5);
|
|
compare_and_swap(polygon, 1, 3);
|
|
}
|
|
#endif
|
|
// This comparison is shared by all sorting networks
|
|
compare_and_swap(polygon, 0, 2);
|
|
#if MAX_POLYGON_VERTEX_COUNT >= 4
|
|
if (polygon.vertex_count >= 4) {
|
|
// This comparison is shared by all sorting networks except the one for
|
|
// triangles
|
|
compare_and_swap(polygon, 2, 3);
|
|
}
|
|
#endif
|
|
// This comparison is shared by all sorting networks
|
|
compare_and_swap(polygon, 0, 1);
|
|
}
|
|
|
|
|
|
/*! Prepares all intermediate values to sample a convex polygon proportional to
|
|
projected solid angle.
|
|
\param vertex_count Number of vertices forming the polygon (at least 3).
|
|
\param vertices List of vertex locations in a coordinate system where the
|
|
shading position is the origin and the normal is the z-axis. The
|
|
polygon should be already clipped against the plane z=0. If
|
|
vertex_count < MAX_POLYGON_VERTEX_COUNT, the first vertex has to be
|
|
repeated at vertex_count. They need not be normalized but if you
|
|
encounter issues with under- or overflow (e.g. NaN or INF outputs),
|
|
normalization may help. The polygon must be convex, and the winding of
|
|
the vertices as seen from the origin must be clockwise. No three
|
|
vertices should be collinear.
|
|
\return Intermediate values for sampling.*/
|
|
projected_solid_angle_polygon_t prepare_projected_solid_angle_polygon_sampling(uint vertex_count, vec3 vertices[MAX_POLYGON_VERTEX_COUNT]) {
|
|
projected_solid_angle_polygon_t polygon;
|
|
// Copy vertices and assign ellipses
|
|
polygon.vertex_count = vertex_count;
|
|
polygon.inner_ellipse_0 = vec2(1.0f, 0.0f);
|
|
polygon.vertices[0] = vertices[0].xy;
|
|
polygon.ellipses[0] = ellipse_from_edge(vertices[0], vertices[1]);
|
|
vec2 previous_ellipse = polygon.ellipses[0];
|
|
[[unroll]]
|
|
for (uint i = 1; i != MAX_POLYGON_VERTEX_COUNT; ++i) {
|
|
polygon.vertices[i] = vertices[i].xy;
|
|
if (i > 2 && i == polygon.vertex_count) break;
|
|
vec2 ellipse = ellipse_from_edge(vertices[i], vertices[(i + 1) % MAX_POLYGON_VERTEX_COUNT]);
|
|
bool ellipse_inner = is_inner_ellipse(ellipse);
|
|
// If the edge is an inner edge, the order is going to flip
|
|
polygon.ellipses[i] = ellipse_inner ? previous_ellipse : ellipse;
|
|
// In doing so, we drop one ellipse, unless we store it explicitly
|
|
polygon.inner_ellipse_0 = (is_inner_ellipse(previous_ellipse) && !ellipse_inner) ? previous_ellipse : polygon.inner_ellipse_0;
|
|
previous_ellipse = ellipse;
|
|
}
|
|
// Same thing for the first vertex (i.e. here we close the loop)
|
|
vec2 ellipse = polygon.ellipses[0];
|
|
bool ellipse_inner = is_inner_ellipse(ellipse);
|
|
polygon.ellipses[0] = ellipse_inner ? previous_ellipse : ellipse;
|
|
polygon.inner_ellipse_0 = (is_inner_ellipse(previous_ellipse) && !ellipse_inner) ? previous_ellipse : polygon.inner_ellipse_0;
|
|
// Compute projected solid angles per sector and in total
|
|
polygon.projected_solid_angle = 0.0f;
|
|
if (is_central_case(polygon)) {
|
|
// In the central case, we have polygon.vertex_count sectors, each
|
|
// bounded by a single ellipse
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT; ++i) {
|
|
if (i > 2 && i == polygon.vertex_count) break;
|
|
polygon.sector_projected_solid_angles[i] = get_ellipse_area_in_sector(polygon.ellipses[i], polygon.vertices[i], polygon.vertices[(i + 1) % MAX_POLYGON_VERTEX_COUNT]);
|
|
polygon.projected_solid_angle += polygon.sector_projected_solid_angles[i];
|
|
}
|
|
}
|
|
else {
|
|
// Sort vertices counter clockwise
|
|
sort_convex_polygon_vertices(polygon);
|
|
// There are polygon.vertex_count - 1 sectors, each bounded by an inner
|
|
// and an outer ellipse
|
|
vec2 inner_ellipse = polygon.inner_ellipse_0;
|
|
float inner_rsqrt_det = get_ellipse_rsqrt_det(inner_ellipse);
|
|
vec2 outer_ellipse;
|
|
float outer_rsqrt_det;
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT - 1; ++i) {
|
|
if (i > 1 && i + 1 == polygon.vertex_count) break;
|
|
vec2 vertex_ellipse = polygon.ellipses[i];
|
|
bool vertex_inner = is_inner_ellipse(vertex_ellipse);
|
|
float vertex_rsqrt_det = get_ellipse_rsqrt_det(vertex_ellipse);
|
|
if (i == 0) {
|
|
outer_ellipse = vertex_ellipse;
|
|
outer_rsqrt_det = vertex_rsqrt_det;
|
|
}
|
|
else {
|
|
inner_ellipse = vertex_inner ? vertex_ellipse : inner_ellipse;
|
|
inner_rsqrt_det = vertex_inner ? vertex_rsqrt_det : inner_rsqrt_det;
|
|
outer_ellipse = vertex_inner ? outer_ellipse : vertex_ellipse;
|
|
outer_rsqrt_det = vertex_inner ? outer_rsqrt_det : vertex_rsqrt_det;
|
|
}
|
|
polygon.sector_projected_solid_angles[i] = get_area_between_ellipses_in_sector(
|
|
inner_ellipse, inner_rsqrt_det, outer_ellipse, outer_rsqrt_det, polygon.vertices[i], polygon.vertices[i + 1]);
|
|
polygon.projected_solid_angle += polygon.sector_projected_solid_angles[i];
|
|
}
|
|
}
|
|
return polygon;
|
|
}
|
|
|
|
|
|
/*! \return A scalar multiple of rhs that is not too far from being normalized.
|
|
For the result, length() returns something between sqrt(2.0f) and 8.0f.
|
|
The sign gets flipped such that the dot product of semi_circle and the
|
|
result is non-negative.
|
|
\note Introduces less latency than normalize() and does not use special
|
|
functions. Useful to avoid under- and overflow when working with
|
|
homogeneous coordinates. The result is undefined if rhs is zero.*/
|
|
vec2 normalize_approx_and_flip(vec2 rhs, vec2 semi_circle) {
|
|
float scaling = abs(rhs.x) + abs(rhs.y);
|
|
// By flipping each bit on the exponent E, we turn it into 1 - E, which is
|
|
// close enough to a reciprocal.
|
|
scaling = uintBitsToFloat(floatBitsToUint(scaling) ^ 0x7F800000u);
|
|
// If the line above causes you any sort of trouble (e.g. because you want
|
|
// to port the code to another language or you are doing differentiable
|
|
// rendering), just use this one instead:
|
|
// scaling = 1.0f / scaling;
|
|
// Flip the sign as needed
|
|
scaling = (dot(rhs, semi_circle) >= 0.0f) ? scaling : -scaling;
|
|
return scaling * rhs;
|
|
}
|
|
|
|
|
|
/*! Returns a solution to the given homogeneous quadratic equation, i.e. a
|
|
non-zero vector root such that dot(root, quadratic * root) == 0.0f. The
|
|
returned root depends continuously on quadratic. Pass -quadratic if you
|
|
want the other root.
|
|
\note The implementation is as proposed by Blinn, except that we do not
|
|
have a special case for quadratic[0][1] + quadratic[1][0] == 0.0f. Unlike
|
|
the standard quadratic formula, it allows us to postpone a division and is
|
|
stable in all cases.
|
|
James F. Blinn 2006, How to Solve a Quadratic Equation, Part 2, IEEE
|
|
Computer Graphics and Applications 26:2 https://doi.org/10.1109/MCG.2006.35
|
|
*/
|
|
vec2 solve_homogeneous_quadratic(mat2 quadratic) {
|
|
float coeff_xy = 0.5f * (quadratic[0][1] + quadratic[1][0]);
|
|
float sqrt_discriminant = sqrt(max(0.0f, coeff_xy * coeff_xy - quadratic[0][0] * quadratic[1][1]));
|
|
float scaled_root = abs(coeff_xy) + sqrt_discriminant;
|
|
return (coeff_xy >= 0.0f) ? vec2(scaled_root, -quadratic[0][0]) : vec2(quadratic[1][1], scaled_root);
|
|
}
|
|
|
|
|
|
/*! Generates a sample between two ellipses and in a specified sector. The
|
|
sample is distributed uniformly with respect to the area measure.
|
|
\param random_numbers A pair of independent uniform random numbers on [0,1]
|
|
\param target_area random_numbers[0] multiplied by the projected solid
|
|
angle of the area to be sampled.
|
|
\param inner_ellipse, outer_ellipse The inner and outer ellipse, as
|
|
produced by ellipse_from_edge().
|
|
\param dir_0, dir_1 Two direction vectors bounding the sector. They need
|
|
not be normalized.
|
|
\param iteration_count The number of iterations to perform. Lower values
|
|
trade speed for bias. Two iterations give practically no bias.
|
|
\return The sample in Cartesian coordinates.*/
|
|
vec2 sample_sector_between_ellipses(vec2 random_numbers, float target_area, vec2 inner_ellipse, vec2 outer_ellipse, vec2 dir_0, vec2 dir_1, uint iteration_count) {
|
|
// For the initialization, split the sector in half
|
|
vec2 quad_dirs[3];
|
|
quad_dirs[0] = normalize(dir_0);
|
|
quad_dirs[2] = normalize(dir_1);
|
|
quad_dirs[1] = quad_dirs[0] + quad_dirs[2];
|
|
// Compute where these lines intersect the ellipses. The six intersection
|
|
// points define two adjacent quads.
|
|
float normalization_factor[2][3] = {
|
|
{
|
|
get_ellipse_normalized_direction_factor(inner_ellipse, quad_dirs[0]),
|
|
get_ellipse_direction_factor(inner_ellipse, quad_dirs[1]),
|
|
get_ellipse_normalized_direction_factor(inner_ellipse, quad_dirs[2])
|
|
},
|
|
{
|
|
get_ellipse_normalized_direction_factor(outer_ellipse, quad_dirs[0]),
|
|
get_ellipse_direction_factor(outer_ellipse, quad_dirs[1]),
|
|
get_ellipse_normalized_direction_factor(outer_ellipse, quad_dirs[2])
|
|
}
|
|
};
|
|
// Compute the relative size of the areas inside these quads
|
|
float sector_areas[2] = {
|
|
normalization_factor[1][0] * normalization_factor[1][1] - normalization_factor[0][0] * normalization_factor[0][1],
|
|
normalization_factor[1][1] * normalization_factor[1][2] - normalization_factor[0][1] * normalization_factor[0][2]
|
|
};
|
|
// Now pick which of the two quads should be sampled for the
|
|
// initialization. If it is not the second, we move data such that the
|
|
// relevant array indices are 1 and 2 anyway.
|
|
float target_quad_area = mix_fma(-sector_areas[0], sector_areas[1], random_numbers[0]);
|
|
quad_dirs[2] = (target_quad_area <= 0.0f) ? quad_dirs[0] : quad_dirs[2];
|
|
normalization_factor[0][2] = (target_quad_area <= 0.0f) ? normalization_factor[0][0] : normalization_factor[0][2];
|
|
normalization_factor[1][2] = (target_quad_area <= 0.0f) ? normalization_factor[1][0] : normalization_factor[1][2];
|
|
target_quad_area += (target_quad_area <= 0.0f) ? sector_areas[0] : -sector_areas[1];
|
|
// We have been a bit lazy about area computation before but now we need
|
|
// all the factors (except for a factor of 0.5 that cancels with a 2 later)
|
|
target_quad_area *= abs(determinant(mat2(quad_dirs[1], quad_dirs[2])));
|
|
// Construct normal vectors for the inner and outer edge of the selected
|
|
// quad. We construct the normal like a half vector (i.e. by addition)
|
|
// because it is less prone to cancellation than an approach using the edge
|
|
// direction (i.e. subtraction of sometimes nearly identical vectors)
|
|
vec2 quad_normals[2] = {
|
|
quad_dirs[1] * normalization_factor[0][1] + quad_dirs[2] * normalization_factor[0][2],
|
|
quad_dirs[1] * normalization_factor[1][1] + quad_dirs[2] * normalization_factor[1][2]
|
|
};
|
|
quad_normals[0] = ellipse_transform(inner_ellipse, quad_normals[0]);
|
|
quad_normals[1] = ellipse_transform(outer_ellipse, quad_normals[1]);
|
|
// Construct complete line equations
|
|
float quad_offsets[2] = {
|
|
dot(quad_normals[0], quad_dirs[1]) * normalization_factor[0][1],
|
|
dot(quad_normals[1], quad_dirs[1]) * normalization_factor[1][1]
|
|
};
|
|
// Now sample the direction within the selected quad by constructing a
|
|
// quadratic equation. This is the initialization for the iteration.
|
|
mat2 quadratic = outerProduct((quad_offsets[1] * normalization_factor[1][2]) * rotate_90(quad_dirs[2]), quad_normals[0]);
|
|
quadratic -= outerProduct((quad_offsets[0] * normalization_factor[0][2]) * rotate_90(quad_dirs[2]) + target_quad_area * quad_normals[0], quad_normals[1]);
|
|
vec2 current_dir = solve_homogeneous_quadratic(quadratic);
|
|
|
|
#ifndef USE_BIASED_PROJECTED_SOLID_ANGLE_SAMPLING
|
|
// For boundary values, the initialization is perfect but the iteration may
|
|
// be unstable, so we disable it
|
|
float acceptable_error = 1.0e-5f;
|
|
iteration_count = (abs(random_numbers.x - 0.5f) <= 0.5f - acceptable_error) ? iteration_count : 0;
|
|
|
|
// Now refine this initialization iteratively
|
|
float inner_rsqrt_det = get_ellipse_rsqrt_det(inner_ellipse);
|
|
float outer_rsqrt_det = get_ellipse_rsqrt_det(outer_ellipse);
|
|
[[dont_unroll]]
|
|
for (uint i = 0; i != iteration_count; ++i) {
|
|
// Avoid under- or overflow and flip the sign so that the clamping to
|
|
// zero below makes sense
|
|
current_dir = normalize_approx_and_flip(current_dir, quad_dirs[1]);
|
|
// Transform current_dir using both ellipses
|
|
vec2 inner_dir = ellipse_transform(inner_ellipse, current_dir);
|
|
vec2 outer_dir = ellipse_transform(outer_ellipse, current_dir);
|
|
// Evaluate the objective function (reusing inner_dir and outer_dir)
|
|
float det_dirs = max(+0.0f, dot(current_dir, rotate_90(quad_dirs[0])));
|
|
float error = target_area - get_area_between_ellipses_in_sector_from_tangents(
|
|
inner_rsqrt_det, det_dirs / (inner_rsqrt_det * dot(quad_dirs[0], inner_dir)),
|
|
outer_rsqrt_det, det_dirs / (outer_rsqrt_det * dot(quad_dirs[0], outer_dir)));
|
|
// Construct a homogeneous quadratic whose solutions include the next
|
|
// step of the iteration
|
|
quadratic = outerProduct(inner_dir - outer_dir, rotate_90(current_dir)) - outerProduct((2.0f * error) * inner_dir, outer_dir);
|
|
current_dir = solve_homogeneous_quadratic(quadratic);
|
|
}
|
|
#endif
|
|
|
|
// The halved sector is at most 90 degrees large, so the dot product with
|
|
// the half vector has to be positive
|
|
current_dir = (dot(current_dir, quad_dirs[1]) >= 0.0f) ? current_dir : -current_dir;
|
|
// Sample a squared radius uniformly between the two ellipses
|
|
float inner_factor = 1.0f / get_ellipse_direction_factor_rsq(inner_ellipse, current_dir);
|
|
float outer_factor = 1.0f / get_ellipse_direction_factor_rsq(outer_ellipse, current_dir);
|
|
current_dir *= sqrt(mix_fma(inner_factor, outer_factor, random_numbers[1]));
|
|
return current_dir;
|
|
}
|
|
|
|
|
|
/*! Produces a sample in the solid angle of the given polygon. If the random
|
|
numbers are uniform in [0,1]^2, the sample is uniform in the projected
|
|
solid angle of the polygon.
|
|
\param polygon Output of prepare_projected_solid_angle_polygon_sampling().
|
|
\param random_numbers A uniform point in [0,1]^2.
|
|
\return A sample on the upper hemisphere (i.e. z>=0) in Cartesian
|
|
coordinates.*/
|
|
vec3 sample_projected_solid_angle_polygon(projected_solid_angle_polygon_t polygon, vec2 random_numbers) {
|
|
float target_projected_solid_angle = random_numbers[0] * polygon.projected_solid_angle;
|
|
// Distinguish between the central case
|
|
vec3 sampled_dir;
|
|
vec2 outer_ellipse;
|
|
vec2 dir_0;
|
|
if (is_central_case(polygon)) {
|
|
// Select a sector and copy the relevant attributes
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT; ++i) {
|
|
if (i > 0) {
|
|
target_projected_solid_angle -= polygon.sector_projected_solid_angles[i - 1];
|
|
}
|
|
outer_ellipse = polygon.ellipses[i];
|
|
dir_0 = polygon.vertices[i];
|
|
if ((i >= 2 && i + 1 == polygon.vertex_count) || target_projected_solid_angle < polygon.sector_projected_solid_angles[i])
|
|
break;
|
|
}
|
|
// Sample a direction within the sector
|
|
float sqrt_det = sqrt(get_ellipse_det(outer_ellipse));
|
|
float angle = 2.0f * target_projected_solid_angle * sqrt_det;
|
|
sampled_dir.xy = (cos(angle) * sqrt_det) * dir_0 + sin(angle) * rotate_90(ellipse_transform(outer_ellipse, dir_0));
|
|
// Sample a squared radius uniformly within the ellipse
|
|
sampled_dir.xy *= sqrt(random_numbers[1] / get_ellipse_direction_factor_rsq(outer_ellipse, sampled_dir.xy));
|
|
}
|
|
// And the decentral case
|
|
else {
|
|
// Select a sector and copy the relevant attributes
|
|
float sector_projected_solid_angle;
|
|
vec2 inner_ellipse = polygon.inner_ellipse_0;
|
|
vec2 dir_1;
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT - 1; ++i) {
|
|
vec2 vertex_ellipse = polygon.ellipses[i];
|
|
if (i == 0) {
|
|
outer_ellipse = vertex_ellipse;
|
|
}
|
|
else {
|
|
target_projected_solid_angle -= polygon.sector_projected_solid_angles[i - 1];
|
|
bool vertex_inner = is_inner_ellipse(vertex_ellipse);
|
|
inner_ellipse = vertex_inner ? vertex_ellipse : inner_ellipse;
|
|
outer_ellipse = vertex_inner ? outer_ellipse : vertex_ellipse;
|
|
}
|
|
dir_0 = polygon.vertices[i];
|
|
dir_1 = polygon.vertices[i + 1];
|
|
sector_projected_solid_angle = polygon.sector_projected_solid_angles[i];
|
|
if ((i >= 1 && i + 2 == polygon.vertex_count) || target_projected_solid_angle < sector_projected_solid_angle)
|
|
break;
|
|
}
|
|
// Sample it
|
|
random_numbers[0] = target_projected_solid_angle / sector_projected_solid_angle;
|
|
sampled_dir.xy = sample_sector_between_ellipses(random_numbers, target_projected_solid_angle, inner_ellipse, outer_ellipse, dir_0, dir_1, 2);
|
|
}
|
|
// Construct the sample
|
|
sampled_dir.z = sqrt(max(0.0f, fma(-sampled_dir.x, sampled_dir.x, fma(-sampled_dir.y, sampled_dir.y, 1.0f))));
|
|
return sampled_dir;
|
|
}
|
|
|
|
|
|
/*! Determines the error of a sample from the projected solid angle of a
|
|
polygon due to the iterative procedure.
|
|
\param polygon Output of prepare_projected_solid_angle_polygon_sampling().
|
|
\param random_numbers The value passed for random_numbers in
|
|
sample_projected_solid_angle_polygon().
|
|
\param sampled_dir The direction returned by
|
|
sample_projected_solid_angle_polygon().
|
|
\return The first component is the signed backward error: The difference
|
|
between random_number_0 and the random number that would yield
|
|
sampled_dir with perfect computation. The second component is the
|
|
backward error, multiplied by the projected solid angle of the polygon.
|
|
The third component is the signed forward error (at least a first-order
|
|
estimate of it): The difference between the exact result of sampling
|
|
and the actual result in radians. Note that all of these are subject
|
|
to rounding error themselves. In the central case, zero is returned.*/
|
|
vec3 compute_projected_solid_angle_polygon_sampling_error(projected_solid_angle_polygon_t polygon, vec2 random_numbers, vec3 sampled_dir) {
|
|
float target_projected_solid_angle = random_numbers[0] * polygon.projected_solid_angle;
|
|
// In the central case, the sampling procedure is exact except for rounding
|
|
// error
|
|
if (is_central_case(polygon)) {
|
|
return vec3(0.0f);
|
|
}
|
|
// In the other case, we repeat some computations to find the error
|
|
else {
|
|
// Select a sector and copy the relevant attributes
|
|
float sector_projected_solid_angle;
|
|
vec2 outer_ellipse;
|
|
vec2 inner_ellipse = polygon.inner_ellipse_0;
|
|
vec2 dir_0;
|
|
[[unroll]]
|
|
for (uint i = 0; i != MAX_POLYGON_VERTEX_COUNT - 1; ++i) {
|
|
if ((i > 1 && i + 1 == polygon.vertex_count) || (i > 0 && target_projected_solid_angle < 0.0f))
|
|
break;
|
|
sector_projected_solid_angle = polygon.sector_projected_solid_angles[i];
|
|
target_projected_solid_angle -= sector_projected_solid_angle;
|
|
vec2 vertex_ellipse = polygon.ellipses[i];
|
|
bool vertex_inner = is_inner_ellipse(vertex_ellipse);
|
|
if (i == 0) {
|
|
outer_ellipse = vertex_ellipse;
|
|
}
|
|
else {
|
|
inner_ellipse = vertex_inner ? vertex_ellipse : inner_ellipse;
|
|
outer_ellipse = vertex_inner ? outer_ellipse : vertex_ellipse;
|
|
}
|
|
dir_0 = polygon.vertices[i];
|
|
}
|
|
target_projected_solid_angle += sector_projected_solid_angle;
|
|
// Compute the area in the sector from sampled_dir and dir_0 between
|
|
// the two ellipses
|
|
float sampled_projected_solid_angle = get_area_between_ellipses_in_sector(
|
|
inner_ellipse, get_ellipse_rsqrt_det(inner_ellipse),
|
|
outer_ellipse, get_ellipse_rsqrt_det(outer_ellipse),
|
|
dir_0, sampled_dir.xy);
|
|
// Compute error in the projected solid angle
|
|
float scaled_backward_error = target_projected_solid_angle - sampled_projected_solid_angle;
|
|
float backward_error = scaled_backward_error / polygon.projected_solid_angle;
|
|
// Evaluate the derivative of the sampled direction with respect to the
|
|
// projected solid angle
|
|
mat2 constraint_matrix;
|
|
vec2 inner_dir = ellipse_transform(inner_ellipse, sampled_dir.xy);
|
|
vec2 outer_dir = ellipse_transform(outer_ellipse, sampled_dir.xy);
|
|
float inner_factor = 1.0f / dot(sampled_dir.xy, inner_dir);
|
|
float outer_factor = 1.0f / dot(sampled_dir.xy, outer_dir);
|
|
constraint_matrix[0] = 0.5f * (inner_factor - outer_factor) * rotate_90(sampled_dir.xy);
|
|
constraint_matrix[1] = ((1.0f - random_numbers[1]) / (inner_factor * inner_factor)) * inner_dir;
|
|
constraint_matrix[1] += (random_numbers[1] / (outer_factor * outer_factor)) * outer_dir;
|
|
constraint_matrix = transpose(constraint_matrix);
|
|
vec3 sample_derivative;
|
|
sample_derivative.xy = (1.0f / determinant(constraint_matrix)) * vec2(constraint_matrix[1][1], -constraint_matrix[0][1]);
|
|
sample_derivative.z = -dot(sampled_dir.xy, sample_derivative.xy) / sampled_dir.z;
|
|
// Evaluate the forward error in radians
|
|
float forward_error = length(sample_derivative) * scaled_backward_error;
|
|
// Return both errors
|
|
return vec3(backward_error, scaled_backward_error, forward_error);
|
|
}
|
|
}
|