Merge branch 'main' into enh-port-edit-gcode-dlg

src/libslic3r/CMakeLists.txt

@@ -157,6 +157,8 @@ set(lisbslic3r_sources
     GCode/PrintExtents.hpp
     GCode/RetractWhenCrossingPerimeters.cpp
     GCode/RetractWhenCrossingPerimeters.hpp
+    GCode/SmallAreaInfillFlowCompensator.cpp
+    GCode/SmallAreaInfillFlowCompensator.hpp
     GCode/SpiralVase.cpp
     GCode/SpiralVase.hpp
     GCode/SeamPlacer.cpp

src/libslic3r/GCode.cpp

@@ -1849,6 +1849,15 @@ static inline std::vector<const PrintInstance*> sort_object_instances_by_max_z(c
 //BBS: add sort logic for seq-print
 std::vector<const PrintInstance*> sort_object_instances_by_model_order(const Print& print, bool init_order)
 {
+    auto find_object_index = [](const Model& model, const ModelObject* obj) {
+        for (int index = 0; index < model.objects.size(); index++)
+        {
+            if (model.objects[index] == obj)
+                return index;
+        }
+        return -1;
+    };
+
     // Build up map from ModelInstance* to PrintInstance*
     std::vector<std::pair<const ModelInstance*, const PrintInstance*>> model_instance_to_print_instance;
     model_instance_to_print_instance.reserve(print.num_object_instances());
@@ -1856,10 +1865,16 @@ std::vector<const PrintInstance*> sort_object_instances_by_model_order(const Pri
         for (const PrintInstance &print_instance : print_object->instances())
         {
             if (init_order)
-                const_cast<ModelInstance*>(print_instance.model_instance)->arrange_order = print_instance.model_instance->id().id;
+                const_cast<ModelInstance*>(print_instance.model_instance)->arrange_order = find_object_index(print.model(), print_object->model_object());
             model_instance_to_print_instance.emplace_back(print_instance.model_instance, &print_instance);
         }
     std::sort(model_instance_to_print_instance.begin(), model_instance_to_print_instance.end(), [](auto &l, auto &r) { return l.first->arrange_order < r.first->arrange_order; });
+    if (init_order) {
+        // Re-assign the arrange_order so each instance has a unique order number
+        for (int k = 0; k < model_instance_to_print_instance.size(); k++) {
+            const_cast<ModelInstance*>(model_instance_to_print_instance[k].first)->arrange_order = k + 1;
+        }
+    }
 
     std::vector<const PrintInstance*> instances;
     instances.reserve(model_instance_to_print_instance.size());
@@ -1980,6 +1995,9 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
     } else
 	    m_enable_extrusion_role_markers = false;
 
+    if (!print.config().small_area_infill_flow_compensation_model.empty())
+        m_small_area_infill_flow_compensator = make_unique<SmallAreaInfillFlowCompensator>(print.config());
+
     // if thumbnail type of BTT_TFT, insert above header
     // if not, it is inserted under the header in its normal spot
     const GCodeThumbnailsFormat m_gcode_thumbnail_format = print.full_print_config().opt_enum<GCodeThumbnailsFormat>("thumbnails_format");
@@ -2215,8 +2233,11 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
         // In non-sequential print, the printing extruders may have been modified by the extruder switches stored in Model::custom_gcode_per_print_z.
         // Therefore initialize the printing extruders from there.
         this->set_extruders(tool_ordering.all_extruders());
-        // Order object instances using a nearest neighbor search.
-        print_object_instances_ordering = chain_print_object_instances(print);
+        print_object_instances_ordering = 
+            // By default, order object instances using a nearest neighbor search.
+            print.config().print_order == PrintOrder::Default ? chain_print_object_instances(print)
+            // Otherwise same order as the object list
+            : sort_object_instances_by_model_order(print);
     }
     if (initial_extruder_id == (unsigned int)-1) {
         // Nothing to print!
@@ -2540,9 +2561,7 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato
                     // another one, set first layer temperatures. This happens before the Z move
                     // is triggered, so machine has more time to reach such temperatures.
                     this->placeholder_parser().set("current_object_idx", int(finished_objects));
-                    //BBS: remove printing_by_object_gcode
-                    //std::string printing_by_object_gcode = this->placeholder_parser_process("printing_by_object_gcode", print.config().printing_by_object_gcode.value, initial_extruder_id);
-                    std::string printing_by_object_gcode;
+                    std::string printing_by_object_gcode = this->placeholder_parser_process("printing_by_object_gcode", print.config().printing_by_object_gcode.value, initial_extruder_id);
                     // Set first layer bed and extruder temperatures, don't wait for it to reach the temperature.
                     this->_print_first_layer_bed_temperature(file, print, printing_by_object_gcode, initial_extruder_id, false);
                     this->_print_first_layer_extruder_temperatures(file, print, printing_by_object_gcode, initial_extruder_id, false);
@@ -4020,7 +4039,7 @@ LayerResult GCode::process_layer(
         std::vector<InstanceToPrint> instances_to_print;
         bool has_prime_tower = print.config().enable_prime_tower
             && print.extruders().size() > 1
-            && (print.config().print_sequence == PrintSequence::ByLayer
+            && ((print.config().print_sequence == PrintSequence::ByLayer && print.config().print_order == PrintOrder::Default)
                 || (print.config().print_sequence == PrintSequence::ByObject && print.objects().size() == 1));
         if (has_prime_tower) {
             int plate_idx = print.get_plate_index();
@@ -5268,15 +5287,25 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
                 for (const Line& line : path.polyline.lines()) {
                     const double line_length = line.length() * SCALING_FACTOR;
                     path_length += line_length;
+                    auto dE = e_per_mm * line_length;
+                    if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) {
+                        auto oldE = dE;
+                        dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role());
+
+                        if (m_config.gcode_comments && oldE > 0 && oldE != dE) {
+                            description += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length);
+                        }
+                    }
                     gcode += m_writer.extrude_to_xy(
                         this->point_to_gcode(line.b),
-                        e_per_mm * line_length,
+                        dE,
                         GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion());
                 }
             } else {
                 // BBS: start to generate gcode from arc fitting data which includes line and arc
                 const std::vector<PathFittingData>& fitting_result = path.polyline.fitting_result;
                 for (size_t fitting_index = 0; fitting_index < fitting_result.size(); fitting_index++) {
+                    std::string tempDescription = description;
                     switch (fitting_result[fitting_index].path_type) {
                     case EMovePathType::Linear_move: {
                         size_t start_index = fitting_result[fitting_index].start_point_index;
@@ -5285,10 +5314,19 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
                             const Line line = Line(path.polyline.points[point_index - 1], path.polyline.points[point_index]);
                             const double line_length = line.length() * SCALING_FACTOR;
                             path_length += line_length;
+                            auto dE = e_per_mm * line_length;
+                            if (m_small_area_infill_flow_compensator  && m_config.small_area_infill_flow_compensation.value) {
+                                auto oldE = dE;
+                                dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role());
+
+                                if (m_config.gcode_comments && oldE > 0 && oldE != dE) {
+                                    tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length);
+                                }
+                            }
                             gcode += m_writer.extrude_to_xy(
                                 this->point_to_gcode(line.b),
-                                e_per_mm * line_length,
-                                GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion());
+                                dE,
+                                GCodeWriter::full_gcode_comment ? tempDescription : "", path.is_force_no_extrusion());
                         }
                         break;
                     }
@@ -5298,12 +5336,21 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
                         const double arc_length = fitting_result[fitting_index].arc_data.length * SCALING_FACTOR;
                         const Vec2d center_offset = this->point_to_gcode(arc.center) - this->point_to_gcode(arc.start_point);
                         path_length += arc_length;
+                        auto dE = e_per_mm * arc_length;
+                        if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) {
+                            auto oldE = dE;
+                            dE = m_small_area_infill_flow_compensator->modify_flow(arc_length, dE, path.role());
+
+                            if (m_config.gcode_comments && oldE > 0 && oldE != dE) {
+                                tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, arc_length);
+                            }
+                        }
                         gcode += m_writer.extrude_arc_to_xy(
                             this->point_to_gcode(arc.end_point),
                             center_offset,
-                            e_per_mm * arc_length,
+                            dE,
                             arc.direction == ArcDirection::Arc_Dir_CCW,
-                            GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion());
+                            GCodeWriter::full_gcode_comment ? tempDescription : "", path.is_force_no_extrusion());
                         break;
                     }
                     default:
@@ -5325,6 +5372,7 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
             pre_fan_enabled = check_overhang_fan(new_points[0].overlap, path.role());
 
         for (size_t i = 1; i < new_points.size(); i++) {
+            std::string tempDescription = description;
             const ProcessedPoint &processed_point = new_points[i];
             const ProcessedPoint &pre_processed_point = new_points[i-1];
             Vec2d p = this->point_to_gcode_quantized(processed_point.p);
@@ -5363,8 +5411,17 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description,
                 gcode += m_writer.set_speed(new_speed, "", comment);
                 last_set_speed = new_speed;
             }
+            auto dE = e_per_mm * line_length;
+            if (m_small_area_infill_flow_compensator  && m_config.small_area_infill_flow_compensation.value) {
+                auto oldE = dE;
+                dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role());
+
+                if (m_config.gcode_comments && oldE > 0 && oldE != dE) {
+                    tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length);
+                }
+            }
             gcode +=
-                m_writer.extrude_to_xy(p, e_per_mm * line_length, GCodeWriter::full_gcode_comment ? description : "");
+                m_writer.extrude_to_xy(p, dE, GCodeWriter::full_gcode_comment ? tempDescription : "");
 
             prev = p;
 

src/libslic3r/GCode.hpp

@@ -23,6 +23,7 @@
 #include "GCode/ExtrusionProcessor.hpp"
 
 #include "GCode/PressureEqualizer.hpp"
+#include "GCode/SmallAreaInfillFlowCompensator.hpp"
 
 #include <memory>
 #include <map>
@@ -536,6 +537,8 @@ private:
 
     std::unique_ptr<WipeTowerIntegration> m_wipe_tower;
 
+    std::unique_ptr<SmallAreaInfillFlowCompensator> m_small_area_infill_flow_compensator;
+
     // Heights (print_z) at which the skirt has already been extruded.
     std::vector<coordf_t>               m_skirt_done;
     // Has the brim been extruded already? Brim is being extruded only for the first object of a multi-object print.
@@ -598,6 +601,7 @@ private:
     friend class WipeTowerIntegration;
     friend class PressureEqualizer;
     friend class Print;
+    friend class SmallAreaInfillFlowCompensator;
 };
 
 std::vector<const PrintInstance*> sort_object_instances_by_model_order(const Print& print, bool init_order = false);

src/libslic3r/GCode/SmallAreaInfillFlowCompensator.cpp

@@ -0,0 +1,88 @@
+// Modify the flow of extrusion lines inversely proportional to the length of
+// the extrusion line. When infill lines get shorter the flow rate will auto-
+// matically be reduced to mitigate the effect of small infill areas being
+// over-extruded.
+
+// Based on original work by Alexander Þór licensed under the GPLv3:
+// https://github.com/Alexander-T-Moss/Small-Area-Flow-Comp
+
+#include <math.h>
+#include <cstring>
+#include <cfloat>
+
+#include "../libslic3r.h"
+#include "../PrintConfig.hpp"
+
+#include "SmallAreaInfillFlowCompensator.hpp"
+
+namespace Slic3r {
+
+bool nearly_equal(double a, double b)
+{
+    return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b && std::nextafter(a, std::numeric_limits<double>::max()) >= b;
+}
+
+SmallAreaInfillFlowCompensator::SmallAreaInfillFlowCompensator(const Slic3r::GCodeConfig& config)
+{
+    for (auto& line : config.small_area_infill_flow_compensation_model.values) {
+        std::istringstream iss(line);
+        std::string        value_str;
+        double             eLength = 0.0;
+
+        if (std::getline(iss, value_str, ',')) {
+            try {
+                eLength = std::stod(value_str);
+                if (std::getline(iss, value_str, ',')) {
+                    eLengths.push_back(eLength);
+                    flowComps.push_back(std::stod(value_str));
+                }
+            } catch (...) {
+                std::stringstream ss;
+                ss << "Error parsing data point in small area infill compensation model:" << line << std::endl;
+
+                throw Slic3r::InvalidArgument(ss.str());
+            }
+        }
+    }
+
+    for (int i = 0; i < eLengths.size(); i++) {
+        if (i == 0) {
+            if (!nearly_equal(eLengths[i], 0.0)) {
+                throw Slic3r::InvalidArgument("First extrusion length for small area infill compensation model must be 0");
+            }
+        } else {
+            if (nearly_equal(eLengths[i], 0.0)) {
+                throw Slic3r::InvalidArgument("Only the first extrusion length for small area infill compensation model can be 0");
+            }
+            if (eLengths[i] <= eLengths[i - 1]) {
+                throw Slic3r::InvalidArgument("Extrusion lengths for subsequent points must be increasing");
+            }
+        }
+    }
+
+    if (!flowComps.empty() && !nearly_equal(flowComps.back(), 1.0)) {
+        throw Slic3r::InvalidArgument("Final compensation factor for small area infill flow compensation model must be 1.0");
+    }
+
+    flowModel.set_points(eLengths, flowComps);
+}
+
+double SmallAreaInfillFlowCompensator::flow_comp_model(const double line_length)
+{
+    if (line_length == 0 || line_length > max_modified_length()) {
+        return 1.0;
+    }
+
+    return flowModel(line_length);
+}
+
+double SmallAreaInfillFlowCompensator::modify_flow(const double line_length, const double dE, const ExtrusionRole role)
+{
+    if (role == ExtrusionRole::erSolidInfill || role == ExtrusionRole::erTopSolidInfill || role == ExtrusionRole::erBottomSurface) {
+        return dE * flow_comp_model(line_length);
+    }
+
+    return dE;
+}
+
+} // namespace Slic3r

