diff options
Diffstat (limited to 'source/blender')
| -rw-r--r-- | source/blender/blenkernel/BKE_curves.hh | 297 | ||||
| -rw-r--r-- | source/blender/blenkernel/CMakeLists.txt | 3 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/curve_bezier.cc | 138 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/curve_catmull_rom.cc | 114 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/curve_nurbs.cc | 252 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/curves_geometry.cc | 367 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/curves_geometry_test.cc | 320 |
7 files changed, 1476 insertions, 15 deletions
diff --git a/source/blender/blenkernel/BKE_curves.hh b/source/blender/blenkernel/BKE_curves.hh index 875cf48ed43..1354fc9ed25 100644 --- a/source/blender/blenkernel/BKE_curves.hh +++ b/source/blender/blenkernel/BKE_curves.hh @@ -25,6 +25,34 @@ namespace blender::bke { +template<typename T, BLI_ENABLE_IF(std::is_integral_v<T>)> +constexpr IndexRange offsets_to_range(Span<T> offsets, int64_t index) +{ + BLI_assert(index >= 0); + BLI_assert(index < offsets.size()); + + const int offset = offsets[index]; + const int offset_next = offsets[index + 1]; + return {offset, offset_next - offset}; +} + +namespace curves::nurbs { + +struct BasisCache { + /** + * For each evaluated point, the weight for alls control points that influences it. + * The vector's size is the evaluated point count multiplied by the spline's order. + */ + Vector<float> weights; + /** + * For each evaluated point, an offset into the curve's control points for the start of #weights. + * In other words, the index of the first control point that influences this evaluated point. + */ + Vector<int> start_indices; +}; + +} // namespace curves::nurbs + /** * Contains derived data, caches, and other information not saved in files, besides a few pointers * to arrays that are kept in the non-runtime struct to avoid dereferencing this whenever they are @@ -32,6 +60,19 @@ namespace blender::bke { */ class CurvesGeometryRuntime { public: + /** + * Cache of offsets into the evaluated array for each curve, accounting for all previous + * evaluated points, Bezier curve vector segments, different resolutions per spline, etc. + */ + mutable Vector<int> evaluated_offsets_cache; + mutable Vector<int> bezier_evaluated_offsets; + mutable std::mutex offsets_cache_mutex; + mutable bool offsets_cache_dirty = true; + + mutable Vector<curves::nurbs::BasisCache> nurbs_basis_cache; + mutable std::mutex nurbs_basis_cache_mutex; + mutable bool nurbs_basis_cache_dirty = true; + /** Cache of evaluated positions. */ mutable Vector<float3> evaluated_position_cache; mutable std::mutex position_cache_mutex; @@ -88,10 +129,11 @@ class CurvesGeometry : public ::CurvesGeometry { IndexRange curves_range() const; /** - * The total number of points in the evaluated poly curve. - * This can depend on the resolution attribute if it exists. + * The index of the first point in every curve. The size of this span is one larger than the + * number of curves. Consider using #range_for_curve rather than using the offsets directly. */ - int evaluated_points_size() const; + Span<int> offsets() const; + MutableSpan<int> offsets(); /** * Access a range of indices of point data for a specific curve. @@ -104,9 +146,63 @@ class CurvesGeometry : public ::CurvesGeometry { /** Mutable access to curve types. Call #tag_topology_changed after changing any type. */ MutableSpan<int8_t> curve_types(); + bool has_curve_with_type(const CurveType type) const; + MutableSpan<float3> positions(); Span<float3> positions() const; + /** Whether the curve loops around to connect to itself, on the curve domain. */ + VArray<bool> cyclic() const; + /** Mutable access to curve cyclic values. Call #tag_topology_changed after changes. */ + MutableSpan<bool> cyclic(); + + /** + * How many evaluated points to create for each segment when evaluating Bezier, + * Catmull Rom, and NURBS curves. On the curve domain. + */ + VArray<int> resolution() const; + /** Mutable access to curve resolution. Call #tag_topology_changed after changes. */ + MutableSpan<int> resolution(); + + /** + * Handle types for Bezier control points. Call #tag_topology_changed after changes. + */ + VArray<int8_t> handle_types_left() const; + MutableSpan<int8_t> handle_types_left(); + VArray<int8_t> handle_types_right() const; + MutableSpan<int8_t> handle_types_right(); + + /** + * The positions of Bezier curve handles. Though these are really control points for the Bezier + * segments, they are stored in separate arrays to better reflect user expectations. Note that + * values may be generated automatically based on the handle types. Call #tag_positions_changed + * after changes. + */ + Span<float3> handle_positions_left() const; + MutableSpan<float3> handle_positions_left(); + Span<float3> handle_positions_right() const; + MutableSpan<float3> handle_positions_right(); + + /** + * The order (degree plus one) of each NURBS curve, on the curve domain. + * Call #tag_topology_changed after changes. + */ + VArray<int8_t> nurbs_orders() const; + MutableSpan<int8_t> nurbs_orders(); + + /** + * The automatic generation mode for each NURBS curve's knots vector, on the curve domain. + * Call #tag_topology_changed after changes. + */ + VArray<int8_t> nurbs_knots_modes() const; + MutableSpan<int8_t> nurbs_knots_modes(); + + /** + * The weight for each control point for NURBS curves. Call #tag_positions_changed after changes. + */ + Span<float> nurbs_weights() const; + MutableSpan<float> nurbs_weights(); + /** * Calculate the largest and smallest position values, only including control points * (rather than evaluated points). The existing values of `min` and `max` are taken into account. @@ -115,20 +211,49 @@ class CurvesGeometry : public ::CurvesGeometry { */ bool bounds_min_max(float3 &min, float3 &max) const; + private: /** - * The index of the first point in every curve. The size of this span is one larger than the - * number of curves. Consider using #range_for_curve rather than using the offsets directly. + * All of the curve indices for curves with a specific type. */ - Span<int> offsets() const; - MutableSpan<int> offsets(); + IndexMask indices_for_curve_type(CurveType type, Vector<int64_t> &r_indices) const; - VArray<bool> cyclic() const; - MutableSpan<bool> cyclic(); + /* -------------------------------------------------------------------- + * Evaluation. + */ + + public: + /** + * The total number of points in the evaluated poly curve. + * This can depend on the resolution attribute if it exists. + */ + int evaluated_points_size() const; + + /** + * Access a range of indices of point data for a specific curve. + * Call #evaluated_offsets() first to ensure that the evaluated offsets cache is current. + */ + IndexRange evaluated_range_for_curve(int index) const; + + /** + * The index of the first evaluated point for every curve. The size of this span is one larger + * than the number of curves. Consider using #evaluated_range_for_curve rather than using the + * offsets directly. + */ + Span<int> evaluated_offsets() const; + + Span<float3> evaluated_positions() const; + + private: + /** + * Make sure the basis weights for NURBS curve's evaluated points are calculated. + */ + void ensure_nurbs_basis_cache() const; /* -------------------------------------------------------------------- * Operations. */ + public: /** * Change the number of elements. New values for existing attributes should be properly * initialized afterwards. @@ -161,6 +286,160 @@ class CurvesGeometry : public ::CurvesGeometry { AttributeDomain to) const; }; +namespace curves { + +/** + * The number of segments between control points, accounting for the last segment of cyclic curves, + * and the fact that curves with two points cannot be cyclic. The logic is simple, but this + * function should be used to make intentions clearer. + */ +inline int curve_segment_size(const int size, const bool cyclic) +{ + return (cyclic && size > 2) ? size : size - 1; +} + +namespace bezier { + +/** + * Return true if the handles that make up a segment both have a vector type. Vector segments for + * Bezier curves have special behavior because they aren't divided into many evaluated points. + */ +bool segment_is_vector(Span<int8_t> handle_types_left, + Span<int8_t> handle_types_right, + int segment_index); + +/** + * Return true if the curve's last cylic segment has a vector type. + * This only makes a difference in the shape of cyclic curves. + */ +bool last_cylic_segment_is_vector(Span<int8_t> handle_types_left, Span<int8_t> handle_types_right); + +/** + * Calculate offsets into the curve's evaluated points for each control point. While most control + * point edges generate the number of edges specified by the resolution, vector segments only + * generate one edge. + * + * The size of the offsets array must be the same as the number of points. The value at each index + * is the evaluated point offset including the following segment. + */ +void calculate_evaluated_offsets(Span<int8_t> handle_types_left, + Span<int8_t> handle_types_right, + bool cyclic, + int resolution, + MutableSpan<int> evaluated_offsets); + +/** + * Evaluate a cubic Bezier segment, using the "forward differencing" method. + * A generic Bezier curve is made up by four points, but in many cases the first and last points + * are referred to as the control points, and the middle points are the corresponding handles. + */ +void evaluate_segment(const float3 &point_0, + const float3 &point_1, + const float3 &point_2, + const float3 &point_3, + MutableSpan<float3> result); + +/** + * Calculate all evaluated points for the Bezier curve. + * + * \param evaluated_offsets: The index in the evaluated points array for each control point, + * including the points from the corresponding segment. Used to varry the number of evaluated + * points per segment, i.e. to make vector segment only have one edge. This is expected to be + * calculated by #calculate_evaluated_offsets, and is the reason why this function doesn't need + * arguments like "cyclic" and "resolution". + */ +void calculate_evaluated_positions(Span<float3> positions, + Span<float3> handles_left, + Span<float3> handles_right, + Span<int> evaluated_offsets, + MutableSpan<float3> evaluated_positions); + +} // namespace bezier + +namespace catmull_rom { + +/** + * Calculate the number of evaluated points that #interpolate_to_evaluated is expected to produce. + * \param size: The number of points in the curve. + * \param resolution: The resolution for each segment. + */ +int calculate_evaluated_size(int size, bool cyclic, int resolution); + +/** + * Evaluate the Catmull Rom curve. The length of the #dst span should be calculated with + * #calculate_evaluated_size and is expected to divide evenly by the #src span's segment size. + */ +void interpolate_to_evaluated(fn::GSpan src, bool cyclic, int resolution, fn::GMutableSpan dst); + +} // namespace catmull_rom + +namespace nurbs { + +/** + * Checks the conditions that a NURBS curve needs to evaluate. + */ +bool check_valid_size_and_order(int size, int8_t order, bool cyclic, KnotsMode knots_mode); + +/** + * Calculate the standard evaluated size for a NURBS curve, using the standard that + * the resolution is multiplied by the number of segments between the control points. + * + * \note Though the number of evaluated points is rather arbitrary, it's useful to have a standard + * for predictability and so that cached basis weights of NURBS curves with these properties can be + * shared. + */ +int calculate_evaluated_size( + int size, int8_t order, bool cyclic, int resolution, KnotsMode knots_mode); + +/** + * Calculate the length of the knot vector for a NURBS curve with the given properties. + * The knots must be longer for a cyclic curve, for example, in order to provide weights for the + * last evaluated points that are also influenced by the first control points. + */ +int knots_size(int size, int8_t order, bool cyclic); + +/** + * Calculate the knots for a spline given its properties, based on built-in standards defined by + * #KnotsMode. + * + * \note Theoretically any sorted values can be used for NURBS knots, but calculating based + * on standard modes allows useful presets, automatic recalculation when the number of points + * changes, and is generally more intuitive than defining the knot vector manually. + */ +void calculate_knots( + int size, KnotsMode mode, int8_t order, bool cyclic, MutableSpan<float> knots); + +/** + * Based on the knots, the order, and other properties of a NURBS curve, calculate a cache that can + * be used to more simply interpolate attributes to the evaluated points later. The cache includes + * two pieces of information for every evaluated point: the first control point that influences it, + * and a weight for each control point. + */ +void calculate_basis_cache(int size, + int evaluated_size, + int8_t order, + bool cyclic, + Span<float> knots, + BasisCache &basis_cache); + +/** + * Using a "basis cache" generated by #BasisCache, interpolate attribute values to the evaluated + * points. The number of evaluated points is determined by the #basis_cache argument. + * + * \param control_weights: An optional span of control point weights, which must have the same size + * as the number of control points in the curve if provided. Using this argument gives a NURBS + * curve the "Rational" behavior that's part of its acryonym; otherwise it is a NUBS. + */ +void interpolate_to_evaluated(const BasisCache &basis_cache, + int8_t order, + Span<float> control_weights, + fn::GSpan src, + fn::GMutableSpan dst); + +} // namespace nurbs + +} // namespace curves + Curves *curves_new_nomain(int point_size, int curves_size); /** diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt index 2e32652647c..2aedd26616e 100644 --- a/source/blender/blenkernel/CMakeLists.txt +++ b/source/blender/blenkernel/CMakeLists.txt @@ -107,10 +107,13 @@ set(SRC intern/curves.cc intern/curves_geometry.cc intern/curve_bevel.c + intern/curve_bezier.cc + intern/curve_catmull_rom.cc intern/curve_convert.c intern/curve_decimate.c intern/curve_deform.c intern/curve_eval.cc + intern/curve_nurbs.cc intern/curve_to_mesh_convert.cc intern/curveprofile.cc intern/customdata.cc diff --git a/source/blender/blenkernel/intern/curve_bezier.cc b/source/blender/blenkernel/intern/curve_bezier.cc new file mode 100644 index 00000000000..84009ad51c0 --- /dev/null +++ b/source/blender/blenkernel/intern/curve_bezier.cc @@ -0,0 +1,138 @@ +/* SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup bke + */ + +#include "BKE_attribute_math.hh" +#include "BKE_curves.hh" + +namespace blender::bke::curves::bezier { + +bool segment_is_vector(const Span<int8_t> handle_types_left, + const Span<int8_t> handle_types_right, + const int segment_index) +{ + BLI_assert(handle_types_left.index_range().drop_back(1).contains(segment_index)); + return handle_types_right[segment_index] == BEZIER_HANDLE_VECTOR && + handle_types_left[segment_index + 1] == BEZIER_HANDLE_VECTOR; +} + +bool last_cylic_segment_is_vector(const Span<int8_t> handle_types_left, + const Span<int8_t> handle_types_right) +{ + return handle_types_right.last() == BEZIER_HANDLE_VECTOR && + handle_types_left.first() == BEZIER_HANDLE_VECTOR; +} + +void calculate_evaluated_offsets(const Span<int8_t> handle_types_left, + const Span<int8_t> handle_types_right, + const bool cyclic, + const int resolution, + MutableSpan<int> evaluated_offsets) +{ + const int size = handle_types_left.size(); + BLI_assert(evaluated_offsets.size() == size); + + if (size == 1) { + evaluated_offsets.first() = 1; + return; + } + + int offset = 0; + + for (const int i : IndexRange(size - 1)) { + offset += segment_is_vector(handle_types_left, handle_types_right, i) ? 1 : resolution; + evaluated_offsets[i] = offset; + } + + if (cyclic) { + offset += last_cylic_segment_is_vector(handle_types_left, handle_types_right) ? 1 : resolution; + evaluated_offsets.last(1) = offset; + } + else { + offset++; + } + + evaluated_offsets.last() = offset; +} + +void evaluate_segment(const float3 &point_0, + const float3 &point_1, + const float3 &point_2, + const float3 &point_3, + MutableSpan<float3> result) +{ + BLI_assert(result.size() > 0); + const float inv_len = 1.0f / static_cast<float>(result.size()); + const float inv_len_squared = inv_len * inv_len; + const float inv_len_cubed = inv_len_squared * inv_len; + + const float3 rt1 = 3.0f * (point_1 - point_0) * inv_len; + const float3 rt2 = 3.0f * (point_0 - 2.0f * point_1 + point_2) * inv_len_squared; + const float3 rt3 = (point_3 - point_0 + 3.0f * (point_1 - point_2)) * inv_len_cubed; + + float3 q0 = point_0; + float3 q1 = rt1 + rt2 + rt3; + float3 q2 = 2.0f * rt2 + 6.0f * rt3; + float3 q3 = 6.0f * rt3; + for (const int i : result.index_range()) { + result[i] = q0; + q0 += q1; + q1 += q2; + q2 += q3; + } +} + +void calculate_evaluated_positions(const Span<float3> positions, + const Span<float3> handles_left, + const Span<float3> handles_right, + const Span<int> evaluated_offsets, + MutableSpan<float3> evaluated_positions) +{ + BLI_assert(evaluated_offsets.last() == evaluated_positions.size()); + BLI_assert(evaluated_offsets.size() == positions.size()); + + /* Evaluate the first segment. */ + evaluate_segment(positions.first(), + handles_right.first(), + handles_left[1], + positions[1], + evaluated_positions.take_front(evaluated_offsets.first())); + + /* Give each task fewer segments as the resolution gets larger. */ + const int grain_size = std::max(evaluated_positions.size() / positions.size() * 32, 1L); + threading::parallel_for( + positions.index_range().drop_back(1).drop_front(1), grain_size, [&](IndexRange range) { + for (const int i : range) { + const IndexRange evaluated_range = offsets_to_range(evaluated_offsets, i - 1); + if (evaluated_range.size() == 1) { + evaluated_positions[evaluated_range.first()] = positions[i]; + } + else { + evaluate_segment(positions[i], + handles_right[i], + handles_left[i + 1], + positions[i + 1], + evaluated_positions.slice(evaluated_range)); + } + } + }); + + /* Evaluate the final cyclic segment if necessary. */ + const IndexRange evaluated_range = offsets_to_range(evaluated_offsets, positions.size() - 2); + if (evaluated_range.size() == 1) { + evaluated_positions.last() = positions.last(); + } + else { + evaluate_segment(positions.last(), + handles_right.last(), + handles_left.first(), + positions.first(), + evaluated_positions.slice(evaluated_range)); + } +} + +/** \} */ + +} // namespace blender::bke::curves::bezier diff --git a/source/blender/blenkernel/intern/curve_catmull_rom.cc b/source/blender/blenkernel/intern/curve_catmull_rom.cc new file mode 100644 index 00000000000..a247c67066b --- /dev/null +++ b/source/blender/blenkernel/intern/curve_catmull_rom.cc @@ -0,0 +1,114 @@ +/* SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup bke + */ + +#include "BKE_attribute_math.hh" +#include "BKE_curves.hh" + +namespace blender::bke::curves::catmull_rom { + +int calculate_evaluated_size(const int size, const bool cyclic, const int resolution) +{ + const int eval_size = resolution * curve_segment_size(size, cyclic); + /* If the curve isn't cyclic, one last point is added to the final point. */ + return (cyclic && size > 2) ? eval_size : eval_size + 1; +} + +/* Adapted from Cycles #catmull_rom_basis_eval function. */ +template<typename T> +static T calculate_basis(const T &a, const T &b, const T &c, const T &d, const float parameter) +{ + const float t = parameter; + const float s = 1.0f - parameter; + const float n0 = -t * s * s; + const float n1 = 2.0f + t * t * (3.0f * t - 5.0f); + const float n2 = 2.0f + s * s * (3.0f * s - 5.0f); + const float n3 = -s * t * t; + return 0.5f * (a * n0 + b * n1 + c * n2 + d * n3); +} + +template<typename T> +static void evaluate_segment(const T &a, const T &b, const T &c, const T &d, MutableSpan<T> dst) +{ + const float step = 1.0f / dst.size(); + dst.first() = b; + for (const int i : dst.index_range().drop_front(1)) { + dst[i] = calculate_basis<T>(a, b, c, d, i * step); + } +} + +template<typename T> +static void interpolate_to_evaluated(const Span<T> src, + const bool cyclic, + const int resolution, + MutableSpan<T> dst) + +{ + BLI_assert(dst.size() == calculate_evaluated_size(src.size(), cyclic, resolution)); + + /* - First deal with one and two point curves need special attention. + * - Then evaluate the first and last segment(s) whose control points need to wrap around + * to the other side of the source array. + * - Finally evaluate all of the segments in the middle in parallel. */ + + if (src.size() == 1) { + dst.first() = src.first(); + return; + } + if (src.size() == 2) { + evaluate_segment(src.first(), src.first(), src.last(), src.last(), dst); + dst.last() = src.last(); + return; + } + + if (cyclic) { + /* The first segment. */ + evaluate_segment(src.last(), src[0], src[1], src[2], dst.take_front(resolution)); + /* The second-to-last segment. */ + evaluate_segment(src.last(2), + src.last(1), + src.last(), + src.first(), + dst.take_back(resolution * 2).drop_back(resolution)); + /* The last segment. */ + evaluate_segment(src.last(1), src.last(), src[0], src[1], dst.take_back(resolution)); + } + else { + /* The first segment. */ + evaluate_segment(src[0], src[0], src[1], src[2], dst.take_front(resolution)); + /* The last segment. */ + evaluate_segment( + src.last(2), src.last(1), src.last(), src.last(), dst.drop_back(1).take_back(resolution)); + /* The final point of the last segment. */ + dst.last() = src.last(); + } + + /* Evaluate every segment that isn't the first or last. */ + const int grain_size = std::max(512 / resolution, 1); + const IndexRange inner_range = src.index_range().drop_back(2).drop_front(1); + threading::parallel_for(inner_range, grain_size, [&](IndexRange range) { + for (const int i : range) { + const IndexRange segment_range(resolution * i, resolution); + evaluate_segment(src[i - 1], src[i], src[i + 1], src[i + 2], dst.slice(segment_range)); + } + }); +} + +void interpolate_to_evaluated(const fn::GSpan src, + const bool cyclic, + const int resolution, + fn::GMutableSpan dst) +{ + attribute_math::convert_to_static_type(src.type(), [&](auto dummy) { + using T = decltype(dummy); + /* TODO: Use DefaultMixer or other generic mixing in the basis evaluation function to simplify + * supporting more types. */ + if constexpr (is_same_any_v<T, float, float2, float3, float4, int8_t, int, int64_t>) { + interpolate_to_evaluated(src.typed<T>(), cyclic, resolution, dst.typed<T>()); + } + }); +} + +} // namespace blender::bke::curves::catmull_rom diff --git a/source/blender/blenkernel/intern/curve_nurbs.cc b/source/blender/blenkernel/intern/curve_nurbs.cc new file mode 100644 index 00000000000..3bec01520e8 --- /dev/null +++ b/source/blender/blenkernel/intern/curve_nurbs.cc @@ -0,0 +1,252 @@ +/* SPDX-License-Identifier: GPL-2.0-or-later */ + +/** \file + * \ingroup bke + */ + +#include "BKE_attribute_math.hh" + +#include "BKE_curves.hh" + +namespace blender::bke::curves::nurbs { + +bool check_valid_size_and_order(const int size, + const int8_t order, + const bool cyclic, + const KnotsMode knots_mode) +{ + if (size < order) { + return false; + } + + if (ELEM(knots_mode, NURBS_KNOT_MODE_BEZIER, NURBS_KNOT_MODE_ENDPOINT_BEZIER)) { + if (knots_mode == NURBS_KNOT_MODE_BEZIER && size <= order) { + return false; + } + return (!cyclic || size % (order - 1) == 0); + } + + return true; +} + +int calculate_evaluated_size(const int size, + const int8_t order, + const bool cyclic, + const int resolution, + const KnotsMode knots_mode) +{ + if (!check_valid_size_and_order(size, order, cyclic, knots_mode)) { + return 0; + } + return resolution * curve_segment_size(size, cyclic); +} + +int knots_size(const int size, const int8_t order, const bool cyclic) +{ + if (cyclic) { + return size + order * 2 - 1; + } + return size + order; +} + +void calculate_knots(const int size, + const KnotsMode mode, + const int8_t order, + const bool cyclic, + MutableSpan<float> knots) +{ + BLI_assert(knots.size() == knots_size(size, order, cyclic)); + UNUSED_VARS_NDEBUG(size); + + const bool is_bezier = ELEM(mode, NURBS_KNOT_MODE_BEZIER, NURBS_KNOT_MODE_ENDPOINT_BEZIER); + const bool is_end_point = ELEM(mode, NURBS_KNOT_MODE_ENDPOINT, NURBS_KNOT_MODE_ENDPOINT_BEZIER); + /* Inner knots are always repeated once except on Bezier case. */ + const int repeat_inner = is_bezier ? order - 1 : 1; + /* How many times to repeat 0.0 at the beginning of knot. */ + const int head = is_end_point ? (order - (cyclic ? 1 : 0)) : + (is_bezier ? min_ii(2, repeat_inner) : 1); + /* Number of knots replicating widths of the starting knots. + * Covers both Cyclic and EndPoint cases. */ + const int tail = cyclic ? 2 * order - 1 : (is_end_point ? order : 0); + + int r = head; + float current = 0.0f; + + const int offset = is_end_point && cyclic ? 1 : 0; + if (offset) { + knots[0] = current; + current += 1.0f; + } + + for (const int i : IndexRange(offset, knots.size() - offset - tail)) { + knots[i] = current; + r--; + if (r == 0) { + current += 1.0; + r = repeat_inner; + } + } + + const int tail_index = knots.size() - tail; + for (const int i : IndexRange(tail)) { + knots[tail_index + i] = current + (knots[i] - knots[0]); + } +} + +static void calculate_basis_for_point(const float parameter, + const int size, + const int degree, + const Span<float> knots, + MutableSpan<float> r_weights, + int &r_start_index) +{ + const int order = degree + 1; + + int start = 0; + int end = 0; + for (const int i : IndexRange(size + degree)) { + const bool knots_equal = knots[i] == knots[i + 1]; + if (knots_equal || parameter < knots[i] || parameter > knots[i + 1]) { + continue; + } + + start = std::max(i - degree, 0); + end = i; + break; + } + + Array<float, 12> buffer(order * 2, 0.0f); + + buffer[end - start] = 1.0f; + + for (const int i_order : IndexRange(2, degree)) { + if (end + i_order >= knots.size()) { + end = size + degree - i_order; + } + for (const int i : IndexRange(end - start + 1)) { + const int knot_index = start + i; + + float new_basis = 0.0f; + if (buffer[i] != 0.0f) { + new_basis += ((parameter - knots[knot_index]) * buffer[i]) / + (knots[knot_index + i_order - 1] - knots[knot_index]); + } + + if (buffer[i + 1] != 0.0f) { + new_basis += ((knots[knot_index + i_order] - parameter) * buffer[i + 1]) / + (knots[knot_index + i_order] - knots[knot_index + 1]); + } + + buffer[i] = new_basis; + } + } + + buffer.as_mutable_span().drop_front(end - start + 1).fill(0.0f); + r_weights.copy_from(buffer.as_span().take_front(order)); + r_start_index = start; +} + +void calculate_basis_cache(const int size, + const int evaluated_size, + const int8_t order, + const bool cyclic, + const Span<float> knots, + BasisCache &basis_cache) +{ + BLI_assert(size > 0); + BLI_assert(evaluated_size > 0); + + const int8_t degree = order - 1; + + basis_cache.weights.resize(evaluated_size * order); + basis_cache.start_indices.resize(evaluated_size); + + if (evaluated_size == 0) { + return; + } + + MutableSpan<float> basis_weights(basis_cache.weights); + MutableSpan<int> basis_start_indices(basis_cache.start_indices); + + const int last_control_point_index = cyclic ? size + degree : size; + const int evaluated_segment_size = curve_segment_size(evaluated_size, cyclic); + + const float start = knots[degree]; + const float end = knots[last_control_point_index]; + const float step = (end - start) / evaluated_segment_size; + for (const int i : IndexRange(evaluated_size)) { + /* Clamp parameter due to floating point inaccuracy. */ + const float parameter = std::clamp(start + step * i, knots[0], knots[size + degree]); + + MutableSpan<float> point_weights = basis_weights.slice(i * order, order); + + calculate_basis_for_point( + parameter, last_control_point_index, degree, knots, point_weights, basis_start_indices[i]); + } +} + +template<typename T> +static void interpolate_to_evaluated(const BasisCache &basis_cache, + const int8_t order, + const Span<T> src, + MutableSpan<T> dst) +{ + attribute_math::DefaultMixer<T> mixer{dst}; + + for (const int i : dst.index_range()) { + Span<float> point_weights = basis_cache.weights.as_span().slice(i * order, order); + + for (const int j : point_weights.index_range()) { + const int point_index = (basis_cache.start_indices[i] + j) % src.size(); + mixer.mix_in(i, src[point_index], point_weights[j]); + } + } + + mixer.finalize(); +} + +template<typename T> +static void interpolate_to_evaluated_rational(const BasisCache &basis_cache, + const int8_t order, + const Span<float> control_weights, + const Span<T> src, + MutableSpan<T> dst) +{ + attribute_math::DefaultMixer<T> mixer{dst}; + + for (const int i : dst.index_range()) { + Span<float> point_weights = basis_cache.weights.as_span().slice(i * order, order); + + for (const int j : point_weights.index_range()) { + const int point_index = (basis_cache.start_indices[i] + j) % src.size(); + const float weight = point_weights[j] * control_weights[point_index]; + mixer.mix_in(i, src[point_index], weight); + } + } + + mixer.finalize(); +} + +void interpolate_to_evaluated(const BasisCache &basis_cache, + const int8_t order, + const Span<float> control_weights, + const fn::GSpan src, + fn::GMutableSpan dst) +{ + BLI_assert(dst.size() == basis_cache.start_indices.size()); + + attribute_math::convert_to_static_type(src.type(), [&](auto dummy) { + using T = decltype(dummy); + if constexpr (!std::is_void_v<attribute_math::DefaultMixer<T>>) { + if (control_weights.is_empty()) { + interpolate_to_evaluated(basis_cache, order, src.typed<T>(), dst.typed<T>()); + } + else { + interpolate_to_evaluated_rational( + basis_cache, order, control_weights, src.typed<T>(), dst.typed<T>()); + } + } + }); +} + +} // namespace blender::bke::curves::nurbs diff --git a/source/blender/blenkernel/intern/curves_geometry.cc b/source/blender/blenkernel/intern/curves_geometry.cc index dd91e788e5a..20d6c4e2f09 100644 --- a/source/blender/blenkernel/intern/curves_geometry.cc +++ b/source/blender/blenkernel/intern/curves_geometry.cc @@ -4,9 +4,12 @@ * \ingroup bke */ +#include <mutex> + #include "MEM_guardedalloc.h" #include "BLI_bounds.hh" +#include "BLI_index_mask_ops.hh" #include "DNA_curves_types.h" @@ -19,6 +22,14 @@ static const std::string ATTR_POSITION = "position"; static const std::string ATTR_RADIUS = "radius"; static const std::string ATTR_CURVE_TYPE = "curve_type"; static const std::string ATTR_CYCLIC = "cyclic"; +static const std::string ATTR_RESOLUTION = "resolution"; +static const std::string ATTR_HANDLE_TYPE_LEFT = "handle_type_left"; +static const std::string ATTR_HANDLE_TYPE_RIGHT = "handle_type_right"; +static const std::string ATTR_HANDLE_POSITION_LEFT = "handle_left"; +static const std::string ATTR_HANDLE_POSITION_RIGHT = "handle_right"; +static const std::string ATTR_NURBS_ORDER = "nurbs_order"; +static const std::string ATTR_NURBS_WEIGHT = "nurbs_weight"; +static const std::string ATTR_NURBS_KNOTS_MODE = "knots_mode"; /* -------------------------------------------------------------------- */ /** \name Constructors/Destructor @@ -152,12 +163,6 @@ IndexRange CurvesGeometry::curves_range() const return IndexRange(this->curves_size()); } -int CurvesGeometry::evaluated_points_size() const -{ - /* TODO: Implement when there are evaluated points. */ - return 0; -} - IndexRange CurvesGeometry::range_for_curve(const int index) const { BLI_assert(this->curve_size > 0); @@ -210,6 +215,22 @@ static VArray<T> get_varray_attribute(const CurvesGeometry &curves, } template<typename T> +static Span<T> get_span_attribute(const CurvesGeometry &curves, + const AttributeDomain domain, + const StringRefNull name) +{ + const int size = domain_size(curves, domain); + const CustomData &custom_data = domain_custom_data(curves, domain); + const CustomDataType type = cpp_type_to_custom_data_type(CPPType::get<T>()); + + T *data = (T *)CustomData_get_layer_named(&custom_data, type, name.c_str()); + if (data == nullptr) { + return {}; + } + return {data, size}; +} + +template<typename T> static MutableSpan<T> get_mutable_attribute(CurvesGeometry &curves, const AttributeDomain domain, const StringRefNull name) @@ -239,6 +260,20 @@ MutableSpan<int8_t> CurvesGeometry::curve_types() return get_mutable_attribute<int8_t>(*this, ATTR_DOMAIN_CURVE, ATTR_CURVE_TYPE); } +bool CurvesGeometry::has_curve_with_type(const CurveType type) const +{ + const VArray<int8_t> curve_types = this->curve_types(); + if (curve_types.is_single()) { + return curve_types.get_internal_single() == type; + } + if (curve_types.is_span()) { + return curve_types.get_internal_span().contains(type); + } + /* The curves types array should be a single value or a span. */ + BLI_assert_unreachable(); + return false; +} + MutableSpan<float3> CurvesGeometry::positions() { this->position = (float(*)[3])CustomData_duplicate_referenced_layer_named( @@ -269,6 +304,324 @@ MutableSpan<bool> CurvesGeometry::cyclic() return get_mutable_attribute<bool>(*this, ATTR_DOMAIN_CURVE, ATTR_CYCLIC); } +VArray<int> CurvesGeometry::resolution() const +{ + return get_varray_attribute<int>(*this, ATTR_DOMAIN_CURVE, ATTR_RESOLUTION, 12); +} + +MutableSpan<int> CurvesGeometry::resolution() +{ + return get_mutable_attribute<int>(*this, ATTR_DOMAIN_CURVE, ATTR_RESOLUTION); +} + +VArray<int8_t> CurvesGeometry::handle_types_left() const +{ + return get_varray_attribute<int8_t>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_TYPE_LEFT, 0); +} +MutableSpan<int8_t> CurvesGeometry::handle_types_left() +{ + return get_mutable_attribute<int8_t>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_TYPE_LEFT); +} + +VArray<int8_t> CurvesGeometry::handle_types_right() const +{ + return get_varray_attribute<int8_t>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_TYPE_RIGHT, 0); +} +MutableSpan<int8_t> CurvesGeometry::handle_types_right() +{ + return get_mutable_attribute<int8_t>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_TYPE_RIGHT); +} + +Span<float3> CurvesGeometry::handle_positions_left() const +{ + return get_span_attribute<float3>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_POSITION_LEFT); +} +MutableSpan<float3> CurvesGeometry::handle_positions_left() +{ + return get_mutable_attribute<float3>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_POSITION_LEFT); +} + +Span<float3> CurvesGeometry::handle_positions_right() const +{ + return get_span_attribute<float3>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_POSITION_RIGHT); +} +MutableSpan<float3> CurvesGeometry::handle_positions_right() +{ + return get_mutable_attribute<float3>(*this, ATTR_DOMAIN_POINT, ATTR_HANDLE_POSITION_RIGHT); +} + +VArray<int8_t> CurvesGeometry::nurbs_orders() const +{ + return get_varray_attribute<int8_t>(*this, ATTR_DOMAIN_CURVE, ATTR_NURBS_ORDER, 4); +} +MutableSpan<int8_t> CurvesGeometry::nurbs_orders() +{ + return get_mutable_attribute<int8_t>(*this, ATTR_DOMAIN_CURVE, ATTR_NURBS_ORDER); +} + +Span<float> CurvesGeometry::nurbs_weights() const +{ + return get_span_attribute<float>(*this, ATTR_DOMAIN_POINT, ATTR_NURBS_WEIGHT); +} +MutableSpan<float> CurvesGeometry::nurbs_weights() +{ + return get_mutable_attribute<float>(*this, ATTR_DOMAIN_POINT, ATTR_NURBS_WEIGHT); +} + +VArray<int8_t> CurvesGeometry::nurbs_knots_modes() const +{ + return get_varray_attribute<int8_t>(*this, ATTR_DOMAIN_CURVE, ATTR_NURBS_KNOTS_MODE, 0); +} +MutableSpan<int8_t> CurvesGeometry::nurbs_knots_modes() +{ + return get_mutable_attribute<int8_t>(*this, ATTR_DOMAIN_CURVE, ATTR_NURBS_KNOTS_MODE); +} + +/** \} */ + +/* -------------------------------------------------------------------- */ +/** \name Evaluation + * \{ */ + +template<typename SizeFn> void build_offsets(MutableSpan<int> offsets, const SizeFn &size_fn) +{ + int offset = 0; + for (const int i : offsets.drop_back(1).index_range()) { + offsets[i] = offset; + offset += size_fn(i); + } + offsets.last() = offset; +} + +static void calculate_evaluated_offsets(const CurvesGeometry &curves, + MutableSpan<int> offsets, + MutableSpan<int> bezier_evaluated_offsets) +{ + VArray<int8_t> types = curves.curve_types(); + VArray<int> resolution = curves.resolution(); + VArray<bool> cyclic = curves.cyclic(); + + VArray_Span<int8_t> handle_types_left{curves.handle_types_left()}; + VArray_Span<int8_t> handle_types_right{curves.handle_types_right()}; + + VArray<int8_t> nurbs_orders = curves.nurbs_orders(); + VArray<int8_t> nurbs_knots_modes = curves.nurbs_knots_modes(); + + build_offsets(offsets, [&](const int curve_index) -> int { + const IndexRange points = curves.range_for_curve(curve_index); + switch (types[curve_index]) { + case CURVE_TYPE_CATMULL_ROM: + return curves::catmull_rom::calculate_evaluated_size( + points.size(), cyclic[curve_index], resolution[curve_index]); + case CURVE_TYPE_POLY: + return points.size(); + case CURVE_TYPE_BEZIER: + curves::bezier::calculate_evaluated_offsets(handle_types_left.slice(points), + handle_types_right.slice(points), + cyclic[curve_index], + resolution[curve_index], + bezier_evaluated_offsets.slice(points)); + return bezier_evaluated_offsets[points.last()]; + case CURVE_TYPE_NURBS: + return curves::nurbs::calculate_evaluated_size(points.size(), + nurbs_orders[curve_index], + cyclic[curve_index], + resolution[curve_index], + KnotsMode(nurbs_knots_modes[curve_index])); + } + BLI_assert_unreachable(); + return 0; + }); +} + +int CurvesGeometry::evaluated_points_size() const +{ + /* This could avoid calculating offsets in the future in simple circumstances. */ + return this->evaluated_offsets().last(); +} + +IndexRange CurvesGeometry::evaluated_range_for_curve(int index) const +{ + BLI_assert(!this->runtime->offsets_cache_dirty); + return offsets_to_range(this->runtime->evaluated_offsets_cache.as_span(), index); +} + +Span<int> CurvesGeometry::evaluated_offsets() const +{ + if (!this->runtime->offsets_cache_dirty) { + return this->runtime->evaluated_offsets_cache; + } + + /* A double checked lock. */ + std::scoped_lock lock{this->runtime->offsets_cache_mutex}; + if (!this->runtime->offsets_cache_dirty) { + return this->runtime->evaluated_offsets_cache; + } + + threading::isolate_task([&]() { + this->runtime->evaluated_offsets_cache.resize(this->curves_size() + 1); + + if (this->has_curve_with_type(CURVE_TYPE_BEZIER)) { + this->runtime->bezier_evaluated_offsets.resize(this->points_size()); + } + else { + this->runtime->bezier_evaluated_offsets.clear_and_make_inline(); + } + + calculate_evaluated_offsets( + *this, this->runtime->evaluated_offsets_cache, this->runtime->bezier_evaluated_offsets); + }); + + this->runtime->offsets_cache_dirty = false; + return this->runtime->evaluated_offsets_cache; +} + +IndexMask CurvesGeometry::indices_for_curve_type(const CurveType type, + Vector<int64_t> &r_indices) const +{ + + VArray<int8_t> types = this->curve_types(); + if (types.is_single()) { + if (types.get_internal_single() == type) { + return IndexMask(types.size()); + } + return {}; + } + Span<int8_t> types_span = types.get_internal_span(); + return index_mask_ops::find_indices_based_on_predicate( + IndexMask(types.size()), 1024, r_indices, [&](const int index) { + return types_span[index] == type; + }); +} + +void CurvesGeometry::ensure_nurbs_basis_cache() const +{ + if (!this->runtime->nurbs_basis_cache_dirty) { + return; + } + + /* A double checked lock. */ + std::scoped_lock lock{this->runtime->nurbs_basis_cache_mutex}; + if (!this->runtime->nurbs_basis_cache_dirty) { + return; + } + + threading::isolate_task([&]() { + Vector<int64_t> nurbs_indices; + const IndexMask nurbs_mask = this->indices_for_curve_type(CURVE_TYPE_NURBS, nurbs_indices); + if (nurbs_mask.is_empty()) { + return; + } + + this->runtime->nurbs_basis_cache.resize(this->curves_size()); + MutableSpan<curves::nurbs::BasisCache> basis_caches(this->runtime->nurbs_basis_cache); + + VArray<bool> cyclic = this->cyclic(); + VArray<int8_t> orders = this->nurbs_orders(); + VArray<int8_t> knots_modes = this->nurbs_knots_modes(); + + threading::parallel_for(nurbs_mask.index_range(), 64, [&](const IndexRange range) { + for (const int curve_index : nurbs_mask.slice(range)) { + const IndexRange points = this->range_for_curve(curve_index); + const IndexRange evaluated_points = this->evaluated_range_for_curve(curve_index); + + const int8_t order = orders[curve_index]; + const bool is_cyclic = cyclic[curve_index]; + const KnotsMode mode = KnotsMode(knots_modes[curve_index]); + + const int knots_size = curves::nurbs::knots_size(points.size(), order, is_cyclic); + Array<float> knots(knots_size); + curves::nurbs::calculate_knots(points.size(), mode, order, is_cyclic, knots); + curves::nurbs::calculate_basis_cache(points.size(), + evaluated_points.size(), + order, + is_cyclic, + knots, + basis_caches[curve_index]); + } + }); + }); +} + +Span<float3> CurvesGeometry::evaluated_positions() const +{ + if (!this->runtime->position_cache_dirty) { + return this->runtime->evaluated_position_cache; + } + + /* A double checked lock. */ + std::scoped_lock lock{this->runtime->position_cache_mutex}; + if (!this->runtime->position_cache_dirty) { + return this->runtime->evaluated_position_cache; + } + + threading::isolate_task([&]() { + this->runtime->evaluated_position_cache.resize(this->evaluated_points_size()); + MutableSpan<float3> evaluated_positions = this->runtime->evaluated_position_cache; + + VArray<int8_t> types = this->curve_types(); + VArray<bool> cyclic = this->cyclic(); + VArray<int> resolution = this->resolution(); + Span<float3> positions = this->positions(); + + Span<float3> handle_positions_left = this->handle_positions_left(); + Span<float3> handle_positions_right = this->handle_positions_right(); + Span<int> bezier_evaluated_offsets = this->runtime->bezier_evaluated_offsets; + + VArray<int8_t> nurbs_orders = this->nurbs_orders(); + Span<float> nurbs_weights = this->nurbs_weights(); + + this->ensure_nurbs_basis_cache(); + + threading::parallel_for(this->curves_range(), 128, [&](IndexRange curves_range) { + for (const int curve_index : curves_range) { + const IndexRange points = this->range_for_curve(curve_index); + const IndexRange evaluated_points = this->evaluated_range_for_curve(curve_index); + + switch (types[curve_index]) { + case CURVE_TYPE_CATMULL_ROM: + curves::catmull_rom::interpolate_to_evaluated( + positions.slice(points), + cyclic[curve_index], + resolution[curve_index], + evaluated_positions.slice(evaluated_points)); + break; + case CURVE_TYPE_POLY: + evaluated_positions.slice(evaluated_points).copy_from(positions.slice(points)); + break; + case CURVE_TYPE_BEZIER: + curves::bezier::calculate_evaluated_positions( + positions.slice(points), + handle_positions_left.slice(points), + handle_positions_right.slice(points), + bezier_evaluated_offsets.slice(points), + evaluated_positions.slice(evaluated_points)); + break; + case CURVE_TYPE_NURBS: { + curves::nurbs::interpolate_to_evaluated(this->runtime->nurbs_basis_cache[curve_index], + nurbs_orders[curve_index], + nurbs_weights.slice(points), + positions.slice(points), + evaluated_positions.slice(evaluated_points)); + break; + } + default: + BLI_assert_unreachable(); + break; + } + } + }); + }); + + return this->runtime->evaluated_position_cache; +} + +/** \} */ + +/* -------------------------------------------------------------------- */ +/** \name Operations + * \{ */ + void CurvesGeometry::resize(const int point_size, const int curve_size) { if (point_size != this->point_size) { @@ -295,6 +648,8 @@ void CurvesGeometry::tag_topology_changed() this->runtime->position_cache_dirty = true; this->runtime->tangent_cache_dirty = true; this->runtime->normal_cache_dirty = true; + this->runtime->offsets_cache_dirty = true; + this->runtime->nurbs_basis_cache_dirty = true; } void CurvesGeometry::tag_normals_changed() { diff --git a/source/blender/blenkernel/intern/curves_geometry_test.cc b/source/blender/blenkernel/intern/curves_geometry_test.cc index 3a43c0c8102..1b3b3f54451 100644 --- a/source/blender/blenkernel/intern/curves_geometry_test.cc +++ b/source/blender/blenkernel/intern/curves_geometry_test.cc @@ -63,4 +63,324 @@ TEST(curves_geometry, Move) EXPECT_EQ(second_other.offsets().data(), offsets_data); } +TEST(curves_geometry, CatmullRomEvaluation) +{ + CurvesGeometry curves(4, 1); + curves.curve_types().fill(CURVE_TYPE_CATMULL_ROM); + curves.resolution().fill(12); + curves.offsets().last() = 4; + curves.cyclic().fill(false); + + MutableSpan<float3> positions = curves.positions(); + positions[0] = {1, 1, 0}; + positions[1] = {0, 1, 0}; + positions[2] = {0, 0, 0}; + positions[3] = {-1, 0, 0}; + + Span<float3> evaluated_positions = curves.evaluated_positions(); + static const Array<float3> result_1{{ + {1, 1, 0}, + {0.948495, 1.00318, 0}, + {0.87963, 1.01157, 0}, + {0.796875, 1.02344, 0}, + {0.703704, 1.03704, 0}, + {0.603588, 1.05064, 0}, + {0.5, 1.0625, 0}, + {0.396412, 1.07089, 0}, + {0.296296, 1.07407, 0}, + {0.203125, 1.07031, 0}, + {0.12037, 1.05787, 0}, + {0.0515046, 1.03501, 0}, + {0, 1, 0}, + {-0.0318287, 0.948495, 0}, + {-0.0462963, 0.87963, 0}, + {-0.046875, 0.796875, 0}, + {-0.037037, 0.703704, 0}, + {-0.0202546, 0.603588, 0}, + {0, 0.5, 0}, + {0.0202546, 0.396412, 0}, + {0.037037, 0.296296, 0}, + {0.046875, 0.203125, 0}, + {0.0462963, 0.12037, 0}, + {0.0318287, 0.0515046, 0}, + {0, 0, 0}, + {-0.0515046, -0.0350116, 0}, + {-0.12037, -0.0578704, 0}, + {-0.203125, -0.0703125, 0}, + {-0.296296, -0.0740741, 0}, + {-0.396412, -0.0708912, 0}, + {-0.5, -0.0625, 0}, + {-0.603588, -0.0506366, 0}, + {-0.703704, -0.037037, 0}, + {-0.796875, -0.0234375, 0}, + {-0.87963, -0.0115741, 0}, + {-0.948495, -0.00318287, 0}, + {-1, 0, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_1[i], 1e-5f); + } + + /* Changing the positions shouldn't cause the evaluated positions array to be reallocated. */ + curves.tag_positions_changed(); + curves.evaluated_positions(); + EXPECT_EQ(curves.evaluated_positions().data(), evaluated_positions.data()); + + /* Call recalculation (which shouldn't happen because low-level accessors don't tag caches). */ + EXPECT_EQ(evaluated_positions[12].x, 0.0f); + EXPECT_EQ(evaluated_positions[12].y, 1.0f); + + positions[0] = {1, 0, 0}; + positions[1] = {1, 1, 0}; + positions[2] = {0, 1, 0}; + positions[3] = {0, 0, 0}; + curves.cyclic().fill(true); + + /* Tag topology changed because the new cyclic value is different. */ + curves.tag_topology_changed(); + + /* Retrieve the data again since the size should be larger than last time (one more segment). */ + evaluated_positions = curves.evaluated_positions(); + static const Array<float3> result_2{{ + {1, 0, 0}, + {1.03819, 0.0515046, 0}, + {1.06944, 0.12037, 0}, + {1.09375, 0.203125, 0}, + {1.11111, 0.296296, 0}, + {1.12153, 0.396412, 0}, + {1.125, 0.5, 0}, + {1.12153, 0.603588, 0}, + {1.11111, 0.703704, 0}, + {1.09375, 0.796875, 0}, + {1.06944, 0.87963, 0}, + {1.03819, 0.948495, 0}, + {1, 1, 0}, + {0.948495, 1.03819, 0}, + {0.87963, 1.06944, 0}, + {0.796875, 1.09375, 0}, + {0.703704, 1.11111, 0}, + {0.603588, 1.12153, 0}, + {0.5, 1.125, 0}, + {0.396412, 1.12153, 0}, + {0.296296, 1.11111, 0}, + {0.203125, 1.09375, 0}, + {0.12037, 1.06944, 0}, + {0.0515046, 1.03819, 0}, + {0, 1, 0}, + {-0.0381944, 0.948495, 0}, + {-0.0694444, 0.87963, 0}, + {-0.09375, 0.796875, 0}, + {-0.111111, 0.703704, 0}, + {-0.121528, 0.603588, 0}, + {-0.125, 0.5, 0}, + {-0.121528, 0.396412, 0}, + {-0.111111, 0.296296, 0}, + {-0.09375, 0.203125, 0}, + {-0.0694444, 0.12037, 0}, + {-0.0381944, 0.0515046, 0}, + {0, 0, 0}, + {0.0515046, -0.0381944, 0}, + {0.12037, -0.0694444, 0}, + {0.203125, -0.09375, 0}, + {0.296296, -0.111111, 0}, + {0.396412, -0.121528, 0}, + {0.5, -0.125, 0}, + {0.603588, -0.121528, 0}, + {0.703704, -0.111111, 0}, + {0.796875, -0.09375, 0}, + {0.87963, -0.0694444, 0}, + {0.948495, -0.0381944, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_2[i], 1e-5f); + } +} + +TEST(curves_geometry, CatmullRomTwoPointCyclic) +{ + CurvesGeometry curves(2, 1); + curves.curve_types().fill(CURVE_TYPE_CATMULL_ROM); + curves.resolution().fill(12); + curves.offsets().last() = 2; + curves.cyclic().fill(true); + + /* The cyclic value should be ignored when there are only two control points. There should + * be 12 evaluated points for the single segment and an extra for the last point. */ + EXPECT_EQ(curves.evaluated_points_size(), 13); +} + +TEST(curves_geometry, BezierPositionEvaluation) +{ + CurvesGeometry curves(2, 1); + curves.curve_types().fill(CURVE_TYPE_BEZIER); + curves.resolution().fill(12); + curves.offsets().last() = 2; + + MutableSpan<float3> handles_left = curves.handle_positions_left(); + MutableSpan<float3> handles_right = curves.handle_positions_right(); + MutableSpan<float3> positions = curves.positions(); + positions.first() = {-1, 0, 0}; + positions.last() = {1, 0, 0}; + handles_right.first() = {-0.5f, 0.5f, 0.0f}; + handles_left.last() = {0, 0, 0}; + + /* Dangling handles shouldn't be used in a non-cyclic curve. */ + handles_left.first() = {100, 100, 100}; + handles_right.last() = {100, 100, 100}; + + Span<float3> evaluated_positions = curves.evaluated_positions(); + static const Array<float3> result_1{{ + {-1, 0, 0}, + {-0.874711, 0.105035, 0}, + {-0.747685, 0.173611, 0}, + {-0.617188, 0.210937, 0}, + {-0.481481, 0.222222, 0}, + {-0.338831, 0.212674, 0}, + {-0.1875, 0.1875, 0}, + {-0.0257524, 0.15191, 0}, + {0.148148, 0.111111, 0}, + {0.335937, 0.0703125, 0}, + {0.539352, 0.0347222, 0}, + {0.760127, 0.00954859, 0}, + {1, 0, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_1[i], 1e-5f); + } + + curves.resize(4, 2); + curves.curve_types().fill(CURVE_TYPE_BEZIER); + curves.resolution().fill(9); + curves.offsets().last() = 4; + handles_left = curves.handle_positions_left(); + handles_right = curves.handle_positions_right(); + positions = curves.positions(); + positions[2] = {-1, 1, 0}; + positions[3] = {1, 1, 0}; + handles_right[2] = {-0.5f, 1.5f, 0.0f}; + handles_left[3] = {0, 1, 0}; + + /* Dangling handles shouldn't be used in a non-cyclic curve. */ + handles_left[2] = {-100, -100, -100}; + handles_right[3] = {-100, -100, -100}; + + evaluated_positions = curves.evaluated_positions(); + EXPECT_EQ(evaluated_positions.size(), 20); + static const Array<float3> result_2{{ + {-1, 0, 0}, + {-0.832647, 0.131687, 0}, + {-0.66118, 0.201646, 0}, + {-0.481481, 0.222222, 0}, + {-0.289438, 0.205761, 0}, + {-0.0809327, 0.164609, 0}, + {0.148148, 0.111111, 0}, + {0.40192, 0.0576133, 0}, + {0.684499, 0.016461, 0}, + {1, 0, 0}, + {-1, 1, 0}, + {-0.832647, 1.13169, 0}, + {-0.66118, 1.20165, 0}, + {-0.481481, 1.22222, 0}, + {-0.289438, 1.20576, 0}, + {-0.0809327, 1.16461, 0}, + {0.148148, 1.11111, 0}, + {0.40192, 1.05761, 0}, + {0.684499, 1.01646, 0}, + {1, 1, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_2[i], 1e-5f); + } +} + +TEST(curves_geometry, NURBSEvaluation) +{ + CurvesGeometry curves(4, 1); + curves.curve_types().fill(CURVE_TYPE_NURBS); + curves.resolution().fill(10); + curves.offsets().last() = 4; + + MutableSpan<float3> positions = curves.positions(); + positions[0] = {1, 1, 0}; + positions[1] = {0, 1, 0}; + positions[2] = {0, 0, 0}; + positions[3] = {-1, 0, 0}; + + Span<float3> evaluated_positions = curves.evaluated_positions(); + static const Array<float3> result_1{{ + {0.166667, 0.833333, 0}, {0.150006, 0.815511, 0}, {0.134453, 0.796582, 0}, + {0.119924, 0.776627, 0}, {0.106339, 0.75573, 0}, {0.0936146, 0.733972, 0}, + {0.0816693, 0.711434, 0}, {0.0704211, 0.6882, 0}, {0.0597879, 0.66435, 0}, + {0.0496877, 0.639968, 0}, {0.0400385, 0.615134, 0}, {0.0307584, 0.589931, 0}, + {0.0217653, 0.564442, 0}, {0.0129772, 0.538747, 0}, {0.00431208, 0.512929, 0}, + {-0.00431208, 0.487071, 0}, {-0.0129772, 0.461253, 0}, {-0.0217653, 0.435558, 0}, + {-0.0307584, 0.410069, 0}, {-0.0400385, 0.384866, 0}, {-0.0496877, 0.360032, 0}, + {-0.0597878, 0.33565, 0}, {-0.0704211, 0.3118, 0}, {-0.0816693, 0.288566, 0}, + {-0.0936146, 0.266028, 0}, {-0.106339, 0.24427, 0}, {-0.119924, 0.223373, 0}, + {-0.134453, 0.203418, 0}, {-0.150006, 0.184489, 0}, {-0.166667, 0.166667, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_1[i], 1e-5f); + } + + /* Test a cyclic curve. */ + curves.cyclic().fill(true); + curves.tag_topology_changed(); + evaluated_positions = curves.evaluated_positions(); + static const Array<float3> result_2{{ + {0.166667, 0.833333, 0}, {0.121333, 0.778667, 0}, + {0.084, 0.716, 0}, {0.0526667, 0.647333, 0}, + {0.0253333, 0.574667, 0}, {0, 0.5, 0}, + {-0.0253333, 0.425333, 0}, {-0.0526667, 0.352667, 0}, + {-0.084, 0.284, 0}, {-0.121333, 0.221333, 0}, + {-0.166667, 0.166667, 0}, {-0.221, 0.121667, 0}, + {-0.281333, 0.0866667, 0}, {-0.343667, 0.0616666, 0}, + {-0.404, 0.0466667, 0}, {-0.458333, 0.0416667, 0}, + {-0.502667, 0.0466667, 0}, {-0.533, 0.0616666, 0}, + {-0.545333, 0.0866667, 0}, {-0.535667, 0.121667, 0}, + {-0.5, 0.166667, 0}, {-0.436, 0.221334, 0}, + {-0.348, 0.284, 0}, {-0.242, 0.352667, 0}, + {-0.124, 0.425333, 0}, {0, 0.5, 0}, + {0.124, 0.574667, 0}, {0.242, 0.647333, 0}, + {0.348, 0.716, 0}, {0.436, 0.778667, 0}, + {0.5, 0.833333, 0}, {0.535667, 0.878334, 0}, + {0.545333, 0.913333, 0}, {0.533, 0.938333, 0}, + {0.502667, 0.953333, 0}, {0.458333, 0.958333, 0}, + {0.404, 0.953333, 0}, {0.343667, 0.938333, 0}, + {0.281333, 0.913333, 0}, {0.221, 0.878333, 0}, + }}; + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_2[i], 1e-5f); + } + + /* Test a circular cyclic curve with weights. */ + positions[0] = {1, 0, 0}; + positions[1] = {1, 1, 0}; + positions[2] = {0, 1, 0}; + positions[3] = {0, 0, 0}; + curves.nurbs_weights().fill(1.0f); + curves.nurbs_weights()[0] = 4.0f; + curves.tag_positions_changed(); + static const Array<float3> result_3{{ + {0.888889, 0.555556, 0}, {0.837792, 0.643703, 0}, {0.773885, 0.727176, 0}, + {0.698961, 0.800967, 0}, {0.616125, 0.860409, 0}, {0.529412, 0.901961, 0}, + {0.443152, 0.923773, 0}, {0.361289, 0.925835, 0}, {0.286853, 0.909695, 0}, + {0.221722, 0.877894, 0}, {0.166667, 0.833333, 0}, {0.122106, 0.778278, 0}, + {0.0903055, 0.713148, 0}, {0.0741654, 0.638711, 0}, {0.0762274, 0.556847, 0}, + {0.0980392, 0.470588, 0}, {0.139591, 0.383875, 0}, {0.199032, 0.301039, 0}, + {0.272824, 0.226114, 0}, {0.356297, 0.162208, 0}, {0.444444, 0.111111, 0}, + {0.531911, 0.0731388, 0}, {0.612554, 0.0468976, 0}, {0.683378, 0.0301622, 0}, + {0.74391, 0.0207962, 0}, {0.794872, 0.017094, 0}, {0.837411, 0.017839, 0}, + {0.872706, 0.0222583, 0}, {0.901798, 0.0299677, 0}, {0.925515, 0.0409445, 0}, + {0.944444, 0.0555556, 0}, {0.959056, 0.0744855, 0}, {0.970032, 0.0982019, 0}, + {0.977742, 0.127294, 0}, {0.982161, 0.162589, 0}, {0.982906, 0.205128, 0}, + {0.979204, 0.256091, 0}, {0.969838, 0.316622, 0}, {0.953102, 0.387446, 0}, + {0.926861, 0.468089, 0}, + }}; + evaluated_positions = curves.evaluated_positions(); + for (const int i : evaluated_positions.index_range()) { + EXPECT_V3_NEAR(evaluated_positions[i], result_3[i], 1e-5f); + } +} + } // namespace blender::bke::tests |