diff --git a/lib/Slic3r.pm b/lib/Slic3r.pm
index aa9aea28cf..4df83f99a0 100644
--- a/lib/Slic3r.pm
+++ b/lib/Slic3r.pm
@@ -193,7 +193,9 @@ sub thread_cleanup {
     *Slic3r::ExtrusionLoop::DESTROY         = sub {};
     *Slic3r::ExtrusionPath::DESTROY         = sub {};
     *Slic3r::ExtrusionPath::Collection::DESTROY = sub {};
+    *Slic3r::ExtrusionSimulator::DESTROY    = sub {};
     *Slic3r::Flow::DESTROY                  = sub {};
+    *Slic3r::Filler::Destroy                = sub {};
     *Slic3r::GCode::DESTROY                 = sub {};
     *Slic3r::GCode::AvoidCrossingPerimeters::DESTROY = sub {};
     *Slic3r::GCode::OozePrevention::DESTROY = sub {};
diff --git a/lib/Slic3r/Fill.pm b/lib/Slic3r/Fill.pm
index 8c63fde596..ec9038b12c 100644
--- a/lib/Slic3r/Fill.pm
+++ b/lib/Slic3r/Fill.pm
@@ -3,45 +3,25 @@ use Moo;
 
 use List::Util qw(max);
 use Slic3r::ExtrusionPath ':roles';
-use Slic3r::Fill::3DHoneycomb;
-use Slic3r::Fill::Base;
-use Slic3r::Fill::Concentric;
-use Slic3r::Fill::Honeycomb;
-use Slic3r::Fill::PlanePath;
-use Slic3r::Fill::Rectilinear;
+
 use Slic3r::Flow ':roles';
 use Slic3r::Geometry qw(X Y PI scale chained_path deg2rad);
 use Slic3r::Geometry::Clipper qw(union union_ex diff diff_ex intersection_ex offset offset2);
 use Slic3r::Surface ':types';
 
-
 has 'bounding_box' => (is => 'ro', required => 0);
 has 'fillers'   => (is => 'rw', default => sub { {} });
 
-our %FillTypes = (
-    archimedeanchords   => 'Slic3r::Fill::ArchimedeanChords',
-    rectilinear         => 'Slic3r::Fill::Rectilinear',
-    grid                => 'Slic3r::Fill::Grid',
-    flowsnake           => 'Slic3r::Fill::Flowsnake',
-    octagramspiral      => 'Slic3r::Fill::OctagramSpiral',
-    hilbertcurve        => 'Slic3r::Fill::HilbertCurve',
-    line                => 'Slic3r::Fill::Line',
-    concentric          => 'Slic3r::Fill::Concentric',
-    honeycomb           => 'Slic3r::Fill::Honeycomb',
-    '3dhoneycomb'       => 'Slic3r::Fill::3DHoneycomb',
-);
-
 sub filler {
     my $self = shift;
     my ($filler) = @_;
     
     if (!ref $self) {
-        return $FillTypes{$filler}->new;
+       return Slic3r::Filler->new_from_type($filler);
     }
-    
-    $self->fillers->{$filler} ||= $FillTypes{$filler}->new(
-        bounding_box => $self->bounding_box,
-    );
+
+    $self->fillers->{$filler} ||= Slic3r::Filler->new_from_type($filler);
+    $self->fillers->{$filler}->set_bounding_box($self->bounding_box);
     return $self->fillers->{$filler};
 }
 
@@ -227,25 +207,25 @@ sub make_fill {
                 -1, # auto width
                 $layerm->layer->object,
             );
-            $f->spacing($internal_flow->spacing);
+            $f->set_spacing($internal_flow->spacing);
             $using_internal_flow = 1;
         } else {
-            $f->spacing($flow->spacing);
+            $f->set_spacing($flow->spacing);
         }
         
-        $f->layer_id($layerm->layer->id);
-        $f->z($layerm->layer->print_z);
-        $f->angle(deg2rad($layerm->region->config->fill_angle));
-        $f->loop_clipping(scale($flow->nozzle_diameter) * &Slic3r::LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER);
+        $f->set_layer_id($layerm->layer->id);
+        $f->set_z($layerm->layer->print_z);
+        $f->set_angle(deg2rad($layerm->region->config->fill_angle));
+        $f->set_loop_clipping(scale($flow->nozzle_diameter) * &Slic3r::LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER);
         
         # apply half spacing using this flow's own spacing and generate infill
-        my @polylines = map $f->fill_surface(
-            $_,
+        my @polylines = $f->fill_surface(
+            $surface,
             density         => $density/100,
             layer_height    => $h,
-        ), @{ $surface->offset(-scale($f->spacing)/2) };
-        
+        );
         next unless @polylines;
+
         
         # calculate actual flow from spacing (which might have been adjusted by the infill
         # pattern generator)
@@ -278,7 +258,7 @@ sub make_fill {
                     mm3_per_mm  => $mm3_per_mm,
                     width       => $flow->width,
                     height      => $flow->height,
-                ), @polylines,
+                ), map @$_, @polylines,
             );
         }
     }
diff --git a/lib/Slic3r/Fill/3DHoneycomb.pm b/lib/Slic3r/Fill/3DHoneycomb.pm
deleted file mode 100644
index 3bf7e547fd..0000000000
--- a/lib/Slic3r/Fill/3DHoneycomb.pm
+++ /dev/null
@@ -1,230 +0,0 @@
-package Slic3r::Fill::3DHoneycomb;
-use Moo;
-
-extends 'Slic3r::Fill::Base';
-
-use POSIX qw(ceil fmod);
-use Slic3r::Geometry qw(scale scaled_epsilon);
-use Slic3r::Geometry::Clipper qw(intersection_pl);
-
-# require bridge flow since most of this pattern hangs in air
-sub use_bridge_flow { 1 }
-
-sub fill_surface {
-    my ($self, $surface, %params) = @_;
-    
-    my $expolygon = $surface->expolygon;
-    my $bb = $expolygon->bounding_box;
-    my $size = $bb->size;
-    
-    my $distance = scale($self->spacing) / $params{density};
-    
-    # align bounding box to a multiple of our honeycomb grid module
-    # (a module is 2*$distance since one $distance half-module is 
-    # growing while the other $distance half-module is shrinking)
-    {
-        my $min = $bb->min_point;
-        $min->translate(
-            -($bb->x_min % (2*$distance)),
-            -($bb->y_min % (2*$distance)),
-        );
-        $bb->merge_point($min);
-    }
-    
-    # generate pattern
-    my @polylines = map Slic3r::Polyline->new(@$_),
-        makeGrid(
-            scale($self->z),
-            $distance,
-            ceil($size->x / $distance) + 1,
-            ceil($size->y / $distance) + 1,  #//
-            (($self->layer_id / $surface->thickness_layers) % 2) + 1,
-        );
-    
-    # move pattern in place
-    $_->translate($bb->x_min, $bb->y_min) for @polylines;
-    
-    # clip pattern to boundaries
-    @polylines = @{intersection_pl(\@polylines, \@$expolygon)};
-    
-    # connect lines
-    unless ($params{dont_connect} || !@polylines) {  # prevent calling leftmost_point() on empty collections
-        my ($expolygon_off) = @{$expolygon->offset_ex(scaled_epsilon)};
-        my $collection = Slic3r::Polyline::Collection->new(@polylines);
-        @polylines = ();
-        foreach my $polyline (@{$collection->chained_path_from($collection->leftmost_point, 0)}) {
-            # try to append this polyline to previous one if any
-            if (@polylines) {
-                my $line = Slic3r::Line->new($polylines[-1]->last_point, $polyline->first_point);
-                if ($line->length <= 1.5*$distance && $expolygon_off->contains_line($line)) {
-                    $polylines[-1]->append_polyline($polyline);
-                    next;
-                }
-            }
-            
-            # make a clone before $collection goes out of scope
-            push @polylines, $polyline->clone;
-        }
-    }
-    
-    # TODO: return ExtrusionLoop objects to get better chained paths
-    return @polylines;
-}
-
-
-=head1 DESCRIPTION
-
-Creates a contiguous sequence of points at a specified height that make
-up a horizontal slice of the edges of a space filling truncated
-octahedron tesselation. The octahedrons are oriented so that the
-square faces are in the horizontal plane with edges parallel to the X
-and Y axes.
-
-Credits: David Eccles (gringer).
-
-=head2 makeGrid(z, gridSize, gridWidth, gridHeight, curveType)
-
-Generate a set of curves (array of array of 2d points) that describe a
-horizontal slice of a truncated regular octahedron with a specified
-grid square size.
-
-=cut
-
-sub makeGrid {
-    my ($z, $gridSize, $gridWidth, $gridHeight, $curveType) = @_;
-    my $scaleFactor = $gridSize;
-    my $normalisedZ = $z / $scaleFactor;
-    my @points = makeNormalisedGrid($normalisedZ, $gridWidth, $gridHeight, $curveType);
-    foreach my $lineRef (@points) {
-        foreach my $pointRef (@$lineRef) {
-            $pointRef->[0] *= $scaleFactor;
-            $pointRef->[1] *= $scaleFactor;
-        }
-    }
-    return @points;
-}
-
-=head1 FUNCTIONS
-=cut
-
-=head2 colinearPoints(offset, gridLength)
-
-Generate an array of points that are in the same direction as the
-basic printing line (i.e. Y points for columns, X points for rows)
-
-Note: a negative offset only causes a change in the perpendicular
-direction
-
-=cut
-
-sub colinearPoints {
-    my ($offset, $baseLocation, $gridLength) = @_;
-    
-    my @points = ();
-    push @points, $baseLocation - abs($offset/2);
-    for (my $i = 0; $i < $gridLength; $i++) {
-        push @points, $baseLocation + $i + abs($offset/2);
-        push @points, $baseLocation + ($i+1) - abs($offset/2);
-    }
-    push @points, $baseLocation + $gridLength + abs($offset/2);
-    return @points;
-}
-
-=head2 colinearPoints(offset, baseLocation, gridLength)
-
-Generate an array of points for the dimension that is perpendicular to
-the basic printing line (i.e. X points for columns, Y points for rows)
-
-=cut
-
-sub perpendPoints {
-    my ($offset, $baseLocation, $gridLength) = @_;
-    
-    my @points = ();
-    my $side = 2*(($baseLocation) % 2) - 1;
-    push @points, $baseLocation - $offset/2 * $side;
-    for (my $i = 0; $i < $gridLength; $i++) {
-        $side = 2*(($i+$baseLocation) % 2) - 1;
-        push @points, $baseLocation + $offset/2 * $side;
-        push @points, $baseLocation + $offset/2 * $side;
-    }
-    push @points, $baseLocation - $offset/2 * $side;
-    
-    return @points;
-}
-
-=head2 trim(pointArrayRef, minX, minY, maxX, maxY)
-
-Trims an array of points to specified rectangular limits. Point
-components that are outside these limits are set to the limits.
-
-=cut
-
-sub trim {
-    my ($pointArrayRef, $minX, $minY, $maxX, $maxY) = @_;
-  
-    foreach (@$pointArrayRef) {
-        $_->[0] = ($_->[0] < $minX) ? $minX : (($_->[0] > $maxX) ? $maxX : $_->[0]);
-        $_->[1] = ($_->[1] < $minY) ? $minY : (($_->[1] > $maxY) ? $maxY : $_->[1]);
-    }
-}
-
-=head2 makeNormalisedGrid(z, gridWidth, gridHeight, curveType)
-
-Generate a set of curves (array of array of 2d points) that describe a
-horizontal slice of a truncated regular octahedron with edge length 1.
-
-curveType specifies which lines to print, 1 for vertical lines
-(columns), 2 for horizontal lines (rows), and 3 for both.
-
-=cut
-
-sub makeNormalisedGrid {
-    my ($z, $gridWidth, $gridHeight, $curveType) = @_;
-    
-    ## offset required to create a regular octagram
-    my $octagramGap = 0.5;
-    
-    # sawtooth wave function for range f($z) = [-$octagramGap .. $octagramGap]
-    my $a = sqrt(2);  # period
-    my $wave = abs(fmod($z, $a) - $a/2)/$a*4 - 1;
-    my $offset = $wave * $octagramGap;
-    
-    my @points = ();
-    if (($curveType & 1) != 0) {
-        for (my $x = 0; $x <= $gridWidth; $x++) {
-            my @xPoints = perpendPoints($offset, $x, $gridHeight);
-            my @yPoints = colinearPoints($offset, 0, $gridHeight);
-            # This is essentially @newPoints = zip(@xPoints, @yPoints)
-            my @newPoints = map [ $xPoints[$_], $yPoints[$_] ], 0..$#xPoints;
-            
-            # trim points to grid edges
-            #trim(\@newPoints, 0, 0, $gridWidth, $gridHeight);
-            
-            if ($x % 2 == 0){
-                push @points, [ @newPoints ];
-            } else {
-                push @points, [ reverse @newPoints ];
-            }
-        }
-    }
-    if (($curveType & 2) != 0) {
-        for (my $y = 0; $y <= $gridHeight; $y++) {
-            my @xPoints = colinearPoints($offset, 0, $gridWidth);
-            my @yPoints = perpendPoints($offset, $y, $gridWidth);
-            my @newPoints = map [ $xPoints[$_], $yPoints[$_] ], 0..$#xPoints;
-            
-            # trim points to grid edges
-            #trim(\@newPoints, 0, 0, $gridWidth, $gridHeight);
-            
-            if ($y % 2 == 0) {
-                push @points, [ @newPoints ];
-            } else {
-                push @points, [ reverse @newPoints ];
-            }
-        }
-    }
-    return @points;
-}
-
-1;
diff --git a/lib/Slic3r/Fill/Base.pm b/lib/Slic3r/Fill/Base.pm
deleted file mode 100644
index 75c8e03e66..0000000000
--- a/lib/Slic3r/Fill/Base.pm
+++ /dev/null
@@ -1,91 +0,0 @@
-package Slic3r::Fill::Base;
-use Moo;
-
-has 'layer_id'            => (is => 'rw');
-has 'z'                   => (is => 'rw'); # in unscaled coordinates
-has 'angle'               => (is => 'rw'); # in radians, ccw, 0 = East
-has 'spacing'             => (is => 'rw'); # in unscaled coordinates
-has 'loop_clipping'       => (is => 'rw', default => sub { 0 }); # in scaled coordinates
-has 'bounding_box'        => (is => 'ro', required => 0);  # Slic3r::Geometry::BoundingBox object
-
-sub adjust_solid_spacing {
-    my $self = shift;
-    my %params = @_;
-    
-    my $number_of_lines = int($params{width} / $params{distance}) + 1;
-    return $params{distance} if $number_of_lines <= 1;
-    
-    my $extra_space = $params{width} % $params{distance};
-    return $params{distance} + $extra_space / ($number_of_lines - 1);
-}
-
-sub no_sort { 0 }
-sub use_bridge_flow { 0 }
-
-
-package Slic3r::Fill::WithDirection;
-use Moo::Role;
-
-use Slic3r::Geometry qw(PI rad2deg);
-
-sub angles () { [0, PI/2] }
-
-sub infill_direction {
-    my $self = shift;
-    my ($surface) = @_;
-    
-    if (!defined $self->angle) {
-        warn "Using undefined infill angle";
-        $self->angle(0);
-    }
-    
-    # set infill angle
-    my (@rotate);
-    $rotate[0] = $self->angle;
-    $rotate[1] = $self->bounding_box
-        ? $self->bounding_box->center
-        : $surface->expolygon->bounding_box->center;
-    my $shift = $rotate[1]->clone;
-    
-    if (defined $self->layer_id) {
-        # alternate fill direction
-        my $layer_num = $self->layer_id / $surface->thickness_layers;
-        my $angle = $self->angles->[$layer_num % @{$self->angles}];
-        $rotate[0] = $self->angle + $angle if $angle;
-    }
-        
-    # use bridge angle
-    if ($surface->bridge_angle >= 0) {
-        Slic3r::debugf "Filling bridge with angle %d\n", rad2deg($surface->bridge_angle);
-        $rotate[0] = $surface->bridge_angle;
-    }
-    
-    $rotate[0] += PI/2;
-    $shift->rotate(@rotate);
-    return [\@rotate, $shift];
-}
-
-# this method accepts any object that implements rotate() and translate()
-sub rotate_points {
-    my $self = shift;
-    my ($expolygon, $rotate_vector) = @_;
-    
-    # rotate points
-    my ($rotate, $shift) = @$rotate_vector;
-    $rotate = [ -$rotate->[0], $rotate->[1] ];
-    $expolygon->rotate(@$rotate);
-    $expolygon->translate(@$shift);
-}
-
-sub rotate_points_back {
-    my $self = shift;
-    my ($paths, $rotate_vector) = @_;
-    
-    my ($rotate, $shift) = @$rotate_vector;
-    $shift = [ map -$_, @$shift ];
-    
-    $_->translate(@$shift) for @$paths;
-    $_->rotate(@$rotate) for @$paths;
-}
-
-1;
diff --git a/lib/Slic3r/Fill/Concentric.pm b/lib/Slic3r/Fill/Concentric.pm
deleted file mode 100644
index ca1837c4ea..0000000000
--- a/lib/Slic3r/Fill/Concentric.pm
+++ /dev/null
@@ -1,57 +0,0 @@
-package Slic3r::Fill::Concentric;
-use Moo;
-
-extends 'Slic3r::Fill::Base';
-
-use Slic3r::Geometry qw(scale unscale X);
-use Slic3r::Geometry::Clipper qw(offset offset2 union_pt_chained);
-
-sub no_sort { 1 }
-
-sub fill_surface {
-    my $self = shift;
-    my ($surface, %params) = @_;
-    
-    # no rotation is supported for this infill pattern
-    
-    my $expolygon = $surface->expolygon;
-    my $bounding_box = $expolygon->bounding_box;
-    
-    my $min_spacing = scale($self->spacing);
-    my $distance = $min_spacing / $params{density};
-    
-    if ($params{density} == 1 && !$params{dont_adjust}) {
-        $distance = $self->adjust_solid_spacing(
-            width       => $bounding_box->size->[X],
-            distance    => $distance,
-        );
-        $self->spacing(unscale $distance);
-    }
-    
-    my @loops = my @last = map $_->clone, @$expolygon;
-    while (@last) {
-        push @loops, @last = @{offset2(\@last, -($distance + 0.5*$min_spacing), +0.5*$min_spacing)};
-    }
-    
-    # generate paths from the outermost to the innermost, to avoid 
-    # adhesion problems of the first central tiny loops
-    @loops = map Slic3r::Polygon->new(@$_),
-        reverse @{union_pt_chained(\@loops)};
-    
-    # split paths using a nearest neighbor search
-    my @paths = ();
-    my $last_pos = Slic3r::Point->new(0,0);
-    foreach my $loop (@loops) {
-        push @paths, $loop->split_at_index($last_pos->nearest_point_index(\@$loop));
-        $last_pos = $paths[-1]->last_point;
-    }
-    
-    # clip the paths to prevent the extruder from getting exactly on the first point of the loop
-    $_->clip_end($self->loop_clipping) for @paths;
-    @paths = grep $_->is_valid, @paths;  # remove empty paths (too short, thus eaten by clipping)
-    
-    # TODO: return ExtrusionLoop objects to get better chained paths
-    return @paths;
-}
-
-1;
diff --git a/lib/Slic3r/Fill/Honeycomb.pm b/lib/Slic3r/Fill/Honeycomb.pm
deleted file mode 100644
index b0fbd65ff2..0000000000
--- a/lib/Slic3r/Fill/Honeycomb.pm
+++ /dev/null
@@ -1,129 +0,0 @@
-package Slic3r::Fill::Honeycomb;
-use Moo;
-
-extends 'Slic3r::Fill::Base';
-with qw(Slic3r::Fill::WithDirection);
-
-has 'cache'         => (is => 'rw', default => sub {{}});
-
-use Slic3r::Geometry qw(PI X Y MIN MAX scale scaled_epsilon);
-use Slic3r::Geometry::Clipper qw(intersection intersection_pl);
-
-sub angles () { [0, PI/3, PI/3*2] }
-
-sub fill_surface {
-    my $self = shift;
-    my ($surface, %params) = @_;
-    
-    my $rotate_vector = $self->infill_direction($surface);
-    
-    # cache hexagons math
-    my $cache_id = sprintf "d%s_s%s", $params{density}, $self->spacing;
-    my $m;
-    if (!($m = $self->cache->{$cache_id})) {
-        $m = $self->cache->{$cache_id} = {};
-        my $min_spacing = scale($self->spacing);
-        $m->{distance} = $min_spacing / $params{density};
-        $m->{hex_side} = $m->{distance} / (sqrt(3)/2);
-        $m->{hex_width} = $m->{distance} * 2;  # $m->{hex_width} == $m->{hex_side} * sqrt(3);
-        my $hex_height = $m->{hex_side} * 2;
-        $m->{pattern_height} = $hex_height + $m->{hex_side};
-        $m->{y_short} = $m->{distance} * sqrt(3)/3;
-        $m->{x_offset} = $min_spacing / 2;
-        $m->{y_offset} = $m->{x_offset} * sqrt(3)/3;
-        $m->{hex_center} = Slic3r::Point->new($m->{hex_width}/2, $m->{hex_side});
-    }
-    
-    my @polygons = ();
-    {
-        # adjust actual bounding box to the nearest multiple of our hex pattern
-        # and align it so that it matches across layers
-        
-        my $bounding_box = $surface->expolygon->bounding_box;
-        {
-            # rotate bounding box according to infill direction
-            my $bb_polygon = $bounding_box->polygon;
-            $bb_polygon->rotate($rotate_vector->[0][0], $m->{hex_center});
-            $bounding_box = $bb_polygon->bounding_box;
-            
-            # extend bounding box so that our pattern will be aligned with other layers
-            # $bounding_box->[X1] and [Y1] represent the displacement between new bounding box offset and old one
-            $bounding_box->merge_point(Slic3r::Point->new(
-                $bounding_box->x_min - ($bounding_box->x_min % $m->{hex_width}),
-                $bounding_box->y_min - ($bounding_box->y_min % $m->{pattern_height}),
-            ));
-        }
-        
-        my $x = $bounding_box->x_min;
-        while ($x <= $bounding_box->x_max) {
-            my $p = [];
-            
-            my @x = ($x + $m->{x_offset}, $x + $m->{distance} - $m->{x_offset});
-            for (1..2) {
-                @$p = reverse @$p; # turn first half upside down
-                my @p = ();
-                for (my $y = $bounding_box->y_min; $y <= $bounding_box->y_max; $y += $m->{y_short} + $m->{hex_side} + $m->{y_short} + $m->{hex_side}) {
-                    push @$p,
-                        [ $x[1], $y + $m->{y_offset} ],
-                        [ $x[0], $y + $m->{y_short} - $m->{y_offset} ],
-                        [ $x[0], $y + $m->{y_short} + $m->{hex_side} + $m->{y_offset} ],
-                        [ $x[1], $y + $m->{y_short} + $m->{hex_side} + $m->{y_short} - $m->{y_offset} ],
-                        [ $x[1], $y + $m->{y_short} + $m->{hex_side} + $m->{y_short} + $m->{hex_side} + $m->{y_offset} ];
-                }
-                @x = map $_ + $m->{distance}, reverse @x; # draw symmetrical pattern
-                $x += $m->{distance};
-            }
-            
-            push @polygons, Slic3r::Polygon->new(@$p);
-        }
-        
-        $_->rotate(-$rotate_vector->[0][0], $m->{hex_center}) for @polygons;
-    }
-    
-    my @paths;
-    if ($params{complete} || 1) {
-        # we were requested to complete each loop;
-        # in this case we don't try to make more continuous paths
-        @paths = map $_->split_at_first_point,
-            @{intersection([ $surface->p ], \@polygons)};
-        
-    } else {
-        # consider polygons as polylines without re-appending the initial point:
-        # this cuts the last segment on purpose, so that the jump to the next 
-        # path is more straight
-        @paths = @{intersection_pl(
-            [ map Slic3r::Polyline->new(@$_), @polygons ],
-            [ @{$surface->expolygon} ],
-        )};
-        
-        # connect paths
-        if (@paths) {  # prevent calling leftmost_point() on empty collections
-            my $collection = Slic3r::Polyline::Collection->new(@paths);
-            @paths = ();
-            foreach my $path (@{$collection->chained_path_from($collection->leftmost_point, 0)}) {
-                if (@paths) {
-                    # distance between first point of this path and last point of last path
-                    my $distance = $paths[-1]->last_point->distance_to($path->first_point);
-                    
-                    if ($distance <= $m->{hex_width}) {
-                        $paths[-1]->append_polyline($path);
-                        next;
-                    }
-                }
-                
-                # make a clone before $collection goes out of scope
-                push @paths, $path->clone;
-            }
-        }
-        
-        # clip paths again to prevent connection segments from crossing the expolygon boundaries
-        @paths = @{intersection_pl(
-            \@paths,
-            [ map @$_, @{$surface->expolygon->offset_ex(scaled_epsilon)} ],
-        )};
-    }
-    
-    return @paths;
-}
-
-1;
diff --git a/lib/Slic3r/Fill/PlanePath.pm b/lib/Slic3r/Fill/PlanePath.pm
deleted file mode 100644
index 556835ec4b..0000000000
--- a/lib/Slic3r/Fill/PlanePath.pm
+++ /dev/null
@@ -1,118 +0,0 @@
-package Slic3r::Fill::PlanePath;
-use Moo;
-
-extends 'Slic3r::Fill::Base';
-with qw(Slic3r::Fill::WithDirection);
-
-use Slic3r::Geometry qw(scale X1 Y1 X2 Y2);
-use Slic3r::Geometry::Clipper qw(intersection_pl);
-
-sub angles () { [0] }
-sub multiplier () { 1 }
-
-sub process_polyline {}
-
-sub fill_surface {
-    my $self = shift;
-    my ($surface, %params) = @_;
-    
-    # rotate polygons
-    my $expolygon = $surface->expolygon->clone;
-    my $rotate_vector = $self->infill_direction($surface);
-    $self->rotate_points($expolygon, $rotate_vector);
-    
-    my $distance_between_lines = scale($self->spacing) / $params{density} * $self->multiplier;
-    
-    # align infill across layers using the object's bounding box
-    my $bb_polygon = $self->bounding_box->polygon;
-    $self->rotate_points($bb_polygon, $rotate_vector);
-    my $bounding_box = $bb_polygon->bounding_box;
-    
-    (ref $self) =~ /::([^:]+)$/;
-    my $path = "Math::PlanePath::$1"->new;
-    
-    my $translate = Slic3r::Point->new(0,0);  # vector
-    if ($path->x_negative || $path->y_negative) {
-        # if the curve extends on both positive and negative coordinate space,
-        # center our expolygon around origin
-        $translate = $bounding_box->center->negative;
-    } else {
-        # if the curve does not extend in negative coordinate space,
-        # move expolygon entirely in positive coordinate space
-        $translate = $bounding_box->min_point->negative;
-    }
-    $expolygon->translate(@$translate);
-    $bounding_box->translate(@$translate);
-    
-    my ($n_lo, $n_hi) = $path->rect_to_n_range(
-        map { $_ / $distance_between_lines }
-            @{$bounding_box->min_point},
-            @{$bounding_box->max_point},
-    );
-    
-    my $polyline = Slic3r::Polyline->new(
-        map [ map { $_ * $distance_between_lines } $path->n_to_xy($_) ], ($n_lo..$n_hi)
-    );
-    return {} if @$polyline <= 1;
-    
-    $self->process_polyline($polyline, $bounding_box);
-    
-    my @paths = @{intersection_pl([$polyline], \@$expolygon)};
-    
-    if (0) {
-        require "Slic3r/SVG.pm";
-        Slic3r::SVG::output("fill.svg",
-            no_arrows       => 1,
-            polygons        => \@$expolygon,
-            green_polygons  => [ $bounding_box->polygon ],
-            polylines       => [ $polyline ],
-            red_polylines   => \@paths,
-        );
-    }
-    
-    # paths must be repositioned and rotated back
-    $_->translate(@{$translate->negative}) for @paths;
-    $self->rotate_points_back(\@paths, $rotate_vector);
-    
-    return @paths;
-}
-
-
-package Slic3r::Fill::ArchimedeanChords;
-use Moo;
-extends 'Slic3r::Fill::PlanePath';
-use Math::PlanePath::ArchimedeanChords;
-
-
-package Slic3r::Fill::Flowsnake;
-use Moo;
-extends 'Slic3r::Fill::PlanePath';
-use Math::PlanePath::Flowsnake;
-use Slic3r::Geometry qw(X);
-
-# Sorry, this fill is currently broken.
-
-sub process_polyline {
-    my $self = shift;
-    my ($polyline, $bounding_box) = @_;
-    
-    $_->[X] += $bounding_box->center->[X] for @$polyline;
-}
-
-
-package Slic3r::Fill::HilbertCurve;
-use Moo;
-extends 'Slic3r::Fill::PlanePath';
-use Math::PlanePath::HilbertCurve;
-
-
-package Slic3r::Fill::OctagramSpiral;
-use Moo;
-extends 'Slic3r::Fill::PlanePath';
-use Math::PlanePath::OctagramSpiral;
-
-sub multiplier () { sqrt(2) }
-
-
-
-1;
diff --git a/lib/Slic3r/Fill/Rectilinear.pm b/lib/Slic3r/Fill/Rectilinear.pm
deleted file mode 100644
index 0922ff7714..0000000000
--- a/lib/Slic3r/Fill/Rectilinear.pm
+++ /dev/null
@@ -1,168 +0,0 @@
-package Slic3r::Fill::Rectilinear;
-use Moo;
-
-extends 'Slic3r::Fill::Base';
-with qw(Slic3r::Fill::WithDirection);
-
-has '_min_spacing'          => (is => 'rw');
-has '_line_spacing'         => (is => 'rw');
-has '_diagonal_distance'    => (is => 'rw');
-has '_line_oscillation'     => (is => 'rw');
-
-use Slic3r::Geometry qw(scale unscale scaled_epsilon);
-use Slic3r::Geometry::Clipper qw(intersection_pl);
-
-sub horizontal_lines { 0 }
-
-sub fill_surface {
-    my $self = shift;
-    my ($surface, %params) = @_;
-    
-    # rotate polygons so that we can work with vertical lines here
-    my $expolygon = $surface->expolygon->clone;
-    my $rotate_vector = $self->infill_direction($surface);
-    $self->rotate_points($expolygon, $rotate_vector);
-    
-    $self->_min_spacing(scale $self->spacing);
-    $self->_line_spacing($self->_min_spacing / $params{density});
-    $self->_diagonal_distance($self->_line_spacing * 2);
-    $self->_line_oscillation($self->_line_spacing - $self->_min_spacing);  # only for Line infill
-    my $bounding_box = $expolygon->bounding_box;
-    
-    # define flow spacing according to requested density
-    if ($params{density} == 1 && !$params{dont_adjust}) {
-        $self->_line_spacing($self->adjust_solid_spacing(
-            width       => $bounding_box->size->x,
-            distance    => $self->_line_spacing,
-        ));
-        $self->spacing(unscale $self->_line_spacing);
-    } else {
-        # extend bounding box so that our pattern will be aligned with other layers
-        $bounding_box->merge_point(Slic3r::Point->new(
-            $bounding_box->x_min - ($bounding_box->x_min % $self->_line_spacing),
-            $bounding_box->y_min - ($bounding_box->y_min % $self->_line_spacing),
-        ));
-    }
-    
-    # generate the basic pattern
-    my $x_max = $bounding_box->x_max + scaled_epsilon;
-    my @lines  = ();
-    for (my $x = $bounding_box->x_min; $x <= $x_max; $x += $self->_line_spacing) {
-        push @lines, $self->_line($#lines, $x, $bounding_box->y_min, $bounding_box->y_max);
-    }
-    if ($self->horizontal_lines) {
-        my $y_max = $bounding_box->y_max + scaled_epsilon;
-        for (my $y = $bounding_box->y_min; $y <= $y_max; $y += $self->_line_spacing) {
-            push @lines, Slic3r::Polyline->new(
-                [$bounding_box->x_min, $y],
-                [$bounding_box->x_max, $y],
-            );
-        }
-    }
-    
-    # clip paths against a slightly larger expolygon, so that the first and last paths
-    # are kept even if the expolygon has vertical sides
-    # the minimum offset for preventing edge lines from being clipped is scaled_epsilon;
-    # however we use a larger offset to support expolygons with slightly skewed sides and 
-    # not perfectly straight
-    my @polylines = @{intersection_pl(\@lines, $expolygon->offset(+scale 0.02))};
-    
-    my $extra = $self->_min_spacing * &Slic3r::INFILL_OVERLAP_OVER_SPACING;
-    foreach my $polyline (@polylines) {
-        my ($first_point, $last_point) = @$polyline[0,-1];
-        if ($first_point->y > $last_point->y) { #>
-            ($first_point, $last_point) = ($last_point, $first_point);
-        }
-        $first_point->set_y($first_point->y - $extra);  #--
-        $last_point->set_y($last_point->y + $extra);    #++
-    }
-    
-    # connect lines
-    unless ($params{dont_connect} || !@polylines) {  # prevent calling leftmost_point() on empty collections
-        # offset the expolygon by max(min_spacing/2, extra)
-        my ($expolygon_off) = @{$expolygon->offset_ex($self->_min_spacing/2)};
-        my $collection = Slic3r::Polyline::Collection->new(@polylines);
-        @polylines = ();
-        
-        foreach my $polyline (@{$collection->chained_path_from($collection->leftmost_point, 0)}) {
-            if (@polylines) {
-                my $first_point = $polyline->first_point;
-                my $last_point = $polylines[-1]->last_point;
-                my @distance = map abs($first_point->$_ - $last_point->$_), qw(x y);
-                
-                # TODO: we should also check that both points are on a fill_boundary to avoid 
-                # connecting paths on the boundaries of internal regions
-                if ($self->_can_connect(@distance) && $expolygon_off->contains_line(Slic3r::Line->new($last_point, $first_point))) {
-                    $polylines[-1]->append_polyline($polyline);
-                    next;
-                }
-            }
-            
-            # make a clone before $collection goes out of scope
-            push @polylines, $polyline->clone;
-        }
-    }
-    
-    # paths must be rotated back
-    $self->rotate_points_back(\@polylines, $rotate_vector);
-    
-    return @polylines;
-}
-
-sub _line {
-    my ($self, $i, $x, $y_min, $y_max) = @_;
-    
-    return Slic3r::Polyline->new(
-        [$x, $y_min],
-        [$x, $y_max],
-    );
-}
-
-sub _can_connect {
-    my ($self, $dist_X, $dist_Y) = @_;
-    
-    return $dist_X <= $self->_diagonal_distance
-        && $dist_Y <= $self->_diagonal_distance;
-}
-
-
-package Slic3r::Fill::Line;
-use Moo;
-extends 'Slic3r::Fill::Rectilinear';
-
-use Slic3r::Geometry qw(scaled_epsilon);
-
-sub _line {
-    my ($self, $i, $x, $y_min, $y_max) = @_;
-    
-    if ($i % 2) {
-        return Slic3r::Polyline->new(
-            [$x - $self->_line_oscillation, $y_min],
-            [$x + $self->_line_oscillation, $y_max],
-        );
-    } else {
-        return Slic3r::Polyline->new(
-            [$x, $y_min],
-            [$x, $y_max],
-        );
-    }
-}
-
-sub _can_connect {
-    my ($self, $dist_X, $dist_Y) = @_;
-    
-    my $TOLERANCE = 10 * scaled_epsilon;
-    return ($dist_X >= ($self->_line_spacing - $self->_line_oscillation) - $TOLERANCE)
-        && ($dist_X <= ($self->_line_spacing + $self->_line_oscillation) + $TOLERANCE)
-        && $dist_Y <= $self->_diagonal_distance;
-}
-
-
-package Slic3r::Fill::Grid;
-use Moo;
-extends 'Slic3r::Fill::Rectilinear';
-
-sub angles () { [0] }
-sub horizontal_lines { 1 }
-
-1;
diff --git a/lib/Slic3r/GUI.pm b/lib/Slic3r/GUI.pm
index 7f4a028a5f..1c811fa789 100644
--- a/lib/Slic3r/GUI.pm
+++ b/lib/Slic3r/GUI.pm
@@ -40,6 +40,7 @@ use Wx 0.9901 qw(:bitmap :dialog :icon :id :misc :systemsettings :toplevelwindow
     :filedialog :font);
 use Wx::Event qw(EVT_IDLE EVT_COMMAND);
 use base 'Wx::App';
+#use base 'Wx::AppConsole';
 
 use constant FILE_WILDCARDS => {
     known   => 'Known files (*.stl, *.obj, *.amf, *.xml)|*.stl;*.STL;*.obj;*.OBJ;*.amf;*.AMF;*.xml;*.XML',
diff --git a/lib/Slic3r/GUI/MainFrame.pm b/lib/Slic3r/GUI/MainFrame.pm
index d97ad22efc..4c04066bb6 100644
--- a/lib/Slic3r/GUI/MainFrame.pm
+++ b/lib/Slic3r/GUI/MainFrame.pm
@@ -294,8 +294,13 @@ sub _init_menubar {
         $self->_append_menu_item($helpMenu, "&About Slic3r", 'Show about dialog', sub {
             wxTheApp->about;
         });
+        if (Slic3r::GUI::debugged()) {
+            $self->_append_menu_item($helpMenu, "&Debug", 'Break to debugger', sub {
+                Slic3r::GUI::break_to_debugger();
+            });
+        }
     }
-    
+
     # menubar
     # assign menubar to frame after appending items, otherwise special items
     # will not be handled correctly
diff --git a/lib/Slic3r/GUI/Plater/2DToolpaths.pm b/lib/Slic3r/GUI/Plater/2DToolpaths.pm
index ee59c2980b..66d669b11c 100644
--- a/lib/Slic3r/GUI/Plater/2DToolpaths.pm
+++ b/lib/Slic3r/GUI/Plater/2DToolpaths.pm
@@ -1,10 +1,13 @@
+# 2D preview of the tool paths of a single layer, using a thin line.
+# OpenGL is used to render the paths.
+
 package Slic3r::GUI::Plater::2DToolpaths;
 use strict;
 use warnings;
 use utf8;
 
 use Slic3r::Print::State ':steps';
-use Wx qw(:misc :sizer :slider :statictext wxWHITE);
+use Wx qw(:misc :sizer :slider :statictext :keycode wxWHITE wxWANTS_CHARS);
 use Wx::Event qw(EVT_SLIDER EVT_KEY_DOWN);
 use base qw(Wx::Panel Class::Accessor);
 
@@ -14,7 +17,7 @@ sub new {
     my $class = shift;
     my ($parent, $print) = @_;
     
-    my $self = $class->SUPER::new($parent, -1, wxDefaultPosition);
+    my $self = $class->SUPER::new($parent, -1, wxDefaultPosition, wxDefaultSize, wxWANTS_CHARS);
     $self->SetBackgroundColour(wxWHITE);
     
     # init GUI elements
@@ -48,17 +51,21 @@ sub new {
     });
     EVT_KEY_DOWN($canvas, sub {
         my ($s, $event) = @_;
-        
         my $key = $event->GetKeyCode;
-        if ($key == 85 || $key == 315) {
+        if ($key == 85 || $key == WXK_LEFT) {
+            # Keys: 'D' or WXK_LEFT
             $slider->SetValue($slider->GetValue + 1);
             $self->set_z($self->{layers_z}[$slider->GetValue]);
-        } elsif ($key == 68 || $key == 317) {
+        } elsif ($key == 68 || $key == WXK_RIGHT) {
+            # Keys: 'U' or WXK_RIGHT
             $slider->SetValue($slider->GetValue - 1);
             $self->set_z($self->{layers_z}[$slider->GetValue]);
+        } elsif ($key >= 49 && $key <= 55) {
+            # Keys: '1' to '3'
+            $canvas->set_simulation_mode($key - 49);
         }
     });
-    
+
     $self->SetSizer($sizer);
     $self->SetMinSize($self->GetSize);
     $sizer->SetSizeHints($self);
@@ -132,6 +139,10 @@ __PACKAGE__->mk_accessors(qw(
     _zoom
     _camera_target
     _drag_start_xy
+    _texture_name
+    _texture_size
+    _extrusion_simulator
+    _simulation_mode
 ));
 
 # make OpenGL::Array thread-safe
@@ -155,7 +166,13 @@ sub new {
 
     # 2D point in model space
     $self->_camera_target(Slic3r::Pointf->new(0,0));
-    
+
+    # Texture for the extrusion simulator. The texture will be allocated / reallocated on Resize.
+    $self->_texture_name(0);
+    $self->_texture_size(Slic3r::Point->new(0,0));
+    $self->_extrusion_simulator(Slic3r::ExtrusionSimulator->new());
+    $self->_simulation_mode(0);
+
     EVT_PAINT($self, sub {
         my $dc = Wx::PaintDC->new($self);
         $self->Render($dc);
@@ -210,6 +227,21 @@ sub new {
     return $self;
 }
 
+sub Destroy {
+    my ($self) = @_;
+    
+    # Deallocate the OpenGL resources.
+    my $context = $self->GetContext;
+    if ($context and $self->texture_id) {
+        $self->SetCurrent($context);
+        glDeleteTextures(1, ($self->texture_id));
+        $self->SetCurrent(0);
+        $self->texture_id(0);
+        $self->texture_size(new Slic3r::Point(0, 0));
+    }
+    return $self->SUPER::Destroy;
+}
+
 sub mouse_event {
     my ($self, $e) = @_;
     
@@ -278,6 +310,14 @@ sub set_z {
     $self->Refresh;
 }
 
+sub set_simulation_mode
+{
+    my ($self, $mode) = @_;
+    $self->_simulation_mode($mode);
+    $self->_dirty(1);
+    $self->Refresh;
+}
+
 sub Render {
     my ($self, $dc) = @_;
     
@@ -299,7 +339,43 @@ sub Render {
     glDisable(GL_DEPTH_TEST);
     glMatrixMode(GL_MODELVIEW);
     glLoadIdentity();
-    
+
+    if ($self->_simulation_mode and $self->_texture_name and $self->_texture_size->x() > 0 and $self->_texture_size->y() > 0) {
+        $self->_simulate_extrusion();
+        my ($x, $y) = $self->GetSizeWH;
+        glEnable(GL_TEXTURE_2D);
+        glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,GL_REPLACE);
+        glBindTexture(GL_TEXTURE_2D, $self->_texture_name);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+        glTexImage2D_c(GL_TEXTURE_2D,
+            0,                        # level (0 normal, heighr is form mip-mapping)
+            GL_RGBA,                  # internal format
+            $self->_texture_size->x(), $self->_texture_size->y(),
+            0,                        # border 
+            GL_RGBA,                  # format RGBA color data
+            GL_UNSIGNED_BYTE,         # unsigned byte data
+            $self->_extrusion_simulator->image_ptr()); # ptr to texture data
+        glMatrixMode(GL_PROJECTION);
+        glPushMatrix();
+        glLoadIdentity();
+        glOrtho(0, 1, 0, 1, 0, 1);
+        glBegin(GL_QUADS);
+        glTexCoord2f(0, 0);
+        glVertex2f(0, 0);
+        glTexCoord2f($x/$self->_texture_size->x(), 0);
+        glVertex2f(1, 0);
+        glTexCoord2f($x/$self->_texture_size->x(), $y/$self->_texture_size->y());
+        glVertex2f(1, 1);
+        glTexCoord2f(0, $y/$self->_texture_size->y());
+        glVertex2f(0, 1);
+        glEnd();
+        glPopMatrix();
+        glBindTexture(GL_TEXTURE_2D, 0);
+    }
+
     # anti-alias
     if (0) {
         glEnable(GL_LINE_SMOOTH);
@@ -308,8 +384,9 @@ sub Render {
         glHint(GL_POLYGON_SMOOTH_HINT, GL_DONT_CARE);
     }
     
+    # Tesselator triangulates polygons with holes on the fly for the rendering purposes only.
     my $tess;
-    if (!(&Wx::wxMSW && $OpenGL::VERSION < 0.6704)) {
+    if ($self->_simulation_mode() == 0 and !(&Wx::wxMSW && $OpenGL::VERSION < 0.6704)) {
         # We can't use the GLU tesselator on MSW with older OpenGL versions
         # because of an upstream bug:
         # http://sourceforge.net/p/pogl/bugs/16/
@@ -448,8 +525,42 @@ sub _draw_path {
             glVertex2f(@{$line->a});
             glVertex2f(@{$line->b});
             glEnd();
+		}
+	}
+}
+
+sub _simulate_extrusion {
+    my ($self) = @_;
+    $self->_extrusion_simulator->reset_accumulator();
+    foreach my $layer (@{$self->layers}) {
+        if (abs($layer->print_z - $self->z) < epsilon) {
+            my $object = $layer->object;
+            my @shifts = (defined $object) ? @{$object->_shifted_copies} : (Slic3r::Point->new(0, 0));
+            foreach my $layerm (@{$layer->regions}) {
+                my @extrusions = ();
+                if ($object->step_done(STEP_PERIMETERS)) {
+                    push @extrusions, @$_ for @{$layerm->perimeters};
+                }
+                if ($object->step_done(STEP_INFILL)) {
+                    push @extrusions, @$_ for @{$layerm->fills};
+                }
+                foreach my $extrusion_entity (@extrusions) {
+                    my @paths = $extrusion_entity->isa('Slic3r::ExtrusionLoop')
+                        ? @$extrusion_entity
+                        : ($extrusion_entity);
+                    foreach my $path (@paths) {
+                        print "width: ", $path->width, 
+                              " height: ", $path->height,
+                              " mm3_per_mm: ", $path->mm3_per_mm,
+                              " height2: ", $path->mm3_per_mm / $path->height,
+                              "\n";
+                        $self->_extrusion_simulator->extrude_to_accumulator($path, $_, $self->_simulation_mode()) for @shifts;
+                    }
+                }
+            }
         }
     }
+    $self->_extrusion_simulator->evaluate_accumulator($self->_simulation_mode());
 }
 
 sub InitGL {
@@ -457,6 +568,10 @@ sub InitGL {
  
     return if $self->init;
     return unless $self->GetContext;
+
+    my $texture_id = 0;
+    ($texture_id) = glGenTextures_p(1);
+    $self->_texture_name($texture_id);
     $self->init(1);
 }
 
@@ -489,6 +604,33 @@ sub Resize {
     
     $self->SetCurrent($self->GetContext);
     my ($x, $y) = $self->GetSizeWH;
+
+    if ($self->_texture_size->x() < $x or $self->_texture_size->y() < $y) {
+        # Allocate a large enough OpenGL texture with power of 2 dimensions.
+        $self->_texture_size->set_x(1) if ($self->_texture_size->x() == 0);
+        $self->_texture_size->set_y(1) if ($self->_texture_size->y() == 0);
+        $self->_texture_size->set_x($self->_texture_size->x() * 2) while ($self->_texture_size->x() < $x);
+        $self->_texture_size->set_y($self->_texture_size->y() * 2) while ($self->_texture_size->y() < $y);
+        #print "screen size ", $x, "x", $y;
+        #print "texture size ", $self->_texture_size->x(), "x", $self->_texture_size->y();
+        # Initialize an empty texture.
+        glBindTexture(GL_TEXTURE_2D, $self->_texture_name);
+        if (1) {
+        glTexImage2D_c(GL_TEXTURE_2D,
+            0,                        # level (0 normal, heighr is form mip-mapping)
+            GL_RGBA,                  # internal format
+            $self->_texture_size->x(), $self->_texture_size->y(),
+            0,                        # border 
+            GL_RGBA,                  # format RGBA color data
+            GL_UNSIGNED_BYTE,         # unsigned byte data
+            0);                       # ptr to texture data
+        }
+        glBindTexture(GL_TEXTURE_2D, 0);
+        $self->_extrusion_simulator->set_image_size($self->_texture_size);
+    }
+    $self->_extrusion_simulator->set_viewport(Slic3r::Geometry::BoundingBox->new_from_points(
+        [Slic3r::Point->new(0, 0), Slic3r::Point->new($x, $y)]));
+
     glViewport(0, 0, $x, $y);
     
     glMatrixMode(GL_PROJECTION);
@@ -544,10 +686,18 @@ sub Resize {
         $x2 = $x1 + $new_x;
     }
     glOrtho($x1, $x2, $y1, $y2, 0, 1);
-    
+
+    # Set the adjusted bounding box at the extrusion simulator.
+    #print "Scene bbox ", $bb->x_min, ",", $bb->y_min, " ", $bb->x_max, ",", $bb->y_max, "\n";
+    #print "Setting simulator bbox ", $x1, ",", $y1, " ", $x2, ",", $y2, "\n";
+    $self->_extrusion_simulator->set_bounding_box(
+        Slic3r::Geometry::BoundingBox->new_from_points(
+            [Slic3r::Point->new($x1, $y1), Slic3r::Point->new($x2, $y2)]));
+
     glMatrixMode(GL_MODELVIEW);
 }
 
+# Thick line drawing is not used anywhere. Probably not tested?
 sub line {
     my (
         $x1, $y1, $x2, $y2,     # coordinates of the line
diff --git a/lib/Slic3r/Print/SupportMaterial.pm b/lib/Slic3r/Print/SupportMaterial.pm
index eea6397afa..acfa5efee6 100644
--- a/lib/Slic3r/Print/SupportMaterial.pm
+++ b/lib/Slic3r/Print/SupportMaterial.pm
@@ -675,8 +675,8 @@ sub generate_toolpaths {
         
         # interface and contact infill
         if (@$interface || @$contact_infill) {
-            $fillers{interface}->angle($interface_angle);
-            $fillers{interface}->spacing($_interface_flow->spacing);
+            $fillers{interface}->set_angle($interface_angle);
+            $fillers{interface}->set_spacing($_interface_flow->spacing);
             
             # find centerline of the external loop
             $interface = offset2($interface, +scaled_epsilon, -(scaled_epsilon + $_interface_flow->scaled_width/2));
@@ -702,7 +702,7 @@ sub generate_toolpaths {
             
             my @paths = ();
             foreach my $expolygon (@{union_ex($interface)}) {
-                my @p = $fillers{interface}->fill_surface(
+                my $polylines = $fillers{interface}->fill_surface(
                     Slic3r::Surface->new(expolygon => $expolygon, surface_type => S_TYPE_INTERNAL),
                     density      => $interface_density,
                     layer_height => $layer->height,
@@ -711,12 +711,12 @@ sub generate_toolpaths {
                 my $mm3_per_mm = $_interface_flow->mm3_per_mm;
                 
                 push @paths, map Slic3r::ExtrusionPath->new(
-                    polyline    => Slic3r::Polyline->new(@$_),
+                    polyline    => $_,
                     role        => EXTR_ROLE_SUPPORTMATERIAL_INTERFACE,
                     mm3_per_mm  => $mm3_per_mm,
                     width       => $_interface_flow->width,
                     height      => $layer->height,
-                ), @p;
+                ), @$polylines,
             }
             
             $layer->support_interface_fills->append(@paths);
@@ -725,11 +725,11 @@ sub generate_toolpaths {
         # support or flange
         if (@$base) {
             my $filler = $fillers{support};
-            $filler->angle($angles[ ($layer_id) % @angles ]);
+            $filler->set_angle($angles[ ($layer_id) % @angles ]);
             
             # We don't use $base_flow->spacing because we need a constant spacing
             # value that guarantees that all layers are correctly aligned.
-            $filler->spacing($flow->spacing);
+            $filler->set_spacing($flow->spacing);
             
             my $density     = $support_density;
             my $base_flow   = $_flow;
@@ -742,13 +742,13 @@ sub generate_toolpaths {
             # base flange
             if ($layer_id == 0) {
                 $filler = $fillers{interface};
-                $filler->angle($self->object_config->support_material_angle + 90);
+                $filler->set_angle($self->object_config->support_material_angle + 90);
                 $density        = 0.5;
                 $base_flow      = $self->first_layer_flow;
                 
                 # use the proper spacing for first layer as we don't need to align
                 # its pattern to the other layers
-                $filler->spacing($base_flow->spacing);
+                $filler->set_spacing($base_flow->spacing);
             } else {
                 # draw a perimeter all around support infill
                 # TODO: use brim ordering algorithm
@@ -767,7 +767,7 @@ sub generate_toolpaths {
             
             my $mm3_per_mm = $base_flow->mm3_per_mm;
             foreach my $expolygon (@$to_infill) {
-                my @p = $filler->fill_surface(
+                my $polylines = $filler->fill_surface(
                     Slic3r::Surface->new(expolygon => $expolygon, surface_type => S_TYPE_INTERNAL),
                     density     => $density,
                     layer_height => $layer->height,
@@ -775,12 +775,12 @@ sub generate_toolpaths {
                 );
                 
                 push @paths, map Slic3r::ExtrusionPath->new(
-                    polyline    => Slic3r::Polyline->new(@$_),
+                    polyline    => $_,
                     role        => EXTR_ROLE_SUPPORTMATERIAL,
                     mm3_per_mm  => $mm3_per_mm,
                     width       => $base_flow->width,
                     height      => $layer->height,
-                ), @p;
+                ), @$polylines;
             }
             
             $layer->support_fills->append(@paths);
diff --git a/slic3r.sublime-project b/slic3r.sublime-project
new file mode 100644
index 0000000000..9deb0e273b
--- /dev/null
+++ b/slic3r.sublime-project
@@ -0,0 +1,61 @@
+{
+	"build_systems":
+	[
+		{
+			"name": "List",
+   			//"file_regex": " at ([^-\\s]*) line ([0-9]*)",
+//   			"file_regex": " at (D\\:\\/src\\/Slic3r\\/.*?) line ([0-9]*)",
+			"shell_cmd": "ls -l"
+		},
+		{
+			"name": "Run",
+			"working_dir": "$project_path",
+   			"file_regex": " at (.*?) line ([0-9]*)",
+			"shell_cmd": "perl slic3r.pl --gui \"..\\Slic3r-tests\\gap fill torture 20 -rt.stl\""
+		},
+		{
+			"name": "full",
+   			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
+			"shell_cmd": "perl Build.pl"
+		},
+		{
+			"name": "xs",
+			"working_dir": "$project_path/xs",
+  			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
+			"shell_cmd": "perl Build install"
+		},
+		{
+			"name": "xs & run",
+			"working_dir": "$project_path/xs",
+  			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
+			"shell_cmd": "perl Build install & cd .. & perl slic3r.pl --gui \"..\\Slic3r-tests\\star3-big2.stl\""
+		}
+	],
+	"folders":
+	[
+		{
+			"path": "."
+		}
+	],
+
+    "settings":
+    {
+        "sublimegdb_workingdir": "${folder:${project_path:run}}",
+  		// NOTE: You MUST provide --interpreter=mi for the plugin to work
+//    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi -ex 'target localhost:2345'",
+//    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi perl --args perl slic3r.pl",
+//    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi perl --args slic3r.pl ",
+//    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi -e C:\\Strawberry\\perl\\bin\\perl.exe -s C:\\Strawberry\\perl\\site\\lib\\auto\\Slic3r\\XS\\XS.xs.dll --args perl slic3r.pl -j 1 --gui D:\\src\\Slic3r-tests\\star3-big.stl",
+    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi perl.exe --args perl slic3r.pl -j 1 --gui", //  D:\\src\\Slic3r-tests\\star3-big.stl",
+//    	"sublimegdb_commandline": "D:\\Qt\\Tools\\mingw492_32\\bin\\gdb.exe --interpreter=mi -x slic3r.gdb",
+//    	"arguments": "slic3r -j 1 --gui ../Slic3r-tests/star3-big.stl",
+//	  	"arguments": "../slic3r.pl -j 1 --gui",
+//        "sublimegdb_exec_cmd": "-exec-continue",
+    	
+    	// Add "pending breakpoints" for symbols that are dynamically loaded from
+    	// external shared libraries
+    	"debug_ext" : true,
+    	"run_after_init": false,
+    	"close_views": false
+    }
+}
diff --git a/xs/Build.PL b/xs/Build.PL
index 5413952bb2..bdf5abbd35 100644
--- a/xs/Build.PL
+++ b/xs/Build.PL
@@ -40,10 +40,15 @@ if (defined $ENV{BOOST_INCLUDEDIR}) {
         push @boost_include, $ENV{BOOST_DIR};
     }
 } else {
-    # Boost library was not defined by the environment.
-    # Try to guess at some default paths.
-    if ($mswin) {
-        for my $path (glob('C:\dev\boost*\include'), glob ('C:\boost*\include')) {
+    push @boost_include, grep { -d $_ }
+        qw(/opt/local/include /usr/local/include /opt/include),
+        qw(/usr/include C:\Boost\include);
+    push @boost_libs, grep { -d $_ }
+        qw(/opt/local/lib /usr/local/lib /opt/lib /usr/lib),
+        qw(C:\Boost\lib /lib);
+    
+    if ($^O eq 'MSWin32') {
+        for my $path (glob('C:\dev\boost*'), glob ('C:\boost*'), glob ('d:\src\boost*')) {
             push @boost_include, $path;
         }
         if (! @boost_include) {
diff --git a/xs/MANIFEST b/xs/MANIFEST
index 0f541515da..2fbd9f0cea 100644
--- a/xs/MANIFEST
+++ b/xs/MANIFEST
@@ -28,6 +28,22 @@ src/libslic3r/ExtrusionEntity.cpp
 src/libslic3r/ExtrusionEntity.hpp
 src/libslic3r/ExtrusionEntityCollection.cpp
 src/libslic3r/ExtrusionEntityCollection.hpp
+src/libslic3r/ExtrusionSimulator.cpp
+src/libslic3r/ExtrusionSimulator.hpp
+src/libslic3r/Fill/FillBase.cpp
+src/libslic3r/Fill/FillBase.hpp
+src/libslic3r/Fill/FillConcentric.cpp
+src/libslic3r/Fill/FillConcentric.hpp
+src/libslic3r/Fill/FillHoneycomb.cpp
+src/libslic3r/Fill/FillHoneycomb.hpp
+src/libslic3r/Fill/Fill3DHoneycomb.cpp
+src/libslic3r/Fill/Fill3DHoneycomb.hpp
+src/libslic3r/Fill/FillPlanePath.cpp
+src/libslic3r/Fill/FillPlanePath.hpp
+src/libslic3r/Fill/FillRectilinear.cpp
+src/libslic3r/Fill/FillRectilinear.hpp
+src/libslic3r/Fill/FillRectilinear2.cpp
+src/libslic3r/Fill/FillRectilinear2.hpp
 src/libslic3r/Flow.cpp
 src/libslic3r/Flow.hpp
 src/libslic3r/GCode.cpp
@@ -129,6 +145,8 @@ xsp/Extruder.xsp
 xsp/ExtrusionEntityCollection.xsp
 xsp/ExtrusionLoop.xsp
 xsp/ExtrusionPath.xsp
+xsp/ExtrusionSimulator.xsp
+xsp/Filler.xsp
 xsp/Flow.xsp
 xsp/GCode.xsp
 xsp/GCodeSender.xsp
diff --git a/xs/lib/Slic3r/XS.pm b/xs/lib/Slic3r/XS.pm
index 4f53966deb..44ac57cdb8 100644
--- a/xs/lib/Slic3r/XS.pm
+++ b/xs/lib/Slic3r/XS.pm
@@ -123,6 +123,26 @@ sub clone {
     );
 }
 
+package Slic3r::ExtrusionSimulator;
+
+sub new {
+    my ($class, %args) = @_;
+    return $class->_new();
+}
+
+package Slic3r::Filler;
+
+sub fill_surface {
+    my ($self, $surface, %args) = @_;
+    $self->set_width($args{width}) if defined($args{width});
+    $self->set_density($args{density}) if defined($args{density});
+    $self->set_distance($args{distance}) if defined($args{distance});
+    $self->set_dont_connect($args{dont_connect}) if defined($args{dont_connect});
+    $self->set_dont_adjust($args{dont_adjust}) if defined($args{dont_adjust});
+    $self->set_complete($args{complete}) if defined($args{complete});
+    return $self->_fill_surface($surface);
+}
+
 package Slic3r::Flow;
 
 sub new {
@@ -215,6 +235,8 @@ for my $class (qw(
         Slic3r::ExtrusionLoop
         Slic3r::ExtrusionPath
         Slic3r::ExtrusionPath::Collection
+        Slic3r::ExtrusionSimulator
+        Slic3r::Filler
         Slic3r::Flow
         Slic3r::GCode
         Slic3r::GCode::AvoidCrossingPerimeters
diff --git a/xs/src/libslic3r/BoundingBox.hpp b/xs/src/libslic3r/BoundingBox.hpp
index 533a19b56e..4d0c1fe361 100644
--- a/xs/src/libslic3r/BoundingBox.hpp
+++ b/xs/src/libslic3r/BoundingBox.hpp
@@ -21,6 +21,7 @@ class BoundingBoxBase
     bool defined;
     
     BoundingBoxBase() : defined(false) {};
+    BoundingBoxBase(const PointClass &pmin, const PointClass &pmax) : min(pmin), max(pmax) {}
     BoundingBoxBase(const std::vector<PointClass> &points);
     void merge(const PointClass &point);
     void merge(const std::vector<PointClass> &points);
@@ -37,6 +38,7 @@ class BoundingBox3Base : public BoundingBoxBase<PointClass>
 {
     public:
     BoundingBox3Base() : BoundingBoxBase<PointClass>() {};
+    BoundingBox3Base(const PointClass &pmin, const PointClass &pmax) : BoundingBoxBase<PointClass>(pmin, pmax) {}
     BoundingBox3Base(const std::vector<PointClass> &points);
     void merge(const PointClass &point);
     void merge(const std::vector<PointClass> &points);
@@ -54,6 +56,7 @@ class BoundingBox : public BoundingBoxBase<Point>
     Polygon polygon() const;
     
     BoundingBox() : BoundingBoxBase<Point>() {};
+    BoundingBox(const Point &pmin, const Point &pmax) : BoundingBoxBase<Point>(pmin, pmax) {};
     BoundingBox(const Points &points) : BoundingBoxBase<Point>(points) {};
     BoundingBox(const Lines &lines);
 };
@@ -65,15 +68,39 @@ class BoundingBox3  : public BoundingBox3Base<Point3> {};
 class BoundingBoxf : public BoundingBoxBase<Pointf> {
     public:
     BoundingBoxf() : BoundingBoxBase<Pointf>() {};
+    BoundingBoxf(const Pointf &pmin, const Pointf &pmax) : BoundingBoxBase<Pointf>(pmin, pmax) {};
     BoundingBoxf(const std::vector<Pointf> &points) : BoundingBoxBase<Pointf>(points) {};
 };
 
 class BoundingBoxf3 : public BoundingBox3Base<Pointf3> {
     public:
     BoundingBoxf3() : BoundingBox3Base<Pointf3>() {};
+    BoundingBoxf3(const Pointf3 &pmin, const Pointf3 &pmax) : BoundingBox3Base<Pointf3>(pmin, pmax) {};
     BoundingBoxf3(const std::vector<Pointf3> &points) : BoundingBox3Base<Pointf3>(points) {};
 };
 
+template<typename VT>
+inline bool operator==(const BoundingBoxBase<VT> &bb1, const BoundingBoxBase<VT> &bb2)
+{
+    return bb1.min == bb2.min && bb1.max == bb2.max;
+}
+
+template<typename VT>
+inline bool operator!=(const BoundingBoxBase<VT> &bb1, const BoundingBoxBase<VT> &bb2)
+{
+    return !(bb1 == bb2);
+}
+
+template<typename VT>
+inline bool empty(const BoundingBoxBase<VT> &bb)
+{
+    return bb.min.x > bb.max.y || bb.min.y > bb.max.y;
+}
+
+template<typename VT>
+inline bool empty(const BoundingBox3Base<VT> &bb)
+{
+    return bb.min.x > bb.max.x || bb.min.y > bb.max.y || bb.min.z > bb.max.z;}
 }
 
 #endif
diff --git a/xs/src/libslic3r/ClipperUtils.cpp b/xs/src/libslic3r/ClipperUtils.cpp
index aea0ef5959..a1297e574c 100644
--- a/xs/src/libslic3r/ClipperUtils.cpp
+++ b/xs/src/libslic3r/ClipperUtils.cpp
@@ -594,6 +594,10 @@ void union_pt_chained(const Slic3r::Polygons &subject, Slic3r::Polygons* retval,
 {
     ClipperLib::PolyTree pt;
     union_pt(subject, &pt, safety_offset_);
+    if (&subject == retval)
+        // It is safe to use the same variable for input and output, because this function makes
+        // a temporary copy of the results.
+        retval->clear();
     traverse_pt(pt.Childs, retval);
 }
 
diff --git a/xs/src/libslic3r/ClipperUtils.hpp b/xs/src/libslic3r/ClipperUtils.hpp
index f3ab6158e4..e7396af76e 100644
--- a/xs/src/libslic3r/ClipperUtils.hpp
+++ b/xs/src/libslic3r/ClipperUtils.hpp
@@ -35,47 +35,47 @@ void scaleClipperPolygons(ClipperLib::Paths &polygons, const double scale);
 
 // offset Polygons
 void offset(const Slic3r::Polygons &polygons, ClipperLib::Paths* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 void offset(const Slic3r::Polygons &polygons, Slic3r::Polygons* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 Slic3r::Polygons offset(const Slic3r::Polygons &polygons, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 
 // offset Polylines
 void offset(const Slic3r::Polylines &polylines, ClipperLib::Paths* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
     double miterLimit = 3);
 void offset(const Slic3r::Polylines &polylines, Slic3r::Polygons* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
     double miterLimit = 3);
 void offset(const Slic3r::Surface &surface, Slic3r::Surfaces* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtSquare, 
     double miterLimit = 3);
 
 void offset(const Slic3r::Polygons &polygons, Slic3r::ExPolygons* retval, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 Slic3r::ExPolygons offset_ex(const Slic3r::Polygons &polygons, const float delta,
-    double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 
 void offset2(const Slic3r::Polygons &polygons, ClipperLib::Paths* retval, const float delta1,
-    const float delta2, double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    const float delta2, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 void offset2(const Slic3r::Polygons &polygons, Slic3r::Polygons* retval, const float delta1,
-    const float delta2, double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    const float delta2, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 Slic3r::Polygons offset2(const Slic3r::Polygons &polygons, const float delta1,
-    const float delta2, double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    const float delta2, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 void offset2(const Slic3r::Polygons &polygons, Slic3r::ExPolygons* retval, const float delta1,
-    const float delta2, double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    const float delta2, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 Slic3r::ExPolygons offset2_ex(const Slic3r::Polygons &polygons, const float delta1,
-    const float delta2, double scale = 100000, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
+    const float delta2, double scale = CLIPPER_OFFSET_SCALE, ClipperLib::JoinType joinType = ClipperLib::jtMiter, 
     double miterLimit = 3);
 
 template <class T>
diff --git a/xs/src/libslic3r/ExPolygon.cpp b/xs/src/libslic3r/ExPolygon.cpp
index 3640d39188..1347f4a947 100644
--- a/xs/src/libslic3r/ExPolygon.cpp
+++ b/xs/src/libslic3r/ExPolygon.cpp
@@ -52,6 +52,15 @@ ExPolygon::translate(double x, double y)
     }
 }
 
+void
+ExPolygon::rotate(double angle)
+{
+    contour.rotate(angle);
+    for (Polygons::iterator it = holes.begin(); it != holes.end(); ++it) {
+        (*it).rotate(angle);
+    }
+}
+
 void
 ExPolygon::rotate(double angle, const Point &center)
 {
diff --git a/xs/src/libslic3r/ExPolygon.hpp b/xs/src/libslic3r/ExPolygon.hpp
index 69fe9fb0ee..53a3fce1db 100644
--- a/xs/src/libslic3r/ExPolygon.hpp
+++ b/xs/src/libslic3r/ExPolygon.hpp
@@ -20,6 +20,7 @@ class ExPolygon
     operator Polygons() const;
     void scale(double factor);
     void translate(double x, double y);
+    void rotate(double angle);
     void rotate(double angle, const Point &center);
     double area() const;
     bool is_valid() const;
@@ -45,6 +46,32 @@ class ExPolygon
     std::string dump_perl() const;
 };
 
+inline Polygons to_polygons(const ExPolygons &src)
+{
+    Polygons polygons;
+    for (ExPolygons::const_iterator it = src.begin(); it != src.end(); ++it) {
+        polygons.push_back(it->contour);
+        for (Polygons::const_iterator ith = it->holes.begin(); ith != it->holes.end(); ++ith) {
+            polygons.push_back(*ith);
+        }
+    }
+    return polygons;
+}
+
+#if SLIC3R_CPPVER > 11
+inline Polygons to_polygons(ExPolygons &&src)
+{
+    Polygons polygons;
+    for (ExPolygons::const_iterator it = src.begin(); it != src.end(); ++it) {
+        polygons.push_back(std::move(it->contour));
+        for (Polygons::const_iterator ith = it->holes.begin(); ith != it->holes.end(); ++ith) {
+            polygons.push_back(std::move(*ith));
+        }
+    }
+    return polygons;
+}
+#endif
+
 }
 
 // start Boost
diff --git a/xs/src/libslic3r/ExtrusionSimulator.cpp b/xs/src/libslic3r/ExtrusionSimulator.cpp
new file mode 100644
index 0000000000..2dbef9a25e
--- /dev/null
+++ b/xs/src/libslic3r/ExtrusionSimulator.cpp
@@ -0,0 +1,1047 @@
+// Optimize the extrusion simulator to the bones.
+#pragma GCC optimize ("O3")
+#undef SLIC3R_DEBUG
+#define NDEBUG
+
+#include <cmath>
+#include <cassert>
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/box.hpp>
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+
+#include <boost/multi_array.hpp>
+
+#include "ExtrusionSimulator.hpp"
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626433832795
+#endif
+
+// Replacement for a template alias.
+// Shorthand for the point_xy.
+template<typename T>
+struct V2
+{
+	typedef boost::geometry::model::d2::point_xy<T> Type;
+};
+
+// Replacement for a template alias.
+// Shorthand for the point with a cartesian coordinate system.
+template<typename T>
+struct V3
+{
+	typedef boost::geometry::model::point<T, 3, boost::geometry::cs::cartesian> Type;
+};
+
+// Replacement for a template alias.
+// Shorthand for the point with a cartesian coordinate system.
+template<typename T>
+struct V4
+{
+	typedef boost::geometry::model::point<T, 4, boost::geometry::cs::cartesian> Type;
+};
+
+typedef V2<int   >::Type V2i;
+typedef V2<float >::Type V2f;
+typedef V2<double>::Type V2d;
+
+// Used for an RGB color.
+typedef V3<unsigned char>::Type V3uc;
+// Used for an RGBA color.
+typedef V4<unsigned char>::Type V4uc;
+
+typedef boost::geometry::model::box<V2i> B2i;
+typedef boost::geometry::model::box<V2f> B2f;
+typedef boost::geometry::model::box<V2d> B2d;
+
+typedef boost::multi_array<unsigned char, 2> 	 A2uc;
+typedef boost::multi_array<int   		, 2> 	 A2i;
+typedef boost::multi_array<float 		, 2> 	 A2f;
+typedef boost::multi_array<double		, 2> 	 A2d;
+
+template<typename T>
+inline void operator+=(
+	      boost::geometry::model::d2::point_xy<T> &v1, 
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	boost::geometry::add_point(v1, v2);
+}
+
+template<typename T>
+inline void operator-=(
+	      boost::geometry::model::d2::point_xy<T> &v1, 
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	boost::geometry::subtract_point(v1, v2);
+}
+
+template<typename T>
+inline void operator*=(boost::geometry::model::d2::point_xy<T> &v, const T c)
+{
+	boost::geometry::multiply_value(v, c);
+}
+
+template<typename T>
+inline void operator/=(boost::geometry::model::d2::point_xy<T> &v, const T c)
+{
+	boost::geometry::divide_value(v, c);
+}
+
+template<typename T>
+inline typename boost::geometry::model::d2::point_xy<T> operator+(
+	const boost::geometry::model::d2::point_xy<T> &v1,
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	boost::geometry::model::d2::point_xy<T> out(v1);
+	out += v2;
+	return out;
+}
+
+template<typename T>
+inline boost::geometry::model::d2::point_xy<T> operator-(
+	const boost::geometry::model::d2::point_xy<T> &v1,
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	boost::geometry::model::d2::point_xy<T> out(v1);
+	out -= v2;
+	return out;
+}
+
+template<typename T>
+inline boost::geometry::model::d2::point_xy<T> operator*(
+	const boost::geometry::model::d2::point_xy<T> &v, const T c)
+{
+	boost::geometry::model::d2::point_xy<T> out(v);
+	out *= c;
+	return out;
+}
+
+template<typename T>
+inline typename boost::geometry::model::d2::point_xy<T> operator*(
+	const T c, const boost::geometry::model::d2::point_xy<T> &v)
+{
+	boost::geometry::model::d2::point_xy<T> out(v);
+	out *= c;
+	return out;
+}
+
+template<typename T>
+inline typename boost::geometry::model::d2::point_xy<T> operator/(
+	const boost::geometry::model::d2::point_xy<T> &v, const T c)
+{
+	boost::geometry::model::d2::point_xy<T> out(v);
+	out /= c;
+	return out;
+}
+
+template<typename T>
+inline T dot(
+	const boost::geometry::model::d2::point_xy<T> &v1,
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	return boost::geometry::dot_product(v1, v2);
+}
+
+template<typename T>
+inline T dot(const boost::geometry::model::d2::point_xy<T> &v)
+{
+	return boost::geometry::dot_product(v, v);
+}
+
+template <typename T>
+inline T cross(
+	const boost::geometry::model::d2::point_xy<T> &v1, 
+	const boost::geometry::model::d2::point_xy<T> &v2)
+{
+	return v1.x() * v2.y() - v2.x() * v1.y();
+}
+
+// Euclidian measure
+template<typename T>
+inline T l2(const boost::geometry::model::d2::point_xy<T> &v)
+{
+	return std::sqrt(dot(v));
+}
+
+// Euclidian measure
+template<typename T>
+inline T mag(const boost::geometry::model::d2::point_xy<T> &v)
+{
+	return l2(v);
+}
+
+template<typename T>
+inline T lerp(T start, T end, T alpha)
+{
+	return start * (T(1.) - alpha) + end * alpha;
+}
+
+template<typename T>
+inline T clamp(T low, T high, T x)
+{
+	return std::max<T>(low, std::min<T>(high, x));
+}
+
+template<typename T>
+inline T sqr(T x)
+{
+	return x * x;
+}
+
+template<typename T>
+inline T dist2_to_line(
+	const boost::geometry::model::d2::point_xy<T> &p0,
+	const boost::geometry::model::d2::point_xy<T> &p1,
+	const boost::geometry::model::d2::point_xy<T> &px)
+{
+	boost::geometry::model::d2::point_xy<T> v  = p1 - p0;
+	boost::geometry::model::d2::point_xy<T> vx = px - p0;
+	T 									    l  = dot(v);
+	T 									    t  = dot(v, vx);
+	if (l != T(0) && t > T(0.)) {
+		t /= l;
+		vx = px - ((t > T(1.)) ? p1 : (p0 + t * v));
+	}
+	return dot(vx);
+}
+
+// Intersect a circle with a line segment.
+// Returns number of intersection points.
+template<typename T>
+int line_circle_intersection(
+	const boost::geometry::model::d2::point_xy<T>	&p0,
+	const boost::geometry::model::d2::point_xy<T>	&p1,
+	const boost::geometry::model::d2::point_xy<T>	&center,
+	const T 										 radius,
+	boost::geometry::model::d2::point_xy<T>			 intersection[2])
+{
+	typedef typename V2<T>::Type V2T;
+	V2T v  = p1 - p0;
+	V2T vc = p0 - center;
+	T   a = dot(v);
+	T   b = T(2.) * dot(vc, v);
+	T   c = dot(vc) - radius * radius;
+	T   d = b * b - T(4.) * a * c;
+
+	if (d < T(0))
+		// The circle misses the ray.
+		return 0;
+
+	int n = 0;
+	if (d == T(0)) {
+		// The circle touches the ray at a single tangent point.
+		T t = - b / (T(2.) * a);
+		if (t >= T(0.) && t <= T(1.))
+			intersection[n ++] = p0 + t * v;
+	} else {
+		// The circle intersects the ray in two points.
+		d = sqrt(d);
+		T t = (- b - d) / (T(2.) * a);
+		if (t >= T(0.) && t <= T(1.))
+			intersection[n ++] = p0 + t * v;
+		t = (- b + d) / (T(2.) * a);
+		if (t >= T(0.) && t <= T(1.))
+			intersection[n ++] = p0 + t * v;
+	}
+	return n;
+}
+
+// Sutherland–Hodgman clipping of a rectangle against an AABB.
+// Expects the first 4 points of rect to be filled at the beginning.
+// The clipping may produce up to 8 points.
+// Returns the number of resulting points.
+template<typename T>
+int clip_rect_by_AABB(
+	boost::geometry::model::d2::point_xy<T>			   					         rect[8], 
+	const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T> > &aabb)
+{
+	typedef typename V2<T>::Type V2T;
+	V2T  result[8];
+	int  nin  = 4;
+	int  nout = 0;
+	V2T *in   = rect;
+	V2T *out  = result;
+	// Clip left
+	{
+		const V2T *S    = in + nin - 1;
+		T          left = aabb.min_corner().x();
+		for (int i = 0; i < nin; ++i) {
+			const V2T &E = in[i];
+			if (E.x() == left) {
+				out[nout++] = E;
+			}
+			else if (E.x() > left) {
+				// E is inside the AABB.
+				if (S->x() < left) {
+					// S is outside the AABB. Calculate an intersection point.
+					T t = (left - S->x()) / (E.x() - S->x());
+					out[nout++] = V2T(left, S->y() + t * (E.y() - S->y()));
+				}
+				out[nout++] = E;
+			}
+			else if (S->x() > left) {
+				// S is inside the AABB, E is outside the AABB.
+				T t = (left - S->x()) / (E.x() - S->x());
+				out[nout++] = V2T(left, S->y() + t * (E.y() - S->y()));
+			}
+			S = &E;
+		}
+		assert(nout <= 8);
+	}
+	// Clip bottom
+	{
+		std::swap(in, out);
+		nin = nout;
+		nout = 0;
+		const V2T *S      = in + nin - 1;
+		T          bottom = aabb.min_corner().y();
+		for (int i = 0; i < nin; ++i) {
+			const V2T &E = in[i];
+			if (E.y() == bottom) {
+				out[nout++] = E;
+			}
+			else if (E.y() > bottom) {
+				// E is inside the AABB.
+				if (S->y() < bottom) {
+					// S is outside the AABB. Calculate an intersection point.
+					T t = (bottom - S->y()) / (E.y() - S->y());
+					out[nout++] = V2T(S->x() + t * (E.x() - S->x()), bottom);
+				}
+				out[nout++] = E;
+			}
+			else if (S->y() > bottom) {
+				// S is inside the AABB, E is outside the AABB.
+				T t = (bottom - S->y()) / (E.y() - S->y());
+				out[nout++] = V2T(S->x() + t * (E.x() - S->x()), bottom);
+			}
+			S = &E;
+		}
+		assert(nout <= 8);
+	}
+	// Clip right
+	{
+		std::swap(in, out);
+		nin = nout;
+		nout = 0;
+		const V2T *S = in + nin - 1;
+		T right = aabb.max_corner().x();
+		for (int i = 0; i < nin; ++i) {
+			const V2T &E = in[i];
+			if (E.x() == right) {
+				out[nout++] = E;
+			}
+			else if (E.x() < right) {
+				// E is inside the AABB.
+				if (S->x() > right) {
+					// S is outside the AABB. Calculate an intersection point.
+					T t = (right - S->x()) / (E.x() - S->x());
+					out[nout++] = V2T(right, S->y() + t * (E.y() - S->y()));
+				}
+				out[nout++] = E;
+			}
+			else if (S->x() < right) {
+				// S is inside the AABB, E is outside the AABB.
+				T t = (right - S->x()) / (E.x() - S->x());
+				out[nout++] = V2T(right, S->y() + t * (E.y() - S->y()));
+			}
+			S = &E;
+		}
+		assert(nout <= 8);
+	}
+	// Clip top
+	{
+		std::swap(in, out);
+		nin = nout;
+		nout = 0;
+		const V2T *S = in + nin - 1;
+		T top = aabb.max_corner().y();
+		for (int i = 0; i < nin; ++i) {
+			const V2T &E = in[i];
+			if (E.y() == top) {
+				out[nout++] = E;
+			}
+			else if (E.y() < top) {
+				// E is inside the AABB.
+				if (S->y() > top) {
+					// S is outside the AABB. Calculate an intersection point.
+					T t = (top - S->y()) / (E.y() - S->y());
+					out[nout++] = V2T(S->x() + t * (E.x() - S->x()), top);
+				}
+				out[nout++] = E;
+			}
+			else if (S->y() < top) {
+				// S is inside the AABB, E is outside the AABB.
+				T t = (top - S->y()) / (E.y() - S->y());
+				out[nout++] = V2T(S->x() + t * (E.x() - S->x()), top);
+			}
+			S = &E;
+		}
+		assert(nout <= 8);
+	}
+
+	assert(nout <= 8);
+	return nout;
+}
+
+// Calculate area of the circle x AABB intersection.
+// The calculation is approximate in a way, that the circular segment
+// intersecting the cell is approximated by its chord (a linear segment).
+template<typename T>
+int clip_circle_by_AABB(
+	const boost::geometry::model::d2::point_xy<T>							    &center,
+	const T 																	 radius,
+	const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T> > &aabb,
+	boost::geometry::model::d2::point_xy<T>			   					         result[8],
+	bool											   					         result_arc[8])
+{
+	typedef typename V2<T>::Type V2T;
+
+	V2T rect[4] = {
+		aabb.min_corner(),
+		V2T(aabb.max_corner().x(), aabb.min_corner().y()),
+		aabb.max_corner(),
+		V2T(aabb.min_corner().x(), aabb.max_corner().y())
+	};
+
+	int  bits_corners = 0;
+	T    r2 = sqr(radius);
+	for (int i = 0; i < 4; ++ i, bits_corners <<= 1)
+		bits_corners |= dot(rect[i] - center) >= r2;
+	bits_corners >>= 1;
+
+	if (bits_corners == 0) {
+		// all inside
+		memcpy(result, rect, sizeof(rect));
+		memset(result_arc, true, 4);
+		return 4;
+	}
+
+	if (bits_corners == 0x0f)
+		// all outside
+		return 0;
+
+	// Some corners are outside, some are inside. Trim the rectangle.
+	int n = 0;
+	for (int i = 0; i < 4; ++ i) {
+		bool inside = (bits_corners & 0x08) == 0;
+		bits_corners <<= 1;
+		V2T chordal_points[2];
+		int n_chordal_points = line_circle_intersection(rect[i], rect[(i + 1)%4], center, radius, chordal_points);
+		if (n_chordal_points == 2) {
+			result_arc[n] = true;
+			result[n ++] = chordal_points[0];
+			result_arc[n] = true;
+			result[n ++] = chordal_points[1];
+		} else {
+			if (inside) {
+				result_arc[n] = false;
+				result[n ++] = rect[i];
+			}
+			if (n_chordal_points == 1) {
+				result_arc[n] = false;
+				result[n ++] = chordal_points[0];
+			}
+		}
+	}
+	return n;
+}
+/*
+// Calculate area of the circle x AABB intersection.
+// The calculation is approximate in a way, that the circular segment
+// intersecting the cell is approximated by its chord (a linear segment).
+template<typename T>
+T circle_AABB_intersection_area(
+	const boost::geometry::model::d2::point_xy<T>							    &center,
+	const T 																	 radius,
+	const boost::geometry::model::box<boost::geometry::model::d2::point_xy<T> > &aabb)
+{
+	typedef typename V2<T>::Type V2T;
+	typedef typename boost::geometry::model::box<V2T> B2T;
+	T radius2 = radius * radius;
+
+	bool intersectionLeft   = sqr(aabb.min_corner().x() - center.x()) < radius2;
+	bool intersectionRight  = sqr(aabb.max_corner().x() - center.x()) < radius2;
+	bool intersectionBottom = sqr(aabb.min_corner().y() - center.y()) < radius2;
+	bool intersectionTop    = sqr(aabb.max_corner().y() - center.y()) < radius2;
+
+	if (! (intersectionLeft || intersectionRight || intersectionTop || intersectionBottom))
+		// No intersection between the aabb and the center.
+		return boost::geometry::point_in_box<V2T, B2T>()::apply(center, aabb) ? 1.f : 0.f;
+
+
+
+	V2T rect[4] = {
+		aabb.min_corner(),
+		V2T(aabb.max_corner().x(), aabb.min_corner().y()),
+		aabb.max_corner(),
+		V2T(aabb.min_corner().x(), aabb.max_corner().y())
+	};
+
+	int  bits_corners = 0;
+	T    r2 = sqr(radius);
+	for (int i = 0; i < 4; ++ i, bits_corners <<= 1)
+		bits_corners |= dot(rect[i] - center) >= r2;
+	bits_corners >>= 1;
+
+	if (bits_corners == 0) {
+		// all inside
+		memcpy(result, rect, sizeof(rect));
+		memset(result_arc, true, 4);
+		return 4;
+	}
+
+	if (bits_corners == 0x0f)
+		// all outside
+		return 0;
+
+	// Some corners are outside, some are inside. Trim the rectangle.
+	int n = 0;
+	for (int i = 0; i < 4; ++ i) {
+		bool inside = (bits_corners & 0x08) == 0;
+		bits_corners <<= 1;
+		V2T chordal_points[2];
+		int n_chordal_points = line_circle_intersection(rect[i], rect[(i + 1)%4], center, radius, chordal_points);
+		if (n_chordal_points == 2) {
+			result_arc[n] = true;
+			result[n ++] = chordal_points[0];
+			result_arc[n] = true;
+			result[n ++] = chordal_points[1];
+		} else {
+			if (inside) {
+				result_arc[n] = false;
+				result[n ++] = rect[i];
+			}
+			if (n_chordal_points == 1) {
+				result_arc[n] = false;
+				result[n ++] = chordal_points[0];
+			}
+		}
+	}
+	return n;
+}
+*/
+
+template<typename T>
+inline T polyArea(const boost::geometry::model::d2::point_xy<T> *poly, int n)
+{
+	T area = T(0);
+	for (int i = 1; i + 1 < n; ++i)
+		area += cross(poly[i] - poly[0], poly[i + 1] - poly[0]);
+	return T(0.5) * area;
+}
+
+template<typename T>
+boost::geometry::model::d2::point_xy<T> polyCentroid(const boost::geometry::model::d2::point_xy<T> *poly, int n)
+{
+	boost::geometry::model::d2::point_xy<T> centroid(T(0), T(0));
+	for (int i = 0; i < n; ++i)
+		centroid += poly[i];
+	return (n == 0) ? centroid : (centroid / float(n));
+}
+
+void gcode_paint_layer(
+	const std::vector<V2f> 	&polyline,
+	float					 width,
+	float					 thickness,
+	A2f 					&acc)
+{
+	int nc = acc.shape()[1];
+	int nr = acc.shape()[0];
+//	printf("gcode_paint_layer %d,%d\n", nc, nr);
+	for (size_t iLine = 1; iLine != polyline.size(); ++iLine) {
+		const V2f &p1 = polyline[iLine - 1];
+		const V2f &p2 = polyline[iLine];
+		// printf("p1, p2:  %f,%f %f,%f\n", p1.x(), p1.y(), p2.x(), p2.y());
+		const V2f  dir = p2 - p1;
+		V2f vperp(- dir.y(), dir.x());
+		vperp = vperp * 0.5f * width / l2(vperp);
+		// Rectangle of the extrusion.
+		V2f rect[4] = { p1 + vperp, p1 - vperp, p2 - vperp, p2 + vperp };
+		// Bounding box of the extrusion.
+		B2f bboxLine(rect[0], rect[0]);
+		boost::geometry::expand(bboxLine, rect[1]);
+		boost::geometry::expand(bboxLine, rect[2]);
+		boost::geometry::expand(bboxLine, rect[3]);
+		B2i bboxLinei(
+			V2i(clamp(0, nc-1, int(floor(bboxLine.min_corner().x()))),
+				clamp(0, nr-1, int(floor(bboxLine.min_corner().y())))),
+			V2i(clamp(0, nc-1, int(ceil (bboxLine.max_corner().x()))),
+				clamp(0, nr-1, int(ceil (bboxLine.max_corner().y())))));
+		// printf("bboxLinei %d,%d %d,%d\n", bboxLinei.min_corner().x(), bboxLinei.min_corner().y(), bboxLinei.max_corner().x(), bboxLinei.max_corner().y());
+#ifdef _DEBUG
+		float area = polyArea(rect, 4);
+		assert(area > 0.f);
+#endif /* _DEBUG */
+		for (int j = bboxLinei.min_corner().y(); j + 1 < bboxLinei.max_corner().y(); ++ j) {
+			for (int i = bboxLinei.min_corner().x(); i + 1 < bboxLinei.max_corner().x(); ++i) {
+				V2f rect2[8];
+				memcpy(rect2, rect, sizeof(rect));
+				int n = clip_rect_by_AABB(rect2, B2f(V2f(float(i), float(j)), V2f(float(i + 1), float(j + 1))));
+				float area = polyArea(rect2, n);
+				assert(area >= 0.f && area <= 1.000001f);
+				acc[j][i] += area * thickness;
+			}
+		}
+	}
+}
+
+void gcode_paint_bitmap(
+	const std::vector<V2f> 	&polyline,
+	float					 width,
+	A2uc 					&bitmap,
+	float					 scale)
+{
+	int nc = bitmap.shape()[1];
+	int nr = bitmap.shape()[0];
+	float r2 = width * width * 0.25f;
+//	printf("gcode_paint_layer %d,%d\n", nc, nr);
+	for (size_t iLine = 1; iLine != polyline.size(); ++iLine) {
+		const V2f &p1 = polyline[iLine - 1];
+		const V2f &p2 = polyline[iLine];
+		// printf("p1, p2:  %f,%f %f,%f\n", p1.x(), p1.y(), p2.x(), p2.y());
+		V2f dir = p2 - p1;
+		dir = dir * 0.5f * width / l2(dir);
+		V2f vperp(- dir.y(), dir.x());
+		// Rectangle of the extrusion.
+		V2f rect[4] = { (p1 + vperp - dir) * scale, (p1 - vperp - dir) * scale, (p2 - vperp + dir) * scale, (p2 + vperp + dir) * scale };
+		// Bounding box of the extrusion.
+		B2f bboxLine(rect[0], rect[0]);
+		boost::geometry::expand(bboxLine, rect[1]);
+		boost::geometry::expand(bboxLine, rect[2]);
+		boost::geometry::expand(bboxLine, rect[3]);
+		B2i bboxLinei(
+			V2i(clamp(0, nc-1, int(floor(bboxLine.min_corner().x()))),
+				clamp(0, nr-1, int(floor(bboxLine.min_corner().y())))),
+			V2i(clamp(0, nc-1, int(ceil (bboxLine.max_corner().x()))),
+				clamp(0, nr-1, int(ceil (bboxLine.max_corner().y())))));
+		// printf("bboxLinei %d,%d %d,%d\n", bboxLinei.min_corner().x(), bboxLinei.min_corner().y(), bboxLinei.max_corner().x(), bboxLinei.max_corner().y());
+		for (int j = bboxLinei.min_corner().y(); j + 1 < bboxLinei.max_corner().y(); ++ j) {
+			for (int i = bboxLinei.min_corner().x(); i + 1 < bboxLinei.max_corner().x(); ++i) {
+				float d2 = dist2_to_line(p1, p2, V2f(float(i) + 0.5f, float(j) + 0.5f) / scale);
+				if (d2 < r2)
+					bitmap[j][i] = 1;
+			}
+		}
+	}
+}
+
+struct Cell
+{
+	// Cell index in the grid.
+	V2i   idx;
+	// Total volume of the material stored in this cell.
+	float volume;
+	// Area covered inside this cell, <0,1>.
+	float area;
+	// Fraction of the area covered by the print head. <0,1>
+	float fraction_covered;
+	// Height of the covered part in excess to the expected layer height.
+	float excess_height;
+
+	bool operator<(const Cell &c2) {
+		return this->excess_height < c2.excess_height;
+	}
+};
+
+struct ExtrusionPoint {
+	V2f   center;
+	float radius;
+	float height;
+};
+
+typedef std::vector<ExtrusionPoint> ExtrusionPoints;
+
+void gcode_spread_points(
+	A2f 					&acc,
+	const A2f				&mask,
+	const ExtrusionPoints   &points, 
+	Slic3r::ExtrusionSimulationType simulationType)
+{
+	int nc = acc.shape()[1];
+	int nr = acc.shape()[0];
+
+	// Maximum radius of the spreading points, to allocate a large enough cell array.
+	float rmax = 0.f;
+	for (ExtrusionPoints::const_iterator it = points.begin(); it != points.end(); ++ it)
+		rmax = std::max(rmax, it->radius);
+	size_t n_rows_max  = size_t(ceil(rmax * 2.f + 2.f));
+	size_t n_cells_max = sqr(n_rows_max);
+	std::vector<std::pair<float, float> > spans;
+	std::vector<Cell>  cells(n_cells_max, Cell());
+	std::vector<float> areas_sum(n_cells_max, 0.f);
+
+	for (ExtrusionPoints::const_iterator it = points.begin(); it != points.end(); ++ it) {
+		const V2f  &center = it->center;
+		const float radius = it->radius;
+		const float radius2 = radius * radius;
+		const float height_target = it->height;
+		B2f bbox(center - V2f(radius, radius), center + V2f(radius, radius));
+		B2i bboxi(
+			V2i(clamp(0, nc-1, int(floor(bbox.min_corner().x()))),
+				clamp(0, nr-1, int(floor(bbox.min_corner().y())))),
+			V2i(clamp(0, nc-1, int(ceil (bbox.max_corner().x()))),
+				clamp(0, nr-1, int(ceil (bbox.max_corner().y())))));
+		/*
+		// Fill in the spans, at which the circle intersects the rows.
+		int row_first = bboxi.min_corner().y();
+		int row_last  = bboxi.max_corner().y();
+		for (; row_first <= row_last; ++ row_first) {
+			float y     = float(j) - center.y();
+			float discr = radius2 - sqr(y);
+			if (discr > 0) {
+				// Circle intersects the row j at 2 points.
+				float d = sqrt(discr);
+				spans.push_back(std.pair<float, float>(center.x() - d, center.x() + d)));
+				break;
+			}
+		}
+		for (int j = row_first + 1; j <= row_last; ++ j) {
+			float y     = float(j) - center.y();
+			float discr = radius2 - sqr(y);
+			if (discr > 0) {
+				// Circle intersects the row j at 2 points.
+				float d = sqrt(discr);
+				spans.push_back(std.pair<float, float>(center.x() - d, center.x() + d)));
+			} else {
+				row_last = j - 1;
+				break;
+			}
+		}
+		*/
+		float area_total     = 0;
+		float volume_total   = 0;
+		float volume_excess  = 0;
+		float volume_deficit = 0;
+		size_t n_cells = 0;
+		float area_circle_total = 0; 
+#if 0
+		// The intermediate lines.
+		for (int j = row_first; j < row_last; ++ j) {
+			const std::pair<float, float> &span1 = spans[j];
+			const std::pair<float, float> &span2 = spans[j+1];
+			float l1 = span1.first;
+			float l2 = span2.first;
+			float r1 = span1.second;
+			float r2 = span2.second;
+			if (l2 < l1)
+				std::swap(l1, l2);
+			if (r1 > r2)
+				std::swap(r1, r2);
+			int il1 = int(floor(l1));
+			int il2 = int(ceil(l2));
+			int ir1 = int(floor(r1));
+			int ir2 = int(floor(r2));
+			assert(il2 <= ir1);
+			for (int i = il1; i < il2; ++ i) {
+				Cell &cell = cells[n_cells ++];
+				cell.idx.x(i);
+				cell.idx.y(j);
+				cell.area = area;
+			}
+			for (int i = il2; i < ir1; ++ i) {
+				Cell &cell = cells[n_cells ++];
+				cell.idx.x(i);
+				cell.idx.y(j);
+				cell.area = 1.f;
+			}
+			for (int i = ir1; i < ir2; ++ i) {
+				Cell &cell = cells[n_cells ++];
+				cell.idx.x(i);
+				cell.idx.y(j);
+				cell.area = area;
+			}
+		}
+#else
+		for (int j = bboxi.min_corner().y(); j < bboxi.max_corner().y(); ++ j) {
+			for (int i = bboxi.min_corner().x(); i < bboxi.max_corner().x(); ++i) {
+				B2f bb(V2f(float(i), float(j)), V2f(float(i + 1), float(j + 1)));
+				V2f poly[8];
+				bool poly_arc[8];
+				int n = clip_circle_by_AABB(center, radius, bb, poly, poly_arc);
+				float area = polyArea(poly, n);
+				assert(area >= 0.f && area <= 1.000001f);
+				if (area == 0.f)
+					continue;
+				Cell &cell = cells[n_cells ++];
+				cell.idx.x(i);
+				cell.idx.y(j);
+				cell.volume  = acc[j][i];
+				cell.area    = mask[j][i];
+				assert(cell.area >= 0.f && cell.area <= 1.000001f);
+				area_circle_total += area;
+				if (cell.area < area)
+					cell.area = area;
+				cell.fraction_covered = clamp(0.f, 1.f, (cell.area > 0) ? (area / cell.area) : 0);
+				if (cell.fraction_covered == 0) {
+					-- n_cells;
+					continue;
+				}
+				float cell_height = cell.volume / cell.area;
+				cell.excess_height = cell_height - height_target;
+				if (cell.excess_height > 0.f)
+					volume_excess  += cell.excess_height * cell.area * cell.fraction_covered;
+				else
+					volume_deficit -= cell.excess_height * cell.area * cell.fraction_covered;
+				volume_total += cell.volume * cell.fraction_covered;
+				area_total   += cell.area * cell.fraction_covered;
+			}
+		}
+#endif
+		float area_circle_total2 = float(M_PI) * sqr(radius);
+		float area_err = fabs(area_circle_total2 - area_circle_total) / area_circle_total2;
+		printf("area_circle_total: %f, %f, %f\n", area_circle_total, area_circle_total2, area_err);
+		float volume_full = float(M_PI) * sqr(radius) * height_target;
+//		if (true) {
+//		printf("volume_total: %f, volume_full: %f, fill factor: %f\n", volume_total, volume_full, 100.f - 100.f * volume_total / volume_full);
+//		printf("volume_full: %f, volume_excess+deficit: %f, volume_excess: %f, volume_deficit: %f\n", volume_full, volume_excess+volume_deficit, volume_excess, volume_deficit);
+		if (simulationType == Slic3r::ExtrusionSimulationSpreadFull || volume_total <= volume_full) {
+			// The volume under the circle is spreaded fully.
+			float height_avg = volume_total / area_total;
+			for (size_t i = 0; i < n_cells; ++ i) {
+				const Cell &cell = cells[i];
+				acc[cell.idx.y()][cell.idx.x()] = (1.f - cell.fraction_covered) * cell.volume + cell.fraction_covered * cell.area * height_avg;
+			}
+		} else if (simulationType == ExtrusionSimulationSpreadExcess) {
+			// The volume under the circle does not fit.
+			// 1) Fill the underfilled cells and remove them from the list.
+			float volume_borrowed_total = 0.;
+			for (size_t i = 0; i < n_cells;) {
+				Cell &cell = cells[i];
+				if (cell.excess_height <= 0) {
+					// Fill in the part of the cell below the circle.
+					float volume_borrowed = - cell.excess_height * cell.area * cell.fraction_covered;
+					assert(volume_borrowed >= 0.f);
+					acc[cell.idx.y()][cell.idx.x()] = cell.volume + volume_borrowed;
+					volume_borrowed_total += volume_borrowed;
+					cell = cells[-- n_cells];
+				} else
+					++ i;
+			}
+			// 2) Sort the remaining cells by their excess height.
+			std::sort(cells.begin(), cells.begin() + n_cells);
+			// 3) Prefix sum the areas per excess height.
+			// The excess height is discrete with the number of excess cells.
+			areas_sum[n_cells-1] = cells[n_cells-1].area * cells[n_cells-1].fraction_covered;
+			for (int i = n_cells - 2; i >= 0; -- i) {
+				const Cell &cell = cells[i];
+				areas_sum[i] = areas_sum[i + 1] + cell.area * cell.fraction_covered;
+			}
+			// 4) Find the excess height, where the volume_excess is over the volume_borrowed_total.
+			float volume_current = 0.f;
+			float excess_height_prev = 0.f;
+			size_t i_top = n_cells;
+			for (size_t i = 0; i < n_cells; ++ i) {
+				const Cell &cell = cells[i];
+				volume_current += (cell.excess_height - excess_height_prev) * areas_sum[i];
+				excess_height_prev = cell.excess_height;
+				if (volume_current > volume_borrowed_total) {
+					i_top = i;
+					break;
+				}
+			}
+			// 5) Remove material from the cells with deficit.
+			// First remove all the excess material from the cells, where the deficit is low.
+			for (size_t i = 0; i < i_top; ++ i) {
+				const Cell &cell = cells[i];
+				float volume_removed = cell.excess_height * cell.area * cell.fraction_covered;
+				acc[cell.idx.y()][cell.idx.x()] = cell.volume - volume_removed;
+				volume_borrowed_total -= volume_removed;
+			}
+			// Second remove some excess material from the cells, where the deficit is high.
+			if (i_top < n_cells) {
+				float height_diff = volume_borrowed_total / areas_sum[i_top];
+				for (size_t i = i_top; i < n_cells; ++ i) {
+					const Cell &cell = cells[i];
+					acc[cell.idx.y()][cell.idx.x()] = cell.volume - height_diff * cell.area * cell.fraction_covered;
+				}
+			}
+		}
+	}
+}
+
+inline std::vector<V3uc> CreatePowerColorGradient24bit()
+{
+	int i;
+	int iColor = 0;
+	std::vector<V3uc> out(6 * 255 + 1, V3uc(0, 0, 0));
+	for (i = 0; i < 256; ++i)
+		out[iColor++] = V3uc(0, 0, i);
+	for (i = 1; i < 256; ++i)
+		out[iColor++] = V3uc(0, i, 255);
+	for (i = 1; i < 256; ++i)
+		out[iColor++] = V3uc(0, 255, 256 - i);
+	for (i = 1; i < 256; ++i)
+		out[iColor++] = V3uc(i, 255, 0);
+	for (i = 1; i < 256; ++i)
+		out[iColor++] = V3uc(255, 256 - i, 0);
+	for (i = 1; i < 256; ++i)
+		out[iColor++] = V3uc(255, 0, i);
+	return out;
+}
+
+namespace Slic3r {
+
+class ExtrusionSimulatorImpl {
+public:
+	std::vector<unsigned char>  image_data;
+	A2f							accumulator;
+	A2uc						bitmap;
+	unsigned int 				bitmap_oversampled;
+	ExtrusionPoints 			extrusion_points;
+	// RGB gradient to color map the fullness of an accumulator bucket into the output image.
+	std::vector<boost::geometry::model::point<unsigned char, 3, boost::geometry::cs::cartesian> > color_gradient;
+};
+
+ExtrusionSimulator::ExtrusionSimulator() :
+	pimpl(new ExtrusionSimulatorImpl)
+{
+	pimpl->color_gradient = CreatePowerColorGradient24bit();
+	pimpl->bitmap_oversampled = 4;
+}
+
+ExtrusionSimulator::~ExtrusionSimulator()
+{
+	delete pimpl;
+	pimpl = NULL;
+}
+
+void ExtrusionSimulator::set_image_size(const Point &image_size)
+{
+	// printf("ExtrusionSimulator::set_image_size()\n");
+	if (this->image_size.x == image_size.x &&
+		this->image_size.y == image_size.y)
+		return;
+
+	// printf("Setting image size: %d, %d\n", image_size.x, image_size.y);
+	this->image_size = image_size;
+	// Allocate the image data in an RGBA format.
+	// printf("Allocating image data, size %d\n", image_size.x * image_size.y * 4);
+	pimpl->image_data.assign(image_size.x * image_size.y * 4, 0);
+	// printf("Allocating image data, allocated\n");
+
+	//FIXME fill the image with red vertical lines.
+	for (size_t r = 0; r < image_size.y; ++ r) {
+		for (size_t c = 0; c < image_size.x; c += 2) {
+			// Color red
+			pimpl->image_data[r * image_size.x * 4 + c * 4] = 255;
+			// Opacity full
+			pimpl->image_data[r * image_size.x * 4 + c * 4 + 3] = 255;
+		}
+	}
+	// printf("Allocating image data, set\n");
+}
+
+void ExtrusionSimulator::set_viewport(const BoundingBox &viewport)
+{
+	// printf("ExtrusionSimulator::set_viewport(%d, %d, %d, %d)\n", viewport.min.x, viewport.min.y, viewport.max.x, viewport.max.y);
+	if (this->viewport != viewport) {
+		this->viewport = viewport;
+		Point sz = viewport.size();
+		pimpl->accumulator.resize(boost::extents[sz.y][sz.x]);
+		pimpl->bitmap.resize(boost::extents[sz.y*pimpl->bitmap_oversampled][sz.x*pimpl->bitmap_oversampled]);
+		// printf("Accumulator size: %d, %d\n", sz.y, sz.x);
+	}
+}
+
+void ExtrusionSimulator::set_bounding_box(const BoundingBox &bbox)
+{
+	this->bbox = bbox;
+}
+
+const void* ExtrusionSimulator::image_ptr() const 
+{
+	return (pimpl->image_data.empty()) ? NULL : (void*)&pimpl->image_data.front();
+}
+
+void ExtrusionSimulator::reset_accumulator()
+{
+	// printf("ExtrusionSimulator::reset_accumulator()\n");
+	Point sz = viewport.size();
+	// printf("Reset accumulator, Accumulator size: %d, %d\n", sz.y, sz.x);
+	memset(&pimpl->accumulator[0][0], 0, sizeof(float) * sz.x * sz.y);
+	memset(&pimpl->bitmap[0][0], 0, sz.x * sz.y * pimpl->bitmap_oversampled * pimpl->bitmap_oversampled);
+	pimpl->extrusion_points.clear();
+	// printf("Reset accumulator, done.\n");
+}
+
+void ExtrusionSimulator::extrude_to_accumulator(const ExtrusionPath &path, const Point &shift, ExtrusionSimulationType simulationType)
+{
+	// printf("Extruding a path. Nr points: %d, width: %f, height: %f\r\n", path.polyline.points.size(), path.width, path.height);
+	// Convert the path to V2f points, shift and scale them to the viewport.
+	std::vector<V2f> polyline;
+	polyline.reserve(path.polyline.points.size());
+	float scalex  = float(viewport.size().x) / float(bbox.size().x);
+	float scaley  = float(viewport.size().y) / float(bbox.size().y);
+	float w = scale_(path.width) * scalex;
+	float h = scale_(path.height) * scalex;
+	w = scale_(path.mm3_per_mm / path.height) * scalex;
+	// printf("scalex: %f, scaley: %f\n", scalex, scaley);
+	// printf("bbox: %d,%d %d,%d\n", bbox.min.x, bbox.min.y, bbox.max.x, bbox.max.y);
+	for (Points::const_iterator it = path.polyline.points.begin(); it != path.polyline.points.end(); ++ it) {
+		// printf("point %d,%d\n", it->x+shift.x, it->y+shift.y);
+		ExtrusionPoint ept;
+		ept.center = V2f(float(it->x+shift.x-bbox.min.x) * scalex, float(it->y+shift.y-bbox.min.y) * scaley);
+		ept.radius = w/2.f;
+		ept.height = 0.5f;
+		polyline.push_back(ept.center);
+		pimpl->extrusion_points.push_back(ept);
+	}
+	// Extrude the polyline into an accumulator.
+	// printf("width scaled: %f, height scaled: %f\n", w, h);
+	gcode_paint_layer(polyline, w, 0.5f, pimpl->accumulator);
+
+ 	if (simulationType > ExtrusionSimulationDontSpread)
+		gcode_paint_bitmap(polyline, w, pimpl->bitmap, pimpl->bitmap_oversampled);
+    // double path.mm3_per_mm;  // mm^3 of plastic per mm of linear head motion
+    // float path.width;
+    // float path.height;
+}
+
+void ExtrusionSimulator::evaluate_accumulator(ExtrusionSimulationType simulationType)
+{
+	// printf("ExtrusionSimulator::evaluate_accumulator()\n");
+	Point sz = viewport.size();
+
+	if (simulationType > ExtrusionSimulationDontSpread) {
+		// Average the cells of a bitmap into a lower resolution floating point mask.
+		A2f mask(boost::extents[sz.y][sz.x]);
+		for (int r = 0; r < sz.y; ++r) {
+			for (int c = 0; c < sz.x; ++c) {
+				float p = 0;
+				for (int j = 0; j < pimpl->bitmap_oversampled; ++ j) {
+					for (int i = 0; i < pimpl->bitmap_oversampled; ++ i) {
+						if (pimpl->bitmap[r * pimpl->bitmap_oversampled + j][c * pimpl->bitmap_oversampled + i])
+							p += 1.f;
+					}
+				}
+				p /= float(pimpl->bitmap_oversampled * pimpl->bitmap_oversampled * 2);
+				mask[r][c] = p;
+			}
+		}
+
+		// Spread the excess of the material.
+		gcode_spread_points(pimpl->accumulator, mask, pimpl->extrusion_points, simulationType);
+	}
+
+	// Color map the accumulator.
+	for (int r = 0; r < sz.y; ++r) {
+		unsigned char *ptr = &pimpl->image_data[(image_size.x * (viewport.min.y + r) + viewport.min.x) * 4];
+		for (int c = 0; c < sz.x; ++c) {
+			#if 1
+			float p   = pimpl->accumulator[r][c];
+			#else
+			float p = mask[r][c];
+			#endif
+			int   idx = int(floor(p * float(pimpl->color_gradient.size()) + 0.5f));
+			V3uc  clr = pimpl->color_gradient[clamp(0, int(pimpl->color_gradient.size()-1), idx)];
+			*ptr ++ = clr.get<0>();
+			*ptr ++ = clr.get<1>();
+			*ptr ++ = clr.get<2>();
+			*ptr ++ = (idx == 0) ? 0 : 255;
+		}
+	}
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/ExtrusionSimulator.hpp b/xs/src/libslic3r/ExtrusionSimulator.hpp
new file mode 100644
index 0000000000..f64a3b7d31
--- /dev/null
+++ b/xs/src/libslic3r/ExtrusionSimulator.hpp
@@ -0,0 +1,60 @@
+#ifndef slic3r_ExtrusionSimulator_hpp_
+#define slic3r_ExtrusionSimulator_hpp_
+
+#include "libslic3r.h"
+#include "Config.hpp"
+#include "ExtrusionEntity.hpp"
+#include "BoundingBox.hpp"
+
+namespace Slic3r {
+
+enum ExtrusionSimulationType
+{
+    ExtrusionSimulationSimple,
+    ExtrusionSimulationDontSpread,
+    ExtrisopmSimulationSpreadNotOverfilled,
+    ExtrusionSimulationSpreadFull,
+    ExtrusionSimulationSpreadExcess,
+};
+
+// An opaque class, to keep the boost stuff away from the header.
+class ExtrusionSimulatorImpl;
+
+class ExtrusionSimulator
+{
+public:
+    ExtrusionSimulator();
+    ~ExtrusionSimulator();
+
+    // Size of the image, that will be returned by image_ptr().
+    // The image may be bigger than the viewport as many graphics drivers 
+    // expect the size of a texture to be rounded to a power of two.
+    void  		set_image_size(const Point &image_size);
+    // Which part of the image shall be rendered to?
+    void  		set_viewport(const BoundingBox &viewport);
+    // Shift and scale of the rendered extrusion paths into the viewport.
+    void		set_bounding_box(const BoundingBox &bbox);
+
+    // Reset the extrusion accumulator to zero for all buckets.
+    void		reset_accumulator();
+    // Paint a thick path into an extrusion buffer.
+    // A simple implementation is provided now, splatting a rectangular extrusion for each linear segment.
+    // In the future, spreading and suqashing of a material will be simulated.
+    void		extrude_to_accumulator(const ExtrusionPath &path, const Point &shift, ExtrusionSimulationType simulationType);
+    // Evaluate the content of the accumulator and paint it into the viewport.
+    // After this call the image_ptr() call will return a valid image.
+    void		evaluate_accumulator(ExtrusionSimulationType simulationType);
+    // An RGBA image of image_size, to be loaded into a GPU texture.
+    const void* image_ptr() const;
+
+private:
+	Point					    	image_size;
+	BoundingBox				    	viewport;
+	BoundingBox 					bbox;
+
+	ExtrusionSimulatorImpl		   *pimpl;
+};
+
+}
+
+#endif /* slic3r_ExtrusionSimulator_hpp_ */
diff --git a/xs/src/libslic3r/Fill/Fill3DHoneycomb.cpp b/xs/src/libslic3r/Fill/Fill3DHoneycomb.cpp
new file mode 100644
index 0000000000..ce3fa3c575
--- /dev/null
+++ b/xs/src/libslic3r/Fill/Fill3DHoneycomb.cpp
@@ -0,0 +1,223 @@
+#include "../ClipperUtils.hpp"
+#include "../PolylineCollection.hpp"
+#include "../Surface.hpp"
+
+#include "Fill3DHoneycomb.hpp"
+
+namespace Slic3r {
+
+/*
+Creates a contiguous sequence of points at a specified height that make
+up a horizontal slice of the edges of a space filling truncated
+octahedron tesselation. The octahedrons are oriented so that the
+square faces are in the horizontal plane with edges parallel to the X
+and Y axes.
+
+Credits: David Eccles (gringer).
+*/
+
+// Generate an array of points that are in the same direction as the
+// basic printing line (i.e. Y points for columns, X points for rows)
+// Note: a negative offset only causes a change in the perpendicular
+// direction
+static std::vector<coordf_t> colinearPoints(const coordf_t offset, const size_t baseLocation, size_t gridLength)
+{
+    const coordf_t offset2 = std::abs(offset / coordf_t(2.));
+    std::vector<coordf_t> points;
+    points.push_back(baseLocation - offset2);
+    for (size_t i = 0; i < gridLength; ++i) {
+        points.push_back(baseLocation + i + offset2);
+        points.push_back(baseLocation + i + 1 - offset2);
+    }
+    points.push_back(baseLocation + gridLength + offset2);
+    return points;
+}
+
+// Generate an array of points for the dimension that is perpendicular to
+// the basic printing line (i.e. X points for columns, Y points for rows)
+static std::vector<coordf_t> perpendPoints(const coordf_t offset, const size_t baseLocation, size_t gridLength)
+{
+    coordf_t offset2 = offset / coordf_t(2.);
+    coord_t  side    = 2 * (baseLocation & 1) - 1;
+    std::vector<coordf_t> points;
+    points.push_back(baseLocation - offset2 * side);
+    for (size_t i = 0; i < gridLength; ++i) {
+        side = 2*((i+baseLocation) & 1) - 1;
+        points.push_back(baseLocation + offset2 * side);
+        points.push_back(baseLocation + offset2 * side);
+    }
+    points.push_back(baseLocation - offset2 * side);
+    return points;
+}
+
+template<typename T>
+static inline T clamp(T low, T high, T x)
+{
+    return std::max<T>(low, std::min<T>(high, x));
+}
+
+// Trims an array of points to specified rectangular limits. Point
+// components that are outside these limits are set to the limits.
+static inline void trim(Pointfs &pts, coordf_t minX, coordf_t minY, coordf_t maxX, coordf_t maxY)
+{
+    for (Pointfs::iterator it = pts.begin(); it != pts.end(); ++ it) {
+        it->x = clamp(minX, maxX, it->x);
+        it->y = clamp(minY, maxY, it->y);
+    }
+}
+
+static inline Pointfs zip(const std::vector<coordf_t> &x, const std::vector<coordf_t> &y)
+{
+    assert(x.size() == y.size());
+    Pointfs out;
+    out.reserve(x.size());
+    for (size_t i = 0; i < x.size(); ++ i)
+        out.push_back(Pointf(x[i], y[i]));
+    return out;
+}
+
+// Generate a set of curves (array of array of 2d points) that describe a
+// horizontal slice of a truncated regular octahedron with edge length 1.
+// curveType specifies which lines to print, 1 for vertical lines
+// (columns), 2 for horizontal lines (rows), and 3 for both.
+static std::vector<Pointfs> makeNormalisedGrid(coordf_t z, size_t gridWidth, size_t gridHeight, size_t curveType)
+{
+    // offset required to create a regular octagram
+    coordf_t octagramGap = coordf_t(0.5);
+    
+    // sawtooth wave function for range f($z) = [-$octagramGap .. $octagramGap]
+    coordf_t a = std::sqrt(coordf_t(2.));  // period
+    coordf_t wave = abs(fmod(z, a) - a/2.)/a*4. - 1.;
+    coordf_t offset = wave * octagramGap;
+    
+    std::vector<Pointfs> points;
+    if ((curveType & 1) != 0) {
+        for (size_t x = 0; x <= gridWidth; ++x) {
+            points.push_back(Pointfs());
+            Pointfs &newPoints = points.back();
+            newPoints = zip(
+                perpendPoints(offset, x, gridHeight), 
+                colinearPoints(offset, 0, gridHeight));
+            // trim points to grid edges
+            trim(newPoints, coordf_t(0.), coordf_t(0.), coordf_t(gridWidth), coordf_t(gridHeight));
+            if (x & 1)
+                std::reverse(newPoints.begin(), newPoints.end());
+        }
+    }
+    if ((curveType & 2) != 0) {
+        for (size_t y = 0; y <= gridHeight; ++y) {
+            points.push_back(Pointfs());
+            Pointfs &newPoints = points.back();
+            newPoints = zip(
+                colinearPoints(offset, 0, gridWidth),
+                perpendPoints(offset, y, gridWidth));
+            // trim points to grid edges
+            trim(newPoints, coordf_t(0.), coordf_t(0.), coordf_t(gridWidth), coordf_t(gridHeight));
+            if (y & 1)
+                std::reverse(newPoints.begin(), newPoints.end());
+        }
+    }
+    return points;
+}
+
+// Generate a set of curves (array of array of 2d points) that describe a
+// horizontal slice of a truncated regular octahedron with a specified
+// grid square size.
+static Polylines makeGrid(coord_t z, coord_t gridSize, size_t gridWidth, size_t gridHeight, size_t curveType)
+{
+    coord_t  scaleFactor = gridSize;
+    coordf_t normalisedZ = coordf_t(z) / coordf_t(scaleFactor);
+    std::vector<Pointfs> polylines = makeNormalisedGrid(normalisedZ, gridWidth, gridHeight, curveType);
+    Polylines result;
+    result.reserve(polylines.size());
+    for (std::vector<Pointfs>::const_iterator it_polylines = polylines.begin(); it_polylines != polylines.end(); ++ it_polylines) {
+        result.push_back(Polyline());
+        Polyline &polyline = result.back();
+        for (Pointfs::const_iterator it = it_polylines->begin(); it != it_polylines->end(); ++ it)
+            polyline.points.push_back(Point(coord_t(it->x * scaleFactor), coord_t(it->y * scaleFactor)));
+    }
+    return result;
+}
+
+Polylines Fill3DHoneycomb::fill_surface(const Surface *surface, const FillParams &params)
+{
+    ExPolygon   expolygon = surface->expolygon;
+    BoundingBox bb = expolygon.contour.bounding_box();
+    Point       size = bb.size();
+    coord_t     distance = coord_t(scale_(this->spacing) / params.density);
+
+    // align bounding box to a multiple of our honeycomb grid module
+    // (a module is 2*$distance since one $distance half-module is 
+    // growing while the other $distance half-module is shrinking)
+    bb.merge(Point(
+        bb.min.x - (bb.min.x % (2*distance)),
+        bb.min.y - (bb.min.y % (2*distance))));
+    
+    // generate pattern
+    Polylines polylines = makeGrid(
+        scale_(this->z),
+        distance,
+        ceil(size.x / distance) + 1,
+        ceil(size.y / distance) + 1,
+        ((this->layer_id / surface->thickness_layers) % 2) + 1);
+    
+    // move pattern in place
+    for (Polylines::iterator it = polylines.begin(); it != polylines.end(); ++ it)
+        it->translate(bb.min.x, bb.min.y);
+
+    // clip pattern to boundaries
+    intersection(polylines, (Polygons)expolygon, &polylines);
+
+    // connect lines
+    if (! params.dont_connect && ! polylines.empty()) { // prevent calling leftmost_point() on empty collections
+        ExPolygon expolygon_off;
+        {
+            ExPolygons expolygons_off = offset_ex(expolygon, SCALED_EPSILON);
+            if (! expolygons_off.empty()) {
+                // When expanding a polygon, the number of islands could only shrink. Therefore the offset_ex shall generate exactly one expanded island for one input island.
+                assert(expolygons_off.size() == 1);
+                std::swap(expolygon_off, expolygons_off.front());
+            }
+        }
+        Polylines chained = PolylineCollection::chained_path_from(
+#if SLIC3R_CPPVER >= 11
+            std::move(polylines), 
+#else
+            polylines,
+#endif
+            PolylineCollection::leftmost_point(polylines), false); // reverse allowed
+#if SLIC3R_CPPVER >= 11
+            assert(polylines.empty());
+#else
+            polylines.clear();
+#endif
+        for (Polylines::iterator it_polyline = chained.begin(); it_polyline != chained.end(); ++ it_polyline) {
+            if (! polylines.empty()) {
+                // Try to connect the lines.
+                Points &pts_end = polylines.back().points;
+                const Point &first_point = it_polyline->points.front();
+                const Point &last_point = pts_end.back();
+                // TODO: we should also check that both points are on a fill_boundary to avoid 
+                // connecting paths on the boundaries of internal regions
+                if (first_point.distance_to(last_point) <= 1.5 * distance && 
+                    expolygon_off.contains(Line(last_point, first_point))) {
+                    // Append the polyline.
+                    pts_end.insert(pts_end.end(), it_polyline->points.begin(), it_polyline->points.end());
+                    continue;
+                }
+            }
+            // The lines cannot be connected.
+#if SLIC3R_CPPVER >= 11
+            polylines.push_back(std::move(*it_polyline));
+#else
+            polylines.push_back(Polyline());
+            std::swap(polylines.back(), *it_polyline);
+#endif
+        }
+    }
+    
+    // TODO: return ExtrusionLoop objects to get better chained paths
+    return polylines;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/Fill3DHoneycomb.hpp b/xs/src/libslic3r/Fill/Fill3DHoneycomb.hpp
new file mode 100644
index 0000000000..1410f5b55a
--- /dev/null
+++ b/xs/src/libslic3r/Fill/Fill3DHoneycomb.hpp
@@ -0,0 +1,24 @@
+#ifndef slic3r_Fill3DHoneycomb_hpp_
+#define slic3r_Fill3DHoneycomb_hpp_
+
+#include <map>
+
+#include "../libslic3r.h"
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+class Fill3DHoneycomb : public FillWithDirection
+{
+public:
+    virtual ~Fill3DHoneycomb() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+	
+	// require bridge flow since most of this pattern hangs in air
+    virtual bool use_bridge_flow() const { return true; }
+};
+
+} // namespace Slic3r
+
+#endif // slic3r_Fill3DHoneycomb_hpp_
diff --git a/xs/src/libslic3r/Fill/FillBase.cpp b/xs/src/libslic3r/Fill/FillBase.cpp
new file mode 100644
index 0000000000..c17f57185d
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillBase.cpp
@@ -0,0 +1,81 @@
+#include "../Surface.hpp"
+
+#include "FillBase.hpp"
+#include "FillConcentric.hpp"
+#include "FillHoneycomb.hpp"
+#include "Fill3DHoneycomb.hpp"
+#include "FillPlanePath.hpp"
+#include "FillRectilinear.hpp"
+#include "FillRectilinear2.hpp"
+
+namespace Slic3r {
+
+Fill* Fill::new_from_type(const std::string &type)
+{
+	if (type == "concentric")
+		return new FillConcentric();
+	if (type == "honeycomb")
+		return new FillHoneycomb();
+	if (type == "3dhoneycomb")
+		return new Fill3DHoneycomb();
+	if (type == "rectilinear")
+//		return new FillRectilinear();
+        return new FillRectilinear2();
+	if (type == "line")
+		return new FillLine();
+	if (type == "grid")
+		return new FillGrid();
+	if (type == "archimedeanchords")
+		return new FillArchimedeanChords();
+	if (type == "hilbertcurve")
+		return new FillHilbertCurve();
+	if (type == "octagramspiral")
+		return new FillOctagramSpiral();
+	CONFESS("unknown type");
+	return NULL;
+}
+
+coord_t Fill::adjust_solid_spacing(const coord_t width, const coord_t distance)
+{
+    coord_t number_of_lines = coord_t(coordf_t(width) / coordf_t(distance)) + 1;
+    coord_t extra_space     = width % distance;
+    return (number_of_lines <= 1) ? 
+    	distance : 
+    	distance + extra_space / (number_of_lines - 1);
+}
+
+std::pair<float, Point> FillWithDirection::infill_direction(const Surface *surface) const
+{
+    // set infill angle
+    float out_angle = this->angle;
+
+	if (out_angle == FLT_MAX) {
+		//FIXME Vojtech: Add a warning?
+        printf("Using undefined infill angle\n");
+        out_angle = 0.f;
+    }
+
+    Point out_shift = empty(this->bounding_box) ? 
+    	surface->expolygon.contour.bounding_box().center() : 
+        this->bounding_box.center();
+
+    if (surface->bridge_angle >= 0) {
+	    // use bridge angle
+		//FIXME Vojtech: Add a debugf?
+        // Slic3r::debugf "Filling bridge with angle %d\n", rad2deg($surface->bridge_angle);
+#ifdef SLIC3R_DEBUG
+        printf("Filling bridge with angle %f\n", surface->bridge_angle);
+#endif /* SLIC3R_DEBUG */
+        out_angle = surface->bridge_angle;
+    } else if (this->layer_id != size_t(-1)) {
+        // alternate fill direction
+        out_angle += this->_layer_angle(this->layer_id / surface->thickness_layers);
+    } else {
+    	printf("Layer_ID undefined!\n");
+    }
+
+    out_angle += float(M_PI/2.);
+    return std::pair<float, Point>(out_angle, out_shift);
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillBase.hpp b/xs/src/libslic3r/Fill/FillBase.hpp
new file mode 100644
index 0000000000..8fa2bee423
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillBase.hpp
@@ -0,0 +1,103 @@
+#ifndef slic3r_FillBase_hpp_
+#define slic3r_FillBase_hpp_
+
+#include <float.h>
+
+#include "../libslic3r.h"
+#include "../BoundingBox.hpp"
+
+namespace Slic3r {
+
+class Surface;
+
+struct FillParams
+{
+    FillParams() { memset(this, 0, sizeof(FillParams)); }
+
+    coordf_t    width;
+    // Fraction in <0, 1>
+    float       density;
+    coordf_t    distance;
+
+    // Don't connect the fill lines around the inner perimeter.
+    bool        dont_connect;
+
+    // Don't adjust spacing to fill the space evenly.
+    bool        dont_adjust;
+
+    // For Honeycomb.
+    // we were requested to complete each loop;
+    // in this case we don't try to make more continuous paths
+    bool        complete;
+};
+
+class Fill
+{
+public:
+    // Index of the layer.
+    size_t      layer_id;
+    // Height of the layer, in unscaled coordinates
+    coordf_t    z;
+    // in unscaled coordinates
+    coordf_t    spacing;
+    // in radians, ccw, 0 = East
+    float       angle;
+    // in scaled coordinates
+    coord_t     loop_clipping;
+    // in scaled coordinates
+    BoundingBox bounding_box;
+
+public:
+    virtual ~Fill() {}
+
+    static Fill* new_from_type(const std::string &type);
+
+    void         set_bounding_box(const Slic3r::BoundingBox &bbox) { bounding_box = bbox; }
+
+    // Use bridge flow for the fill?
+    virtual bool use_bridge_flow() const { return false; }
+
+    // Do not sort the fill lines to optimize the print head path?
+    virtual bool no_sort() const { return false; }
+
+    // Perform the fill.
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params) = 0;
+
+protected:
+    Fill() :
+        layer_id(size_t(-1)),
+        z(0.f),
+        spacing(0.f),
+        // Initial angle is undefined.
+        angle(FLT_MAX),
+        loop_clipping(0),
+        // The initial bounding box is empty, therefore undefined.
+        bounding_box(Point(0, 0), Point(-1, -1))
+        {}
+
+    static coord_t adjust_solid_spacing(const coord_t width, const coord_t distance);
+};
+
+// An interface class to Perl, aggregating an instance of a Fill and a FillData.
+class Filler
+{
+public:
+    Filler() : fill(NULL) {}
+    ~Filler() { delete fill; fill = NULL; }
+    Fill        *fill;
+    FillParams   params;
+};
+
+class FillWithDirection : public Fill
+{
+public:
+    virtual float _layer_angle(size_t idx) const { 
+        bool odd = idx & 1;
+        return (idx & 1) ? float(M_PI/2.) : 0;
+    }
+    virtual std::pair<float, Point> infill_direction(const Surface *surface) const ;
+};
+
+} // namespace Slic3r
+
+#endif // slic3r_FillBase_hpp_
diff --git a/xs/src/libslic3r/Fill/FillConcentric.cpp b/xs/src/libslic3r/Fill/FillConcentric.cpp
new file mode 100644
index 0000000000..30b34456dd
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillConcentric.cpp
@@ -0,0 +1,60 @@
+#include "../ClipperUtils.hpp"
+#include "../ExPolygon.hpp"
+#include "../Surface.hpp"
+
+#include "FillConcentric.hpp"
+
+namespace Slic3r {
+
+Polylines FillConcentric::fill_surface(const Surface *surface, const FillParams &params)
+{
+    // no rotation is supported for this infill pattern
+    ExPolygon expolygon = surface->expolygon;
+    BoundingBox bounding_box = expolygon.contour.bounding_box();
+    
+    coord_t min_spacing = scale_(this->spacing);
+    coord_t distance = coord_t(min_spacing / params.density);
+    
+    if (params.density > 0.9999f && !params.dont_adjust) {
+        distance = this->adjust_solid_spacing(bounding_box.size().x, distance);
+        this->spacing = unscale(distance);
+    }
+
+    Polygons loops = (Polygons)expolygon;
+    Polygons last  = loops;
+    while (! last.empty()) {
+        last = offset2(last, -(distance + min_spacing/2), +min_spacing/2);
+        loops.insert(loops.end(), last.begin(), last.end());
+    }
+
+    // generate paths from the outermost to the innermost, to avoid
+    // adhesion problems of the first central tiny loops
+    union_pt_chained(loops, &loops, false);
+    
+    // split paths using a nearest neighbor search
+    Polylines paths;
+    Point last_pos(0, 0);
+    for (Polygons::const_iterator it_loop = loops.begin(); it_loop != loops.end(); ++ it_loop) {
+        paths.push_back(it_loop->split_at_index(last_pos.nearest_point_index(*it_loop)));
+        last_pos = paths.back().last_point();
+    }
+
+    // clip the paths to prevent the extruder from getting exactly on the first point of the loop
+    // Keep valid paths only.
+    size_t j = 0;
+    for (size_t i = 0; i < paths.size(); ++ i) {
+        paths[i].clip_end(this->loop_clipping);
+        if (paths[i].is_valid()) {
+            if (j < i)
+                std::swap(paths[j], paths[i]);
+            ++ j;
+        }
+    }
+    if (j < paths.size())
+        paths.erase(paths.begin() + j, paths.end());
+
+    // TODO: return ExtrusionLoop objects to get better chained paths
+    return paths;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillConcentric.hpp b/xs/src/libslic3r/Fill/FillConcentric.hpp
new file mode 100644
index 0000000000..6d58ec0894
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillConcentric.hpp
@@ -0,0 +1,20 @@
+#ifndef slic3r_FillConcentric_hpp_
+#define slic3r_FillConcentric_hpp_
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+class FillConcentric : public Fill
+{
+public:
+    virtual ~FillConcentric() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+
+protected:
+	virtual bool no_sort() const { return true; }
+};
+
+} // namespace Slic3r
+
+#endif // slic3r_FillConcentric_hpp_
diff --git a/xs/src/libslic3r/Fill/FillHoneycomb.cpp b/xs/src/libslic3r/Fill/FillHoneycomb.cpp
new file mode 100644
index 0000000000..408a5fa81b
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillHoneycomb.cpp
@@ -0,0 +1,133 @@
+#include "../ClipperUtils.hpp"
+#include "../PolylineCollection.hpp"
+#include "../Surface.hpp"
+
+#include "FillHoneycomb.hpp"
+
+namespace Slic3r {
+
+Polylines FillHoneycomb::fill_surface(const Surface *surface, const FillParams &params)
+{
+    std::pair<float, Point> rotate_vector = this->infill_direction(surface);
+    
+    // cache hexagons math
+    CacheID cache_id(params.density, this->spacing);
+    Cache::iterator it_m = this->cache.find(cache_id);
+    if (it_m == this->cache.end()) {
+#if SLIC3R_CPPVER > 11
+        it_m = this->cache.emplace_hint(it_m);
+#else
+        it_m = this->cache.insert(it_m, std::pair<CacheID, CacheData>(cache_id, CacheData()));
+#endif
+        CacheData &m = it_m->second;
+        coord_t min_spacing = scale_(this->spacing);
+        m.distance = min_spacing / params.density;
+        m.hex_side = m.distance / (sqrt(3)/2);
+        m.hex_width = m.distance * 2; // $m->{hex_width} == $m->{hex_side} * sqrt(3);
+        coord_t hex_height = m.hex_side * 2;
+        m.pattern_height = hex_height + m.hex_side;
+        m.y_short = m.distance * sqrt(3)/3;
+        m.x_offset = min_spacing / 2;
+        m.y_offset = m.x_offset * sqrt(3)/3;
+        m.hex_center = Point(m.hex_width/2, m.hex_side);
+    }
+    CacheData &m = it_m->second;
+
+    Polygons polygons;
+    {
+        // adjust actual bounding box to the nearest multiple of our hex pattern
+        // and align it so that it matches across layers
+        
+        BoundingBox bounding_box = surface->expolygon.contour.bounding_box();
+        {
+            // rotate bounding box according to infill direction
+            Polygon bb_polygon = bounding_box.polygon();
+            bb_polygon.rotate(rotate_vector.first, m.hex_center);
+            bounding_box = bb_polygon.bounding_box();
+            
+            // extend bounding box so that our pattern will be aligned with other layers
+            // $bounding_box->[X1] and [Y1] represent the displacement between new bounding box offset and old one
+            bounding_box.merge(Point(
+                bounding_box.min.x - (bounding_box.min.x % m.hex_width),
+                bounding_box.min.y - (bounding_box.min.y % m.pattern_height)));
+        }
+        
+        coord_t x = bounding_box.min.x;
+        while (x <= bounding_box.max.x) {
+            Polygon p;
+            coord_t ax[2] = { x + m.x_offset, x + m.distance - m.x_offset };
+            for (size_t i = 0; i < 2; ++ i) {
+                std::reverse(p.points.begin(), p.points.end()); // turn first half upside down
+                for (coord_t y = bounding_box.min.y; y <= bounding_box.max.y; y += m.y_short + m.hex_side + m.y_short + m.hex_side) {
+                    p.points.push_back(Point(ax[1], y + m.y_offset));
+                    p.points.push_back(Point(ax[0], y + m.y_short - m.y_offset));
+                    p.points.push_back(Point(ax[0], y + m.y_short + m.hex_side + m.y_offset));
+                    p.points.push_back(Point(ax[1], y + m.y_short + m.hex_side + m.y_short - m.y_offset));
+                    p.points.push_back(Point(ax[1], y + m.y_short + m.hex_side + m.y_short + m.hex_side + m.y_offset));
+                }
+                ax[0] = ax[0] + m.distance;
+                ax[1] = ax[1] + m.distance;
+                std::swap(ax[0], ax[1]); // draw symmetrical pattern
+                x += m.distance;
+            }
+            p.rotate(-rotate_vector.first, m.hex_center);
+            polygons.push_back(p);
+        }
+    }
+    
+    Polylines paths;
+    if (params.complete || true) {
+        // we were requested to complete each loop;
+        // in this case we don't try to make more continuous paths
+        Polygons polygons_trimmed = intersection((Polygons)*surface, polygons);
+        for (Polygons::iterator it = polygons_trimmed.begin(); it != polygons_trimmed.end(); ++ it)
+            paths.push_back(it->split_at_first_point());
+    } else {
+        // consider polygons as polylines without re-appending the initial point:
+        // this cuts the last segment on purpose, so that the jump to the next 
+        // path is more straight
+        {
+            Polylines p;
+            for (Polygons::iterator it = polygons.begin(); it != polygons.end(); ++ it)
+                p.push_back((Polyline)(*it));
+            paths = intersection(p, (Polygons)*surface);
+        }
+
+        // connect paths
+        if (! paths.empty()) { // prevent calling leftmost_point() on empty collections
+            Polylines chained = PolylineCollection::chained_path_from(
+#if SLIC3R_CPPVER >= 11
+                std::move(paths), 
+#else
+                paths,
+#endif
+                PolylineCollection::leftmost_point(paths), false);
+            assert(paths.empty());
+            paths.clear();
+            for (Polylines::iterator it_path = chained.begin(); it_path != chained.end(); ++ it_path) {
+                if (! paths.empty()) {
+                    // distance between first point of this path and last point of last path
+                    double distance = paths.back().last_point().distance_to(it_path->first_point());
+                    if (distance <= m.hex_width) {
+                        paths.back().points.insert(paths.back().points.end(), it_path->points.begin(), it_path->points.end());
+                        continue;
+                    }
+                }
+                // Don't connect the paths.
+                paths.push_back(*it_path);
+            }
+        }
+        
+        // clip paths again to prevent connection segments from crossing the expolygon boundaries
+        Polylines paths_trimmed = intersection(paths, to_polygons(offset_ex(surface->expolygon, SCALED_EPSILON)));
+#if SLIC3R_CPPVER >= 11
+        paths = std::move(paths_trimmed);
+#else
+        std::swap(paths, paths_trimmed);
+#endif
+    }
+    
+    return paths;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillHoneycomb.hpp b/xs/src/libslic3r/Fill/FillHoneycomb.hpp
new file mode 100644
index 0000000000..63fe1f6f4c
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillHoneycomb.hpp
@@ -0,0 +1,50 @@
+#ifndef slic3r_FillHoneycomb_hpp_
+#define slic3r_FillHoneycomb_hpp_
+
+#include <map>
+
+#include "../libslic3r.h"
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+class FillHoneycomb : public FillWithDirection
+{
+public:
+    virtual ~FillHoneycomb() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+
+protected:
+	// Caching the 
+	struct CacheID 
+	{
+		CacheID(float adensity, coordf_t aspacing) : 
+			density(adensity), spacing(aspacing) {}
+		float		density;
+		coordf_t	spacing;
+		bool operator<(const CacheID &other) const 
+			{ return (density < other.density) || (density == other.density && spacing < other.spacing); }
+		bool operator==(const CacheID &other) const 
+			{ return density == other.density && spacing == other.spacing; }
+	};
+	struct CacheData
+	{
+		coord_t	distance;
+        coord_t hex_side;
+        coord_t hex_width;
+        coord_t	pattern_height;
+        coord_t y_short;
+        coord_t x_offset;
+        coord_t	y_offset;
+        Point	hex_center;
+    };
+    typedef std::map<CacheID, CacheData> Cache;
+	Cache cache;
+
+    virtual float _layer_angle(size_t idx) const { return 0.5f * float(M_PI) * (idx % 3); }
+};
+
+} // namespace Slic3r
+
+#endif // slic3r_FillHoneycomb_hpp_
diff --git a/xs/src/libslic3r/Fill/FillPlanePath.cpp b/xs/src/libslic3r/Fill/FillPlanePath.cpp
new file mode 100644
index 0000000000..0a271cbd8b
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillPlanePath.cpp
@@ -0,0 +1,197 @@
+#include "../ClipperUtils.hpp"
+#include "../PolylineCollection.hpp"
+#include "../Surface.hpp"
+
+#include "FillPlanePath.hpp"
+
+namespace Slic3r {
+
+Polylines FillPlanePath::fill_surface(const Surface *surface, const FillParams &params)
+{
+    ExPolygon expolygon = surface->expolygon;
+    std::pair<float, Point> rotate_vector = this->infill_direction(surface);
+    expolygon.rotate(- rotate_vector.first);
+
+    coord_t distance_between_lines = scale_(this->spacing) / params.density;
+    
+    // align infill across layers using the object's bounding box
+    Polygon bb_polygon = this->bounding_box.polygon();
+    bb_polygon.rotate(- rotate_vector.first);
+    BoundingBox bounding_box = bb_polygon.bounding_box();
+    
+    Point shift = this->_centered() ? 
+        bounding_box.center() :
+        bounding_box.min;
+    expolygon.translate(-shift.x, -shift.y);
+    bounding_box.translate(-shift.x, -shift.y);
+
+    Pointfs pts = _generate(
+        coord_t(ceil(coordf_t(bounding_box.min.x) / distance_between_lines)),
+        coord_t(ceil(coordf_t(bounding_box.min.y) / distance_between_lines)),
+        coord_t(ceil(coordf_t(bounding_box.max.x) / distance_between_lines)),
+        coord_t(ceil(coordf_t(bounding_box.max.y) / distance_between_lines)));
+
+    Polylines polylines;
+    if (pts.size() >= 2) {
+        // Convert points to a polyline, upscale.
+        polylines.push_back(Polyline());
+        Polyline &polyline = polylines.back();
+        polyline.points.reserve(pts.size());
+        for (Pointfs::iterator it = pts.begin(); it != pts.end(); ++ it)
+            polyline.points.push_back(Point(
+                coord_t(floor(it->x * distance_between_lines + 0.5)), 
+                coord_t(floor(it->y * distance_between_lines + 0.5))));
+//      intersection(polylines_src, offset((Polygons)expolygon, scale_(0.02)), &polylines);
+        intersection(polylines, (Polygons)expolygon, &polylines);
+
+/*        
+        if (1) {
+            require "Slic3r/SVG.pm";
+            print "Writing fill.svg\n";
+            Slic3r::SVG::output("fill.svg",
+                no_arrows       => 1,
+                polygons        => \@$expolygon,
+                green_polygons  => [ $bounding_box->polygon ],
+                polylines       => [ $polyline ],
+                red_polylines   => \@paths,
+            );
+        }
+*/
+        
+        // paths must be repositioned and rotated back
+        for (Polylines::iterator it = polylines.begin(); it != polylines.end(); ++ it) {
+            it->translate(shift.x, shift.y);
+            it->rotate(rotate_vector.first);
+        }
+    }
+    
+    return polylines;
+}
+
+// Follow an Archimedean spiral, in polar coordinates: r=a+b\theta
+Pointfs FillArchimedeanChords::_generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y)
+{
+    // Radius to achieve.
+    coordf_t rmax = std::sqrt(coordf_t(max_x)*coordf_t(max_x)+coordf_t(max_y)*coordf_t(max_y)) * std::sqrt(2.) + 1.5;
+    // Now unwind the spiral.
+    coordf_t a = 1.;
+    coordf_t b = 1./(2.*M_PI);
+    coordf_t theta = 0.;
+    coordf_t r = 1;
+    Pointfs out;
+    //FIXME Vojtech: If used as a solid infill, there is a gap left at the center.
+    out.push_back(Pointf(0, 0));
+    out.push_back(Pointf(1, 0));
+    while (r < rmax) {
+        theta += 1. / r;
+        r = a + b * theta;
+        out.push_back(Pointf(r * cos(theta), r * sin(theta)));
+    }
+    return out;
+}
+
+// Adapted from 
+// http://cpansearch.perl.org/src/KRYDE/Math-PlanePath-122/lib/Math/PlanePath/HilbertCurve.pm
+//
+// state=0    3--2   plain
+//               |
+//            0--1
+//
+// state=4    1--2  transpose
+//            |  |
+//            0  3
+//
+// state=8
+//
+// state=12   3  0  rot180 + transpose
+//            |  |
+//            2--1
+//
+static inline Point hilbert_n_to_xy(const size_t n)
+{
+    static const int next_state[16] = { 4,0,0,12, 0,4,4,8, 12,8,8,4, 8,12,12,0 };
+    static const int digit_to_x[16] = { 0,1,1,0, 0,0,1,1, 1,0,0,1, 1,1,0,0 };
+    static const int digit_to_y[16] = { 0,0,1,1, 0,1,1,0, 1,1,0,0, 1,0,0,1 };
+
+    // Number of 2 bit digits.
+    size_t ndigits = 0;
+    {
+        size_t nc = n;
+        while(nc > 0) {
+            nc >>= 2;
+            ++ ndigits;
+        }
+    }
+    int state    = (ndigits & 1) ? 4 : 0;
+    int dirstate = (ndigits & 1) ? 0 : 4;
+    coord_t x = 0;
+    coord_t y = 0;
+    for (int i = (int)ndigits - 1; i >= 0; -- i) {
+        int digit = (n >> (i * 2)) & 3;
+        state += digit;
+        if (digit != 3)
+            dirstate = state; // lowest non-3 digit
+        x |= digit_to_x[state] << i;
+        y |= digit_to_y[state] << i;
+        state = next_state[state];
+    }
+    return Point(x, y);
+}
+
+Pointfs FillHilbertCurve::_generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y)
+{
+    // Minimum power of two square to fit the domain.
+    size_t sz = 2;
+    size_t pw = 1;
+    {
+        size_t sz0 = std::max(max_x + 1 - min_x, max_y + 1 - min_y);
+        while (sz < sz0) {
+            sz = sz << 1;
+            ++ pw;
+        }
+    }
+
+    size_t sz2 = sz * sz;
+    Pointfs line;
+    line.reserve(sz2);
+    for (size_t i = 0; i < sz2; ++ i) {
+        Point p = hilbert_n_to_xy(i);
+        line.push_back(Pointf(p.x + min_x, p.y + min_y));
+    }
+    return line;
+}
+
+Pointfs FillOctagramSpiral::_generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y)
+{
+    // Radius to achieve.
+    coordf_t rmax = std::sqrt(coordf_t(max_x)*coordf_t(max_x)+coordf_t(max_y)*coordf_t(max_y)) * std::sqrt(2.) + 1.5;
+    // Now unwind the spiral.
+    coordf_t r = 0;
+    coordf_t r_inc = sqrt(2.);
+    Pointfs out;
+    out.push_back(Pointf(0, 0));
+    while (r < rmax) {
+        r += r_inc;
+        coordf_t rx = r / sqrt(2.);
+        coordf_t r2 = r + rx;
+        out.push_back(Pointf( r,  0.));
+        out.push_back(Pointf( r2, rx));
+        out.push_back(Pointf( rx, rx));
+        out.push_back(Pointf( rx, r2));
+        out.push_back(Pointf(0.,  r));
+        out.push_back(Pointf(-rx, r2));
+        out.push_back(Pointf(-rx, rx));
+        out.push_back(Pointf(-r2, rx));
+        out.push_back(Pointf(-r,  0.));
+        out.push_back(Pointf(-r2, -rx));
+        out.push_back(Pointf(-rx, -rx));
+        out.push_back(Pointf(-rx, -r2));
+        out.push_back(Pointf(0., -r));
+        out.push_back(Pointf( rx, -r2));
+        out.push_back(Pointf( rx, -rx));
+        out.push_back(Pointf( r2+r_inc, -rx));
+    }
+    return out;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillPlanePath.hpp b/xs/src/libslic3r/Fill/FillPlanePath.hpp
new file mode 100644
index 0000000000..79e5522a95
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillPlanePath.hpp
@@ -0,0 +1,60 @@
+#ifndef slic3r_FillPlanePath_hpp_
+#define slic3r_FillPlanePath_hpp_
+
+#include <map>
+
+#include "../libslic3r.h"
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+// The original Perl code used path generators from Math::PlanePath library:
+// http://user42.tuxfamily.org/math-planepath/
+// http://user42.tuxfamily.org/math-planepath/gallery.html
+
+class FillPlanePath : public FillWithDirection
+{
+public:
+    virtual ~FillPlanePath() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+
+protected:
+    virtual float _layer_angle(size_t idx) const { return 0.f; }
+    virtual bool  _centered() const = 0;
+    virtual Pointfs _generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y) = 0;
+};
+
+class FillArchimedeanChords : public FillPlanePath
+{
+public:
+    virtual ~FillArchimedeanChords() {}
+
+protected:
+    virtual bool  _centered() const { return true; }
+    virtual Pointfs _generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y);
+};
+
+class FillHilbertCurve : public FillPlanePath
+{
+public:
+    virtual ~FillHilbertCurve() {}
+
+protected:
+    virtual bool  _centered() const { return false; }
+    virtual Pointfs _generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y);
+};
+
+class FillOctagramSpiral : public FillPlanePath
+{
+public:
+    virtual ~FillOctagramSpiral() {}
+
+protected:
+    virtual bool  _centered() const { return true; }
+    virtual Pointfs _generate(coord_t min_x, coord_t min_y, coord_t max_x, coord_t max_y);
+};
+
+} // namespace Slic3r
+
+#endif // slic3r_FillPlanePath_hpp_
diff --git a/xs/src/libslic3r/Fill/FillRectilinear.cpp b/xs/src/libslic3r/Fill/FillRectilinear.cpp
new file mode 100644
index 0000000000..e2f996e2a2
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillRectilinear.cpp
@@ -0,0 +1,138 @@
+#include "../ClipperUtils.hpp"
+#include "../ExPolygon.hpp"
+#include "../PolylineCollection.hpp"
+#include "../Surface.hpp"
+
+#include "FillRectilinear.hpp"
+
+namespace Slic3r {
+
+Polylines FillRectilinear::fill_surface(const Surface *surface, const FillParams &params)
+{
+    // rotate polygons so that we can work with vertical lines here
+    ExPolygon expolygon = surface->expolygon;
+    std::pair<float, Point> rotate_vector = this->infill_direction(surface);
+    expolygon.rotate(- rotate_vector.first);
+    // No need to translate the polygon anyhow for the infill.
+    // The infill will be performed inside a bounding box of the expolygon and its absolute position does not matter.
+//    expolygon.translate(rotate_vector.second.x, rotate_vector.second.y);
+
+    this->_min_spacing = scale_(this->spacing);
+    assert(params.density > 0.0001f && params.density <= 1.f);
+    this->_line_spacing = coord_t(coordf_t(this->_min_spacing) / params.density);
+    this->_diagonal_distance = this->_line_spacing * 2;
+    this->_line_oscillation = this->_line_spacing - this->_min_spacing; // only for Line infill
+    BoundingBox bounding_box = expolygon.contour.bounding_box();
+    
+    // define flow spacing according to requested density
+    if (params.density > 0.9999f && !params.dont_adjust) {
+        this->_line_spacing = this->adjust_solid_spacing(bounding_box.size().x, this->_line_spacing);
+        this->spacing = unscale(this->_line_spacing);
+    } else {
+        // extend bounding box so that our pattern will be aligned with other layers
+        bounding_box.merge(Point(
+            bounding_box.min.x - (bounding_box.min.x % this->_line_spacing),
+            bounding_box.min.y - (bounding_box.min.y % this->_line_spacing)));
+    }
+
+    // generate the basic pattern
+    coord_t x_max = bounding_box.max.x + SCALED_EPSILON;
+    Lines lines;
+    for (coord_t x = bounding_box.min.x; x <= x_max; x += this->_line_spacing)
+        lines.push_back(this->_line(lines.size(), x, bounding_box.min.y, bounding_box.max.y));
+    if (this->_horizontal_lines()) {
+        coord_t y_max = bounding_box.max.y + SCALED_EPSILON;
+        for (coord_t y = bounding_box.min.y; y <= y_max; y += this->_line_spacing)
+            lines.push_back(Line(Point(bounding_box.min.x, y), Point(bounding_box.max.x, y)));
+    }
+
+    // clip paths against a slightly larger expolygon, so that the first and last paths
+    // are kept even if the expolygon has vertical sides
+    // the minimum offset for preventing edge lines from being clipped is SCALED_EPSILON;
+    // however we use a larger offset to support expolygons with slightly skewed sides and 
+    // not perfectly straight
+    //FIXME Vojtech: Update the intersecton function to work directly with lines.
+    Polylines polylines_src;
+    polylines_src.reserve(lines.size());
+    for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it) {
+        polylines_src.push_back(Polyline());
+        Points &pts = polylines_src.back().points;
+        pts.reserve(2);
+        pts.push_back(it->a);
+        pts.push_back(it->b);
+    }
+    Polylines polylines = intersection(polylines_src, offset((Polygons)expolygon, scale_(0.02)), false);
+
+    // FIXME Vojtech: This is only performed for horizontal lines, not for the vertical lines!
+    const float INFILL_OVERLAP_OVER_SPACING = 0.3f;
+    coord_t extra = coord_t(floor(this->_min_spacing * INFILL_OVERLAP_OVER_SPACING + 0.5f));
+    for (Polylines::iterator it_polyline = polylines.begin(); it_polyline != polylines.end(); ++ it_polyline) {
+        Point *first_point = &it_polyline->points.front();
+        Point *last_point  = &it_polyline->points.back();
+        if (first_point->y > last_point->y)
+            std::swap(first_point, last_point);
+        first_point->y -= extra;
+        last_point->y += extra;
+    }
+
+    // connect lines
+    if (! params.dont_connect && ! polylines.empty()) { // prevent calling leftmost_point() on empty collections
+        // offset the expolygon by max(min_spacing/2, extra)
+        ExPolygon expolygon_off;
+        {
+            ExPolygons expolygons_off = offset_ex(expolygon, this->_min_spacing/2);
+            if (! expolygons_off.empty()) {
+                // When expanding a polygon, the number of islands could only shrink. Therefore the offset_ex shall generate exactly one expanded island for one input island.
+                assert(expolygons_off.size() == 1);
+                std::swap(expolygon_off, expolygons_off.front());
+            }
+        }
+        Polylines chained = PolylineCollection::chained_path_from(
+#if SLIC3R_CPPVER >= 11
+            std::move(polylines), 
+#else
+            polylines,
+#endif
+            PolylineCollection::leftmost_point(polylines), false); // reverse allowed
+#if SLIC3R_CPPVER >= 11
+            assert(polylines.empty());
+#else
+            polylines.clear();
+#endif
+        for (Polylines::iterator it_polyline = chained.begin(); it_polyline != chained.end(); ++ it_polyline) {
+            if (! polylines.empty()) {
+                // Try to connect the lines.
+                Points &pts_end = polylines.back().points;
+                const Point &first_point = it_polyline->points.front();
+                const Point &last_point = pts_end.back();
+                // Distance in X, Y.
+                const Vector distance = first_point.vector_to(last_point);
+                // TODO: we should also check that both points are on a fill_boundary to avoid 
+                // connecting paths on the boundaries of internal regions
+                if (this->_can_connect(std::abs(distance.x), std::abs(distance.y)) && 
+                    expolygon_off.contains(Line(last_point, first_point))) {
+                    // Append the polyline.
+                    pts_end.insert(pts_end.end(), it_polyline->points.begin(), it_polyline->points.end());
+                    continue;
+                }
+            }
+            // The lines cannot be connected.
+#if SLIC3R_CPPVER >= 11
+            polylines.push_back(std::move(*it_polyline));
+#else
+            polylines.push_back(Polyline());
+            std::swap(polylines.back(), *it_polyline);
+#endif
+        }
+    }
+
+    // paths must be rotated back
+    for (Polylines::iterator it = polylines.begin(); it != polylines.end(); ++ it) {
+        // No need to translate, the absolute position is irrelevant.
+        // it->translate(- rotate_vector.second.x, - rotate_vector.second.y);
+        it->rotate(rotate_vector.first);
+    }
+    return polylines;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillRectilinear.hpp b/xs/src/libslic3r/Fill/FillRectilinear.hpp
new file mode 100644
index 0000000000..d10c4ea65d
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillRectilinear.hpp
@@ -0,0 +1,72 @@
+#ifndef slic3r_FillRectilinear_hpp_
+#define slic3r_FillRectilinear_hpp_
+
+#include "../libslic3r.h"
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+class Surface;
+
+class FillRectilinear : public FillWithDirection
+{
+public:
+    virtual ~FillRectilinear() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+
+protected:
+	coord_t _min_spacing;
+	coord_t _line_spacing;
+	// distance threshold for allowing the horizontal infill lines to be connected into a continuous path
+	coord_t _diagonal_distance;
+	// only for line infill
+	coord_t _line_oscillation;
+
+	// Enabled for the grid infill, disabled for the rectilinear and line infill.
+	virtual bool _horizontal_lines() const { return false; }
+
+	virtual Line _line(int i, coord_t x, coord_t y_min, coord_t y_max) const 
+		{ return Line(Point(x, y_min), Point(x, y_max)); }
+	
+	virtual bool _can_connect(coord_t dist_X, coord_t dist_Y) {
+	    return dist_X <= this->_diagonal_distance
+	        && dist_Y <= this->_diagonal_distance;
+    }
+};
+
+class FillLine : public FillRectilinear
+{
+public:
+    virtual ~FillLine() {}
+
+protected:
+	virtual Line _line(int i, coord_t x, coord_t y_min, coord_t y_max) const {
+		coord_t osc = (i & 1) ? this->_line_oscillation : 0;
+		return Line(Point(x - osc, y_min), Point(x + osc, y_max));
+	}
+
+	virtual bool _can_connect(coord_t dist_X, coord_t dist_Y)
+	{
+	    coord_t TOLERANCE = 10 * SCALED_EPSILON;
+    	return (dist_X >= (this->_line_spacing - this->_line_oscillation) - TOLERANCE)
+        	&& (dist_X <= (this->_line_spacing + this->_line_oscillation) + TOLERANCE)
+        	&& (dist_Y <= this->_diagonal_distance);
+    }
+};
+
+class FillGrid : public FillRectilinear
+{
+public:
+    virtual ~FillGrid() {}
+
+protected:
+	// The grid fill will keep the angle constant between the layers, see the implementation of Slic3r::Fill::Base.
+    virtual float _layer_angle(size_t idx) const { return 0.f; }
+	// Flag for Slic3r::Fill::Rectilinear to fill both directions.
+	virtual bool _horizontal_lines() const { return true; }
+};
+
+}; // namespace Slic3r
+
+#endif // slic3r_FillRectilinear_hpp_
diff --git a/xs/src/libslic3r/Fill/FillRectilinear2.cpp b/xs/src/libslic3r/Fill/FillRectilinear2.cpp
new file mode 100644
index 0000000000..1b90193054
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillRectilinear2.cpp
@@ -0,0 +1,956 @@
+#include <assert.h>
+#include <stdint.h>
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include <boost/static_assert.hpp>
+
+#include "../ClipperUtils.hpp"
+#include "../ExPolygon.hpp"
+#include "../Surface.hpp"
+
+#include "FillRectilinear2.hpp"
+
+#ifdef SLIC3R_DEBUG
+#include "SVG.hpp"
+#endif
+
+#if defined(SLIC3R_DEBUG) and defined(_WIN32)
+#include <Windows.h>
+#pragma comment(lib, "user32.lib")
+    static inline void assert_fail(const char *assertion, const char *file, unsigned line, const char *function)
+    {
+        printf("Assert: %s in function %s\nfile %s:%d\n", assertion, function, file, line);
+        if (IsDebuggerPresent()) {
+            DebugBreak();
+        } else {
+            ExitProcess(-1);
+        }
+    }
+    #undef assert
+    #define assert(expr)       \
+        ((expr)                \
+        ? static_cast<void>(0) \
+        : assert_fail (#expr, __FILE__, __LINE__, __FUNCTION__))
+#endif /* SLIC3R_DEBUG */
+
+namespace Slic3r {
+
+#ifndef clamp
+template<typename T>
+static inline T clamp(T low, T high, T x)
+{
+    return std::max<T>(low, std::min<T>(high, x));
+}
+#endif /* clamp */
+
+#ifndef sqr
+template<typename T>
+static inline T sqr(T x)
+{
+    return x * x;
+}
+#endif /* sqr */
+
+#ifndef mag2
+static inline coordf_t mag2(const Point &p)
+{
+    return sqr(coordf_t(p.x)) + sqr(coordf_t(p.y));
+}
+#endif /* mag2 */
+
+#ifndef mag
+static inline coordf_t mag(const Point &p)
+{
+    return std::sqrt(mag2(p));
+}
+#endif /* mag */
+
+enum Orientation
+{
+    ORIENTATION_CCW = 1,
+    ORIENTATION_CW = -1,
+    ORIENTATION_COLINEAR = 0
+};
+
+// Return orientation of the three points (clockwise, counter-clockwise, colinear)
+// The predicate is exact for the coord_t type, using 64bit signed integers for the temporaries.
+//FIXME Make sure the temporaries do not overflow,
+// which means, the coord_t types must not have some of the topmost bits utilized.
+static inline Orientation orient(const Point &a, const Point &b, const Point &c)
+{
+    BOOST_STATIC_ASSERT(sizeof(coord_t) * 2 == sizeof(int64_t));
+    int64_t u = int64_t(b.x) * int64_t(c.y) - int64_t(b.y) * int64_t(c.x);
+    int64_t v = int64_t(a.x) * int64_t(c.y) - int64_t(a.y) * int64_t(c.x);
+    int64_t w = int64_t(a.x) * int64_t(b.y) - int64_t(a.y) * int64_t(b.x);
+    int64_t d = u - v + w;
+    return (d > 0) ? ORIENTATION_CCW : ((d == 0) ? ORIENTATION_COLINEAR : ORIENTATION_CW);
+}
+
+// Return orientation of the polygon.
+// The input polygon must not contain duplicate points.
+static inline bool is_ccw(const Polygon &poly)
+{
+    // The polygon shall be at least a triangle.
+    assert(poly.points.size() >= 3);
+    if (poly.points.size() < 3)
+        return true;
+
+    // 1) Find the lowest lexicographical point.
+    int     imin = 0;
+    for (size_t i = 1; i < poly.points.size(); ++ i) {
+        const Point &pmin = poly.points[imin];
+        const Point &p    = poly.points[i];
+        if (p.x < pmin.x || (p.x == pmin.x && p.y < pmin.y))
+            imin = i;
+    }
+
+    // 2) Detect the orientation of the corner imin.
+    size_t iPrev = ((imin == 0) ? poly.points.size() : imin) - 1;
+    size_t iNext = ((imin + 1 == poly.points.size()) ? 0 : imin + 1);
+    Orientation o = orient(poly.points[iPrev], poly.points[imin], poly.points[iNext]);
+    // The lowest bottom point must not be collinear if the polygon does not contain duplicate points
+    // or overlapping segments.
+    assert(o != ORIENTATION_COLINEAR);
+    return o == ORIENTATION_CCW;
+}
+
+// Having a segment of a closed polygon, calculate its Euclidian length.
+// The segment indices seg1 and seg2 signify an end point of an edge in the forward direction of the loop,
+// therefore the point p1 lies on poly.points[seg1-1], poly.points[seg1] etc.
+static inline coordf_t segment_length(const Polygon &poly, size_t seg1, const Point &p1, size_t seg2, const Point &p2)
+{
+#ifdef SLIC3R_DEBUG
+    // Verify that p1 lies on seg1. This is difficult to verify precisely,
+    // but at least verify, that p1 lies in the bounding box of seg1.
+    for (size_t i = 0; i < 2; ++ i) {
+        size_t seg = (i == 0) ? seg1 : seg2;
+        Point  px  = (i == 0) ? p1   : p2;
+        Point  pa  = poly.points[((seg == 0) ? poly.points.size() : seg) - 1];
+        Point  pb  = poly.points[seg];
+        if (pa.x > pb.x)
+            std::swap(pa.x, pb.x);
+        if (pa.y > pb.y)
+            std::swap(pa.y, pb.y);
+        assert(px.x >= pa.x && px.x <= pb.x);
+        assert(px.y >= pa.y && px.y <= pb.y);
+    }
+#endif /* SLIC3R_DEBUG */
+    const Point *pPrev = &p1;
+    const Point *pThis = NULL;
+    coordf_t len = 0;
+    if (seg1 <= seg2) {
+        for (size_t i = seg1; i < seg2; ++ i, pPrev = pThis)
+           len += pPrev->distance_to(*(pThis = &poly.points[i]));
+    } else {
+        for (size_t i = seg1; i < poly.points.size(); ++ i, pPrev = pThis)
+           len += pPrev->distance_to(*(pThis = &poly.points[i]));
+        for (size_t i = 0; i < seg2; ++ i, pPrev = pThis)
+           len += pPrev->distance_to(*(pThis = &poly.points[i]));
+    }
+    len += pPrev->distance_to(p2);
+    return len;
+}
+
+// Append a segment of a closed polygon to a polyline.
+// The segment indices seg1 and seg2 signify an end point of an edge in the forward direction of the loop.
+// Only insert intermediate points between seg1 and seg2.
+static inline void polygon_segment_append(Points &out, const Polygon &polygon, size_t seg1, size_t seg2)
+{
+    if (seg1 == seg2) {
+        // Nothing to append from this segment.
+    } else if (seg1 < seg2) {
+        // Do not append a point pointed to by seg2.
+        out.insert(out.end(), polygon.points.begin() + seg1, polygon.points.begin() + seg2);
+    } else {
+        out.reserve(out.size() + seg2 + polygon.points.size() - seg1);
+        out.insert(out.end(), polygon.points.begin() + seg1, polygon.points.end());
+        // Do not append a point pointed to by seg2.
+        out.insert(out.end(), polygon.points.begin(), polygon.points.begin() + seg2);
+    }
+}
+
+// Append a segment of a closed polygon to a polyline.
+// The segment indices seg1 and seg2 signify an end point of an edge in the forward direction of the loop,
+// but this time the segment is traversed backward.
+// Only insert intermediate points between seg1 and seg2.
+static inline void polygon_segment_append_reversed(Points &out, const Polygon &polygon, size_t seg1, size_t seg2)
+{
+    if (seg1 >= seg2) {
+        out.reserve(seg1 - seg2);
+        for (size_t i = seg1; i > seg2; -- i)
+            out.push_back(polygon.points[i - 1]);
+    } else {
+        // it could be, that seg1 == seg2. In that case, append the complete loop.
+        out.reserve(out.size() + seg2 + polygon.points.size() - seg1);
+        for (size_t i = seg1; i > 0; -- i)
+            out.push_back(polygon.points[i - 1]);
+        for (size_t i = polygon.points.size(); i > seg2; -- i)
+            out.push_back(polygon.points[i - 1]);
+    }
+}
+
+// Intersection point of a vertical line with a polygon segment.
+class SegmentIntersection
+{
+public:
+    SegmentIntersection() : 
+        iContour(0),
+        iSegment(0),
+        pos(0),
+        type(UNKNOWN),
+        consumed_vertical_up(false),
+        consumed_perimeter_right(false)
+        {}
+
+    // Index of a contour in ExPolygonWithOffset, with which this vertical line intersects.
+    size_t      iContour;
+    // Index of a segment in iContour, with which this vertical line intersects.
+    size_t      iSegment;
+    // y position of the intersection.
+    coord_t     pos;
+
+    // Kind of intersection. With the original contour, or with the inner offestted contour?
+    // A vertical segment will be at least intersected by OUTER_LOW, OUTER_HIGH,
+    // but it could be intersected with OUTER_LOW, INNER_LOW, INNER_HIGH, OUTER_HIGH,
+    // and there may be more than one pair of INNER_LOW, INNER_HIGH between OUTER_LOW, OUTER_HIGH.
+    enum SegmentIntersectionType {
+        OUTER_LOW   = 0,
+        OUTER_HIGH  = 1,
+        INNER_LOW   = 2,
+        INNER_HIGH  = 3,
+        UNKNOWN     = -1
+    };
+    SegmentIntersectionType type;
+
+    // Was this segment along the y axis consumed?
+    // Up means up along the vertical segment.
+    bool consumed_vertical_up;
+    // Was a segment of the inner perimeter contour consumed?
+    // Right means right from the vertical segment.
+    bool consumed_perimeter_right;
+
+    // For the INNER_LOW type, this point may be connected to another INNER_LOW point following a perimeter contour.
+    // For the INNER_HIGH type, this point may be connected to another INNER_HIGH point following a perimeter contour.
+    // If INNER_LOW is connected to INNER_HIGH or vice versa,
+    // one has to make sure the vertical infill line does not overlap with the connecting perimeter line.
+    bool is_inner() const { return type == INNER_LOW  || type == INNER_HIGH; }
+    bool is_outer() const { return type == OUTER_LOW  || type == OUTER_HIGH; }
+    bool is_low  () const { return type == INNER_LOW  || type == OUTER_LOW; }
+    bool is_high () const { return type == INNER_HIGH || type == OUTER_HIGH; }
+
+    bool operator<(const SegmentIntersection &other) const
+        { return pos < other.pos; }
+};
+
+// A vertical line with intersection points with polygons.
+class SegmentedIntersectionLine
+{
+public:
+    // Index of this vertical intersection line.
+    size_t                              idx;
+    // x position of this vertical intersection line.
+    coord_t                             pos;
+    // List of intersection points with polygons, sorted increasingly by the y axis.
+    std::vector<SegmentIntersection>    intersections;
+};
+
+// A container maintaining an expolygon with its inner offsetted polygon.
+// The purpose of the inner offsetted polygon is to provide segments to connect the infill lines.
+struct ExPolygonWithOffset
+{
+public:
+    ExPolygonWithOffset(const ExPolygon &aexpolygon, coord_t aoffset) : expolygon(aexpolygon)
+    {
+        polygons_inner = offset((Polygons)expolygon, aoffset, 
+            CLIPPER_OFFSET_SCALE, 
+            ClipperLib::jtMiter, 
+            // for the infill pattern, don't cut the corners.
+            // default miterLimt = 3
+            10.);
+        n_contours_outer = 1 + expolygon.holes.size();
+        n_contours_inner = polygons_inner.size();
+        n_contours = n_contours_outer + n_contours_inner;
+        polygons_inner_ccw.assign(polygons_inner.size(), false);
+        for (size_t i = 0; i < polygons_inner.size(); ++ i)
+            polygons_inner_ccw[i] = is_ccw(polygons_inner[i]);
+        #ifdef SLIC3R_DEBUG
+        // Verify orientation of the expolygon.
+        assert(is_ccw(expolygon.contour));
+        for (size_t i = 0; i < expolygon.holes.size(); ++ i)
+            assert(is_ccw(expolygon.holes[i]));
+        #endif /* SLIC3R_DEBUG */
+    }
+
+    // Outer contour of the expolygon.
+    bool             is_contour_external(size_t idx) const { return idx == 0; }
+    // Any contour of the expolygon.
+    bool             is_contour_outer(size_t idx) const { return idx < n_contours_inner; }
+    // Contour of the shrunk expolygon.
+    bool             is_contour_inner(size_t idx) const { return idx >= n_contours_inner; }
+
+    const Polygon&   contour(size_t idx) const {
+        return  is_contour_external(idx) ? expolygon.contour :
+               (is_contour_outer(idx) ? expolygon.holes[idx - 1] : polygons_inner[idx - n_contours_inner]);
+    }
+
+    bool             is_contour_ccw(size_t idx) const {
+        return  is_contour_external(idx) || (is_contour_inner(idx) && polygons_inner_ccw[idx - n_contours_inner]);
+    }
+
+    const ExPolygon &expolygon;
+    Polygons         polygons_inner;
+
+    size_t           n_contours_outer;
+    size_t           n_contours_inner;
+    size_t           n_contours;
+
+protected:
+
+    // For each polygon of polygons_inner, remember its orientation.
+    std::vector<unsigned char> polygons_inner_ccw;
+};
+
+// For a vertical line, an inner contour and an intersection point,
+// find an intersection point on the previous resp. next vertical line.
+// The intersection point is connected with the prev resp. next intersection point with iInnerContour.
+// Return -1 if there is no such point on the previous resp. next vertical line.
+static inline int intersection_on_prev_next_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset,
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    bool                                           dir_is_next)
+{
+    size_t iVerticalLineOther = iVerticalLine;
+    if (dir_is_next) {
+        if (++ iVerticalLineOther == segs.size())
+            // No successive vertical line.
+            return -1;
+    } else if (iVerticalLineOther -- == 0) {
+        // No preceding vertical line.
+        return -1;
+    }
+
+    const SegmentedIntersectionLine &il    = segs[iVerticalLine];
+    const SegmentIntersection       &itsct = il.intersections[iIntersection];
+    const SegmentedIntersectionLine &il2   = segs[iVerticalLineOther];
+    const Polygon                   &poly  = poly_with_offset.contour(iInnerContour);
+    const bool                       ccw   = poly_with_offset.is_contour_ccw(iInnerContour);
+    // Resulting index of an intersection point on il2.
+    int                              out   = -1;
+    // Find an intersection point on iVerticalLineOther, intersecting iInnerContour
+    // at the same orientation as iIntersection, and being closest to iIntersection
+    // in the number of contour segments, when following the direction of the contour.
+    int                              dmin  = std::numeric_limits<int>::max();
+    for (size_t i = 0; i < il2.intersections.size(); ++ i) {
+        const SegmentIntersection &itsct2 = il2.intersections[i];
+        if (itsct.iContour == itsct2.iContour && itsct.type == itsct2.type) {
+            // The intersection points lie on the same contour and have the same orientation.
+            // Find the intersection point with a shortest path in the direction of the contour.
+            int d = int(itsct2.iSegment) - int(itsct.iSegment);
+            if (ccw != dir_is_next)
+                d = - d;
+            if (d < 0)
+                d += int(poly.points.size());
+            if (d < dmin) {
+                out = i;
+                dmin = d;
+            }
+        }
+    }
+    //FIXME this routine is not asymptotic optimal, it will be slow if there are many intersection points along the line.
+    return out;
+}
+
+static inline int intersection_on_prev_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset, 
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection)
+{
+    return intersection_on_prev_next_vertical_line(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, false);
+}
+
+static inline intersection_on_next_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset, 
+    const std::vector<SegmentedIntersectionLine>  &segs, 
+    size_t                                         iVerticalLine, 
+    size_t                                         iInnerContour, 
+    size_t                                         iIntersection)
+{
+    return intersection_on_prev_next_vertical_line(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, true);
+}
+
+// Find an intersection on a previous line, but return -1, if the connecting segment of a perimeter was already extruded.
+static inline int intersection_unused_on_prev_next_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset,
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    bool                                           dir_is_next)
+{
+    int iIntersectionOther = intersection_on_prev_next_vertical_line(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, dir_is_next);
+    if (iIntersectionOther == -1)
+        return -1;
+    //FIXME this routine will propose a connecting line even if the connecting perimeter segment intersects iVertical line multiple times before reaching iIntersectionOther.
+    assert(dir_is_next ? (iVerticalLine + 1 < segs.size()) : (iVerticalLine > 0));
+    const SegmentedIntersectionLine &il_this      = segs[iVerticalLine];
+    const SegmentIntersection       &itsct_this   = il_this.intersections[iIntersection];
+    const SegmentedIntersectionLine &il_other     = segs[dir_is_next ? (iVerticalLine+1) : (iVerticalLine-1)];
+    const SegmentIntersection       &itsct_other  = il_other.intersections[iIntersectionOther];
+    assert(itsct_other.is_inner());
+    assert(itsct_other.is_low() || iIntersectionOther > 1);
+    if (dir_is_next ? itsct_this.consumed_perimeter_right : itsct_other.consumed_perimeter_right)
+        // This perimeter segment was already consumed.
+        return -1;
+    if (itsct_other.is_low() ? itsct_other.consumed_vertical_up : il_other.intersections[iIntersectionOther-1].consumed_vertical_up)
+        // This vertical segment was already consumed.
+        return -1;
+    return iIntersectionOther;
+}
+
+static inline intersection_unused_on_prev_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset, 
+    const std::vector<SegmentedIntersectionLine>  &segs, 
+    size_t                                         iVerticalLine, 
+    size_t                                         iInnerContour, 
+    size_t                                         iIntersection)
+{
+    return intersection_unused_on_prev_next_vertical_line(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, false);
+}
+
+static inline intersection_unused_on_next_vertical_line(
+    const ExPolygonWithOffset                     &poly_with_offset, 
+    const std::vector<SegmentedIntersectionLine>  &segs, 
+    size_t                                         iVerticalLine, 
+    size_t                                         iInnerContour, 
+    size_t                                         iIntersection)
+{
+    return intersection_unused_on_prev_next_vertical_line(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, true);
+}
+
+// Measure an Euclidian length of a perimeter segment when going from iIntersection to iIntersection2.
+static inline coordf_t measure_perimeter_prev_next_segment_length(
+    const ExPolygonWithOffset                     &poly_with_offset, 
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    size_t                                         iIntersection2,
+    bool                                           dir_is_next)
+{
+    size_t iVerticalLineOther = iVerticalLine;
+    if (dir_is_next) {
+        if (++ iVerticalLineOther == segs.size())
+            // No successive vertical line.
+            return coordf_t(-1);
+    } else if (iVerticalLineOther -- == 0) {
+        // No preceding vertical line.
+        return coordf_t(-1);
+    }
+
+    const SegmentedIntersectionLine &il     = segs[iVerticalLine];
+    const SegmentIntersection       &itsct  = il.intersections[iIntersection];
+    const SegmentedIntersectionLine &il2    = segs[iVerticalLineOther];
+    const SegmentIntersection       &itsct2 = il2.intersections[iIntersection2];
+    const Polygon                   &poly   = poly_with_offset.contour(iInnerContour);
+    const bool                       ccw    = poly_with_offset.is_contour_ccw(iInnerContour);
+    assert(itsct.type == itsct2.type);
+    assert(itsct.iContour == itsct2.iContour);
+    assert(itsct.is_inner());
+    const bool                       forward = (itsct.is_low() == ccw) == dir_is_next;
+
+    Point p1(il.pos, itsct.pos);
+    Point p2(il2.pos, itsct2.pos);
+    return forward ?
+        segment_length(poly, itsct .iSegment, p1, itsct2.iSegment, p2) :
+        segment_length(poly, itsct2.iSegment, p2, itsct .iSegment, p1);
+}
+
+static inline coordf_t measure_perimeter_prev_segment_length(
+    const ExPolygonWithOffset                     &poly_with_offset,
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    size_t                                         iIntersection2)
+{
+    return measure_perimeter_prev_next_segment_length(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, iIntersection2, false);
+}
+
+static inline coordf_t measure_perimeter_next_segment_length(
+    const ExPolygonWithOffset                     &poly_with_offset,
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    size_t                                         iIntersection2)
+{
+    return measure_perimeter_prev_next_segment_length(poly_with_offset, segs, iVerticalLine, iInnerContour, iIntersection, iIntersection2, true);
+}
+
+// Append the points of a perimeter segment when going from iIntersection to iIntersection2.
+// The first point (the point of iIntersection) will not be inserted,
+// the last point will be inserted.
+static inline void emit_perimeter_prev_next_segment(
+    const ExPolygonWithOffset                     &poly_with_offset,
+    const std::vector<SegmentedIntersectionLine>  &segs,
+    size_t                                         iVerticalLine,
+    size_t                                         iInnerContour,
+    size_t                                         iIntersection,
+    size_t                                         iIntersection2,
+    Polyline                                      &out,
+    bool                                           dir_is_next)
+{
+    size_t iVerticalLineOther = iVerticalLine;
+    if (dir_is_next) {
+        ++ iVerticalLineOther;
+        assert(iVerticalLineOther < segs.size());
+    } else {
+        assert(iVerticalLineOther > 0);
+        -- iVerticalLineOther;
+    }
+
+    const SegmentedIntersectionLine &il     = segs[iVerticalLine];
+    const SegmentIntersection       &itsct  = il.intersections[iIntersection];
+    const SegmentedIntersectionLine &il2    = segs[iVerticalLineOther];
+    const SegmentIntersection       &itsct2 = il2.intersections[iIntersection2];
+    const Polygon                   &poly   = poly_with_offset.contour(iInnerContour);
+    const bool                       ccw    = poly_with_offset.is_contour_ccw(iInnerContour);
+    assert(itsct.type == itsct2.type);
+    assert(itsct.iContour == itsct2.iContour);
+    assert(itsct.is_inner());
+    const bool                       forward = (itsct.is_low() == ccw) == dir_is_next;
+    // Do not append the first point.
+    // out.points.push_back(Point(il.pos, itsct.pos));
+    if (forward)
+        polygon_segment_append(out.points, poly, itsct.iSegment, itsct2.iSegment);
+    else
+        polygon_segment_append_reversed(out.points, poly, itsct.iSegment, itsct2.iSegment);
+    // Append the last point.
+    out.points.push_back(Point(il2.pos, itsct2.pos));
+}
+
+Polylines FillRectilinear2::fill_surface(const Surface *surface, const FillParams &params)
+{
+    // rotate polygons so that we can work with vertical lines here
+    ExPolygon expolygon = surface->expolygon;
+    std::pair<float, Point> rotate_vector = this->infill_direction(surface);
+    expolygon.rotate(- rotate_vector.first);
+    // No need to translate the polygon anyhow for the infill.
+    // The infill will be performed inside a bounding box of the expolygon and its absolute position does not matter.
+//    expolygon.translate(rotate_vector.second.x, rotate_vector.second.y);
+
+    this->_min_spacing = scale_(this->spacing);
+    assert(params.density > 0.0001f && params.density <= 1.f);
+    this->_line_spacing = coord_t(coordf_t(this->_min_spacing) / params.density);
+    this->_diagonal_distance = this->_line_spacing * 2;
+    BoundingBox bounding_box = expolygon.contour.bounding_box();
+    
+    // define flow spacing according to requested density
+    if (params.density > 0.9999f && !params.dont_adjust) {
+        this->_line_spacing = this->adjust_solid_spacing(bounding_box.size().x, this->_line_spacing);
+        this->spacing = unscale(this->_line_spacing);
+    } else {
+        // extend bounding box so that our pattern will be aligned with other layers
+        bounding_box.merge(Point(
+            bounding_box.min.x - (bounding_box.min.x % this->_line_spacing),
+            bounding_box.min.y - (bounding_box.min.y % this->_line_spacing)));
+    }
+
+    // Intersect a set of euqally spaced vertical lines wiht expolygon.
+    size_t  n_vlines = (bounding_box.max.x - bounding_box.min.x + SCALED_EPSILON) / this->_line_spacing;
+    coord_t x0 = bounding_box.min.x + this->_line_spacing;
+    // On these polygons the infill lines will be connected.
+    ExPolygonWithOffset poly_with_offset(expolygon, - _min_spacing / 2);
+
+#ifdef SLIC3R_DEBUG
+    char path[2048];
+    static int iRun = 0;
+    sprintf(path, "out/FillRectilinear2-%d.svg", iRun);
+    BoundingBox bbox_svg = expolygon.contour.bounding_box();
+    bbox_svg.min.x -= coord_t(1. / SCALING_FACTOR);
+    bbox_svg.min.y -= coord_t(1. / SCALING_FACTOR);
+    bbox_svg.max.x += coord_t(1. / SCALING_FACTOR);
+    bbox_svg.max.y += coord_t(1. / SCALING_FACTOR);
+    ::Slic3r::SVG svg(path, bbox_svg);
+    svg.draw(expolygon.lines());
+    svg.draw(poly_with_offset.polygons_inner);
+    {
+        char path2[2048];
+        sprintf(path2, "out/FillRectilinear2-initial-%d.svg", iRun);
+        ::Slic3r::SVG svg(path2, bbox_svg);
+        svg.draw(expolygon.lines());
+        svg.draw(poly_with_offset.polygons_inner);
+        svg.Close();
+    }
+    iRun ++;
+#endif /* SLIC3R_DEBUG */
+
+    // For each contour
+    // Allocate the storage for the segments.
+    std::vector<SegmentedIntersectionLine> segs(n_vlines, SegmentedIntersectionLine());
+    for (size_t i = 0; i < n_vlines; ++ i) {
+        segs[i].idx = i;
+        segs[i].pos = x0 + i * this->_line_spacing;
+    }
+    for (size_t iContour = 0; iContour < poly_with_offset.n_contours; ++ iContour) {
+        const Points &contour = poly_with_offset.contour(iContour);
+        if (contour.size() < 2)
+            continue;
+        // For each segment
+        for (size_t iSegment = 0; iSegment < contour.size(); ++ iSegment) {
+            size_t iPrev = ((iSegment == 0) ? contour.size() : iSegment) - 1;
+            const Point &p1 = contour[iPrev];
+            const Point &p2 = contour[iSegment];
+            // Which of the equally spaced vertical lines is intersected by this segment?
+            coord_t l = p1.x;
+            coord_t r = p2.x;
+            if (l > r)
+                std::swap(l, r);
+            // il, ir are the left / right indices of vertical lines intersecting a segment
+            int il = (l - x0) / this->_line_spacing;
+            while (il * this->_line_spacing + x0 < l)
+                ++ il;
+            il = std::max(int(0), il);
+            int ir = (r - x0 + this->_line_spacing) / this->_line_spacing;
+            while (ir * this->_line_spacing + x0 > r)
+                -- ir;
+            ir = std::min(int(segs.size()) - 1, ir);
+            if (il > ir)
+                // No vertical line intersects this segment.
+                continue;
+            assert(il >= 0 && il < segs.size());
+            assert(ir >= 0 && ir < segs.size());
+            if (l == r) {
+                // The segment is vertical.
+                SegmentIntersection is;
+                is.iContour = iContour;
+                is.iSegment = iSegment;
+                is.pos = p1.y;
+                segs[il].intersections.push_back(is);
+                is.pos = p2.y;
+                segs[il].intersections.push_back(is);
+                continue;
+            }
+            for (int i = il; i <= ir; ++ i) {
+                SegmentIntersection is;
+                is.iContour = iContour;
+                is.iSegment = iSegment;
+                assert(l <= segs[i].pos);
+                assert(r >= segs[i].pos);
+                // Calculate the intersection position in y axis. x is known.
+                double t = double(segs[i].pos - p1.x) / double(p2.x - p1.x);
+                assert(t > -0.000001 && t < 1.000001);
+                t = clamp(0., 1., t);
+                coord_t lo = p1.y;
+                coord_t hi = p2.y;
+                if (lo > hi)
+                    std::swap(lo, hi);
+                is.pos = p1.y + coord_t(t * double(p2.y - p1.y));
+                assert(is.pos > lo - 0.000001 && is.pos < hi + 0.000001);
+                is.pos = clamp(lo, hi, is.pos);
+                segs[i].intersections.push_back(is);
+            }
+        }
+    }
+
+    // Sort the intersections along their segments, specify the intersection types.
+    for (size_t i_seg = 0; i_seg < segs.size(); ++ i_seg) {
+        SegmentedIntersectionLine &sil = segs[i_seg];
+        // Sort the intersection points. This needs to be verified, because the intersection points were calculated
+        // using imprecise arithmetics.
+        std::sort(sil.intersections.begin(), sil.intersections.end());
+        // Verify the order, bubble sort the intersections until sorted.
+        bool modified = false;
+        do {
+            modified = false;
+            for (size_t i = 1; i < sil.intersections.size(); ++ i) {
+                size_t iContour1 = sil.intersections[i-1].iContour;
+                size_t iContour2 = sil.intersections[i].iContour;
+                const Points &contour1 = poly_with_offset.contour(iContour1);
+                const Points &contour2 = poly_with_offset.contour(iContour2);
+                size_t iSegment1 = sil.intersections[i-1].iSegment;
+                size_t iPrev1    = ((iSegment1 == 0) ? contour1.size() : iSegment1) - 1;
+                size_t iSegment2 = sil.intersections[i].iSegment;
+                size_t iPrev2    = ((iSegment2 == 0) ? contour2.size() : iSegment2) - 1;
+                bool   swap = false;
+                if (iContour1 == iContour2 && iSegment1 == iSegment2) {
+                    // The same segment, it has to be vertical.
+                    assert(iPrev1 == iPrev2);
+                    swap = contour1[iPrev1].y > contour1[iContour1].y;
+                    #ifdef SLIC3R_DEBUG
+                    if (swap)
+                        printf("Swapping when single vertical segment\n");
+                    #endif
+                } else {
+                    // Segments are in a general position. Here an exact airthmetics may come into play.
+                    coord_t y1max = std::max(contour1[iPrev1].y, contour1[iSegment1].y);
+                    coord_t y2min = std::min(contour2[iPrev2].y, contour2[iSegment2].y);
+                    if (y1max < y2min) {
+                        // The segments are separated, nothing to do.
+                    } else {
+                        // Use an exact predicate to verify, that segment1 is below segment2.
+                        const Point *a = &contour1[iPrev1];
+                        const Point *b = &contour1[iSegment1];
+                        const Point *c = &contour2[iPrev2];
+                        const Point *d = &contour2[iSegment2];
+#ifdef SLIC3R_DEBUG
+                        const Point  x1(sil.pos, sil.intersections[i-1].pos);
+                        const Point  x2(sil.pos, sil.intersections[i  ].pos);
+                        bool successive = false;
+#endif /* SLIC3R_DEBUG */
+                        if (a->x > b->x)
+                        	std::swap(a, b);
+                        if (c->x > d->x)
+                            std::swap(c, d);
+                        bool upper_more_left = false;
+                        if (a->x > c->x) {
+                            upper_more_left = true;
+                            std::swap(a, c);
+                            std::swap(b, d);
+                        }
+                        if (a == c || b == c) {
+                            assert(iContour1 == iContour2);
+                            assert(iSegment1 == iPrev2 || iPrev1 == iSegment2);
+                            std::swap(c, d);
+                            assert(a != c && b != c);
+#ifdef SLIC3R_DEBUG
+                            successive = true;
+#endif /* SLIC3R_DEBUG */
+                        }
+                        Orientation o = orient(*a, *b, *c);
+                        assert(! ORIENTATION_COLINEAR);
+                        swap = upper_more_left != (o == ORIENTATION_CW);
+#ifdef SLIC3R_DEBUG
+                        if (swap)
+                            printf(successive ? 
+                                "Swapping when iContour1 == iContour2 and successive segments\n" :
+                                "Swapping when exact predicate\n");
+#endif
+                    }
+                }
+                if (swap) {
+                    // Swap the intersection points, but keep the original positions, so they are sorted.
+                    std::swap(sil.intersections[i-1], sil.intersections[i]);
+                    std::swap(sil.intersections[i-1].pos, sil.intersections[i].pos);
+                    modified = true;
+                }
+            }
+        } while (modified);
+        // Assign the intersection types.
+        for (size_t i = 0; i < sil.intersections.size(); ++ i) {
+            // What is the orientation of the segment at the intersection point?
+            size_t iContour = sil.intersections[i].iContour;
+            const Points &contour = poly_with_offset.contour(iContour);
+            size_t iSegment = sil.intersections[i].iSegment;
+            size_t iPrev    = ((iSegment == 0) ? contour.size() : iSegment) - 1;
+            coord_t dir = contour[iSegment].x - contour[iPrev].x;
+            bool ccw = poly_with_offset.is_contour_ccw(iContour);
+            bool low = (dir > 0) == ccw;
+            sil.intersections[i].type = poly_with_offset.is_contour_outer(iContour) ? 
+                (low ? SegmentIntersection::OUTER_LOW : SegmentIntersection::OUTER_HIGH) :
+                (low ? SegmentIntersection::INNER_LOW : SegmentIntersection::INNER_HIGH);
+        }
+    }
+
+#ifdef SLIC3R_DEBUG
+    // Verify the segments & paint them.
+    for (size_t i_seg = 0; i_seg < segs.size(); ++ i_seg) {
+        SegmentedIntersectionLine &sil = segs[i_seg];
+        // The intersection points have to be even.
+        assert((sil.intersections.size() & 1) == 0);
+        for (size_t i = 0; i < sil.intersections.size();) {
+            // An intersection segment crossing the bigger contour may cross the inner offsetted contour even number of times.
+            assert(sil.intersections[i].type == SegmentIntersection::OUTER_LOW);
+            size_t j = i + 1;
+            assert(j < sil.intersections.size());
+            assert(sil.intersections[j].type == SegmentIntersection::INNER_LOW || sil.intersections[j].type == SegmentIntersection::OUTER_HIGH);
+            for (; j < sil.intersections.size() && sil.intersections[j].is_inner(); ++ j) ;
+            assert(j < sil.intersections.size());
+            assert((j & 1) == 1);
+            assert(sil.intersections[j].type == SegmentIntersection::OUTER_HIGH);
+            assert(i + 1 == j || sil.intersections[j - 1].type == SegmentIntersection::INNER_HIGH);
+            if (i + 1 == j) {
+                svg.draw(Line(Point(sil.pos, sil.intersections[i].pos), Point(sil.pos, sil.intersections[j].pos)), "blue");
+            } else {
+                svg.draw(Line(Point(sil.pos, sil.intersections[i].pos), Point(sil.pos, sil.intersections[i+1].pos)), "green");
+                svg.draw(Line(Point(sil.pos, sil.intersections[i+1].pos), Point(sil.pos, sil.intersections[j-1].pos)), (j - i + 1 > 4) ? "yellow" : "magenta");
+                svg.draw(Line(Point(sil.pos, sil.intersections[j-1].pos), Point(sil.pos, sil.intersections[j].pos)), "green");
+            }
+            i = j + 1;
+        }
+    }
+    svg.Close();
+#endif /* SLIC3R_DEBUG */
+
+    // Now construct a graph.
+    // Find the first point.
+    //FIXME ideally one would plan the initial point to be closest to the current print head position.
+    size_t    i_vline = 0;
+    size_t    i_intersection = size_t(-1);
+    // Follow the line, connect the lines into a graph.
+    // Until no new line could be added to the output path:
+    Point     pointLast;
+    Polylines polylines_out;
+    Polyline *polyline_current = NULL;
+    for (;;) {
+        if (i_intersection == size_t(-1)) {
+            // The path has been interrupted. Find a next starting point, closest to the previous extruder position.
+            coordf_t dist2min = std::numeric_limits<coordf_t>().max();
+            for (size_t i_vline2 = 0; i_vline2 < segs.size(); ++ i_vline2) {
+                const SegmentedIntersectionLine &seg = segs[i_vline2];
+                if (! seg.intersections.empty()) {
+                    assert(seg.intersections.size() > 1);
+                    // Even number of intersections with the loops.
+                    assert((seg.intersections.size() & 1) == 0);
+                    assert(seg.intersections.front().type == SegmentIntersection::OUTER_LOW);
+                    for (size_t i = 0; i < seg.intersections.size(); ++ i) {
+                        const SegmentIntersection &intrsctn = seg.intersections[i];
+                        if (intrsctn.is_outer()) {
+                            assert(intrsctn.is_low() || i > 0);
+                            bool consumed = intrsctn.is_low() ? 
+                                intrsctn.consumed_vertical_up : 
+                                seg.intersections[i-1].consumed_vertical_up;
+                            if (! consumed) {
+                                coordf_t dist2 = sqr(coordf_t(pointLast.x - seg.pos)) + sqr(coordf_t(pointLast.y - intrsctn.pos));
+                                if (dist2 < dist2min) {
+                                    dist2min = dist2;
+                                    i_vline = i_vline2;
+                                    i_intersection = i;
+                                    if (polylines_out.empty()) {
+                                        // Initial state, take the first line, which is the first from the left.
+                                        goto found;
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+            if (i_intersection == size_t(-1))
+                // We are finished.
+                break;
+        found:
+            // Start a new path.
+            polylines_out.push_back(Polyline());
+            polyline_current = &polylines_out.back();
+            // Emit the first point of a path.
+            pointLast = Point(segs[i_vline].pos, segs[i_vline].intersections[i_intersection].pos);
+            polyline_current->points.push_back(pointLast);
+        }
+
+        // From the initial point (i_vline, i_intersection), follow a path.
+        SegmentedIntersectionLine &seg      = segs[i_vline];
+        SegmentIntersection       *intrsctn = &seg.intersections[i_intersection];
+        bool going_up = intrsctn->is_low();
+        bool try_connect = false;
+        if (going_up) {
+            assert(! intrsctn->consumed_vertical_up);
+            assert(i_intersection + 1 < seg.intersections.size());
+            // Step back to the beginning of the vertical segment to mark it as consumed.
+            if (intrsctn->is_inner()) {
+                assert(i_intersection > 0);
+                -- intrsctn;
+                -- i_intersection;
+            }
+            // Consume the complete vertical segment up to the outer contour.
+            do {
+                intrsctn->consumed_vertical_up = true;
+                ++ intrsctn;
+                ++ i_intersection;
+                assert(i_intersection < seg.intersections.size());
+            } while (intrsctn->type != SegmentIntersection::OUTER_HIGH);
+            if ((intrsctn - 1)->is_inner()) {
+                // Step back.
+                -- intrsctn;
+                -- i_intersection;
+                assert(intrsctn->type == SegmentIntersection::INNER_HIGH);
+                try_connect = true;
+            }
+        } else {
+            // Going down.
+            assert(intrsctn->is_high());
+            assert(i_intersection > 0);
+            assert(! (intrsctn - 1)->consumed_vertical_up);
+            // Consume the complete vertical segment up to the outer contour.
+            if (intrsctn->is_inner())
+                intrsctn->consumed_vertical_up = true;
+            do {
+                assert(i_intersection > 0);
+                -- intrsctn;
+                -- i_intersection;
+                intrsctn->consumed_vertical_up = true;
+            } while (intrsctn->type != SegmentIntersection::OUTER_LOW);
+            if ((intrsctn + 1)->is_inner()) {
+                // Step back.
+                ++ intrsctn;
+                ++ i_intersection;
+                assert(intrsctn->type == SegmentIntersection::INNER_LOW);
+                try_connect = true;
+            }
+        }
+        if (try_connect) {
+            // Decide, whether to finish the segment, or whether to follow the perimeter.
+            int iPrev = intersection_unused_on_prev_vertical_line(poly_with_offset, segs, i_vline, intrsctn->iContour, i_intersection);
+            int iNext = intersection_unused_on_next_vertical_line(poly_with_offset, segs, i_vline, intrsctn->iContour, i_intersection);
+            if (iPrev != -1 || iNext != -1) {
+                // Zig zag
+                coord_t distPrev = (iPrev == -1) ? std::numeric_limits<coord_t>::max() :
+                    measure_perimeter_prev_segment_length(poly_with_offset, segs, i_vline, intrsctn->iContour, i_intersection, iPrev);
+                coord_t distNext = (iNext == -1) ? std::numeric_limits<coord_t>::max() :
+                    measure_perimeter_next_segment_length(poly_with_offset, segs, i_vline, intrsctn->iContour, i_intersection, iNext);
+                // Take the shorter path.
+                bool take_next = (iPrev != -1 && iNext != -1) ? (distNext < distPrev) : distNext != -1;
+                assert(intrsctn->is_inner());
+                polyline_current->points.push_back(Point(seg.pos, intrsctn->pos));
+                emit_perimeter_prev_next_segment(poly_with_offset, segs, i_vline, intrsctn->iContour, i_intersection, take_next ? iNext : iPrev, *polyline_current, take_next);
+                // Mark both the left and right connecting segment as consumed, because one cannot go to this intersection point as it has been consumed.
+                if (iPrev != -1)
+                    segs[i_vline-1].intersections[iPrev].consumed_perimeter_right = true;
+                if (iNext != -1)
+                    intrsctn->consumed_perimeter_right = true;
+                //FIXME consume the left / right connecting segments at the other end of this line? Currently it is not critical because a perimeter segment is not followed if the vertical segment at the other side has already been consumed.
+                // Advance to the neighbor line.
+                if (take_next) {
+                    ++ i_vline;
+                    i_intersection = iNext;
+                } else {
+                    -- i_vline;
+                    i_intersection = iPrev;
+                }
+                continue;
+            }
+            // Take the complete line up to the outer contour.
+            if (going_up)
+                ++ intrsctn;
+            else
+                -- intrsctn;
+        }
+        // Finish the current vertical line,
+        // reset the current vertical line to pick a new starting point in the next round.
+        assert(intrsctn->is_outer());
+        assert(intrsctn->is_high() == going_up);
+        pointLast = Point(seg.pos, intrsctn->pos);
+        polyline_current->points.push_back(pointLast);
+        intrsctn = NULL;
+        i_intersection = -1;
+        polyline_current = NULL;
+    }
+
+    // paths must be rotated back
+    for (Polylines::iterator it = polylines_out.begin(); it != polylines_out.end(); ++ it) {
+        // No need to translate, the absolute position is irrelevant.
+        // it->translate(- rotate_vector.second.x, - rotate_vector.second.y);
+        it->rotate(rotate_vector.first);
+    }
+    return polylines_out;
+}
+
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/Fill/FillRectilinear2.hpp b/xs/src/libslic3r/Fill/FillRectilinear2.hpp
new file mode 100644
index 0000000000..0c4473e4a8
--- /dev/null
+++ b/xs/src/libslic3r/Fill/FillRectilinear2.hpp
@@ -0,0 +1,37 @@
+#ifndef slic3r_FillRectilinear2_hpp_
+#define slic3r_FillRectilinear2_hpp_
+
+#include "../libslic3r.h"
+
+#include "FillBase.hpp"
+
+namespace Slic3r {
+
+class Surface;
+
+class FillRectilinear2 : public FillWithDirection
+{
+public:
+    virtual ~FillRectilinear2() {}
+    virtual Polylines fill_surface(const Surface *surface, const FillParams &params);
+
+protected:
+	coord_t _min_spacing;
+	coord_t _line_spacing;
+	// distance threshold for allowing the horizontal infill lines to be connected into a continuous path
+	coord_t _diagonal_distance;
+};
+
+class FillGrid2 : public FillRectilinear2
+{
+public:
+    virtual ~FillGrid2() {}
+
+protected:
+	// The grid fill will keep the angle constant between the layers, see the implementation of Slic3r::Fill::Base.
+    virtual float _layer_angle(size_t idx) const { return 0.f; }
+};
+
+}; // namespace Slic3r
+
+#endif // slic3r_FillRectilinear2_hpp_
diff --git a/xs/src/libslic3r/Geometry.cpp b/xs/src/libslic3r/Geometry.cpp
index 6d864f6313..cdf913a277 100644
--- a/xs/src/libslic3r/Geometry.cpp
+++ b/xs/src/libslic3r/Geometry.cpp
@@ -11,12 +11,186 @@
 #include <map>
 #include <set>
 #include <utility>
+#include <stack>
 #include <vector>
 
 #ifdef SLIC3R_DEBUG
 #include "SVG.hpp"
 #endif
 
+#ifdef SLIC3R_DEBUG
+namespace boost { namespace polygon {
+
+// The following code for the visualization of the boost Voronoi diagram is based on:
+//
+// Boost.Polygon library voronoi_graphic_utils.hpp header file
+//          Copyright Andrii Sydorchuk 2010-2012.
+// Distributed under the Boost Software License, Version 1.0.
+//    (See accompanying file LICENSE_1_0.txt or copy at
+//          http://www.boost.org/LICENSE_1_0.txt)
+template <typename CT>
+class voronoi_visual_utils {
+ public:
+  // Discretize parabolic Voronoi edge.
+  // Parabolic Voronoi edges are always formed by one point and one segment
+  // from the initial input set.
+  //
+  // Args:
+  //   point: input point.
+  //   segment: input segment.
+  //   max_dist: maximum discretization distance.
+  //   discretization: point discretization of the given Voronoi edge.
+  //
+  // Template arguments:
+  //   InCT: coordinate type of the input geometries (usually integer).
+  //   Point: point type, should model point concept.
+  //   Segment: segment type, should model segment concept.
+  //
+  // Important:
+  //   discretization should contain both edge endpoints initially.
+  template <class InCT1, class InCT2,
+            template<class> class Point,
+            template<class> class Segment>
+  static
+  typename enable_if<
+    typename gtl_and<
+      typename gtl_if<
+        typename is_point_concept<
+          typename geometry_concept< Point<InCT1> >::type
+        >::type
+      >::type,
+      typename gtl_if<
+        typename is_segment_concept<
+          typename geometry_concept< Segment<InCT2> >::type
+        >::type
+      >::type
+    >::type,
+    void
+  >::type discretize(
+      const Point<InCT1>& point,
+      const Segment<InCT2>& segment,
+      const CT max_dist,
+      std::vector< Point<CT> >* discretization) {
+    // Apply the linear transformation to move start point of the segment to
+    // the point with coordinates (0, 0) and the direction of the segment to
+    // coincide the positive direction of the x-axis.
+    CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
+    CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
+    CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y;
+
+    // Compute x-coordinates of the endpoints of the edge
+    // in the transformed space.
+    CT projection_start = sqr_segment_length *
+        get_point_projection((*discretization)[0], segment);
+    CT projection_end = sqr_segment_length *
+        get_point_projection((*discretization)[1], segment);
+
+    // Compute parabola parameters in the transformed space.
+    // Parabola has next representation:
+    // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
+    CT point_vec_x = cast(x(point)) - cast(x(low(segment)));
+    CT point_vec_y = cast(y(point)) - cast(y(low(segment)));
+    CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y;
+    CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x;
+
+    // Save the last point.
+    Point<CT> last_point = (*discretization)[1];
+    discretization->pop_back();
+
+    // Use stack to avoid recursion.
+    std::stack<CT> point_stack;
+    point_stack.push(projection_end);
+    CT cur_x = projection_start;
+    CT cur_y = parabola_y(cur_x, rot_x, rot_y);
+
+    // Adjust max_dist parameter in the transformed space.
+    const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length;
+    while (!point_stack.empty()) {
+      CT new_x = point_stack.top();
+      CT new_y = parabola_y(new_x, rot_x, rot_y);
+
+      // Compute coordinates of the point of the parabola that is
+      // furthest from the current line segment.
+      CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x;
+      CT mid_y = parabola_y(mid_x, rot_x, rot_y);
+
+      // Compute maximum distance between the given parabolic arc
+      // and line segment that discretize it.
+      CT dist = (new_y - cur_y) * (mid_x - cur_x) -
+          (new_x - cur_x) * (mid_y - cur_y);
+      dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
+          (new_x - cur_x) * (new_x - cur_x));
+      if (dist <= max_dist_transformed) {
+        // Distance between parabola and line segment is less than max_dist.
+        point_stack.pop();
+        CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) /
+            sqr_segment_length + cast(x(low(segment)));
+        CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) /
+            sqr_segment_length + cast(y(low(segment)));
+        discretization->push_back(Point<CT>(inter_x, inter_y));
+        cur_x = new_x;
+        cur_y = new_y;
+      } else {
+        point_stack.push(mid_x);
+      }
+    }
+
+    // Update last point.
+    discretization->back() = last_point;
+  }
+
+ private:
+  // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
+  static CT parabola_y(CT x, CT a, CT b) {
+    return ((x - a) * (x - a) + b * b) / (b + b);
+  }
+
+  // Get normalized length of the distance between:
+  //   1) point projection onto the segment
+  //   2) start point of the segment
+  // Return this length divided by the segment length. This is made to avoid
+  // sqrt computation during transformation from the initial space to the
+  // transformed one and vice versa. The assumption is made that projection of
+  // the point lies between the start-point and endpoint of the segment.
+  template <class InCT,
+            template<class> class Point,
+            template<class> class Segment>
+  static
+  typename enable_if<
+    typename gtl_and<
+      typename gtl_if<
+        typename is_point_concept<
+          typename geometry_concept< Point<int> >::type
+        >::type
+      >::type,
+      typename gtl_if<
+        typename is_segment_concept<
+          typename geometry_concept< Segment<long> >::type
+        >::type
+      >::type
+    >::type,
+    CT
+  >::type get_point_projection(
+      const Point<CT>& point, const Segment<InCT>& segment) {
+    CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
+    CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
+    CT point_vec_x = x(point) - cast(x(low(segment)));
+    CT point_vec_y = y(point) - cast(y(low(segment)));
+    CT sqr_segment_length =
+        segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y;
+    CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y;
+    return vec_dot / sqr_segment_length;
+  }
+
+  template <typename InCT>
+  static CT cast(const InCT& value) {
+    return static_cast<CT>(value);
+  }
+};
+
+} } // namespace boost::polygon
+#endif
+
 using namespace boost::polygon;  // provides also high() and low()
 
 namespace Slic3r { namespace Geometry {
@@ -290,6 +464,294 @@ arrange(size_t total_parts, Pointf part, coordf_t dist, const BoundingBoxf* bb)
     return positions;
 }
 
+#ifdef SLIC3R_DEBUG
+// The following code for the visualization of the boost Voronoi diagram is based on:
+//
+// Boost.Polygon library voronoi_visualizer.cpp file
+//          Copyright Andrii Sydorchuk 2010-2012.
+// Distributed under the Boost Software License, Version 1.0.
+//    (See accompanying file LICENSE_1_0.txt or copy at
+//          http://www.boost.org/LICENSE_1_0.txt)
+namespace Voronoi { namespace Internal {
+
+    typedef double coordinate_type;
+    typedef boost::polygon::point_data<coordinate_type> point_type;
+    typedef boost::polygon::segment_data<coordinate_type> segment_type;
+    typedef boost::polygon::rectangle_data<coordinate_type> rect_type;
+//    typedef voronoi_builder<int> VB;
+    typedef boost::polygon::voronoi_diagram<coordinate_type> VD;
+    typedef VD::cell_type cell_type;
+    typedef VD::cell_type::source_index_type source_index_type;
+    typedef VD::cell_type::source_category_type source_category_type;
+    typedef VD::edge_type edge_type;
+    typedef VD::cell_container_type cell_container_type;
+    typedef VD::cell_container_type vertex_container_type;
+    typedef VD::edge_container_type edge_container_type;
+    typedef VD::const_cell_iterator const_cell_iterator;
+    typedef VD::const_vertex_iterator const_vertex_iterator;
+    typedef VD::const_edge_iterator const_edge_iterator;
+
+    static const std::size_t EXTERNAL_COLOR = 1;
+
+    inline void color_exterior(const VD::edge_type* edge) 
+    {
+        if (edge->color() == EXTERNAL_COLOR)
+            return;
+        edge->color(EXTERNAL_COLOR);
+        edge->twin()->color(EXTERNAL_COLOR);
+        const VD::vertex_type* v = edge->vertex1();
+        if (v == NULL || !edge->is_primary())
+            return;
+        v->color(EXTERNAL_COLOR);
+        const VD::edge_type* e = v->incident_edge();
+        do {
+            color_exterior(e);
+            e = e->rot_next();
+        } while (e != v->incident_edge());
+    }
+
+    inline point_type retrieve_point(const std::vector<segment_type> &segments, const cell_type& cell) 
+    {
+        assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT);
+        return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]);
+    }
+
+    inline void clip_infinite_edge(const std::vector<segment_type> &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector<point_type>* clipped_edge) 
+    {
+        const cell_type& cell1 = *edge.cell();
+        const cell_type& cell2 = *edge.twin()->cell();
+        point_type origin, direction;
+        // Infinite edges could not be created by two segment sites.
+        if (cell1.contains_point() && cell2.contains_point()) {
+            point_type p1 = retrieve_point(segments, cell1);
+            point_type p2 = retrieve_point(segments, cell2);
+            origin.x((p1.x() + p2.x()) * 0.5);
+            origin.y((p1.y() + p2.y()) * 0.5);
+            direction.x(p1.y() - p2.y());
+            direction.y(p2.x() - p1.x());
+        } else {
+            origin = cell1.contains_segment() ? retrieve_point(segments, cell2) : retrieve_point(segments, cell1);
+            segment_type segment = cell1.contains_segment() ? segments[cell1.source_index()] : segments[cell2.source_index()];
+            coordinate_type dx = high(segment).x() - low(segment).x();
+            coordinate_type dy = high(segment).y() - low(segment).y();
+            if ((low(segment) == origin) ^ cell1.contains_point()) {
+                direction.x(dy);
+                direction.y(-dx);
+            } else {
+                direction.x(-dy);
+                direction.y(dx);
+            }
+        }
+        coordinate_type koef = bbox_max_size / (std::max)(fabs(direction.x()), fabs(direction.y()));
+        if (edge.vertex0() == NULL) {
+            clipped_edge->push_back(point_type(
+                origin.x() - direction.x() * koef,
+                origin.y() - direction.y() * koef));
+        } else {
+            clipped_edge->push_back(
+                point_type(edge.vertex0()->x(), edge.vertex0()->y()));
+        }
+        if (edge.vertex1() == NULL) {
+            clipped_edge->push_back(point_type(
+                origin.x() + direction.x() * koef,
+                origin.y() + direction.y() * koef));
+        } else {
+            clipped_edge->push_back(
+                point_type(edge.vertex1()->x(), edge.vertex1()->y()));
+        }
+    }
+
+    inline void sample_curved_edge(const std::vector<segment_type> &segments, const edge_type& edge, std::vector<point_type> &sampled_edge, coordinate_type max_dist) 
+    {
+        point_type point = edge.cell()->contains_point() ?
+            retrieve_point(segments, *edge.cell()) :
+            retrieve_point(segments, *edge.twin()->cell());
+        segment_type segment = edge.cell()->contains_point() ?
+            segments[edge.twin()->cell()->source_index()] :
+            segments[edge.cell()->source_index()];
+        ::boost::polygon::voronoi_visual_utils<coordinate_type>::discretize(point, segment, max_dist, &sampled_edge);
+    }
+
+} /* namespace Internal */ } // namespace Voronoi
+
+static inline void dump_voronoi_to_svg(const Lines &lines, /* const */ voronoi_diagram<double> &vd, const ThickPolylines *polylines, const char *path)
+{
+    const double        scale                       = 0.2;
+    const std::string   inputSegmentPointColor      = "lightseagreen";
+    const coord_t       inputSegmentPointRadius     = coord_t(0.09 * scale / SCALING_FACTOR); 
+    const std::string   inputSegmentColor           = "lightseagreen";
+    const coord_t       inputSegmentLineWidth       = coord_t(0.03 * scale / SCALING_FACTOR);
+
+    const std::string   voronoiPointColor           = "black";
+    const coord_t       voronoiPointRadius          = coord_t(0.06 * scale / SCALING_FACTOR);
+    const std::string   voronoiLineColorPrimary     = "black";
+    const std::string   voronoiLineColorSecondary   = "green";
+    const std::string   voronoiArcColor             = "red";
+    const coord_t       voronoiLineWidth            = coord_t(0.02 * scale / SCALING_FACTOR);
+
+    const bool          internalEdgesOnly           = false;
+    const bool          primaryEdgesOnly            = false;
+
+    BoundingBox bbox = BoundingBox(lines);
+    bbox.min.x -= coord_t(1. / SCALING_FACTOR);
+    bbox.min.y -= coord_t(1. / SCALING_FACTOR);
+    bbox.max.x += coord_t(1. / SCALING_FACTOR);
+    bbox.max.y += coord_t(1. / SCALING_FACTOR);
+
+    ::Slic3r::SVG svg(path, bbox);
+
+    if (polylines != NULL)
+        svg.draw(*polylines, "lime", "lime", voronoiLineWidth);
+
+//    bbox.scale(1.2);
+    // For clipping of half-lines to some reasonable value.
+    // The line will then be clipped by the SVG viewer anyway.
+    const double bbox_dim_max = double(bbox.max.x - bbox.min.x) + double(bbox.max.y - bbox.min.y);
+    // For the discretization of the Voronoi parabolic segments.
+    const double discretization_step = 0.0005 * bbox_dim_max;
+
+    // Make a copy of the input segments with the double type.
+    std::vector<Voronoi::Internal::segment_type> segments;
+    for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it)
+        segments.push_back(Voronoi::Internal::segment_type(
+            Voronoi::Internal::point_type(double(it->a.x), double(it->a.y)), 
+            Voronoi::Internal::point_type(double(it->b.x), double(it->b.y))));
+    
+    // Color exterior edges.
+    for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it)
+        if (!it->is_finite())
+            Voronoi::Internal::color_exterior(&(*it));
+
+    // Draw the end points of the input polygon.
+    for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) {
+        svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius);
+        svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius);
+    }
+    // Draw the input polygon.
+    for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it)
+        svg.draw(Line(Point(coord_t(it->a.x), coord_t(it->a.y)), Point(coord_t(it->b.x), coord_t(it->b.y))), inputSegmentColor, inputSegmentLineWidth);
+
+#if 1
+    // Draw voronoi vertices.
+    for (voronoi_diagram<double>::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it)
+        if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR)
+            svg.draw(Point(coord_t(it->x()), coord_t(it->y())), voronoiPointColor, voronoiPointRadius);
+
+    for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) {
+        if (primaryEdgesOnly && !it->is_primary())
+            continue;
+        if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR))
+            continue;
+        std::vector<Voronoi::Internal::point_type> samples;
+        std::string color = voronoiLineColorPrimary;
+        if (!it->is_finite()) {
+            Voronoi::Internal::clip_infinite_edge(segments, *it, bbox_dim_max, &samples);
+            if (! it->is_primary())
+                color = voronoiLineColorSecondary;
+        } else {
+            // Store both points of the segment into samples. sample_curved_edge will split the initial line
+            // until the discretization_step is reached.
+            samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y()));
+            samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y()));
+            if (it->is_curved()) {
+                Voronoi::Internal::sample_curved_edge(segments, *it, samples, discretization_step);
+                color = voronoiArcColor;
+            } else if (! it->is_primary())
+                color = voronoiLineColorSecondary;
+        }
+        for (std::size_t i = 0; i + 1 < samples.size(); ++i)
+            svg.draw(Line(Point(coord_t(samples[i].x()), coord_t(samples[i].y())), Point(coord_t(samples[i+1].x()), coord_t(samples[i+1].y()))), color, voronoiLineWidth);
+    }
+#endif
+
+    if (polylines != NULL)
+        svg.draw(*polylines, "blue", voronoiLineWidth);
+
+    svg.Close();
+}
+#endif /* SLIC3R_DEBUG */
+
+// Euclidian distance of two boost::polygon points.
+template<typename T>
+T dist(const boost::polygon::point_data<T> &p1,const boost::polygon::point_data<T> &p2)
+{
+	T dx = p2.x() - p1.x();
+	T dy = p2.y() - p1.y();
+	return sqrt(dx*dx+dy*dy);
+}
+
+// Find a foot point of "px" on a segment "seg".
+template<typename segment_type, typename point_type>
+inline point_type project_point_to_segment(segment_type &seg, point_type &px)
+{
+    typedef typename point_type::coordinate_type T;
+    const point_type &p0 = low(seg);
+    const point_type &p1 = high(seg);
+    const point_type  dir(p1.x()-p0.x(), p1.y()-p0.y());
+    const point_type  dproj(px.x()-p0.x(), px.y()-p0.y());
+    const T           t = (dir.x()*dproj.x() + dir.y()*dproj.y()) / (dir.x()*dir.x() + dir.y()*dir.y());
+    assert(t >= T(-1e-6) && t <= T(1. + 1e-6));
+    return point_type(p0.x() + t*dir.x(), p0.y() + t*dir.y());
+}
+
+template<typename VD, typename SEGMENTS>
+inline const typename VD::point_type retrieve_cell_point(const typename VD::cell_type& cell, const SEGMENTS &segments)
+{
+    assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT);
+    return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]);
+}
+
+template<typename VD, typename SEGMENTS>
+inline std::pair<typename VD::coord_type, typename VD::coord_type>
+measure_edge_thickness(const VD &vd, const typename VD::edge_type& edge, const SEGMENTS &segments)
+{
+	typedef typename VD::coord_type T;
+    const typename VD::point_type  pa(edge.vertex0()->x(), edge.vertex0()->y());
+    const typename VD::point_type  pb(edge.vertex1()->x(), edge.vertex1()->y());
+    const typename VD::cell_type  &cell1 = *edge.cell();
+    const typename VD::cell_type  &cell2 = *edge.twin()->cell();
+    if (cell1.contains_segment()) {
+        if (cell2.contains_segment()) {
+            // Both cells contain a linear segment, the left / right cells are symmetric.
+            // Project pa, pb to the left segment.
+            const typename VD::segment_type segment1 = segments[cell1.source_index()];
+            const typename VD::point_type p1a = project_point_to_segment(segment1, pa);
+            const typename VD::point_type p1b = project_point_to_segment(segment1, pb);
+            return std::pair<T, T>(T(2.)*dist(pa, p1a), T(2.)*dist(pb, p1b));
+        } else {
+            // 1st cell contains a linear segment, 2nd cell contains a point.
+            // The medial axis between the cells is a parabolic arc.
+            // Project pa, pb to the left segment.
+            const typename  VD::point_type p2 = retrieve_cell_point<VD>(cell2, segments);
+            return std::pair<T, T>(T(2.)*dist(pa, p2), T(2.)*dist(pb, p2));
+        }
+    } else if (cell2.contains_segment()) {
+        // 1st cell contains a point, 2nd cell contains a linear segment.
+        // The medial axis between the cells is a parabolic arc.
+        const typename VD::point_type p1 = retrieve_cell_point<VD>(cell1, segments);
+        return std::pair<T, T>(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1));
+    } else {
+        // Both cells contain a point. The left / right regions are triangular and symmetric.
+        const typename VD::point_type p1 = retrieve_cell_point<VD>(cell1, segments);
+        return std::pair<T, T>(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1));
+    }
+}
+
+// Converts the Line instances of Lines vector to VD::segment_type.
+template<typename VD>
+class Lines2VDSegments
+{
+public:
+    Lines2VDSegments(const Lines &alines) : lines(alines) {}
+    typename VD::segment_type operator[](size_t idx) const {
+        return typename VD::segment_type(
+            typename VD::point_type(typename VD::coord_type(lines[idx].a.x), typename VD::coord_type(lines[idx].a.y)),
+            typename VD::point_type(typename VD::coord_type(lines[idx].b.x), typename VD::coord_type(lines[idx].b.y)));
+    }
+private:
+    const Lines &lines;
+};
+
 void
 MedialAxis::build(ThickPolylines* polylines)
 {
@@ -373,6 +835,25 @@ MedialAxis::build(ThickPolylines* polylines)
         // append polyline to result
         polylines->push_back(polyline);
     }
+
+    #ifdef SLIC3R_DEBUG
+    {
+        char path[2048];
+        static int iRun = 0;
+        sprintf(path, "out/MedialAxis-%d.svg", iRun ++);
+        dump_voronoi_to_svg(this->lines, this->vd, polylines, path);
+
+
+        printf("Thick lines: ");
+        for (ThickPolylines::const_iterator it = polylines->begin(); it != polylines->end(); ++ it) {
+            ThickLines lines = it->thicklines();
+            for (ThickLines::const_iterator it2 = lines.begin(); it2 != lines.end(); ++ it2) {
+                printf("%f,%f ", it2->a_width, it2->b_width);
+            }
+        }
+        printf("\n");
+    }
+    #endif /* SLIC3R_DEBUG */
 }
 
 void
diff --git a/xs/src/libslic3r/Geometry.hpp b/xs/src/libslic3r/Geometry.hpp
index 14bc3f0ca8..75762a8c87 100644
--- a/xs/src/libslic3r/Geometry.hpp
+++ b/xs/src/libslic3r/Geometry.hpp
@@ -52,7 +52,13 @@ class MedialAxis {
     void build(Polylines* polylines);
     
     private:
-    typedef voronoi_diagram<double> VD;
+    class VD : public voronoi_diagram<double> {
+    public:
+        typedef double                                          coord_type;
+        typedef boost::polygon::point_data<coordinate_type>     point_type;
+        typedef boost::polygon::segment_data<coordinate_type>   segment_type;
+        typedef boost::polygon::rectangle_data<coordinate_type> rect_type;
+    };
     VD vd;
     std::set<const VD::edge_type*> edges, valid_edges;
     std::map<const VD::edge_type*, std::pair<coordf_t,coordf_t> > thickness;
diff --git a/xs/src/libslic3r/LayerRegion.cpp b/xs/src/libslic3r/LayerRegion.cpp
index f459c58788..50f6ebbea4 100644
--- a/xs/src/libslic3r/LayerRegion.cpp
+++ b/xs/src/libslic3r/LayerRegion.cpp
@@ -112,7 +112,7 @@ LayerRegion::process_external_surfaces(const Layer* lower_layer)
             );
             
             #ifdef SLIC3R_DEBUG
-            printf("Processing bridge at layer %zu:\n", this->layer()->id();
+            printf("Processing bridge at layer %zu:\n", this->layer()->id());
             #endif
             
             if (bd.detect_angle()) {
diff --git a/xs/src/libslic3r/MultiPoint.cpp b/xs/src/libslic3r/MultiPoint.cpp
index 6857d6393e..0b7f8f3c21 100644
--- a/xs/src/libslic3r/MultiPoint.cpp
+++ b/xs/src/libslic3r/MultiPoint.cpp
@@ -30,6 +30,19 @@ MultiPoint::translate(const Point &vector)
     this->translate(vector.x, vector.y);
 }
 
+void
+MultiPoint::rotate(double angle)
+{
+    double s     = sin(angle);
+    double c     = cos(angle);
+    for (Points::iterator it = points.begin(); it != points.end(); ++it) {
+	    double cur_x = (double)it->x;
+	    double cur_y = (double)it->y;
+	    it->x = (coord_t)round(c * cur_x - s * cur_y);
+	    it->y = (coord_t)round(c * cur_y + s * cur_x);
+    }
+}
+
 void
 MultiPoint::rotate(double angle, const Point &center)
 {
@@ -89,15 +102,33 @@ MultiPoint::bounding_box() const
     return BoundingBox(this->points);
 }
 
-void
+bool 
+MultiPoint::has_duplicate_points() const
+{
+    for (size_t i = 1; i < points.size(); ++i)
+        if (points[i-1].coincides_with(points[i]))
+            return true;
+    return false;
+}
+
+bool
 MultiPoint::remove_duplicate_points()
 {
-    for (size_t i = 1; i < this->points.size(); ++i) {
-        if (this->points.at(i).coincides_with(this->points.at(i-1))) {
-            this->points.erase(this->points.begin() + i);
-            --i;
+    size_t j = 0;
+    for (size_t i = 1; i < points.size(); ++i) {
+        if (points[j].coincides_with(points[i])) {
+            // Just increase index i.
+        } else {
+            ++ j;
+            if (j < i)
+                points[j] = points[i];
         }
     }
+    if (++ j < points.size()) {
+        points.erase(points.begin() + j, points.end());
+        return true;
+    }
+    return false;
 }
 
 void
diff --git a/xs/src/libslic3r/MultiPoint.hpp b/xs/src/libslic3r/MultiPoint.hpp
index d057b14f12..08d2b86701 100644
--- a/xs/src/libslic3r/MultiPoint.hpp
+++ b/xs/src/libslic3r/MultiPoint.hpp
@@ -22,6 +22,7 @@ class MultiPoint
     void scale(double factor);
     void translate(double x, double y);
     void translate(const Point &vector);
+    void rotate(double angle);
     void rotate(double angle, const Point &center);
     void reverse();
     Point first_point() const;
@@ -32,7 +33,10 @@ class MultiPoint
     int find_point(const Point &point) const;
     bool has_boundary_point(const Point &point) const;
     BoundingBox bounding_box() const;
-    void remove_duplicate_points();
+    // Return true if there are exact duplicates.
+    bool has_duplicate_points() const;
+    // Remove exact duplicates, return true if any duplicate has been removed.
+    bool remove_duplicate_points();
     void append(const Point &point);
     void append(const Points &points);
     void append(const Points::const_iterator &begin, const Points::const_iterator &end);
diff --git a/xs/src/libslic3r/Point.cpp b/xs/src/libslic3r/Point.cpp
index 29ce0b1532..d7819bce28 100644
--- a/xs/src/libslic3r/Point.cpp
+++ b/xs/src/libslic3r/Point.cpp
@@ -54,13 +54,26 @@ Point::translate(const Vector &vector)
     this->translate(vector.x, vector.y);
 }
 
+void
+Point::rotate(double angle)
+{
+    double cur_x = (double)this->x;
+    double cur_y = (double)this->y;
+    double s     = sin(angle);
+    double c     = cos(angle);
+    this->x = (coord_t)round(c * cur_x - s * cur_y);
+    this->y = (coord_t)round(c * cur_y + s * cur_x);
+}
+
 void
 Point::rotate(double angle, const Point &center)
 {
     double cur_x = (double)this->x;
     double cur_y = (double)this->y;
-    this->x = (coord_t)round( (double)center.x + cos(angle) * (cur_x - (double)center.x) - sin(angle) * (cur_y - (double)center.y) );
-    this->y = (coord_t)round( (double)center.y + cos(angle) * (cur_y - (double)center.y) + sin(angle) * (cur_x - (double)center.x) );
+    double s     = sin(angle);
+    double c     = cos(angle);
+    this->x = (coord_t)round( (double)center.x + c * (cur_x - (double)center.x) - s * (cur_y - (double)center.y) );
+    this->y = (coord_t)round( (double)center.y + c * (cur_y - (double)center.y) + s * (cur_x - (double)center.x) );
 }
 
 bool
@@ -301,6 +314,12 @@ operator+(const Point& point1, const Point& point2)
     return Point(point1.x + point2.x, point1.y + point2.y);
 }
 
+Point
+operator-(const Point& point1, const Point& point2)
+{
+    return Point(point1.x - point2.x, point1.y - point2.y);
+}
+
 Point
 operator*(double scalar, const Point& point2)
 {
@@ -349,13 +368,26 @@ Pointf::translate(const Vectorf &vector)
     this->translate(vector.x, vector.y);
 }
 
+void
+Pointf::rotate(double angle)
+{
+    double cur_x = this->x;
+    double cur_y = this->y;
+    double s     = sin(angle);
+    double c     = cos(angle);
+    this->x = c * cur_x - s * cur_y;
+    this->y = c * cur_y + s * cur_x;
+}
+
 void
 Pointf::rotate(double angle, const Pointf &center)
 {
     double cur_x = this->x;
     double cur_y = this->y;
-    this->x = center.x + cos(angle) * (cur_x - center.x) - sin(angle) * (cur_y - center.y);
-    this->y = center.y + cos(angle) * (cur_y - center.y) + sin(angle) * (cur_x - center.x);
+    double s     = sin(angle);
+    double c     = cos(angle);
+    this->x = center.x + c * (cur_x - center.x) - s * (cur_y - center.y);
+    this->y = center.y + c * (cur_y - center.y) + s * (cur_x - center.x);
 }
 
 Pointf
diff --git a/xs/src/libslic3r/Point.hpp b/xs/src/libslic3r/Point.hpp
index 6638e0eb3f..45331aa714 100644
--- a/xs/src/libslic3r/Point.hpp
+++ b/xs/src/libslic3r/Point.hpp
@@ -42,6 +42,7 @@ class Point
     void scale(double factor);
     void translate(double x, double y);
     void translate(const Vector &vector);
+    void rotate(double angle);
     void rotate(double angle, const Point &center);
     bool coincides_with(const Point &point) const;
     bool coincides_with_epsilon(const Point &point) const;
@@ -64,6 +65,7 @@ class Point
 };
 
 Point operator+(const Point& point1, const Point& point2);
+Point operator-(const Point& point1, const Point& point2);
 Point operator*(double scalar, const Point& point2);
 
 class Point3 : public Point
@@ -92,6 +94,7 @@ class Pointf
     void scale(double factor);
     void translate(double x, double y);
     void translate(const Vectorf &vector);
+    void rotate(double angle);
     void rotate(double angle, const Pointf &center);
     Pointf negative() const;
     Vectorf vector_to(const Pointf &point) const;
diff --git a/xs/src/libslic3r/PolylineCollection.cpp b/xs/src/libslic3r/PolylineCollection.cpp
index f3077b0e1e..fcf95e6c1f 100644
--- a/xs/src/libslic3r/PolylineCollection.cpp
+++ b/xs/src/libslic3r/PolylineCollection.cpp
@@ -2,50 +2,122 @@
 
 namespace Slic3r {
 
-void
-PolylineCollection::chained_path(PolylineCollection* retval, bool no_reverse) const
+struct Chaining
 {
-    if (this->polylines.empty()) return;
-    this->chained_path_from(this->polylines.front().first_point(), retval, no_reverse);
-}
+    Point first;
+    Point last;
+    size_t idx;
+};
 
-void
-PolylineCollection::chained_path_from(Point start_near, PolylineCollection* retval, bool no_reverse) const
+#ifndef sqr
+template<typename T>
+inline T sqr(T x) { return x * x; }
+#endif /* sqr */
+
+template<typename T>
+inline int nearest_point_index(const std::vector<Chaining> &pairs, const Point &start_near, bool no_reverse)
 {
-    Polylines my_paths = this->polylines;
-    
-    Points endpoints;
-    for (Polylines::const_iterator it = my_paths.begin(); it != my_paths.end(); ++it) {
-        endpoints.push_back(it->first_point());
-        if (no_reverse) {
-            endpoints.push_back(it->first_point());
-        } else {
-            endpoints.push_back(it->last_point());
+    T dmin = std::numeric_limits<T>::max();
+    int idx = 0;
+    for (std::vector<Chaining>::const_iterator it = pairs.begin(); it != pairs.end(); ++it) {
+        T d = sqr(T(start_near.x - it->first.x));
+        if (d <= dmin) {
+            d += sqr(T(start_near.y - it->first.y));
+            if (d < dmin) {
+                idx = (it - pairs.begin()) * 2;
+                dmin = d;
+                if (dmin < EPSILON)
+                    break;
+            }
+        }
+        if (! no_reverse) {
+            d = sqr(T(start_near.x - it->last.x));
+            if (d <= dmin) {
+                d += sqr(T(start_near.y - it->last.y));
+                if (d < dmin) {
+                    idx = (it - pairs.begin()) * 2 + 1;
+                    dmin = d;
+                    if (dmin < EPSILON)
+                        break;
+                }
+            }
         }
     }
-    
-    while (!my_paths.empty()) {
+    return idx;
+}
+
+Polylines PolylineCollection::chained_path_from(
+#if SLIC3R_CPPVER > 11
+    Polylines &&src,
+#else
+    const Polylines &src,
+#endif
+    Point start_near, 
+    bool  no_reverse)
+{
+    std::vector<Chaining> endpoints;
+    endpoints.reserve(src.size());
+    for (size_t i = 0; i < src.size(); ++ i) {
+        Chaining c;
+        c.first = src[i].first_point();
+        if (! no_reverse)
+            c.last = src[i].last_point();
+        c.idx = i;
+        endpoints.push_back(c);
+    }
+    Polylines retval;
+    while (! endpoints.empty()) {
         // find nearest point
-        int start_index = start_near.nearest_point_index(endpoints);
-        int path_index = start_index/2;
-        if (start_index % 2 && !no_reverse) {
-            my_paths.at(path_index).reverse();
-        }
-        retval->polylines.push_back(my_paths.at(path_index));
-        my_paths.erase(my_paths.begin() + path_index);
-        endpoints.erase(endpoints.begin() + 2*path_index, endpoints.begin() + 2*path_index + 2);
-        start_near = retval->polylines.back().last_point();
+        int endpoint_index = nearest_point_index<double>(endpoints, start_near, no_reverse);
+        assert(endpoint_index >= 0 && endpoint_index < endpoints.size() * 2);
+#if SLIC3R_CPPVER > 11
+        retval.push_back(std::move(src[endpoints[endpoint_index/2].idx]));
+#else
+        retval.push_back(src[endpoints[endpoint_index/2].idx]);
+#endif
+        if (endpoint_index & 1)
+            retval.back().reverse();
+        endpoints.erase(endpoints.begin() + endpoint_index/2);
+        start_near = retval.back().last_point();
     }
+    return retval;
 }
 
-Point
-PolylineCollection::leftmost_point() const
+#if SLIC3R_CPPVER > 11
+Polylines PolylineCollection::chained_path(Polylines &&src, bool no_reverse)
 {
-    if (this->polylines.empty()) CONFESS("leftmost_point() called on empty PolylineCollection");
-    Point p = this->polylines.front().leftmost_point();
-    for (Polylines::const_iterator it = this->polylines.begin() + 1; it != this->polylines.end(); ++it) {
+    return (src.empty() || src.front().empty()) ?
+        Polylines() :
+        chained_path_from(std::move(src), src.front().first_point(), no_reverse);
+}
+Polylines PolylineCollection::chained_path_from(Polylines src, Point start_near, bool no_reverse)
+{
+    return chained_path_from(std::move(src), start_near, no_reverse);
+}
+Polylines PolylineCollection::chained_path(Polylines src, bool no_reverse)
+{
+    return (src.empty() || src.front().empty()) ?
+        Polylines() :
+        chained_path_from(std::move(src), src.front().first_point(), no_reverse);
+}
+#else
+Polylines PolylineCollection::chained_path(const Polylines &src, bool no_reverse)
+{
+    return (src.empty() || src.front().points.empty()) ?
+        Polylines() :
+        chained_path_from(src, src.front().first_point(), no_reverse);
+}
+#endif
+
+Point PolylineCollection::leftmost_point(const Polylines &polylines)
+{
+    if (polylines.empty()) CONFESS("leftmost_point() called on empty PolylineCollection");
+    Polylines::const_iterator it = polylines.begin();
+    Point p = it->leftmost_point();
+    for (++ it; it != polylines.end(); ++it) {
         Point p2 = it->leftmost_point();
-        if (p2.x < p.x) p = p2;
+        if (p2.x < p.x) 
+            p = p2;
     }
     return p;
 }
@@ -56,4 +128,4 @@ PolylineCollection::append(const Polylines &pp)
     this->polylines.insert(this->polylines.end(), pp.begin(), pp.end());
 }
 
-}
+} // namespace Slic3r
diff --git a/xs/src/libslic3r/PolylineCollection.hpp b/xs/src/libslic3r/PolylineCollection.hpp
index d41e8c465c..492f0318fb 100644
--- a/xs/src/libslic3r/PolylineCollection.hpp
+++ b/xs/src/libslic3r/PolylineCollection.hpp
@@ -8,12 +8,26 @@ namespace Slic3r {
 
 class PolylineCollection
 {
-    public:
+public:
     Polylines polylines;
-    void chained_path(PolylineCollection* retval, bool no_reverse = false) const;
-    void chained_path_from(Point start_near, PolylineCollection* retval, bool no_reverse = false) const;
-    Point leftmost_point() const;
+    void chained_path(PolylineCollection* retval, bool no_reverse = false) const
+    	{ retval->polylines = chained_path(this->polylines, no_reverse); }
+    void chained_path_from(Point start_near, PolylineCollection* retval, bool no_reverse = false) const 
+    	{ retval->polylines = chained_path_from(this->polylines, start_near, no_reverse); }
+    Point leftmost_point() const
+    	{ return leftmost_point(polylines); }
     void append(const Polylines &polylines);
+
+	static Point     leftmost_point(const Polylines &polylines);
+#if SLIC3R_CPPVER > 11
+	static Polylines chained_path(Polylines &&src, bool no_reverse = false);
+	static Polylines chained_path_from(Polylines &&src, Point start_near, bool no_reverse = false);
+	static Polylines chained_path(Polylines src, bool no_reverse = false);
+	static Polylines chained_path_from(Polylines src, Point start_near, bool no_reverse = false);
+#else
+	static Polylines chained_path(const Polylines &src, bool no_reverse = false);
+	static Polylines chained_path_from(const Polylines &src, Point start_near, bool no_reverse = false);
+#endif
 };
 
 }
diff --git a/xs/src/libslic3r/SVG.cpp b/xs/src/libslic3r/SVG.cpp
index af367ae76e..259c6eacfc 100644
--- a/xs/src/libslic3r/SVG.cpp
+++ b/xs/src/libslic3r/SVG.cpp
@@ -19,18 +19,54 @@ SVG::SVG(const char* filename)
 	    );
 }
 
+SVG::SVG(const char* filename, const BoundingBox &bbox)
+    : arrows(false), fill("grey"), stroke("black"), filename(filename), origin(bbox.min)
+{
+    this->f = fopen(filename, "w");
+    float w = COORD(bbox.max.x - bbox.min.x);
+    float h = COORD(bbox.max.y - bbox.min.y);
+    fprintf(this->f,
+        "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"yes\"?>\n"
+        "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.0//EN\" \"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n"
+        "<svg height=\"%f\" width=\"%f\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:svg=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n"
+        "   <marker id=\"endArrow\" markerHeight=\"8\" markerUnits=\"strokeWidth\" markerWidth=\"10\" orient=\"auto\" refX=\"1\" refY=\"5\" viewBox=\"0 0 10 10\">\n"
+        "      <polyline fill=\"darkblue\" points=\"0,0 10,5 0,10 1,5\" />\n"
+        "   </marker>\n",
+        h, w);
+}
+
 void
-SVG::draw(const Line &line, std::string stroke)
+SVG::draw(const Line &line, std::string stroke, coord_t stroke_width)
 {
     fprintf(this->f,
-        "   <line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke: %s; stroke-width: 1\"",
-        COORD(line.a.x), COORD(line.a.y), COORD(line.b.x), COORD(line.b.y), stroke.c_str()
-        );
+        "   <line x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\" style=\"stroke: %s; stroke-width: %f\"",
+        COORD(line.a.x - origin.x), COORD(line.a.y - origin.y), COORD(line.b.x - origin.x), COORD(line.b.y - origin.y), stroke.c_str(), (stroke_width == 0) ? 1.f : COORD(stroke_width));
     if (this->arrows)
         fprintf(this->f, " marker-end=\"url(#endArrow)\"");
     fprintf(this->f, "/>\n");
 }
 
+void SVG::draw(const ThickLine &line, const std::string &fill, const std::string &stroke, coord_t stroke_width)
+{
+    Pointf dir(line.b.x-line.a.x, line.b.y-line.a.y);
+    Pointf perp(-dir.y, dir.x);
+    coordf_t len = sqrt(perp.x*perp.x + perp.y*perp.y);
+    coordf_t da  = coordf_t(0.5)*line.a_width/len;
+    coordf_t db  = coordf_t(0.5)*line.b_width/len;
+    fprintf(this->f,
+        "   <polygon points=\"%f,%f %f,%f %f,%f %f,%f\" style=\"fill:%s; stroke: %s; stroke-width: %f\"/>\n",
+        COORD(line.a.x-da*perp.x-origin.x),
+        COORD(line.a.y-da*perp.y-origin.y),
+        COORD(line.b.x-db*perp.x-origin.x),
+        COORD(line.b.y-db*perp.y-origin.y),
+        COORD(line.b.x+db*perp.x-origin.x),
+        COORD(line.b.y+db*perp.y-origin.y),
+        COORD(line.a.x+da*perp.x-origin.x),
+        COORD(line.a.y+da*perp.y-origin.y),
+        fill.c_str(), stroke.c_str(),
+        (stroke_width == 0) ? 1.f : COORD(stroke_width));
+}
+
 void
 SVG::draw(const Lines &lines, std::string stroke)
 {
@@ -80,31 +116,45 @@ SVG::draw(const Polygons &polygons, std::string fill)
 }
 
 void
-SVG::draw(const Polyline &polyline, std::string stroke)
+SVG::draw(const Polyline &polyline, std::string stroke, coord_t stroke_width)
 {
     this->stroke = stroke;
-    this->path(this->get_path_d(polyline, false), false);
+    this->path(this->get_path_d(polyline, false), false, stroke_width);
 }
 
 void
-SVG::draw(const Polylines &polylines, std::string stroke)
+SVG::draw(const Polylines &polylines, std::string stroke, coord_t stroke_width)
 {
     for (Polylines::const_iterator it = polylines.begin(); it != polylines.end(); ++it)
-        this->draw(*it, stroke);
+        this->draw(*it, fill, stroke_width);
+}
+
+void SVG::draw(const ThickLines &thicklines, const std::string &fill, const std::string &stroke, coord_t stroke_width)
+{
+    for (ThickLines::const_iterator it = thicklines.begin(); it != thicklines.end(); ++it)
+        this->draw(*it, fill, stroke, stroke_width);
 }
 
 void
-SVG::draw(const ThickPolylines &polylines, std::string stroke)
+SVG::draw(const ThickPolylines &polylines, const std::string &stroke, coord_t stroke_width)
 {
     for (ThickPolylines::const_iterator it = polylines.begin(); it != polylines.end(); ++it)
-        this->draw((Polyline)*it, stroke);
+        this->draw((Polyline)*it, stroke, stroke_width);
+}
+
+void 
+SVG::draw(const ThickPolylines &thickpolylines, const std::string &fill, const std::string &stroke, coord_t stroke_width)
+{
+    for (ThickPolylines::const_iterator it = thickpolylines.begin(); it != thickpolylines.end(); ++ it)
+        draw(it->thicklines(), fill, stroke, stroke_width);
 }
 
 void
-SVG::draw(const Point &point, std::string fill, unsigned int radius)
+SVG::draw(const Point &point, std::string fill, coord_t iradius)
 {
+    float radius = (iradius == 0) ? 3.f : COORD(iradius);
     std::ostringstream svg;
-    svg << "   <circle cx=\"" << COORD(point.x) << "\" cy=\"" << COORD(point.y)
+    svg << "   <circle cx=\"" << COORD(point.x - origin.x) << "\" cy=\"" << COORD(point.y - origin.y)
         << "\" r=\"" << radius << "\" "
         << "style=\"stroke: none; fill: " << fill << "\" />";
     
@@ -112,22 +162,26 @@ SVG::draw(const Point &point, std::string fill, unsigned int radius)
 }
 
 void
-SVG::draw(const Points &points, std::string fill, unsigned int radius)
+SVG::draw(const Points &points, std::string fill, coord_t radius)
 {
     for (Points::const_iterator it = points.begin(); it != points.end(); ++it)
         this->draw(*it, fill, radius);
 }
 
 void
-SVG::path(const std::string &d, bool fill)
+SVG::path(const std::string &d, bool fill, coord_t stroke_width)
 {
+    float lineWidth = 0.f;
+    if (! fill)
+        lineWidth = (stroke_width == 0) ? 2.f : COORD(stroke_width);
+
     fprintf(
         this->f,
-        "   <path d=\"%s\" style=\"fill: %s; stroke: %s; stroke-width: %s; fill-type: evenodd\" %s />\n",
+        "   <path d=\"%s\" style=\"fill: %s; stroke: %s; stroke-width: %f; fill-type: evenodd\" %s />\n",
         d.c_str(),
         fill ? this->fill.c_str() : "none",
         this->stroke.c_str(),
-        fill ? "0" : "2",
+        lineWidth, 
         (this->arrows && !fill) ? " marker-end=\"url(#endArrow)\"" : ""
     );
 }
@@ -138,8 +192,8 @@ SVG::get_path_d(const MultiPoint &mp, bool closed) const
     std::ostringstream d;
     d << "M ";
     for (Points::const_iterator p = mp.points.begin(); p != mp.points.end(); ++p) {
-        d << COORD(p->x) << " ";
-        d << COORD(p->y) << " ";
+        d << COORD(p->x - origin.x) << " ";
+        d << COORD(p->y - origin.y) << " ";
     }
     if (closed) d << "z";
     return d.str();
@@ -150,6 +204,7 @@ SVG::Close()
 {
     fprintf(this->f, "</svg>\n");
     fclose(this->f);
+    this->f = NULL;
     printf("SVG written to %s\n", this->filename.c_str());
 }
 
diff --git a/xs/src/libslic3r/SVG.hpp b/xs/src/libslic3r/SVG.hpp
index f0ae06d042..aa08c0e203 100644
--- a/xs/src/libslic3r/SVG.hpp
+++ b/xs/src/libslic3r/SVG.hpp
@@ -13,27 +13,34 @@ class SVG
     public:
     bool arrows;
     std::string fill, stroke;
-    
+    Point origin;
+
     SVG(const char* filename);
-    void draw(const Line &line, std::string stroke = "black");
+    SVG(const char* filename, const BoundingBox &bbox);
+    ~SVG() { if (f != NULL) Close(); }
+
+    void draw(const Line &line, std::string stroke = "black", coord_t stroke_width = 0);
+    void draw(const ThickLine &line, const std::string &fill, const std::string &stroke, coord_t stroke_width = 0);
     void draw(const Lines &lines, std::string stroke = "black");
     void draw(const IntersectionLines &lines, std::string stroke = "black");
     void draw(const ExPolygon &expolygon, std::string fill = "grey");
     void draw(const ExPolygons &expolygons, std::string fill = "grey");
     void draw(const Polygon &polygon, std::string fill = "grey");
     void draw(const Polygons &polygons, std::string fill = "grey");
-    void draw(const Polyline &polyline, std::string stroke = "black");
-    void draw(const Polylines &polylines, std::string stroke = "black");
-    void draw(const ThickPolylines &polylines, std::string stroke = "black");
-    void draw(const Point &point, std::string fill = "black", unsigned int radius = 3);
-    void draw(const Points &points, std::string fill = "black", unsigned int radius = 3);
+    void draw(const Polyline &polyline, std::string stroke = "black", coord_t stroke_width = 0);
+    void draw(const Polylines &polylines, std::string stroke = "black", coord_t stroke_width = 0);
+    void draw(const ThickLines &thicklines, const std::string &fill = "lime", const std::string &stroke = "black", coord_t stroke_width = 0);
+    void draw(const ThickPolylines &polylines, const std::string &stroke = "black", coord_t stroke_width = 0);
+    void draw(const ThickPolylines &thickpolylines, const std::string &fill, const std::string &stroke, coord_t stroke_width);
+    void draw(const Point &point, std::string fill = "black", coord_t radius = 0);
+    void draw(const Points &points, std::string fill = "black", coord_t radius = 0);
     void Close();
     
     private:
     std::string filename;
     FILE* f;
     
-    void path(const std::string &d, bool fill);
+    void path(const std::string &d, bool fill, coord_t stroke_width = 0);
     std::string get_path_d(const MultiPoint &mp, bool closed = false) const;
 };
 
diff --git a/xs/src/libslic3r/TriangleMesh.cpp b/xs/src/libslic3r/TriangleMesh.cpp
index 484fee9fc1..91e7982aa8 100644
--- a/xs/src/libslic3r/TriangleMesh.cpp
+++ b/xs/src/libslic3r/TriangleMesh.cpp
@@ -13,6 +13,7 @@
 #include <assert.h>
 
 #ifdef SLIC3R_DEBUG
+// #define SLIC3R_TRIANGLEMESH_DEBUG
 #include "SVG.hpp"
 #endif
 
@@ -445,7 +446,7 @@ TriangleMeshSlicer::slice(const std::vector<float> &z, std::vector<Polygons>* la
         float min_z = fminf(facet->vertex[0].z, fminf(facet->vertex[1].z, facet->vertex[2].z));
         float max_z = fmaxf(facet->vertex[0].z, fmaxf(facet->vertex[1].z, facet->vertex[2].z));
         
-        #ifdef SLIC3R_DEBUG
+        #ifdef SLIC3R_TRIANGLEMESH_DEBUG
         printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx,
             facet->vertex[0].x, facet->vertex[0].y, facet->vertex[0].z,
             facet->vertex[1].x, facet->vertex[1].y, facet->vertex[1].z,
@@ -457,7 +458,7 @@ TriangleMeshSlicer::slice(const std::vector<float> &z, std::vector<Polygons>* la
         std::vector<float>::const_iterator min_layer, max_layer;
         min_layer = std::lower_bound(z.begin(), z.end(), min_z); // first layer whose slice_z is >= min_z
         max_layer = std::upper_bound(z.begin() + (min_layer - z.begin()), z.end(), max_z) - 1; // last layer whose slice_z is <= max_z
-        #ifdef SLIC3R_DEBUG
+        #ifdef SLIC3R_TRIANGLEMESH_DEBUG
         printf("layers: min = %d, max = %d\n", (int)(min_layer - z.begin()), (int)(max_layer - z.begin()));
         #endif
         
@@ -473,7 +474,7 @@ TriangleMeshSlicer::slice(const std::vector<float> &z, std::vector<Polygons>* la
     layers->resize(z.size());
     for (std::vector<IntersectionLines>::iterator it = lines.begin(); it != lines.end(); ++it) {
         size_t layer_idx = it - lines.begin();
-        #ifdef SLIC3R_DEBUG
+        #ifdef SLIC3R_TRIANGLEMESH_DEBUG
         printf("Layer %zu:\n", layer_idx);
         #endif
         this->make_loops(*it, &(*layers)[layer_idx]);
@@ -488,7 +489,7 @@ TriangleMeshSlicer::slice(const std::vector<float> &z, std::vector<ExPolygons>*
     
     layers->resize(z.size());
     for (std::vector<Polygons>::const_iterator loops = layers_p.begin(); loops != layers_p.end(); ++loops) {
-        #ifdef SLIC3R_DEBUG
+        #ifdef SLIC3R_TRIANGLEMESH_DEBUG
         size_t layer_id = loops - layers_p.begin();
         printf("Layer %zu (slice_z = %.2f):\n", layer_id, z[layer_id]);
         #endif
@@ -712,7 +713,7 @@ TriangleMeshSlicer::make_loops(std::vector<IntersectionLine> &lines, Polygons* l
                     }
                     loops->push_back(p);
                     
-                    #ifdef SLIC3R_DEBUG
+                    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
                     printf("  Discovered %s polygon of %d points\n", (p.is_counter_clockwise() ? "ccw" : "cw"), (int)p.points.size());
                     #endif
                     
@@ -721,7 +722,7 @@ TriangleMeshSlicer::make_loops(std::vector<IntersectionLine> &lines, Polygons* l
                 
                 // we can't close this loop!
                 //// push @failed_loops, [@loop];
-                //#ifdef SLIC3R_DEBUG
+                //#ifdef SLIC3R_TRIANGLEMESH_DEBUG
                 printf("  Unable to close this loop having %d points\n", (int)loop.size());
                 //#endif
                 goto CYCLE;
@@ -833,7 +834,7 @@ TriangleMeshSlicer::make_expolygons(const Polygons &loops, ExPolygons* slices)
     ExPolygons ex_slices;
     offset2(p_slices, &ex_slices, +safety_offset, -safety_offset);
     
-    #ifdef SLIC3R_DEBUG
+    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
     size_t holes_count = 0;
     for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++e) {
         holes_count += e->holes.size();
@@ -1052,7 +1053,7 @@ TriangleMeshSlicer::TriangleMeshSlicer(TriangleMesh* _mesh) : mesh(_mesh), v_sca
                 }
                 this->facets_edges[facet_idx][i] = edge_idx;
                 
-                #ifdef SLIC3R_DEBUG
+                #ifdef SLIC3R_TRIANGLEMESH_DEBUG
                 printf("  [facet %d, edge %d] a_id = %d, b_id = %d   --> edge %d\n", facet_idx, i, a_id, b_id, edge_idx);
                 #endif
             }
diff --git a/xs/src/libslic3r/libslic3r.h b/xs/src/libslic3r/libslic3r.h
index 168e3f49cd..a087455cff 100644
--- a/xs/src/libslic3r/libslic3r.h
+++ b/xs/src/libslic3r/libslic3r.h
@@ -39,4 +39,14 @@ using namespace Slic3r;
 void confess_at(const char *file, int line, const char *func, const char *pat, ...);
 /* End implementation of CONFESS("foo"): */
 
+// Which C++ version is supported?
+// For example, could optimized functions with move semantics be used?
+#if __cplusplus==201402L
+	#define SLIC3R_CPPVER 14
+#elif __cplusplus==201103L
+	#define SLIC3R_CPPVER 11
+#else
+	#define SLIC3R_CPPVER 0
+#endif
+
 #endif
diff --git a/xs/src/perlglue.cpp b/xs/src/perlglue.cpp
index 55e00bc3a6..d3f6c48e50 100644
--- a/xs/src/perlglue.cpp
+++ b/xs/src/perlglue.cpp
@@ -10,6 +10,8 @@ REGISTER_CLASS(ExtrusionPath, "ExtrusionPath");
 REGISTER_CLASS(ExtrusionLoop, "ExtrusionLoop");
 // there is no ExtrusionLoop::Collection or ExtrusionEntity::Collection
 REGISTER_CLASS(ExtrusionEntityCollection, "ExtrusionPath::Collection");
+REGISTER_CLASS(ExtrusionSimulator, "ExtrusionSimulator");
+REGISTER_CLASS(Filler, "Filler");
 REGISTER_CLASS(Flow, "Flow");
 REGISTER_CLASS(AvoidCrossingPerimeters, "GCode::AvoidCrossingPerimeters");
 REGISTER_CLASS(OozePrevention, "GCode::OozePrevention");
diff --git a/xs/src/slic3r/GUI/3DScene.cpp b/xs/src/slic3r/GUI/3DScene.cpp
index e398a1c37c..323b4240b3 100644
--- a/xs/src/slic3r/GUI/3DScene.cpp
+++ b/xs/src/slic3r/GUI/3DScene.cpp
@@ -1,5 +1,7 @@
 #include "3DScene.hpp"
 
+#include <assert.h>
+
 namespace Slic3r {
 
 // caller is responsible for supplying NO lines with zero length
@@ -15,6 +17,10 @@ _3DScene::_extrusionentity_to_verts_do(const Lines &lines, const std::vector<dou
     // two triangles for each corner
     tverts->reserve_more(3 * 3 * 2 * (lines.size() + 1));
     */
+
+    assert(! lines.empty());
+    if (lines.empty())
+        return;
     
     Line prev_line;
     Pointf prev_b1, prev_b2;
diff --git a/xs/src/slic3r/GUI/GUI.cpp b/xs/src/slic3r/GUI/GUI.cpp
index fd5135f0a8..934d07ebd7 100644
--- a/xs/src/slic3r/GUI/GUI.cpp
+++ b/xs/src/slic3r/GUI/GUI.cpp
@@ -36,4 +36,23 @@ enable_screensaver()
     #endif
 }
 
+bool 
+debugged()
+{
+    #ifdef _WIN32
+    return IsDebuggerPresent();
+	#else
+	return false;
+    #endif /* _WIN32 */
+}
+
+void
+break_to_debugger()
+{
+    #ifdef _WIN32
+    if (IsDebuggerPresent())
+        DebugBreak();
+    #endif /* _WIN32 */
+}
+
 } }
diff --git a/xs/src/slic3r/GUI/GUI.hpp b/xs/src/slic3r/GUI/GUI.hpp
index 3ce72fd61b..85656c0a3a 100644
--- a/xs/src/slic3r/GUI/GUI.hpp
+++ b/xs/src/slic3r/GUI/GUI.hpp
@@ -5,6 +5,8 @@ namespace Slic3r { namespace GUI {
 
 void disable_screensaver();
 void enable_screensaver();
+bool debugged();
+void break_to_debugger();
 
 } }
 
diff --git a/xs/xsp/ExtrusionSimulator.xsp b/xs/xsp/ExtrusionSimulator.xsp
new file mode 100644
index 0000000000..9395913b41
--- /dev/null
+++ b/xs/xsp/ExtrusionSimulator.xsp
@@ -0,0 +1,50 @@
+%module{Slic3r::XS};
+
+%{
+#include <xsinit.h>
+#include "libslic3r/ExtrusionSimulator.hpp"
+%}
+
+%name{Slic3r::ExtrusionSimulator} class ExtrusionSimulator {
+    ~ExtrusionSimulator();
+    %name{_new} ExtrusionSimulator();
+
+    Clone<ExtrusionSimulator> clone()
+        %code{% RETVAL = THIS; %};
+    
+    void  set_image_size(Point *image_size)
+        %code{% THIS->set_image_size(*image_size); %};
+    void  set_viewport(BoundingBox *viewport)
+        %code{% THIS->set_viewport(*viewport); %};
+    void  set_bounding_box(BoundingBox *bbox)
+        %code{% THIS->set_bounding_box(*bbox); %};
+
+    void  reset_accumulator();
+    void  extrude_to_accumulator(ExtrusionPath *path, Point *shift, ExtrusionSimulationType simulationType)
+        %code{% THIS->extrude_to_accumulator(*path, *shift, simulationType); %};
+    void  evaluate_accumulator(ExtrusionSimulationType simulationType);
+    void* image_ptr()
+        %code{% RETVAL = const_cast<void*>(const_cast<Slic3r::ExtrusionSimulator*>(THIS)->image_ptr()); %};
+
+%{
+
+%}
+};
+
+%package{Slic3r::ExtrusionSimulator};
+%{
+
+IV
+_constant()
+  ALIAS:
+    EXTRSIM_SIMPLE          = ExtrusionSimulationSimple
+    EXTRSIM_DONT_SPREAD     = ExtrusionSimulationDontSpread
+    EXTRSIM_SPREAD_NFULL     = ExtrisopmSimulationSpreadNotOverfilled
+    EXTRSIM_SPREAD_FULL     = ExtrusionSimulationSpreadFull
+    EXTRSIM_SPREAD_EXCESS   = ExtrusionSimulationSpreadExcess
+  PROTOTYPE:
+  CODE:
+    RETVAL = ix;
+  OUTPUT: RETVAL
+
+%}
diff --git a/xs/xsp/Filler.xsp b/xs/xsp/Filler.xsp
new file mode 100644
index 0000000000..b7871664f3
--- /dev/null
+++ b/xs/xsp/Filler.xsp
@@ -0,0 +1,68 @@
+%module{Slic3r::XS};
+
+%{
+#include <xsinit.h>
+#include "libslic3r/Fill/FillBase.hpp"
+%}
+
+%name{Slic3r::Filler} class Filler {
+    ~Filler();
+
+    void set_bounding_box(BoundingBox *bbox)
+        %code{% THIS->fill->set_bounding_box(*bbox); %};
+    void set_spacing(coordf_t spacing)
+        %code{% THIS->fill->spacing = spacing; %};
+    coordf_t spacing()
+        %code{% RETVAL = THIS->fill->spacing; %};
+    void set_layer_id(size_t layer_id)
+        %code{% THIS->fill->layer_id = layer_id; %};
+    void set_z(coordf_t z)
+        %code{% THIS->fill->z = z; %};
+    void set_angle(float angle)
+        %code{% THIS->fill->angle = angle; %};
+    void set_loop_clipping(coordf_t clipping)
+        %code{% THIS->fill->loop_clipping = clipping; %};
+
+    bool use_bridge_flow()
+        %code{% RETVAL = THIS->fill->use_bridge_flow(); %};
+    bool no_sort()
+        %code{% RETVAL = THIS->fill->no_sort(); %};
+
+    void set_width(float width)
+        %code{% THIS->params.width = width; %};
+    void set_density(float density)
+        %code{% THIS->params.density = density; %};
+    void set_distance(float distance)
+        %code{% THIS->params.distance = distance; %};
+    void set_dont_connect(bool dont_connect)
+        %code{% THIS->params.dont_connect = dont_connect; %};
+    void set_dont_adjust(bool dont_adjust)
+        %code{% THIS->params.dont_adjust = dont_adjust; %};
+    void set_complete(bool complete)
+        %code{% THIS->params.complete = complete; %};
+
+    PolylineCollection* _fill_surface(Surface *surface)
+        %code{% 
+            PolylineCollection *pc = NULL;
+            if (THIS->fill != NULL) {
+                pc = new PolylineCollection();
+                pc->polylines = THIS->fill->fill_surface(surface, THIS->params);
+            }
+            RETVAL =  pc;
+        %};
+
+%{
+
+Filler*
+new_from_type(CLASS, type)
+    char*               CLASS;
+    std::string         type;
+    CODE:
+        Filler *filler = new Filler();
+        filler->fill = Fill::new_from_type(type);
+        RETVAL = filler;
+    OUTPUT:
+        RETVAL
+
+%}
+};
diff --git a/xs/xsp/GUI.xsp b/xs/xsp/GUI.xsp
index a422040e16..875c49e269 100644
--- a/xs/xsp/GUI.xsp
+++ b/xs/xsp/GUI.xsp
@@ -13,3 +13,9 @@ void disable_screensaver()
 
 void enable_screensaver()
     %code{% Slic3r::GUI::enable_screensaver(); %};
+
+bool debugged()
+    %code{% RETVAL=Slic3r::GUI::debugged(); %};
+
+void break_to_debugger()
+    %code{% Slic3r::GUI::break_to_debugger(); %};
diff --git a/xs/xsp/my.map b/xs/xsp/my.map
index 5f87f51c2c..4d4d5d997b 100644
--- a/xs/xsp/my.map
+++ b/xs/xsp/my.map
@@ -112,6 +112,14 @@ ExtrusionLoop*             O_OBJECT_SLIC3R
 Ref<ExtrusionLoop>         O_OBJECT_SLIC3R_T
 Clone<ExtrusionLoop>       O_OBJECT_SLIC3R_T
 
+ExtrusionSimulator*        O_OBJECT_SLIC3R
+Ref<ExtrusionSimulator>    O_OBJECT_SLIC3R_T
+Clone<ExtrusionSimulator>  O_OBJECT_SLIC3R_T
+
+Filler*                    O_OBJECT_SLIC3R
+Ref<Filler>                O_OBJECT_SLIC3R_T
+Clone<Filler>              O_OBJECT_SLIC3R_T
+
 Flow*                      O_OBJECT_SLIC3R
 Ref<Flow>                  O_OBJECT_SLIC3R_T
 Clone<Flow>                O_OBJECT_SLIC3R_T
@@ -214,6 +222,7 @@ GLVertexArray*             O_OBJECT_SLIC3R
 Axis                  T_UV
 ExtrusionLoopRole     T_UV
 ExtrusionRole     T_UV
+ExtrusionSimulationType T_UV
 FlowRole     T_UV
 PrintStep     T_UV
 PrintObjectStep     T_UV
diff --git a/xs/xsp/typemap.xspt b/xs/xsp/typemap.xspt
index 37d7a9620c..3453f0067a 100644
--- a/xs/xsp/typemap.xspt
+++ b/xs/xsp/typemap.xspt
@@ -57,6 +57,9 @@
 %typemap{ExPolygonCollection*};
 %typemap{Ref<ExPolygonCollection>}{simple};
 %typemap{Clone<ExPolygonCollection>}{simple};
+%typemap{Filler*};
+%typemap{Ref<Filler>}{simple};
+%typemap{Clone<Filler>}{simple};
 %typemap{Flow*};
 %typemap{Ref<Flow>}{simple};
 %typemap{Clone<Flow>}{simple};
@@ -81,6 +84,9 @@
 %typemap{ExtrusionLoop*};
 %typemap{Ref<ExtrusionLoop>}{simple};
 %typemap{Clone<ExtrusionLoop>}{simple};
+%typemap{ExtrusionSimulator*};
+%typemap{Ref<ExtrusionSimulator>}{simple};
+%typemap{Clone<ExtrusionSimulator>}{simple};
 %typemap{TriangleMesh*};
 %typemap{Ref<TriangleMesh>}{simple};
 %typemap{Clone<TriangleMesh>}{simple};
@@ -223,6 +229,12 @@
     $CVar = (ExtrusionRole)SvUV($PerlVar);
   %};
 };
+%typemap{ExtrusionSimulationType}{parsed}{
+  %cpp_type{ExtrusionSimulationType};
+  %precall_code{%
+    $CVar = (ExtrusionSimulationType)SvUV($PerlVar);
+  %};
+};
 %typemap{FlowRole}{parsed}{
   %cpp_type{FlowRole};
   %precall_code{%