src/libslic3r/GCode/SmallAreaInfillFlowCompensator.hpp

@@ -0,0 +1,35 @@
+#ifndef slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_
+#define slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_
+
+#include "../libslic3r.h"
+#include "../PrintConfig.hpp"
+#include "../ExtrusionEntity.hpp"
+#include "spline/spline.h"
+
+namespace Slic3r {
+
+class SmallAreaInfillFlowCompensator
+{
+public:
+    SmallAreaInfillFlowCompensator() = delete;
+    explicit SmallAreaInfillFlowCompensator(const Slic3r::GCodeConfig& config);
+    ~SmallAreaInfillFlowCompensator() = default;
+
+    double modify_flow(const double line_length, const double dE, const ExtrusionRole role);
+
+private:
+    // Model points
+    std::vector<double> eLengths;
+    std::vector<double> flowComps;
+
+    // TODO: Cubic Spline
+    tk::spline flowModel;
+
+    double flow_comp_model(const double line_length);
+
+    double max_modified_length() { return eLengths.back(); }
+};
+
+} // namespace Slic3r
+
+#endif /* slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_ */

src/libslic3r/Model.cpp

@@ -836,7 +836,7 @@ end:
 // BBS: backup all in one dir
 std::string Model::get_auxiliary_file_temp_path()
 {
-    return get_backup_path("/Auxiliaries");
+    return get_backup_path("Auxiliaries");
 }
 
 // BBS: backup dir

src/libslic3r/Preset.cpp

@@ -786,7 +786,7 @@ static std::vector<std::string> s_Preset_print_options {
     "independent_support_layer_height",
     "support_angle", "support_interface_top_layers", "support_interface_bottom_layers",
     "support_interface_pattern", "support_interface_spacing", "support_interface_loop_pattern",
-    "support_top_z_distance", "support_on_build_plate_only","support_critical_regions_only", "bridge_no_support", "thick_bridges", "thick_internal_bridges","dont_filter_internal_bridges", "max_bridge_length", "print_sequence", "support_remove_small_overhang",
+    "support_top_z_distance", "support_on_build_plate_only","support_critical_regions_only", "bridge_no_support", "thick_bridges", "thick_internal_bridges","dont_filter_internal_bridges", "max_bridge_length", "print_sequence", "print_order", "support_remove_small_overhang",
     "filename_format", "wall_filament", "support_bottom_z_distance",
     "sparse_infill_filament", "solid_infill_filament", "support_filament", "support_interface_filament","support_interface_not_for_body",
     "ooze_prevention", "standby_temperature_delta", "interface_shells", "line_width", "initial_layer_line_width",
@@ -817,6 +817,7 @@ static std::vector<std::string> s_Preset_print_options {
      "wipe_tower_cone_angle", "wipe_tower_extra_spacing", "wipe_tower_extruder", "wiping_volumes_extruders","wipe_tower_bridging", "single_extruder_multi_material_priming",
      "wipe_tower_rotation_angle", "tree_support_branch_distance_organic", "tree_support_branch_diameter_organic", "tree_support_branch_angle_organic",
      "hole_to_polyhole", "hole_to_polyhole_threshold", "hole_to_polyhole_twisted", "mmu_segmented_region_max_width", "mmu_segmented_region_interlocking_depth",
+     "small_area_infill_flow_compensation", "small_area_infill_flow_compensation_model",
 };
 
 static std::vector<std::string> s_Preset_filament_options {

src/libslic3r/Print.cpp

@@ -314,6 +314,7 @@ bool Print::invalidate_state_by_config_options(const ConfigOptionResolver & /* n
             //|| opt_key == "resolution"
             //BBS: when enable arc fitting, we must re-generate perimeter
             || opt_key == "enable_arc_fitting"
+            || opt_key == "print_order"
             || opt_key == "wall_sequence") {
             osteps.emplace_back(posPerimeters);
             osteps.emplace_back(posEstimateCurledExtrusions);
@@ -1043,6 +1044,7 @@ boost::regex regex_g92e0 { "^[ \\t]*[gG]92[ \\t]*[eE](0(\\.0*)?|\\.0+)[ \\t]*(;.
 StringObjectException Print::validate(StringObjectException *warning, Polygons* collison_polygons, std::vector<std::pair<Polygon, float>>* height_polygons) const
 {
     std::vector<unsigned int> extruders = this->extruders();
+    unsigned int nozzles = m_config.nozzle_diameter.size();
 
     if (m_objects.empty())
         return {std::string()};
@@ -1050,7 +1052,7 @@ StringObjectException Print::validate(StringObjectException *warning, Polygons*
     if (extruders.empty())
         return { L("No extrusions under current settings.") };
 
-    if (extruders.size() > 1 && m_config.print_sequence != PrintSequence::ByObject) {
+    if (nozzles < 2 && extruders.size() > 1 && m_config.print_sequence != PrintSequence::ByObject) {
         auto ret = check_multi_filament_valid(*this);
         if (!ret.string.empty())
         {

src/libslic3r/PrintConfig.cpp

@@ -195,6 +195,12 @@ static t_config_enum_values s_keys_map_PrintSequence {
 };
 CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(PrintSequence)
 
+static t_config_enum_values s_keys_map_PrintOrder{
+    { "default",     int(PrintOrder::Default) },
+    { "as_obj_list", int(PrintOrder::AsObjectList)},
+};
+CONFIG_OPTION_ENUM_DEFINE_STATIC_MAPS(PrintOrder)
+
 static t_config_enum_values s_keys_map_SlicingMode {
     { "regular",        int(SlicingMode::Regular) },
     { "even_odd",       int(SlicingMode::EvenOdd) },
@@ -1171,6 +1177,17 @@ void PrintConfigDef::init_fff_params()
     def->mode = comSimple;
     def->set_default_value(new ConfigOptionEnum<PrintSequence>(PrintSequence::ByLayer));
 
+    def = this->add("print_order", coEnum);
+    def->label = L("Layer order");
+    def->tooltip = L("Print order within a single layer");
+    def->enum_keys_map = &ConfigOptionEnum<PrintOrder>::get_enum_values();
+    def->enum_values.push_back("default");
+    def->enum_values.push_back("as_obj_list");
+    def->enum_labels.push_back(L("Default"));
+    def->enum_labels.push_back(L("As object list"));
+    def->mode = comAdvanced;
+    def->set_default_value(new ConfigOptionEnum<PrintOrder>(PrintOrder::Default));
+
     def = this->add("slow_down_for_layer_cooling", coBools);
     def->label = L("Slow printing down for better layer cooling");
     def->tooltip = L("Enable this option to slow printing speed down to make the final layer time not shorter than "
@@ -2668,6 +2685,26 @@ def = this->add("filament_loading_speed", coFloats);
     def->mode = comAdvanced;
     def->set_default_value(new ConfigOptionString(""));
 
+    def = this->add("small_area_infill_flow_compensation", coBool);
+    def->label = L("Enable Flow Compensation");
+    def->tooltip = L("Enable flow compensation for small infill areas");
+    def->mode = comAdvanced;
+    def->set_default_value(new ConfigOptionBool(false));
+
+    def = this->add("small_area_infill_flow_compensation_model", coStrings);
+    def->label = L("Flow Compensation Model");
+    def->tooltip = L(
+        "Flow Compensation Model, used to adjust the flow for small infill "
+        "areas. The model is expressed as a comma separated pair of values for "
+        "extrusion length and flow correction factors, one per line, in the "
+        "following format: \"1.234,5.678\"");
+    def->mode = comAdvanced;
+    def->gui_flags = "serialized";
+    def->multiline = true;
+    def->full_width = true;
+    def->height = 15;
+    def->set_default_value(new ConfigOptionStrings{"0,0", "\n0.2,0.4444", "\n0.4,0.6145", "\n0.6,0.7059", "\n0.8,0.7619", "\n1.5,0.8571", "\n2,0.8889", "\n3,0.9231", "\n5,0.9520", "\n10,1"});
+
     {
         struct AxisDefault {
             std::string         name;

src/libslic3r/PrintConfig.hpp

@@ -106,6 +106,13 @@ enum class PrintSequence {
     Count,
 };
 
+enum class PrintOrder
+{
+    Default,
+    AsObjectList,
+    Count,
+};
+
 enum class SlicingMode
 {
     // Regular, applying ClipperLib::pftNonZero rule when creating ExPolygons.
@@ -920,6 +927,7 @@ PRINT_CONFIG_CLASS_DEFINE(
     
     ((ConfigOptionEnum<WallSequence>,  wall_sequence))
     ((ConfigOptionBool,                is_infill_first))
+    ((ConfigOptionBool,                small_area_infill_flow_compensation))
 )
 
 PRINT_CONFIG_CLASS_DEFINE(
@@ -1067,6 +1075,8 @@ PRINT_CONFIG_CLASS_DEFINE(
     ((ConfigOptionBool,                enable_filament_ramming))
     ((ConfigOptionBool,                support_multi_bed_types))
 
+    // Small Area Infill Flow Compensation
+    ((ConfigOptionStrings,              small_area_infill_flow_compensation_model))
 )
 
 // This object is mapped to Perl as Slic3r::Config::Print.
@@ -1098,6 +1108,7 @@ PRINT_CONFIG_CLASS_DERIVED_DEFINE(
     ((ConfigOptionInts,               overhang_fan_speed))
     ((ConfigOptionEnumsGeneric,       overhang_fan_threshold))
     ((ConfigOptionEnum<PrintSequence>,print_sequence))
+    ((ConfigOptionEnum<PrintOrder>,   print_order))
     ((ConfigOptionInts,               first_layer_print_sequence))
     ((ConfigOptionBools,              slow_down_for_layer_cooling))
     ((ConfigOptionInts,               close_fan_the_first_x_layers))

src/libslic3r/PrintObject.cpp

@@ -928,6 +928,10 @@ bool PrintObject::invalidate_state_by_config_options(
             || opt_key == "wipe_on_loops"
             || opt_key == "wipe_speed") {
             steps.emplace_back(posPerimeters);
+        } else if (
+               opt_key == "small_area_infill_flow_compensation"
+            || opt_key == "small_area_infill_flow_compensation_model") {
+            steps.emplace_back(posSlice);
         } else if (opt_key == "gap_infill_speed"
             || opt_key == "filter_out_gap_fill" ) {
             // Return true if gap-fill speed has changed from zero value to non-zero or from non-zero value to zero.

src/slic3r/GUI/Auxiliary.cpp

@@ -33,17 +33,29 @@ wxDEFINE_EVENT(EVT_AUXILIARY_IMPORT, wxCommandEvent);
 wxDEFINE_EVENT(EVT_AUXILIARY_UPDATE_COVER, wxCommandEvent);
 wxDEFINE_EVENT(EVT_AUXILIARY_UPDATE_DELETE, wxCommandEvent);
 wxDEFINE_EVENT(EVT_AUXILIARY_UPDATE_RENAME, wxCommandEvent);
+wxDEFINE_EVENT(EVT_AUXILIARY_DONE, wxCommandEvent);
 
 
 const std::vector<std::string> license_list = {
-    "BSD License",
-    "Apache License",
-    "GPL License",
-    "LGPL License",
-    "MIT License",
-    "CC License"
+    "",
+    "CC0",
+    "BY",
+    "BY-SA",
+    "BY-ND",
+    "BY-NC",
+    "BY-NC-SA",
+    "BY-NC-ND",
 };
 
+static std::shared_ptr<ModelInfo> ensure_model_info()
+{
+    auto& model = wxGetApp().plater()->model();
+    if (model.model_info == nullptr) {
+        model.model_info = std::make_shared<ModelInfo>();
+    }
+    return model.model_info;
+}
+
 AuFile::AuFile(wxWindow *parent, fs::path file_path, wxString file_name, AuxiliaryFolderType type, wxWindowID id, const wxPoint &pos, const wxSize &size, long style)
 {
     m_type      = type;
@@ -344,7 +356,7 @@ void AuFile::on_input_enter(wxCommandEvent &evt)
     }
 
     auto     existing  = false;
-    auto     dir       = m_file_path.branch_path();
+    auto     dir       = m_file_path.parent_path();
     auto     new_fullname = new_file_name + m_file_path.extension().string();
 
     
@@ -454,14 +466,12 @@ void AuFile::on_mouse_left_up(wxMouseEvent &evt)
 
 void AuFile::on_set_cover()
 {
-    if (wxGetApp().plater()->model().model_info == nullptr) { wxGetApp().plater()->model().model_info = std::make_shared<ModelInfo>(); }
-
     fs::path path(into_path(m_file_name));
-    wxGetApp().plater()->model().model_info->cover_file = path.string();
+    ensure_model_info()->cover_file = path.string();
     //wxGetApp().plater()->model().model_info->cover_file = m_file_name.ToStdString();
 
-    auto full_path          = m_file_path.branch_path();
-    auto full_root_path         = full_path.branch_path();
+    auto full_path          = m_file_path.parent_path();
+    auto full_root_path         = full_path.parent_path();
     auto full_root_path_str = encode_path(full_root_path.string().c_str());
     auto dir       = wxString::Format("%s/.thumbnails", full_root_path_str);
 
@@ -505,8 +515,8 @@ void AuFile::on_set_delete()
     auto     is_fine = fs::remove(bfs_path);
 
     if (m_cover) {
-        auto full_path          = m_file_path.branch_path();
-        auto full_root_path     = full_path.branch_path();
+        auto full_path          = m_file_path.parent_path();
+        auto full_root_path     = full_path.parent_path();
         auto full_root_path_str = encode_path(full_root_path.string().c_str());
         auto dir                = wxString::Format("%s/.thumbnails", full_root_path_str);
         fs::path dir_path(dir.ToStdWstring());
@@ -520,8 +530,11 @@ void AuFile::on_set_delete()
         if (fs::exists(fs::path(middle_img_path))) { fs::remove(fs::path(middle_img_path)); }
     }
 
-    if (wxGetApp().plater()->model().model_info == nullptr) { wxGetApp().plater()->model().model_info = std::make_shared<ModelInfo>(); }
-    if (wxGetApp().plater()->model().model_info->cover_file == m_file_name) { wxGetApp().plater()->model().model_info->cover_file = ""; }
+    if (wxGetApp().plater()->model().model_info != nullptr) {
+        if (wxGetApp().plater()->model().model_info->cover_file == m_file_name) {
+            wxGetApp().plater()->model().model_info->cover_file = "";
+        }
+    }
 
     if (is_fine) {
         auto evt = wxCommandEvent(EVT_AUXILIARY_UPDATE_DELETE);
@@ -669,6 +682,7 @@ void AuFolderPanel::update(std::vector<fs::path> paths)
     }
     m_gsizer_content->Layout();
     Layout();
+    Refresh();
 }
 
 void AuFolderPanel::msw_rescale() 
@@ -820,9 +834,22 @@ void AuxiliaryPanel::init_bitmap()
 
 void AuxiliaryPanel::init_tabpanel()
 {
-    auto        m_side_tools     = new wxPanel(this, wxID_ANY, wxDefaultPosition, wxSize(FromDIP(220), FromDIP(18)));
+    StateColor btn_bg_green(std::pair<wxColour, int>(wxColour(206, 206, 206), StateColor::Disabled),
+                            std::pair<wxColour, int>(wxColour(0, 137, 123), StateColor::Pressed),
+                            std::pair<wxColour, int>(wxColour(38, 166, 154), StateColor::Hovered),
+                            std::pair<wxColour, int>(wxColour(0, 150, 136), StateColor::Normal));
+    auto back_btn = new Button(this, _L("Back"), "assemble_return", wxBORDER_NONE | wxBU_LEFT | wxBU_EXACTFIT);
+    back_btn->SetSize(wxSize(FromDIP(220), FromDIP(18)));
+    back_btn->SetBackgroundColor(btn_bg_green);
+    back_btn->SetCornerRadius(0);
+    back_btn->Bind(wxEVT_COMMAND_BUTTON_CLICKED, [this](wxEvent& e) {
+        auto event = wxCommandEvent(EVT_AUXILIARY_DONE);
+        event.SetEventObject(m_parent);
+        wxPostEvent(m_parent, event);
+    });
+
     wxBoxSizer *sizer_side_tools = new wxBoxSizer(wxVERTICAL);
-    sizer_side_tools->Add(m_side_tools, 1, wxEXPAND, 0);
+    sizer_side_tools->Add(back_btn, 1, wxEXPAND, 0);
     m_tabpanel = new Tabbook(this, wxID_ANY, wxDefaultPosition, wxDefaultSize, sizer_side_tools, wxNB_LEFT | wxTAB_TRAVERSAL | wxNB_NOPAGETHEME);
     m_tabpanel->SetBackgroundColour(wxColour("#FEFFFF"));
     m_tabpanel->Bind(wxEVT_BOOKCTRL_PAGE_CHANGED, [this](wxBookCtrlEvent &e) { ; });
@@ -872,20 +899,7 @@ bool AuxiliaryPanel::Show(bool show) { return wxPanel::Show(show); }
 void AuxiliaryPanel::init_auxiliary()
 {
     Model &model = wxGetApp().plater()->model();
-    m_root_dir   = encode_path(model.get_auxiliary_file_temp_path().c_str());
-    if (wxDirExists(m_root_dir)) {
-        fs::path path_to_del(m_root_dir.ToStdWstring());
-        try {
-            fs::remove_all(path_to_del);
-        } catch (...) {
-            BOOST_LOG_TRIVIAL(error) << "Failed  removing the auxiliary directory " << m_root_dir.c_str();
-        }
-    }
-
-    fs::path top_dir_path(m_root_dir.ToStdWstring());
-    fs::create_directory(top_dir_path);
-
-    for (auto folder : s_default_folders) create_folder(folder);
+    Reload(encode_path(model.get_auxiliary_file_temp_path().c_str()), {});
 }
 
 void AuxiliaryPanel::on_import_file(wxCommandEvent &event)
@@ -947,7 +961,7 @@ void AuxiliaryPanel::on_import_file(wxCommandEvent &event)
            
 
             boost::system::error_code ec;
-            if (!fs::copy_file(src_bfs_path, fs::path(dir_path.ToStdWstring()), fs::copy_option::overwrite_if_exists, ec)) continue;
+            if (!fs::copy_file(src_bfs_path, fs::path(dir_path.ToStdWstring()), fs::copy_options::overwrite_existing, ec)) continue;
             Slic3r::put_other_changes();
 
             // add in file list
@@ -987,76 +1001,22 @@ std::string AuxiliaryPanel::replaceSpace(std::string s, std::string ts, std::str
     return s;
 }
 
-void AuxiliaryPanel::Reload(wxString aux_path)
+void AuxiliaryPanel::Reload(wxString aux_path, std::map<std::string, std::vector<json>> paths)
 {
-    fs::path new_aux_path(aux_path.ToStdWstring());
-
-    try {
-        fs::remove_all(fs::path(m_root_dir.ToStdWstring()));
-    } catch (...) {
-        BOOST_LOG_TRIVIAL(error) << "Failed  removing the auxiliary directory " << m_root_dir.c_str();
-    }
-
     m_root_dir = aux_path;
     m_paths_list.clear();
-    // Check new path. If not exist, create a new one.
-    if (!fs::exists(new_aux_path)) {
-        fs::create_directory(new_aux_path);
-        // Create default folders if they are not loaded
-        for (auto folder : s_default_folders) {
-            wxString folder_path = aux_path + "/" + folder;
-            if (fs::exists(folder_path.ToStdWstring())) continue;
-            fs::create_directory(folder_path.ToStdWstring());
+
+    for (const auto & path : paths) {
+        m_paths_list[path.first] = std::vector<fs::path>{};
+        for (const auto & j : path.second) {
+            m_paths_list[path.first].push_back(j["_filepath"]);
         }
-        update_all_panel();
-        m_designer_panel->update_info();
-        return;
-    }
-
-    // Load from new path
-    std::vector<fs::path>  dir_cache;
-    fs::directory_iterator iter_end;
-
-    for (fs::directory_iterator iter(new_aux_path); iter != iter_end; iter++) {
-        wxString path = iter->path().generic_wstring();
-        dir_cache.push_back(iter->path());
-    }
-
-    for (auto dir : dir_cache) {
-        for (fs::directory_iterator iter(dir); iter != iter_end; iter++) {
-            if (fs::is_directory(iter->path())) continue;
-            wxString file_path     = iter->path().generic_wstring();
-            //auto     file_path_str = encode_path(file_path.c_str());
-
-            for (auto folder : s_default_folders) {
-                auto idx = file_path.find(folder.ToStdString());
-                if (idx != std::string::npos) {
-                    auto iter = m_paths_list.find(folder.ToStdString());
-                    auto     file_path_str = fs::path(file_path.ToStdWstring());
-
-                    if (iter != m_paths_list.end()) {
-                        m_paths_list[folder.ToStdString()].push_back(file_path_str);
-                        break;
-                    } else {
-                        m_paths_list[folder.ToStdString()] = std::vector<fs::path>{file_path_str};
-                        break;
-                    }
-                }
-            }
-        }
-    }
-
-    // Create default folders if they are not loaded
-    wxDataViewItemArray default_items;
-    for (auto folder : s_default_folders) {
-        wxString folder_path = aux_path + "/" + folder;
-        if (fs::exists(folder_path.ToStdWstring())) continue;
-        fs::create_directory(folder_path.ToStdWstring());
     }
 
     update_all_panel();
     update_all_cover();
     m_designer_panel->update_info();
+    m_tabpanel->SetSelection(0);
 }
 
 void AuxiliaryPanel::update_all_panel()
@@ -1121,22 +1081,21 @@ void AuxiliaryPanel::update_all_cover()
      m_imput_model_name->GetTextCtrl()->SetSize(wxSize(FromDIP(450), -1));
      m_sizer_model_name->Add(m_imput_model_name, 0, wxALIGN_CENTER, 0);
 
-     /*
      wxBoxSizer *m_sizer_license = new wxBoxSizer(wxHORIZONTAL);
-     auto m_text_license = new wxStaticText(this, wxID_ANY, _L("License"), wxDefaultPosition, wxSize(120, -1), 0);
+     auto m_text_license = new wxStaticText(this, wxID_ANY, _L("License"), wxDefaultPosition, wxSize(180, -1), 0);
      m_text_license->Wrap(-1);
      m_sizer_license->Add(m_text_license, 0, wxALIGN_CENTER, 0);
 
-     m_combo_license = new wxComboBox(this, wxID_ANY, wxEmptyString, wxDefaultPosition, wxSize(450, -1), 0, NULL, wxCB_READONLY);
+     m_combo_license = new ComboBox(this, wxID_ANY, wxEmptyString, wxDefaultPosition, wxSize(FromDIP(450), -1), 0, NULL, wxCB_READONLY);
      m_sizer_license->Add(m_combo_license, 0, wxALIGN_CENTER, 0);
-     */
+
      m_sizer_body->Add( 0, 0, 0, wxTOP, FromDIP(50) );
      m_sizer_body->Add(m_sizer_designer, 0, wxLEFT, FromDIP(50));
      m_sizer_body->Add( 0, 0, 0, wxTOP, FromDIP(20));
      m_sizer_body->Add(m_sizer_model_name, 0, wxLEFT, FromDIP(50));
-     //m_sizer_body->Add(0, 0, 0, wxTOP, FromDIP(20));
-     //m_sizer_body->Add(m_sizer_license, 0, wxLEFT, FromDIP(50));
-     //init_license_list();
+     m_sizer_body->Add(0, 0, 0, wxTOP, FromDIP(20));
+     m_sizer_body->Add(m_sizer_license, 0, wxLEFT, FromDIP(50));
+     init_license_list();
 
      SetSizer(m_sizer_body);
      Layout();
@@ -1144,52 +1103,35 @@ void AuxiliaryPanel::update_all_cover()
 
      m_input_designer->Bind(wxEVT_TEXT, &DesignerPanel::on_input_enter_designer, this);
      m_imput_model_name->Bind(wxEVT_TEXT, &DesignerPanel::on_input_enter_model, this);
-     //m_combo_license->Connect(wxEVT_COMMAND_COMBOBOX_SELECTED, wxCommandEventHandler(DesignerPanel::on_select_license), NULL, this);
+     m_combo_license->Bind(wxEVT_COMMAND_COMBOBOX_SELECTED, &DesignerPanel::on_select_license, this);
 }
 
  DesignerPanel::~DesignerPanel()
  {
-     //m_combo_license->Disconnect(wxEVT_COMMAND_COMBOBOX_SELECTED, wxCommandEventHandler(DesignerPanel::on_select_license), NULL, this);
  }
 
  void DesignerPanel::init_license_list()
  {
-     /*
      wxArrayString text_licese;
      for (int i = 0; i < license_list.size(); i++) {
          text_licese.Add(license_list[i]);
      }
      m_combo_license->Set(text_licese);
-     */
  }
 
  void DesignerPanel::on_select_license(wxCommandEvent&evt)
  {
      int selected = evt.GetInt();
      if (selected >= 0 && selected < license_list.size()) {
-         if (wxGetApp().plater()->model().model_info == nullptr) {
-             wxGetApp().plater()->model().model_info = std::make_shared<ModelInfo>();
-         }
-         if (wxGetApp().plater()->model().model_info != nullptr) {
-             wxGetApp().plater()->model().model_info->license = license_list[selected];
-         }
+         ensure_model_info()->license = license_list[selected];
      }
  }
 
-bool DesignerPanel::Show(bool show) 
-{
-    if ( wxGetApp().plater()->model().design_info != nullptr) {
-        wxString text = wxString::FromUTF8(wxGetApp().plater()->model().design_info->Designer);
-        m_input_designer->GetTextCtrl()->SetValue(text);
-    }
-
-     if (wxGetApp().plater()->model().model_info != nullptr) { 
-         wxString text = wxString::FromUTF8(wxGetApp().plater()->model().model_info->model_name);
-         m_imput_model_name->GetTextCtrl()->SetValue(text);
-     }
-    
-    return wxPanel::Show(show);
-}
+bool DesignerPanel::Show(bool show)
+ {
+     if (show) update_info();
+     return wxPanel::Show(show);
+ }
 
 void DesignerPanel::on_input_enter_designer(wxCommandEvent &evt) 
 { 
@@ -1200,9 +1142,7 @@ void DesignerPanel::on_input_enter_designer(wxCommandEvent &evt)
 void DesignerPanel::on_input_enter_model(wxCommandEvent &evt) 
 {
     auto text   = evt.GetString();
-    if (wxGetApp().plater()->model().model_info) {
-        wxGetApp().plater()->model().model_info->model_name = std::string(text.ToUTF8().data());
-    }
+    ensure_model_info()->model_name = std::string(text.ToUTF8().data());
 }
 
 void DesignerPanel::update_info() 
@@ -1215,10 +1155,13 @@ void DesignerPanel::update_info()
     }
 
     if (wxGetApp().plater()->model().model_info != nullptr) {
-        wxString text = wxString::FromUTF8(wxGetApp().plater()->model().model_info->model_name);
-        m_imput_model_name->GetTextCtrl()->SetValue(text);
+        m_imput_model_name->GetTextCtrl()->SetValue(wxString::FromUTF8(wxGetApp().plater()->model().model_info->model_name));
+        if (!m_combo_license->SetStringSelection(wxString::FromUTF8(wxGetApp().plater()->model().model_info->license))) {
+            m_combo_license->SetSelection(0);
+        }
     } else {
-         m_imput_model_name->GetTextCtrl()->SetValue(wxEmptyString);
+        m_imput_model_name->GetTextCtrl()->SetValue(wxEmptyString);
+        m_combo_license->SetSelection(0);
     }
 }
 
@@ -1226,6 +1169,7 @@ void DesignerPanel::msw_rescale()
 {
     m_input_designer->GetTextCtrl()->SetSize(wxSize(FromDIP(450), -1));
     m_imput_model_name->GetTextCtrl()->SetSize(wxSize(FromDIP(450), -1));
+    m_combo_license->SetSize(wxSize(FromDIP(450), -1));
 }
 
 }} // namespace Slic3r::GUI

