[Bf-blender-cvs] [5015c8219b6] microfacet_hair: 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}
Weizhen Huang
noreply at git.blender.org
Thu Jan 5 18:04:51 CET 2023
Commit: 5015c8219b69e7f632ec966634b2b8e7749e6c90
Author: Weizhen Huang
Date: Thu Jan 5 17:46:59 2023 +0100
Branches: microfacet_hair
https://developer.blender.org/rB5015c8219b69e7f632ec966634b2b8e7749e6c90
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}
===================================================================
M source/blender/blenkernel/BKE_curves.hh
M source/blender/blenkernel/intern/curve_catmull_rom.cc
M source/blender/blenkernel/intern/curves_geometry.cc
===================================================================
diff --git a/source/blender/blenkernel/BKE_curves.hh b/source/blender/blenkernel/BKE_curves.hh
index 44d0dd132aa..49d0ac96a35 100644
--- a/source/blender/blenkernel/BKE_curves.hh
+++ b/source/blender/blenkernel/BKE_curves.hh
@@ -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);
/**
diff --git a/source/blender/blenkernel/intern/curve_catmull_rom.cc b/source/blender/blenkernel/intern/curve_catmull_rom.cc
index d16b82ed1d4..9bf5a5f1708 100644
--- a/source/blender/blenkernel/intern/curve_catmull_rom.cc
+++ b/source/blender/blenkernel/intern/curve_catmull_rom.cc
@@ -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](f
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list