Merge branch 'pm_seam_curling_fix_250' into master_250

3 files changed, 119 insertions, 107 deletions

diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp
index 4185ae41b..2a716a300 100644
--- a/src/libslic3r/GCode/SeamPlacer.cpp
+++ b/src/libslic3r/GCode/SeamPlacer.cpp
@@ -64,6 +64,16 @@ float gauss(float value, float mean_x_coord, float mean_value, float falloff_spe
     return mean_value * (std::exp(exponent) - 1.0f) / (std::exp(1.0f) - 1.0f);
 }
 
+float compute_angle_penalty(float ccw_angle) {
+    // This function is used:
+    // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x)))
+    // looks scary, but it is gaussian combined with sigmoid,
+    // so that concave points have much smaller penalty over convex ones
+    // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png
+    return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) +
+            1.0f / (2 + std::exp(-ccw_angle));
+}
+
 /// Coordinate frame
 class Frame {
 public:
@@ -457,22 +467,41 @@ Polygons extract_perimeter_polygons(const Layer *layer, const SeamPosition confi
 //each SeamCandidate also contains pointer to shared Perimeter structure representing the polygon
 // if Custom Seam modifiers are present, oversamples the polygon if necessary to better fit user intentions
 void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const LayerRegion *region,
-        const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) {
+        bool arachne_generated, const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) {
     if (orig_polygon.size() == 0) {
         return;
     }
-
     Polygon polygon = orig_polygon;
+    bool was_clockwise = polygon.make_counter_clockwise();
+
     std::vector<float> lengths { };
     for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) {
         lengths.push_back((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm());
     }
     lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.1));
-
-    bool was_clockwise = polygon.make_counter_clockwise();
-    std::vector<float> local_angles = calculate_polygon_angles_at_vertices(polygon, lengths,
+    std::vector<float> polygon_angles = calculate_polygon_angles_at_vertices(polygon, lengths,
             SeamPlacer::polygon_local_angles_arm_distance);
 
+    // resample smooth surfaces from arachne, so that alignment finds short path down, and does not create unnecesary curves
+    if (arachne_generated && std::all_of(polygon_angles.begin(), polygon_angles.end(), [](float angle) {
+    	return compute_angle_penalty(angle) > SeamPlacer::sharp_angle_penalty_snapping_threshold;
+    })) {
+    	float total_dist = std::accumulate(lengths.begin(), lengths.end(), 0.0f);
+    	float avg_dist = total_dist / float(lengths.size());
+    	if (avg_dist < SeamPlacer::seam_align_tolerable_dist * 2.0f){
+			coord_t sampling_dist = scaled(avg_dist*0.2f);
+
+			polygon.points = polygon.equally_spaced_points(sampling_dist);
+			lengths.clear();
+			for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) {
+				lengths.push_back((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm());
+			}
+			lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.1));
+			polygon_angles = calculate_polygon_angles_at_vertices(polygon, lengths, avg_dist);
+		}
+    }
+
+
     result.perimeters.push_back( { });
     Perimeter &perimeter = result.perimeters.back();
 
@@ -498,16 +527,16 @@ void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const
         } else {
             position = orig_polygon_points.front();
             orig_polygon_points.pop();
-            local_ccw_angle = was_clockwise ? -local_angles[orig_angle_index] : local_angles[orig_angle_index];
+            local_ccw_angle = was_clockwise ? -polygon_angles[orig_angle_index] : polygon_angles[orig_angle_index];
             orig_angle_index++;
             orig_point = true;
         }
 
-        if (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) {
+        if (global_model_info.is_enforced(position, perimeter.flow_width)) {
             type = EnforcedBlockedSeamPoint::Enforced;
         }
 
-        if (global_model_info.is_blocked(position, SeamPlacer::enforcer_blocker_distance_tolerance)) {
+        if (global_model_info.is_blocked(position, perimeter.flow_width)) {
             type = EnforcedBlockedSeamPoint::Blocked;
         }
         some_point_enforced = some_point_enforced || type == EnforcedBlockedSeamPoint::Enforced;
@@ -560,7 +589,8 @@ void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const
                 break;
             }
             viable_points_indices.push_back(last_enforced_idx);
-            if (abs(result.points[last_enforced_idx].local_ccw_angle) > SeamPlacer::sharp_angle_snapping_threshold) {
+            if (compute_angle_penalty(result.points[last_enforced_idx].local_ccw_angle)
+            		< SeamPlacer::sharp_angle_penalty_snapping_threshold) {
                 orig_large_angle_points_indices.push_back(last_enforced_idx);
             }
             last_enforced_idx = next_index(last_enforced_idx);
@@ -827,17 +857,10 @@ struct SeamComparator {
         return is_first_not_much_worse(a, b) && is_first_not_much_worse(b, a);
     }
 
-    float compute_angle_penalty(float ccw_angle) const {
-        // This function is used:
-        // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x)))
-        // looks scary, but it is gaussian combined with sigmoid,
-        // so that concave points have much smaller penalty over convex ones
-        // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png
-        return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) +
-                1.0f / (2 + std::exp(-ccw_angle));
-    }
-
     float weight(const SeamCandidate &a) const {
+    	if (setup == SeamPosition::spAligned && a.central_enforcer) {
+    		return 2.0f;
+    	}
         return a.visibility + angle_importance * compute_angle_penalty(a.local_ccw_angle) / (1.0f + angle_importance);
     }
 };
