Browse Source

Initial commit of the new Pressure Equalizer, the EdgeGrid
signed distance field structure.
The EdgeGrid is used to avoid placing the seams on overhangs.

bubnikv 8 years ago
parent
commit
f518e0675c

+ 15 - 3
lib/Slic3r/Print/GCode.pm

@@ -10,6 +10,7 @@ has '_spiral_vase'                   => (is => 'rw');
 has '_vibration_limit'               => (is => 'rw');
 has '_arc_fitting'                   => (is => 'rw');
 has '_pressure_regulator'            => (is => 'rw');
+has '_pressure_equalizer'            => (is => 'rw');
 has '_skirt_done'                    => (is => 'rw', default => sub { {} });  # print_z => 1
 has '_brim_done'                     => (is => 'rw');
 has '_second_layer_things_done'      => (is => 'rw');
@@ -114,6 +115,11 @@ sub BUILD {
     
     $self->_pressure_regulator(Slic3r::GCode::PressureRegulator->new(config => $self->config))
         if $self->config->pressure_advance > 0;
+
+    $self->_pressure_equalizer(Slic3r::GCode::PressureEqualizer->new($self->config))
+        if defined($ENV{"SLIC3R_PRESSURE_EQUALIZER"}) && $ENV{'SLIC3R_PRESSURE_EQUALIZER'} == 1;
+
+    $self->_gcodegen->set_enable_extrusion_role_markers(defined $self->_pressure_equalizer);
 }
 
 # Export a G-code for the complete print.
@@ -560,7 +566,7 @@ sub process_layer {
                     }
                 }
             }
-        }
+        } # for regions
         
         # tweak extruder ordering to save toolchanges
         my @extruders = sort keys %by_extruder;
@@ -586,7 +592,7 @@ sub process_layer {
                 }
             }
         }
-    }
+    } # for object copies
     
     # apply spiral vase post-processing if this layer contains suitable geometry
     # (we must feed all the G-code into the post-processor, including the first 
@@ -658,7 +664,13 @@ sub filter {
     # this depends on actual speeds
     $gcode = $self->_pressure_regulator->process($gcode, $flush)
         if defined $self->_pressure_regulator;
-    
+
+    # apply pressure equalization if enabled;
+#    print "G-code before filter:\n", $gcode;
+    $gcode = $self->_pressure_equalizer->process($gcode, $flush)
+        if defined $self->_pressure_equalizer;
+#    print "G-code after filter:\n", $gcode;
+
     # apply arc fitting if enabled;
     # this does not depend on speeds but changes G1 XY commands into G2/G2 IJ
     $gcode = $self->_arc_fitting->process($gcode)

+ 20 - 6
slic3r.sublime-project

@@ -11,30 +11,44 @@
 			"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\""
+			"shell_cmd": "chdir & perl slic3r.pl  --DataDir \"C:\\Users\\Public\\Documents\\Prusa3D\\Slic3r settings MK2\" --gui \"..\\Slic3r-tests\\gap fill torture 20 -rt.stl\""
 		},
 		{
 			"name": "full",
    			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
-			"shell_cmd": "perl Build.pl"
+			"shell_cmd": "chdir & perl Build.pl"
 		},
 		{
 			"name": "xs",
 			"working_dir": "$project_path/xs",
-  			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
-			"shell_cmd": "perl Build install"
+			// for Visual Studio:
+  			"file_regex": "^(..[^:]*)\\(([0-9]+)\\)(.*)$",
+  			// For GCC:
+// 			"file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$",
+			"shell_cmd": "chdir & perl Build install",
+			"env": {
+//				"PATH": "C:\\Program Files (x86)\\MSBuild\\12.0\\bin\\amd64;C:\\Program Files (x86)\\Microsoft Visual Studio 12.0\\VC\\BIN\\amd64;C:\\Program Files (x86)\\Microsoft Visual Studio 12.0\\Common7\\IDE;C:\\Program Files (x86)\\Microsoft Visual Studio 12.0\\Common7\\Tools;%PATH%;c:\\wperl64d\\site\\bin;c:\\wperl64d\\bin",
+//				"PERL_CPANM_HOME": "c:\\wperl64d\\cpanm",
+//				"WXDIR": "D:\\src-perl\\wxWidgets-3.0.3-beta1",
+//				"BOOST_DIR": "D:\\src-perl\\boost_1_61_0",
+//				"BOOST_INCLUDEDIR": "D:\\src-perl\\boost_1_61_0",
+//				"BOOST_LIBRARYDIR": "D:\\src-perl\\boost_1_61_0\\stage\\x64\\lib",
+//				"SLIC3R_STATIC": "1"
+			}
 		},
 		{
 			"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\""
+			"shell_cmd": "chdir & perl Build install & cd .. & perl slic3r.pl --gui \"..\\Slic3r-tests\\star3-big2.stl\""
 		}
 	],
 	"folders":
 	[
 		{
-			"path": "."
+			"path": ".",
+//	        "folder_exclude_patterns": [".svn", "._d", ".metadata", ".settings"],
+	        "file_exclude_patterns": ["XS.c"]
 		}
 	],
 

+ 15 - 1
xs/Build.PL

@@ -25,8 +25,22 @@ if ($mswin) {
     push @cflags, qw(-DNOMINMAX -D_USE_MATH_DEFINES);
 }
 
+my @early_includes = ();
 my @INC  = qw(-Isrc/libslic3r);
 my @LIBS = $cpp_guess->is_msvc ? qw(-LIBPATH:src/libslic3r) : qw(-Lsrc/libslic3r);
+   
+#if ($ENV{SLIC3R_GUI}) 
+{
+    print "Slic3r will be built with GUI support\n";
+    require Alien::wxWidgets;
+    Alien::wxWidgets->load;
+    push @INC, Alien::wxWidgets->include_path;
+    push @cflags, qw(-DSLIC3R_GUI -DUNICODE), Alien::wxWidgets->defines, Alien::wxWidgets->c_flags;
+    my $alienwx_libraries = Alien::wxWidgets->libraries(qw(gl html));
+    $alienwx_libraries =~ s/-L/-LIBPATH:/g if ($cpp_guess->is_msvc);
+    push @ldflags, Alien::wxWidgets->link_flags, $alienwx_libraries;
+#    push @early_includes, qw(slic3r/GUI/wxinit.h);
+}
 
 # search for Boost in a number of places
 my @boost_include = ();
@@ -195,7 +209,7 @@ my $build = Module::Build::WithXSpp->new(
         ostream
         sstream
         libslic3r/GCodeSender.hpp
-    )]
+    ), @early_includes]
 );
 
 $build->create_build_script;

+ 3 - 0
xs/MANIFEST

@@ -52,6 +52,8 @@ src/libslic3r/GCodeSender.cpp
 src/libslic3r/GCodeSender.hpp
 src/libslic3r/GCodeWriter.cpp
 src/libslic3r/GCodeWriter.hpp
+src/libslic3r/GCode/PressureEqualizer.cpp
+src/libslic3r/GCode/PressureEqualizer.hpp
 src/libslic3r/Geometry.cpp
 src/libslic3r/Geometry.hpp
 src/libslic3r/Layer.cpp
@@ -151,6 +153,7 @@ xsp/Flow.xsp
 xsp/GCode.xsp
 xsp/GCodeSender.xsp
 xsp/GCodeWriter.xsp
+xsp/GCodePressureEqualizer.xsp
 xsp/Geometry.xsp
 xsp/GUI.xsp
 xsp/GUI_3DScene.xsp

+ 13 - 0
xs/lib/Slic3r/XS.pm

@@ -4,6 +4,19 @@ use strict;
 
 our $VERSION = '0.01';
 
+# We have to load these modules in order to have Wx.pm find the correct paths
+# for wxWidgets dlls on MSW.
+# We avoid loading these on OS X because Wx::Load() initializes a Wx App
+# automatically and it steals focus even when we're not running Slic3r in GUI mode.
+# TODO: only load these when compiling with GUI support
+BEGIN {
+    if ($^O eq 'MSWin32') {
+        eval "use Wx";
+#        eval "use Wx::Html";
+        eval "use Wx::Print";  # because of some Wx bug, thread creation fails if we don't have this (looks like Wx::Printout is hard-coded in some thread cleanup code)
+    }
+}
+
 use Carp qw();
 use XSLoader;
 XSLoader::load(__PACKAGE__, $VERSION);

+ 1025 - 0
xs/src/libslic3r/EdgeGrid.cpp

@@ -0,0 +1,1025 @@
+#include <algorithm>
+#include <vector>
+
+#include <wx/image.h>
+
+#include "libslic3r.h"
+#include "EdgeGrid.hpp"
+
+#if 0
+// Enable debugging and assert in this file.
+#define DEBUG
+#define _DEBUG
+#undef NDEBUG
+#endif
+
+#include <assert.h>
+
+namespace Slic3r {
+
+EdgeGrid::Grid::Grid() : 
+	m_rows(0), m_cols(0) 
+{
+}
+
+EdgeGrid::Grid::~Grid() 
+{
+	m_contours.clear();
+	m_cell_data.clear();
+	m_cells.clear();
+}
+
+void EdgeGrid::Grid::create(const Polygons &polygons, coord_t resolution)
+{
+	// Count the contours.
+	size_t ncontours = 0;
+	for (size_t j = 0; j < polygons.size(); ++ j)
+		if (! polygons[j].points.empty())
+			++ ncontours;
+
+	// Collect the contours.
+	m_contours.assign(ncontours, NULL);
+	ncontours = 0;
+	for (size_t j = 0; j < polygons.size(); ++ j)
+		if (! polygons[j].points.empty())
+			m_contours[ncontours++] = &polygons[j].points;
+
+	create_from_m_contours(resolution);
+}
+
+void EdgeGrid::Grid::create(const ExPolygon &expoly, coord_t resolution)
+{
+	// Count the contours.
+	size_t ncontours = 0;
+	if (! expoly.contour.points.empty())
+		++ ncontours;
+	for (size_t j = 0; j < expoly.holes.size(); ++ j)
+		if (! expoly.holes[j].points.empty())
+			++ ncontours;
+
+	// Collect the contours.
+	m_contours.assign(ncontours, NULL);
+	ncontours = 0;
+	if (! expoly.contour.points.empty())
+		m_contours[ncontours++] = &expoly.contour.points;
+	for (size_t j = 0; j < expoly.holes.size(); ++ j)
+		if (! expoly.holes[j].points.empty())
+			m_contours[ncontours++] = &expoly.holes[j].points;
+
+	create_from_m_contours(resolution);
+}
+
+void EdgeGrid::Grid::create(const ExPolygons &expolygons, coord_t resolution)
+{
+	// Count the contours.
+	size_t ncontours = 0;
+	for (size_t i = 0; i < expolygons.size(); ++ i) {
+		const ExPolygon &expoly = expolygons[i];
+		if (! expoly.contour.points.empty())
+			++ ncontours;
+		for (size_t j = 0; j < expoly.holes.size(); ++ j)
+			if (! expoly.holes[j].points.empty())
+				++ ncontours;
+	}
+
+	// Collect the contours.
+	m_contours.assign(ncontours, NULL);
+	ncontours = 0;
+	for (size_t i = 0; i < expolygons.size(); ++ i) {
+		const ExPolygon &expoly = expolygons[i];
+		if (! expoly.contour.points.empty())
+			m_contours[ncontours++] = &expoly.contour.points;
+		for (size_t j = 0; j < expoly.holes.size(); ++ j)
+			if (! expoly.holes[j].points.empty())
+				m_contours[ncontours++] = &expoly.holes[j].points;
+	}
+
+	create_from_m_contours(resolution);
+}
+
+void EdgeGrid::Grid::create(const ExPolygonCollection &expolygons, coord_t resolution)
+{
+	create(expolygons.expolygons, resolution);
+}
+
+// m_contours has been initialized. Now fill in the edge grid.
+void EdgeGrid::Grid::create_from_m_contours(coord_t resolution)
+{
+	// 1) Measure the bounding box.
+	m_bbox.defined = false;
+	for (size_t i = 0; i < m_contours.size(); ++ i) {
+		const Slic3r::Points &pts = *m_contours[i];
+		for (size_t j = 0; j < pts.size(); ++ j)
+			m_bbox.merge(pts[j]);
+	}
+	coord_t eps = 16;
+	m_bbox.min.x -= eps;
+	m_bbox.min.y -= eps;
+	m_bbox.max.x += eps;
+	m_bbox.max.y += eps;
+
+	// 2) Initialize the edge grid.
+	m_resolution = resolution;
+	m_cols = (m_bbox.max.x - m_bbox.min.x + m_resolution - 1) / m_resolution;
+	m_rows = (m_bbox.max.y - m_bbox.min.y + m_resolution - 1) / m_resolution;
+	m_cells.assign(m_rows * m_cols, Cell());
+
+	// 3) First round of contour rasterization, count the edges per grid cell.
+	for (size_t i = 0; i < m_contours.size(); ++ i) {
+		const Slic3r::Points &pts = *m_contours[i];
+		for (size_t j = 0; j < pts.size(); ++ j) {
+			// End points of the line segment.
+			Slic3r::Point p1(pts[j]);
+			Slic3r::Point p2 = pts[(j + 1 == pts.size()) ? 0 : j + 1];
+			p1.x -= m_bbox.min.x;
+			p1.y -= m_bbox.min.y;
+			p2.x -= m_bbox.min.x;
+			p2.y -= m_bbox.min.y;
+		   	// Get the cells of the end points.
+		    coord_t ix    = p1.x / m_resolution;
+		    coord_t iy    = p1.y / m_resolution;
+		    coord_t ixb   = p2.x / m_resolution;
+		    coord_t iyb   = p2.y / m_resolution;
+			assert(ix >= 0 && ix < m_cols);
+			assert(iy >= 0 && iy < m_rows);
+			assert(ixb >= 0 && ixb < m_cols);
+			assert(iyb >= 0 && iyb < m_rows);
+			// Account for the end points.
+			++ m_cells[iy*m_cols+ix].end;
+			if (ix == ixb && iy == iyb)
+				// Both ends fall into the same cell.
+				continue;
+		    // Raster the centeral part of the line.
+			coord_t dx = std::abs(p2.x - p1.x);
+			coord_t dy = std::abs(p2.y - p1.y);
+			if (p1.x < p2.x) {
+				int64_t ex = int64_t((ix + 1)*m_resolution - p1.x) * int64_t(dy);
+				if (p1.y < p2.y) {
+					int64_t ey = int64_t((iy + 1)*m_resolution - p1.y) * int64_t(dx);
+					do {
+						assert(ix <= ixb && iy <= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix += 1;
+						}
+						else if (ex == ey) {
+							ex = int64_t(dy) * m_resolution;
+							ey = int64_t(dx) * m_resolution;
+							ix += 1;
+							iy += 1;
+						}
+						else {
+							assert(ex > ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy += 1;
+						}
+						++m_cells[iy*m_cols + ix].end;
+					} while (ix != ixb || iy != iyb);
+				}
+				else {
+					int64_t ey = int64_t(p1.y - iy*m_resolution) * int64_t(dx);
+					do {
+						assert(ix <= ixb && iy >= iyb);
+						if (ex <= ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix += 1;
+						}
+						else {
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy -= 1;
+						}
+						++m_cells[iy*m_cols + ix].end;
+					} while (ix != ixb || iy != iyb);
+				}
+			}
+			else {
+				int64_t ex = int64_t(p1.x - ix*m_resolution) * int64_t(dy);
+				if (p1.y < p2.y) {
+					int64_t ey = int64_t((iy + 1)*m_resolution - p1.y) * int64_t(dx);
+					do {
+						assert(ix >= ixb && iy <= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix -= 1;
+						}
+						else {
+							assert(ex >= ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy += 1;
+						}
+						++m_cells[iy*m_cols + ix].end;
+					} while (ix != ixb || iy != iyb);
+				}
+				else {
+					int64_t ey = int64_t(p1.y - iy*m_resolution) * int64_t(dx);
+					do {
+						assert(ix >= ixb && iy >= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix -= 1;
+						}
+						else if (ex == ey) {
+							ex = int64_t(dy) * m_resolution;
+							ey = int64_t(dx) * m_resolution;
+							ix -= 1;
+							iy -= 1;
+						}
+						else {
+							assert(ex > ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy -= 1;
+						}
+						++m_cells[iy*m_cols + ix].end;
+					} while (ix != ixb || iy != iyb);
+				}
+			}
+		}
+	}
+
+	// 4) Prefix sum the numbers of hits per cells to get an index into m_cell_data.
+	size_t cnt = m_cells.front().end;
+	for (size_t i = 1; i < m_cells.size(); ++ i) {
+		m_cells[i].begin = cnt;
+		cnt += m_cells[i].end;
+		m_cells[i].end = cnt;
+	}
+
+	// 5) Allocate the cell data.
+	m_cell_data.assign(cnt, std::pair<size_t, size_t>(size_t(-1), size_t(-1)));
+
+	// 6) Finally fill in m_cell_data by rasterizing the lines once again.
+	for (size_t i = 0; i < m_cells.size(); ++i)
+		m_cells[i].end = m_cells[i].begin;
+	for (size_t i = 0; i < m_contours.size(); ++i) {
+		const Slic3r::Points &pts = *m_contours[i];
+		for (size_t j = 0; j < pts.size(); ++j) {
+			// End points of the line segment.
+			Slic3r::Point p1(pts[j]);
+			Slic3r::Point p2 = pts[(j + 1 == pts.size()) ? 0 : j + 1];
+			p1.x -= m_bbox.min.x;
+			p1.y -= m_bbox.min.y;
+			p2.x -= m_bbox.min.x;
+			p2.y -= m_bbox.min.y;
+			// Get the cells of the end points.
+			coord_t ix = p1.x / m_resolution;
+			coord_t iy = p1.y / m_resolution;
+			coord_t ixb = p2.x / m_resolution;
+			coord_t iyb = p2.y / m_resolution;
+			assert(ix >= 0 && ix < m_cols);
+			assert(iy >= 0 && iy < m_rows);
+			assert(ixb >= 0 && ixb < m_cols);
+			assert(iyb >= 0 && iyb < m_rows);
+			// Account for the end points.
+			m_cell_data[m_cells[iy*m_cols + ix].end++] = std::pair<size_t, size_t>(i, j);
+			if (ix == ixb && iy == iyb)
+				// Both ends fall into the same cell.
+				continue;
+			// Raster the centeral part of the line.
+			coord_t dx = std::abs(p2.x - p1.x);
+			coord_t dy = std::abs(p2.y - p1.y);
+			if (p1.x < p2.x) {
+				int64_t ex = int64_t((ix + 1)*m_resolution - p1.x) * int64_t(dy);
+				if (p1.y < p2.y) {
+					int64_t ey = int64_t((iy + 1)*m_resolution - p1.y) * int64_t(dx);
+					do {
+						assert(ix <= ixb && iy <= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix += 1;
+						}
+						else if (ex == ey) {
+							ex = int64_t(dy) * m_resolution;
+							ey = int64_t(dx) * m_resolution;
+							ix += 1;
+							iy += 1;
+						}
+						else {
+							assert(ex > ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy += 1;
+						}
+						m_cell_data[m_cells[iy*m_cols + ix].end++] = std::pair<size_t, size_t>(i, j);
+					} while (ix != ixb || iy != iyb);
+				}
+				else {
+					int64_t ey = int64_t(p1.y - iy*m_resolution) * int64_t(dx);
+					do {
+						assert(ix <= ixb && iy >= iyb);
+						if (ex <= ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix += 1;
+						}
+						else {
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy -= 1;
+						}
+						m_cell_data[m_cells[iy*m_cols + ix].end++] = std::pair<size_t, size_t>(i, j);
+					} while (ix != ixb || iy != iyb);
+				}
+			}
+			else {
+				int64_t ex = int64_t(p1.x - ix*m_resolution) * int64_t(dy);
+				if (p1.y < p2.y) {
+					int64_t ey = int64_t((iy + 1)*m_resolution - p1.y) * int64_t(dx);
+					do {
+						assert(ix >= ixb && iy <= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix -= 1;
+						}
+						else {
+							assert(ex >= ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy += 1;
+						}
+						m_cell_data[m_cells[iy*m_cols + ix].end++] = std::pair<size_t, size_t>(i, j);
+					} while (ix != ixb || iy != iyb);
+				}
+				else {
+					int64_t ey = int64_t(p1.y - iy*m_resolution) * int64_t(dx);
+					do {
+						assert(ix >= ixb && iy >= iyb);
+						if (ex < ey) {
+							ey -= ex;
+							ex = int64_t(dy) * m_resolution;
+							ix -= 1;
+						}
+						else if (ex == ey) {
+							ex = int64_t(dy) * m_resolution;
+							ey = int64_t(dx) * m_resolution;
+							ix -= 1;
+							iy -= 1;
+						}
+						else {
+							assert(ex > ey);
+							ex -= ey;
+							ey = int64_t(dx) * m_resolution;
+							iy -= 1;
+						}
+						m_cell_data[m_cells[iy*m_cols + ix].end++] = std::pair<size_t, size_t>(i, j);
+					} while (ix != ixb || iy != iyb);
+				}
+			}
+		}
+	}
+}
+
+template<const int INCX, const int INCY>
+struct PropagateDanielssonSingleStep {
+	PropagateDanielssonSingleStep(float *aL, unsigned char *asigns, size_t astride, coord_t aresolution) :
+		L(aL), signs(asigns), stride(astride), resolution(aresolution) {}
+	inline void operator()(int r, int c, int addr_delta) {
+		size_t  addr = r * stride + c;
+		if ((signs[addr] & 2) == 0) {
+			float  *v = &L[addr << 1];
+			float   l = v[0] * v[0] + v[1] * v[1];
+			float  *v2s = v + (addr_delta << 1);
+			float	v2[2] = {
+				v2s[0] + INCX * resolution,
+				v2s[1] + INCY * resolution
+			};
+			float l2 = v2[0] * v2[0] + v2[1] * v2[1];
+			if (l2 < l) {
+				v[0] = v2[0];
+				v[1] = v2[1];
+			}
+		}
+	}
+	float		   *L;
+	unsigned char  *signs;
+	size_t			stride;
+	coord_t			resolution;
+};
+
+struct PropagateDanielssonSingleVStep3 {
+	PropagateDanielssonSingleVStep3(float *aL, unsigned char *asigns, size_t astride, coord_t aresolution) :
+		L(aL), signs(asigns), stride(astride), resolution(aresolution) {}
+	inline void operator()(int r, int c, int addr_delta, bool has_l, bool has_r) {
+		size_t  addr = r * stride + c;
+		if ((signs[addr] & 2) == 0) {
+			float  *v    = &L[addr<<1];
+			float   l    = v[0]*v[0]+v[1]*v[1];
+			float  *v2s   = v+(addr_delta<<1);
+			float	v2[2] = {
+				v2s[0],
+				v2s[1] + resolution
+			};
+			float   l2   = v2[0]*v2[0]+v2[1]*v2[1];
+			if (l2 < l) {
+				v[0] = v2[0];
+				v[1] = v2[1];
+			}
+			if (has_l) {
+				float  *v2sl = v2s - 1;
+				v2[0] = v2sl[0] + resolution;
+				v2[1] = v2sl[1] + resolution;
+				l2   = v2[0]*v2[0]+v2[1]*v2[1];
+				if (l2 < l) {
+					v[0] = v2[0];
+					v[1] = v2[1];
+				}
+			}
+			if (has_r) {
+				float  *v2sr = v2s + 1;
+				v2[0] = v2sr[0] + resolution;
+				v2[1] = v2sr[1] + resolution;
+				l2   = v2[0]*v2[0]+v2[1]*v2[1];
+				if (l2 < l) {
+					v[0] = v2[0];
+					v[1] = v2[1];
+				}
+			}
+		}
+	}
+	float		   *L;
+	unsigned char  *signs;
+	size_t			stride;
+	coord_t			resolution;
+};
+
+void EdgeGrid::Grid::calculate_sdf()
+{
+	// 1) Initialize a signum and an unsigned vector to a zero iso surface.
+	size_t nrows = m_rows + 1;
+	size_t ncols = m_cols + 1;
+	// Unsigned vectors towards the closest point on the surface.
+	std::vector<float> L(nrows * ncols * 2, FLT_MAX);
+	// Bit 0 set - negative.
+	// Bit 1 set - original value, the distance value shall not be changed by the Danielsson propagation.
+	// Bit 2 set - signum not propagated yet.
+	std::vector<unsigned char> signs(nrows * ncols, 4);
+	// SDF will be initially filled with unsigned DF.
+//	m_signed_distance_field.assign(nrows * ncols, FLT_MAX);
+	float search_radius = float(m_resolution<<1);
+	m_signed_distance_field.assign(nrows * ncols, search_radius);
+	// For each cell:
+	for (size_t r = 0; r < m_rows; ++ r) {
+		for (size_t c = 0; c < m_cols; ++ c) {
+			const Cell &cell = m_cells[r * m_cols + c];
+			// For each segment in the cell:
+			for (size_t i = cell.begin; i != cell.end; ++ i) {
+				const Slic3r::Points &pts = *m_contours[m_cell_data[i].first];
+				size_t ipt = m_cell_data[i].second;
+				// End points of the line segment.
+				const Slic3r::Point &p1 = pts[ipt];
+				const Slic3r::Point &p2 = pts[(ipt + 1 == pts.size()) ? 0 : ipt + 1];
+				// Segment vector
+				const Slic3r::Point v_seg = p1.vector_to(p2);
+				// l2 of v_seg
+				const int64_t l2_seg = int64_t(v_seg.x) * int64_t(v_seg.x) + int64_t(v_seg.y) * int64_t(v_seg.y);
+				// For each corner of this cell and its 1 ring neighbours:
+				for (int corner_y = -1; corner_y < 3; ++ corner_y) {
+					coord_t corner_r = r + corner_y;
+					if (corner_r < 0 || corner_r >= nrows)
+						continue;
+					for (int corner_x = -1; corner_x < 3; ++ corner_x) {
+						coord_t corner_c = c + corner_x;
+						if (corner_c < 0 || corner_c >= ncols)
+							continue;
+						float  &d_min = m_signed_distance_field[corner_r * ncols + corner_c];
+						Slic3r::Point pt(m_bbox.min.x + corner_c * m_resolution, m_bbox.min.y + corner_r * m_resolution);
+						Slic3r::Point v_pt = p1.vector_to(pt);
+						// dot(p2-p1, pt-p1)
+						int64_t t_pt = int64_t(v_seg.x) * int64_t(v_pt.x) + int64_t(v_seg.y) * int64_t(v_pt.y);
+						if (t_pt < 0) {
+							// Closest to p1.
+							double dabs = sqrt(int64_t(v_pt.x) * int64_t(v_pt.x) + int64_t(v_pt.y) * int64_t(v_pt.y));
+							if (dabs < d_min) {
+								// Previous point.
+								const Slic3r::Point &p0 = pts[(ipt == 0) ? (pts.size() - 1) : ipt - 1];
+								Slic3r::Point v_seg_prev = p0.vector_to(p1);
+								int64_t t2_pt = int64_t(v_seg_prev.x) * int64_t(v_pt.x) + int64_t(v_seg_prev.y) * int64_t(v_pt.y);
+								if (t2_pt > 0) {
+									// Inside the wedge between the previous and the next segment.
+									// Set the signum depending on whether the vertex is convex or reflex.
+									int64_t det = int64_t(v_seg_prev.x) * int64_t(v_seg.y) - int64_t(v_seg_prev.y) * int64_t(v_seg.x);
+									assert(det != 0);
+									d_min = dabs;
+									// Fill in an unsigned vector towards the zero iso surface.
+									float *l = &L[(corner_r * ncols + corner_c) << 1];
+									l[0] = std::abs(v_pt.x);
+									l[1] = std::abs(v_pt.y);
+								#ifdef _DEBUG
+									double dabs2 = sqrt(l[0]*l[0]+l[1]*l[1]);
+									assert(std::abs(dabs-dabs2) < 1e-4 * std::max(dabs, dabs2));
+								#endif /* _DEBUG */
+									signs[corner_r * ncols + corner_c] = ((det < 0) ? 1 : 0) | 2;
+								}
+							}
+						}
+						else if (t_pt > l2_seg) {
+							// Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the same cell.
+							continue;
+						} else {
+							// Closest to the segment.
+							assert(t_pt >= 0 && t_pt <= l2_seg);
+							int64_t d_seg = int64_t(v_seg.y) * int64_t(v_pt.x) - int64_t(v_seg.x) * int64_t(v_pt.y);
+							double d = double(d_seg) / sqrt(double(l2_seg));
+							double dabs = std::abs(d);
+							if (dabs < d_min) {
+								d_min = dabs;
+								// Fill in an unsigned vector towards the zero iso surface.
+								float *l = &L[(corner_r * ncols + corner_c) << 1];
+								float linv = float(d_seg) / float(l2_seg);
+								l[0] = std::abs(float(v_seg.y) * linv);
+								l[1] = std::abs(float(v_seg.x) * linv);
+								#ifdef _DEBUG
+									double dabs2 = sqrt(l[0]*l[0]+l[1]*l[1]);
+									assert(std::abs(dabs-dabs2) <= 1e-4 * std::max(dabs, dabs2));
+								#endif /* _DEBUG */
+								signs[corner_r * ncols + corner_c] = ((d_seg < 0) ? 1 : 0) | 2;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+#if 1
+	{ 
+		wxImage img(ncols, nrows);
+		unsigned char *data = img.GetData();
+		memset(data, 0, ncols * nrows * 3);
+		for (coord_t r = 0; r < nrows; ++r) {
+			for (coord_t c = 0; c < ncols; ++c) {
+				unsigned char *pxl = data + (((nrows - r - 1) * ncols) + c) * 3;
+				float d = m_signed_distance_field[r * ncols + c];
+				if (d != search_radius) {
+					float s = 255 * d / search_radius;
+					int is = std::max(0, std::min(255, int(floor(s + 0.5f))));
+					pxl[0] = 255;
+					pxl[1] = 255 - is;
+					pxl[2] = 255 - is;
+				}
+				else {
+					pxl[0] = 0;
+					pxl[1] = 255;
+					pxl[2] = 0;
+				}
+			}
+		}
+		img.SaveFile("out\\unsigned_df.png", wxBITMAP_TYPE_PNG);
+	}
+	{
+		wxImage img(ncols, nrows);
+		unsigned char *data = img.GetData();
+		memset(data, 0, ncols * nrows * 3);
+		for (coord_t r = 0; r < nrows; ++r) {
+			for (coord_t c = 0; c < ncols; ++c) {
+				unsigned char *pxl = data + (((nrows - r - 1) * ncols) + c) * 3;
+				float d = m_signed_distance_field[r * ncols + c];
+				if (d != search_radius) {
+					float s = 255 * d / search_radius;
+					int is = std::max(0, std::min(255, int(floor(s + 0.5f))));
+					if ((signs[r * ncols + c] & 1) == 0) {
+						// Positive
+						pxl[0] = 255;
+						pxl[1] = 255 - is;
+						pxl[2] = 255 - is;
+					}
+					else {
+						// Negative
+						pxl[0] = 255 - is;
+						pxl[1] = 255 - is;
+						pxl[2] = 255;
+					}
+				}
+				else {
+					pxl[0] = 0;
+					pxl[1] = 255;
+					pxl[2] = 0;
+				}
+			}
+		}
+		img.SaveFile("out\\signed_df.png", wxBITMAP_TYPE_PNG);
+	}
+#endif
+
+	// 2) Propagate the signum.
+	#define PROPAGATE_SIGNUM_SINGLE_STEP(DELTA) do { \
+		size_t 		   addr    = r * ncols + c;				\
+		unsigned char &cur_val = signs[addr]; 				\
+		if (cur_val & 4) {	 								\
+			unsigned char old_val = signs[addr + (DELTA)];  \
+			if ((old_val & 4) == 0)							\
+				cur_val = old_val & 1; 						\
+		} 													\
+	} while (0);
+	// Top to bottom propagation.
+	for (size_t r = 0; r < nrows; ++ r) {
+		if (r > 0)
+			for (size_t c = 0; c < ncols; ++ c)
+				PROPAGATE_SIGNUM_SINGLE_STEP(- int(ncols));
+		for (size_t c = 1; c < ncols; ++ c)
+			PROPAGATE_SIGNUM_SINGLE_STEP(- 1);
+		for (int c = int(ncols) - 2; c >= 0; -- c)
+			PROPAGATE_SIGNUM_SINGLE_STEP(+ 1);
+	}
+	// Bottom to top propagation.
+	for (int r = int(nrows) - 2; r >= 0; -- r) {
+		for (size_t c = 0; c < ncols; ++ c)
+			PROPAGATE_SIGNUM_SINGLE_STEP(+ ncols);
+		for (size_t c = 1; c < ncols; ++ c)
+			PROPAGATE_SIGNUM_SINGLE_STEP(- 1);
+		for (int c = int(ncols) - 2; c >= 0; -- c)
+			PROPAGATE_SIGNUM_SINGLE_STEP(+ 1);
+	}
+	#undef PROPAGATE_SIGNUM_SINGLE_STEP
+
+	// 3) Propagate the distance by the Danielsson chamfer metric.
+	// Top to bottom propagation.
+	PropagateDanielssonSingleStep<1, 0> danielsson_hstep(L.data(), signs.data(), ncols, m_resolution);
+	PropagateDanielssonSingleStep<0, 1> danielsson_vstep(L.data(), signs.data(), ncols, m_resolution);
+	PropagateDanielssonSingleVStep3 	danielsson_vstep3(L.data(), signs.data(), ncols, m_resolution);
+	// Top to bottom propagation.
+	for (size_t r = 0; r < nrows; ++ r) {
+		if (r > 0)
+			for (size_t c = 0; c < ncols; ++ c)
+				danielsson_vstep(r, c, -int(ncols));
+//				PROPAGATE_DANIELSSON_SINGLE_VSTEP3(-int(ncols), c != 0, c + 1 != ncols);
+		for (size_t c = 1; c < ncols; ++ c)
+			danielsson_hstep(r, c, -1);
+		for (int c = int(ncols) - 2; c >= 0; -- c)
+			danielsson_hstep(r, c, +1);
+	}
+	// Bottom to top propagation.
+	for (int r = int(nrows) - 2; r >= 0; -- r) {
+		for (size_t c = 0; c < ncols; ++ c)
+			danielsson_vstep(r, c, +ncols);
+//			PROPAGATE_DANIELSSON_SINGLE_VSTEP3(+int(ncols), c != 0, c + 1 != ncols);
+		for (size_t c = 1; c < ncols; ++ c)
+			danielsson_hstep(r, c, -1);
+		for (int c = int(ncols) - 2; c >= 0; -- c)
+			danielsson_hstep(r, c, +1);
+	}
+
+	// Update signed distance field from absolte vectors to the iso-surface.
+	for (size_t r = 0; r < nrows; ++ r) {
+		for (size_t c = 0; c < ncols; ++ c) {
+			size_t  addr = r * ncols + c;
+			float  *v    = &L[addr<<1];
+			float   d    = sqrt(v[0]*v[0]+v[1]*v[1]);
+			if (signs[addr] & 1)
+				d = -d;
+			m_signed_distance_field[addr] = d;
+		}
+	}
+
+#if 1
+	{
+		wxImage img(ncols, nrows);
+		unsigned char *data = img.GetData();
+		memset(data, 0, ncols * nrows * 3);
+		float search_radius = float(m_resolution * 5);
+		for (coord_t r = 0; r < nrows; ++r) {
+			for (coord_t c = 0; c < ncols; ++c) {
+				unsigned char *pxl = data + (((nrows - r - 1) * ncols) + c) * 3;
+				unsigned char sign = signs[r * ncols + c];
+				switch (sign) {
+				case 0:
+					// Positive, outside of a narrow band.
+					pxl[0] = 0;
+					pxl[1] = 0;
+					pxl[2] = 255;
+					break;
+				case 1:
+					// Negative, outside of a narrow band.
+					pxl[0] = 255;
+					pxl[1] = 0;
+					pxl[2] = 0;
+					break;
+				case 2:
+					// Positive, outside of a narrow band.
+					pxl[0] = 100;
+					pxl[1] = 100;
+					pxl[2] = 255;
+					break;
+				case 3:
+					// Negative, outside of a narrow band.
+					pxl[0] = 255;
+					pxl[1] = 100; 
+					pxl[2] = 100;
+					break;
+				case 4:
+					// This shall not happen. Undefined signum.
+					pxl[0] = 0;
+					pxl[1] = 255;
+					pxl[2] = 0;
+					break;
+				default:
+					// This shall not happen. Invalid signum value.
+					pxl[0] = 255;
+					pxl[1] = 255;
+					pxl[2] = 255;
+					break;
+				}
+			}
+		}
+		img.SaveFile("out\\signed_df-signs.png", wxBITMAP_TYPE_PNG);
+	}
+#endif
+
+#if 1
+	{
+		wxImage img(ncols, nrows);
+		unsigned char *data = img.GetData();
+		memset(data, 0, ncols * nrows * 3);
+		float search_radius = float(m_resolution * 5);
+		for (coord_t r = 0; r < nrows; ++r) {
+			for (coord_t c = 0; c < ncols; ++c) {
+				unsigned char *pxl = data + (((nrows - r - 1) * ncols) + c) * 3;
+				float d = m_signed_distance_field[r * ncols + c];
+				float s = 255.f * fabs(d) / search_radius;
+				int is = std::max(0, std::min(255, int(floor(s + 0.5f))));
+				if (d < 0.f) {
+					pxl[0] = 255;
+					pxl[1] = 255 - is;
+					pxl[2] = 255 - is;
+				}
+				else {
+					pxl[0] = 255 - is;
+					pxl[1] = 255 - is;
+					pxl[2] = 255;
+				}
+			}
+		}
+		img.SaveFile("out\\signed_df2.png", wxBITMAP_TYPE_PNG);
+	}
+#endif
+}
+
+float EdgeGrid::Grid::signed_distance_bilinear(const Point &pt) const
+{
+	coord_t x = pt.x - m_bbox.min.x;
+	coord_t y = pt.y - m_bbox.min.y;
+	coord_t w = m_resolution * m_cols;
+	coord_t h = m_resolution * m_rows;
+	bool    clamped = false;
+	coord_t xcl = x;
+	coord_t ycl = y;
+	if (x < 0) {
+		xcl = 0;
+		clamped = true;
+	} else if (x >= w) {
+		xcl = w - 1;
+		clamped = true;
+	}
+	if (y < 0) {
+		ycl = 0;
+		clamped = true;
+	} else if (y >= h) {
+		ycl = h - 1;
+		clamped = true;
+	}
+ 
+	coord_t cell_c = coord_t(floor(xcl / m_resolution));
+	coord_t cell_r = coord_t(floor(ycl / m_resolution));
+	float   tx = float(xcl - cell_c * m_resolution) / float(m_resolution);
+	assert(tx >= -1e-5 && tx < 1.f + 1e-5);
+	float   ty = float(ycl - cell_r * m_resolution) / float(m_resolution);
+	assert(ty >= -1e-5 && ty < 1.f + 1e-5);
+	size_t  addr = cell_r * (m_cols + 1) + cell_c;
+	float   f00 = m_signed_distance_field[addr];
+	float   f01 = m_signed_distance_field[addr+1];
+	addr += m_cols + 1;
+	float   f10 = m_signed_distance_field[addr];
+	float   f11 = m_signed_distance_field[addr+1];
+	float   f0  = (1.f - tx) * f00 + tx * f01;
+	float   f1  = (1.f - tx) * f10 + tx * f11;
+	float	f   = (1.f - ty) * f0 + ty * f1;
+
+	if (clamped) {
+		if (f > 0) {
+			if (x < 0)
+				f += -x;
+			else if (x >= w)
+				f += x - w + 1;
+			if (y < 0)
+				f += -y;
+			else if (y >= h)
+				f += y - h + 1;
+		} else {			
+			if (x < 0)
+				f -= -x;
+			else if (x >= w)
+				f -= x - w + 1;
+			if (y < 0)
+				f -= -y;
+			else if (y >= h)
+				f -= y - h + 1;
+		}
+	}
+
+	return f;
+}
+
+bool EdgeGrid::Grid::signed_distance_edges(const Point &pt, coord_t search_radius, coordf_t &result_min_dist, bool *pon_segment) const {
+	BoundingBox bbox;
+	bbox.min = bbox.max = Point(pt.x - m_bbox.min.x, pt.y - m_bbox.min.y);
+	bbox.defined = true;
+	// Upper boundary, round to grid and test validity.
+	bbox.max.x += search_radius;
+	bbox.max.y += search_radius;
+	if (bbox.max.x < 0 || bbox.max.y < 0)
+		return false;
+	bbox.max.x /= m_resolution;
+	bbox.max.y /= m_resolution;
+	if (bbox.max.x >= m_cols)
+		bbox.max.x = m_cols - 1;
+	if (bbox.max.y >= m_rows)
+		bbox.max.y = m_rows - 1;
+	// Lower boundary, round to grid and test validity.
+	bbox.min.x -= search_radius;
+	bbox.min.y -= search_radius;
+	if (bbox.min.x < 0)
+		bbox.min.x = 0;
+	if (bbox.min.y < 0)
+		bbox.min.y = 0;
+	bbox.min.x /= m_resolution;
+	bbox.min.y /= m_resolution;
+	// Is the interval empty?
+	if (bbox.min.x > bbox.max.x ||
+		bbox.min.y > bbox.max.y)
+		return false;
+	// Traverse all cells in the bounding box.
+	float d_min = search_radius;
+	// Signum of the distance field at pt.
+	int sign_min = 0;
+	bool on_segment = false;
+	for (int r = bbox.min.y; r <= bbox.max.y; ++ r) {
+		for (int c = bbox.min.x; c <= bbox.max.x; ++ c) {
+			const Cell &cell = m_cells[r * m_cols + c];
+			for (size_t i = cell.begin; i < cell.end; ++ i) {
+				const Slic3r::Points &pts = *m_contours[m_cell_data[i].first];
+				size_t ipt = m_cell_data[i].second;
+				// End points of the line segment.
+				const Slic3r::Point &p1 = pts[ipt];
+				const Slic3r::Point &p2 = pts[(ipt + 1 == pts.size()) ? 0 : ipt + 1];
+				Slic3r::Point v_seg = p1.vector_to(p2);
+				Slic3r::Point v_pt = p1.vector_to(pt);
+				// dot(p2-p1, pt-p1)
+				int64_t t_pt = int64_t(v_seg.x) * int64_t(v_pt.x) + int64_t(v_seg.y) * int64_t(v_pt.y);
+				// l2 of seg
+				int64_t l2_seg = int64_t(v_seg.x) * int64_t(v_seg.x) + int64_t(v_seg.y) * int64_t(v_seg.y);
+				if (t_pt < 0) {
+					// Closest to p1.
+					double dabs = sqrt(int64_t(v_pt.x) * int64_t(v_pt.x) + int64_t(v_pt.y) * int64_t(v_pt.y));
+					if (dabs < d_min) {
+						// Previous point.
+						const Slic3r::Point &p0 = pts[(ipt == 0) ? (pts.size() - 1) : ipt - 1];
+						Slic3r::Point v_seg_prev = p0.vector_to(p1);
+						int64_t t2_pt = int64_t(v_seg_prev.x) * int64_t(v_pt.x) + int64_t(v_seg_prev.y) * int64_t(v_pt.y);
+						if (t2_pt > 0) {
+							// Inside the wedge between the previous and the next segment.
+							d_min = dabs;
+							// Set the signum depending on whether the vertex is convex or reflex.
+							int64_t det = int64_t(v_seg_prev.x) * int64_t(v_seg.y) - int64_t(v_seg_prev.y) * int64_t(v_seg.x);
+							assert(det != 0);
+							sign_min = (det > 0) ? 1 : -1;
+							on_segment = false;
+						}
+					}
+				}
+				else if (t_pt > l2_seg) {
+					// Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the same cell.
+					continue;
+				} else {
+					// Closest to the segment.
+					assert(t_pt >= 0 && t_pt <= l2_seg);
+					int64_t d_seg = int64_t(v_seg.y) * int64_t(v_pt.x) - int64_t(v_seg.x) * int64_t(v_pt.y);
+					double d = double(d_seg) / sqrt(double(l2_seg));
+					double dabs = std::abs(d);
+					if (dabs < d_min) {
+						d_min = dabs;
+						sign_min = (d_seg < 0) ? -1 : ((d_seg == 0) ? 0 : 1);
+						on_segment = true;
+					}
+				}
+			}
+		}
+	}
+	if (d_min >= search_radius)
+		return false;
+	result_min_dist = d_min * sign_min;
+	if (pon_segment != NULL)
+		*pon_segment = on_segment;
+	return true;
+}
+
+bool EdgeGrid::Grid::signed_distance(const Point &pt, coord_t search_radius, coordf_t &result_min_dist) const
+{
+	if (signed_distance_edges(pt, search_radius, result_min_dist))
+		return true;
+	if (m_signed_distance_field.empty())
+		return false;
+	result_min_dist = signed_distance_bilinear(pt);
+	return true;
+}
+
+void EdgeGrid::save_png(const EdgeGrid::Grid &grid, const BoundingBox &bbox, coord_t resolution, const char *path)
+{
+	unsigned int w = (bbox.max.x - bbox.min.x + resolution - 1) / resolution;
+	unsigned int h = (bbox.max.y - bbox.min.y + resolution - 1) / resolution;
+	wxImage img(w, h);
+    unsigned char *data = img.GetData();
+    memset(data, 0, w * h * 3);
+
+	static int iRun = 0;
+	++iRun;
+    
+    const coord_t search_radius = grid.resolution() * 2;
+	const coord_t display_blend_radius = grid.resolution() * 5;
+	for (coord_t r = 0; r < h; ++r) {
+    	for (coord_t c = 0; c < w; ++ c) {
+			unsigned char *pxl = data + (((h - r - 1) * w) + c) * 3;
+			Point pt(c * resolution + bbox.min.x, r * resolution + bbox.min.y);
+			coordf_t min_dist;
+			bool on_segment;
+//			if (grid.signed_distance_edges(pt, search_radius, min_dist, &on_segment)) {
+			if (grid.signed_distance(pt, search_radius, min_dist)) {
+				//FIXME
+				on_segment = true;
+				float s = 255 * std::abs(min_dist) / float(display_blend_radius);
+				int is = std::max(0, std::min(255, int(floor(s + 0.5f))));
+				if (min_dist < 0) {
+					if (on_segment) {
+						pxl[0] = 255;
+						pxl[1] = 255 - is;
+						pxl[2] = 255 - is;
+					} else {
+						pxl[0] = 128;
+						pxl[1] = 128;
+						pxl[2] = 255 - is;						
+					}
+				}
+				else {
+					if (on_segment) {
+						pxl[0] = 255 - is;
+						pxl[1] = 255 - is;
+						pxl[2] = 255;
+					} else {
+						pxl[0] = 255 - is;
+						pxl[1] = 0;
+						pxl[2] = 255;
+					}
+				}
+			} else {
+				pxl[0] = 0;
+				pxl[1] = 255;
+				pxl[2] = 0;
+			}
+
+			float gridx = float(pt.x - grid.bbox().min.x) / float(grid.resolution());
+			float gridy = float(pt.y - grid.bbox().min.y) / float(grid.resolution());
+			if (gridx >= -0.4f && gridy >= -0.4f && gridx <= grid.cols() + 0.4f && gridy <= grid.rows() + 0.4f) {
+				int ix = int(floor(gridx + 0.5f));
+				int iy = int(floor(gridy + 0.5f));
+				float dx = gridx - float(ix);
+				float dy = gridy - float(iy);
+				float d = sqrt(dx*dx + dy*dy) * float(grid.resolution()) / float(resolution);
+				if (d < 1.f) {
+					// Less than 1 pixel from the grid point.
+					float t = 0.5f + 0.5f * d;
+					pxl[0] = (unsigned char)(t * pxl[0]);
+					pxl[1] = (unsigned char)(t * pxl[1]);
+					pxl[2] = (unsigned char)(t * pxl[2]);
+				}
+			}
+
+			float dgrid = fabs(min_dist) / float(grid.resolution());
+			float igrid = floor(dgrid + 0.5f);
+			dgrid = std::abs(dgrid - igrid) * float(grid.resolution()) / float(resolution);
+			if (dgrid < 1.f) {
+				// Less than 1 pixel from the grid point.
+				float t = 0.5f + 0.5f * dgrid;
+				pxl[0] = (unsigned char)(t * pxl[0]);
+				pxl[1] = (unsigned char)(t * pxl[1]);
+				pxl[2] = (unsigned char)(t * pxl[2]);
+				if (igrid > 0.f) {
+					// Other than zero iso contour.
+					int g = pxl[1] + 255.f * (1.f - t);
+					pxl[1] = std::min(g, 255);
+				}
+			}
+		}
+    }
+
+    img.SaveFile(path, wxBITMAP_TYPE_PNG);
+}
+
+} // namespace Slic3r

+ 81 - 0
xs/src/libslic3r/EdgeGrid.hpp

@@ -0,0 +1,81 @@
+#ifndef slic3r_EdgeGrid_hpp_
+#define slic3r_EdgeGrid_hpp_
+
+#include <stdint.h>
+#include <math.h>
+
+#include "Point.hpp"
+#include "BoundingBox.hpp"
+#include "ExPolygon.hpp"
+#include "ExPolygonCollection.hpp"
+
+namespace Slic3r {
+namespace EdgeGrid {
+
+struct Grid
+{
+	Grid();
+	~Grid();
+
+	void create(const Polygons &polygons, coord_t resolution);
+	void create(const ExPolygon &expoly, coord_t resolution);
+	void create(const ExPolygons &expolygons, coord_t resolution);
+	void create(const ExPolygonCollection &expolygons, coord_t resolution);
+
+	// Fill in a rough m_signed_distance_field from the edge grid.
+	// The rough SDF is used by signed_distance() for distances outside of the search_radius.
+	void calculate_sdf();
+
+	// Return an estimate of the signed distance based on m_signed_distance_field grid.
+	float signed_distance_bilinear(const Point &pt) const;
+
+	// Calculate a signed distance to the contours in search_radius from the point.
+	bool signed_distance_edges(const Point &pt, coord_t search_radius, coordf_t &result_min_dist, bool *pon_segment = NULL) const;
+
+	// Calculate a signed distance to the contours in search_radius from the point. If no edge is found in search_radius,
+	// return an interpolated value from m_signed_distance_field, if it exists.
+	bool signed_distance(const Point &pt, coord_t search_radius, coordf_t &result_min_dist) const;
+
+	const BoundingBox& 	bbox() const { return m_bbox; }
+	const coord_t 		resolution() const { return m_resolution; }
+	const size_t		rows() const { return m_rows; }
+	const size_t		cols() const { return m_cols; }
+
+protected:
+	void create_from_m_contours(coord_t resolution);
+
+	// Bounding box around the contours.
+	BoundingBox 								m_bbox;
+	// Grid dimensions.
+	coord_t										m_resolution;
+	size_t										m_rows;
+	size_t										m_cols;
+
+	// Referencing the source contours.
+	// This format allows one to work with any Slic3r fixed point contour format
+	// (Polygon, ExPolygon, ExPolygonCollection etc).
+	std::vector<const Slic3r::Points*>			m_contours;
+
+	// Referencing a contour and a line segment of m_contours.
+	std::vector<std::pair<size_t, size_t>>		m_cell_data;
+
+	struct Cell {
+		Cell() : begin(0), end(0) {}
+		size_t begin;
+		size_t end;
+	};
+	// Full grid of cells.
+	std::vector<Cell> 							m_cells;
+
+	// Distance field derived from the edge grid, seed filled by the Danielsson chamfer metric.
+	// May be empty.
+	std::vector<float>							m_signed_distance_field;
+};
+
+// Debugging utility. Save the signed distance field.
+extern void save_png(const Grid &grid, const BoundingBox &bbox, coord_t resolution, const char *path);
+
+} // namespace EdgeGrid
+} // namespace Slic3r
+
+#endif /* slic3r_EdgeGrid_hpp_ */

+ 369 - 58
xs/src/libslic3r/GCode.cpp

@@ -1,9 +1,21 @@
 #include "GCode.hpp"
 #include "ExtrusionEntity.hpp"
+#include "EdgeGrid.hpp"
 #include <algorithm>
 #include <cstdlib>
 #include <math.h>
 
+#include "SVG.hpp"
+
+#if 0
+// Enable debugging and asserts, even in the release build.
+#define DEBUG
+#define _DEBUG
+#undef NDEBUG
+#endif
+
+#include <assert.h>
+
 namespace Slic3r {
 
 AvoidCrossingPerimeters::AvoidCrossingPerimeters()
@@ -198,10 +210,20 @@ Wipe::wipe(GCode &gcodegen, bool toolchange)
 #define EXTRUDER_CONFIG(OPT) this->config.OPT.get_at(this->writer.extruder()->id)
 
 GCode::GCode()
-    : placeholder_parser(NULL), enable_loop_clipping(true), enable_cooling_markers(false), layer_count(0),
+    : placeholder_parser(NULL), enable_loop_clipping(true), 
+        enable_cooling_markers(false), enable_extrusion_role_markers(false), 
+        layer_count(0),
         layer_index(-1), layer(NULL), first_layer(false), elapsed_time(0.0), volumetric_speed(0),
-        _last_pos_defined(false)
+        _last_pos_defined(false),
+        _lower_layer_edge_grid(NULL),
+        _last_extrusion_role(erNone)
+{
+}
+
+GCode::~GCode()
 {
+    delete _lower_layer_edge_grid;
+    _lower_layer_edge_grid = NULL;
 }
 
 const Point&
@@ -279,6 +301,8 @@ GCode::change_layer(const Layer &layer)
     this->layer = &layer;
     this->layer_index++;
     this->first_layer = (layer.id() == 0);
+    delete this->_lower_layer_edge_grid;
+    this->_lower_layer_edge_grid = NULL;
     
     // avoid computing islands and overhangs if they're not needed
     if (this->config.avoid_crossing_perimeters) {
@@ -308,12 +332,230 @@ GCode::change_layer(const Layer &layer)
     return gcode;
 }
 
+static inline const char* ExtrusionRole2String(const ExtrusionRole role)
+{
+    switch (role) {
+    case erNone:                        return "erNone";
+    case erPerimeter:                   return "erPerimeter";
+    case erExternalPerimeter:           return "erExternalPerimeter";
+    case erOverhangPerimeter:           return "erOverhangPerimeter";
+    case erInternalInfill:              return "erInternalInfill";
+    case erSolidInfill:                 return "erSolidInfill";
+    case erTopSolidInfill:              return "erTopSolidInfill";
+    case erBridgeInfill:                return "erBridgeInfill";
+    case erGapFill:                     return "erGapFill";
+    case erSkirt:                       return "erSkirt";
+    case erSupportMaterial:             return "erSupportMaterial";
+    case erSupportMaterialInterface:    return "erSupportMaterialInterface";
+    default:                            return "erInvalid";
+    };
+}
+
+static inline const char* ExtrusionLoopRole2String(const ExtrusionLoopRole role)
+{
+    switch (role) {
+    case elrDefault:                    return "elrDefault";
+    case elrContourInternalPerimeter:   return "elrContourInternalPerimeter";
+    case elrSkirt:                      return "elrSkirt";
+    default:                            return "elrInvalid";
+    }
+};
+
+// Return a value in <0, 1> of a cubic B-spline kernel centered around zero.
+// The B-spline is re-scaled so it has value 1 at zero.
+static inline float bspline_kernel(float x)
+{
+    x = std::abs(x);
+	if (x < 1.f) {
+		return 1.f - (3. / 2.) * x * x + (3.f / 4.f) * x * x * x;
+	}
+	else if (x < 2.f) {
+		x -= 1.f;
+		float x2 = x * x;
+		float x3 = x2 * x;
+		return (1.f / 4.f) - (3.f / 4.f) * x + (3.f / 4.f) * x2 - (1.f / 4.f) * x3;
+	}
+	else
+        return 0;
+}
+
+static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance)
+{
+    // The extrudate is not fully supported by the lower layer. Fit a polynomial penalty curve.
+    // Solved by sympy package:
+/*
+from sympy import *
+(x,a,b,c,d,r,z)=symbols('x a b c d r z')
+p = a + b*x + c*x*x + d*x*x*x
+p2 = p.subs(solve([p.subs(x, -r), p.diff(x).subs(x, -r), p.diff(x,x).subs(x, -r), p.subs(x, 0)-z], [a, b, c, d]))
+from sympy.plotting import plot
+plot(p2.subs(r,0.2).subs(z,1.), (x, -1, 3), adaptive=False, nb_of_points=400)
+*/
+    if (overlap_distance < - nozzle_r) {
+        // The extrudate is fully supported by the lower layer. This is the ideal case, therefore zero penalty.
+        return 0.f;
+    } else {
+        float x  = overlap_distance / nozzle_r;
+        float x2 = x * x;
+        float x3 = x2 * x;
+        return weight_zero * (1.f + 3.f * x + 3.f * x2 + x3);
+    }
+}
+
+static Points::iterator project_point_to_polygon_and_insert(Polygon &polygon, const Point &pt, double eps)
+{
+    assert(polygon.points.size() >= 2);
+    if (polygon.points.size() <= 1)
+    if (polygon.points.size() == 1)
+        return polygon.points.begin();
+
+    Point  pt_min;
+    double d_min = std::numeric_limits<double>::max();
+    size_t i_min = size_t(-1);
+
+    for (size_t i = 0; i < polygon.points.size(); ++ i) {
+        size_t j = i + 1;
+        if (j == polygon.points.size())
+            j = 0;
+        const Point &p1 = polygon.points[i];
+        const Point &p2 = polygon.points[j];
+        const Slic3r::Point v_seg = p1.vector_to(p2);
+        const Slic3r::Point v_pt  = p1.vector_to(pt);
+        const int64_t l2_seg = int64_t(v_seg.x) * int64_t(v_seg.x) + int64_t(v_seg.y) * int64_t(v_seg.y);
+        int64_t t_pt = int64_t(v_seg.x) * int64_t(v_pt.x) + int64_t(v_seg.y) * int64_t(v_pt.y);
+        if (t_pt < 0) {
+            // Closest to p1.
+            double dabs = sqrt(int64_t(v_pt.x) * int64_t(v_pt.x) + int64_t(v_pt.y) * int64_t(v_pt.y));
+            if (dabs < d_min) {
+                d_min  = dabs;
+                i_min  = i;
+                pt_min = p1;
+            }
+        }
+        else if (t_pt > l2_seg) {
+            // Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the next step.
+            continue;
+        } else {
+            // Closest to the segment.
+            assert(t_pt >= 0 && t_pt <= l2_seg);
+            int64_t d_seg = int64_t(v_seg.y) * int64_t(v_pt.x) - int64_t(v_seg.x) * int64_t(v_pt.y);
+            double d = double(d_seg) / sqrt(double(l2_seg));
+            double dabs = std::abs(d);
+            if (dabs < d_min) {
+                d_min  = dabs;
+                i_min  = i;
+                // Evaluate the foot point.
+                pt_min = p1;
+                double linv = double(d_seg) / double(l2_seg);
+                pt_min.x = pt.x - coord_t(floor(double(v_seg.y) * linv + 0.5));
+				pt_min.y = pt.y + coord_t(floor(double(v_seg.x) * linv + 0.5));
+				assert(Line(p1, p2).distance_to(pt_min) < scale_(1e-5));
+            }
+        }
+    }
+
+	assert(i_min != size_t(-1));
+    if (pt_min.distance_to(polygon.points[i_min]) > eps) {
+        // Insert a new point on the segment i_min, i_min+1.
+        return polygon.points.insert(polygon.points.begin() + (i_min + 1), pt_min);
+    }
+    return polygon.points.begin() + i_min;
+}
+
+std::vector<float> polygon_parameter_by_length(const Polygon &polygon)
+{
+    // Parametrize the polygon by its length.
+    std::vector<float> lengths(polygon.points.size()+1, 0.);
+    for (size_t i = 1; i < polygon.points.size(); ++ i)
+        lengths[i] = lengths[i-1] + polygon.points[i].distance_to(polygon.points[i-1]);
+    lengths.back() = lengths[lengths.size()-2] + polygon.points.front().distance_to(polygon.points.back());
+    return lengths;
+}
+
+std::vector<float> polygon_angles_at_vertices(const Polygon &polygon, const std::vector<float> &lengths, float min_arm_length)
+{
+    assert(polygon.points.size() + 1 == lengths.size());
+    if (min_arm_length > 0.25f * lengths.back())
+        min_arm_length = 0.25f * lengths.back();
+
+    // Find the initial prev / next point span.
+    size_t idx_prev = polygon.points.size();
+    size_t idx_curr = 0;
+    size_t idx_next = 1;
+    while (idx_prev > idx_curr && lengths.back() - lengths[idx_prev] < min_arm_length)
+        -- idx_prev;
+    while (idx_next < idx_prev && lengths[idx_next] < min_arm_length)
+        ++ idx_next;
+
+    std::vector<float> angles(polygon.points.size(), 0.f);
+    for (; idx_curr < polygon.points.size(); ++ idx_curr) {
+        // Move idx_prev up until the distance between idx_prev and idx_curr is lower than min_arm_length.
+        if (idx_prev >= idx_curr) {
+            while (idx_prev < polygon.points.size() && lengths.back() - lengths[idx_prev] + lengths[idx_curr] > min_arm_length)
+                ++ idx_prev;
+            if (idx_prev == polygon.points.size())
+                idx_prev = 0;
+        }
+        while (idx_prev < idx_curr && lengths[idx_curr] - lengths[idx_prev] > min_arm_length)
+            ++ idx_prev;
+        // Move idx_prev one step back.
+        if (idx_prev == 0)
+            idx_prev = polygon.points.size() - 1;
+        else
+            -- idx_prev;
+        // Move idx_next up until the distance between idx_curr and idx_next is greater than min_arm_length.
+        if (idx_curr <= idx_next) {
+            while (idx_next < polygon.points.size() && lengths[idx_next] - lengths[idx_curr] < min_arm_length)
+                ++ idx_next;
+            if (idx_next == polygon.points.size())
+                idx_next = 0;
+        }
+        while (idx_next < idx_curr && lengths.back() - lengths[idx_curr] + lengths[idx_next] < min_arm_length)
+            ++ idx_next;
+        // Calculate angle between idx_prev, idx_curr, idx_next.
+        const Point &p0 = polygon.points[idx_prev];
+        const Point &p1 = polygon.points[idx_curr];
+        const Point &p2 = polygon.points[idx_next];
+        const Point  v1 = p0.vector_to(p1);
+        const Point  v2 = p1.vector_to(p2);
+		int64_t dot   = int64_t(v1.x)*int64_t(v2.x) + int64_t(v1.y)*int64_t(v2.y);
+		int64_t cross = int64_t(v1.x)*int64_t(v2.y) - int64_t(v1.y)*int64_t(v2.x);
+		float angle = float(atan2(double(cross), double(dot)));
+        angles[idx_curr] = angle;
+    }
+
+    return angles;
+}
+
 std::string
 GCode::extrude(ExtrusionLoop loop, std::string description, double speed)
 {
     // get a copy; don't modify the orientation of the original loop object otherwise
     // next copies (if any) would not detect the correct orientation
-    
+
+    if (this->layer->lower_layer != NULL) {
+        if (this->_lower_layer_edge_grid == NULL) {
+            // Create the distance field for a layer below.
+            const coord_t distance_field_resolution = scale_(1.f);
+            this->_lower_layer_edge_grid = new EdgeGrid::Grid();
+            this->_lower_layer_edge_grid->create(this->layer->lower_layer->slices, distance_field_resolution);
+            this->_lower_layer_edge_grid->calculate_sdf();
+            #if 0
+            {
+                static int iRun = 0;
+                char path[2048];
+                sprintf(path, "out\\GCode_extrude_loop_edge_grid-%d.png", iRun++);
+                BoundingBox bbox = this->_lower_layer_edge_grid->bbox();
+                bbox.min.x -= scale_(5.f);
+                bbox.min.y -= scale_(5.f);
+                bbox.max.x += scale_(5.f);
+                bbox.max.y += scale_(5.f);
+                EdgeGrid::save_png(*this->_lower_layer_edge_grid, bbox, scale_(0.1f), path);
+            }
+            #endif
+        }
+    }
+  
     // extrude all loops ccw
     bool was_clockwise = loop.make_counter_clockwise();
     
@@ -326,68 +568,126 @@ GCode::extrude(ExtrusionLoop loop, std::string description, double speed)
     if (this->config.spiral_vase) {
         loop.split_at(last_pos);
     } else if (seam_position == spNearest || seam_position == spAligned) {
-        const Polygon polygon = loop.polygon();
-        
-        // simplify polygon in order to skip false positives in concave/convex detection
-        // (loop is always ccw as polygon.simplify() only works on ccw polygons)
-        Polygons simplified = polygon.simplify(scale_(EXTRUDER_CONFIG(nozzle_diameter))/2);
-        
-        // restore original winding order so that concave and convex detection always happens
-        // on the right/outer side of the polygon
-        if (was_clockwise) {
-            for (Polygons::iterator p = simplified.begin(); p != simplified.end(); ++p)
-                p->reverse();
-        }
-        
-        // concave vertices have priority
-        Points candidates;
-        for (Polygons::const_iterator p = simplified.begin(); p != simplified.end(); ++p) {
-            Points concave = p->concave_points(PI*4/3);
-            candidates.insert(candidates.end(), concave.begin(), concave.end());
+        Polygon        polygon    = loop.polygon();
+        const coordf_t nozzle_dmr = EXTRUDER_CONFIG(nozzle_diameter);
+        const coord_t  nozzle_r   = scale_(0.5*nozzle_dmr);
+
+        // Retrieve the last start position for this object.
+        float last_pos_weight = 1.f;
+        if (seam_position == spAligned && this->layer != NULL && this->_seam_position.count(this->layer->object()) > 0) {
+            last_pos = this->_seam_position[this->layer->object()];
+            last_pos_weight = 5.f;
         }
-        
-        // if no concave points were found, look for convex vertices
-        if (candidates.empty()) {
-            for (Polygons::const_iterator p = simplified.begin(); p != simplified.end(); ++p) {
-                Points convex = p->convex_points(PI*2/3);
-                candidates.insert(candidates.end(), convex.begin(), convex.end());
+
+        // Insert a projection of last_pos into the polygon.
+		size_t last_pos_proj_idx;
+		{
+			auto it = project_point_to_polygon_and_insert(polygon, last_pos, 0.1 * nozzle_r);
+			last_pos_proj_idx = it - polygon.points.begin();
+		}
+        Point last_pos_proj = polygon.points[last_pos_proj_idx];
+        // Parametrize the polygon by its length.
+        std::vector<float> lengths = polygon_parameter_by_length(polygon);
+
+        // For each polygon point, store a penalty.
+        // First calculate the angles, store them as penalties. The angles are caluculated over a minimum arm length of nozzle_r.
+        std::vector<float> penalties = polygon_angles_at_vertices(polygon, lengths, nozzle_r);
+        // No penalty for reflex points, slight penalty for convex points, high penalty for flat surfaces.
+        const float penaltyConvexVertex = 1.f;
+        const float penaltyFlatSurface  = 5.f;
+        const float penaltySeam         = 1.3f;
+        const float penaltyOverhangHalf = 10.f;
+        // Penalty for visible seams.
+        for (size_t i = 0; i < polygon.points.size(); ++ i) {
+            float ccwAngle = penalties[i];
+            if (was_clockwise)
+                ccwAngle = - ccwAngle;
+            float penalty = 0;
+//            if (ccwAngle <- float(PI/3.))
+            if (ccwAngle <- float(0.6 * PI))
+                // Sharp reflex vertex. We love that, it hides the seam perfectly.
+                penalty = 0.f;
+//            else if (ccwAngle > float(PI/3.))
+            else if (ccwAngle > float(0.6 * PI))
+                // Seams on sharp convex vertices are more visible than on reflex vertices.
+                penalty = penaltyConvexVertex;
+            else if (ccwAngle < 0.f) {
+                // Interpolate penalty between maximum and zero.
+                penalty = penaltyFlatSurface * bspline_kernel(ccwAngle * (PI * 2. / 3.));
+            } else {
+                assert(ccwAngle >= 0.f);
+                // Interpolate penalty between maximum and the penalty for a convex vertex.
+                penalty = penaltyConvexVertex + (penaltyFlatSurface - penaltyConvexVertex) * bspline_kernel(ccwAngle * (PI * 2. / 3.));
             }
+            // Give a negative penalty for points close to the last point or the prefered seam location.
+            //float dist_to_last_pos_proj = last_pos_proj.distance_to(polygon.points[i]);
+            float dist_to_last_pos_proj = (i < last_pos_proj_idx) ? 
+                std::min(lengths[last_pos_proj_idx] - lengths[i], lengths.back() - lengths[last_pos_proj_idx] + lengths[i]) : 
+                std::min(lengths[i] - lengths[last_pos_proj_idx], lengths.back() - lengths[i] + lengths[last_pos_proj_idx]);
+            float dist_max = 0.1f * lengths.back(); // 5.f * nozzle_dmr
+            penalty -= last_pos_weight * bspline_kernel(dist_to_last_pos_proj / dist_max);
+            penalties[i] = std::max(0.f, penalty);
         }
-        
-        // retrieve the last start position for this object
-        if (this->layer != NULL && this->_seam_position.count(this->layer->object()) > 0) {
-            last_pos = this->_seam_position[this->layer->object()];
-        }
-        
-        Point point;
-        if (seam_position == spNearest) {
-            if (candidates.empty()) candidates = polygon.points;
-            last_pos.nearest_point(candidates, &point);
-            
-            // On 32-bit Linux, Clipper will change some point coordinates by 1 unit
-            // while performing simplify_polygons(), thus split_at_vertex() won't 
-            // find them anymore.
-            if (!loop.split_at_vertex(point)) loop.split_at(point);
-        } else if (!candidates.empty()) {
-            Points non_overhang;
-            for (Points::const_iterator p = candidates.begin(); p != candidates.end(); ++p) {
-                if (!loop.has_overhang_point(*p))
-                    non_overhang.push_back(*p);
+
+        // Penalty for overhangs.
+        if (this->_lower_layer_edge_grid) {
+            // Use the edge grid distance field structure over the lower layer to calculate overhangs.
+            coord_t nozzle_r = scale_(0.5*nozzle_dmr);
+            coord_t search_r = scale_(0.8*nozzle_dmr);
+            for (size_t i = 0; i < polygon.points.size(); ++ i) {
+                const Point &p = polygon.points[i];
+                coordf_t dist;
+                // Signed distance is positive outside the object, negative inside the object.
+                // The point is considered at an overhang, if it is more than nozzle radius
+                // outside of the lower layer contour.
+                bool found = this->_lower_layer_edge_grid->signed_distance(p, search_r, dist);
+                // If the approximate Signed Distance Field was initialized over this->_lower_layer_edge_grid,
+                // then the signed distnace shall always be known.
+                assert(found);
+                penalties[i] += extrudate_overlap_penalty(nozzle_r, penaltyOverhangHalf, dist);
             }
-            
-            if (!non_overhang.empty())
-                candidates = non_overhang;
-            
-            last_pos.nearest_point(candidates, &point);
-            if (!loop.split_at_vertex(point)) loop.split_at(point);  // see note above
-        } else {
-            point = last_pos.projection_onto(polygon);
-            loop.split_at(point);
         }
-        if (this->layer != NULL)
-            this->_seam_position[this->layer->object()] = point;
+
+        // Find a point with a minimum penalty.
+        size_t idx_min = std::min_element(penalties.begin(), penalties.end()) - penalties.begin();
+
+        // Export the contour into a SVG file.
+        #if 0
+        {
+            static int iRun = 0;
+            char path[2048];
+            sprintf(path, "out\\GCode_extrude_loop-%d.svg", iRun ++);
+            SVG svg(path);
+			if (this->layer->lower_layer != NULL)
+				svg.draw(this->layer->lower_layer->slices.expolygons);
+            for (size_t i = 0; i < loop.paths.size(); ++ i)
+                svg.draw(loop.paths[i].as_polyline(), "red");
+            Polylines polylines;
+            for (size_t i = 0; i < loop.paths.size(); ++ i)
+                polylines.push_back(loop.paths[i].as_polyline());
+            Slic3r::Polygons polygons;
+            coordf_t nozzle_dmr = EXTRUDER_CONFIG(nozzle_diameter);
+            coord_t delta = scale_(0.5*nozzle_dmr);
+            Slic3r::offset(polylines, &polygons, delta);
+//            for (size_t i = 0; i < polygons.size(); ++ i) svg.draw((Polyline)polygons[i], "blue");
+			svg.draw(last_pos, "green", 3);
+			svg.draw(polygon.points[idx_min], "yellow", 3);
+			svg.Close();
+        }
+        #endif
+
+        // Split the loop at the point with a minium penalty.
+        if (!loop.split_at_vertex(polygon.points[idx_min]))
+            // The point is not in the original loop. Insert it.
+            loop.split_at(polygon.points[idx_min]);
+
     } else if (seam_position == spRandom) {
         if (loop.role == elrContourInternalPerimeter) {
+            // This loop does not contain any other loop. Set a random position.
+            // The other loops will get a seam close to the random point chosen
+            // on the inner most contour.
+            //FIXME This works correctly for inner contours first only.
+            //FIXME Better parametrize the loop by its length.
             Polygon polygon = loop.polygon();
             Point centroid = polygon.centroid();
             last_pos = Point(polygon.bounding_box().max.x, centroid.y);
@@ -416,6 +716,8 @@ GCode::extrude(ExtrusionLoop loop, std::string description, double speed)
     // extrude along the path
     std::string gcode;
     for (ExtrusionPaths::const_iterator path = paths.begin(); path != paths.end(); ++path)
+//    description += ExtrusionLoopRole2String(loop.role);
+//    description += ExtrusionRole2String(path->role);
         gcode += this->_extrude(*path, description, speed);
     
     // reset acceleration
@@ -477,6 +779,7 @@ GCode::extrude(const ExtrusionEntity &entity, std::string description, double sp
 std::string
 GCode::extrude(const ExtrusionPath &path, std::string description, double speed)
 {
+//    description += ExtrusionRole2String(path.role);
     std::string gcode = this->_extrude(path, description, speed);
     
     // reset acceleration
@@ -561,6 +864,14 @@ GCode::_extrude(ExtrusionPath path, std::string description, double speed)
     double F = speed * 60;  // convert mm/sec to mm/min
     
     // extrude arc or line
+    if (this->enable_extrusion_role_markers) {
+        if (path.role != this->_last_extrusion_role) {
+            this->_last_extrusion_role = path.role;
+            char buf[32];
+            sprintf(buf, ";_EXTRUSION_ROLE:%d\n", int(path.role));
+            gcode += buf;
+        }
+    }
     if (path.is_bridge() && this->enable_cooling_markers)
         gcode += ";_BRIDGE_FAN_START\n";
     gcode += this->writer.set_speed(F, "", this->enable_cooling_markers ? ";_EXTRUDE_SET_SPEED" : "");

+ 10 - 0
xs/src/libslic3r/GCode.hpp

@@ -14,7 +14,9 @@
 
 namespace Slic3r {
 
+// Forward declarations.
 class GCode;
+namespace EdgeGrid { class Grid; }
 
 class AvoidCrossingPerimeters {
     public:
@@ -77,15 +79,23 @@ class GCode {
     AvoidCrossingPerimeters avoid_crossing_perimeters;
     bool enable_loop_clipping;
     bool enable_cooling_markers;
+    // Markers for the Pressure Equalizer to recognize the extrusion type.
+    // The Pressure Equalizer removes the markers from the final G-code.
+    bool enable_extrusion_role_markers;
     size_t layer_count;
     int layer_index; // just a counter
     const Layer* layer;
     std::map<const PrintObject*,Point> _seam_position;
+    // Distance Field structure to 
+    EdgeGrid::Grid *_lower_layer_edge_grid;
     bool first_layer; // this flag triggers first layer speeds
     float elapsed_time; // seconds
     double volumetric_speed;
+    // Support for the extrusion role markers. Which marker is active?
+    ExtrusionRole _last_extrusion_role;
     
     GCode();
+    ~GCode();
     const Point& last_pos() const;
     void set_last_pos(const Point &pos);
     bool last_pos_defined() const;

+ 608 - 0
xs/src/libslic3r/GCode/PressureEqualizer.cpp

@@ -0,0 +1,608 @@
+#include <memory.h>
+#include <string.h>
+
+#include "../libslic3r.h"
+#include "../PrintConfig.hpp"
+
+#include "PressureEqualizer.hpp"
+
+namespace Slic3r {
+
+GCodePressureEqualizer::GCodePressureEqualizer(const Slic3r::GCodeConfig *config) : 
+    m_config(config)
+{
+    reset();
+}
+
+GCodePressureEqualizer::~GCodePressureEqualizer()
+{
+}
+
+void GCodePressureEqualizer::reset()
+{
+    circular_buffer_pos     = 0;
+    circular_buffer_size    = 100;
+    circular_buffer_items   = 0;
+    circular_buffer.assign(circular_buffer_size, GCodeLine());
+
+    output_buffer.clear();
+    output_buffer_length    = 0;
+
+    m_current_extruder = 0;
+    // Zero the position of the XYZE axes + the current feed
+    memset(m_current_pos, 0, sizeof(float) * 5);
+    m_current_extrusion_role = erNone;
+    // Expect the first command to fill the nozzle (deretract).
+    m_retracted = true;
+
+    // Calculate filamet crossections for the multiple extruders.
+    m_filament_crossections.clear();
+    for (size_t i = 0; i < m_config->filament_diameter.values.size(); ++ i) {
+        double r = m_config->filament_diameter.values[i];
+        double a = 0.25f*M_PI*r*r;
+        m_filament_crossections.push_back(float(a));
+    }
+
+    m_max_segment_length = 20.f;
+    // Volumetric rate of a 0.45mm x 0.2mm extrusion at 60mm/min XY movement: 0.45*0.2*60*60=5.4*60 = 324 mm^3/min
+    // Volumetric rate of a 0.45mm x 0.2mm extrusion at 20mm/min XY movement: 0.45*0.2*20*60=1.8*60 = 108 mm^3/min
+    // Slope of the volumetric rate, changing from 20mm^3/min to 60mm^3/min over 2 seconds: (5.4-1.8)*60*60/2=60*60*1.8 = 6480 mm^3/min^2
+    m_max_volumetric_extrusion_rate_slope_positive = 6480.f;
+    m_max_volumetric_extrusion_rate_slope_negative = 6480.f;
+
+    for (size_t i = 0; i < numExtrusionRoles; ++ i) {
+        m_max_volumetric_extrusion_rate_slopes[i].negative = m_max_volumetric_extrusion_rate_slope_negative;
+        m_max_volumetric_extrusion_rate_slopes[i].positive = m_max_volumetric_extrusion_rate_slope_positive;
+    }
+
+    // Don't regulate the pressure in infill.
+    m_max_volumetric_extrusion_rate_slopes[erBridgeInfill].negative = 0;
+    m_max_volumetric_extrusion_rate_slopes[erBridgeInfill].positive = 0;
+    // Don't regulate the pressure in gap fill.
+    m_max_volumetric_extrusion_rate_slopes[erGapFill].negative = 0;
+    m_max_volumetric_extrusion_rate_slopes[erGapFill].positive = 0;
+
+    m_stat.reset();
+}
+
+const char* GCodePressureEqualizer::process(const char *szGCode, bool flush)
+{
+    // Reset length of the output_buffer.
+    output_buffer_length = 0;
+
+    if (szGCode != 0) {
+        const char *p = szGCode;
+        while (*p != 0) {
+            // Find end of the line.
+            const char *endl = p;
+            // Slic3r always generates end of lines in a Unix style.
+            for (; *endl != 0 && *endl != '\n'; ++ endl) ;
+            if (circular_buffer_items == circular_buffer_size)
+                // Buffer is full. Push out the oldest line.
+                output_gcode_line(circular_buffer[circular_buffer_pos]);
+            else
+                ++ circular_buffer_items;
+            // Process a G-code line, store it into the provided GCodeLine object.
+            size_t idx_tail = circular_buffer_pos;
+            circular_buffer_pos = circular_buffer_idx_next(circular_buffer_pos);
+            if (! process_line(p, endl - p, circular_buffer[idx_tail])) {
+                // The line has to be forgotten. It contains comment marks, which shall be
+                // filtered out of the target g-code.
+                circular_buffer_pos = idx_tail;
+                -- circular_buffer_items;
+            }
+            p = endl;
+            if (*p == '\n') 
+                ++ p;
+        }
+    }
+
+    if (flush) {
+        // Flush the remaining valid lines of the circular buffer.
+        for (size_t idx = circular_buffer_idx_head(); circular_buffer_items > 0; -- circular_buffer_items) {
+            output_gcode_line(circular_buffer[idx]);
+            if (++ idx == circular_buffer_size)
+                idx = 0;
+        }
+        // Reset the index pointer.
+        assert(circular_buffer_items == 0);
+        circular_buffer_pos = 0;
+
+#if 1 
+        printf("Statistics: \n");
+        printf("Minimum volumetric extrusion rate: %f\n", m_stat.volumetric_extrusion_rate_min);
+        printf("Maximum volumetric extrusion rate: %f\n", m_stat.volumetric_extrusion_rate_max);
+        if (m_stat.extrusion_length > 0)
+            m_stat.volumetric_extrusion_rate_avg /= m_stat.extrusion_length;
+        printf("Average volumetric extrusion rate: %f\n", m_stat.volumetric_extrusion_rate_avg);
+        m_stat.reset();
+#endif
+    } 
+
+    return output_buffer.data();
+}
+
+// Is a white space?
+static inline bool is_ws(const char c) { return c == ' ' || c == '\t'; }
+// Is it an end of line? Consider a comment to be an end of line as well.
+static inline bool is_eol(const char c) { return c == 0 || c == '\r' || c == '\n' || c == ';'; };
+// Is it a white space or end of line?
+static inline bool is_ws_or_eol(const char c) { return is_ws(c) || is_eol(c); };
+
+// Eat whitespaces.
+static void eatws(const char *&line)
+{
+    while (is_ws(*line)) 
+        ++ line;
+}
+
+// Parse an int starting at the current position of a line.
+// If succeeded, the line pointer is advanced.
+static inline int parse_int(const char *&line)
+{
+    char *endptr = NULL;
+    long result = strtol(line, &endptr, 10);
+    if (endptr == NULL || !is_ws_or_eol(*endptr))
+        throw std::runtime_error("GCodePressureEqualizer: Error parsing an int");
+    line = endptr;
+    return int(result);
+};
+
+// Parse an int starting at the current position of a line.
+// If succeeded, the line pointer is advanced.
+static inline float parse_float(const char *&line)
+{
+    char *endptr = NULL;
+    float result = strtof(line, &endptr);
+    if (endptr == NULL || !is_ws_or_eol(*endptr))
+        throw std::runtime_error("GCodePressureEqualizer: Error parsing a float");
+    line = endptr;
+    return result;
+};
+
+#define EXTRUSION_ROLE_TAG ";_EXTRUSION_ROLE:"
+bool GCodePressureEqualizer::process_line(const char *line, const size_t len, GCodeLine &buf)
+{
+    if (strncmp(line, EXTRUSION_ROLE_TAG, strlen(EXTRUSION_ROLE_TAG)) == 0) {
+        line += strlen(EXTRUSION_ROLE_TAG);
+        int role = atoi(line);
+        this->m_current_extrusion_role = ExtrusionRole(role);
+        return false;
+    }
+
+    // Set the type, copy the line to the buffer.
+    buf.type = GCODELINETYPE_OTHER;
+    buf.modified = false;
+    if (buf.raw.size() < len + 1)
+        buf.raw.assign(line, line + len + 1);
+    else
+        memcpy(buf.raw.data(), line, len);
+    buf.raw[len] = 0;
+    buf.raw_length = len;
+
+    memcpy(buf.pos_start, m_current_pos, sizeof(float)*5);
+    memcpy(buf.pos_end, m_current_pos, sizeof(float)*5);
+    memset(buf.pos_provided, 0, 5);
+
+    buf.volumetric_extrusion_rate = 0.f;
+    buf.volumetric_extrusion_rate_start = 0.f;
+    buf.volumetric_extrusion_rate_end = 0.f;
+    buf.max_volumetric_extrusion_rate_slope_positive = 0.f;
+    buf.max_volumetric_extrusion_rate_slope_negative = 0.f;
+	buf.extrusion_role = m_current_extrusion_role;
+
+    // Parse the G-code line, store the result into the buf.
+    switch (toupper(*line ++)) {
+    case 'G': {
+        int gcode = parse_int(line);
+        eatws(line);
+        switch (gcode) {
+        case 0:
+        case 1:
+        {
+            // G0, G1: A FFF 3D printer does not make a difference between the two.
+            float new_pos[5];
+            memcpy(new_pos, m_current_pos, sizeof(float)*5);
+            bool  changed[5] = { false, false, false, false, false };
+            while (!is_eol(*line)) {
+                char axis = toupper(*line++);
+                int  i = -1;
+                switch (axis) {
+                case 'X':
+                case 'Y':
+                case 'Z':
+                    i = axis - 'X';
+                    break;
+                case 'E':
+                    i = 3;
+                    break;
+                case 'F':
+                    i = 4;
+                    break;
+                default:
+                    assert(false);
+                }
+                if (i == -1)
+                    throw std::runtime_error(std::string("GCodePressureEqualizer: Invalid axis for G0/G1: ") + axis);
+                buf.pos_provided[i] = true;
+                new_pos[i] = parse_float(line);
+                if (i == 3 && m_config->use_relative_e_distances.value)
+                    new_pos[i] += m_current_pos[i];
+                changed[i] = new_pos[i] != m_current_pos[i];
+                eatws(line);
+            }
+            if (changed[3]) {
+                // Extrusion, retract or unretract.
+                float diff = new_pos[3] - m_current_pos[3];
+                if (diff < 0) {
+                    buf.type = GCODELINETYPE_RETRACT;
+                    m_retracted = true;
+                } else if (! changed[0] && ! changed[1] && ! changed[2]) {
+                    // assert(m_retracted);
+                    buf.type = GCODELINETYPE_UNRETRACT;
+                    m_retracted = false;
+                } else {
+                    assert(changed[0] || changed[1]);
+                    // Moving in XY plane.
+                    buf.type = GCODELINETYPE_EXTRUDE;
+                    // Calculate the volumetric extrusion rate.
+                    float diff[4];
+                    for (size_t i = 0; i < 4; ++ i)
+                        diff[i] = new_pos[i] - m_current_pos[i];
+                    // volumetric extrusion rate = A_filament * F_xyz * L_e / L_xyz [mm^3/min]
+                    float len2 = diff[0]*diff[0]+diff[1]*diff[1]+diff[2]*diff[2];
+                    float rate = m_filament_crossections[m_current_extruder] * new_pos[4] * sqrt((diff[3]*diff[3])/len2);
+                    buf.volumetric_extrusion_rate       = rate;
+                    buf.volumetric_extrusion_rate_start = rate;
+                    buf.volumetric_extrusion_rate_end   = rate;
+                    m_stat.update(rate, sqrt(len2));
+                    if (rate < 10.f) {
+                    	printf("Extremely low flow rate: %f\n", rate);
+                    }
+                }
+            } else if (changed[0] || changed[1] || changed[2]) {
+                // Moving without extrusion.
+                buf.type = GCODELINETYPE_MOVE;
+            }
+            memcpy(m_current_pos, new_pos, sizeof(float) * 5);
+            break;
+        }
+        case 92:
+        {
+            // G92 : Set Position
+            // Set a logical coordinate position to a new value without actually moving the machine motors.
+            // Which axes to set?
+            bool set = false;
+            while (!is_eol(*line)) {
+                char axis = toupper(*line++);
+                switch (axis) {
+                case 'X':
+                case 'Y':
+                case 'Z':
+                    m_current_pos[axis - 'X'] = (!is_ws_or_eol(*line)) ? parse_float(line) : 0.f;
+                    set = true;
+                    break;
+                case 'E':
+                    m_current_pos[3] = (!is_ws_or_eol(*line)) ? parse_float(line) : 0.f;
+                    set = true;
+                    break;
+                default:
+                    throw std::runtime_error(std::string("GCodePressureEqualizer: Incorrect axis in a G92 G-code: ") + axis);
+                }
+                eatws(line);
+            }
+            assert(set);
+            break;
+        }
+        case 10:
+        case 22:
+            // Firmware retract.
+            buf.type = GCODELINETYPE_RETRACT;
+            m_retracted = true;
+            break;
+        case 11:
+        case 23:
+            // Firmware unretract.
+            buf.type = GCODELINETYPE_UNRETRACT;
+            m_retracted = false;
+            break;
+        default:
+            // Ignore the rest.
+        break;
+        }
+        break;
+    }
+    case 'M': {
+        int mcode = parse_int(line);
+        eatws(line);
+        switch (mcode) {
+        default:
+            // Ignore the rest of the M-codes.
+        break;
+        }
+        break;
+    }
+    case 'T':
+    {
+        // Activate an extruder head.
+        int new_extruder = parse_int(line);
+        if (new_extruder != m_current_extruder) {
+            m_current_extruder = new_extruder;
+            m_retracted = true;
+            buf.type = GCODELINETYPE_TOOL_CHANGE;
+        } else {
+            buf.type = GCODELINETYPE_NOOP;
+        }
+        break;
+    }
+    }
+
+    buf.extruder_id = m_current_extruder;
+    memcpy(buf.pos_end, m_current_pos, sizeof(float)*5);
+
+    adjust_volumetric_rate();
+	return true;
+}
+
+void GCodePressureEqualizer::output_gcode_line(GCodeLine &line)
+{
+    if (! line.modified) {
+        push_to_output(line.raw.data(), line.raw_length, true);
+        return;
+    }
+
+    // The line was modified.
+    // Find the comment.
+    const char *comment = line.raw.data();
+    while (*comment != ';' && *comment != 0) ++comment;
+    if (*comment != ';')
+        comment = NULL;
+    
+    // Emit the line with lowered extrusion rates.
+    float l2 = line.dist_xyz2();
+    float l = sqrt(l2);
+    size_t nSegments = size_t(ceil(l / m_max_segment_length));
+    char text[2048];
+    if (nSegments == 1) {
+        // Just update this segment.
+        push_line_to_output(line, line.feedrate() * line.volumetric_correction_avg(), comment);
+    } else {
+        bool accelerating = line.volumetric_extrusion_rate_start < line.volumetric_extrusion_rate_end;
+        // Update the initial and final feed rate values.
+        line.pos_start[4] = line.volumetric_extrusion_rate_start * line.pos_end[4] / line.volumetric_extrusion_rate;
+        line.pos_end  [4] = line.volumetric_extrusion_rate_end   * line.pos_end[4] / line.volumetric_extrusion_rate;
+        float feed_avg = 0.5f * (line.pos_start[4] + line.pos_end[4]);
+        // Limiting volumetric extrusion rate slope for this segment.
+        float max_volumetric_extrusion_rate_slope = accelerating ?
+            line.max_volumetric_extrusion_rate_slope_positive : line.max_volumetric_extrusion_rate_slope_negative; 
+        // Total time for the segment, corrected for the possibly lowered volumetric feed rate,
+        // if accelerating / decelerating over the complete segment.
+        float t_total = line.dist_xyz() / feed_avg;
+        // Time of the acceleration / deceleration part of the segment, if accelerating / decelerating
+        // with the maximum volumetric extrusion rate slope.
+        float t_acc    = 0.5f * (line.volumetric_extrusion_rate_start + line.volumetric_extrusion_rate_end) / max_volumetric_extrusion_rate_slope;
+        float l_acc    = l;
+        float l_steady = 0.f;
+        if (t_acc < t_total) {
+            // One may achieve higher print speeds if part of the segment is not speed limited.
+            float l_acc    = t_acc * feed_avg;
+            float l_steady = l - l_acc;
+            if (l_steady < 0.5f * m_max_segment_length) {
+                l_acc    = l;
+                l_steady = 0.f;
+            } else
+                nSegments = size_t(ceil(l_acc / m_max_segment_length));
+        }
+        float pos_start[5];
+        float pos_end  [5];
+        float pos_end2 [4];
+        memcpy(pos_start, line.pos_start, sizeof(float)*5);
+        memcpy(pos_end  , line.pos_end  , sizeof(float)*5);
+        if (l_steady > 0.f) {
+            // There will be a steady feed segment emitted.
+            if (accelerating) {
+                // Prepare the final steady feed rate segment.
+                memcpy(pos_end2, pos_end, sizeof(float)*4);
+                float t = l_acc / l;
+                for (int i = 0; i < 4; ++ i) {
+                    pos_end[i] = pos_start[i] + (pos_end[i] - pos_start[i]) * t;
+                    line.pos_provided[i] = true;
+                }
+            } else {
+                // Emit the steady feed rate segment.
+                float t = l_steady / l;
+                for (int i = 0; i < 4; ++ i) {
+                    line.pos_end[i] = pos_start[i] + (pos_end[i] - pos_start[i]) * t;
+                    line.pos_provided[i] = true;
+                }
+                push_line_to_output(line, pos_start[4], comment);
+                comment = NULL;
+                memcpy(line.pos_start, line.pos_end, sizeof(float)*5);
+                memcpy(pos_start, line.pos_end, sizeof(float)*5);
+            }
+        }
+        // Split the segment into pieces.
+        for (size_t i = 1; i < nSegments; ++ i) {
+            float t = float(i) / float(nSegments);
+            for (size_t j = 0; j < 4; ++ j) {
+                line.pos_end[j] = pos_start[j] + (pos_end[j] - pos_start[j]) * t;
+                line.pos_provided[j] = true;
+            } 
+            // Interpolate the feed rate at the center of the segment.
+            push_line_to_output(line, pos_start[4] + (pos_end[4] - pos_start[4]) * (float(i) - 0.5f) / float(nSegments), comment);
+            comment = NULL;
+            memcpy(line.pos_start, line.pos_end, sizeof(float)*5);
+        }
+		if (l_steady > 0.f && accelerating) {
+            for (int i = 0; i < 4; ++ i) {
+                line.pos_end[i] = pos_end2[i];
+                line.pos_provided[i] = true;
+            }
+            push_line_to_output(line, pos_end[4], comment);
+        }
+    }
+}
+
+void GCodePressureEqualizer::adjust_volumetric_rate()
+{
+    if (circular_buffer_items < 2)
+        return;
+
+    // Go back from the current circular_buffer_pos and lower the feedtrate to decrease the slope of the extrusion rate changes.
+    const size_t idx_head = circular_buffer_idx_head();
+    const size_t idx_tail = circular_buffer_idx_prev(circular_buffer_idx_tail());
+    size_t idx = idx_tail;
+    if (idx == idx_head || ! circular_buffer[idx].extruding())
+        // Nothing to do, the last move is not extruding.
+        return;
+
+    float feedrate_per_extrusion_role[numExtrusionRoles];
+    for (size_t i = 0; i < numExtrusionRoles; ++ i)
+        feedrate_per_extrusion_role[i] = FLT_MAX;
+    feedrate_per_extrusion_role[circular_buffer[idx].extrusion_role] = circular_buffer[idx].volumetric_extrusion_rate_start;
+
+    bool modified = true;
+    while (modified && idx != idx_head) {
+        size_t idx_prev = circular_buffer_idx_prev(idx);
+        for (; ! circular_buffer[idx_prev].extruding() && idx_prev != idx_head; idx_prev = circular_buffer_idx_prev(idx_prev)) ;
+        if (! circular_buffer[idx_prev].extruding())
+        	break;
+        float rate_succ = circular_buffer[idx].volumetric_extrusion_rate_start;
+        // What is the gradient of the extrusion rate between idx_prev and idx?
+        idx = idx_prev;
+        GCodeLine &line = circular_buffer[idx];
+        for (size_t iRole = 1; iRole < numExtrusionRoles; ++ iRole) {
+            float rate_slope = m_max_volumetric_extrusion_rate_slopes[iRole].negative;
+            if (rate_slope == 0)
+                // The negative rate is unlimited.
+                continue;
+            float rate_end = feedrate_per_extrusion_role[iRole];
+            if (iRole == line.extrusion_role && rate_succ < rate_end)
+                rate_end = rate_succ;
+            if (line.volumetric_extrusion_rate_end > rate_end) {
+                line.volumetric_extrusion_rate_end = rate_end;
+                line.modified = true;
+            } else if (iRole == line.extrusion_role) {
+                rate_end = line.volumetric_extrusion_rate_end;
+            } else if (rate_end == FLT_MAX) {
+                // The rate for ExtrusionRole iRole is unlimited.
+                continue;
+            } else {
+                // Use the original, 'floating' extrusion rate as a starting point for the limiter.
+            }
+//            modified = false;
+            float rate_start = rate_end + rate_slope * line.time_corrected();
+            if (rate_start < line.volumetric_extrusion_rate_start) {
+                // Limit the volumetric extrusion rate at the start of this segment due to a segment 
+                // of ExtrusionType iRole, which will be extruded in the future.
+                line.volumetric_extrusion_rate_start = rate_start;
+                line.max_volumetric_extrusion_rate_slope_negative = rate_slope;
+                line.modified = true;
+//              modified = true;
+            }
+            feedrate_per_extrusion_role[iRole] = (iRole == line.extrusion_role) ? line.volumetric_extrusion_rate_start : rate_start;
+        }
+    }
+
+    // Go forward and adjust the feedrate to decrease the slope of the extrusion rate changes.
+    for (size_t i = 0; i < numExtrusionRoles; ++ i)
+        feedrate_per_extrusion_role[i] = FLT_MAX;
+    feedrate_per_extrusion_role[circular_buffer[idx].extrusion_role] = circular_buffer[idx].volumetric_extrusion_rate_end;
+
+    assert(circular_buffer[idx].extruding());
+    while (idx != idx_tail) {
+        size_t idx_next = circular_buffer_idx_next(idx);
+        for (; ! circular_buffer[idx_next].extruding() && idx_next != idx_tail; idx_next = circular_buffer_idx_next(idx_next)) ;
+        if (! circular_buffer[idx_next].extruding())
+        	break;
+        float rate_prec = circular_buffer[idx].volumetric_extrusion_rate_end;
+        // What is the gradient of the extrusion rate between idx_prev and idx?
+        idx = idx_next;
+        GCodeLine &line = circular_buffer[idx];
+        for (size_t iRole = 1; iRole < numExtrusionRoles; ++ iRole) {
+            float rate_slope = m_max_volumetric_extrusion_rate_slopes[iRole].positive;
+            if (rate_slope == 0)
+                // The positive rate is unlimited.
+                continue;
+            float rate_start = feedrate_per_extrusion_role[iRole];
+            if (iRole == line.extrusion_role && rate_prec < rate_start)
+                rate_start = rate_prec;
+            if (line.volumetric_extrusion_rate_start > rate_start) {
+                line.volumetric_extrusion_rate_start = rate_start;
+                line.modified = true;
+            } else if (iRole == line.extrusion_role) {
+                rate_start = line.volumetric_extrusion_rate_start;
+            } else if (rate_start == FLT_MAX) {
+                // The rate for ExtrusionRole iRole is unlimited.
+                continue;
+            } else {
+                // Use the original, 'floating' extrusion rate as a starting point for the limiter.
+            }
+            float rate_end = (rate_slope == 0) ? FLT_MAX : rate_start + rate_slope * line.time_corrected();
+            if (rate_end < line.volumetric_extrusion_rate_end) {
+                // Limit the volumetric extrusion rate at the start of this segment due to a segment 
+                // of ExtrusionType iRole, which was extruded before.
+                line.volumetric_extrusion_rate_end = rate_end;
+                line.max_volumetric_extrusion_rate_slope_positive = rate_slope;
+                line.modified = true;
+            }
+            feedrate_per_extrusion_role[iRole] = (iRole == line.extrusion_role) ? line.volumetric_extrusion_rate_end : rate_end;
+        }
+    }
+}
+
+void GCodePressureEqualizer::push_axis_to_output(const char axis, const float value, bool add_eol)
+{
+    char buf[2048];
+    int len = sprintf(buf, 
+        (axis == 'E') ? " %c%.3f" : " %c%.5f",
+        axis, value);
+    push_to_output(buf, len, add_eol);
+}
+
+void GCodePressureEqualizer::push_to_output(const char *text, const size_t len, bool add_eol)
+{
+    // New length of the output buffer content.
+    size_t len_new = output_buffer_length + len + 1;
+    if (add_eol)
+        ++ len_new;
+
+    // Resize the output buffer to a power of 2 higher than the required memory.
+    if (output_buffer.size() < len_new) {
+        size_t v = len_new;
+        // Compute the next highest power of 2 of 32-bit v
+        // http://graphics.stanford.edu/~seander/bithacks.html
+        v--;
+        v |= v >> 1;
+        v |= v >> 2;
+        v |= v >> 4;
+        v |= v >> 8;
+        v |= v >> 16;
+        v++;
+        output_buffer.resize(v);
+    }
+
+    // Copy the text to the output.
+    if (len != 0) {
+        memcpy(output_buffer.data() + output_buffer_length, text, len);
+        output_buffer_length += len;
+    }
+    if (add_eol)
+        output_buffer[output_buffer_length ++] = '\n';
+    output_buffer[output_buffer_length] = 0;
+}
+
+void GCodePressureEqualizer::push_line_to_output(const GCodeLine &line, const float new_feedrate, const char *comment)
+{
+    push_to_output("G1", 2, false);
+    for (size_t i = 0; i < 3; ++ i)
+        if (line.pos_provided[i])
+            push_axis_to_output('X'+i, line.pos_end[i]);
+    push_axis_to_output('E', m_config->use_relative_e_distances.value ? (line.pos_end[3] - line.pos_start[3]) : line.pos_end[3]);
+//    if (line.pos_provided[4] || fabs(line.feedrate() - new_feedrate) > 1e-5)
+        push_axis_to_output('F', new_feedrate);
+    // output comment and EOL
+    push_to_output(comment, (comment == NULL) ? 0 : strlen(comment), true);
+} 
+
+} // namespace Slic3r

Some files were not shown because too many files changed in this diff