src/slic3r/GUI/Auxiliary.hpp

@@ -181,7 +181,7 @@ public:
 
     ::TextInput*        m_input_designer {nullptr};
     ::TextInput*        m_imput_model_name {nullptr};
-    //wxComboBox*         m_combo_license {nullptr};
+    ComboBox*           m_combo_license {nullptr};
     bool Show(bool show) override;
     void                init_license_list();
     void                on_input_enter_designer(wxCommandEvent &evt);
@@ -232,7 +232,7 @@ public:
     void                                            create_folder(wxString name = wxEmptyString);
     std::string                                     replaceSpace(std::string s, std::string ts, std::string ns);
     void                                            on_import_file(wxCommandEvent &event);
-    void                                            Reload(wxString aux_path);
+    void                                            Reload(wxString aux_path, std::map<std::string, std::vector<json>> paths);
 
     void update_all_panel();
     void update_all_cover();
@@ -242,6 +242,7 @@ wxDECLARE_EVENT(EVT_AUXILIARY_IMPORT, wxCommandEvent);
 wxDECLARE_EVENT(EVT_AUXILIARY_UPDATE_COVER, wxCommandEvent);
 wxDECLARE_EVENT(EVT_AUXILIARY_UPDATE_DELETE, wxCommandEvent);
 wxDECLARE_EVENT(EVT_AUXILIARY_UPDATE_RENAME, wxCommandEvent);
+wxDECLARE_EVENT(EVT_AUXILIARY_DONE, wxCommandEvent);
 }} // namespace Slic3r::GUI
 
 #endif

src/slic3r/GUI/ConfigManipulation.cpp

@@ -661,9 +661,10 @@ void ConfigManipulation::toggle_print_fff_options(DynamicPrintConfig *config, co
     for (auto el : { "ironing_pattern", "ironing_flow", "ironing_spacing", "ironing_speed", "ironing_angle" })
         toggle_line(el, has_ironing);
     
-    // bool have_sequential_printing = (config->opt_enum<PrintSequence>("print_sequence") == PrintSequence::ByObject);
+    bool have_sequential_printing = (config->opt_enum<PrintSequence>("print_sequence") == PrintSequence::ByObject);
     // for (auto el : { "extruder_clearance_radius", "extruder_clearance_height_to_rod", "extruder_clearance_height_to_lid" })
     //     toggle_field(el, have_sequential_printing);
+    toggle_field("print_order", !have_sequential_printing);
     
     bool have_ooze_prevention = config->opt_bool("ooze_prevention");
     toggle_field("standby_temperature_delta", have_ooze_prevention);
@@ -742,6 +743,10 @@ void ConfigManipulation::toggle_print_fff_options(DynamicPrintConfig *config, co
         apply(config, &new_conf);
     }
     toggle_line("timelapse_type", is_BBL_Printer);
+
+
+    bool have_small_area_infill_flow_compensation = config->opt_bool("small_area_infill_flow_compensation");
+    toggle_line("small_area_infill_flow_compensation_model", have_small_area_infill_flow_compensation);
 }
 
 void ConfigManipulation::update_print_sla_config(DynamicPrintConfig* config, const bool is_global_config/* = false*/)

src/slic3r/GUI/GLCanvas3D.cpp

@@ -7304,10 +7304,10 @@ void GLCanvas3D::_render_overlays()
 
 	auto curr_plate = wxGetApp().plater()->get_partplate_list().get_curr_plate();
     auto curr_print_seq = curr_plate->get_real_print_seq();
-    bool sequential_print = (curr_print_seq == PrintSequence::ByObject);
+    const Print* print = fff_print();
+    bool sequential_print = (curr_print_seq == PrintSequence::ByObject) || print->config().print_order == PrintOrder::AsObjectList;
     std::vector<const ModelInstance*> sorted_instances;
     if (sequential_print) {
-        const Print* print = fff_print();
         if (print) {
             for (const PrintObject *print_object : print->objects())
             {

src/slic3r/GUI/GUI_App.cpp

@@ -2656,7 +2656,7 @@ bool GUI_App::on_init_inner()
 
     sidebar().obj_list()->init();
     //sidebar().aux_list()->init_auxiliary();
-    //mainframe->m_auxiliary->init_auxiliary();
+    mainframe->m_project->init_auxiliary();
 
 //     update_mode(); // !!! do that later
     SetTopWindow(mainframe);

src/slic3r/GUI/GUI_ObjectList.cpp

@@ -1545,26 +1545,6 @@ void ObjectList::OnBeginDrag(wxDataViewEvent &event)
     }
 
     if (type & itObject) {
-        int curr_obj_id = m_objects_model->GetIdByItem(event.GetItem());
-        PartPlateList& partplate_list = wxGetApp().plater()->get_partplate_list();
-        int from_plate = partplate_list.find_instance(curr_obj_id, 0);
-        if (from_plate == -1) {
-            event.Veto();
-            return;
-        }
-        auto curr_plate_seq = partplate_list.get_plate(from_plate)->get_print_seq();
-        if (curr_plate_seq == PrintSequence::ByDefault) {
-            auto curr_preset_config = wxGetApp().preset_bundle->prints.get_edited_preset().config;
-            if (curr_preset_config.has("print_sequence"))
-                curr_plate_seq = curr_preset_config.option<ConfigOptionEnum<PrintSequence>>("print_sequence")->value;
-        }
-
-        if (curr_plate_seq != PrintSequence::ByObject) {
-            //drag forbidden under bylayer mode
-            event.Veto();
-            return;
-        }
-
         m_dragged_data.init(m_objects_model->GetIdByItem(item), type);
     }
     else if (type & itVolume){

src/slic3r/GUI/Gizmos/GLGizmoCut.cpp

@@ -717,7 +717,12 @@ bool GLGizmoCut3D::render_reset_button(const std::string& label_id, const std::s
 
 static double get_grabber_mean_size(const BoundingBoxf3& bb)
 {
+#if ENABLE_FIXED_GRABBER
+    // Orca: make grabber larger
+    return 32. * GLGizmoBase::INV_ZOOM;
+#else
     return (bb.size().x() + bb.size().y() + bb.size().z()) / 30.;
+#endif
 }
 
 indexed_triangle_set GLGizmoCut3D::its_make_groove_plane()
@@ -2503,13 +2508,13 @@ void GLGizmoCut3D::add_horizontal_shift(float shift)
 
 void GLGizmoCut3D::render_color_marker(float size, const ImU32& color)
 {
-    ImGui::SameLine();
     const float radius = 0.5f * size;
     ImVec2 pos = ImGui::GetCurrentWindow()->DC.CursorPos;
-    pos.x += size;
-    pos.y += 1.25f * radius;
+    pos.x += radius;
+    pos.y += 1.4f * radius;
     ImGui::GetCurrentWindow()->DrawList->AddNgonFilled(pos, radius, color, 6);
-    m_imgui->text(" ");
+    m_imgui->text("  ");
+    ImGui::SameLine();
 }
 
 void GLGizmoCut3D::render_groove_float_input(const std::string& label, float& in_val, const float& init_val, float& in_tolerance)
@@ -2730,20 +2735,27 @@ void GLGizmoCut3D::render_cut_plane_input_window(CutConnectors &connectors, floa
 
         // render "After Cut" section
 
-        float label_width = 0;
+        ImVec2 label_size;
         for (const wxString &label : {_L("Upper part"), _L("Lower part")}) {
-            const float width = m_imgui->calc_text_size(label).x + m_imgui->scaled(1.5f);
-            if (label_width < width)
-                    label_width = width;
+            const ImVec2 text_size = ImGuiWrapper::calc_text_size(label);
+            if (label_size.x < text_size.x)
+                label_size.x = text_size.x;
+            if (label_size.y < text_size.y)
+                label_size.y = text_size.y;
         }
 
-        auto render_part_action_line = [this, label_width, &connectors](const wxString &label, const wxString &suffix, bool &keep_part,
+        const float marker_size = label_size.y;
+        const float h_shift     = marker_size + label_size.x + m_imgui->scaled(2.f);
+
+        auto render_part_action_line = [this, h_shift, marker_size, &connectors](const wxString &label, const wxString &suffix, bool &keep_part,
                                                                         bool &place_on_cut_part, bool &rotate_part) {
             bool keep = true;
+
             ImGui::AlignTextToFramePadding();
+            render_color_marker(marker_size, ImGuiWrapper::to_ImU32(suffix == "##upper" ? UPPER_PART_COLOR : LOWER_PART_COLOR));
             m_imgui->text(label);
 
-            ImGui::SameLine(label_width);
+            ImGui::SameLine(h_shift);
 
             m_imgui->disabled_begin(!connectors.empty() || m_keep_as_parts);
             m_imgui->bbl_checkbox(_L("Keep") + suffix, connectors.empty() ? keep_part : keep);

src/slic3r/GUI/PartPlate.cpp

@@ -51,6 +51,7 @@ static unsigned int GLOBAL_PLATE_INDEX = 0;
 
 static const double LOGICAL_PART_PLATE_GAP = 1. / 5.;
 static const int PARTPLATE_ICON_SIZE = 16;
+static const int PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE = 12;
 static const int PARTPLATE_ICON_GAP_TOP = 3;
 static const int PARTPLATE_ICON_GAP_LEFT = 3;
 static const int PARTPLATE_ICON_GAP_Y = 5;
@@ -571,6 +572,42 @@ void PartPlate::calc_vertex_for_number(int index, bool one_number, GLModel &buff
 		BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << "Unable to generate geometry buffers for icons\n";
 }
 
+void PartPlate::calc_vertex_for_plate_name_edit_icon(GLTexture *texture, int index, PickingModel &model) {
+    model.reset();
+
+	auto    bed_ext = get_extents(m_shape);
+	auto    factor  = bed_ext.size()(1) / 200.0;
+	wxCoord w, h;
+	h = int(factor * 16);
+	ExPolygon poly;
+	Vec2d     p        = bed_ext[3];
+	float     offset_x = 1;
+	h = PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE;
+    p += Vec2d(0, PARTPLATE_TEXT_OFFSET_Y + h);
+	if (texture && texture->get_width() > 0 && texture->get_height()) {
+		w    = int(factor * (texture->get_original_width() * 16) / texture->get_height()) + 1;
+
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + w), scale_(p(1) - h )});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + w + PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE), scale_(p(1) - h)});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + w + PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE), scale_(p(1))});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + w), scale_(p(1) )});
+
+		if (!init_model_from_poly(model.model, poly, GROUND_Z))
+			BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << "Unable to generate geometry buffers for icons\n";
+	} else {
+
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + offset_x ), scale_(p(1) - h )});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + offset_x + PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE), scale_(p(1) - h)});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + offset_x + PARTPLATE_EDIT_PLATE_NAME_ICON_SIZE), scale_(p(1))});
+		poly.contour.append({scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + offset_x), scale_(p(1) )});
+
+		if (!init_model_from_poly(model.model, poly, GROUND_Z))
+			BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << "Unable to generate geometry buffers for icons\n";
+    }
+
+    init_raycaster_from_model(model);
+}
+
 void PartPlate::calc_vertex_for_icons(int index, PickingModel &model)
 {
     model.reset();
@@ -975,6 +1012,13 @@ void PartPlate::render_icons(bool bottom, bool only_name, int hover_id)
                     render_icon_texture(m_lock_icon.model, m_partplate_list->m_lockopen_texture);
             }
 