@@ -1046,12 +1069,14 @@ public:
 void SeamPlacer::gather_seam_candidates(const PrintObject *po,
         const SeamPlacerImpl::GlobalModelInfo &global_model_info, const SeamPosition configured_seam_preference) {
     using namespace SeamPlacerImpl;
+    bool arachne_generated = po->config().perimeter_generator == PerimeterGeneratorType::Arachne;
 
     PrintObjectSeamData &seam_data = m_seam_per_object.emplace(po, PrintObjectSeamData { }).first->second;
     seam_data.layers.resize(po->layer_count());
 
     tbb::parallel_for(tbb::blocked_range<size_t>(0, po->layers().size()),
-            [po, configured_seam_preference, &global_model_info, &seam_data](tbb::blocked_range<size_t> r) {
+            [po, configured_seam_preference, arachne_generated, &global_model_info, &seam_data]
+			 (tbb::blocked_range<size_t> r) {
                 for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) {
                     PrintObjectSeamData::LayerSeams &layer_seams = seam_data.layers[layer_idx];
                     const Layer *layer = po->get_layer(layer_idx);
@@ -1061,7 +1086,7 @@ void SeamPlacer::gather_seam_candidates(const PrintObject *po,
                     Polygons polygons = extract_perimeter_polygons(layer, configured_seam_preference, regions);
                     for (size_t poly_index = 0; poly_index < polygons.size(); ++poly_index) {
                         process_perimeter_polygon(polygons[poly_index], unscaled_z,
-                                regions[poly_index], global_model_info, layer_seams);
+                                regions[poly_index], arachne_generated, global_model_info, layer_seams);
                     }
                     auto functor = SeamCandidateCoordinateFunctor { layer_seams.points };
                     seam_data.layers[layer_idx].points_tree =
@@ -1135,16 +1160,12 @@ void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po)
 // Used by align_seam_points().
 std::optional<std::pair<size_t, size_t>> SeamPlacer::find_next_seam_in_layer(
         const std::vector<PrintObjectSeamData::LayerSeams> &layers,
-        const std::pair<size_t, size_t> &prev_point_index,
-        const size_t layer_idx, const float slice_z,
+        const Vec3f& projected_position,
+        const size_t layer_idx, const float max_distance,
         const SeamPlacerImpl::SeamComparator &comparator) const {
     using namespace SeamPlacerImpl;
-
-    const SeamCandidate &last_point = layers[prev_point_index.first].points[prev_point_index.second];
-
-    Vec3f projected_position { last_point.position.x(), last_point.position.y(), slice_z };
     std::vector<size_t> nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position,
-            SeamPlacer::seam_align_tolerable_dist);
+            max_distance);
 
     if (nearby_points_indices.empty()) {
         return {};
@@ -1185,7 +1206,7 @@ std::optional<std::pair<size_t, size_t>> SeamPlacer::find_next_seam_in_layer(
     // First try to pick central enforcer if any present
     if (next_layer_seam.central_enforcer
             && (next_layer_seam.position - projected_position).squaredNorm()
-                    < sqr(3 * SeamPlacer::seam_align_tolerable_dist)) {
+                    < sqr(3 * max_distance)) {
         return {std::pair<size_t, size_t> {layer_idx, nearest_point.perimeter.seam_index}};
     }
 
@@ -1203,69 +1224,60 @@ std::optional<std::pair<size_t, size_t>> SeamPlacer::find_next_seam_in_layer(
 
 std::vector<std::pair<size_t, size_t>> SeamPlacer::find_seam_string(const PrintObject *po,
         std::pair<size_t, size_t> start_seam, const SeamPlacerImpl::SeamComparator &comparator,
-        std::optional<std::pair<size_t, size_t>> &out_best_moved_seam, size_t &out_moved_seams_count) const {
-    out_best_moved_seam.reset();
-    out_moved_seams_count = 0;
+        float& string_weight) const {
+    string_weight = 0.0f;
     const std::vector<PrintObjectSeamData::LayerSeams> &layers = m_seam_per_object.find(po)->second.layers;
     int layer_idx = start_seam.first;
-    int seam_index = start_seam.second;
 
     //initialize searching for seam string - cluster of nearby seams on previous and next layers
     int next_layer = layer_idx + 1;
+    int step = 1;
     std::pair<size_t, size_t> prev_point_index = start_seam;
     std::vector<std::pair<size_t, size_t>> seam_string { start_seam };
 
-    //find seams or potential seams in forward direction; there is a budget of skips allowed
-    while (next_layer < int(layers.size())) {
-        auto maybe_next_seam = find_next_seam_in_layer(layers, prev_point_index, next_layer,
-                float(po->get_layer(next_layer)->slice_z), comparator);
-        if (maybe_next_seam.has_value()) {
+    auto reverse_lookup_direction = [&]() {
+        step = -1;
+        prev_point_index = start_seam;
+        next_layer = layer_idx - 1;
+    };
 
-            // For old macOS (pre 10.14), std::optional does not have .value() method, so the code is using operator*() instead.
-            std::pair<size_t, size_t> next_seam_coords = maybe_next_seam.operator*();
-            const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second];
-            bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second;
-            out_moved_seams_count += is_moved;
-            if (is_moved && (!out_best_moved_seam.has_value() ||
-                    comparator.is_first_better(next_seam,
-                            layers[out_best_moved_seam.operator*().first].points[out_best_moved_seam.operator*().second]))) {
-                out_best_moved_seam = { next_seam_coords };
+    while (next_layer >= 0) {
+        if (next_layer >= int(layers.size())) {
+            reverse_lookup_direction();
+            if (next_layer < 0) {
+                break;
             }
-
-            seam_string.push_back(maybe_next_seam.operator*());
-            prev_point_index = seam_string.back();
-            //String added, prev_point_index updated
-        } else {
-            break;
         }
-        next_layer++;
-    }
+        float max_distance = SeamPlacer::seam_align_tolerable_dist;
+        Vec3f prev_position = layers[prev_point_index.first].points[prev_point_index.second].position;
+        Vec3f projected_position = prev_position;
+        projected_position.z() = float(po->get_layer(next_layer)->slice_z);
 
-    //do additional check in back direction
-    next_layer = layer_idx - 1;
-    prev_point_index = std::pair<size_t, size_t>(layer_idx, seam_index);
-    while (next_layer >= 0) {
-        auto maybe_next_seam = find_next_seam_in_layer(layers, prev_point_index, next_layer,
-                float(po->get_layer(next_layer)->slice_z), comparator);
-        if (maybe_next_seam.has_value()) {
+        std::optional<std::pair<size_t, size_t>> maybe_next_seam = find_next_seam_in_layer(layers, projected_position,
+                next_layer,
+                max_distance, comparator);
 
+        if (maybe_next_seam.has_value()) {
+            // For old macOS (pre 10.14), std::optional does not have .value() method, so the code is using operator*() instead.
             std::pair<size_t, size_t> next_seam_coords = maybe_next_seam.operator*();
             const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second];
             bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second;
-            out_moved_seams_count += is_moved;
-            if (is_moved && (!out_best_moved_seam.has_value() ||
-                    comparator.is_first_better(next_seam,
-                            layers[out_best_moved_seam.operator*().first].points[out_best_moved_seam.operator*().second]))) {
-                out_best_moved_seam = { next_seam_coords };
-            }
-
+            string_weight += comparator.weight(next_seam) -
+                    is_moved ? comparator.weight(layers[next_seam_coords.first].points[next_seam.perimeter.seam_index]) : 0.0f;
             seam_string.push_back(maybe_next_seam.operator*());
             prev_point_index = seam_string.back();
             //String added, prev_point_index updated
         } else {
-            break;
+            if (step == 1) {
+                reverse_lookup_direction();
+                if (next_layer < 0) {
+                     break;
+                 }
+            } else {
+                break;
+            }
         }
-        next_layer--;
+        next_layer += step;
     }
 
     return seam_string;
@@ -1311,7 +1323,7 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::
     }
 
     //sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice
-    std::sort(seams.begin(), seams.end(),
+    std::stable_sort(seams.begin(), seams.end(),
             [&comparator, &layers](const std::pair<size_t, size_t> &left,
                     const std::pair<size_t, size_t> &right) {
                 return comparator.is_first_better(layers[left.first].points[left.second],
@@ -1337,22 +1349,21 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::
             // This perimeter is already aligned, skip seam
             continue;
         } else {
-            std::optional<std::pair<size_t, size_t>> best_moved_seam;
-            size_t moved_seams_count;
-            seam_string = this->find_seam_string(po, { layer_idx, seam_index }, comparator, best_moved_seam,
-                    moved_seams_count);
-            if (best_moved_seam.has_value()) {
-                size_t alternative_moved_seams_count;
-                alternative_seam_string = this->find_seam_string(po, best_moved_seam.operator*(), comparator,
-                        best_moved_seam, alternative_moved_seams_count);
+            float seam_string_weight;
+            seam_string = this->find_seam_string(po, { layer_idx, seam_index }, comparator, seam_string_weight);
+            size_t step_size = 1 + seam_string.size() / 20;
+            for (size_t alternative_start = 0; alternative_start < seam_string.size(); alternative_start+=step_size) {
+                float alternative_seam_string_weight = 0;
+                size_t start_layer_idx = seam_string[alternative_start].first;
+                size_t seam_idx = layers[start_layer_idx].points[seam_string[alternative_start].second].perimeter.seam_index;
+                alternative_seam_string = this->find_seam_string(po, std::pair<size_t,size_t>(start_layer_idx, seam_idx), comparator,
+                        alternative_seam_string_weight);
                 if (alternative_seam_string.size() >= SeamPlacer::seam_align_minimum_string_seams &&
-                        alternative_moved_seams_count < moved_seams_count) {
+                        alternative_seam_string_weight > seam_string_weight) {
+                    seam_string_weight = alternative_seam_string_weight;
                     seam_string = std::move(alternative_seam_string);
-                    // finish loop. but repeat the alignment for the current seam, since it could be skipped due to alternative path being aligned.
-                    global_index--;
                 }
             }
-
             if (seam_string.size() < seam_align_minimum_string_seams) {
                 //string NOT long enough to be worth aligning, skip
                 continue;
@@ -1365,26 +1376,29 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::
                         return left.first < right.first;
                     });
 
+            //repeat the alignment for the current seam, since it could be skipped due to alternative path being aligned.
+            global_index--;
+
             // gather all positions of seams and their weights (weights are derived as negative penalty, they are made positive in next step)
             observations.resize(seam_string.size());
             observation_points.resize(seam_string.size());
             weights.resize(seam_string.size());
 
             //gather points positions and weights
-            // The algorithm uses only angle to compute penalty, to enforce snapping to sharp corners, if they are present
-            // after several experiments approach that gives best results is to snap the weight to one for sharp corners, and
-            //  leave it small for others. However, this can result in non-smooth line over area with a lot of unaligned sharp corners.
+            float total_length = 0.0f;
+            Vec3f last_point_pos = layers[seam_string[0].first].points[seam_string[0].second].position;
             for (size_t index = 0; index < seam_string.size(); ++index) {
                 Vec3f pos = layers[seam_string[index].first].points[seam_string[index].second].position;
+                total_length += (last_point_pos - pos).norm();
+                last_point_pos = pos;
                 observations[index] = pos.head<2>();
                 observation_points[index] = pos.z();
-                weights[index] = std::min(1.0f,
-                        comparator.weight(layers[seam_string[index].first].points[seam_string[index].second]));
+                weights[index] = comparator.weight(layers[seam_string[index].first].points[seam_string[index].second]);
             }
 
             // Curve Fitting
             size_t number_of_segments = std::max(size_t(1),
-                    size_t(observations.size() / SeamPlacer::seam_align_seams_per_segment));
+                    size_t(total_length / SeamPlacer::seam_align_mm_per_segment));
             auto curve = Geometry::fit_cubic_bspline(observations, observation_points, weights, number_of_segments);
 
             // Do alignment - compute fitted point for each point in the string from its Z coord, and store the position into
@@ -1392,14 +1406,14 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::
             for (size_t index = 0; index < seam_string.size(); ++index) {
                 const auto &pair = seam_string[index];
                 const float t =
-                        abs(layers[pair.first].points[pair.second].local_ccw_angle)
-                                > SeamPlacer::sharp_angle_snapping_threshold
-                                ? 1.0 : 0.0f;
+                        compute_angle_penalty(layers[pair.first].points[pair.second].local_ccw_angle)
+                                < SeamPlacer::sharp_angle_penalty_snapping_threshold
+                                ? 0.8f : 0.0f;
                 Vec3f current_pos = layers[pair.first].points[pair.second].position;
                 Vec2f fitted_pos = curve.get_fitted_value(current_pos.z());
 
                 //interpolate between current and fitted position, prefer current pos for large weights.
-                Vec3f final_position = t * current_pos + (1 - t) * to_3d(fitted_pos, current_pos.z());
+                Vec3f final_position = t * current_pos + (1.0f - t) * to_3d(fitted_pos, current_pos.z());
 
                 Perimeter &perimeter = layers[pair.first].points[pair.second].perimeter;
                 perimeter.seam_index = pair.second;
diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp
index cfad8c1ba..1f1643d0d 100644
--- a/src/libslic3r/GCode/SeamPlacer.hpp
+++ b/src/libslic3r/GCode/SeamPlacer.hpp
@@ -128,17 +128,16 @@ public:
 
     // arm length used during angles computation
     static constexpr float polygon_local_angles_arm_distance = 0.3f;
-    static constexpr float sharp_angle_snapping_threshold = 0.3f * float(PI);
+    // value for angles with penalty lower than this threshold - such angles will be snapped to their original position instead of spline interpolated position
+    static constexpr float sharp_angle_penalty_snapping_threshold = 0.6f;
 
     // max tolerable distance from the previous layer is overhang_distance_tolerance_factor * flow_width
     static constexpr float overhang_distance_tolerance_factor = 0.5f;
 
     // determines angle importance compared to visibility ( neutral value is 1.0f. )
     static constexpr float angle_importance_aligned = 0.6f;
-    static constexpr float angle_importance_nearest = 1.0f; // use much higher angle importance for nearest mode, to combat the visiblity info noise
+    static constexpr float angle_importance_nearest = 1.0f; // use much higher angle importance for nearest mode, to combat the visibility info noise
 
-    // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer
-    static constexpr float enforcer_blocker_distance_tolerance = 0.35f;
     // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size
     static constexpr float enforcer_oversampling_distance = 0.2f;
 
@@ -148,9 +147,9 @@ public:
     // seam_align_tolerable_dist - if next layer closest point is too far away, break aligned string
     static constexpr float seam_align_tolerable_dist = 1.0f;
     // minimum number of seams needed in cluster to make alignment happen
-    static constexpr size_t seam_align_minimum_string_seams = 10;
-    // points covered by spline; determines number of splines for the given string
-    static constexpr size_t seam_align_seams_per_segment = 16;
+    static constexpr size_t seam_align_minimum_string_seams = 6;
+    // millimeters covered by spline; determines number of splines for the given string
+    static constexpr size_t seam_align_mm_per_segment = 4.0f;
 
     //The following data structures hold all perimeter points for all PrintObject.
     std::unordered_map<const PrintObject*, PrintObjectSeamData> m_seam_per_object;
@@ -169,12 +168,11 @@ private:
     std::vector<std::pair<size_t, size_t>> find_seam_string(const PrintObject *po,
             std::pair<size_t, size_t> start_seam,
             const SeamPlacerImpl::SeamComparator &comparator,
-            std::optional<std::pair<size_t, size_t>> &out_best_moved_seam,
-            size_t& out_moved_seams_count) const;
+            float& string_weight) const;
     std::optional<std::pair<size_t, size_t>> find_next_seam_in_layer(
             const std::vector<PrintObjectSeamData::LayerSeams> &layers,
-            const std::pair<size_t, size_t> &prev_point_index,
-            const size_t layer_idx, const float slice_z,
+            const Vec3f& projected_position,
+            const size_t layer_idx, const float max_distance,
             const SeamPlacerImpl::SeamComparator &comparator) const;
 };
 
diff --git a/src/libslic3r/Geometry/Curves.hpp b/src/libslic3r/Geometry/Curves.hpp
index 5a0d5481f..6ccdea366 100644
--- a/src/libslic3r/Geometry/Curves.hpp
+++ b/src/libslic3r/Geometry/Curves.hpp
@@ -15,8 +15,8 @@ template<int Dimension, typename NumberType>
 struct PolynomialCurve {
     Eigen::MatrixXf coefficients;
 
-    Vec3f get_fitted_value(const NumberType value) const {
-        auto result = Vec<Dimension, NumberType>::Zero();
+    Vec<Dimension, NumberType> get_fitted_value(const NumberType& value) const {
+        Vec<Dimension, NumberType> result = Vec<Dimension, NumberType>::Zero();
         size_t order = this->coefficients.rows() - 1;
         auto x = NumberType(1.);
         for (size_t index = 0; index < order + 1; ++index, x *= value)