Cycles: new Microfacet-based Hair BSDF with elliptical cross-section support #105600

Merged
Weizhen Huang merged 114 commits from weizhen/blender:microfacet_hair into main 2023-08-18 12:46:20 +02:00
3 changed files with 122 additions and 51 deletions
Showing only changes of commit 5015c8219b - Show all commits

View File

@ -710,7 +710,6 @@ void calculate_tangents(Span<float3> positions,
void calculate_normals(Span<float3> positions,
bool is_cyclic,
int resolution,
Span<float3> tangents,
MutableSpan<float3> normals);
/**

View File

@ -6,6 +6,8 @@
#include "BKE_attribute_math.hh"
#include "BKE_curves.hh"
#include "BLI_math_rotation_legacy.hh"
#include "BLI_math_vector.hh"
namespace blender::bke::curves::catmull_rom {
@ -139,7 +141,8 @@ static void interpolate_to_evaluated(const int order,
}
template<typename T>
static void interpolate_to_evaluated(const Span<T> src,
static void interpolate_to_evaluated(const int order,
const Span<T> src,
const bool cyclic,
const int resolution,
MutableSpan<T> dst)
@ -147,7 +150,7 @@ static void interpolate_to_evaluated(const Span<T> src,
{
BLI_assert(dst.size() == calculate_evaluated_num(src.size(), cyclic, resolution));
interpolate_to_evaluated(
0,
order,
src,
cyclic,
[resolution](const int segment_i) -> IndexRange {
@ -180,7 +183,7 @@ void interpolate_to_evaluated(const GSpan src,
{
attribute_math::convert_to_static_type(src.type(), [&](auto dummy) {
using T = decltype(dummy);
interpolate_to_evaluated(src.typed<T>(), cyclic, resolution, dst.typed<T>());
interpolate_to_evaluated(0, src.typed<T>(), cyclic, resolution, dst.typed<T>());
});
}
@ -206,25 +209,49 @@ void calculate_tangents(const Span<float3> positions,
return;
}
interpolate_to_evaluated(
1,
positions,
is_cyclic,
[resolution](const int segment_i) -> IndexRange {
return {segment_i * resolution, resolution};
},
evaluated_tangents);
interpolate_to_evaluated(1, positions, is_cyclic, resolution, evaluated_tangents);
for (const int i : evaluated_tangents.index_range()) {
evaluated_tangents[i] = math::normalize(evaluated_tangents[i]);
}
}
/* Calculate analytical normals from the first and the second derivative. */
/* Use a combination of Newton iterations and bisection to find the root in the monotonic interval
* [x_begin, x_end], given precision. Idea taken from the paper High-Performance Polynomial Root
* Finding for Graphics by Cem Yuksel, HPG 2022. */
/* TODO: move this function elsewhere as a utility function? Also test corner cases. */
static float find_root_newton_bisection(float x_begin,
float x_end,
std::function<float(float)> eval,
std::function<float(float)> eval_derivative,
const float precision)
{
float x_mid = 0.5f * (x_begin + x_end);
float y_begin = eval(x_begin);
BLI_assert(signf(y_begin) != signf(eval(x_end)));
while (x_mid - x_begin > precision) {
const float y_mid = eval(x_mid);
if (signf(y_begin) == signf(y_mid)) {
x_begin = x_mid;
y_begin = y_mid;
}
else {
x_end = x_mid;
}
const float x_newton = x_mid - y_mid / eval_derivative(x_mid);
x_mid = (x_newton > x_begin && x_newton < x_end) ? x_newton : 0.5f * (x_begin + x_end);
}
return x_mid;
}
/* Calculate normals by minimizing the potential energy due to twist and bending. Global
* estimation involves integration and is thus too costly. Instead, we start with the first
* curvature vector and propogate it along the curve. */
void calculate_normals(const Span<float3> positions,
const bool is_cyclic,
const int resolution,
const Span<float3> evaluated_tangents,
MutableSpan<float3> evaluated_normals)
{
if (evaluated_normals.size() == 1) {
@ -232,48 +259,94 @@ void calculate_normals(const Span<float3> positions,
return;
}
/* Compute the second derivative (r'') of the control points, evaluated at t = 0. */
interpolate_to_evaluated(
2,
positions,
is_cyclic,
[resolution](const int segment_i) -> IndexRange {
return {segment_i * resolution, 1};
},
evaluated_normals);
/* TODO: check if derivatives == 0.*/
Vector<float3> first_derivatives_(evaluated_normals.size());
MutableSpan<float3> first_derivatives = first_derivatives_;
interpolate_to_evaluated(1, positions, is_cyclic, resolution, first_derivatives);
/* The normals along the segment is interpolated between the control points. */
const float step = 1.0f / resolution;
for (const int i : positions.index_range()) {
if (i == positions.size() - 1 && !is_cyclic) {
break;
}
const float3 last_normal = evaluated_normals[i * resolution];
float3 next_normal = (i == positions.size() - 1) ? evaluated_normals[0] :
evaluated_normals[(i + 1) * resolution];
if (!is_cyclic && math::dot(last_normal, next_normal) < 0) {
/* Prevent sudden normal changes. */
next_normal = -next_normal;
evaluated_normals[(i + 1) * resolution] = next_normal;
}
for (int j = 1; j < resolution; j++) {
const float t = j * step;
evaluated_normals[i * resolution + j] = (1.0f - t) * last_normal + t * next_normal;
}
}
Vector<float3> second_derivatives_(evaluated_normals.size());
MutableSpan<float3> second_derivatives = second_derivatives_;
interpolate_to_evaluated(2, positions, is_cyclic, resolution, second_derivatives);
MutableSpan<float3> binormals(evaluated_normals);
/* TODO: below are hair-specific values, assuming elliptical cross-section. Maybe make them more
* general by specifying user-defined weight? */
/* Values taken from Table 9.5 in book Chemical and Physical Behavior of Human Hair by Clarence
* R. Robbins, 5th edition. Unit: GPa. */
const float youngs_modulus = 3.89f;
const float shear_modulus = 0.89f;
/* Assuming major axis a = 1. */
const float aspect_ratio = 0.5f; /* Minimal realistic value. */
const float aspect_ratio_squared = aspect_ratio * aspect_ratio;
const float torsion_constant = M_PI * aspect_ratio_squared * aspect_ratio /
(1.0f + aspect_ratio_squared);
const float moment_of_inertia_coefficient = M_PI * aspect_ratio * 0.25f *
(1.0f - aspect_ratio_squared);
/* Compute evaluated_normals of the control points. */
const float dt = 1.0f / resolution;
/* Compute evaluated_normals. */
for (const int i : evaluated_normals.index_range()) {
/* B = normalize(r' x r'') */
binormals[i] = math::cross(evaluated_tangents[i], evaluated_normals[i]);
/* B = r' x r'' */
const float3 binormal = math::cross(first_derivatives[i], second_derivatives[i]);
/* N = B x T */
evaluated_normals[i] = math::normalize(math::cross(binormals[i], evaluated_tangents[i]));
}
/* TODO: Catmull-Rom is not C2 continuous. Some smoothening maybe? */
float3 curvature_vector = math::normalize(math::cross(binormal, first_derivatives[i]));
if (positions.size() == 1) {
return;
if (i == 0) {
/* TODO: Does it make sense to propogate from where the curvature is the largest? */
evaluated_normals[i] = curvature_vector;
continue;
}
const float first_derivative_norm = math::length(first_derivatives[i]);
/* Arc length depends on the unit, assuming meter. */
const float arc_length = first_derivative_norm * dt;
const float curvature = math::length(binormal) / pow3f(first_derivative_norm);
const float weight_twist = 0.5f * shear_modulus * torsion_constant / arc_length;
const float weight_bending = 0.5f * arc_length * curvature * curvature * youngs_modulus *
moment_of_inertia_coefficient;
const float3 last_tangent = math::normalize(first_derivatives[i - 1]);
const float3 current_tangent = first_derivatives[i] / first_derivative_norm;
const float angle = angle_normalized_v3v3(last_tangent, current_tangent);
const float3 last_normal = evaluated_normals[i - 1];
const float3 minimal_twist_normal =
angle > 0 ?
math::rotate_direction_around_axis(
last_normal, math::normalize(math::cross(last_tangent, current_tangent)), angle) :
last_normal;
/* Angle between the curvature vector and the minimal twist normal. */
float theta_t = angle_normalized_v3v3(curvature_vector, minimal_twist_normal);
if (theta_t > M_PI_2) {
curvature_vector = -curvature_vector;
theta_t = M_PI - theta_t;
}
/* Minimizing the total potential energy U(θ) = wt * (θ - θt)^2 + wb * sin^2(θ) by solving the
* equation: 2 * wt * (θ - θt) + wb * sin(2θ) = 0, with θ being the angle between the computed
* normal and the curvature vector. */
const float theta_begin = 0.0f;
const float theta_end = weight_twist + weight_bending * cosf(2.0f * theta_t) < 0 ?
0.5f * acosf(-weight_twist / weight_bending) :
theta_t;
const float theta = find_root_newton_bisection(
theta_begin,
theta_end,
[weight_bending, weight_twist, theta_t](float t) -> float {
return 2.0f * weight_twist * (t - theta_t) + weight_bending * sinf(2.0f * t);
},
[weight_bending, weight_twist](float t) -> float {
return 2.0f * (weight_twist + weight_bending * cosf(2.0f * t));
},
1e-5f);
const float3 axis = math::dot(current_tangent,
math::cross(curvature_vector, minimal_twist_normal)) > 0 ?
current_tangent :
-current_tangent;
evaluated_normals[i] = math::rotate_direction_around_axis(curvature_vector, axis, theta);
}
}

View File

@ -784,7 +784,6 @@ Span<float3> CurvesGeometry::evaluated_normals() const
curves::catmull_rom::calculate_normals(positions.slice(points),
cyclic[curve_index],
resolution[curve_index],
evaluated_tangents.slice(evaluated_points),
evaluated_normals.slice(evaluated_points));
break;
default: