diff --git a/lib/Slic3r.pm b/lib/Slic3r.pm index 8b1f2a7a2..47ed095f3 100644 --- a/lib/Slic3r.pm +++ b/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 {}; diff --git a/lib/Slic3r/GUI.pm b/lib/Slic3r/GUI.pm index 93ae53d6c..c828c9576 100644 --- a/lib/Slic3r/GUI.pm +++ b/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; diff --git a/lib/Slic3r/GUI/OptionsGroup/Field.pm b/lib/Slic3r/GUI/OptionsGroup/Field.pm index 9d3fb7e0e..65cea5f08 100644 --- a/lib/Slic3r/GUI/OptionsGroup/Field.pm +++ b/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) = @_; diff --git a/lib/Slic3r/GUI/Plater.pm b/lib/Slic3r/GUI/Plater.pm index 8bc43b7fe..4c5e7d003 100644 --- a/lib/Slic3r/GUI/Plater.pm +++ b/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,9 +2628,16 @@ 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; diff --git a/lib/Slic3r/GUI/Plater/3DPreview.pm b/lib/Slic3r/GUI/Plater/3DPreview.pm index 8fa9188ef..ce7a5f579 100644 --- a/lib/Slic3r/GUI/Plater/3DPreview.pm +++ b/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}) { - foreach my $layer (@{$object->layers}, @{$object->support_layers}) { + if(defined $obj_idx) { # Load only given object + 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 - $self->canvas->load_print_toolpaths($self->print); - - foreach my $object (@{$self->print->objects}) { - $self->canvas->load_print_object_toolpaths($object); + if(defined $obj_idx) { # Load only one object + $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); - #my @volume_ids = $self->canvas->load_object($object->model_object); - #$self->canvas->volumes->[$_]->color->[3] = 0.2 for @volume_ids; + foreach my $object (@{$self->print->objects}) { + $self->canvas->load_print_object_toolpaths($object); + + #my @volume_ids = $self->canvas->load_object($object->model_object); + #$self->canvas->volumes->[$_]->color->[3] = 0.2 for @volume_ids; + } } $self->_loaded(1); } diff --git a/lib/Slic3r/GUI/Plater/ObjectSettingsDialog.pm b/lib/Slic3r/GUI/Plater/ObjectSettingsDialog.pm index 752b4a32c..b6d911c93 100644 --- a/lib/Slic3r/GUI/Plater/ObjectSettingsDialog.pm +++ b/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; diff --git a/lib/Slic3r/GUI/Plater/SplineControl.pm b/lib/Slic3r/GUI/Plater/SplineControl.pm new file mode 100644 index 000000000..469842a40 --- /dev/null +++ b/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; diff --git a/lib/Slic3r/GUI/PresetEditor.pm b/lib/Slic3r/GUI/PresetEditor.pm index 1eeb3c031..759a71442 100644 --- a/lib/Slic3r/GUI/PresetEditor.pm +++ b/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'); @@ -864,6 +869,12 @@ sub _update { for qw(extra_perimeters thin_walls overhangs seam_position external_perimeters_first 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() @@ -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); diff --git a/lib/Slic3r/Print/Object.pm b/lib/Slic3r/Print/Object.pm index fa15cc528..6eb9993aa 100644 --- a/lib/Slic3r/Print/Object.pm +++ b/lib/Slic3r/Print/Object.pm @@ -48,9 +48,9 @@ sub slice { return if $self->step_done(STEP_SLICE); $self->set_step_started(STEP_SLICE); $self->print->status_cb->(10, "Processing triangulated mesh"); - + $self->_slice; - + # detect slicing errors my $warning_thrown = 0; for my $i (0 .. ($self->layer_count - 1)) { diff --git a/lib/Slic3r/Print/State.pm b/lib/Slic3r/Print/State.pm index d242e3760..91d4a6ef9 100644 --- a/lib/Slic3r/Print/State.pm +++ b/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); diff --git a/lib/Slic3r/Test.pm b/lib/Slic3r/Test.pm index 8f0b203e1..f167e30aa 100644 --- a/lib/Slic3r/Test.pm +++ b/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; } diff --git a/t/adaptive_slicing.t b/t/adaptive_slicing.t new file mode 100644 index 000000000..64c3769fc --- /dev/null +++ b/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__ diff --git a/var/variable_layer_height.png b/var/variable_layer_height.png new file mode 100644 index 000000000..2fbdd6920 Binary files /dev/null and b/var/variable_layer_height.png differ diff --git a/xs/MANIFEST b/xs/MANIFEST index 3d36de159..519004743 100644 --- a/xs/MANIFEST +++ b/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 diff --git a/xs/lib/Slic3r/XS.pm b/xs/lib/Slic3r/XS.pm index 364264d67..7f269c1a8 100644 --- a/xs/lib/Slic3r/XS.pm +++ b/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 diff --git a/xs/src/BSpline/BSpline.cpp b/xs/src/BSpline/BSpline.cpp new file mode 100644 index 000000000..6cce57319 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + + + + + +/* + * 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 static inline + T abs(const T t) + { + return (t < 0) ? -t : t; + } + + template static inline + const T& min(const T& a, + const T& b) + { + return (a < b) ? a : b; + } + + template static inline + const T& max(const T& a, + const T& b) + { + return (a > b) ? a : b; + } +}; + +////////////////////////////////////////////////////////////////////// +template class Matrix : public BandedMatrix +{ + 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::operator= (e); + return *this; + } + + }; + ////////////////////////////////////////////////////////////////////// + // Our private state structure, which hides our use of some matrix + // template classes. + +template struct BSplineBaseP +{ + typedef Matrix MatrixT; + + MatrixT Q; // Holds P+Q and its factorization + std::vector X; + std::vector 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 const double BSplineBase::BoundaryConditions[3][4] = + { + // 0 1 M-1 M + { + -4, + -1, + -1, + -4 }, + { + 0, + 1, + 1, + 0 }, + { + 2, + -1, + -1, + 2 } }; +////////////////////////////////////////////////////////////////////// +template inline bool BSplineBase::Debug(int on) +{ + static bool debug = false; + if (on >= 0) + debug = (on > 0); + return debug; +} + +////////////////////////////////////////////////////////////////////// +template const double BSplineBase::BS_PI = 3.1415927; +////////////////////////////////////////////////////////////////////// +template const char * BSplineBase::ImplVersion() +{ + return ("$Id: BSpline.cpp 6352 2008-05-05 04:40:39Z martinc $"); +} + +////////////////////////////////////////////////////////////////////// +template const char * BSplineBase::IfaceVersion() +{ + return (_BSPLINEBASE_IFACE_ID); +} + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +template BSplineBase::~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 BSplineBase::BSplineBase(const BSplineBase &bb) : + K(bb.K), BC(bb.BC), OK(bb.OK), base(new BSplineBaseP(*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 BSplineBase::BSplineBase(const T *x, + int nx, + double wl, + int bc, + int num_nodes) : + NX(0), K(2), OK(false), base(new BSplineBaseP) +{ + setDomain(x, nx, wl, bc, num_nodes); +} + +////////////////////////////////////////////////////////////////////// +// Methods +template bool BSplineBase::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 double BSplineBase::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 inline double BSplineBase::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 BSpline * BSplineBase::apply(const T *y) +{ + return new BSpline (*this, y); +} +////////////////////////////////////////////////////////////////////// +/* + * Evaluate the closed basis function at node m for value x, + * using the parameters for the current boundary conditions. + */ +template double BSplineBase::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 double BSplineBase::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 double BSplineBase::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 void BSplineBase::calculateQ() +{ + Matrix &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 void BSplineBase::addP() +{ + // Add directly to Q's elements + Matrix &P = base->Q; + std::vector &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 bool BSplineBase::factor() +{ + Matrix &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 inline double BSplineBase::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 bool BSplineBase::Setup(int num_nodes) +{ + std::vector &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 const T * BSplineBase::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 std::ostream &operator<<(std::ostream &out, + const std::vector &c) +{ + for (typename std::vector::const_iterator it = c.begin(); it < c.end(); ++it) + out << *it << ", "; + out << std::endl; + return out; +} + + + + + +/* + * Original BSpline.cpp start here + */ + + +////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////// +// BSpline Class +////////////////////////////////////////////////////////////////////// + +template struct BSplineP { + std::vector spline; + std::vector A; +}; + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +/* + * This BSpline constructor constructs and sets up a new base and + * solves for the spline curve coeffiecients all at once. + */ +template BSpline::BSpline(const T *x, + int nx, + const T *y, + double wl, + int bc_type, + int num_nodes) : + BSplineBase(x, nx, wl, bc_type, num_nodes), s(new BSplineP) { + solve(y); +} +////////////////////////////////////////////////////////////////////// +/* + * Create a new spline given a BSplineBase. + */ +template BSpline::BSpline(BSplineBase &bb, + const T *y) : + BSplineBase(bb), s(new BSplineP) { + solve(y); +} +////////////////////////////////////////////////////////////////////// +/* + * (Re)calculate the spline for the given set of y values. + */ +template bool BSpline::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 &B = s->A; + std::vector &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 BSpline::~BSpline() { + delete s; +} +////////////////////////////////////////////////////////////////////// +template T BSpline::coefficient(int n) { + if (OK) + if (0 <= n && n <= M) + return s->A[n]; + return 0; +} +////////////////////////////////////////////////////////////////////// +template T BSpline::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 T BSpline::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; +//template class BSplineBase; + +/// Instantiate BSpline +template class BSpline; +//template class BSpline; diff --git a/xs/src/BSpline/BSpline.h b/xs/src/BSpline/BSpline.h new file mode 100644 index 000000000..9114385c1 --- /dev/null +++ b/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 + + +/* + * 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 BSpline; + +/* + * Opaque member structure to hide the matrix implementation. + */ +template 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 x; + vector y; + { ... } + int bc = BSplineBase::BC_ZERO_SECOND; + BSpline::Debug = true; + BSpline spline (x.begin(), x.size(), y.begin(), wl, bc); + if (spline.ok()) + { + ostream_iterator 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 and BSpline. + * 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 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 *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 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 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 BSpline : public BSplineBase +{ +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::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 &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::Debug; + using BSplineBase::Basis; + using BSplineBase::DBasis; + +protected: + + using BSplineBase::OK; + using BSplineBase::M; + using BSplineBase::NX; + using BSplineBase::DX; + using BSplineBase::base; + using BSplineBase::xmin; + using BSplineBase::xmax; + + // Our hidden state structure + BSplineP *s; + T mean; // Fit without mean and add it in later + +}; + +#endif diff --git a/xs/src/BSpline/BandedMatrix.h b/xs/src/BSpline/BandedMatrix.h new file mode 100644 index 000000000..bbe526a8b --- /dev/null +++ b/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 +#include +#include + +template class BandedMatrixRow; + + +template 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[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 & operator= (const BandedMatrix &b) + { + return Copy (*this, b); + } + + BandedMatrix & 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[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 operator[] (int row) const + { + return BandedMatrixRow(*this, row); + } + + BandedMatrixRow operator[] (int row) + { + return BandedMatrixRow(*this, row); + } + + +private: + + int top; + int bot; + int nbands; + std::vector *bands; + int N; + T out_of_bounds; + +}; + + +template +std::ostream &operator<< (std::ostream &out, const BandedMatrix &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 BandedMatrixRow +{ +public: + BandedMatrixRow (const BandedMatrix &_m, int _row) : bm(_m), i(_row) + { } + + BandedMatrixRow (BandedMatrix &_m, int _row) : bm(_m), i(_row) + { } + + ~BandedMatrixRow () {} + + typename BandedMatrix::element_type & operator[] (int j) + { + return const_cast &>(bm).element (i, j); + } + + const typename BandedMatrix::element_type & operator[] (int j) const + { + return bm.element (i, j); + } + +private: + const BandedMatrix &bm; + int i; +}; + + +/* + * Vector multiplication + */ + +template +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 +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 +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 */ + diff --git a/xs/src/BSpline/COPYRIGHT b/xs/src/BSpline/COPYRIGHT new file mode 100644 index 000000000..50c77bb13 --- /dev/null +++ b/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. diff --git a/xs/src/libslic3r/LayerHeightSpline.cpp b/xs/src/libslic3r/LayerHeightSpline.cpp new file mode 100644 index 000000000..5fe16e2f9 --- /dev/null +++ b/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 // 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 layers) +{ + this->_layers = layers; + + // generate updated layer height list from layers + this->_layer_heights.clear(); + coordf_t last_z = 0; + for (std::vector::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 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 LayerHeightSpline::getInterpolatedLayers() const +{ + std::vector 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(&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; +} + + +} diff --git a/xs/src/libslic3r/LayerHeightSpline.hpp b/xs/src/libslic3r/LayerHeightSpline.hpp new file mode 100644 index 000000000..214f2466c --- /dev/null +++ b/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 layers); + bool updateLayerHeights(std::vector 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 getOriginalLayers() const { return this->_layers; }; + std::vector 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 _layers; + std::vector _layer_heights; + std::vector _spline_layers; + std::vector _spline_layer_heights; + std::unique_ptr> _layer_height_spline; +}; + +} + +#endif diff --git a/xs/src/libslic3r/Model.cpp b/xs/src/libslic3r/Model.cpp index 287306d85..aa29858af 100644 --- a/xs/src/libslic3r/Model.cpp +++ b/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 diff --git a/xs/src/libslic3r/Model.hpp b/xs/src/libslic3r/Model.hpp index 8d594417b..db6299ee4 100644 --- a/xs/src/libslic3r/Model.hpp +++ b/xs/src/libslic3r/Model.hpp @@ -7,6 +7,7 @@ #include "Layer.hpp" #include "Point.hpp" #include "TriangleMesh.hpp" +#include "LayerHeightSpline.hpp" #include #include #include @@ -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 diff --git a/xs/src/libslic3r/Print.cpp b/xs/src/libslic3r/Print.cpp index a8119e288..5531bd4c7 100644 --- a/xs/src/libslic3r/Print.cpp +++ b/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" - || opt_key == "resolution" + } else if (opt_key == "nozzle_diameter") { + osteps.insert(posLayers); + } else if (opt_key == "resolution" || opt_key == "z_steps_per_mm") { osteps.insert(posSlice); } else if (opt_key == "avoid_crossing_perimeters" diff --git a/xs/src/libslic3r/Print.hpp b/xs/src/libslic3r/Print.hpp index 5e050ccd7..196565b57 100644 --- a/xs/src/libslic3r/Print.hpp +++ b/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; diff --git a/xs/src/libslic3r/PrintConfig.cpp b/xs/src/libslic3r/PrintConfig.cpp index 565f4ced9..431a5ea66 100644 --- a/xs/src/libslic3r/PrintConfig.cpp +++ b/xs/src/libslic3r/PrintConfig.cpp @@ -22,7 +22,26 @@ PrintConfigDef::PrintConfigDef() external_fill_pattern.enum_labels.push_back("Octagram Spiral"); 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."; diff --git a/xs/src/libslic3r/PrintConfig.hpp b/xs/src/libslic3r/PrintConfig.hpp index 1b7c435ae..b837bb735 100644 --- a/xs/src/libslic3r/PrintConfig.hpp +++ b/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 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); diff --git a/xs/src/libslic3r/PrintObject.cpp b/xs/src/libslic3r/PrintObject.cpp index 0c72d7638..10b56c37f 100644 --- a/xs/src/libslic3r/PrintObject.cpp +++ b/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 == "raft_layers" + || opt_key == "adaptive_slicing" + || 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; - if(this->_print->config.z_steps_per_mm > 0) { - coordf_t min_dz = 1 / this->_print->config.z_steps_per_mm * 4; - result = int(layer_height / min_dz + 0.5) * min_dz; - } + coordf_t result = layer_height; + if(this->_print->config.z_steps_per_mm > 0) { + coordf_t min_dz = 1 / this->_print->config.z_steps_per_mm * 4; + result = int(layer_height / min_dz + 0.5) * min_dz; + } return result > 0 ? result : layer_height; } @@ -553,64 +560,154 @@ std::vector PrintObject::generate_object_layers(coordf_t first_layer_h std::vector 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 object_extruders = this->_print->object_extruders(); - for (std::set::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)); - } - coordf_t layer_height = std::min(min_nozzle_diameter, this->config.layer_height.getFloat()); - layer_height = this->adjust_layer_height(layer_height); + for (std::set::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_layer_height = std::max(min_layer_height, this->_print->config.min_layer_height.get_at(*it_extruder)); + max_layer_height = std::min(max_layer_height, this->_print->config.max_layer_height.get_at(*it_extruder)); + + } + 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 - 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; + // Update object size at the spline object to define upper border + this->layer_height_spline.setObjectHeight(unscale(this->size.z)); + if (this->state.is_done(posLayers)) { + // layer heights are already generated, just update layers from spline + // 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. - // 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) { - coordf_t diff = result.back() - unscale(this->size.z); - int last_layer = result.size()-1; + // push modified spline object back to model + this->_model_object->layer_height_spline = this->layer_height_spline; - 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 - std::cout << new_h << std::endl; - 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 - std::cout << new_h << std::endl; - result[last_layer] = result[last_layer-1] + new_h; + // apply z-gradation (this is redundant for static layer height...) + coordf_t gradation = 1 / this->_print->config.z_steps_per_mm * 4; + if(this->_print->config.z_steps_per_mm > 0) { + coordf_t last_z = 0; + coordf_t height; + for(std::vector::iterator l = result.begin(); l != result.end(); ++l) { + height = *l - last_z; + + coordf_t gradation_effect = unscale((scale_(height)) % (scale_(gradation))); + if(gradation_effect > gradation/2 && (height + (gradation-gradation_effect)) <= max_layer_height) { // round up + height = height + (gradation-gradation_effect); + }else{ // round down + 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 PrintObject::generate_object_layers(coordf_t first_layer_h void PrintObject::_slice() { - 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 raft_height = 0; + coordf_t first_layer_height = this->config.first_layer_height.get_abs_value(this->config.layer_height.value); // take raft layers into account diff --git a/xs/src/libslic3r/SlicingAdaptive.cpp b/xs/src/libslic3r/SlicingAdaptive.cpp new file mode 100644 index 000000000..7cd25a0f4 --- /dev/null +++ b/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 face_z_span(const stl_facet *f) +{ + return std::pair( + 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_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_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 span1 = face_z_span(f1); + std::pair 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(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 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(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 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 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(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 diff --git a/xs/src/libslic3r/SlicingAdaptive.hpp b/xs/src/libslic3r/SlicingAdaptive.hpp new file mode 100644 index 000000000..be2f8f1c9 --- /dev/null +++ b/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 m_meshes; + // Collected faces of all meshes, sorted by raising Z of the bottom most face. + std::vector m_faces; + // Z component of face normals, normalized. + std::vector m_face_normal_z; +}; + +}; // namespace Slic3r + +#endif /* slic3r_SlicingAdaptive_hpp_ */ diff --git a/xs/src/perlglue.cpp b/xs/src/perlglue.cpp index a60de7252..0a194cbba 100644 --- a/xs/src/perlglue.cpp +++ b/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"); diff --git a/xs/xsp/LayerHeightSpline.xsp b/xs/xsp/LayerHeightSpline.xsp new file mode 100644 index 000000000..67d7f1e11 --- /dev/null +++ b/xs/xsp/LayerHeightSpline.xsp @@ -0,0 +1,25 @@ +%module{Slic3r::XS}; + +%{ +#include +#include "libslic3r/LayerHeightSpline.hpp" +%} + +%name{Slic3r::LayerHeightSpline} class LayerHeightSpline { + LayerHeightSpline(); + Clone clone() + %code%{ RETVAL = THIS; %}; + + void setObjectHeight(coordf_t object_height); + bool hasData(); + bool setLayers(std::vector layers) + %code%{ RETVAL = THIS->setLayers(layers); %}; + bool updateLayerHeights(std::vector heights) + %code%{ RETVAL = THIS->updateLayerHeights(heights); %}; + bool layersUpdated(); + bool layerHeightsUpdated(); + void clear(); + std::vector getOriginalLayers(); + std::vector getInterpolatedLayers(); + coordf_t getLayerHeightAt(coordf_t height); +}; diff --git a/xs/xsp/Model.xsp b/xs/xsp/Model.xsp index 272db104e..6772aee07 100644 --- a/xs/xsp/Model.xsp +++ b/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 layer_height_spline() + %code%{ RETVAL = &THIS->layer_height_spline; %}; + void set_layer_height_spline(LayerHeightSpline* spline) + %code%{ THIS->layer_height_spline = *spline; %}; + Ref origin_translation() %code%{ RETVAL = &THIS->origin_translation; %}; void set_origin_translation(Pointf3* point) diff --git a/xs/xsp/Print.xsp b/xs/xsp/Print.xsp index 59169b325..9397961dd 100644 --- a/xs/xsp/Print.xsp +++ b/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 layer_height_spline() + %code%{ RETVAL = &THIS->layer_height_spline; %}; Ref size() %code%{ RETVAL = &THIS->size; %}; Clone bounding_box(); diff --git a/xs/xsp/SlicingAdaptive.xsp b/xs/xsp/SlicingAdaptive.xsp new file mode 100644 index 000000000..09526a938 --- /dev/null +++ b/xs/xsp/SlicingAdaptive.xsp @@ -0,0 +1,16 @@ +%module{Slic3r::XS}; + +%{ +#include +#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); +}; diff --git a/xs/xsp/my.map b/xs/xsp/my.map index 89025f078..9a8bd2e3e 100644 --- a/xs/xsp/my.map +++ b/xs/xsp/my.map @@ -178,6 +178,14 @@ Ref O_OBJECT_SLIC3R_T SupportLayer* O_OBJECT_SLIC3R Ref O_OBJECT_SLIC3R_T +SlicingAdaptive* O_OBJECT_SLIC3R +Ref O_OBJECT_SLIC3R_T +Clone O_OBJECT_SLIC3R_T + +LayerHeightSpline* O_OBJECT_SLIC3R +Ref O_OBJECT_SLIC3R_T +Clone O_OBJECT_SLIC3R_T + PlaceholderParser* O_OBJECT_SLIC3R Ref O_OBJECT_SLIC3R_T Clone O_OBJECT_SLIC3R_T diff --git a/xs/xsp/typemap.xspt b/xs/xsp/typemap.xspt index 1433b828f..6a1847b94 100644 --- a/xs/xsp/typemap.xspt +++ b/xs/xsp/typemap.xspt @@ -133,6 +133,14 @@ %typemap{Layer*}; %typemap{Ref}{simple}; +%typemap{SlicingAdaptive*}; +%typemap{Ref}{simple}; +%typemap{Clone}{simple}; + +%typemap{LayerHeightSpline*}; +%typemap{Ref}{simple}; +%typemap{Clone}{simple}; + %typemap{SupportLayer*}; %typemap{Ref}{simple};