+			if (hover_id == 6) {
+                render_icon_texture(m_plate_name_edit_icon.model, m_partplate_list->m_plate_name_edit_hovered_texture);
+                show_tooltip(_u8L("Edit current plate name"));
+			}
+			else
+                render_icon_texture(m_plate_name_edit_icon.model, m_partplate_list->m_plate_name_edit_texture);
+
             if (m_partplate_list->render_plate_settings) {
                 if (hover_id == 5) {
                     if (get_bed_type() == BedType::btDefault && get_print_seq() == PrintSequence::ByDefault && get_first_layer_print_sequence().empty())
@@ -1272,6 +1316,9 @@ void PartPlate::register_raycasters_for_picking(GLCanvas3D &canvas)
     register_model_for_picking(canvas, m_lock_icon, picking_id_component(4));
     if (m_partplate_list->render_plate_settings)
         register_model_for_picking(canvas, m_plate_settings_icon, picking_id_component(5));
+
+    canvas.remove_raycasters_for_picking(SceneRaycaster::EType::Bed, picking_id_component(6));
+    register_model_for_picking(canvas, m_plate_name_edit_icon, picking_id_component(6));
 }
 
 int PartPlate::picking_id_component(int idx) const
@@ -1781,7 +1828,12 @@ void PartPlate::generate_plate_name_texture()
 	poly.contour.append({ scale_(p(0) + PARTPLATE_ICON_GAP_LEFT + offset_x), scale_(p(1) - PARTPLATE_TEXT_OFFSET_Y) });
 
     if (!init_model_from_poly(m_plate_name_icon, poly, GROUND_Z))
-		BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << "Unable to generate geometry buffers for icons\n";
+        BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << "Unable to generate geometry buffers for icons\n";
+
+	auto canvas = this->m_partplate_list->m_plater->get_view3D_canvas3D();
+    canvas->remove_raycasters_for_picking(SceneRaycaster::EType::Bed, picking_id_component(6));
+    calc_vertex_for_plate_name_edit_icon(&m_name_texture, 0, m_plate_name_edit_icon);
+    register_model_for_picking(*canvas, m_plate_name_edit_icon, picking_id_component(6));
 }
 void PartPlate::set_plate_name(const std::string& name) 
 { 
@@ -3124,6 +3176,20 @@ void PartPlateList::generate_icon_textures()
 		}
 	}
 
+	// if (m_plate_name_edit_texture.get_id() == 0)
+	{
+		file_name = path + (m_is_dark ? "plate_name_edit_dark.svg" : "plate_name_edit.svg");
+		if (!m_plate_name_edit_texture.load_from_svg_file(file_name, true, false, false, icon_size)) {
+			BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << boost::format(":load file %1% failed") % file_name;
+		 }
+	}
+    // if (m_plate_name_edit_hovered_texture.get_id() == 0)
+	{
+		file_name = path + (m_is_dark ? "plate_name_edit_hover_dark.svg" : "plate_name_edit_hover.svg");
+		if (!m_plate_name_edit_hovered_texture.load_from_svg_file(file_name, true, false, false, icon_size)) {
+		BOOST_LOG_TRIVIAL(error) << __FUNCTION__ << boost::format(":load file %1% failed") % file_name;
+		}
+	}
 
 	std::string text_str = "01";
 	wxFont* font = find_font(text_str,32);
@@ -3161,7 +3227,8 @@ void PartPlateList::release_icon_textures()
 	m_plate_settings_texture.reset();
 	m_plate_settings_texture.reset();
 	m_plate_settings_hovered_texture.reset();
-
+	m_plate_name_edit_texture.reset();
+	m_plate_name_edit_hovered_texture.reset();
 	for (int i = 0;i < MAX_PLATE_COUNT; i++) {
 		m_idx_textures[i].reset();
 	}

src/slic3r/GUI/PartPlate.hpp

@@ -138,6 +138,7 @@ private:
     PickingModel m_orient_icon;
     PickingModel m_lock_icon;
     PickingModel m_plate_settings_icon;
+    PickingModel m_plate_name_edit_icon;
     GLModel m_plate_idx_icon;
     GLTexture m_texture;
 
@@ -169,6 +170,7 @@ private:
     void calc_gridlines(const ExPolygon& poly, const BoundingBox& pp_bbox);
     void calc_height_limit();
     void calc_vertex_for_number(int index, bool one_number, GLModel &buffer);
+    void calc_vertex_for_plate_name_edit_icon(GLTexture *texture, int index, PickingModel &model);
     void calc_vertex_for_icons(int index, PickingModel &model);
     // void calc_vertex_for_icons_background(int icon_count, GLModel &buffer);
     void render_background(bool force_default_color = false);

src/slic3r/GUI/Project.cpp

@@ -62,6 +62,11 @@ ProjectPanel::ProjectPanel(wxWindow *parent, wxWindowID id, const wxPoint &pos,
 
     Bind(EVT_PROJECT_RELOAD, &ProjectPanel::on_reload, this);
 
+    m_auxiliary = new AuxiliaryPanel(this);
+    m_auxiliary->Hide();
+    main_sizer->Add(m_auxiliary, wxSizerFlags().Expand().Proportion(1));
+    Bind(EVT_AUXILIARY_DONE, [this](wxCommandEvent& e) { update_model_data();});
+
     SetSizer(main_sizer);
     Layout();
     Fit();
@@ -91,38 +96,39 @@ void ProjectPanel::on_reload(wxCommandEvent& evt)
         std::string model_author;
         std::string cover_file;
         std::string description;
-        std::map<std::string, std::vector<json>> files;
-
         std::string p_name;
         std::string p_author;
         std::string p_description;
         std::string p_cover_file;
 
+        std::map<std::string, std::vector<json>> files;
+
         Model model = wxGetApp().plater()->model();
 
-        license = model.model_info->license;
-        model_name = model.model_info->model_name;
-        cover_file = model.model_info->cover_file;
-        description = model.model_info->description;
-        update_type = model.model_info->origin;
+        auto model_info = model.model_info;
+        if (model_info != nullptr) {
+            license     = model_info->license;
+            model_name  = model_info->model_name;
+            cover_file  = model_info->cover_file;
+            description = model_info->description;
+            update_type = model_info->origin;
 
+            try {
+                if (!model_info->copyright.empty()) {
+                    json copy_right = json::parse(model_info->copyright);
 
-        try {
-            if (!model.model_info->copyright.empty()) {
-                json copy_right = json::parse(model.model_info->copyright);
-
-                if (copy_right.is_array()) {
-                    for (auto it = copy_right.begin(); it != copy_right.end(); it++) {
-                        if ((*it).contains("author")) {
-                            model_author = (*it)["author"].get<std::string>();
+                    if (copy_right.is_array()) {
+                        for (auto it = copy_right.begin(); it != copy_right.end(); it++) {
+                            if ((*it).contains("author")) {
+                                model_author = (*it)["author"].get<std::string>();
+                            }
                         }
                     }
                 }
+            } catch (...) {
+                ;
             }
         }
-        catch (...) {
-            ;
-        }
 
         if (model_author.empty() && model.design_info != nullptr)
             model_author = model.design_info->Designer;
@@ -134,12 +140,44 @@ void ProjectPanel::on_reload(wxCommandEvent& evt)
             p_author = model.profile_info->ProfileUserName;
         }
 
-        //file info
+        // file info
         std::string file_path = encode_path(wxGetApp().plater()->model().get_auxiliary_file_temp_path().c_str());
         if (!file_path.empty()) {
             files = Reload(file_path);
+            wxGetApp().CallAfter([this, file_path, files] { m_auxiliary->Reload(file_path, files); });
+        } else {
+            clear_model_info();
+            return;
         }
-        else {
+
+        bool has_content = false;
+        for (const string& v : {
+                 update_type,
+                 license,
+                 model_name,
+                 model_author,
+                 cover_file,
+                 description,
+                 p_name,
+                 p_author,
+                 p_description,
+                 p_cover_file,
+             }) {
+            if (!v.empty()) {
+                has_content = true;
+                break;
+            }
+        }
+        if (!has_content) {
+            for (const auto & file : files) {
+                if (!file.second.empty()) {
+                    has_content = true;
+                    break;
+                }
+            }
+        }
+        if (!has_content) {
+            // Nothing to show, just return
             clear_model_info();
             return;
         }
@@ -180,6 +218,7 @@ void ProjectPanel::on_reload(wxCommandEvent& evt)
 
 void ProjectPanel::msw_rescale() 
 {
+    m_auxiliary->msw_rescale();
 }
 
 void ProjectPanel::on_size(wxSizeEvent &event)
@@ -215,6 +254,9 @@ void ProjectPanel::OnScriptMessage(wxWebViewEvent& evt)
         else if (strCmd == "request_3mf_info") {
             m_web_init_completed = true;
         }
+        else if (strCmd == "edit_project_info") {
+            show_info_editor(true);
+        }
         else if (strCmd == "debug_info") {
             //wxString msg =  j["msg"];
             //OutputDebugString(wxString::Format("Model_Web: msg = %s \r\n", msg));
@@ -227,14 +269,24 @@ void ProjectPanel::OnScriptMessage(wxWebViewEvent& evt)
     }
 }
 
+void ProjectPanel::show_info_editor(bool show)
+{
+    m_browser->Show(!show);
+    m_auxiliary->Show(show);
+    Layout();
+}
+
 void ProjectPanel::update_model_data()
 {
     Model model = wxGetApp().plater()->model();
+    show_info_editor(false);
     clear_model_info();
 
+    m_auxiliary->init_auxiliary();
+
     //basics info
-    if (model.model_info == nullptr)
-        return;
+    //if (model.model_info == nullptr)
+    //    return;
     
     auto event = wxCommandEvent(EVT_PROJECT_RELOAD);
     event.SetEventObject(this);
@@ -258,7 +310,6 @@ std::map<std::string, std::vector<json>> ProjectPanel::Reload(wxString aux_path)
 {
     std::vector<fs::path>                           dir_cache;
     fs::directory_iterator                          iter_end;
-    wxString                                        m_root_dir;
     std::map<std::string, std::vector<json>> m_paths_list;
 
     const static std::array<wxString, 5> s_default_folders = {
@@ -276,24 +327,16 @@ std::map<std::string, std::vector<json>> ProjectPanel::Reload(wxString aux_path)
 
     fs::path new_aux_path(aux_path.ToStdWstring());
 
-    try {
-        fs::remove_all(fs::path(m_root_dir.ToStdWstring()));
-    }
-    catch (...) {
-        BOOST_LOG_TRIVIAL(error) << "Failed  removing the auxiliary directory" << m_root_dir.c_str();
-    }
-
-    m_root_dir = aux_path;
     // Check new path. If not exist, create a new one.
     if (!fs::exists(new_aux_path)) {
         fs::create_directory(new_aux_path);
-        // Create default folders if they are not loaded
-        for (auto folder : s_default_folders) {
-            wxString folder_path = aux_path + "/" + folder;
-            if (fs::exists(folder_path.ToStdWstring())) continue;
-            fs::create_directory(folder_path.ToStdWstring());
-        }
-        return m_paths_list;
+    }
+
+    // Create default folders if they are not loaded
+    for (auto folder : s_default_folders) {
+        wxString folder_path = aux_path + "/" + folder;
+        if (fs::exists(folder_path.ToStdWstring())) continue;
+        fs::create_directory(folder_path.ToStdWstring());
     }
 
     // Load from new path
@@ -320,7 +363,8 @@ std::map<std::string, std::vector<json>> ProjectPanel::Reload(wxString aux_path)
                     wxString file_name = encode_path(file_path.c_str());
                     wxStat(file_name, &strucStat);
                     wxFileOffset filelen = strucStat.st_size;
-     
+
+                    pfile_obj["_filepath"] = file_path;
                     pfile_obj["filename"] = wxGetApp().url_encode(file_path_obj.filename().string().c_str());
                     pfile_obj["size"] = formatBytes((unsigned long)filelen);
 

src/slic3r/GUI/Project.hpp

@@ -32,6 +32,7 @@
 #include "Event.hpp"
 #include "libslic3r/ProjectTask.hpp"
 #include "wxExtensions.hpp"
+#include "Auxiliary.hpp"
 
 #define AUFILE_GREY700 wxColour(107, 107, 107)
 #define AUFILE_GREY500 wxColour(158, 158, 158)
@@ -63,9 +64,12 @@ private:
     bool       m_reload_already = {false};
 
     wxWebView* m_browser = {nullptr};
+    AuxiliaryPanel*   m_auxiliary{nullptr};
     wxString   m_project_home_url;
     wxString   m_root_dir;
     static inline int m_sequence_id = 8000;
+
+    void show_info_editor(bool show);
     
 
 public:
@@ -81,6 +85,7 @@ public:
     void msw_rescale();
     void update_model_data();
     void clear_model_info();
+    void init_auxiliary() { m_auxiliary->init_auxiliary(); }
 
     bool Show(bool show);
     void OnScriptMessage(wxWebViewEvent& evt);

src/slic3r/GUI/Selection.cpp

@@ -1324,8 +1324,11 @@ void Selection::scale_legacy(const Vec3d& scale, TransformationType transformati
                     v.set_instance_offset(m_cache.dragging_center + m * (m_cache.volumes_data[i].get_instance_position() - m_cache.dragging_center));
 
                 v.set_instance_scaling_factor(new_scale);
+                // Restore mirror state
+                v.set_instance_mirror(m_cache.volumes_data[i].get_instance_transform().get_mirror());
             }
             else {
+                const auto mirror = v.get_instance_mirror();
                 if (transformation_type.world() && (std::abs(scale.x() - scale.y()) > EPSILON || std::abs(scale.x() - scale.z()) > EPSILON)) {
                     // Non-uniform scaling. Transform the scaling factors into the local coordinate system.
                     // This is only possible, if the instance rotation is mulitples of ninety degrees.
@@ -1334,16 +1337,24 @@ void Selection::scale_legacy(const Vec3d& scale, TransformationType transformati
                 }
                 else
                     v.set_instance_scaling_factor(scale);
+                // Restore mirror state
+                v.set_instance_mirror(mirror);
             }
 
             // update the instance assemble transform
             ModelObject* object = m_model->objects[v.object_idx()];
             Geometry::Transformation assemble_transform = object->instances[v.instance_idx()]->get_assemble_transformation();
+            const auto               mirror             = assemble_transform.get_mirror();
             assemble_transform.set_scaling_factor(v.get_instance_scaling_factor());
+            assemble_transform.set_mirror(mirror);
             object->instances[v.instance_idx()]->set_assemble_transformation(assemble_transform);
         }
-        else if (is_single_volume() || is_single_modifier())
+        else if (is_single_volume() || is_single_modifier()) {
+            const auto mirror = v.get_volume_transformation().get_mirror();
             v.set_volume_scaling_factor(scale);
+            // Restore mirror state
+            v.set_volume_mirror(mirror);
+        }
         else {
             Transform3d m = Geometry::assemble_transform(Vec3d::Zero(), Vec3d::Zero(), scale);
             if (m_mode == Instance) {
@@ -1354,6 +1365,8 @@ void Selection::scale_legacy(const Vec3d& scale, TransformationType transformati
                     v.set_instance_offset(m_cache.dragging_center + m * (m_cache.volumes_data[i].get_instance_position() - m_cache.dragging_center));
 
                 v.set_instance_scaling_factor(new_scale);
+                // Restore mirror state
+                v.set_instance_mirror(m_cache.volumes_data[i].get_instance_transform().get_mirror());
             }
             else if (m_mode == Volume) {
                 Eigen::Matrix<double, 3, 3, Eigen::DontAlign> new_matrix = (m * m_cache.volumes_data[i].get_volume_scale_matrix()).matrix().block(0, 0, 3, 3);
@@ -1364,6 +1377,8 @@ void Selection::scale_legacy(const Vec3d& scale, TransformationType transformati
                     v.set_volume_offset(m_cache.dragging_center - m_cache.volumes_data[i].get_instance_position() + offset);
                 }
                 v.set_volume_scaling_factor(new_scale);
+                // Restore mirror state
+                v.set_volume_mirror(m_cache.volumes_data[i].get_volume_transform().get_mirror());
             }
         }
     }

src/slic3r/GUI/Tab.cpp

@@ -1988,6 +1988,14 @@ void TabPrint::build()
         optgroup->append_single_option_line("only_one_wall_first_layer");
         optgroup->append_single_option_line("reduce_crossing_wall");
         optgroup->append_single_option_line("max_travel_detour_distance");
+
+        optgroup = page->new_optgroup(L("Small Area Infill Flow Compensation (experimental)"), L"param_advanced");
+        optgroup->append_single_option_line("small_area_infill_flow_compensation");
+        Option option = optgroup->get_option("small_area_infill_flow_compensation_model");
+        option.opt.full_width = true;
+        option.opt.is_code = true;
+        option.opt.height = 15;
+        optgroup->append_single_option_line(option);
         
         optgroup = page->new_optgroup(L("Bridging"), L"param_advanced");
         optgroup->append_single_option_line("bridge_flow");
@@ -2198,6 +2206,7 @@ void TabPrint::build()
         optgroup = page->new_optgroup(L("Special mode"), L"param_special");
         optgroup->append_single_option_line("slicing_mode");
         optgroup->append_single_option_line("print_sequence", "sequent-print");
+        optgroup->append_single_option_line("print_order");
         optgroup->append_single_option_line("spiral_mode", "spiral-vase");
         optgroup->append_single_option_line("spiral_mode_smooth", "spiral-vase#smooth");
         optgroup->append_single_option_line("spiral_mode_max_xy_smoothing", "spiral-vase#max-xy-smoothing");
@@ -2218,7 +2227,7 @@ void TabPrint::build()
         optgroup->append_single_option_line("gcode_comments");
         optgroup->append_single_option_line("gcode_label_objects");
         optgroup->append_single_option_line("exclude_object");
-        Option option = optgroup->get_option("filename_format");
+        option = optgroup->get_option("filename_format");
         // option.opt.full_width = true;
         option.opt.is_code = true;
         option.opt.multiline = true;

src/spline/spline.h

@@ -0,0 +1,951 @@
+/*
+ * spline.h
+ *
+ * simple cubic spline interpolation library without external
+ * dependencies
+ *
+ * ---------------------------------------------------------------------
+ * Copyright (C) 2011, 2014, 2016, 2021 Tino Kluge (ttk448 at gmail.com)
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * ---------------------------------------------------------------------
+ *
+ */
+
+
+#ifndef TK_SPLINE_H
+#define TK_SPLINE_H
+
+#include <cstdio>
+#include <cassert>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#ifdef HAVE_SSTREAM
+#include <sstream>
+#include <string>
+#endif // HAVE_SSTREAM
+
+// not ideal but disable unused-function warnings
+// (we get them because we have implementations in the header file,
+// and this is because we want to be able to quickly separate them
+// into a cpp file if necessary)
+#if !defined(_MSC_VER)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-function"
+#endif
+
+// unnamed namespace only because the implementation is in this
+// header file and we don't want to export symbols to the obj files
+namespace
+{
+
+namespace tk
+{
+
+// spline interpolation
+class spline
+{
+public:
+    // spline types
+    enum spline_type {
+        linear = 10,            // linear interpolation
+        cspline = 30,           // cubic splines (classical C^2)
+        cspline_hermite = 31    // cubic hermite splines (local, only C^1)
+    };
+
+    // boundary condition type for the spline end-points
+    enum bd_type {
+        first_deriv = 1,
+        second_deriv = 2,
+        not_a_knot = 3
+    };
+
+protected:
+    std::vector<double> m_x,m_y;            // x,y coordinates of points
+    // interpolation parameters
+    // f(x) = a_i + b_i*(x-x_i) + c_i*(x-x_i)^2 + d_i*(x-x_i)^3
+    // where a_i = y_i, or else it won't go through grid points
+    std::vector<double> m_b,m_c,m_d;        // spline coefficients
+    double m_c0;                            // for left extrapolation
+    spline_type m_type;
+    bd_type m_left, m_right;
+    double  m_left_value, m_right_value;
+    bool m_made_monotonic;
+    void set_coeffs_from_b();               // calculate c_i, d_i from b_i
+    size_t find_closest(double x) const;    // closest idx so that m_x[idx]<=x
+
+public:
+    // default constructor: set boundary condition to be zero curvature
+    // at both ends, i.e. natural splines
+    spline(): m_type(cspline),
+        m_left(second_deriv), m_right(second_deriv),
+        m_left_value(0.0), m_right_value(0.0), m_made_monotonic(false)
+    {
+        ;
+    }
+    spline(const std::vector<double>& X, const std::vector<double>& Y,
+           spline_type type = cspline,
+           bool make_monotonic = false,
+           bd_type left  = second_deriv, double left_value  = 0.0,
+           bd_type right = second_deriv, double right_value = 0.0
+          ):
+        m_type(type),
+        m_left(left), m_right(right),
+        m_left_value(left_value), m_right_value(right_value),
+        m_made_monotonic(false) // false correct here: make_monotonic() sets it
+    {
+        this->set_points(X,Y,m_type);
+        if(make_monotonic) {
+            this->make_monotonic();
+        }
+    }
+
+
+    // modify boundary conditions: if called it must be before set_points()
+    void set_boundary(bd_type left, double left_value,
+                      bd_type right, double right_value);
+
+    // set all data points (cubic_spline=false means linear interpolation)
+    void set_points(const std::vector<double>& x,
+                    const std::vector<double>& y,
+                    spline_type type=cspline);
+
+    // adjust coefficients so that the spline becomes piecewise monotonic
+    // where possible
+    //   this is done by adjusting slopes at grid points by a non-negative
+    //   factor and this will break C^2
+    //   this can also break boundary conditions if adjustments need to
+    //   be made at the boundary points
+    // returns false if no adjustments have been made, true otherwise
+    bool make_monotonic();
+
+    // evaluates the spline at point x
+    double operator() (double x) const;
+    double deriv(int order, double x) const;
+
+    // solves for all x so that: spline(x) = y
+    std::vector<double> solve(double y, bool ignore_extrapolation=true) const;
+
+    // returns the input data points
+    std::vector<double> get_x() const { return m_x; }
+    std::vector<double> get_y() const { return m_y; }
+    double get_x_min() const { assert(!m_x.empty()); return m_x.front(); }
+    double get_x_max() const { assert(!m_x.empty()); return m_x.back(); }
+
+#ifdef HAVE_SSTREAM
+    // spline info string, i.e. spline type, boundary conditions etc.
+    std::string info() const;
+#endif // HAVE_SSTREAM
+
+};
+
+
+
+namespace internal
+{
+
+// band matrix solver
+class band_matrix
+{
+private:
+    std::vector< std::vector<double> > m_upper;  // upper band
+    std::vector< std::vector<double> > m_lower;  // lower band
+public:
+    band_matrix() {};                             // constructor
+    band_matrix(int dim, int n_u, int n_l);       // constructor
+    ~band_matrix() {};                            // destructor
+    void resize(int dim, int n_u, int n_l);      // init with dim,n_u,n_l
+    int dim() const;                             // matrix dimension
+    int num_upper() const
+    {
+        return (int)m_upper.size()-1;
+    }
+    int num_lower() const
+    {
+        return (int)m_lower.size()-1;
+    }
+    // access operator
+    double & operator () (int i, int j);            // write
+    double   operator () (int i, int j) const;      // read
+    // we can store an additional diagonal (in m_lower)
+    double& saved_diag(int i);
+    double  saved_diag(int i) const;
+    void lu_decompose();
+    std::vector<double> r_solve(const std::vector<double>& b) const;
+    std::vector<double> l_solve(const std::vector<double>& b) const;
+    std::vector<double> lu_solve(const std::vector<double>& b,
+                                 bool is_lu_decomposed=false);
+
+};
+
+double get_eps();
+
+std::vector<double> solve_cubic(double a, double b, double c, double d,
+                                int newton_iter=0);
+
+} // namespace internal
+
+
+
+
+// ---------------------------------------------------------------------
+// implementation part, which could be separated into a cpp file
+// ---------------------------------------------------------------------
+
+// spline implementation
+// -----------------------
+
+void spline::set_boundary(spline::bd_type left, double left_value,
+                          spline::bd_type right, double right_value)
+{
+    assert(m_x.size()==0);          // set_points() must not have happened yet
+    m_left=left;
+    m_right=right;
+    m_left_value=left_value;
+    m_right_value=right_value;
+}
+
+
+void spline::set_coeffs_from_b()
+{
+    assert(m_x.size()==m_y.size());
+    assert(m_x.size()==m_b.size());
+    assert(m_x.size()>2);
+    size_t n=m_b.size();
+    if(m_c.size()!=n)
+        m_c.resize(n);
+    if(m_d.size()!=n)
+        m_d.resize(n);
+
+    for(size_t i=0; i<n-1; i++) {
+        const double h  = m_x[i+1]-m_x[i];
+        // from continuity and differentiability condition
+        m_c[i] = ( 3.0*(m_y[i+1]-m_y[i])/h - (2.0*m_b[i]+m_b[i+1]) ) / h;
+        // from differentiability condition
+        m_d[i] = ( (m_b[i+1]-m_b[i])/(3.0*h) - 2.0/3.0*m_c[i] ) / h;
+    }
+
+    // for left extrapolation coefficients
+    m_c0 = (m_left==first_deriv) ? 0.0 : m_c[0];
+}
+
+void spline::set_points(const std::vector<double>& x,
+                        const std::vector<double>& y,
+                        spline_type type)
+{
+    assert(x.size()==y.size());
+    assert(x.size()>=3);
+    // not-a-knot with 3 points has many solutions
+    if(m_left==not_a_knot || m_right==not_a_knot)
+        assert(x.size()>=4);
+    m_type=type;
+    m_made_monotonic=false;
+    m_x=x;
+    m_y=y;
+    int n = (int) x.size();
+    // check strict monotonicity of input vector x
+    for(int i=0; i<n-1; i++) {
+        assert(m_x[i]<m_x[i+1]);
+    }
+
+
+    if(type==linear) {
+        // linear interpolation
+        m_d.resize(n);
+        m_c.resize(n);
+        m_b.resize(n);
+        for(int i=0; i<n-1; i++) {
+            m_d[i]=0.0;
+            m_c[i]=0.0;
+            m_b[i]=(m_y[i+1]-m_y[i])/(m_x[i+1]-m_x[i]);
+        }
+        // ignore boundary conditions, set slope equal to the last segment
+        m_b[n-1]=m_b[n-2];
+        m_c[n-1]=0.0;
+        m_d[n-1]=0.0;
+    } else if(type==cspline) {
+        // classical cubic splines which are C^2 (twice cont differentiable)
+        // this requires solving an equation system
+
+        // setting up the matrix and right hand side of the equation system
+        // for the parameters b[]
+        int n_upper = (m_left  == spline::not_a_knot) ? 2 : 1;
+        int n_lower = (m_right == spline::not_a_knot) ? 2 : 1;
+        internal::band_matrix A(n,n_upper,n_lower);
+        std::vector<double>  rhs(n);
+        for(int i=1; i<n-1; i++) {
+            A(i,i-1)=1.0/3.0*(x[i]-x[i-1]);
+            A(i,i)=2.0/3.0*(x[i+1]-x[i-1]);
+            A(i,i+1)=1.0/3.0*(x[i+1]-x[i]);
+            rhs[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
+        }
+        // boundary conditions
+        if(m_left == spline::second_deriv) {
+            // 2*c[0] = f''
+            A(0,0)=2.0;
+            A(0,1)=0.0;
+            rhs[0]=m_left_value;
+        } else if(m_left == spline::first_deriv) {
+            // b[0] = f', needs to be re-expressed in terms of c:
+            // (2c[0]+c[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
+            A(0,0)=2.0*(x[1]-x[0]);
+            A(0,1)=1.0*(x[1]-x[0]);
+            rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value);
+        } else if(m_left == spline::not_a_knot) {
+            // f'''(x[1]) exists, i.e. d[0]=d[1], or re-expressed in c:
+            // -h1*c[0] + (h0+h1)*c[1] - h0*c[2] = 0
+            A(0,0) = -(x[2]-x[1]);
+            A(0,1) = x[2]-x[0];
+            A(0,2) = -(x[1]-x[0]);
+            rhs[0] = 0.0;
+        } else {
+            assert(false);
+        }
+        if(m_right == spline::second_deriv) {
+            // 2*c[n-1] = f''
+            A(n-1,n-1)=2.0;
+            A(n-1,n-2)=0.0;
+            rhs[n-1]=m_right_value;
+        } else if(m_right == spline::first_deriv) {
+            // b[n-1] = f', needs to be re-expressed in terms of c:
+            // (c[n-2]+2c[n-1])(x[n-1]-x[n-2])
+            // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
+            A(n-1,n-1)=2.0*(x[n-1]-x[n-2]);
+            A(n-1,n-2)=1.0*(x[n-1]-x[n-2]);
+            rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
+        } else if(m_right == spline::not_a_knot) {
+            // f'''(x[n-2]) exists, i.e. d[n-3]=d[n-2], or re-expressed in c:
+            // -h_{n-2}*c[n-3] + (h_{n-3}+h_{n-2})*c[n-2] - h_{n-3}*c[n-1] = 0
+            A(n-1,n-3) = -(x[n-1]-x[n-2]);
+            A(n-1,n-2) = x[n-1]-x[n-3];
+            A(n-1,n-1) = -(x[n-2]-x[n-3]);
+            rhs[0] = 0.0;
+        } else {
+            assert(false);
+        }
+
+        // solve the equation system to obtain the parameters c[]
+        m_c=A.lu_solve(rhs);
+
+        // calculate parameters b[] and d[] based on c[]
+        m_d.resize(n);
+        m_b.resize(n);
+        for(int i=0; i<n-1; i++) {
+            m_d[i]=1.0/3.0*(m_c[i+1]-m_c[i])/(x[i+1]-x[i]);
+            m_b[i]=(y[i+1]-y[i])/(x[i+1]-x[i])
+                   - 1.0/3.0*(2.0*m_c[i]+m_c[i+1])*(x[i+1]-x[i]);
+        }
+        // for the right extrapolation coefficients (zero cubic term)
+        // f_{n-1}(x) = y_{n-1} + b*(x-x_{n-1}) + c*(x-x_{n-1})^2
+        double h=x[n-1]-x[n-2];
+        // m_c[n-1] is determined by the boundary condition
+        m_d[n-1]=0.0;
+        m_b[n-1]=3.0*m_d[n-2]*h*h+2.0*m_c[n-2]*h+m_b[n-2];   // = f'_{n-2}(x_{n-1})
+        if(m_right==first_deriv)
+            m_c[n-1]=0.0;   // force linear extrapolation
+
+    } else if(type==cspline_hermite) {
+        // hermite cubic splines which are C^1 (cont. differentiable)
+        // and derivatives are specified on each grid point
+        // (here we use 3-point finite differences)
+        m_b.resize(n);
+        m_c.resize(n);
+        m_d.resize(n);
+        // set b to match 1st order derivative finite difference
+        for(int i=1; i<n-1; i++) {
+            const double h  = m_x[i+1]-m_x[i];
+            const double hl = m_x[i]-m_x[i-1];
+            m_b[i] = -h/(hl*(hl+h))*m_y[i-1] + (h-hl)/(hl*h)*m_y[i]
+                     +  hl/(h*(hl+h))*m_y[i+1];
+        }
+        // boundary conditions determine b[0] and b[n-1]
+        if(m_left==first_deriv) {
+            m_b[0]=m_left_value;
+        } else if(m_left==second_deriv) {
+            const double h = m_x[1]-m_x[0];
+            m_b[0]=0.5*(-m_b[1]-0.5*m_left_value*h+3.0*(m_y[1]-m_y[0])/h);
+        } else if(m_left == not_a_knot) {
+            // f''' continuous at x[1]
+            const double h0 = m_x[1]-m_x[0];
+            const double h1 = m_x[2]-m_x[1];
+            m_b[0]= -m_b[1] + 2.0*(m_y[1]-m_y[0])/h0
+                    + h0*h0/(h1*h1)*(m_b[1]+m_b[2]-2.0*(m_y[2]-m_y[1])/h1);
+        } else {
+            assert(false);
+        }
+        if(m_right==first_deriv) {
+            m_b[n-1]=m_right_value;
+            m_c[n-1]=0.0;
+        } else if(m_right==second_deriv) {
+            const double h = m_x[n-1]-m_x[n-2];
+            m_b[n-1]=0.5*(-m_b[n-2]+0.5*m_right_value*h+3.0*(m_y[n-1]-m_y[n-2])/h);
+            m_c[n-1]=0.5*m_right_value;
+        } else if(m_right == not_a_knot) {
+            // f''' continuous at x[n-2]
+            const double h0 = m_x[n-2]-m_x[n-3];
+            const double h1 = m_x[n-1]-m_x[n-2];
+            m_b[n-1]= -m_b[n-2] + 2.0*(m_y[n-1]-m_y[n-2])/h1 + h1*h1/(h0*h0)
+                      *(m_b[n-3]+m_b[n-2]-2.0*(m_y[n-2]-m_y[n-3])/h0);
+            // f'' continuous at x[n-1]: c[n-1] = 3*d[n-2]*h[n-2] + c[n-1]
+            m_c[n-1]=(m_b[n-2]+2.0*m_b[n-1])/h1-3.0*(m_y[n-1]-m_y[n-2])/(h1*h1);
+        } else {
+            assert(false);
+        }
+        m_d[n-1]=0.0;
+
+        // parameters c and d are determined by continuity and differentiability
+        set_coeffs_from_b();
+
+    } else {
+        assert(false);
+    }
+
+    // for left extrapolation coefficients
+    m_c0 = (m_left==first_deriv) ? 0.0 : m_c[0];
+}
+
+bool spline::make_monotonic()
+{
+    assert(m_x.size()==m_y.size());
+    assert(m_x.size()==m_b.size());
+    assert(m_x.size()>2);
+    bool modified = false;
+    const int n=(int)m_x.size();
+    // make sure: input data monotonic increasing --> b_i>=0
+    //            input data monotonic decreasing --> b_i<=0
+    for(int i=0; i<n; i++) {
+        int im1 = std::max(i-1, 0);
+        int ip1 = std::min(i+1, n-1);
+        if( ((m_y[im1]<=m_y[i]) && (m_y[i]<=m_y[ip1]) && m_b[i]<0.0) ||
+            ((m_y[im1]>=m_y[i]) && (m_y[i]>=m_y[ip1]) && m_b[i]>0.0) ) {
+            modified=true;
+            m_b[i]=0.0;
+        }
+    }
+    // if input data is monotonic (b[i], b[i+1], avg have all the same sign)
+    // ensure a sufficient criteria for monotonicity is satisfied:
+    //     sqrt(b[i]^2+b[i+1]^2) <= 3 |avg|, with avg=(y[i+1]-y[i])/h,
+    for(int i=0; i<n-1; i++) {
+        double h = m_x[i+1]-m_x[i];
+        double avg = (m_y[i+1]-m_y[i])/h;
+        if( avg==0.0 && (m_b[i]!=0.0 || m_b[i+1]!=0.0) ) {
+            modified=true;
+            m_b[i]=0.0;
+            m_b[i+1]=0.0;
+        } else if( (m_b[i]>=0.0 && m_b[i+1]>=0.0 && avg>0.0) ||
+                   (m_b[i]<=0.0 && m_b[i+1]<=0.0 && avg<0.0) ) {
+            // input data is monotonic
+            double r = sqrt(m_b[i]*m_b[i]+m_b[i+1]*m_b[i+1])/std::fabs(avg);
+            if(r>3.0) {
+                // sufficient criteria for monotonicity: r<=3
+                // adjust b[i] and b[i+1]
+                modified=true;
+                m_b[i]   *= (3.0/r);
+                m_b[i+1] *= (3.0/r);
+            }
+        }
+    }
+
+    if(modified==true) {
+        set_coeffs_from_b();
+        m_made_monotonic=true;
+    }
+
+    return modified;
+}
+
+// return the closest idx so that m_x[idx] <= x (return 0 if x<m_x[0])
+size_t spline::find_closest(double x) const
+{
+    std::vector<double>::const_iterator it;
+    it=std::upper_bound(m_x.begin(),m_x.end(),x);       // *it > x
+    size_t idx = std::max( int(it-m_x.begin())-1, 0);   // m_x[idx] <= x
+    return idx;
+}
+
+double spline::operator() (double x) const
+{
+    // polynomial evaluation using Horner's scheme
+    // TODO: consider more numerically accurate algorithms, e.g.:
+    //   - Clenshaw
+    //   - Even-Odd method by A.C.R. Newbery
+    //   - Compensated Horner Scheme
+    size_t n=m_x.size();
+    size_t idx=find_closest(x);
+
+    double h=x-m_x[idx];
+    double interpol;
+    if(x<m_x[0]) {
+        // extrapolation to the left
+        interpol=(m_c0*h + m_b[0])*h + m_y[0];
+    } else if(x>m_x[n-1]) {
+        // extrapolation to the right
+        interpol=(m_c[n-1]*h + m_b[n-1])*h + m_y[n-1];
+    } else {
+        // interpolation
+        interpol=((m_d[idx]*h + m_c[idx])*h + m_b[idx])*h + m_y[idx];
+    }
+    return interpol;
+}
+
+double spline::deriv(int order, double x) const
+{
+    assert(order>0);
+    size_t n=m_x.size();
+    size_t idx = find_closest(x);
+
+    double h=x-m_x[idx];
+    double interpol;
+    if(x<m_x[0]) {
+        // extrapolation to the left
+        switch(order) {
+        case 1:
+            interpol=2.0*m_c0*h + m_b[0];
+            break;
+        case 2:
+            interpol=2.0*m_c0;
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    } else if(x>m_x[n-1]) {
+        // extrapolation to the right
+        switch(order) {
+        case 1:
+            interpol=2.0*m_c[n-1]*h + m_b[n-1];
+            break;
+        case 2:
+            interpol=2.0*m_c[n-1];
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    } else {
+        // interpolation
+        switch(order) {
+        case 1:
+            interpol=(3.0*m_d[idx]*h + 2.0*m_c[idx])*h + m_b[idx];
+            break;
+        case 2:
+            interpol=6.0*m_d[idx]*h + 2.0*m_c[idx];
+            break;
+        case 3:
+            interpol=6.0*m_d[idx];
+            break;
+        default:
+            interpol=0.0;
+            break;
+        }
+    }
+    return interpol;
+}
+
+std::vector<double> spline::solve(double y, bool ignore_extrapolation) const
+{
+    std::vector<double> x;          // roots for the entire spline
+    std::vector<double> root;       // roots for each piecewise cubic
+    const size_t n=m_x.size();
+
+    // left extrapolation
+    if(ignore_extrapolation==false) {
+        root = internal::solve_cubic(m_y[0]-y,m_b[0],m_c0,0.0,1);
+        for(size_t j=0; j<root.size(); j++) {
+            if(root[j]<0.0) {
+                x.push_back(m_x[0]+root[j]);
+            }
+        }
+    }
+
+    // brute force check if piecewise cubic has roots in their resp. segment
+    // TODO: make more efficient
+    for(size_t i=0; i<n-1; i++) {
+        root = internal::solve_cubic(m_y[i]-y,m_b[i],m_c[i],m_d[i],1);
+        for(size_t j=0; j<root.size(); j++) {
+            double h = (i>0) ? (m_x[i]-m_x[i-1]) : 0.0;
+            double eps = internal::get_eps()*512.0*std::min(h,1.0);
+            if( (-eps<=root[j]) && (root[j]<m_x[i+1]-m_x[i]) ) {
+                double new_root = m_x[i]+root[j];
+                if(x.size()>0 && x.back()+eps > new_root) {
+                    x.back()=new_root;      // avoid spurious duplicate roots
+                } else {
+                    x.push_back(new_root);
+                }
+            }
+        }
+    }
+
+    // right extrapolation
+    if(ignore_extrapolation==false) {
+        root = internal::solve_cubic(m_y[n-1]-y,m_b[n-1],m_c[n-1],0.0,1);
+        for(size_t j=0; j<root.size(); j++) {
+            if(0.0<=root[j]) {
+                x.push_back(m_x[n-1]+root[j]);
+            }
+        }
+    }
+
+    return x;
+};
+
+
+#ifdef HAVE_SSTREAM
+std::string spline::info() const
+{
+    std::stringstream ss;
+    ss << "type " << m_type << ", left boundary deriv " << m_left << " = ";
+    ss << m_left_value << ", right boundary deriv " << m_right << " = ";
+    ss << m_right_value << std::endl;
+    if(m_made_monotonic) {
+        ss << "(spline has been adjusted for piece-wise monotonicity)";
+    }
+    return ss.str();
+}
+#endif // HAVE_SSTREAM
+
+
+namespace internal
+{
+
+// band_matrix implementation
+// -------------------------
+
+band_matrix::band_matrix(int dim, int n_u, int n_l)
+{
+    resize(dim, n_u, n_l);
+}
+void band_matrix::resize(int dim, int n_u, int n_l)
+{
+    assert(dim>0);
+    assert(n_u>=0);
+    assert(n_l>=0);
+    m_upper.resize(n_u+1);
+    m_lower.resize(n_l+1);
+    for(size_t i=0; i<m_upper.size(); i++) {
+        m_upper[i].resize(dim);
+    }
+    for(size_t i=0; i<m_lower.size(); i++) {
+        m_lower[i].resize(dim);
+    }
+}
+int band_matrix::dim() const
+{
+    if(m_upper.size()>0) {
+        return (int)m_upper[0].size();
+    } else {
+        return 0;
+    }
+}
+
+
+// defines the new operator (), so that we can access the elements
+// by A(i,j), index going from i=0,...,dim()-1
+double & band_matrix::operator () (int i, int j)
+{
+    int k=j-i;       // what band is the entry
+    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
+    assert( (-num_lower()<=k) && (k<=num_upper()) );
+    // k=0 -> diagonal, k<0 lower left part, k>0 upper right part
+    if(k>=0)    return m_upper[k][i];
+    else        return m_lower[-k][i];
+}
+double band_matrix::operator () (int i, int j) const
+{
+    int k=j-i;       // what band is the entry
+    assert( (i>=0) && (i<dim()) && (j>=0) && (j<dim()) );
+    assert( (-num_lower()<=k) && (k<=num_upper()) );
+    // k=0 -> diagonal, k<0 lower left part, k>0 upper right part
+    if(k>=0)    return m_upper[k][i];
+    else        return m_lower[-k][i];
+}
+// second diag (used in LU decomposition), saved in m_lower
+double band_matrix::saved_diag(int i) const
+{
+    assert( (i>=0) && (i<dim()) );
+    return m_lower[0][i];
+}
+double & band_matrix::saved_diag(int i)
+{
+    assert( (i>=0) && (i<dim()) );
+    return m_lower[0][i];
+}
+
+// LR-Decomposition of a band matrix
+void band_matrix::lu_decompose()
+{
+    int  i_max,j_max;
+    int  j_min;
+    double x;
+
+    // preconditioning
+    // normalize column i so that a_ii=1
+    for(int i=0; i<this->dim(); i++) {
+        assert(this->operator()(i,i)!=0.0);
+        this->saved_diag(i)=1.0/this->operator()(i,i);
+        j_min=std::max(0,i-this->num_lower());
+        j_max=std::min(this->dim()-1,i+this->num_upper());
+        for(int j=j_min; j<=j_max; j++) {
+            this->operator()(i,j) *= this->saved_diag(i);
+        }
+        this->operator()(i,i)=1.0;          // prevents rounding errors
+    }
+
+    // Gauss LR-Decomposition
+    for(int k=0; k<this->dim(); k++) {
+        i_max=std::min(this->dim()-1,k+this->num_lower());  // num_lower not a mistake!
+        for(int i=k+1; i<=i_max; i++) {
+            assert(this->operator()(k,k)!=0.0);
+            x=-this->operator()(i,k)/this->operator()(k,k);
+            this->operator()(i,k)=-x;                         // assembly part of L
+            j_max=std::min(this->dim()-1,k+this->num_upper());
+            for(int j=k+1; j<=j_max; j++) {
+                // assembly part of R
+                this->operator()(i,j)=this->operator()(i,j)+x*this->operator()(k,j);
+            }
+        }
+    }
+}
+// solves Ly=b
+std::vector<double> band_matrix::l_solve(const std::vector<double>& b) const
+{
+    assert( this->dim()==(int)b.size() );
+    std::vector<double> x(this->dim());
+    int j_start;
+    double sum;
+    for(int i=0; i<this->dim(); i++) {
+        sum=0;
+        j_start=std::max(0,i-this->num_lower());
+        for(int j=j_start; j<i; j++) sum += this->operator()(i,j)*x[j];
+        x[i]=(b[i]*this->saved_diag(i)) - sum;
+    }
+    return x;
+}
+// solves Rx=y
+std::vector<double> band_matrix::r_solve(const std::vector<double>& b) const
+{
+    assert( this->dim()==(int)b.size() );
+    std::vector<double> x(this->dim());
+    int j_stop;
+    double sum;
+    for(int i=this->dim()-1; i>=0; i--) {
+        sum=0;
+        j_stop=std::min(this->dim()-1,i+this->num_upper());
+        for(int j=i+1; j<=j_stop; j++) sum += this->operator()(i,j)*x[j];
+        x[i]=( b[i] - sum ) / this->operator()(i,i);
+    }
+    return x;
+}
+
+std::vector<double> band_matrix::lu_solve(const std::vector<double>& b,
+        bool is_lu_decomposed)
+{
+    assert( this->dim()==(int)b.size() );
+    std::vector<double>  x,y;
+    if(is_lu_decomposed==false) {
+        this->lu_decompose();
+    }
+    y=this->l_solve(b);
+    x=this->r_solve(y);
+    return x;
+}
+
+// machine precision of a double, i.e. the successor of 1 is 1+eps
+double get_eps()
+{
+    //return std::numeric_limits<double>::epsilon();    // __DBL_EPSILON__
+    return 2.2204460492503131e-16;                      // 2^-52
+}
+
+// solutions for a + b*x = 0
+std::vector<double> solve_linear(double a, double b)
+{
+    std::vector<double> x;      // roots
+    if(b==0.0) {
+        if(a==0.0) {
+            // 0*x = 0
+            x.resize(1);
+            x[0] = 0.0;   // any x solves it but we need to pick one
+            return x;
+        } else {
+            // 0*x + ... = 0, no solution
+            return x;
+        }
+    } else {
+        x.resize(1);
+        x[0] = -a/b;
+        return x;
+    }
+}
+
+// solutions for a + b*x + c*x^2 = 0
+std::vector<double> solve_quadratic(double a, double b, double c,
+                                    int newton_iter=0)
+{
+    if(c==0.0) {
+        return solve_linear(a,b);
+    }
+    // rescale so that we solve x^2 + 2p x + q = (x+p)^2 + q - p^2 = 0
+    double p=0.5*b/c;
+    double q=a/c;
+    double discr = p*p-q;
+    const double eps=0.5*internal::get_eps();
+    double discr_err = (6.0*(p*p)+3.0*fabs(q)+fabs(discr))*eps;
+
+    std::vector<double> x;      // roots
+    if(fabs(discr)<=discr_err) {
+        // discriminant is zero --> one root
+        x.resize(1);
+        x[0] = -p;
+    } else if(discr<0) {
+        // no root
+    } else {
+        // two roots
+        x.resize(2);
+        x[0] = -p - sqrt(discr);
+        x[1] = -p + sqrt(discr);
+    }
+
+    // improve solution via newton steps
+    for(size_t i=0; i<x.size(); i++) {
+        for(int k=0; k<newton_iter; k++) {
+            double f  = (c*x[i] + b)*x[i] + a;
+            double f1 = 2.0*c*x[i] + b;
+            // only adjust if slope is large enough
+            if(fabs(f1)>1e-8) {
+                x[i] -= f/f1;
+            }
+        }
+    }
+
+    return x;
+}
+
+// solutions for the cubic equation: a + b*x +c*x^2 + d*x^3 = 0
+// this is a naive implementation of the analytic solution without
+// optimisation for speed or numerical accuracy
+// newton_iter: number of newton iterations to improve analytical solution
+// see also
+//   gsl: gsl_poly_solve_cubic() in solve_cubic.c
+//   octave: roots.m - via eigenvalues of the Frobenius companion matrix
+std::vector<double> solve_cubic(double a, double b, double c, double d,
+                                int newton_iter)
+{
+    if(d==0.0) {
+        return solve_quadratic(a,b,c,newton_iter);
+    }
+
+    // convert to normalised form: a + bx + cx^2 + x^3 = 0
+    if(d!=1.0) {
+        a/=d;
+        b/=d;
+        c/=d;
+    }
+
+    // convert to depressed cubic: z^3 - 3pz - 2q = 0
+    // via substitution: z = x + c/3
+    std::vector<double> z;              // roots of the depressed cubic
+    double p = -(1.0/3.0)*b + (1.0/9.0)*(c*c);
+    double r = 2.0*(c*c)-9.0*b;
+    double q = -0.5*a - (1.0/54.0)*(c*r);
+    double discr=p*p*p-q*q;             // discriminant
+    // calculating numerical round-off errors with assumptions:
+    //  - each operation is precise but each intermediate result x
+    //    when stored has max error of x*eps
+    //  - only multiplication with a power of 2 introduces no new error
+    //  - a,b,c,d and some fractions (e.g. 1/3) have rounding errors eps
+    //  - p_err << |p|, q_err << |q|, ... (this is violated in rare cases)
+    // would be more elegant to use boost::numeric::interval<double>
+    const double eps = internal::get_eps();
+    double p_err = eps*((3.0/3.0)*fabs(b)+(4.0/9.0)*(c*c)+fabs(p));
+    double r_err = eps*(6.0*(c*c)+18.0*fabs(b)+fabs(r));
+    double q_err = 0.5*fabs(a)*eps + (1.0/54.0)*fabs(c)*(r_err+fabs(r)*3.0*eps)
+                   + fabs(q)*eps;
+    double discr_err = (p*p) * (3.0*p_err + fabs(p)*2.0*eps)
+                       + fabs(q) * (2.0*q_err + fabs(q)*eps) + fabs(discr)*eps;
+
+    // depending on the discriminant we get different solutions
+    if(fabs(discr)<=discr_err) {
+        // discriminant zero: one or two real roots
+        if(fabs(p)<=p_err) {
+            // p and q are zero: single root
+            z.resize(1);
+            z[0] = 0.0;             // triple root
+        } else {
+            z.resize(2);
+            z[0] = 2.0*q/p;         // single root
+            z[1] = -0.5*z[0];       // double root
+        }
+    } else if(discr>0) {
+        // three real roots: via trigonometric solution
+        z.resize(3);
+        double ac = (1.0/3.0) * acos( q/(p*sqrt(p)) );
+        double sq = 2.0*sqrt(p);
+        z[0] = sq * cos(ac);
+        z[1] = sq * cos(ac-2.0*M_PI/3.0);
+        z[2] = sq * cos(ac-4.0*M_PI/3.0);
+    } else if (discr<0.0) {
+        // single real root: via Cardano's fromula
+        z.resize(1);
+        double sgnq = (q >= 0 ? 1 : -1);
+        double basis = fabs(q) + sqrt(-discr);
+        double C = sgnq * pow(basis, 1.0/3.0); // c++11 has std::cbrt()
+        z[0] = C + p/C;
+    }
+    for(size_t i=0; i<z.size(); i++) {
+        // convert depressed cubic roots to original cubic: x = z - c/3
+        z[i] -= (1.0/3.0)*c;
+        // improve solution via newton steps
+        for(int k=0; k<newton_iter; k++) {
+            double f  = ((z[i] + c)*z[i] + b)*z[i] + a;
+            double f1 = (3.0*z[i] + 2.0*c)*z[i] + b;
+            // only adjust if slope is large enough
+            if(fabs(f1)>1e-8) {
+                z[i] -= f/f1;
+            }
+        }
+    }
+    // ensure if a=0 we get exactly x=0 as root
+    // TODO: remove this fudge
+    if(a==0.0) {
+        assert(z.size()>0);     // cubic should always have at least one root
+        double xmin=fabs(z[0]);
+        size_t imin=0;
+        for(size_t i=1; i<z.size(); i++) {
+            if(xmin>fabs(z[i])) {
+                xmin=fabs(z[i]);
+                imin=i;
+            }
+        }
+        z[imin]=0.0;        // replace the smallest absolute value with 0
+    }
+    std::sort(z.begin(), z.end());
+    return z;
+}
+
+
+} // namespace internal
+
+
+} // namespace tk
+
+
+} // namespace
+
+#if !defined(_MSC_VER)
+#pragma GCC diagnostic pop
+#endif
+
+#endif /* TK_SPLINE_H */