blender
/
blender
128
395
Fork 567
Code Issues 6.2k Pull Requests 750 Projects 19 Wiki Activity

Compute hair normal with minimal total potential energy 

Although computed in the geometry node, this normal is specific for
hairs is only implemented for Catmull-Rom splines. Questions remain
whether this should be generalized, also a few design ideas are obscure.
{F14115834}
This commit is contained in:
Weizhen Huang 2023-01-05 17:46:59 +01:00
parent 8af16fd087
commit 5015c8219b
3 changed files with 122 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
source/blender/blenkernel/BKE_curves.hh
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);
/**

171
source/blender/blenkernel/intern/curve_catmull_rom.cc
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);
}
}

1
source/blender/blenkernel/intern/curves_geometry.cc
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: