Merge branch 'adaptive-slicing' into merge-adaptive

# Conflicts:
#	xs/src/libslic3r/Model.cpp
#	xs/src/libslic3r/Model.hpp

lib/Slic3r.pm

@@ -229,6 +229,7 @@ sub thread_cleanup {
     *Slic3r::Geometry::BoundingBoxf::DESTROY = sub {};
     *Slic3r::Geometry::BoundingBoxf3::DESTROY = sub {};
     *Slic3r::Layer::PerimeterGenerator::DESTROY = sub {};
+    *Slic3r::LayerHeightSpline::DESTROY     = sub {};
     *Slic3r::Line::DESTROY                  = sub {};
     *Slic3r::Linef3::DESTROY                = sub {};
     *Slic3r::Model::DESTROY                 = sub {};

lib/Slic3r/GUI.pm

@@ -42,6 +42,7 @@ use Slic3r::GUI::Plater::ObjectCutDialog;
 use Slic3r::GUI::Plater::ObjectSettingsDialog;
 use Slic3r::GUI::Plater::LambdaObjectDialog;
 use Slic3r::GUI::Plater::OverrideSettingsPanel;
+use Slic3r::GUI::Plater::SplineControl;
 use Slic3r::GUI::Preferences;
 use Slic3r::GUI::ProgressStatusBar;
 use Slic3r::GUI::Projector;

lib/Slic3r/GUI/OptionsGroup/Field.pm

@@ -580,6 +580,17 @@ sub get_value {
     return $self->slider->GetValue/$self->scale;
 }
 
+# Update internal scaling
+sub set_scale {
+    my ($self, $scale) = @_;
+    $self->disable_change_event(1);
+    my $current_value = $self->get_value;
+    $self->slider->SetRange($self->slider->GetMin / $self->scale * $scale, $self->slider->GetMax / $self->scale * $scale);
+    $self->scale($scale);
+    $self->set_value($current_value);
+    $self->disable_change_event(0);
+}
+
 sub _update_textctrl {
     my ($self) = @_;
     

lib/Slic3r/GUI/Plater.pm

@@ -28,7 +28,7 @@ use threads::shared qw(shared_clone);
 use Wx qw(:button :cursor :dialog :filedialog :keycode :icon :font :id :misc 
     :panel :sizer :toolbar :window wxTheApp :notebook :combobox);
 use Wx::Event qw(EVT_BUTTON EVT_COMMAND EVT_KEY_DOWN EVT_MOUSE_EVENTS EVT_PAINT EVT_TOOL 
-    EVT_CHOICE EVT_COMBOBOX EVT_TIMER EVT_NOTEBOOK_PAGE_CHANGED EVT_LEFT_UP);
+    EVT_CHOICE EVT_COMBOBOX EVT_TIMER EVT_NOTEBOOK_PAGE_CHANGED EVT_LEFT_UP EVT_CLOSE);
 use base qw(Wx::Panel Class::Accessor);
 
 __PACKAGE__->mk_accessors(qw(presets));
@@ -46,6 +46,7 @@ use constant TB_45CCW   => &Wx::NewId;
 use constant TB_SCALE   => &Wx::NewId;
 use constant TB_SPLIT   => &Wx::NewId;
 use constant TB_CUT     => &Wx::NewId;
+use constant TB_LAYERS  => &Wx::NewId;
 use constant TB_SETTINGS => &Wx::NewId;
 
 # package variables to avoid passing lexicals to threads
@@ -195,6 +196,7 @@ sub new {
         $self->{htoolbar}->AddTool(TB_CUT, "Cut…", Wx::Bitmap->new($Slic3r::var->("package.png"), wxBITMAP_TYPE_PNG), '');
         $self->{htoolbar}->AddSeparator;
         $self->{htoolbar}->AddTool(TB_SETTINGS, "Settings…", Wx::Bitmap->new($Slic3r::var->("cog.png"), wxBITMAP_TYPE_PNG), '');
+        $self->{htoolbar}->AddTool(TB_LAYERS, "Layer heights…", Wx::Bitmap->new($Slic3r::var->("variable_layer_height.png"), wxBITMAP_TYPE_PNG), '');
     } else {
         my %tbar_buttons = (
             add             => "Add…",
@@ -208,10 +210,11 @@ sub new {
             changescale     => "Scale…",
             split           => "Split",
             cut             => "Cut…",
+            layers          => "Layer heights…",
             settings        => "Settings…",
         );
         $self->{btoolbar} = Wx::BoxSizer->new(wxHORIZONTAL);
-        for (qw(add remove reset arrange increase decrease rotate45ccw rotate45cw changescale split cut settings)) {
+        for (qw(add remove reset arrange increase decrease rotate45ccw rotate45cw changescale split cut layers settings)) {
             $self->{"btn_$_"} = Wx::Button->new($self, -1, $tbar_buttons{$_}, wxDefaultPosition, wxDefaultSize, wxBU_EXACTFIT);
             $self->{btoolbar}->Add($self->{"btn_$_"});
         }
@@ -245,6 +248,7 @@ sub new {
             changescale     arrow_out.png
             split           shape_ungroup.png
             cut             package.png
+            layers          cog.png
             settings        cog.png
         );
         for (grep $self->{"btn_$_"}, keys %icons) {
@@ -281,6 +285,7 @@ sub new {
         EVT_TOOL($self, TB_SCALE, sub { $self->changescale(undef); });
         EVT_TOOL($self, TB_SPLIT, sub { $self->split_object; });
         EVT_TOOL($self, TB_CUT, sub { $_[0]->object_cut_dialog });
+        EVT_TOOL($self, TB_LAYERS, sub { $_[0]->object_layers_dialog });
         EVT_TOOL($self, TB_SETTINGS, sub { $_[0]->object_settings_dialog });
     } else {
         EVT_BUTTON($self, $self->{btn_add}, sub { $self->add; });
@@ -294,6 +299,7 @@ sub new {
         EVT_BUTTON($self, $self->{btn_changescale}, sub { $self->changescale(undef); });
         EVT_BUTTON($self, $self->{btn_split}, sub { $self->split_object; });
         EVT_BUTTON($self, $self->{btn_cut}, sub { $_[0]->object_cut_dialog });
+        EVT_BUTTON($self, $self->{btn_layers}, sub { $_[0]->object_layers_dialog });
         EVT_BUTTON($self, $self->{btn_settings}, sub { $_[0]->object_settings_dialog });
     }
     
@@ -1821,6 +1827,7 @@ sub async_apply_config {
     # reset preview canvases (invalidated contents will be hidden)
     $self->{toolpaths2D}->reload_print if $self->{toolpaths2D};
     $self->{preview3D}->reload_print if $self->{preview3D};
+    $self->{AdaptiveLayersDialog}->reload_preview if $self->{AdaptiveLayersDialog};
     
     if (!$Slic3r::GUI::Settings->{_}{background_processing}) {
         $self->hide_preview if $invalidated;
@@ -1897,6 +1904,7 @@ sub stop_background_process {
     
     $self->{toolpaths2D}->reload_print if $self->{toolpaths2D};
     $self->{preview3D}->reload_print if $self->{preview3D};
+    $self->{AdaptiveLayersDialog}->reload_preview if $self->{AdaptiveLayersDialog};
     
     if ($self->{process_thread}) {
         Slic3r::debugf "Killing background process.\n";
@@ -2040,6 +2048,7 @@ sub on_process_completed {
     return if $error;
     $self->{toolpaths2D}->reload_print if $self->{toolpaths2D};
     $self->{preview3D}->reload_print if $self->{preview3D};
+    $self->{AdaptiveLayersDialog}->reload_preview if $self->{AdaptiveLayersDialog};
     
     # if we have an export filename, start a new thread for exporting G-code
     if ($self->{export_gcode_output_file}) {
@@ -2619,10 +2628,17 @@ sub object_cut_dialog {
 	}
 }
 
-sub object_settings_dialog {
+sub object_layers_dialog {
     my $self = shift;
     my ($obj_idx) = @_;
 
+    $self->object_settings_dialog($obj_idx, adaptive_layers => 1);
+}
+
+sub object_settings_dialog {
+    my $self = shift;
+    my ($obj_idx, %params) = @_;
+    
     if (!defined $obj_idx) {
         ($obj_idx, undef) = $self->selected_object;
     }
@@ -2636,9 +2652,15 @@ sub object_settings_dialog {
     my $dlg = Slic3r::GUI::Plater::ObjectSettingsDialog->new($self,
 		object          => $self->{objects}[$obj_idx],
 		model_object    => $model_object,
+		obj_idx         => $obj_idx,
 	);
+	# store pointer to the adaptive layer tab to push preview updates
+	$self->{AdaptiveLayersDialog} = $dlg->{adaptive_layers};
+	# and jump directly to the tab if called by "promo-button"
+	$dlg->{tabpanel}->SetSelection(1) if $params{adaptive_layers};
 	$self->pause_background_process;
 	$dlg->ShowModal;
+	$self->{AdaptiveLayersDialog} = undef;
 	
     # update thumbnail since parts may have changed
     if ($dlg->PartsChanged) {
@@ -2695,11 +2717,11 @@ sub selection_changed {
     
     my $method = $have_sel ? 'Enable' : 'Disable';
     $self->{"btn_$_"}->$method
-        for grep $self->{"btn_$_"}, qw(remove increase decrease rotate45cw rotate45ccw changescale split cut settings);
+        for grep $self->{"btn_$_"}, qw(remove increase decrease rotate45cw rotate45ccw changescale split cut layers settings);
     
     if ($self->{htoolbar}) {
         $self->{htoolbar}->EnableTool($_, $have_sel)
-            for (TB_REMOVE, TB_MORE, TB_FEWER, TB_45CW, TB_45CCW, TB_SCALE, TB_SPLIT, TB_CUT, TB_SETTINGS);
+            for (TB_REMOVE, TB_MORE, TB_FEWER, TB_45CW, TB_45CCW, TB_SCALE, TB_SPLIT, TB_CUT, TB_LAYERS, TB_SETTINGS);
     }
     
     if ($self->{object_info_size}) { # have we already loaded the info pane?
@@ -2917,6 +2939,9 @@ sub object_menu {
     wxTheApp->append_menu_item($menu, "Cut…", 'Open the 3D cutting tool', sub {
         $self->object_cut_dialog;
     }, undef, 'package.png');
+    wxTheApp->append_menu_item($menu, "Layer heights…", 'Open the dynamic layer height control', sub {
+        $self->object_layers_dialog;
+    }, undef, 'cog.png');
     $menu->AppendSeparator();
     wxTheApp->append_menu_item($menu, "Settings…", 'Open the object editor dialog', sub {
         $self->object_settings_dialog;

lib/Slic3r/GUI/Plater/3DPreview.pm

@@ -75,15 +75,15 @@ sub new {
 }
 
 sub reload_print {
-    my ($self) = @_;
+    my ($self, $obj_idx) = @_;
     
     $self->canvas->reset_objects;
     $self->_loaded(0);
-    $self->load_print;
+    $self->load_print($obj_idx);
 }
 
 sub load_print {
-    my ($self) = @_;
+    my ($self, $obj_idx) = @_;
     
     return if $self->_loaded;
     
@@ -100,10 +100,16 @@ sub load_print {
     my $z_idx;
     {
         my %z = ();  # z => 1
-        foreach my $object (@{$self->{print}->objects}) {
+        if(defined $obj_idx) { # Load only given object
-            foreach my $layer (@{$object->layers}, @{$object->support_layers}) {
+            foreach my $layer (@{$self->{print}->get_object($obj_idx)->layers}) {
                 $z{$layer->print_z} = 1;
             }
+        }else{ # Load all objects on the plater + support material
+            foreach my $object (@{$self->{print}->objects}) {
+                foreach my $layer (@{$object->layers}, @{$object->support_layers}) {
+                    $z{$layer->print_z} = 1;
+                }
+            }
         }
         $self->enabled(1);
         $self->{layers_z} = [ sort { $a <=> $b } keys %z ];
@@ -129,14 +135,18 @@ sub load_print {
             $self->canvas->colors([ $self->canvas->default_colors ]);
         }
         
-        # load skirt and brim
+        if(defined $obj_idx) { # Load only one object
-        $self->canvas->load_print_toolpaths($self->print);
+            $self->canvas->load_print_object_toolpaths($self->{print}->get_object($obj_idx));
+        }else{ # load all objects
+	        # load skirt and brim
+            $self->canvas->load_print_toolpaths($self->print);
             
-        foreach my $object (@{$self->print->objects}) {
+            foreach my $object (@{$self->print->objects}) {
-            $self->canvas->load_print_object_toolpaths($object);
+                $self->canvas->load_print_object_toolpaths($object);
 
-            #my @volume_ids = $self->canvas->load_object($object->model_object);
+                #my @volume_ids = $self->canvas->load_object($object->model_object);
-            #$self->canvas->volumes->[$_]->color->[3] = 0.2 for @volume_ids;
+                #$self->canvas->volumes->[$_]->color->[3] = 0.2 for @volume_ids;
+            }
         }
         $self->_loaded(1);
     }

lib/Slic3r/GUI/Plater/ObjectSettingsDialog.pm

@@ -14,11 +14,16 @@ use base 'Wx::Dialog';
 sub new {
     my $class = shift;
     my ($parent, %params) = @_;
-    my $self = $class->SUPER::new($parent, -1, "Settings for " . $params{object}->name, wxDefaultPosition, [700,500], wxDEFAULT_DIALOG_STYLE | wxRESIZE_BORDER);
+    my $self = $class->SUPER::new($parent, -1, "Settings for " . $params{object}->name, wxDefaultPosition, [1000,700], wxDEFAULT_DIALOG_STYLE | wxRESIZE_BORDER);
     $self->{$_} = $params{$_} for keys %params;
     
     $self->{tabpanel} = Wx::Notebook->new($self, -1, wxDefaultPosition, wxDefaultSize, wxNB_TOP | wxTAB_TRAVERSAL);
     $self->{tabpanel}->AddPage($self->{parts} = Slic3r::GUI::Plater::ObjectPartsPanel->new($self->{tabpanel}, model_object => $params{model_object}), "Parts");
+    $self->{tabpanel}->AddPage($self->{adaptive_layers} = Slic3r::GUI::Plater::ObjectDialog::AdaptiveLayersTab->new( $self->{tabpanel}, 
+        plater         => $parent,
+        model_object   => $params{model_object},
+        obj_idx        => $params{obj_idx}
+        ), "Adaptive Layers");
     $self->{tabpanel}->AddPage($self->{layers} = Slic3r::GUI::Plater::ObjectDialog::LayersTab->new($self->{tabpanel}), "Layers");
     
     my $buttons = $self->CreateStdDialogButtonSizer(wxOK);
@@ -68,6 +73,168 @@ sub model_object {
     return $self->GetParent->GetParent->{model_object};
 }
 
+package Slic3r::GUI::Plater::ObjectDialog::AdaptiveLayersTab;
+use Slic3r::Geometry qw(X Y Z scale unscale);
+use Slic3r::Print::State ':steps';
+use List::Util qw(min max sum first);
+use Wx qw(wxTheApp :dialog :id :misc :sizer :systemsettings :statictext wxTAB_TRAVERSAL);
+use base 'Slic3r::GUI::Plater::ObjectDialog::BaseTab';
+
+sub new {
+    my $class = shift;
+    my ($parent, %params) = @_;
+    my $self = $class->SUPER::new($parent, -1, wxDefaultPosition, wxDefaultSize);
+    my $model_object = $self->{model_object} = $params{model_object};
+    my $obj_idx = $self->{obj_idx} = $params{obj_idx};
+    my $plater = $self->{plater} = $params{plater};
+    my $object = $self->{object} = $self->{plater}->{print}->get_object($self->{obj_idx});
+
+    # store last raft height to correctly draw z-indicator plane during a running background job where the printObject is not valid
+    $self->{last_raft_height} = 0;
+
+    # Initialize 3D toolpaths preview
+    if ($Slic3r::GUI::have_OpenGL) {
+        $self->{preview3D} = Slic3r::GUI::Plater::3DPreview->new($self, $plater->{print});
+        $self->{preview3D}->canvas->set_auto_bed_shape;
+        $self->{preview3D}->canvas->SetSize([500,500]);
+        $self->{preview3D}->canvas->SetMinSize($self->{preview3D}->canvas->GetSize);
+        # object already processed?
+        wxTheApp->CallAfter(sub {
+            if (!$plater->{processed}) {
+                $self->_trigger_slicing(0); # trigger processing without invalidating STEP_SLICE to keep current height distribution
+            }else{
+                $self->{preview3D}->reload_print($obj_idx);
+                $self->{preview3D}->canvas->zoom_to_volumes;
+                $self->{preview_zoomed} = 1;
+            }
+        });
+    }
+
+    $self->{splineControl} = Slic3r::GUI::Plater::SplineControl->new($self, Wx::Size->new(150, 200), $object);
+
+    my $optgroup;
+    $optgroup = $self->{optgroup} = Slic3r::GUI::OptionsGroup->new(
+        parent      => $self,
+        title       => 'Adaptive quality %',
+        on_change   => sub {
+            my ($opt_id) = @_;
+            # There seems to be an issue with wxWidgets 3.0.2/3.0.3, where the slider
+            # genates tens of events for a single value change.
+            # Only trigger the recalculation if the value changes
+            # or a live preview was activated and the mesh cut is not valid yet.
+            if ($self->{adaptive_quality} != $optgroup->get_value($opt_id)) {
+                $self->{adaptive_quality} = $optgroup->get_value($opt_id);
+                $self->{model_object}->config->set('adaptive_slicing_quality', $optgroup->get_value($opt_id));
+                $self->{object}->config->set('adaptive_slicing_quality', $optgroup->get_value($opt_id));
+                $self->{object}->invalidate_step(STEP_LAYERS);
+                # trigger re-slicing
+                $self->_trigger_slicing;
+            }
+        },
+        label_width  => 0,
+    );
+
+    $optgroup->append_single_option_line(Slic3r::GUI::OptionsGroup::Option->new(
+        opt_id      => 'adaptive_slicing_quality',
+        type        => 'slider',
+        label       => '',
+        default     => $object->config->get('adaptive_slicing_quality'),
+        min         => 0,
+        max         => 100,
+        full_width  => 1,
+    ));
+    $optgroup->get_field('adaptive_slicing_quality')->set_scale(1);
+    $self->{adaptive_quality} = $object->config->get('adaptive_slicing_quality');
+    # init quality slider
+    if(!$object->config->get('adaptive_slicing')) {
+        # disable slider
+        $optgroup->get_field('adaptive_slicing_quality')->disable;
+    }
+
+    my $right_sizer = Wx::BoxSizer->new(wxVERTICAL);
+    $right_sizer->Add($self->{splineControl}, 1, wxEXPAND | wxALL, 0);
+    $right_sizer->Add($optgroup->sizer, 0, wxEXPAND | wxALL, 0);
+
+
+    $self->{sizer} = Wx::BoxSizer->new(wxHORIZONTAL);
+    $self->{sizer}->Add($self->{preview3D}, 3, wxEXPAND | wxTOP | wxBOTTOM, 0) if $self->{preview3D};
+    $self->{sizer}->Add($right_sizer, 1, wxEXPAND | wxTOP | wxBOTTOM, 10);
+
+    $self->SetSizerAndFit($self->{sizer});
+
+    # init spline control values
+    # determine min and max layer height from perimeter extruder capabilities.
+    my %extruders;
+    for my $region_id (0 .. ($object->region_count - 1)) {
+        foreach (qw(perimeter_extruder infill_extruder solid_infill_extruder)) {
+            my $extruder_id = $self->{plater}->{print}->get_region($region_id)->config->get($_)-1;
+            $extruders{$extruder_id} = $extruder_id;
+        }
+    }
+    my $min_height = max(map {$self->{plater}->{print}->config->get_at('min_layer_height', $_)} (values %extruders));
+    my $max_height = min(map {$self->{plater}->{print}->config->get_at('max_layer_height', $_)} (values %extruders));
+
+    $self->{splineControl}->set_size_parameters($min_height, $max_height, unscale($object->size->z));
+
+    $self->{splineControl}->on_layer_update(sub {
+        # trigger re-slicing
+        $self->_trigger_slicing;
+    });
+
+    $self->{splineControl}->on_z_indicator(sub {
+        my ($z) = @_;
+
+        if($z) { # compensate raft height
+	        $z += $self->{last_raft_height};
+        }
+        $self->{preview3D}->canvas->SetCuttingPlane(Z, $z, []);
+        $self->{preview3D}->canvas->Render;
+    });
+
+    return $self;
+}
+
+# This is called by the plater after processing to update the preview and spline
+sub reload_preview {
+    my ($self) = @_;
+    $self->{splineControl}->update;
+    $self->{preview3D}->reload_print($self->{obj_idx});
+    my $object = $self->{plater}->{print}->get_object($self->{obj_idx});
+    if($object->layer_count-1 > 0) {
+        my $first_layer = $self->{object}->get_layer(0);
+        $self->{last_raft_height} = max(0, $first_layer->print_z - $first_layer->height);
+        $self->{preview3D}->set_z(unscale($self->{object}->size->z));
+        if(!$self->{preview_zoomed}) {
+            $self->{preview3D}->canvas->set_auto_bed_shape;
+            $self->{preview3D}->canvas->zoom_to_volumes;
+            $self->{preview_zoomed} = 1;
+        }
+    }
+}
+
+# Trigger background slicing at the plater
+sub _trigger_slicing {
+    my ($self, $invalidate) = @_;
+    $invalidate //= 1;
+    my $object = $self->{plater}->{print}->get_object($self->{obj_idx});
+    $self->{model_object}->set_layer_height_spline($self->{object}->layer_height_spline); # push modified spline object to model_object
+    #$self->{plater}->pause_background_process;
+    $self->{plater}->stop_background_process;
+    if (!$Slic3r::GUI::Settings->{_}{background_processing}) {
+        $self->{plater}->statusbar->SetCancelCallback(sub {
+            $self->{plater}->stop_background_process;
+            $self->{plater}->statusbar->SetStatusText("Slicing cancelled");
+            $self->{plater}->preview_notebook->SetSelection(0);
+        });
+        $object->invalidate_step(STEP_SLICE) if($invalidate);
+        $self->{plater}->start_background_process;
+    }else{
+        $object->invalidate_step(STEP_SLICE) if($invalidate);
+        $self->{plater}->schedule_background_process;
+    }
+}
+
+
 package Slic3r::GUI::Plater::ObjectDialog::LayersTab;
 use Wx qw(:dialog :id :misc :sizer :systemsettings);
 use Wx::Grid;

lib/Slic3r/GUI/Plater/SplineControl.pm

@@ -0,0 +1,378 @@
+package Slic3r::GUI::Plater::SplineControl;
+use strict;
+use warnings;
+use utf8;
+
+use List::Util qw(min max first);
+use Slic3r::Geometry qw(X Y scale unscale);
+use Wx qw(:misc :pen :brush :sizer :font :cursor wxTAB_TRAVERSAL);
+use Wx::Event qw(EVT_MOUSE_EVENTS EVT_PAINT EVT_ERASE_BACKGROUND EVT_SIZE);
+use base 'Wx::Panel';
+
+sub new {
+    my $class = shift;
+    my ($parent, $size, $object) = @_;
+    
+    my $self = $class->SUPER::new($parent, -1, wxDefaultPosition, $size, wxTAB_TRAVERSAL);
+    
+    $self->{object} = $object;
+    $self->{is_valid} = 0;
+    
+    # This has only effect on MacOS. On Windows and Linux/GTK, the background is painted by $self->repaint().
+    $self->SetBackgroundColour(Wx::wxWHITE);
+
+    $self->{line_pen}           = Wx::Pen->new(Wx::Colour->new(50,50,50), 1, wxSOLID);
+    $self->{original_pen}       = Wx::Pen->new(Wx::Colour->new(200,200,200), 1, wxSOLID);
+    $self->{interactive_pen}    = Wx::Pen->new(Wx::Colour->new(255,0,0), 1, wxSOLID);
+    $self->{resulting_pen}      = Wx::Pen->new(Wx::Colour->new(5,120,160), 1, wxSOLID);
+
+    $self->{user_drawn_background} = $^O ne 'darwin';
+
+    # scale plot data to actual canvas, documentation in set_size_parameters
+    $self->{scaling_factor_x} = 1;
+    $self->{scaling_factor_y} = 1;
+    $self->{min_layer_height} = 0.1;
+    $self->{max_layer_height} = 0.4;
+    $self->{mousover_layer_height} = undef; # display layer height at mousover position
+    $self->{object_height} = 1.0;
+    
+    # initialize values
+    $self->update;
+    
+    EVT_PAINT($self, \&repaint);
+    EVT_ERASE_BACKGROUND($self, sub {}) if $self->{user_drawn_background};
+    EVT_MOUSE_EVENTS($self, \&mouse_event);
+    EVT_SIZE($self, sub {
+        $self->_update_canvas_size;
+        $self->Refresh;
+    });
+
+    return $self;
+}
+
+sub repaint {
+    my ($self, $event) = @_;
+    
+    my $dc = Wx::AutoBufferedPaintDC->new($self);
+    my $size = $self->GetSize;
+    my @size = ($size->GetWidth, $size->GetHeight);
+    
+
+    if ($self->{user_drawn_background}) {
+        # On all systems the AutoBufferedPaintDC() achieves double buffering.
+        # On MacOS the background is erased, on Windows the background is not erased 
+        # and on Linux/GTK the background is erased to gray color.
+        # Fill DC with the background on Windows & Linux/GTK.
+        my $brush_background = Wx::Brush->new(Wx::wxWHITE, wxSOLID);
+        $dc->SetPen(wxWHITE_PEN);
+        $dc->SetBrush($brush_background);
+        my $rect = $self->GetUpdateRegion()->GetBox();
+        $dc->DrawRectangle($rect->GetLeft(), $rect->GetTop(), $rect->GetWidth(), $rect->GetHeight());
+    }
+    
+    # draw scale (min and max indicator at the bottom)
+    $dc->SetTextForeground(Wx::Colour->new(0,0,0));
+    $dc->SetFont(Wx::Font->new(10, wxDEFAULT, wxNORMAL, wxNORMAL));
+    $dc->DrawLabel(sprintf('%.4g', $self->{min_layer_height}), Wx::Rect->new(0, $size[1]/2, $size[0], $size[1]/2), wxALIGN_LEFT | wxALIGN_BOTTOM);
+    $dc->DrawLabel(sprintf('%.2g', $self->{max_layer_height}), Wx::Rect->new(0, $size[1]/2, $size[0], $size[1]/2), wxALIGN_RIGHT | wxALIGN_BOTTOM);
+    if($self->{mousover_layer_height}){
+        $dc->DrawLabel(sprintf('%4.2fmm', $self->{mousover_layer_height}), Wx::Rect->new(0, 0, $size[0], $size[1]), wxALIGN_RIGHT | wxALIGN_TOP);
+    }
+
+    if($self->{is_valid}) {
+
+        # draw original spline as reference
+        if($self->{original_height_spline}) {
+            #draw spline
+            $self->_draw_layer_height_spline($dc, $self->{original_height_spline}, $self->{original_pen});
+        }
+
+        # draw interactive (currently modified by the user) layers as lines and spline
+        if($self->{interactive_height_spline}) {
+            # draw layer lines
+            my @interpolated_layers = @{$self->{interactive_height_spline}->getInterpolatedLayers};
+            $self->_draw_layers_as_lines($dc, $self->{interactive_pen}, \@interpolated_layers);
+
+            #draw spline
+            $self->_draw_layer_height_spline($dc, $self->{interactive_height_spline}, $self->{interactive_pen});
+        }
+
+        # draw resulting layers as lines
+        unless($self->{interactive_heights}) {
+            $self->_draw_layers_as_lines($dc, $self->{resulting_pen}, $self->{interpolated_layers});
+        }
+
+        # Always draw current BSpline, gives a reference during a modification
+        $self->_draw_layer_height_spline($dc, $self->{object}->layer_height_spline, $self->{line_pen});
+    }
+
+    $event->Skip;
+}
+
+# Set basic parameters for this control.
+# min/max_layer_height are required to define the x-range, object_height is used to scale the y-range.
+# Must be called if object selection changes.
+sub set_size_parameters {
+    my ($self, $min_layer_height, $max_layer_height, $object_height) = @_;
+    
+    $self->{min_layer_height} = $min_layer_height;
+    $self->{max_layer_height} = $max_layer_height;
+    $self->{object_height} = $object_height;
+    
+    $self->_update_canvas_size;
+    $self->Refresh;
+}
+
+# Layers have been modified externally, re-initialize this control with new values
+sub update {
+    my $self = shift;
+
+    if(($self->{object}->layer_height_spline->layersUpdated || !$self->{heights}) && $self->{object}->layer_height_spline->hasData) {
+        $self->{original_height_spline} = $self->{object}->layer_height_spline->clone; # make a copy to display the unmodified original spline
+        $self->{original_layers} = $self->{object}->layer_height_spline->getOriginalLayers;
+        $self->{interpolated_layers} = $self->{object}->layer_height_spline->getInterpolatedLayers; # Initialize to current values
+    
+        # initialize height vector
+        $self->{heights} = ();
+        $self->{interactive_heights} = ();
+        foreach my $z (@{$self->{original_layers}}) {
+            push (@{$self->{heights}}, $self->{object}->layer_height_spline->getLayerHeightAt($z));
+        }
+        $self->{is_valid} = 1;
+    }
+    $self->Refresh;
+}
+
+# Callback to notify parent element if layers have changed and reslicing should be triggered
+sub on_layer_update {
+    my ($self, $cb) = @_;
+    $self->{on_layer_update} = $cb;
+}
+
+# Callback to tell parent element at which z-position the mouse currently hovers to update indicator in 3D-view
+sub on_z_indicator {
+    my ($self, $cb) = @_;
+    $self->{on_z_indicator} = $cb;
+}
+
+
+sub mouse_event {
+    my ($self, $event) = @_;
+    
+    my $pos = $event->GetPosition;
+    my @obj_pos = $self->pixel_to_point($pos);
+    
+    if ($event->ButtonDown) {
+        if ($event->LeftDown) {
+           # start dragging
+           $self->{left_drag_start_pos} = $pos;
+           $self->{interactive_height_spline} = $self->{object}->layer_height_spline->clone;
+        }
+        if ($event->RightDown) {
+           # start dragging
+           $self->{right_drag_start_pos} = $pos;
+           $self->{interactive_height_spline} = $self->{object}->layer_height_spline->clone;
+        }
+    } elsif ($event->LeftUp) {
+        if($self->{left_drag_start_pos}) {
+            $self->_modification_done;
+        }
+        $self->{left_drag_start_pos} = undef;
+    } elsif ($event->RightUp) {
+        if($self->{right_drag_start_pos}) {
+            $self->_modification_done;
+        }
+        $self->{right_drag_start_pos} = undef;
+    } elsif ($event->Dragging) {
+        if($self->{left_drag_start_pos}) {
+        
+            my @start_pos = $self->pixel_to_point($self->{left_drag_start_pos});
+            my $range = abs($start_pos[1] - $obj_pos[1]);
+
+            # compute updated interactive layer heights
+            $self->_interactive_quadratic_curve($start_pos[1], $obj_pos[0], $range);
+
+            unless($self->{interactive_height_spline}->updateLayerHeights($self->{interactive_heights})) {
+                die "Unable to update interactive interpolated layers!\n";
+            }
+            $self->Refresh;
+        } elsif($self->{right_drag_start_pos}) {
+            my @start_pos = $self->pixel_to_point($self->{right_drag_start_pos});
+            my $range = $obj_pos[1] - $start_pos[1];
+
+            # compute updated interactive layer heights
+            $self->_interactive_linear_curve($start_pos[1], $obj_pos[0], $range);
+            unless($self->{interactive_height_spline}->updateLayerHeights($self->{interactive_heights})) {
+                die "Unable to update interactive interpolated layers!\n";
+            }
+            $self->Refresh;
+        }
+    } elsif ($event->Moving) {
+        if($self->{on_z_indicator}) {
+            $self->{on_z_indicator}->($obj_pos[1]);
+            $self->{mousover_layer_height} = $self->{object}->layer_height_spline->getLayerHeightAt($obj_pos[1]);
+            $self->Refresh;
+            $self->Update;
+        }
+    } elsif ($event->Leaving) {
+        if($self->{on_z_indicator} && !$self->{left_drag_start_pos}) {
+            $self->{on_z_indicator}->(undef);
+        }
+        $self->{mousover_layer_height} = undef;
+        $self->Refresh;
+        $self->Update;
+    }
+}
+
+# Push modified heights to the spline object and update after user modification
+sub _modification_done {
+    my $self = shift;
+    
+    if($self->{interactive_heights}) {
+        $self->{heights} = $self->{interactive_heights};
+        $self->{interactive_heights} = ();
+        # update spline database
+        unless($self->{object}->layer_height_spline->updateLayerHeights($self->{heights})) {
+            die "Unable to update interpolated layers!\n";
+        }
+        $self->{interpolated_layers} = $self->{object}->layer_height_spline->getInterpolatedLayers;
+    }
+    $self->Refresh;
+    $self->{on_layer_update}->(@{$self->{interpolated_layers}});
+    $self->{interactive_height_spline} = undef;
+}
+
+# Internal function to cache scaling factors
+sub _update_canvas_size {
+    my $self = shift;
+    
+    # when the canvas is not rendered yet, its GetSize() method returns 0,0
+    my $canvas_size = $self->GetSize;
+    my ($canvas_w, $canvas_h) = ($canvas_size->GetWidth, $canvas_size->GetHeight);
+    return if $canvas_w == 0;
+
+    my $padding = $self->{canvas_padding} = 10; # border size in pixels
+
+    my @size = ($canvas_w - 2*$padding, $canvas_h - 2*$padding);
+    $self->{canvas_size} = [@size];
+    
+    $self->{scaling_factor_x} = $size[0]/($self->{max_layer_height} - $self->{min_layer_height});
+    $self->{scaling_factor_y} = $size[1]/$self->{object_height};
+}
+
+# calculate new set of layers with quadaratic modifier for interactive display
+sub _interactive_quadratic_curve {
+    my ($self, $mod_z, $target_layer_height, $range) = @_;
+    
+    $self->{interactive_heights} = (); # reset interactive curve
+    
+    # iterate over original points provided by spline
+    my $last_z = 0;
+    foreach my $i (0..@{$self->{heights}}-1  ) {
+        my $z = $self->{original_layers}[$i];
+        my $layer_h = $self->{heights}[$i];
+        my $quadratic_factor = $self->_quadratic_factor($mod_z, $range, $z);
+        my $diff = $target_layer_height - $layer_h;
+        $layer_h  += $diff * $quadratic_factor;
+        push (@{$self->{interactive_heights}}, $layer_h);
+    }
+}
+
+# calculate new set of layers with linear modifier for interactive display
+sub _interactive_linear_curve {
+    my ($self, $mod_z, $target_layer_height, $range) = @_;
+
+    $self->{interactive_heights} = (); # reset interactive curve
+    my $from;
+    my $to;
+
+    if($range >= 0) {
+        $from = $mod_z;
+        $to = $mod_z + $range;
+    }else{
+        $from = $mod_z + $range;
+        $to = $mod_z;
+    }
+
+    # iterate over original points provided by spline
+    foreach my $i (0..@{$self->{heights}}-1  ) {
+        if(($self->{original_layers}[$i]) >= $from && ($self->{original_layers}[$i]< $to)) {
+            push (@{$self->{interactive_heights}}, $target_layer_height);
+        }else{
+           push (@{$self->{interactive_heights}}, $self->{heights}[$i]);
+        }
+    }
+}
+
+sub _quadratic_factor {
+    my ($self, $fixpoint, $range, $value) = @_;
+    
+    # avoid division by zero
+    $range = 0.00001 if $range <= 0;
+    
+    my $dist = abs($fixpoint - $value);
+    my $x = $dist/$range; # normalize
+    my $result = 1-($x*$x);
+    
+    return max(0, $result);
+}
+
+# Draw a set of layers as lines
+sub _draw_layers_as_lines {
+    my ($self, $dc, $pen, $layers) = @_;
+    
+    $dc->SetPen($pen);
+    my $last_z = 0.0;    
+    foreach my $z (@$layers) {
+        my $layer_h = $z - $last_z;
+        my $pl = $self->point_to_pixel(0, $z);
+        my $pr = $self->point_to_pixel($layer_h, $z);
+        $dc->DrawLine($pl->x, $pl->y, $pr->x, $pr->y);
+        $last_z = $z;
+    }
+}
+
+# Draw the resulting spline from a LayerHeightSpline object over the full canvas height
+sub _draw_layer_height_spline {
+    my ($self, $dc, $layer_height_spline, $pen) = @_;
+    
+    my @size = @{$self->{canvas_size}};
+    
+    $dc->SetPen($pen);
+    my @points = ();
+    foreach my $pixel (0..$size[1]) {
+        my @z = $self->pixel_to_point(Wx::Point->new(0, $pixel));
+        my $h = $layer_height_spline->getLayerHeightAt($z[1]);
+        my $p = $self->point_to_pixel($h, $z[1]);
+        push (@points, $p);
+    }
+    $dc->DrawLines(\@points);
+}
+
+# Takes a 2-tupel [layer_height (x), height(y)] and converts it
+# into a Wx::Point in scaled canvas coordinates
+sub point_to_pixel {
+    my ($self, @point) = @_;
+    
+    my @size = @{$self->{canvas_size}};
+    
+    my $x = ($point[0] - $self->{min_layer_height})*$self->{scaling_factor_x} + $self->{canvas_padding};
+    my $y = $size[1] - $point[1]*$self->{scaling_factor_y} + $self->{canvas_padding}; # invert y-axis
+
+    return Wx::Point->new(min(max($x, $self->{canvas_padding}), $size[0]+$self->{canvas_padding}), min(max($y, $self->{canvas_padding}), $size[1]+$self->{canvas_padding})); # limit to canvas size
+}
+
+# Takes a Wx::Point in scaled canvas coordinates and converts it
+# into a 2-tupel [layer_height (x), height(y)]
+sub pixel_to_point {
+    my ($self, $point) = @_;
+    
+    my @size = @{$self->{canvas_size}};
+    
+    my $x = ($point->x-$self->{canvas_padding})/$self->{scaling_factor_x} + $self->{min_layer_height};
+    my $y = ($size[1] - $point->y)/$self->{scaling_factor_y}; # invert y-axis
+    
+    return (min(max($x, $self->{min_layer_height}), $self->{max_layer_height}), min(max($y, 0), $self->{object_height})); # limit to object size and layer constraints 
+}
+
+1;

lib/Slic3r/GUI/PresetEditor.pm

@@ -437,6 +437,7 @@ sub title { 'Print Settings' }
 sub options {
     return qw(
         layer_height first_layer_height
+        adaptive_slicing adaptive_slicing_quality match_horizontal_surfaces
         perimeters spiral_vase
         top_solid_layers bottom_solid_layers
         extra_perimeters avoid_crossing_perimeters thin_walls overhangs
@@ -512,6 +513,10 @@ sub build {
             my $optgroup = $page->new_optgroup('Layer height');
             $optgroup->append_single_option_line('layer_height');
             $optgroup->append_single_option_line('first_layer_height');
+            $optgroup->append_single_option_line('adaptive_slicing');
+            $optgroup->append_single_option_line('adaptive_slicing_quality');
+            $optgroup->get_field('adaptive_slicing_quality')->set_scale(1);
+            $optgroup->append_single_option_line('match_horizontal_surfaces');
         }
         {
             my $optgroup = $page->new_optgroup('Vertical shells');
@@ -865,6 +870,12 @@ sub _update {
             external_perimeter_extrusion_width
             perimeter_speed small_perimeter_speed external_perimeter_speed);
 
+    my $have_adaptive_slicing = $config->adaptive_slicing;
+    $self->get_field($_)->toggle($have_adaptive_slicing)
+        for qw(adaptive_slicing_quality match_horizontal_surfaces);
+    $self->get_field($_)->toggle(!$have_adaptive_slicing)
+        for qw(layer_height);
+    
     my $have_infill = $config->fill_density > 0;
     # infill_extruder uses the same logic as in Print::extruders()
     $self->get_field($_)->toggle($have_infill)
@@ -1146,10 +1157,10 @@ sub _update {
     
     $self->_update_description;
     
-    my $cooling = $self->config->cooling;
+    my $cooling = $self->{config}->cooling;
     $self->get_field($_)->toggle($cooling)
         for qw(max_fan_speed fan_below_layer_time slowdown_below_layer_time min_print_speed);
-    $self->get_field($_)->toggle($cooling || $self->config->fan_always_on)
+    $self->get_field($_)->toggle($cooling || $self->{config}->fan_always_on)
         for qw(min_fan_speed disable_fan_first_layers);
 }
 
@@ -1199,7 +1210,7 @@ sub options {
         use_firmware_retraction pressure_advance vibration_limit
         use_volumetric_e
         start_gcode end_gcode before_layer_gcode layer_gcode toolchange_gcode between_objects_gcode
-        nozzle_diameter extruder_offset
+        nozzle_diameter extruder_offset min_layer_height max_layer_height
         retract_length retract_lift retract_speed retract_restart_extra retract_before_travel retract_layer_change wipe
         retract_length_toolchange retract_restart_extra_toolchange retract_lift_above retract_lift_below
         printer_settings_id
@@ -1444,7 +1455,7 @@ sub _extruders_count_changed {
     $self->_update;
 }
 
-sub _extruder_options { qw(nozzle_diameter extruder_offset retract_length retract_lift retract_lift_above retract_lift_below retract_speed retract_restart_extra retract_before_travel wipe
+sub _extruder_options { qw(nozzle_diameter min_layer_height max_layer_height extruder_offset retract_length retract_lift retract_lift_above retract_lift_below retract_speed retract_restart_extra retract_before_travel wipe
     retract_layer_change retract_length_toolchange retract_restart_extra_toolchange) }
 
 sub _build_extruder_pages {
@@ -1473,6 +1484,11 @@ sub _build_extruder_pages {
             my $optgroup = $page->new_optgroup('Size');
             $optgroup->append_single_option_line('nozzle_diameter', $extruder_idx);
         }
+        {
+            my $optgroup = $page->new_optgroup('Limits');
+            $optgroup->append_single_option_line($_, $extruder_idx)
+               for qw(min_layer_height max_layer_height);
+        }
         {
             my $optgroup = $page->new_optgroup('Position (for multi-extruder printers)');
             $optgroup->append_single_option_line('extruder_offset', $extruder_idx);

lib/Slic3r/Print/State.pm

@@ -5,7 +5,7 @@ use warnings;
 
 require Exporter;
 our @ISA = qw(Exporter);
-our @EXPORT_OK   = qw(STEP_SLICE STEP_PERIMETERS STEP_PREPARE_INFILL 
+our @EXPORT_OK   = qw(STEP_LAYERS STEP_SLICE STEP_PERIMETERS STEP_PREPARE_INFILL
                     STEP_INFILL STEP_SUPPORTMATERIAL STEP_SKIRT STEP_BRIM);
 our %EXPORT_TAGS = (steps => \@EXPORT_OK);
 

lib/Slic3r/Test.pm

@@ -147,6 +147,16 @@ sub mesh {
         $facets = [
             [0,1,2],[1,3,2],[3,4,5],[2,3,5],[6,0,2],[6,2,7],[5,8,7],[5,7,2],[9,10,8],[9,8,5],[9,0,6],[9,6,10],[9,11,1],[9,1,0],[3,1,11],[4,3,11],[5,11,9],[5,4,11],[12,10,6],[12,6,13],[6,7,14],[13,6,14],[7,8,15],[14,7,15],[15,8,10],[15,10,12],[12,13,14],[12,14,15]
         ];
+    } elsif ($name eq 'slopy_cube') {
+    	$vertices = [
+	        [-10,-10,0], [-10,-10,20], [-10,10,0], [-10,10,20], [0,-10,10], [10,-10,0], [2.92893,-10,10], [10,-10,2.92893], 
+	        [0,-10,20], [10,10,0], [0,10,10], [0,10,20], [2.92893,10,10], [10,10,2.92893],
+	    ];
+	    $facets = [
+	    	[0,1,2], [2,1,3], [4,0,5], [4,1,0], [6,4,7], [7,4,5], [4,8,1], [0,2,5], [5,2,9], [2,10,9], [10,3,11], 
+	    	[2,3,10], [9,10,12], [13,9,12], [3,1,8], [11,3,8], [10,11,8], [4,10,8], [6,12,10], [4,6,10], 
+	    	[7,13,12], [6,7,12], [7,5,9], [13,7,9],
+	    ];
     } else {
         return undef;
     }

t/adaptive_slicing.t

@@ -0,0 +1,147 @@
+use Test::More tests => 9;
+use strict;
+use warnings;
+
+BEGIN {
+    use FindBin;
+    use lib "$FindBin::Bin/../lib";
+    use local::lib "$FindBin::Bin/../local-lib";
+}
+
+use List::Util qw(first sum);
+use Slic3r;
+use Slic3r::Test qw(_eq);
+use Slic3r::Geometry qw(Z PI scale unscale);
+
+use Devel::Peek;
+
+my $config = Slic3r::Config->new_from_defaults;
+
+my $generate_gcode = sub {
+    my ($conf) = @_;
+    $conf ||= $config;
+    
+    my $print = Slic3r::Test::init_print('slopy_cube', config => $conf);
+    
+    my @z = ();
+    my @increments = ();
+    Slic3r::GCode::Reader->new->parse(Slic3r::Test::gcode($print), sub {
+        my ($self, $cmd, $args, $info) = @_;
+        
+        if ($info->{dist_Z}) {
+            push @z, 1*$args->{Z};
+            push @increments, $info->{dist_Z};
+        }
+    });
+    
+    return (@z);
+};
+
+my $horizontal_feature_test = sub {
+    my (@z) = $generate_gcode->();
+    
+    ok  (_eq($z[0], $config->get_value('first_layer_height') + $config->z_offset), 'first layer height.');
+    
+    ok  (_eq($z[1], $config->get_value('first_layer_height') + $config->get('max_layer_height')->[0] + $config->z_offset), 'second layer height.');
+        
+    cmp_ok((first { _eq($_, 10.0) } @z[1..$#z]), '>', 0, 'horizontal facet matched');
+    
+    1;
+};
+
+my $height_gradation_test = sub {
+    my (@z) = $generate_gcode->();
+
+    my $gradation = 1 / $config->get('z_steps_per_mm');
+    # +1 is a "dirty fix" to avoid rounding issues with the modulo operator...
+    my @results = map {unscale((scale($_)+1) % scale($gradation))} @z;
+
+    ok  (_eq(sum(@results), 0), 'layer z is multiple of gradation ' . $gradation );
+    
+    1;
+};
+
+
+my $adaptive_slicing = Slic3r::SlicingAdaptive->new();
+my $mesh = Slic3r::Test::mesh('slopy_cube');
+$adaptive_slicing->add_mesh($mesh);
+$adaptive_slicing->prepare(20);
+
+
+subtest 'max layer_height limited by extruder capabilities' => sub {
+    plan tests => 3;
+
+    ok  (_eq($adaptive_slicing->next_layer_height(1, 20, 0.1, 0.15), 0.15), 'low');
+    ok  (_eq($adaptive_slicing->next_layer_height(1, 20, 0.1, 0.4), 0.4), 'higher');
+    ok  (_eq($adaptive_slicing->next_layer_height(1, 20, 0.1, 0.65), 0.65), 'highest');
+};
+  
+subtest 'min layer_height limited by extruder capabilities' => sub {
+    plan tests => 3;
+
+    ok  (_eq($adaptive_slicing->next_layer_height(4, 99, 0.1, 0.15), 0.1), 'low');
+    ok  (_eq($adaptive_slicing->next_layer_height(4, 98, 0.2, 0.4), 0.2), 'higher');
+    ok  (_eq($adaptive_slicing->next_layer_height(4, 99, 0.3, 0.65), 0.3), 'highest');
+};
+
+subtest 'correct layer_height depending on the facet normals' => sub {
+    plan tests => 3;
+
+    ok  (_eq($adaptive_slicing->next_layer_height(1, 10, 0.1, 0.5), 0.5), 'limit');
+    ok  (_eq($adaptive_slicing->next_layer_height(4, 80, 0.1, 0.5), 0.1546), '45deg facet, quality_value: 0.2');
+    ok  (_eq($adaptive_slicing->next_layer_height(4, 50, 0.1, 0.5), 0.3352), '45deg facet, quality_value: 0.5');
+};
+
+
+# 2.92893 ist lower slope edge
+# distance to slope must be higher than min extruder cap.
+# slopes layer height must be greater than the distance to the slope
+ok  (_eq($adaptive_slicing->next_layer_height(2.798, 80, 0.1, 0.5), 0.1546), 'reducing layer_height due to higher slopy facet');
+
+# slopes layer height must be smaller than the distance to the slope
+ok  (_eq($adaptive_slicing->next_layer_height(2.6289, 85, 0.1, 0.5), 0.3), 'reducing layer_height to z-diff');
+
+subtest 'horizontal planes' => sub {
+    plan tests => 3;
+
+    ok  (_eq($adaptive_slicing->horizontal_facet_distance(1, 1.2), 1.2), 'max_height limit');
+    ok  (_eq($adaptive_slicing->horizontal_facet_distance(8.5, 4), 1.5), 'normal horizontal facet');
+    ok  (_eq($adaptive_slicing->horizontal_facet_distance(17, 5), 3.0), 'object maximum');
+};
+
+# shrink current layer to fit another layer under horizontal facet
+$config->set('start_gcode', '');  # to avoid dealing with the nozzle lift in start G-code
+$config->set('z_offset', 0);
+
+$config->set('adaptive_slicing', 1);
+$config->set('first_layer_height', 0.42893); # to catch lower slope edge
+$config->set('nozzle_diameter', [0.5]);
+$config->set('min_layer_height', [0.1]);
+$config->set('max_layer_height', [0.5]);
+$config->set('adaptive_slicing_quality', [81]);
+# slope height: 7,07107 (2.92893 to 10)
+
+subtest 'shrink to match horizontal facets' => sub {
+    plan skip_all => 'spline smoothing currently prevents exact horizontal facet matching';
+    plan tests => 3;
+    $horizontal_feature_test->();
+};
+    
+# widen current layer to match horizontal facet
+$config->set('adaptive_slicing_quality', [0.1]);
+
+subtest 'widen to match horizontal facets' => sub {
+    plan skip_all => 'spline smoothing currently prevents exact horizontal facet matching';
+    plan tests => 3;
+    $horizontal_feature_test->();
+};
+
+subtest 'layer height gradation' => sub {
+    plan tests => 5;
+    foreach my $gradation (1/0.001*4, 1/0.01*4, 1/0.02*4, 1/0.05*4, 1/0.08*4) {
+        $config->set('z_steps_per_mm', $gradation);
+        $height_gradation_test->();
+    }
+};
+
+__END__

xs/MANIFEST

@@ -103,6 +103,8 @@ src/libslic3r/Layer.cpp
 src/libslic3r/Layer.hpp
 src/libslic3r/LayerRegion.cpp
 src/libslic3r/LayerRegionFill.cpp
+src/libslic3r/LayerHeightSpline.hpp
+src/libslic3r/LayerHeightSpline.cpp
 src/libslic3r/libslic3r.h
 src/libslic3r/Line.cpp
 src/libslic3r/Line.hpp
@@ -215,6 +217,7 @@ xsp/Geometry.xsp
 xsp/GUI.xsp
 xsp/GUI_3DScene.xsp
 xsp/Layer.xsp
+xsp/LayerHeightSpline.xsp
 xsp/Line.xsp
 xsp/Model.xsp
 xsp/MotionPlanner.xsp

xs/lib/Slic3r/XS.pm

@@ -240,6 +240,7 @@ for my $class (qw(
         Slic3r::Layer
         Slic3r::Layer::Region
         Slic3r::Layer::Support
+        Slic3r::LayerHeightSpline
         Slic3r::Line
         Slic3r::Linef3
         Slic3r::Model
@@ -258,6 +259,7 @@ for my $class (qw(
         Slic3r::Print::Object
         Slic3r::Print::Region
         Slic3r::Print::State
+        Slic3r::SlicingAdaptive
         Slic3r::Surface
         Slic3r::Surface::Collection
         Slic3r::TriangleMesh

xs/src/BSpline/BSpline.cpp

@@ -0,0 +1,879 @@
+// -*- mode: c++; c-basic-offset: 4; -*-
+/************************************************************************
+ * Copyright 2009 University Corporation for Atmospheric Research.
+ * All rights reserved.
+ *
+ * Use of this code is subject to the standard BSD license:
+ *
+ *  http://www.opensource.org/licenses/bsd-license.html
+ *
+ * See the COPYRIGHT file in the source distribution for the license text,
+ * or see this web page:
+ *
+ *  http://www.eol.ucar.edu/homes/granger/bspline/doc/
+ *
+ *************************************************************************/
+
+/**
+ * @file
+ *
+ * This file defines the implementation for the BSpline and BSplineBase
+ * templates.
+ **/
+#include "BSpline.h"
+#include "BandedMatrix.h"
+
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <iostream>
+#include <iomanip>
+#include <map>
+#include <assert.h>
+
+
+
+
+
+/*
+ * Original BSplineBase.cpp start here
+ */
+
+/*
+ * This class simulates a namespace for private symbols used by this template
+ * implementation which should not pollute the global namespace.
+ */
+class my
+{
+    public:
+        template<class T> static inline
+        T abs(const T t)
+        {
+            return (t < 0) ? -t : t;
+        }
+
+        template<class T> static inline
+        const T& min(const T& a,
+                     const T& b)
+        {
+            return (a < b) ? a : b;
+        }
+
+        template<class T> static inline
+        const T& max(const T& a,
+                     const T& b)
+        {
+            return (a > b) ? a : b;
+        }
+};
+
+//////////////////////////////////////////////////////////////////////
+template<class T> class Matrix : public BandedMatrix<T>
+{
+    public:
+        Matrix &operator +=(const Matrix &B)
+        {
+            Matrix &A = *this;
+            typename Matrix::size_type M = A.num_rows();
+            typename Matrix::size_type N = A.num_cols();
+
+            assert(M==B.num_rows());
+            assert(N==B.num_cols());
+
+            typename Matrix::size_type i, j;
+            for (i=0; i<M; i++)
+                for (j=0; j<N; j++)
+                    A[i][j] += B[i][j];
+            return A;
+        }
+
+        inline Matrix & operator=(const Matrix &b)
+        {
+            return Copy(*this, b);
+        }
+
+        inline Matrix & operator=(const T &e)
+        {
+            BandedMatrix<T>::operator= (e);
+            return *this;
+        }
+
+    };
+    //////////////////////////////////////////////////////////////////////
+    // Our private state structure, which hides our use of some matrix
+    // template classes.
+
+template<class T> struct BSplineBaseP
+{
+        typedef Matrix<T> MatrixT;
+
+        MatrixT Q; // Holds P+Q and its factorization
+        std::vector<T> X;
+        std::vector<T> Nodes;
+};
+
+//////////////////////////////////////////////////////////////////////
+
+// This array contains the beta parameter for the boundary condition
+// constraints.  The boundary condition type--0, 1, or 2--is the first
+// index into the array, followed by the index of the endpoints.  See the
+// Beta() method.
+template<class T> const double BSplineBase<T>::BoundaryConditions[3][4] =
+    {
+            //  0       1       M-1     M
+                {
+                        -4,
+                        -1,
+                        -1,
+                        -4 },
+                {
+                        0,
+                        1,
+                        1,
+                        0 },
+                {
+                        2,
+                        -1,
+                        -1,
+                        2 } };
+//////////////////////////////////////////////////////////////////////
+template<class T> inline bool BSplineBase<T>::Debug(int on)
+{
+    static bool debug = false;
+    if (on >= 0)
+        debug = (on > 0);
+    return debug;
+}
+
+//////////////////////////////////////////////////////////////////////
+template<class T> const double BSplineBase<T>::BS_PI = 3.1415927;
+//////////////////////////////////////////////////////////////////////
+template<class T> const char * BSplineBase<T>::ImplVersion()
+{
+    return ("$Id: BSpline.cpp 6352 2008-05-05 04:40:39Z martinc $");
+}
+
+//////////////////////////////////////////////////////////////////////
+template<class T> const char * BSplineBase<T>::IfaceVersion()
+{
+    return (_BSPLINEBASE_IFACE_ID);
+}
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+template<class T> BSplineBase<T>::~BSplineBase()
+{
+    delete base;
+}
+
+// This is a member-wise copy except for replacing our
+// private base structure with the source's, rather than just copying
+// the pointer.  But we use the compiler's default copy constructor for
+// constructing our BSplineBaseP.
+template<class T> BSplineBase<T>::BSplineBase(const BSplineBase<T> &bb) :
+    K(bb.K), BC(bb.BC), OK(bb.OK), base(new BSplineBaseP<T>(*bb.base))
+{
+    xmin = bb.xmin;
+    xmax = bb.xmax;
+    alpha = bb.alpha;
+    waveLength = bb.waveLength;
+    DX = bb.DX;
+    M = bb.M;
+    NX = base->X.size();
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> BSplineBase<T>::BSplineBase(const T *x,
+                                              int nx,
+                                              double wl,
+                                              int bc,
+                                              int num_nodes) :
+    NX(0), K(2), OK(false), base(new BSplineBaseP<T>)
+{
+    setDomain(x, nx, wl, bc, num_nodes);
+}
+
+//////////////////////////////////////////////////////////////////////
+// Methods
+template<class T> bool BSplineBase<T>::setDomain(const T *x,
+                                                 int nx,
+                                                 double wl,
+                                                 int bc,
+                                                 int num_nodes)
+{
+    if ((nx <= 0) || (x == 0) || (wl< 0) || (bc< 0) || (bc> 2)) {
+        return false;
+    }
+    OK = false;
+    waveLength = wl;
+    BC = bc;
+    // Copy the x array into our storage.
+    base->X.resize(nx);
+    std::copy(x, x+nx, base->X.begin());
+    NX = base->X.size();
+
+    // The Setup() method determines the number and size of node intervals.
+    if (Setup(num_nodes)) {
+        if (Debug()) {
+            std::cerr << "Using M node intervals: " << M << " of length DX: "
+                    << DX << std::endl;
+            std::cerr << "X min: " << xmin << " ; X max: " << xmax << std::endl;
+            std::cerr << "Data points per interval: " << (float)NX/(float)M
+                    << std::endl;
+            std::cerr << "Nodes per wavelength: " << (float)waveLength
+                    /(float)DX << std::endl;
+            std::cerr << "Derivative constraint degree: " << K << std::endl;
+        }
+
+        // Now we can calculate alpha and our Q matrix
+        alpha = Alpha(waveLength);
+        if (Debug()) {
+            std::cerr << "Cutoff wavelength: " << waveLength << " ; "
+                    << "Alpha: " << alpha << std::endl;
+            std::cerr << "Calculating Q..." << std::endl;
+        }
+        calculateQ();
+        if (Debug() && M < 30) {
+            std::cerr.fill(' ');
+            std::cerr.precision(2);
+            std::cerr.width(5);
+            std::cerr << base->Q << std::endl;
+        }
+
+        if (Debug())
+            std::cerr << "Calculating P..." << std::endl;
+        addP();
+        if (Debug()) {
+            std::cerr << "Done." << std::endl;
+            if (M < 30) {
+                std::cerr << "Array Q after addition of P." << std::endl;
+                std::cerr << base->Q;
+            }
+        }
+
+        // Now perform the LU factorization on Q
+        if (Debug())
+            std::cerr << "Beginning LU factoring of P+Q..." << std::endl;
+        if (!factor()) {
+            if (Debug())
+                std::cerr << "Factoring failed." << std::endl;
+        } else {
+            if (Debug())
+                std::cerr << "Done." << std::endl;
+            OK = true;
+        }
+    }
+    return OK;
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Calculate the alpha parameter given a wavelength.
+ */
+template<class T> double BSplineBase<T>::Alpha(double wl)
+{
+    // K is the degree of the derivative constraint: 1, 2, or 3
+    double a = (double) (wl / (2 * BS_PI * DX));
+    a *= a; // a^2
+    if (K == 2)
+        a = a * a; // a^4
+    else if (K == 3)
+        a = a * a * a; // a^6
+    return a;
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Return the correct beta value given the node index.  The value depends
+ * on the node index and the current boundary condition type.
+ */
+template<class T> inline double BSplineBase<T>::Beta(int m)
+{
+    if (m > 1 && m < M-1)
+        return 0.0;
+    if (m >= M-1)
+        m -= M-3;
+    assert(0 <= BC && BC <= 2);
+    assert(0 <= m && m <= 3);
+    return BoundaryConditions[BC][m];
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Given an array of y data points defined over the domain
+ * of x data points in this BSplineBase, create a BSpline
+ * object which contains the smoothed curve for the y array.
+ */
+template<class T> BSpline<T> * BSplineBase<T>::apply(const T *y)
+{
+    return new BSpline<T> (*this, y);
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Evaluate the closed basis function at node m for value x,
+ * using the parameters for the current boundary conditions.
+ */
+template<class T> double BSplineBase<T>::Basis(int m,
+                                               T x)
+{
+    double y = 0;
+    double xm = xmin + (m * DX);
+    double z = my::abs((double)(x - xm) / (double)DX);
+    if (z < 2.0) {
+        z = 2 - z;
+        y = 0.25 * (z*z*z);
+        z -= 1.0;
+        if (z > 0)
+            y -= (z*z*z);
+    }
+
+    // Boundary conditions, if any, are an additional addend.
+    if (m == 0 || m == 1)
+        y += Beta(m) * Basis(-1, x);
+    else if (m == M-1 || m == M)
+        y += Beta(m) * Basis(M+1, x);
+
+    return y;
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Evaluate the deriviative of the closed basis function at node m for
+ * value x, using the parameters for the current boundary conditions.
+ */
+template<class T> double BSplineBase<T>::DBasis(int m,
+                                                T x)
+{
+    double dy = 0;
+    double xm = xmin + (m * DX);
+    double delta = (double)(x - xm) / (double)DX;
+    double z = my::abs(delta);
+    if (z < 2.0) {
+        z = 2.0 - z;
+        dy = 0.25 * z * z;
+        z -= 1.0;
+
+        if (z > 0) {
+            dy -= z * z;
+        }
+        dy *= ((delta > 0) ? -1.0 : 1.0) * 3.0 / DX;
+    }
+
+    // Boundary conditions, if any, are an additional addend.
+    if (m == 0 || m == 1)
+        dy += Beta(m) * DBasis(-1, x);
+    else if (m == M-1 || m == M)
+        dy += Beta(m) * DBasis(M+1, x);
+
+    return dy;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> double BSplineBase<T>::qDelta(int m1,
+                                                int m2)
+/*
+ * Return the integral of the product of the basis function derivative
+ * restricted to the node domain, 0 to M.
+ */
+{
+    // These are the products of the Kth derivative of the
+    // normalized basis functions
+    // given a distance m nodes apart, qparts[K-1][m], 0 <= m <= 3
+    // Each column is the integral over each unit domain, -2 to 2
+    static const double qparts[3][4][4] =
+        {
+                    {
+                                {
+                                        0.11250,
+                                        0.63750,
+                                        0.63750,
+                                        0.11250 },
+                                {
+                                        0.00000,
+                                        0.13125,
+                                        -0.54375,
+                                        0.13125 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        -0.22500,
+                                        -0.22500 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        0.00000,
+                                        -0.01875 } },
+                    {
+                                {
+                                        0.75000,
+                                        2.25000,
+                                        2.25000,
+                                        0.75000 },
+                                {
+                                        0.00000,
+                                        -1.12500,
+                                        -1.12500,
+                                        -1.12500 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        0.00000,
+                                        0.00000 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        0.00000,
+                                        0.37500 } },
+                    {
+                                {
+                                        2.25000,
+                                        20.25000,
+                                        20.25000,
+                                        2.25000 },
+                                {
+                                        0.00000,
+                                        -6.75000,
+                                        -20.25000,
+                                        -6.75000 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        6.75000,
+                                        6.75000 },
+                                {
+                                        0.00000,
+                                        0.00000,
+                                        0.00000,
+                                        -2.25000 } } };
+
+    if (m1 > m2)
+        std::swap(m1, m2);
+
+    if (m2 - m1 > 3)
+        return 0.0;
+
+    double q = 0;
+    for (int m = my::max(m1-2, 0); m < my::min(m1+2, M); ++m)
+        q += qparts[K-1][m2-m1][m-m1+2];
+    return q * alpha;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> void BSplineBase<T>::calculateQ()
+{
+    Matrix<T> &Q = base->Q;
+    Q.setup(M+1, 3);
+    Q = 0;
+    if (alpha == 0)
+        return;
+
+    // First fill in the q values without the boundary constraints.
+    int i;
+    for (i = 0; i <= M; ++i) {
+        Q[i][i] = qDelta(i, i);
+        for (int j = 1; j < 4 && i+j <= M; ++j) {
+            Q[i][i+j] = Q[i+j][i] = qDelta(i, i+j);
+        }
+    }
+
+    // Now add the boundary constraints:
+    // First the upper left corner.
+    float b1, b2, q;
+    for (i = 0; i <= 1; ++i) {
+        b1 = Beta(i);
+        for (int j = i; j < i+4; ++j) {
+            b2 = Beta(j);
+            assert(j-i >= 0 && j - i < 4);
+            q = 0.0;
+            if (i+1 < 4)
+                q += b2*qDelta(-1, i);
+            if (j+1 < 4)
+                q += b1*qDelta(-1, j);
+            q += b1*b2*qDelta(-1, -1);
+            Q[j][i] = (Q[i][j] += q);
+        }
+    }
+
+    // Then the lower right
+    for (i = M-1; i <= M; ++i) {
+        b1 = Beta(i);
+        for (int j = i - 3; j <= i; ++j) {
+            b2 = Beta(j);
+            q = 0.0;
+            if (M+1-i < 4)
+                q += b2*qDelta(i, M+1);
+            if (M+1-j < 4)
+                q += b1*qDelta(j, M+1);
+            q += b1*b2*qDelta(M+1, M+1);
+            Q[j][i] = (Q[i][j] += q);
+        }
+    }
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> void BSplineBase<T>::addP()
+{
+    // Add directly to Q's elements
+    Matrix<T> &P = base->Q;
+    std::vector<T> &X = base->X;
+
+    // For each data point, sum the product of the nearest, non-zero Basis
+    // nodes
+    int mx, m, n, i;
+    for (i = 0; i < NX; ++i) {
+        // Which node does this put us in?
+        T &x = X[i];
+        mx = (int)((x - xmin) / DX);
+
+        // Loop over the upper triangle of nonzero basis functions,
+        // and add in the products on each side of the diagonal.
+        for (m = my::max(0, mx-1); m <= my::min(M, mx+2); ++m) {
+            float pn;
+            float pm = Basis(m, x);
+            float sum = pm * pm;
+            P[m][m] += sum;
+            for (n = m+1; n <= my::min(M, mx+2); ++n) {
+                pn = Basis(n, x);
+                sum = pm * pn;
+                P[m][n] += sum;
+                P[n][m] += sum;
+            }
+        }
+    }
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> bool BSplineBase<T>::factor()
+{
+    Matrix<T> &LU = base->Q;
+
+    if (LU_factor_banded(LU, 3) != 0) {
+        if (Debug())
+            std::cerr << "LU_factor_banded() failed." << std::endl;
+        return false;
+    }
+    if (Debug() && M < 30)
+        std::cerr << "LU decomposition: " << std::endl << LU << std::endl;
+    return true;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> inline double BSplineBase<T>::Ratiod(int &ni,
+                                                       double &deltax,
+                                                       double &ratiof)
+{
+    deltax = (xmax - xmin) / ni;
+    ratiof = waveLength / deltax;
+    double ratiod = (double) NX / (double) (ni + 1);
+    return ratiod;
+}
+//////////////////////////////////////////////////////////////////////
+// Setup the number of nodes (and hence deltax) for the given domain and
+// cutoff wavelength.  According to Ooyama, the derivative constraint
+// approximates a lo-pass filter if the cutoff wavelength is about 4*deltax
+// or more, but it should at least be 2*deltax.  We can increase the number
+// of nodes to increase the number of nodes per cutoff wavelength.
+// However, to get a reasonable representation of the data, the setup
+// enforces at least as many nodes as data points in the domain.  (This
+// constraint assumes reasonably even distribution of data points, since
+// its really the density of data points which matters.)
+//
+// Return zero if the setup fails, non-zero otherwise.
+//
+// The algorithm in this routine is mostly taken from the FORTRAN
+// implementation by James Franklin, NOAA/HRD.
+//
+template<class T> bool BSplineBase<T>::Setup(int num_nodes)
+{
+    std::vector<T> &X = base->X;
+
+    // Find the min and max of the x domain
+    xmin = X[0];
+    xmax = X[0];
+
+    int i;
+    for (i = 1; i < NX; ++i) {
+        if (X[i] < xmin)
+            xmin = X[i];
+        else if (X[i] > xmax)
+            xmax = X[i];
+    }
+    if (Debug())
+	std::cerr << "Xmax=" << xmax << ", Xmin=" << xmin << std::endl;
+
+    // Number of node intervals (number of spline nodes - 1).
+    int ni;
+    double deltax;
+
+    if (num_nodes >= 2) {
+        // We've been told explicitly the number of nodes to use.
+        ni = num_nodes - 1;
+        if (waveLength == 0) {
+            waveLength = 1.0;
+        }
+	if (Debug())
+	{
+	    std::cerr << "Num nodes explicitly given as " << num_nodes
+		      << ", wavelength set to " << waveLength << std::endl;
+	}
+    } else if (waveLength == 0) {
+        // Turn off frequency constraint and just set two node intervals per
+        // data point.
+        ni = NX * 2;
+        waveLength = 1;
+	if (Debug())
+	{
+	    std::cerr << "Frequency constraint disabled, using 2 intervals "
+		      << "per node: " << ni << " intervals, wavelength="
+		      << waveLength << std::endl;
+	}
+    } else if (waveLength > xmax - xmin) {
+	if (Debug())
+	{
+	    std::cerr << "Wavelength " << waveLength << " exceeds X span: "
+		      << xmin << " - " << xmax << std::endl;
+	}
+        return (false);
+    } else {
+	if (Debug())
+	{
+	    std::cerr << "Searching for a reasonable number of "
+		      << "intervals for wavelength " << waveLength
+		      << " while keeping at least 2 intervals per "
+		      << "wavelength ..." << std::endl;
+	}
+        // Minimum acceptable number of node intervals per cutoff wavelength.
+        static const double fmin = 2.0;
+
+	// Start out at a minimum number of intervals, meaning the maximum
+	// number of points per interval, then work up to the maximum
+	// number of intervals for which the intervals per wavelength is
+	// still adequate.  I think the minimum must be more than 2 since
+	// the basis function is evaluated on multiple nodes.
+	ni = 5;
+
+        double ratiof; // Nodes per wavelength for current deltax
+        double ratiod; // Points per node interval
+
+        // Increase the number of node intervals until we reach the minimum
+        // number of intervals per cutoff wavelength, but only as long as
+        // we can maintain at least one point per interval.
+        do {
+            if (Ratiod(++ni, deltax, ratiof) < 1.0)
+	    {
+		if (Debug())
+		{
+		    std::cerr << "At " << ni << " intervals, fewer than "
+			      << "one point per interval, and "
+			      << "intervals per wavelength is "
+			      << ratiof << "." << std::endl;
+		}
+                return false;
+	    }
+        } while (ratiof < fmin);
+
+        // Now increase the number of intervals until we have at least 4
+        // intervals per cutoff wavelength, but only as long as we can
+        // maintain at least 2 points per node interval.  There's also no
+        // point to increasing the number of intervals if we already have
+        // 15 or more nodes per cutoff wavelength.
+        //
+        do {
+            if ((ratiod = Ratiod(++ni, deltax, ratiof)) < 1.0 || ratiof > 15.0) {
+                --ni;
+                break;
+            }
+        } while (ratiof < 4 || ratiod > 2.0);
+
+	if (Debug())
+	{
+	    std::cerr << "Found " << ni << " intervals, "
+		      << "length " << deltax << ", "
+		      << ratiof << " nodes per wavelength " << waveLength
+		      << ", "
+		      << ratiod << " data points per interval." << std::endl;
+	}
+    }
+
+    // Store the calculations in our state
+    M = ni;
+    DX = (xmax - xmin) / ni;
+
+    return (true);
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> const T * BSplineBase<T>::nodes(int *nn)
+{
+    if (base->Nodes.size() == 0) {
+        base->Nodes.reserve(M+1);
+        for (int i = 0; i <= M; ++i) {
+            base->Nodes.push_back(xmin + (i * DX));
+        }
+    }
+
+    if (nn)
+        *nn = base->Nodes.size();
+
+    assert(base->Nodes.size() == (unsigned)(M+1));
+    return &base->Nodes[0];
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> std::ostream &operator<<(std::ostream &out,
+                                           const std::vector<T> &c)
+{
+    for (typename std::vector<T>::const_iterator it = c.begin(); it < c.end(); ++it)
+        out << *it << ", ";
+    out << std::endl;
+    return out;
+}
+
+
+
+
+
+/*
+ * Original BSpline.cpp start here
+ */
+
+
+//////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////
+// BSpline Class
+//////////////////////////////////////////////////////////////////////
+
+template<class T> struct BSplineP {
+        std::vector<T> spline;
+        std::vector<T> A;
+};
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+/*
+ * This BSpline constructor constructs and sets up a new base and
+ * solves for the spline curve coeffiecients all at once.
+ */
+template<class T> BSpline<T>::BSpline(const T *x,
+                                      int nx,
+                                      const T *y,
+                                      double wl,
+                                      int bc_type,
+                                      int num_nodes) :
+    BSplineBase<T>(x, nx, wl, bc_type, num_nodes), s(new BSplineP<T>) {
+    solve(y);
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * Create a new spline given a BSplineBase.
+ */
+template<class T> BSpline<T>::BSpline(BSplineBase<T> &bb,
+                                      const T *y) :
+    BSplineBase<T>(bb), s(new BSplineP<T>) {
+    solve(y);
+}
+//////////////////////////////////////////////////////////////////////
+/*
+ * (Re)calculate the spline for the given set of y values.
+ */
+template<class T> bool BSpline<T>::solve(const T *y) {
+    if (!OK)
+        return false;
+
+    // Any previously calculated curve is now invalid.
+    s->spline.clear();
+    OK = false;
+
+    // Given an array of data points over x and its precalculated
+    // P+Q matrix, calculate the b vector and solve for the coefficients.
+    std::vector<T> &B = s->A;
+    std::vector<T> &A = s->A;
+    A.clear();
+    A.resize(M+1);
+
+    if (Debug())
+        std::cerr << "Solving for B..." << std::endl;
+
+    // Find the mean of these data
+    mean = 0.0;
+    int i;
+    for (i = 0; i < NX; ++i) {
+        mean += y[i];
+    }
+    mean = mean / (double)NX;
+    if (Debug())
+        std::cerr << "Mean for y: " << mean << std::endl;
+
+    int mx, m, j;
+    for (j = 0; j < NX; ++j) {
+        // Which node does this put us in?
+        T &xj = base->X[j];
+        T yj = y[j] - mean;
+        mx = (int)((xj - xmin) / DX);
+
+        for (m = my::max(0, mx-1); m <= my::min(mx+2, M); ++m) {
+            B[m] += yj * Basis(m, xj);
+        }
+    }
+
+    if (Debug() && M < 30) {
+        std::cerr << "Solution a for (P+Q)a = b" << std::endl;
+        std::cerr << " b: " << B << std::endl;
+    }
+
+    // Now solve for the A vector in place.
+    if (LU_solve_banded(base->Q, A, 3) != 0) {
+        if (Debug())
+            std::cerr << "LU_solve_banded() failed." << std::endl;
+    } else {
+        OK = true;
+        if (Debug())
+            std::cerr << "Done." << std::endl;
+        if (Debug() && M < 30) {
+            std::cerr << " a: " << A << std::endl;
+            std::cerr << "LU factor of (P+Q) = " << std::endl << base->Q
+                    << std::endl;
+        }
+    }
+    return (OK);
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> BSpline<T>::~BSpline() {
+    delete s;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> T BSpline<T>::coefficient(int n) {
+    if (OK)
+        if (0 <= n && n <= M)
+            return s->A[n];
+    return 0;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> T BSpline<T>::evaluate(T x) {
+    T y = 0;
+    if (OK) {
+        int n = (int)((x - xmin)/DX);
+        for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) {
+            y += s->A[i] * Basis(i, x);
+        }
+        y += mean;
+    }
+    return y;
+}
+//////////////////////////////////////////////////////////////////////
+template<class T> T BSpline<T>::slope(T x) {
+    T dy = 0;
+    if (OK) {
+        int n = (int)((x - xmin)/DX);
+        for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) {
+            dy += s->A[i] * DBasis(i, x);
+        }
+    }
+    return dy;
+}
+
+
+
+
+
+
+
+/// Instantiate BSplineBase
+template class  BSplineBase<double>;
+//template class  BSplineBase<float>;
+
+/// Instantiate BSpline
+template class BSpline<double>;
+//template class BSpline<float>;

xs/src/BSpline/BSpline.h

@@ -0,0 +1,452 @@
+/* -*- mode: c++; c-basic-offset: 4; -*- */
+/************************************************************************
+ * Copyright 2009 University Corporation for Atmospheric Research.
+ * All rights reserved.
+ *
+ * Use of this code is subject to the standard BSD license:
+ *
+ *  http://www.opensource.org/licenses/bsd-license.html
+ *
+ * See the COPYRIGHT file in the source distribution for the license text,
+ * or see this web page:
+ *
+ *  http://www.eol.ucar.edu/homes/granger/bspline/doc/
+ *
+ *************************************************************************/
+
+#ifndef BSPLINE_H
+#define BSPLINE_H
+
+#include <vector>
+
+
+/*
+ * Warning: original BSplineBase.h/cpp merged into BSpline.h/cpp to avoid dependency issues caused by Build::WithXSpp which tries to compile all .cpp files in /src
+ * Original BSplineBase.h starts here!!!
+ */
+
+
+
+#ifndef _BSPLINEBASE_IFACE_ID
+#define _BSPLINEBASE_IFACE_ID "$Id: BSpline.h 6353 2008-05-05 19:30:48Z martinc $"
+#endif
+
+/**
+ * @file
+ *
+ * This file defines the BSpline library interface.
+ *
+ */
+template <class T> class BSpline;
+
+/*
+ * Opaque member structure to hide the matrix implementation.
+ */
+template <class T> struct BSplineBaseP;
+
+/**
+ * @class BSplineBase
+ *
+ * The base class for a spline object containing the nodes for a given
+ * domain, cutoff wavelength, and boundary condition.
+
+ * To smooth a single curve, the BSpline interface contains a constructor
+ * which both sets up the domain and solves for the spline.  Subsequent
+ * curves over the same domain can be created by apply()ing them to the
+ * BSpline object, where apply() is a BSplineBase method.  [See apply().]
+ * New curves can also be smoothed within the same BSpline object by
+ * calling solve() with the new set of y values.  [See BSpline.]  A
+ * BSplineBase can be created on its own, in which case all of the
+ * computations dependent on the x values, boundary conditions, and cutoff
+ * wavelength have already been completed.
+ *
+ * The solution of the cubic b-spline is divided into two parts.  The first
+ * is the setup of the domain given the x values, boundary conditions, and
+ * wavelength.  The second is the solution of the spline for a set of y
+ * values corresponding to the x values in the domain.  The first part is
+ * done in the creation of the BSplineBase object (or when calling the
+ * setDomain method).  The second part is done when creating a BSpline
+ * object (or calling solve() on a BSpline object).
+ *
+ * A BSpline object can be created with either one of its constructors, or
+ * by calling apply() on an existing BSplineBase object.  Once a spline has
+ * been solved, it can be evaluated at any x value.  The following example
+ * creates a spline curve and evaluates it over the domain:
+ *
+@verbatim
+
+    vector<float> x;
+    vector<float> y;
+    { ... }
+    int bc = BSplineBase<float>::BC_ZERO_SECOND;
+    BSpline<float>::Debug = true;
+    BSpline<float> spline (x.begin(), x.size(), y.begin(), wl, bc);
+    if (spline.ok())
+    {
+        ostream_iterator<float> of(cout, "\t ");
+        float xi = spline.Xmin();
+    float xs = (spline.Xmax() - xi) / 2000.0;
+    for (; xi <= spline.Xmax(); xi += xs)
+    {
+        *of++ = spline.evaluate (xi);
+    }
+    }
+
+@endverbatim
+ *
+ * In the usual usage, the BSplineBase can compute a reasonable number of
+ * nodes for the spline, balancing between a few desirable factors.  There
+ * needs to be at least 2 nodes per cutoff wavelength (preferably 4 or
+ * more) for the derivative constraint to reliably approximate a lo-pass
+ * filter.  There should be at least 1 and preferably about 2 data points
+ * per node (measured just by their number and not by any check of the
+ * density of points across the domain).  Lastly, of course, the fewer the
+ * nodes then the faster the computation of the spline.  The computation of
+ * the number of nodes happens in the Setup() method during BSplineBase
+ * construction and when setDomain() is called.  If the setup fails to find
+ * a desirable number of nodes, then the BSplineBase object will return
+ * false from ok().
+ *
+ * The ok() method returns false when a BSplineBase or BSpline could not
+ * complete any operation successfully.  In particular, as mentioned above,
+ * ok() will return false if some problem was detected with the domain
+ * values or if no reasonable number of nodes could be found for the given
+ * cutoff wavelength.  Also, ok() on a BSpline object will return false if
+ * the matrix equation could not be solved, such as after BSpline
+ * construction or after a call to apply().
+ *
+ * If letting Setup() determine the number of nodes is not acceptable, the
+ * constructors and setDomain() accept the parameter num_nodes.  By
+ * default, num_nodes is passed as zero, forcing Setup() to calculate the
+ * number of nodes.  However, if num_nodes is passed as 2 or greater, then
+ * Setup() will bypass its own algorithm and accept the given number of
+ * nodes instead.  Obviously, it's up to the programmer to understand the
+ * affects of the number of nodes on the representation of the data and on
+ * the solution (or non-solution) of the spline.  Remember to check the
+ * ok() method to detect when the spline solution has failed.
+ *
+ * The interface for the BSplineBase and BSpline templates is defined in
+ * the header file BSpline.h.  The implementation is defined in BSpline.cpp.
+ * Source files which will instantiate the template should include the
+ * implementation file and @em not the interface.  If the implementation
+ * for a specific type will be linked from elsewhere, such as a
+ * static library or Windows DLL, source files should only include the
+ * interface file.  On Windows, applications should link with the import
+ * library BSpline.lib and make sure BSpline.dll is on the path.  The DLL
+ * contains an implementation for BSpline<float> and BSpline<double>.
+ * For debugging, an application can include the implementation to get its
+ * own instantiation.
+ *
+ * The algorithm is based on the cubic spline described by Katsuyuki Ooyama
+ * in Montly Weather Review, Vol 115, October 1987.  This implementation
+ * has benefited from comparisons with a previous FORTRAN implementation by
+ * James L. Franklin, NOAA/Hurricane Research Division.  In particular, the
+ * algorithm in the Setup() method is based mostly on his implementation
+ * (VICSETUP).  The Setup() method finds a suitable default for the number
+ * of nodes given a domain and cutoff frequency.  This implementation
+ * adopts most of the same constraints, including a constraint that the
+ * cutoff wavelength not be greater than the span of the domain values: wl
+ * < max(x) - min(x). If this is not an acceptable constraint, then use the
+ * num_nodes parameter to specify the number of nodes explicitly.
+ *
+ * The cubic b-spline is formulated as the sum of some multiple of the
+ * basis function centered at each node in the domain.  The number of nodes
+ * is determined by the desired cutoff wavelength and a desirable number of
+ * x values per node.  The basis function is continuous and differentiable
+ * up to the second degree.  A derivative constraint is included in the
+ * solution to achieve the effect of a low-pass frequency filter with the
+ * given cutoff wavelength.  The derivative constraint can be disabled by
+ * specifying a wavelength value of zero, which reduces the analysis to a
+ * least squares fit to a cubic b-spline.  The domain nodes, boundary
+ * constraints, and wavelength determine a linear system of equations,
+ * Qa=b, where a is the vector of basis function coefficients at each node.
+ * The coefficient vector is solved by first LU factoring along the
+ * diagonally banded matrix Q in BSplineBase.  The BSpline object then
+ * computes the B vector for a set of y values and solves for the
+ * coefficient vector with the LU matrix.  Only the diagonal bands are
+ * stored in memory and calculated during LU factoring and back
+ * substitution, and the basis function is evaluated as few times as
+ * possible in computing the diagonal matrix and B vector.
+ *
+ * @author Gary Granger (http://www.eol.ucar.edu/homes/granger)
+ *
+@verbatim
+Copyright (c) 1998-2009
+University Corporation for Atmospheric Research, UCAR
+All rights reserved.
+@endverbatim
+ **/
+
+template <class T>
+class BSplineBase
+{
+public:
+    // Datum type
+    typedef T datum_type;
+
+    /// Return a string describing the implementation version.
+    static const char *ImplVersion();
+
+    /// Return a string describing the interface version.
+    static const char *IfaceVersion();
+
+    /**
+     * Call this class method with a value greater than zero to enable
+     * debug messages, or with zero to disable messages.  Calling with
+     * no arguments returns true if debugging enabled, else false.
+     */
+    static bool Debug (int on = -1);
+
+    /**
+     * Boundary condition types.
+     */
+    enum BoundaryConditionTypes
+    {
+    /// Set the endpoints of the spline to zero.
+    BC_ZERO_ENDPOINTS = 0,
+    /// Set the first derivative of the spline to zero at the endpoints.
+    BC_ZERO_FIRST = 1,
+    /// Set the second derivative to zero.
+    BC_ZERO_SECOND = 2
+    };
+
+public:
+
+    /**
+     * Construct a spline domain for the given set of x values, cutoff
+     * wavelength, and boundary condition type.  The parameters are the
+     * same as for setDomain().  Call ok() to check whether domain
+     * setup succeeded after construction.
+     */
+    BSplineBase (const T *x, int nx,
+         double wl, int bc_type = BC_ZERO_SECOND,
+         int num_nodes = 0);
+
+    /// Copy constructor
+    BSplineBase (const BSplineBase &);
+
+    /**
+     * Change the domain of this base.  [If this is part of a BSpline
+     * object, this method {\em will not} change the existing curve or
+     * re-apply the smoothing to any set of y values.]
+     *
+     * The x values can be in any order, but they must be of sufficient
+     * density to support the requested cutoff wavelength.  The setup of
+     * the domain may fail because of either inconsistency between the x
+     * density and the cutoff wavelength, or because the resulting matrix
+     * could not be factored.  If setup fails, the method returns false.
+     *
+     * @param x     The array of x values in the domain.
+     * @param nx    The number of values in the @p x array.
+     * @param wl    The cutoff wavelength, in the same units as the
+     *          @p x values.  A wavelength of zero disables
+     *          the derivative constraint.
+     * @param bc_type   The enumerated boundary condition type.  If
+     *          omitted it defaults to BC_ZERO_SECOND.
+     * @param num_nodes The number of nodes to use for the cubic b-spline.
+     *          If less than 2 a reasonable number will be
+     *          calculated automatically, if possible, taking
+     *          into account the given cutoff wavelength.
+     *
+     * @see ok().
+     */
+    bool setDomain (const T *x, int nx, double wl,
+            int bc_type = BC_ZERO_SECOND,
+            int num_nodes = 0);
+
+    /**
+     * Create a BSpline smoothed curve for the given set of NX y values.
+     * The returned object will need to be deleted by the caller.
+     * @param y The array of y values corresponding to each of the nX()
+     *      x values in the domain.
+     * @see ok()
+     */
+    BSpline<T> *apply (const T *y);
+
+    /**
+     * Return array of the node coordinates.  Returns 0 if not ok().  The
+     * array of nodes returned by nodes() belongs to the object and should
+     * not be deleted; it will also be invalid if the object is destroyed.
+     */
+    const T *nodes (int *nnodes);
+
+    /**
+     * Return the number of nodes (one more than the number of intervals).
+     */
+    int nNodes () { return M+1; }
+
+    /**
+     * Number of original x values.
+     */
+    int nX () { return NX; }
+
+    /// Minimum x value found.
+    T Xmin () { return xmin; }
+
+    /// Maximum x value found.
+    T Xmax () { return xmin + (M * DX); }
+
+    /**
+     * Return the Alpha value for a given wavelength.  Note that this
+     * depends on the current node interval length (DX).
+     */
+    double Alpha (double wavelength);
+
+    /**
+     * Return alpha currently in use by this domain.
+     */
+    double Alpha () { return alpha; }
+
+    /**
+     * Return the current state of the object, either ok or not ok.
+     * Use this method to test for valid state after construction or after
+     * a call to setDomain().  ok() will return false if either fail, such
+     * as when an appropriate number of nodes and node interval cannot be
+     * found for a given wavelength, or when the linear equation for the
+     * coefficients cannot be solved.
+     */
+    bool ok () { return OK; }
+
+    virtual ~BSplineBase();
+
+protected:
+
+    typedef BSplineBaseP<T> Base;
+
+    // Provided
+    double waveLength;  // Cutoff wavelength (l sub c)
+    int NX;
+    int K;  // Degree of derivative constraint (currently fixed at 2)
+    int BC;         // Boundary conditions type (0,1,2)
+
+    // Derived
+    T xmax;
+    T xmin;
+    int M;          // Number of intervals (M+1 nodes)
+    double DX;          // Interval length in same units as X
+    double alpha;
+    bool OK;
+    Base *base;         // Hide more complicated state members
+                    // from the public interface.
+
+    bool Setup (int num_nodes = 0);
+    void calculateQ ();
+    double qDelta (int m1, int m2);
+    double Beta (int m);
+    void addP ();
+    bool factor ();
+    double Basis (int m, T x);
+    double DBasis (int m, T x);
+
+    static const double BoundaryConditions[3][4];
+    static const double BS_PI;
+
+    double Ratiod (int&, double &, double &);
+};
+
+
+
+/*
+ * Original BSpline.h start here
+ */
+
+
+
+template <class T> struct BSplineP;
+
+
+/**
+ * Used to evaluate a BSpline.
+ * Inherits the BSplineBase domain information and interface and adds
+ * smoothing.  See the BSplineBase documentation for a summary of the
+ * BSpline interface.
+ */
+template <class T>
+class BSpline : public BSplineBase<T>
+{
+public:
+    /**
+     * Create a single spline with the parameters required to set up
+     * the domain and subsequently smooth the given set of y values.
+     * The y values must correspond to each of the values in the x array.
+     * If either the domain setup fails or the spline cannot be solved,
+     * the state will be set to not ok.
+     *
+     * @see ok().
+     *
+     * @param x		The array of x values in the domain.
+     * @param nx	The number of values in the @p x array.
+     * @param y		The array of y values corresponding to each of the
+     *			nX() x values in the domain.
+     * @param wl	The cutoff wavelength, in the same units as the
+     *			@p x values.  A wavelength of zero disables
+     *			the derivative constraint.
+     * @param bc_type	The enumerated boundary condition type.  If
+     *			omitted it defaults to BC_ZERO_SECOND.
+     * @param num_nodes The number of nodes to use for the cubic b-spline.
+     *			If less than 2 a "reasonable" number will be
+     *			calculated automatically, taking into account
+     *			the given cutoff wavelength.
+     */
+    BSpline (const T *x, int nx, 		/* independent variable */
+	     const T *y,			/* dependent values @ ea X */
+	     double wl,				/* cutoff wavelength */
+	     int bc_type = BSplineBase<T>::BC_ZERO_SECOND,
+	     int num_nodes = 0);
+
+    /**
+     * A BSpline curve can be derived from a separate @p base and a set
+     * of data points @p y over that base.
+     */
+    BSpline (BSplineBase<T> &base, const T *y);
+
+    /**
+     * Solve the spline curve for a new set of y values.  Returns false
+     * if the solution fails.
+     *
+     * @param y The array of y values corresponding to each of the nX()
+     *		x values in the domain.
+     */
+    bool solve (const T *y);
+
+    /**
+     * Return the evaluation of the smoothed curve
+     * at a particular @p x value.  If current state is not ok(), returns 0.
+     */
+    T evaluate (T x);
+
+    /**
+     * Return the first derivative of the spline curve at the given @p x.
+     * Returns zero if the current state is not ok().
+     */
+    T slope (T x);
+
+    /**
+     * Return the @p n-th basis coefficient, from 0 to M.  If the current
+     * state is not ok(), or @p n is out of range, the method returns zero.
+     */
+    T coefficient (int n);
+
+    virtual ~BSpline();
+
+    using BSplineBase<T>::Debug;
+    using BSplineBase<T>::Basis;
+    using BSplineBase<T>::DBasis;
+
+protected:
+
+    using BSplineBase<T>::OK;
+    using BSplineBase<T>::M;
+    using BSplineBase<T>::NX;
+    using BSplineBase<T>::DX;
+    using BSplineBase<T>::base;
+    using BSplineBase<T>::xmin;
+    using BSplineBase<T>::xmax;
+
+    // Our hidden state structure
+    BSplineP<T> *s;
+    T mean;			// Fit without mean and add it in later
+
+};
+
+#endif 

xs/src/BSpline/BandedMatrix.h

@@ -0,0 +1,375 @@
+/* -*- mode: c++; c-basic-offset: 4; -*- */
+/************************************************************************
+ * Copyright 2009 University Corporation for Atmospheric Research.
+ * All rights reserved.
+ *
+ * Use of this code is subject to the standard BSD license:
+ *
+ *  http://www.opensource.org/licenses/bsd-license.html
+ *
+ * See the COPYRIGHT file in the source distribution for the license text,
+ * or see this web page:
+ *
+ *  http://www.eol.ucar.edu/homes/granger/bspline/doc/
+ *
+ *************************************************************************/
+
+/**
+ * Template for a diagonally banded matrix.
+ **/
+#ifndef _BANDEDMATRIX_ID
+#define _BANDEDMATRIX_ID "$Id$"
+
+#include <vector>
+#include <algorithm>
+#include <iostream>
+
+template <class T> class BandedMatrixRow;
+
+
+template <class T> class BandedMatrix
+{
+public:
+    typedef unsigned int size_type;
+    typedef T element_type;
+
+    // Create a banded matrix with the same number of bands above and below
+    // the diagonal.
+    BandedMatrix (int N_ = 1, int nbands_off_diagonal = 0) : bands(0)
+    {
+	if (! setup (N_, nbands_off_diagonal))
+	    setup ();
+    }
+
+    // Create a banded matrix by naming the first and last non-zero bands,
+    // where the diagonal is at zero, and bands below the diagonal are
+    // negative, bands above the diagonal are positive.
+    BandedMatrix (int N_, int first, int last) : bands(0)
+    {
+	if (! setup (N_, first, last))
+	    setup ();
+    }
+
+    // Copy constructor
+    BandedMatrix (const BandedMatrix &b) : bands(0)
+    {
+	Copy (*this, b);
+    }
+
+    inline bool setup (int N_ = 1, int noff = 0)
+    {
+	return setup (N_, -noff, noff);
+    }
+
+    bool setup (int N_, int first, int last)
+    {
+	// Check our limits first and make sure they make sense.
+	// Don't change anything until we know it will work.
+	if (first > last || N_ <= 0)
+	    return false;
+
+	// Need at least as many N_ as columns and as rows in the bands.
+	if (N_ < abs(first) || N_ < abs(last))
+	    return false;
+
+	top = last;
+	bot = first;
+	N = N_;
+	out_of_bounds = T();
+
+	// Finally setup the diagonal vectors
+	nbands = last - first + 1;
+	if (bands) delete[] bands;
+	bands = new std::vector<T>[nbands];
+	int i;
+	for (i = 0; i < nbands; ++i)
+	{
+	    // The length of each array varies with its distance from the
+	    // diagonal
+	    int len = N - (abs(bot + i));
+	    bands[i].clear ();
+	    bands[i].resize (len);
+	}
+	return true;
+    }
+
+    BandedMatrix<T> & operator= (const BandedMatrix<T> &b) 
+    {
+	return Copy (*this, b);
+    }
+
+    BandedMatrix<T> & operator= (const T &e)
+    {
+	int i;
+	for (i = 0; i < nbands; ++i)
+	{
+	    std::fill_n (bands[i].begin(), bands[i].size(), e);
+	}
+	out_of_bounds = e;
+	return (*this);
+    }
+
+    ~BandedMatrix ()
+    {
+	if (bands)
+	    delete[] bands;
+    }
+
+private:
+    // Return false if coordinates are out of bounds
+    inline bool check_bounds (int i, int j, int &v, int &e) const
+    {
+	v = (j - i) - bot;
+	e = (i >= j) ? j : i;
+	return !(v < 0 || v >= nbands || 
+		 e < 0 || (unsigned int)e >= bands[v].size());
+    }
+
+    static BandedMatrix & Copy (BandedMatrix &a, const BandedMatrix &b)
+    {
+	if (a.bands) delete[] a.bands;
+	a.top = b.top;
+	a.bot = b.bot;
+	a.N = b.N;
+	a.out_of_bounds = b.out_of_bounds;
+	a.nbands = a.top - a.bot + 1;
+	a.bands = new std::vector<T>[a.nbands];
+	int i;
+	for (i = 0; i < a.nbands; ++i)
+	{
+	    a.bands[i] = b.bands[i];
+	}
+	return a;
+    }
+
+public:
+    T &element (int i, int j)
+    {
+	int v, e;
+	if (check_bounds(i, j, v, e))
+	    return (bands[v][e]);
+	else
+	    return out_of_bounds;
+    }
+
+    const T &element (int i, int j) const
+    {
+	int v, e;
+	if (check_bounds(i, j, v, e))
+	    return (bands[v][e]);
+	else
+	    return out_of_bounds;
+    }
+
+    inline T & operator() (int i, int j) 
+    {
+	return element (i-1,j-1);
+    }
+
+    inline const T & operator() (int i, int j) const
+    {
+	return element (i-1,j-1);
+    }
+
+    size_type num_rows() const { return N; }
+
+    size_type num_cols() const { return N; }
+
+    const BandedMatrixRow<T> operator[] (int row) const
+    {
+	return BandedMatrixRow<T>(*this, row);
+    }
+
+    BandedMatrixRow<T> operator[] (int row)
+    {
+	return BandedMatrixRow<T>(*this, row);
+    }
+
+
+private:
+
+    int top;
+    int bot;
+    int nbands;
+    std::vector<T> *bands;
+    int N;
+    T out_of_bounds;
+
+};
+
+
+template <class T>
+std::ostream &operator<< (std::ostream &out, const BandedMatrix<T> &m)
+{
+    unsigned int i, j;
+    for (i = 0; i < m.num_rows(); ++i)
+    {
+	for (j = 0; j < m.num_cols(); ++j)
+	{
+	    out << m.element (i, j) << " ";
+	}
+	out << std::endl;
+    }
+    return out;
+}
+
+
+
+/*
+ * Helper class for the intermediate in the [][] operation.
+ */
+template <class T> class BandedMatrixRow
+{
+public:
+    BandedMatrixRow (const BandedMatrix<T> &_m, int _row) : bm(_m), i(_row)
+    { }
+
+    BandedMatrixRow (BandedMatrix<T> &_m, int _row) : bm(_m), i(_row)
+    { }
+
+    ~BandedMatrixRow () {}
+
+    typename BandedMatrix<T>::element_type & operator[] (int j)
+    {
+	return const_cast<BandedMatrix<T> &>(bm).element (i, j);
+    }
+
+    const typename BandedMatrix<T>::element_type & operator[] (int j) const
+    {
+	return bm.element (i, j);
+    }
+
+private:
+    const BandedMatrix<T> &bm;
+    int i;
+};
+
+
+/*
+ * Vector multiplication
+ */
+
+template <class Vector, class Matrix>
+Vector operator* (const Matrix &m, const Vector &v)
+{
+    typename Matrix::size_type M = m.num_rows();
+    typename Matrix::size_type N = m.num_cols();
+
+    assert (N <= v.size());
+    //if (M > v.size())
+    //	return Vector();
+
+    Vector r(N);
+    for (unsigned int i = 0; i < M; ++i)
+    {
+	typename Matrix::element_type sum = 0;
+	for (unsigned int j = 0; j < N; ++j)
+	{
+	    sum += m[i][j] * v[j];
+	}
+	r[i] = sum;
+    }
+    return r;
+}
+
+
+
+/*
+ * LU factor a diagonally banded matrix using Crout's algorithm, but
+ * limiting the trailing sub-matrix multiplication to the non-zero
+ * elements in the diagonal bands.  Return nonzero if a problem occurs.
+ */
+template <class MT>
+int LU_factor_banded (MT &A, unsigned int bands)
+{
+    typename MT::size_type M = A.num_rows();
+    typename MT::size_type N = A.num_cols();
+    if (M != N)
+	return 1;
+
+    typename MT::size_type i,j,k;
+    typename MT::element_type sum;
+
+    for (j = 1; j <= N; ++j)
+    {
+	// Check for zero pivot
+        if ( A(j,j) == 0 )                 
+            return 1;
+
+	// Calculate rows above and on diagonal. A(1,j) remains as A(1,j).
+	for (i = (j > bands) ? j-bands : 1; i <= j; ++i)
+	{	
+	    sum = 0;
+	    for (k = (j > bands) ? j-bands : 1; k < i; ++k)
+	    {
+		sum += A(i,k)*A(k,j);
+	    }
+	    A(i,j) -= sum;
+	}
+
+	// Calculate rows below the diagonal.
+	for (i = j+1; (i <= M) && (i <= j+bands); ++i)
+	{
+	    sum = 0;
+	    for (k = (i > bands) ? i-bands : 1; k < j; ++k)
+	    {
+		sum += A(i,k)*A(k,j);
+	    }
+	    A(i,j) = (A(i,j) - sum) / A(j,j);
+	}
+    }
+
+    return 0;
+}   
+
+
+
+/*
+ * Solving (LU)x = B.  First forward substitute to solve for y, Ly = B.
+ * Then backwards substitute to find x, Ux = y.  Return nonzero if a
+ * problem occurs.  Limit the substitution sums to the elements on the
+ * bands above and below the diagonal.
+ */
+template <class MT, class Vector>
+int LU_solve_banded(const MT &A, Vector &b, unsigned int bands)
+{
+    typename MT::size_type i,j;
+    typename MT::size_type M = A.num_rows();
+    typename MT::size_type N = A.num_cols();
+    typename MT::element_type sum;
+
+    if (M != N || M == 0)
+	return 1;
+
+    // Forward substitution to find y.  The diagonals of the lower
+    // triangular matrix are taken to be 1.
+    for (i = 2; i <= M; ++i)
+    {
+	sum = b[i-1];
+	for (j = (i > bands) ? i-bands : 1; j < i; ++j)
+	{
+	    sum -= A(i,j)*b[j-1];
+	}
+	b[i-1] = sum;
+    }
+
+    // Now for the backward substitution
+    b[M-1] /= A(M,M);
+    for (i = M-1; i >= 1; --i)
+    {
+	if (A(i,i) == 0)	// oops!
+	    return 1;
+	sum = b[i-1];
+	for (j = i+1; (j <= N) && (j <= i+bands); ++j)
+	{
+	    sum -= A(i,j)*b[j-1];
+	}
+	b[i-1] = sum / A(i,i);
+    }
+
+    return 0;
+}
+
+
+#endif /* _BANDEDMATRIX_ID */
+

xs/src/BSpline/COPYRIGHT

@@ -0,0 +1,44 @@
+Copyright (c) 1998-2009,2015
+University Corporation for Atmospheric Research, UCAR
+All rights reserved.
+
+This software is licensed with the standard BSD license:
+
+ http://www.opensource.org/licenses/bsd-license.html
+
+When citing this software, here is a suggested reference:
+
+ This software is written by Gary Granger of the National Center for
+ Atmospheric Research (NCAR), sponsored by the National Science Foundation
+ (NSF).  The algorithm is based on the cubic spline described by Katsuyuki
+ Ooyama in Montly Weather Review, Vol 115, October 1987.  This
+ implementation has benefited from comparisons with a previous FORTRAN
+ implementation by James L. Franklin, NOAA/Hurricane Research Division.
+
+The text of the license is reproduced below:
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+
+    * Neither the name of the UCAR nor the names of its contributors may be
+      used to endorse or promote products derived from this software
+      without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.

xs/src/libslic3r/LayerHeightSpline.cpp

@@ -0,0 +1,203 @@
+/*
+ * This class represents a set of layers and their heights.
+ * It is intended for smoothing the height distribution (avoid very thin
+ * layers next to thick layers) and to correctly interpolate higher layers if
+ * a layer height changes somewhere in a lower position at the object.
+ * Uses http://www.eol.ucar.edu/homes/granger/bspline/doc/ for spline computation.
+ */
+
+#include "LayerHeightSpline.hpp"
+#include <cmath> // std::abs
+
+namespace Slic3r {
+
+LayerHeightSpline::LayerHeightSpline()
+:   _object_height(0)
+{
+    this->_is_valid = false;
+    this->_layers_updated = false;
+    this->_layer_heights_updated = false;
+}
+
+LayerHeightSpline::LayerHeightSpline(const LayerHeightSpline &other)
+{
+    *this = other;
+}
+
+LayerHeightSpline& LayerHeightSpline::operator=(const LayerHeightSpline &other)
+{
+    this->_object_height = other._object_height;
+    this->_layers = other._layers;
+    this->_layer_heights = other._layer_heights;
+    this->_is_valid = other._is_valid;
+    this->_layers_updated = other._layers_updated;
+    this->_layer_heights_updated = other._layer_heights_updated;
+    if(this->_is_valid) {
+        this->_updateBSpline();
+    }
+    return *this;
+}
+
+/*
+ * Indicates whether the object has valid data and the spline was successfully computed or not.
+ */
+bool LayerHeightSpline::hasData()
+{
+    return this->_is_valid;
+}
+
+/*
+ * Set absolute layer positions in object coordinates.
+ * Heights (thickness of each layer) is generated from this list.
+ */
+bool LayerHeightSpline::setLayers(std::vector<coordf_t> layers)
+{
+    this->_layers = layers;
+
+    // generate updated layer height list from layers
+    this->_layer_heights.clear();
+    coordf_t last_z = 0;
+    for (std::vector<coordf_t>::const_iterator l = this->_layers.begin(); l != this->_layers.end(); ++l) {
+        this->_layer_heights.push_back(*l-last_z);
+        last_z = *l;
+    }
+
+    this->_layers_updated = true;
+    this->_layer_heights_updated = false;
+
+    return this->_updateBSpline();
+}
+
+
+/*
+ * Update only the desired thickness of the layers, but not their positions!
+ * This modifies the y-values for the spline computation and only affects
+ * the resulting layers which can be obtained with getInterpolatedLayers.
+ * The argument vector must be of the same size as the layers vector.
+ */
+bool LayerHeightSpline::updateLayerHeights(std::vector<coordf_t> heights)
+{
+    bool result = false;
+
+    // do we receive the correct number of values?
+    if(heights.size() == this->_layers.size()) {
+        this->_layer_heights = heights;
+        result = this->_updateBSpline();
+    }else{
+        std::cerr << "Unable to update layer heights. You provided " << heights.size() << " layers, but " << this->_layers.size()-1 << " expected" << std::endl;
+    }
+
+    this->_layers_updated = false;
+    this->_layer_heights_updated = true;
+
+    return result;
+}
+
+/*
+ * Reset the this object, remove database and interpolated results.
+ */
+void LayerHeightSpline::clear()
+{
+    this->_layers.clear();
+    this->_layer_heights.clear();
+    this->_layer_height_spline.reset();
+    this->_is_valid = false;
+    this->_layers_updated = false;
+    this->_layer_heights_updated = false;
+}
+
+
+/*
+ * Get a full set of layer z-positions by interpolation along the spline.
+ */
+std::vector<coordf_t> LayerHeightSpline::getInterpolatedLayers() const
+{
+    std::vector<coordf_t> layers;
+    if(this->_is_valid) {
+        // preserve first layer for bed contact
+        layers.push_back(this->_layers[0]);
+        coordf_t z = this->_layers[0];
+        coordf_t h;
+        coordf_t h_diff = 0;
+        coordf_t last_h_diff = 0;
+        coordf_t eps = 0.0001;
+        while(z <= this->_object_height) {
+            h = 0;
+            h_diff = 0;
+            // find intersection between layer height and spline
+            do {
+                last_h_diff = h_diff;
+                h += h_diff/2;
+                h = this->_layer_height_spline->evaluate(z+h);
+                h_diff = this->_layer_height_spline->evaluate(z+h) - h;
+            } while(std::abs(h_diff) > eps && std::abs(h_diff - last_h_diff) > eps);
+
+            if(z+h > this->_object_height) {
+                z += this->_layer_height_spline->evaluate(layers.back()); // re-use last layer height if outside of defined range
+            }else{
+                z += h;
+            }
+            layers.push_back(z);
+        }
+        // how to make sure, the last layer is not higher than object while maintaining between min/max layer height?
+    }
+    return layers;
+}
+
+/*
+ * Evaluate interpolated layer height (thickness) at given z-position
+ */
+const coordf_t LayerHeightSpline::getLayerHeightAt(coordf_t height)
+{
+    coordf_t result = 0;
+    if (this->_is_valid) {
+        if(height <= this->_layers[0]) {
+            result = this->_layers[0]; // return first_layer height
+        }else if (height > this->_layers.back()){
+            result = this->_layer_height_spline->evaluate(this->_layers.back()); // repeat last value for height > last layer
+        }else{
+            result = this->_layer_height_spline->evaluate(height); // return interpolated layer height
+        }
+    }
+    return result;
+}
+
+/*
+ * Internal method to re-compute the spline
+ */
+bool LayerHeightSpline::_updateBSpline()
+{
+    bool result = false;
+    //TODO: exception if not enough points?
+
+    // copy layer vectors and duplicate a datapoint at the front / end to achieve correct boundary conditions
+    this->_spline_layers = this->_layers;
+    this->_spline_layers[0] = 0;
+    this->_spline_layers.push_back(this->_spline_layers.back()+1);
+
+    this->_spline_layer_heights = this->_layer_heights;
+    this->_spline_layer_heights[0] = this->_spline_layer_heights[1]; // override fixed first layer height with first "real" layer
+    this->_spline_layer_heights.push_back(this->_spline_layer_heights.back());
+
+    this->_layer_height_spline.reset(new BSpline<double>(&this->_spline_layers[0],
+            this->_spline_layers.size(),
+            &this->_spline_layer_heights[0],
+            0,
+            1,
+            0)
+    );
+
+    if (this->_layer_height_spline->ok()) {
+        result = true;
+    } else {
+        result = false;
+        std::cerr << "Spline setup failed." << std::endl;
+    }
+
+    this->_is_valid = result;
+
+    return result;
+}
+
+
+}

xs/src/libslic3r/LayerHeightSpline.hpp

@@ -0,0 +1,43 @@
+#ifndef slic3r_LayerHeightSpline_hpp_
+#define slic3r_LayerHeightSpline_hpp_
+
+#include "libslic3r.h"
+#include "BSpline/BSpline.h" // Warning: original BSplineBase.h/cpp merged into BSpline.h/cpp to avoid dependency issues caused by Build::WithXSpp which tries to compile all .cpp files in /src
+
+namespace Slic3r {
+
+
+class LayerHeightSpline
+{
+    public:
+    LayerHeightSpline();
+    LayerHeightSpline(const LayerHeightSpline &other);
+    LayerHeightSpline& operator=(const LayerHeightSpline &other);
+    void setObjectHeight(coordf_t object_height) { this->_object_height = object_height; };
+    bool hasData(); // indicate that we have valid data;
+    bool setLayers(std::vector<coordf_t> layers);
+    bool updateLayerHeights(std::vector<coordf_t> heights);
+    bool layersUpdated() const { return this->_layers_updated; }; // true if the basis set of layers was updated (by the slicing algorithm)
+    bool layerHeightsUpdated() const { return this->_layer_heights_updated; }; // true if the heights where updated (by the spline control user interface)
+    void clear();
+    std::vector<coordf_t> getOriginalLayers() const { return this->_layers; };
+    std::vector<coordf_t> getInterpolatedLayers() const;
+    const coordf_t getLayerHeightAt(coordf_t height);
+
+    private:
+    bool _updateBSpline();
+
+    coordf_t _object_height;
+    bool _is_valid;
+    bool _layers_updated;
+    bool _layer_heights_updated;
+    std::vector<coordf_t> _layers;
+    std::vector<coordf_t> _layer_heights;
+    std::vector<coordf_t> _spline_layers;
+    std::vector<coordf_t> _spline_layer_heights;
+    std::unique_ptr<BSpline<double>> _layer_height_spline;
+};
+
+}
+
+#endif

xs/src/libslic3r/Model.cpp

@@ -430,6 +430,7 @@ ModelObject::ModelObject(Model *model, const ModelObject &other, bool copy_volum
     config(other.config),
     layer_height_ranges(other.layer_height_ranges),
     part_number(other.part_number),
+    layer_height_spline(other.layer_height_spline),
     origin_translation(other.origin_translation),
     _bounding_box(other._bounding_box),
     _bounding_box_valid(other._bounding_box_valid),
@@ -460,6 +461,7 @@ ModelObject::swap(ModelObject &other)
     std::swap(this->volumes,                other.volumes);
     std::swap(this->config,                 other.config);
     std::swap(this->layer_height_ranges,    other.layer_height_ranges);
+    std::swap(this->layer_height_spline,    other.layer_height_spline);
     std::swap(this->origin_translation,     other.origin_translation);
     std::swap(this->_bounding_box,          other._bounding_box);
     std::swap(this->_bounding_box_valid,    other._bounding_box_valid);
@@ -511,7 +513,6 @@ ModelObject::add_instance()
 {
     ModelInstance* i = new ModelInstance(this);
     this->instances.push_back(i);
-    this->invalidate_bounding_box();
     return i;
 }
 
@@ -520,7 +521,6 @@ ModelObject::add_instance(const ModelInstance &other)
 {
     ModelInstance* i = new ModelInstance(this, other);
     this->instances.push_back(i);
-    this->invalidate_bounding_box();
     return i;
 }
 
@@ -530,7 +530,6 @@ ModelObject::delete_instance(size_t idx)
     ModelInstancePtrs::iterator i = this->instances.begin() + idx;
     delete *i;
     this->instances.erase(i);
-    this->invalidate_bounding_box();
 }
 
 void
@@ -544,6 +543,7 @@ ModelObject::clear_instances()
 {
     while (!this->instances.empty())
         this->delete_last_instance();
+
 }
 
 // this returns the bounding box of the *transformed* instances

xs/src/libslic3r/Model.hpp

@@ -7,6 +7,7 @@
 #include "Layer.hpp"
 #include "Point.hpp"
 #include "TriangleMesh.hpp"
+#include "LayerHeightSpline.hpp"
 #include <map>
 #include <string>
 #include <utility>
@@ -267,6 +268,7 @@ class ModelObject
     t_layer_height_ranges layer_height_ranges; ///< Variation of a layer thickness for spans of Z coordinates.
 
     int part_number; ///< It's used for the 3MF items part numbers in the build element.
+    LayerHeightSpline layer_height_spline;     ///< Spline based variations of layer thickness for interactive user manipulation
 
     Pointf3 origin_translation;
     ///< This vector accumulates the total translation applied to the object by the

xs/src/libslic3r/Print.cpp

@@ -172,8 +172,9 @@ Print::invalidate_state_by_config(const PrintConfigBase &config)
             || opt_key == "brim_connections_width") {
             steps.insert(psBrim);
             steps.insert(psSkirt);
-        } else if (opt_key == "nozzle_diameter"
+        } else if (opt_key == "nozzle_diameter") {
-            || opt_key == "resolution"
+                osteps.insert(posLayers);
+        } else if (opt_key == "resolution"
             || opt_key == "z_steps_per_mm") {
             osteps.insert(posSlice);
         } else if (opt_key == "avoid_crossing_perimeters"

xs/src/libslic3r/Print.hpp

@@ -13,7 +13,8 @@
 #include "Layer.hpp"
 #include "Model.hpp"
 #include "PlaceholderParser.hpp"
-
+#include "SlicingAdaptive.hpp"
+#include "LayerHeightSpline.hpp"
 
 namespace Slic3r {
 
@@ -26,7 +27,7 @@ enum PrintStep {
     psSkirt, psBrim,
 };
 enum PrintObjectStep {
-    posSlice, posPerimeters, posDetectSurfaces,
+    posLayers, posSlice, posPerimeters, posDetectSurfaces,
     posPrepareInfill, posInfill, posSupportMaterial,
 };
 
@@ -81,6 +82,8 @@ class PrintObject
     PrintObjectConfig config;
     t_layer_height_ranges layer_height_ranges;
     
+    LayerHeightSpline layer_height_spline;
+
     // this is set to true when LayerRegion->slices is split in top/internal/bottom
     // so that next call to make_perimeters() performs a union() before computing loops
     bool typed_slices;

xs/src/libslic3r/PrintConfig.cpp

@@ -23,6 +23,25 @@ PrintConfigDef::PrintConfigDef()
     
     ConfigOptionDef* def;
 
+    def = this->add("adaptive_slicing", coBool);
+    def->label = "Use adaptive slicing";
+    def->category = "Layers and Perimeters";
+    def->tooltip = "Automatically determine layer heights by the objects topology instead of using the static value.";
+    def->cli = "adaptive-slicing!";
+    def->default_value = new ConfigOptionBool(false);
+
+    def = this->add("adaptive_slicing_quality", coPercent);
+    def->label = "Adaptive quality";
+    def->category = "Layers and Perimeters";
+    def->tooltip = "Controls the quality / printing time tradeoff for adaptive layer generation. 0 -> fastest printing with max layer height, 100 -> highest quality, min layer height";
+    def->sidetext = "%";
+    def->cli = "adaptive_slicing_quality=f";
+    def->min = 0;
+    def->max = 100;
+    def->gui_type = "slider";
+    def->width = 200;
+    def->default_value = new ConfigOptionPercent(75);
+
     def = this->add("avoid_crossing_perimeters", coBool);
     def->label = "Avoid crossing perimeters";
     def->category = "Layers and Perimeters";
@@ -752,6 +771,12 @@ PrintConfigDef::PrintConfigDef()
     def->min = 0;
     def->default_value = new ConfigOptionFloat(0.3);
 
+    def = this->add("match_horizontal_surfaces", coBool);
+    def->label = "Match horizontal surfaces";
+    def->tooltip = "Try to match horizontal surfaces during the slicing process. Matching is not guaranteed, very small surfaces and multiple surfaces with low vertical distance might cause bad results.";
+    def->cli = "match-horizontal-surfaces!";
+    def->default_value = new ConfigOptionBool(false);
+
     def = this->add("max_fan_speed", coInt);
     def->label = "Max";
     def->tooltip = "This setting represents the maximum speed of your fan.";
@@ -761,6 +786,18 @@ PrintConfigDef::PrintConfigDef()
     def->max = 100;
     def->default_value = new ConfigOptionInt(100);
 
+    def = this->add("max_layer_height", coFloats);
+	def->label = "Max";
+	def->tooltip = "This is the highest printable layer height for this extruder and limits the resolution for adaptive slicing. Typical values are slightly smaller than nozzle_diameter.";
+	def->sidetext = "mm";
+	def->cli = "max-layer-height=f@";
+	def->min = 0;
+	{
+		ConfigOptionFloats* opt = new ConfigOptionFloats();
+		opt->values.push_back(0.3);
+		def->default_value = opt;
+	}
+
     def = this->add("max_print_speed", coFloat);
     def->label = "Max print speed";
     def->category = "Speed";
@@ -788,6 +825,18 @@ PrintConfigDef::PrintConfigDef()
     def->max = 100;
     def->default_value = new ConfigOptionInt(35);
 
+    def = this->add("min_layer_height", coFloats);
+	def->label = "Min";
+	def->tooltip = "This is the lowest printable layer height for this extruder and limits the resolution for adaptive slicing. Typical values are 0.1 or 0.05.";
+	def->sidetext = "mm";
+	def->cli = "min-layer-height=f@";
+	def->min = 0;
+	{
+		ConfigOptionFloats* opt = new ConfigOptionFloats();
+		opt->values.push_back(0.15);
+		def->default_value = opt;
+	}
+
     def = this->add("min_print_speed", coFloat);
     def->label = "Min print speed";
     def->tooltip = "Slic3r will not scale speed down below this speed.";

xs/src/libslic3r/PrintConfig.hpp

@@ -154,12 +154,15 @@ class StaticPrintConfig : public PrintConfigBase, public StaticConfig
 class PrintObjectConfig : public virtual StaticPrintConfig
 {
     public:
+	ConfigOptionBool                adaptive_slicing;
+	ConfigOptionPercent             adaptive_slicing_quality;
     ConfigOptionBool                dont_support_bridges;
     ConfigOptionFloatOrPercent      extrusion_width;
     ConfigOptionFloatOrPercent      first_layer_height;
     ConfigOptionBool                infill_only_where_needed;
     ConfigOptionBool                interface_shells;
     ConfigOptionFloat               layer_height;
+    ConfigOptionBool                match_horizontal_surfaces;
     ConfigOptionInt                 raft_layers;
     ConfigOptionFloat               regions_overlap;
     ConfigOptionEnum<SeamPosition>  seam_position;
@@ -186,12 +189,15 @@ class PrintObjectConfig : public virtual StaticPrintConfig
     }
     
     virtual ConfigOption* optptr(const t_config_option_key &opt_key, bool create = false) {
+        OPT_PTR(adaptive_slicing);
+        OPT_PTR(adaptive_slicing_quality);
         OPT_PTR(dont_support_bridges);
         OPT_PTR(extrusion_width);
         OPT_PTR(first_layer_height);
         OPT_PTR(infill_only_where_needed);
         OPT_PTR(interface_shells);
         OPT_PTR(layer_height);
+        OPT_PTR(match_horizontal_surfaces);
         OPT_PTR(raft_layers);
         OPT_PTR(regions_overlap);
         OPT_PTR(seam_position);
@@ -428,7 +434,9 @@ class PrintConfig : public GCodeConfig
     ConfigOptionBool                infill_first;
     ConfigOptionFloat               interior_brim_width;
     ConfigOptionInt                 max_fan_speed;
+    ConfigOptionFloats              max_layer_height;
     ConfigOptionInt                 min_fan_speed;
+    ConfigOptionFloats              min_layer_height;
     ConfigOptionFloat               min_print_speed;
     ConfigOptionFloat               min_skirt_length;
     ConfigOptionFloats              nozzle_diameter;
@@ -488,7 +496,9 @@ class PrintConfig : public GCodeConfig
         OPT_PTR(infill_first);
         OPT_PTR(interior_brim_width);
         OPT_PTR(max_fan_speed);
+        OPT_PTR(max_layer_height);
         OPT_PTR(min_fan_speed);
+        OPT_PTR(min_layer_height);
         OPT_PTR(min_print_speed);
         OPT_PTR(min_skirt_length);
         OPT_PTR(nozzle_diameter);

xs/src/libslic3r/PrintObject.cpp

@@ -10,7 +10,8 @@ namespace Slic3r {
 PrintObject::PrintObject(Print* print, ModelObject* model_object, const BoundingBoxf3 &modobj_bbox)
 :   typed_slices(false),
     _print(print),
-    _model_object(model_object)
+    _model_object(model_object),
+    layer_height_spline(model_object->layer_height_spline)
 {
     // Compute the translation to be applied to our meshes so that we work with smaller coordinates
     {
@@ -223,9 +224,13 @@ PrintObject::invalidate_state_by_config(const PrintConfigBase &config)
     for (const t_config_option_key &opt_key : diff) {
         if (opt_key == "layer_height"
             || opt_key == "first_layer_height"
-            || opt_key == "xy_size_compensation"
+            || opt_key == "adaptive_slicing"
-            || opt_key == "raft_layers"
+            || opt_key == "adaptive_slicing_quality"
+            || opt_key == "match_horizontal_surfaces"
             || opt_key == "regions_overlap") {
+            steps.insert(posLayers);
+        } else if (opt_key == "xy_size_compensation"
+            || opt_key == "raft_layers") {
             steps.insert(posSlice);
         } else if (opt_key == "support_material_contact_distance") {
             steps.insert(posSlice);
@@ -296,6 +301,8 @@ PrintObject::invalidate_step(PrintObjectStep step)
         this->invalidate_step(posPerimeters);
         this->invalidate_step(posDetectSurfaces);
         this->invalidate_step(posSupportMaterial);
+    }else if (step == posLayers) {
+        this->invalidate_step(posSlice);
     } else if (step == posSupportMaterial) {
         this->_print->invalidate_step(psSkirt);
         this->_print->invalidate_step(psBrim);
@@ -538,11 +545,11 @@ PrintObject::bridge_over_infill()
 // adjust the layer height to the next multiple of the z full-step resolution
 coordf_t PrintObject::adjust_layer_height(coordf_t layer_height) const
 {
-	coordf_t result = layer_height;
+    coordf_t result = layer_height;
-	if(this->_print->config.z_steps_per_mm > 0) {
+    if(this->_print->config.z_steps_per_mm > 0) {
-		coordf_t min_dz = 1 / this->_print->config.z_steps_per_mm * 4;
+        coordf_t min_dz = 1 / this->_print->config.z_steps_per_mm * 4;
-		result = int(layer_height / min_dz + 0.5) * min_dz;
+        result = int(layer_height / min_dz + 0.5) * min_dz;
-	}
+    }
 
     return result > 0 ? result : layer_height;
 }
@@ -553,64 +560,154 @@ std::vector<coordf_t> PrintObject::generate_object_layers(coordf_t first_layer_h
 
     std::vector<coordf_t> result;
 
-	coordf_t min_nozzle_diameter = 1.0;
+    // collect values from config
+    coordf_t min_nozzle_diameter = 1.0;
+    coordf_t min_layer_height = 0.0;
+    coordf_t max_layer_height = 10.0;
     std::set<size_t> object_extruders = this->_print->object_extruders();
-	for (std::set<size_t>::const_iterator it_extruder = object_extruders.begin(); it_extruder != object_extruders.end(); ++ it_extruder) {
+    for (std::set<size_t>::const_iterator it_extruder = object_extruders.begin(); it_extruder != object_extruders.end(); ++ it_extruder) {
-		min_nozzle_diameter = std::min(min_nozzle_diameter, this->_print->config.nozzle_diameter.get_at(*it_extruder));
+        min_nozzle_diameter = std::min(min_nozzle_diameter, this->_print->config.nozzle_diameter.get_at(*it_extruder));
-	}
+        min_layer_height = std::max(min_layer_height, this->_print->config.min_layer_height.get_at(*it_extruder));
-	coordf_t layer_height = std::min(min_nozzle_diameter, this->config.layer_height.getFloat());
+        max_layer_height = std::min(max_layer_height, this->_print->config.max_layer_height.get_at(*it_extruder));
-	layer_height = this->adjust_layer_height(layer_height);
+
+    }
+    coordf_t layer_height = std::min(min_nozzle_diameter, this->config.layer_height.getFloat());
+    layer_height = this->adjust_layer_height(layer_height);
     this->config.layer_height.value = layer_height;
+
+    // respect first layer height
     if(first_layer_height) {
         result.push_back(first_layer_height);
     }
 
     coordf_t print_z = first_layer_height;
     coordf_t height = first_layer_height;
-    // loop until we have at least one layer and the max slice_z reaches the object height
-    while (print_z < unscale(this->size.z)) {
-        height = layer_height;
 
-        // look for an applicable custom range
+    // Update object size at the spline object to define upper border
-        for (t_layer_height_ranges::const_iterator it_range = this->layer_height_ranges.begin(); it_range != this->layer_height_ranges.end(); ++ it_range) {
+    this->layer_height_spline.setObjectHeight(unscale(this->size.z));
-            if(print_z >= it_range->first.first && print_z <= it_range->first.second) {
+    if (this->state.is_done(posLayers)) {
-                if(it_range->second > 0) {
+        // layer heights are already generated, just update layers from spline
-                    height = it_range->second;
+        // we don't need to respect first layer here, it's correctly provided by the spline object
+        result = this->layer_height_spline.getInterpolatedLayers();
+    }else{ // create new set of layers
+        // create stateful objects and variables for the adaptive slicing process
+        SlicingAdaptive as;
+        coordf_t adaptive_quality = this->config.adaptive_slicing_quality.value;
+        if(this->config.adaptive_slicing.value) {
+            const ModelVolumePtrs volumes = this->model_object()->volumes;
+            for (ModelVolumePtrs::const_iterator it = volumes.begin(); it != volumes.end(); ++ it)
+                if (! (*it)->modifier)
+                    as.add_mesh(&(*it)->mesh);
+            as.prepare(unscale(this->size.z));
+        }
+
+        // loop until we have at least one layer and the max slice_z reaches the object height
+        while (print_z < unscale(this->size.z)) {
+
+            if (this->config.adaptive_slicing.value) {
+                height = 999;
+
+                // determine next layer height
+                height = as.next_layer_height(print_z, adaptive_quality, min_layer_height, max_layer_height);
+
+                // check for horizontal features and object size
+                if(this->config.match_horizontal_surfaces.value) {
+                    coordf_t horizontal_dist = as.horizontal_facet_distance(print_z + height, min_layer_height);
+                    if((horizontal_dist < min_layer_height) && (horizontal_dist > 0)) {
+                        #ifdef SLIC3R_DEBUG
+                        std::cout << "Horizontal feature ahead, distance: " << horizontal_dist << std::endl;
+                        #endif
+                        // can we shrink the current layer a bit?
+                        if(height-(min_layer_height - horizontal_dist) > min_layer_height) {
+                            // yes we can
+                            height -= (min_layer_height - horizontal_dist);
+                            #ifdef SLIC3R_DEBUG
+                            std::cout << "Shrink layer height to " << height << std::endl;
+                            #endif
+                        }else{
+                            // no, current layer would become too thin
+                            height += horizontal_dist;
+                            #ifdef SLIC3R_DEBUG
+                            std::cout << "Widen layer height to " << height << std::endl;
+                            #endif
+                        }
+                    }
+                }
+            }else{
+                height = layer_height;
+            }
+
+            // look for an applicable custom range
+            for (t_layer_height_ranges::const_iterator it_range = this->layer_height_ranges.begin(); it_range != this->layer_height_ranges.end(); ++ it_range) {
+                if(print_z >= it_range->first.first && print_z <= it_range->first.second) {
+                    if(it_range->second > 0) {
+                        height = it_range->second;
+                    }
+                }
+            }
+
+            print_z += height;
+
+            result.push_back(print_z);
+        }
+
+        // Reduce or thicken the top layer in order to match the original object size.
+        // This is not actually related to z_steps_per_mm but we only enable it in case
+        // user provided that value, as it means they really care about the layer height
+        // accuracy and we don't provide unexpected result for people noticing the last
+        // layer has a different layer height.
+        if (this->_print->config.z_steps_per_mm > 0 && result.size() > 1 && !this->config.adaptive_slicing.value) {
+            coordf_t diff = result.back() - unscale(this->size.z);
+            int last_layer = result.size()-1;
+
+            if (diff < 0) {
+                // we need to thicken last layer
+                coordf_t new_h = result[last_layer] - result[last_layer-1];
+                new_h = std::min(min_nozzle_diameter, new_h - diff); // add (negativ) diff value
+                result[last_layer] = result[last_layer-1] + new_h;
+            } else {
+                // we need to reduce last layer
+                coordf_t new_h = result[last_layer] - result[last_layer-1];
+                if(min_nozzle_diameter/2 < new_h) { //prevent generation of a too small layer
+                    new_h = std::max(min_nozzle_diameter/2, new_h - diff); // subtract (positive) diff value
+                    result[last_layer] = result[last_layer-1] + new_h;
                 }
             }
         }
 
-        print_z += height;
+        // Store layer vector for interactive manipulation
+        this->layer_height_spline.setLayers(result);
+        if (this->config.adaptive_slicing.value) { // smoothing after adaptive algorithm
+            result = this->layer_height_spline.getInterpolatedLayers();
+        }
 
-        result.push_back(print_z);
+        this->state.set_done(posLayers);
     }
 
-    // Reduce or thicken the top layer in order to match the original object size.
+    // push modified spline object back to model
-    // This is not actually related to z_steps_per_mm but we only enable it in case
+    this->_model_object->layer_height_spline = this->layer_height_spline;
-    // user provided that value, as it means they really care about the layer height
-    // accuracy and we don't provide unexpected result for people noticing the last
-    // layer has a different layer height.
-    if (this->_print->config.z_steps_per_mm > 0 && result.size() > 1) {
-        coordf_t diff = result.back() - unscale(this->size.z);
-        int last_layer = result.size()-1;
 
-        if (diff < 0) {
+    // apply z-gradation (this is redundant for static layer height...)
-            // we need to thicken last layer
+    coordf_t gradation = 1 / this->_print->config.z_steps_per_mm * 4;
-            coordf_t new_h = result[last_layer] - result[last_layer-1];
+    if(this->_print->config.z_steps_per_mm > 0) {
-            new_h = std::min(min_nozzle_diameter, new_h - diff); // add (negativ) diff value
+        coordf_t last_z = 0;
-            std::cout << new_h << std::endl;
+        coordf_t height;
-            result[last_layer] = result[last_layer-1] + new_h;
+        for(std::vector<coordf_t>::iterator l = result.begin(); l != result.end(); ++l) {
-        } else {
+            height = *l - last_z;
-            // we need to reduce last layer
+
-            coordf_t new_h = result[last_layer] - result[last_layer-1];
+            coordf_t gradation_effect = unscale((scale_(height)) % (scale_(gradation)));
-            if(min_nozzle_diameter/2 < new_h) { //prevent generation of a too small layer
+            if(gradation_effect > gradation/2 && (height + (gradation-gradation_effect)) <= max_layer_height) { // round up
-                new_h = std::max(min_nozzle_diameter/2, new_h - diff); // subtract (positive) diff value
+                height = height + (gradation-gradation_effect);
-                std::cout << new_h << std::endl;
+            }else{ // round down
-                result[last_layer] = result[last_layer-1] + new_h;
+                height = height - gradation_effect;
             }
+            height = std::min(std::max(height, min_layer_height), max_layer_height);
+            *l = last_z + height;
+            last_z = *l;
         }
     }
-	return result;
+
+    return result;
 }
 
 // 1) Decides Z positions of the layers,
@@ -625,8 +722,8 @@ std::vector<coordf_t> PrintObject::generate_object_layers(coordf_t first_layer_h
 void PrintObject::_slice()
 {
 
-	coordf_t raft_height = 0;
+    coordf_t raft_height = 0;
-	coordf_t first_layer_height = this->config.first_layer_height.get_abs_value(this->config.layer_height.value);
+    coordf_t first_layer_height = this->config.first_layer_height.get_abs_value(this->config.layer_height.value);
 
 
     // take raft layers into account

xs/src/libslic3r/SlicingAdaptive.cpp

@@ -0,0 +1,182 @@
+#include "libslic3r.h"
+#include "TriangleMesh.hpp"
+#include "SlicingAdaptive.hpp"
+
+#ifdef SLIC3R_DEBUG
+    #undef NDEBUG
+    #define DEBUG
+    #define _DEBUG
+#endif
+
+/* This constant essentially describes the volumetric error at the surface which is induced
+ * by stacking "elliptic" extrusion threads.
+ * It is empirically determined by
+ * 1. measuring the surface profile of printed parts to find
+ * the ratio between layer height and profile height and then
+ * 2. computing the geometric difference between the model-surface and the elliptic profile.
+ * [Link to detailed description follows]
+ */
+#define SURFACE_CONST 0.18403
+
+namespace Slic3r
+{
+
+void SlicingAdaptive::clear()
+{
+    m_meshes.clear();
+    m_faces.clear();
+    m_face_normal_z.clear();
+}
+
+std::pair<float, float> face_z_span(const stl_facet *f)
+{
+    return std::pair<float, float>(
+        std::min(std::min(f->vertex[0].z, f->vertex[1].z), f->vertex[2].z),
+        std::max(std::max(f->vertex[0].z, f->vertex[1].z), f->vertex[2].z));
+}
+
+void SlicingAdaptive::prepare(coordf_t object_size)
+{
+    this->object_size = object_size;
+
+    // 1) Collect faces of all meshes.
+    int nfaces_total = 0;
+    for (std::vector<const TriangleMesh*>::const_iterator it_mesh = m_meshes.begin(); it_mesh != m_meshes.end(); ++ it_mesh)
+        nfaces_total += (*it_mesh)->stl.stats.number_of_facets;
+    m_faces.reserve(nfaces_total);
+    for (std::vector<const TriangleMesh*>::const_iterator it_mesh = m_meshes.begin(); it_mesh != m_meshes.end(); ++ it_mesh)
+        for (int i = 0; i < (*it_mesh)->stl.stats.number_of_facets; ++ i)
+            m_faces.push_back((*it_mesh)->stl.facet_start + i);
+
+    // 2) Sort faces lexicographically by their Z span.
+    std::sort(m_faces.begin(), m_faces.end(), [](const stl_facet *f1, const stl_facet *f2) {
+        std::pair<float, float> span1 = face_z_span(f1);
+        std::pair<float, float> span2 = face_z_span(f2);
+        return span1 < span2;
+    });
+
+    // 3) Generate Z components of the facet normals.
+    m_face_normal_z.assign(m_faces.size(), 0.f);
+    for (size_t iface = 0; iface < m_faces.size(); ++ iface)
+        m_face_normal_z[iface] = m_faces[iface]->normal.z;
+
+    // 4) Reset current facet pointer
+    this->current_facet = 0;
+}
+
+float SlicingAdaptive::next_layer_height(coordf_t z, coordf_t quality_factor, coordf_t min_layer_height, coordf_t max_layer_height)
+{
+    float height = max_layer_height;
+    // factor must be between 0-1, 0 is highest quality, 1 highest print speed.
+
+    // factor must be between 0-1, 0 is highest quality, 1 highest print speed.
+    // Invert the slider scale (100% should represent a very high quality for the user)
+    quality_factor = std::max(0.f, std::min<float>(1.f, 1 - quality_factor/100.f));
+
+    float delta_min = SURFACE_CONST * min_layer_height;
+    float delta_max = SURFACE_CONST * max_layer_height + 0.5 * max_layer_height;
+    float scaled_quality_factor = quality_factor * (delta_max - delta_min) + delta_min;
+
+    bool first_hit = false;
+
+    // find all facets intersecting the slice-layer
+    int ordered_id = current_facet;
+    for (; ordered_id < int(m_faces.size()); ++ ordered_id) {
+        std::pair<float, float> zspan = face_z_span(m_faces[ordered_id]);
+        // facet's minimum is higher than slice_z -> end loop
+        if (zspan.first >= z)
+            break;
+        // facet's maximum is higher than slice_z -> store the first event for next layer_height call to begin at this point
+        if (zspan.second > z) {
+            // first event?
+            if (! first_hit) {
+                first_hit = true;
+                current_facet = ordered_id;
+            }
+            // skip touching facets which could otherwise cause small height values
+            if (zspan.second <= z + EPSILON)
+                continue;
+            // compute height for this facet and store minimum of all heights
+            height = std::min(height, this->_layer_height_from_facet(ordered_id, scaled_quality_factor));
+        }
+    }
+
+    // lower height limit due to printer capabilities
+    height = std::max<float>(height, min_layer_height);
+
+    // check for sloped facets inside the determined layer and correct height if necessary
+    if (height > min_layer_height) {
+        for (; ordered_id < int(m_faces.size()); ++ ordered_id) {
+            std::pair<float, float> zspan = face_z_span(m_faces[ordered_id]);
+            // facet's minimum is higher than slice_z + height -> end loop
+            if (zspan.first >= z + height)
+                break;
+
+            // skip touching facets which could otherwise cause small cusp values
+            if (zspan.second <= z + EPSILON)
+                continue;
+
+            // Compute new height for this facet and check against height.
+            float reduced_height = this->_layer_height_from_facet(ordered_id, scaled_quality_factor);
+
+            float z_diff = zspan.first - z;
+
+            if (reduced_height > z_diff) {
+                if (reduced_height < height) {
+#ifdef DEBUG
+                    std::cout << "adaptive layer computation: height is reduced from " << height;
+#endif
+                    height = reduced_height;
+#ifdef DEBUG
+                    std::cout << "to " << height << " due to higher facet" << std::endl;
+#endif
+                }
+            } else {
+#ifdef DEBUG
+                std::cout << "cusp computation, height is reduced from " << height;
+#endif
+                height = z_diff;
+#ifdef DEBUG
+                std::cout << "to " << height << " due to z-diff" << std::endl;
+#endif
+            }
+        }
+        // lower height limit due to printer capabilities again
+        height = std::max(height, float(min_layer_height));
+    }
+
+#ifdef DEBUG
+    std::cout << "adaptive layer computation, layer-bottom at z:" << z << ", quality_factor:" << quality_factor << ", resulting layer height:" << height << std::endl;
+#endif
+    return height;
+}
+
+// Returns the distance to the next horizontal facet in Z-dir 
+// to consider horizontal object features in slice thickness
+float SlicingAdaptive::horizontal_facet_distance(coordf_t z, coordf_t max_layer_height)
+{
+    for (size_t i = 0; i < m_faces.size(); ++ i) {
+        std::pair<float, float> zspan = face_z_span(m_faces[i]);
+        // facet's minimum is higher than max forward distance -> end loop
+        if (zspan.first > z + max_layer_height)
+            break;
+        // min_z == max_z -> horizontal facet
+        if (zspan.first > z && zspan.first == zspan.second)
+            return zspan.first - z;
+    }
+
+    // objects maximum?
+    return (z + max_layer_height > this->object_size) ?
+        std::max<float>(this->object_size - z, 0.f) :
+        max_layer_height;
+}
+
+// for a given facet, compute maximum height within the allowed surface roughness / stairstepping deviation
+float SlicingAdaptive::_layer_height_from_facet(int ordered_id, float scaled_quality_factor)
+{
+    float normal_z = std::abs(m_face_normal_z[ordered_id]);
+    float height = scaled_quality_factor/(SURFACE_CONST + normal_z/2);
+    return height;
+}
+
+}; // namespace Slic3r

xs/src/libslic3r/SlicingAdaptive.hpp

@@ -0,0 +1,38 @@
+#ifndef slic3r_SlicingAdaptive_hpp_
+#define slic3r_SlicingAdaptive_hpp_
+
+#include "admesh/stl.h"
+
+namespace Slic3r
+{
+
+class TriangleMesh;
+
+class SlicingAdaptive
+{
+public:
+    SlicingAdaptive() {};
+    ~SlicingAdaptive() {};
+    void clear();
+    void add_mesh(const TriangleMesh *mesh) { m_meshes.push_back(mesh); }
+    void prepare(coordf_t object_size);
+    float next_layer_height(coordf_t z, coordf_t quality_factor, coordf_t min_layer_height, coordf_t max_layer_height);
+    float horizontal_facet_distance(coordf_t z, coordf_t max_layer_height);
+
+private:
+    float _layer_height_from_facet(int ordered_id, float scaled_quality_factor);
+
+protected:
+    // id of the current facet from last iteration
+    coordf_t                            object_size;
+    int                                 current_facet;
+    std::vector<const TriangleMesh*>	m_meshes;
+    // Collected faces of all meshes, sorted by raising Z of the bottom most face.
+    std::vector<const stl_facet*>		m_faces;
+    // Z component of face normals, normalized.
+    std::vector<float>					m_face_normal_z;
+};
+
+}; // namespace Slic3r
+
+#endif /* slic3r_SlicingAdaptive_hpp_ */

xs/src/perlglue.cpp

@@ -24,6 +24,7 @@ REGISTER_CLASS(GCodeWriter, "GCode::Writer");
 REGISTER_CLASS(Layer, "Layer");
 REGISTER_CLASS(SupportLayer, "Layer::Support");
 REGISTER_CLASS(LayerRegion, "Layer::Region");
+REGISTER_CLASS(LayerHeightSpline, "LayerHeightSpline");
 REGISTER_CLASS(Line, "Line");
 REGISTER_CLASS(Linef3, "Linef3");
 REGISTER_CLASS(PerimeterGenerator, "Layer::PerimeterGenerator");
@@ -56,6 +57,7 @@ REGISTER_CLASS(GCodeConfig, "Config::GCode");
 REGISTER_CLASS(PrintConfig, "Config::Print");
 REGISTER_CLASS(FullPrintConfig, "Config::Full");
 REGISTER_CLASS(SLAPrint, "SLAPrint");
+REGISTER_CLASS(SlicingAdaptive, "SlicingAdaptive");
 REGISTER_CLASS(Surface, "Surface");
 REGISTER_CLASS(SurfaceCollection, "Surface::Collection");
 REGISTER_CLASS(TriangleMesh, "TriangleMesh");

xs/xsp/LayerHeightSpline.xsp

@@ -0,0 +1,25 @@
+%module{Slic3r::XS};
+
+%{
+#include <xsinit.h>
+#include "libslic3r/LayerHeightSpline.hpp"
+%}
+
+%name{Slic3r::LayerHeightSpline} class LayerHeightSpline {
+    LayerHeightSpline();
+    Clone<LayerHeightSpline> clone()
+        %code%{ RETVAL = THIS; %};
+    
+    void setObjectHeight(coordf_t object_height);
+    bool hasData();
+    bool setLayers(std::vector<double> layers)
+        %code%{ RETVAL = THIS->setLayers(layers); %};
+    bool updateLayerHeights(std::vector<double> heights)
+        %code%{ RETVAL = THIS->updateLayerHeights(heights); %};
+    bool layersUpdated();
+    bool layerHeightsUpdated();
+    void clear();
+    std::vector<double> getOriginalLayers();
+    std::vector<double> getInterpolatedLayers();
+    coordf_t getLayerHeightAt(coordf_t height);
+};

xs/xsp/Model.xsp

@@ -210,6 +210,11 @@ ModelMaterial::attributes()
     void set_layer_height_ranges(t_layer_height_ranges ranges)
         %code%{ THIS->layer_height_ranges = ranges; %};
 
+    Ref<LayerHeightSpline> layer_height_spline()
+        %code%{ RETVAL = &THIS->layer_height_spline; %};
+    void set_layer_height_spline(LayerHeightSpline* spline)
+        %code%{ THIS->layer_height_spline = *spline; %};
+
     Ref<Pointf3> origin_translation()
         %code%{ RETVAL = &THIS->origin_translation; %};
     void set_origin_translation(Pointf3* point)

xs/xsp/Print.xsp

@@ -12,6 +12,7 @@
 IV
 _constant()
   ALIAS:
+    STEP_LAYERS             = posLayers
     STEP_SLICE              = posSlice
     STEP_PERIMETERS         = posPerimeters
     STEP_DETECT_SURFACES    = posDetectSurfaces
@@ -59,6 +60,8 @@ _constant()
     Points copies();
     t_layer_height_ranges layer_height_ranges()
         %code%{ RETVAL = THIS->layer_height_ranges; %};
+    Ref<LayerHeightSpline> layer_height_spline()
+        %code%{ RETVAL = &THIS->layer_height_spline; %};
     Ref<Point3> size()
         %code%{ RETVAL = &THIS->size; %};
     Clone<BoundingBox> bounding_box();

xs/xsp/SlicingAdaptive.xsp

@@ -0,0 +1,16 @@
+%module{Slic3r::XS};
+
+%{
+#include <xsinit.h>
+#include "libslic3r/SlicingAdaptive.hpp"
+%}
+
+%name{Slic3r::SlicingAdaptive} class SlicingAdaptive {
+    SlicingAdaptive();
+    ~SlicingAdaptive();
+    void clear();
+    void add_mesh(TriangleMesh *mesh);
+    void prepare(coordf_t object_size);
+    float next_layer_height(coordf_t z, coordf_t quality_factor, coordf_t min_layer_height, coordf_t max_layer_height);
+    float horizontal_facet_distance(coordf_t z, coordf_t max_layer_height);
+};

xs/xsp/my.map

@@ -178,6 +178,14 @@ Ref<Layer>                 O_OBJECT_SLIC3R_T
 SupportLayer*              O_OBJECT_SLIC3R
 Ref<SupportLayer>          O_OBJECT_SLIC3R_T
 
+SlicingAdaptive*           O_OBJECT_SLIC3R
+Ref<SlicingAdaptive>       O_OBJECT_SLIC3R_T
+Clone<SlicingAdaptive>     O_OBJECT_SLIC3R_T
+
+LayerHeightSpline*         O_OBJECT_SLIC3R
+Ref<LayerHeightSpline>     O_OBJECT_SLIC3R_T
+Clone<LayerHeightSpline>   O_OBJECT_SLIC3R_T
+
 PlaceholderParser*         O_OBJECT_SLIC3R
 Ref<PlaceholderParser>     O_OBJECT_SLIC3R_T
 Clone<PlaceholderParser>   O_OBJECT_SLIC3R_T

xs/xsp/typemap.xspt

@@ -133,6 +133,14 @@
 %typemap{Layer*};
 %typemap{Ref<Layer>}{simple};
 
+%typemap{SlicingAdaptive*};
+%typemap{Ref<SlicingAdaptive>}{simple};
+%typemap{Clone<SlicingAdaptive>}{simple};
+
+%typemap{LayerHeightSpline*};
+%typemap{Ref<LayerHeightSpline>}{simple};
+%typemap{Clone<LayerHeightSpline>}{simple};
+
 %typemap{SupportLayer*};
 %typemap{Ref<SupportLayer>}{simple};