Update to trunk

(bzr r11950.1.148)

11 files changed, 674 insertions, 106 deletions

diff --git a/src/libdepixelize/CMakeLists.txt b/src/libdepixelize/CMakeLists.txt
index e05849e29..895e16e85 100644
--- a/src/libdepixelize/CMakeLists.txt
+++ b/src/libdepixelize/CMakeLists.txt
@@ -7,14 +7,17 @@ set(libdepixelize_SRC
 	kopftracer2011.h
 	splines.h
 
-	priv/branchless.h
+        priv/branchless.h
 	priv/colorspace.h
+	priv/curvature.h
 	priv/homogeneoussplines.h
+	priv/integral
 	priv/iterator.h
+	priv/optimization-kopf2011.h
 	priv/pixelgraph.h
 	priv/point.h
 	priv/simplifiedvoronoi.h
-	priv/splines.h
+	priv/splines-kopf2011.h
 )
 
 add_inkscape_lib(depixelize_LIB "${libdepixelize_SRC}")
diff --git a/src/libdepixelize/Makefile_insert b/src/libdepixelize/Makefile_insert
index 421d32439..8aed7244f 100644
--- a/src/libdepixelize/Makefile_insert
+++ b/src/libdepixelize/Makefile_insert
@@ -5,15 +5,18 @@ libdepixelize/all: libdepixelize/libdepixelize.a
 libdepixelize/clean:
 	rm -f libdepixelize/libdepixelize.a $(libdepixelize_libdepixelize_a_OBJECTS)
 
-libdepixelize_libdepixelize_a_SOURCES =          \
-	libdepixelize/kopftracer2011.cpp         \
-	libdepixelize/kopftracer2011.h           \
-	libdepixelize/splines.h                  \
-	libdepixelize/priv/branchless.h          \
-	libdepixelize/priv/colorspace.h          \
-	libdepixelize/priv/homogeneoussplines.h  \
-	libdepixelize/priv/iterator.h            \
-	libdepixelize/priv/pixelgraph.h          \
-	libdepixelize/priv/point.h               \
-	libdepixelize/priv/simplifiedvoronoi.h   \
-	libdepixelize/priv/splines.h
+libdepixelize_libdepixelize_a_SOURCES =            \
+	libdepixelize/kopftracer2011.cpp           \
+	libdepixelize/kopftracer2011.h             \
+	libdepixelize/splines.h                    \
+	libdepixelize/priv/branchless.h            \
+	libdepixelize/priv/colorspace.h            \
+	libdepixelize/priv/curvature.h             \
+	libdepixelize/priv/homogeneoussplines.h    \
+	libdepixelize/priv/integral.h              \
+	libdepixelize/priv/iterator.h              \
+	libdepixelize/priv/optimization-kopf2011.h \
+	libdepixelize/priv/pixelgraph.h            \
+	libdepixelize/priv/point.h                 \
+	libdepixelize/priv/simplifiedvoronoi.h     \
+	libdepixelize/priv/splines-kopf2011.h
diff --git a/src/libdepixelize/kopftracer2011.cpp b/src/libdepixelize/kopftracer2011.cpp
index 26ad8863b..5e6e26048 100644
--- a/src/libdepixelize/kopftracer2011.cpp
+++ b/src/libdepixelize/kopftracer2011.cpp
@@ -37,7 +37,7 @@
 #include "priv/colorspace.h"
 #include "priv/homogeneoussplines.h"
 #include "priv/branchless.h"
-#include "priv/splines.h"
+#include "priv/splines-kopf2011.h"
 #include "priv/iterator.h"
 
 namespace Tracer {
@@ -146,6 +146,12 @@ SimplifiedVoronoi<T> Kopf2011::_voronoi(const Glib::RefPtr<Gdk::Pixbuf const> &b
     graph.checkConsistency();
 #endif
 
+    _remove_puzzle_pattern(graph);
+
+#ifndef NDEBUG
+    graph.checkConsistency();
+#endif
+
     return SimplifiedVoronoi<T>(graph);
 }
 
@@ -309,6 +315,40 @@ void Kopf2011::_remove_crossing_edges_unsafe(PixelGraph &graph,
     }
 }
 
+inline
+void Kopf2011::_remove_puzzle_pattern(PixelGraph &graph)
+{
+    if ( graph.width() < 2 || graph.height() < 2 )
+        return;
+
+    PixelGraph::iterator it = graph.begin();
+    for ( int i = 0 ; i + 1 != graph.height() ; ++i ) {
+        PixelGraph::iterator it2 = it;
+        for ( int j = 0 ; j + 1 != graph.width() ; ++j ) {
+            // Evil pattern currently not handled correctly in SimplifiedVoronoi
+            if ( it2->adj.right + it2->adj.bottom
+                 + graph.nodeBottomRight(it2)->adj.left
+                 + graph.nodeBottomRight(it2)->adj.top == 3 ) {
+                // We fake a new connection =)
+                it2->adj.right = true;
+                graph.nodeRight(it2)->adj.left = true;
+
+                it2->adj.bottom = true;
+                graph.nodeBottom(it2)->adj.top = true;
+
+                graph.nodeBottomRight(it2)->adj.left = true;
+                graph.nodeBottom(it2)->adj.right = true;
+
+                graph.nodeBottomRight(it2)->adj.top = true;
+                graph.nodeRight(it2)->adj.bottom = true;
+            }
+
+            it2 = graph.nodeRight(it2);
+        }
+        it = graph.nodeBottom(it);
+    }
+}
+
 inline int Heuristics::curves(const PixelGraph &graph,
                               PixelGraph::const_iterator a,
                               PixelGraph::const_iterator b)
@@ -383,22 +423,25 @@ inline void Heuristics::SparsePixels::operator ()(const PixelGraph &graph,
     {
         unsigned x = graph.toX(diagonals[MAIN_DIAGONAL].first.first);
         unsigned y = graph.toY(diagonals[MAIN_DIAGONAL].first.first);
-        unsigned minor = std::min(x, y);
         unsigned displace = radius - 1;
 
-        if ( displace > minor ) {
-            displace = minor;
-            radius = displace + 1;
+        {
+            unsigned minor = std::min(x, y);
+
+            if ( displace > minor ) {
+                displace = minor;
+                radius = displace + 1;
+            }
         }
 
         displace = radius;
 
-        if ( x + displace >= graph.width() ) {
+        if ( x + displace >= unsigned(graph.width()) ) {
             displace = unsigned(graph.width()) - x - 1;
             radius = displace;
         }
 
-        if ( y + displace >= graph.height() ) {
+        if ( y + displace >= unsigned(graph.height()) ) {
             displace = unsigned(graph.height()) - y - 1;
             radius = displace;
         }
@@ -434,7 +477,6 @@ inline void Heuristics::SparsePixels::operator ()(const PixelGraph &graph,
     int minor = std::min(diagonals[0].second, diagonals[1].second);
     for ( int i = 0 ; i != 2 ; ++i )
         diagonals[i].second -= minor;
-
     std::swap(diagonals[0].second, diagonals[1].second);
 }
 
diff --git a/src/libdepixelize/kopftracer2011.h b/src/libdepixelize/kopftracer2011.h
index 73f3b7274..aff39d3d8 100644
--- a/src/libdepixelize/kopftracer2011.h
+++ b/src/libdepixelize/kopftracer2011.h
@@ -73,12 +73,17 @@ public:
     /**
      * # Exceptions
      *
+     * \p options.optimize and options.nthreads will be ignored
+     *
      * Glib::FileError
      * Gdk::PixbufError
      */
     static Splines to_voronoi(const std::string &filename,
                               const Options &options = Options());
 
+    /*
+     * \p options.optimize and options.nthreads will be ignored
+     */
     static Splines to_voronoi(const Glib::RefPtr<Gdk::Pixbuf const> &buf,
                               const Options &options = Options());
 
@@ -95,7 +100,7 @@ public:
                               const Options &options = Options());
 
 private:
-    typedef double Precision;
+    typedef Geom::Coord Precision;
 
     template<class T>
     static SimplifiedVoronoi<T> _voronoi(const Glib::RefPtr<Gdk::Pixbuf const> &buf,
@@ -105,6 +110,7 @@ private:
     static void _remove_crossing_edges_safe(PixelGraph &graph);
     static void _remove_crossing_edges_unsafe(PixelGraph &graph,
                                               const Options &options);
+    static void _remove_puzzle_pattern(PixelGraph &graph);
 };
 
 } // namespace Tracer
diff --git a/src/libdepixelize/priv/curvature.h b/src/libdepixelize/priv/curvature.h
new file mode 100644
index 000000000..3a1af4197
--- /dev/null
+++ b/src/libdepixelize/priv/curvature.h
@@ -0,0 +1,115 @@
+/*  This file is part of the libdepixelize project
+    Copyright (C) 2013 Vinícius dos Santos Oliveira <vini.ipsmaker@gmail.com>
+
+    GNU Lesser General Public License Usage
+    This library is free software; you can redistribute it and/or modify it
+    under the terms of the GNU Lesser General Public License as published by the
+    Free Software Foundation; either version 2.1 of the License, or (at your
+    option) any later version.
+    You should have received a copy of the GNU Lesser General Public License
+    along with this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    GNU General Public License Usage
+    Alternatively, this library may be used under the terms of the GNU General
+    Public License as published by the Free Software Foundation, either version
+    2 of the License, or (at your option) any later version.
+    You should have received a copy of the GNU General Public License along with
+    this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+*/
+
+#ifndef LIBDEPIXELIZE_TRACER_CURVATURE_H
+#define LIBDEPIXELIZE_TRACER_CURVATURE_H
+
+#include "point.h"
+#include <cmath>
+
+namespace Tracer {
+
+/**
+ * Curvature function for a quadratic Bézier curve where the points are know.
+ * Its main use is to numerically integrate the curvature function and then
+ * smooth the B-Splines generated by the Kopf-Lischinski algorithm.
+ */
+template<class T>
+struct Curvature
+{
+    Curvature(Point<T> p0, Point<T> c1, Point<T> p2) :
+        p0(p0), c1(c1), p2(p2)
+    {}
+
+    T operator()(T t) const;
+
+    /**
+     * The derivative of x
+     */
+    T xPrime(T t) const;
+
+    /**
+     * The derivative of y
+     */
+    T yPrime(T t) const;
+
+    /**
+     * The second derivative of x
+     */
+    T xPrimePrime() const;
+
+    /**
+     * The second derivative of y
+     */
+    T yPrimePrime() const;
+
+    Point<T> p0, c1, p2;
+};
+
+template<class T>
+T Curvature<T>::operator()(T t) const
+{
+    T num = xPrime(t) * yPrimePrime() - yPrime(t) * xPrimePrime();
+    T den = std::pow(xPrime(t) * xPrime(t) + yPrime(t) * yPrime(t), T(3) / 2);
+    return num / den;
+}
+
+template<class T>
+T Curvature<T>::xPrime(T t) const
+{
+    return (1-t)*2*(c1.x-p0.x) + t*2*(p2.x-c1.x);
+}
+
+template<class T>
+T Curvature<T>::yPrime(T t) const
+{
+    return (1-t)*2*(c1.y-p0.y) + t*2*(p2.y-c1.y);
+}
+
+template<class T>
+T Curvature<T>::xPrimePrime() const
+{
+    return 2 * (p2.x - 2*c1.x + p0.x);
+}
+
+template<class T>
+T Curvature<T>::yPrimePrime() const
+{
+    return 2 * (p2.y - 2*c1.y + p0.y);
+}
+
+} // namespace Tracer
+
+#endif // LIBDEPIXELIZE_TRACER_CURVATURE_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/libdepixelize/priv/integral.h b/src/libdepixelize/priv/integral.h
new file mode 100644
index 000000000..fc1a49f27
--- /dev/null
+++ b/src/libdepixelize/priv/integral.h
@@ -0,0 +1,61 @@
+/*  This file is part of the libdepixelize project
+    Copyright (C) 2013 Vinícius dos Santos Oliveira <vini.ipsmaker@gmail.com>
+
+    GNU Lesser General Public License Usage
+    This library is free software; you can redistribute it and/or modify it
+    under the terms of the GNU Lesser General Public License as published by the
+    Free Software Foundation; either version 2.1 of the License, or (at your
+    option) any later version.
+    You should have received a copy of the GNU Lesser General Public License
+    along with this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    GNU General Public License Usage
+    Alternatively, this library may be used under the terms of the GNU General
+    Public License as published by the Free Software Foundation, either version
+    2 of the License, or (at your option) any later version.
+    You should have received a copy of the GNU General Public License along with
+    this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+*/
+
+#ifndef LIBDEPIXELIZE_TRACER_INTEGRAL_H
+#define LIBDEPIXELIZE_TRACER_INTEGRAL_H
+
+#include <2geom/coord.h>
+
+namespace Tracer {
+
+/**
+ * Compute the integral numerically using Gaussian Quadrature rule with
+ * \p samples number of samples.
+ */
+template<class T, class F>
+Geom::Coord integral(F f, T begin, T end, unsigned samples)
+{
+    T ret = 0;
+    const T width = (end - begin) / samples;
+
+    for ( unsigned i = 0 ; i != samples ; ++i )
+        ret += width * f(begin + width * (i + .5));
+
+    return ret;
+}
+
+} // namespace Tracer
+
+#endif // LIBDEPIXELIZE_TRACER_INTEGRAL_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/libdepixelize/priv/optimization-kopf2011.h b/src/libdepixelize/priv/optimization-kopf2011.h
new file mode 100644
index 000000000..0c9011ca7
--- /dev/null
+++ b/src/libdepixelize/priv/optimization-kopf2011.h
@@ -0,0 +1,263 @@
+/*  This file is part of the libdepixelize project
+    Copyright (C) 2013 Vinícius dos Santos Oliveira <vini.ipsmaker@gmail.com>
+
+    GNU Lesser General Public License Usage
+    This library is free software; you can redistribute it and/or modify it
+    under the terms of the GNU Lesser General Public License as published by the
+    Free Software Foundation; either version 2.1 of the License, or (at your
+    option) any later version.
+    You should have received a copy of the GNU Lesser General Public License
+    along with this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    GNU General Public License Usage
+    Alternatively, this library may be used under the terms of the GNU General
+    Public License as published by the Free Software Foundation, either version
+    2 of the License, or (at your option) any later version.
+    You should have received a copy of the GNU General Public License along with
+    this library.  If not, see <http://www.gnu.org/licenses/>.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+*/
+
+#ifndef LIBDEPIXELIZE_TRACER_OPTIMIZATION_KOPF2011_H
+#define LIBDEPIXELIZE_TRACER_OPTIMIZATION_KOPF2011_H
+
+#include "curvature.h"
+#include "integral.h"
+#include <cmath>
+#include <limits>
+
+namespace Tracer {
+
+/**
+ * m is angular coeficient
+ */
+template<class T>
+bool is_valid_border_m(T a)
+{
+    if ( a < 0 )
+        a *= -1;
+
+    // TODO: alternative behaviour if has no infinity
+    return a == std::numeric_limits<T>::infinity()
+        || a == 3 || a == 1;
+}
+
+/**
+ * Return true if the four points are considered a border.
+ */
+template<class T>
+bool is_border(const Point<T> (&points)[4])
+{
+    T dy[2];
+    T dx[2];
+    T m[2];
+    if ( points[1].y == points[2].y ) {
+        dy[0] = points[1].y - points[0].y;
+        dy[1] = points[3].y - points[2].y;
+
+        dx[0] = points[1].x - points[0].x;
+        dx[1] = points[3].x - points[2].x;
+
+        m[0] = dy[0] / dx[0];
+        m[1] = dy[1] / dx[1];
+    } else if ( points[1].x == points[2].x ) {
+        // It's easier to have a unified logic, then we'll fake dx and dy
+        dy[0] = points[1].x - points[0].x;
+        dy[1] = points[3].x - points[2].x;
+
+        dx[0] = points[1].y - points[0].y;
+        dx[1] = points[3].y - points[2].y;
+
+        m[0] = dy[0] / dx[0];
+        m[1] = dy[1] / dx[1];
+    } else {
+        return false;
+    }
+
+    return m[0] == -m[1] && is_valid_border_m(m[0]);
+}
+
+template<class T>
+typename std::vector< Point<T> >::iterator
+skip1visible(typename std::vector< Point<T> >::iterator it,
+             typename std::vector< Point<T> >::iterator end)
+{
+    for ( ++it ; it != end ; ++it ) {
+        if ( it->visible )
+            return it;
+    }
+    return end;
+}
+
+/**
+ * Return how many elements should be skipped.
+ */
+template<class T>
+typename std::vector< Point<T> >::difference_type
+border_detection(typename std::vector< Point<T> >::iterator it,
+                 typename std::vector< Point<T> >::iterator end)
+{
+    typename std::vector< Point<T> >::iterator begin = it;
+
+    if ( end - it < 4 )
+        return 0;
+
+    Point<T> last[4];
+    typename std::vector< Point<T> >::iterator prev = it;
+
+    for ( int i = 0 ; i != 4 ; ++i ) {
+        if ( it == end )
+            return 0;
+        last[i] = *it;
+        prev = it;
+        it = skip1visible<T>(it, end);
+    }
+
+    if ( !is_border(last) )
+        return 0;
+
+    if ( it == end )
+        return prev - begin;
+
+    bool got_another = false;
+    for ( it = skip1visible<T>(it, end) ; it != end
+              ; it = skip1visible<T>(it, end) ) {
+        if ( !got_another ) {
+            last[0] = last[2];
+            last[1] = last[3];
+            last[2] = *it;
+
+            got_another = true;
+            continue;
+        }
+        last[3] = *it;
+
+        if ( !is_border(last) )
+            return prev - begin;
+        prev = it;
+    }
+
+    return prev - begin;
+}
+
+template<class T>
+T smoothness_energy(Point<T> c0, Point<T> c1, Point<T> c2)
+{
+    Point<T> p0 = midpoint(c0, c1);
+    Point<T> p2 = midpoint(c1, c2);
+    Curvature<T> cur(p0, c1, p2);
+
+    return std::abs(integral<T>(cur, 0, 1, 16));
+}
+
+template<class T>
+T positional_energy(Point<T> guess, Point<T> initial)
+{
+    using std::pow;
+
+    return pow(pow(guess.x - initial.x, 2)
+               + pow(guess.y - initial.y, 2), 2);
+}
+
+/**
+ * Kopf-Lischinski simple relaxation procedure: a random new offset position
+ * within a small radius around its current location.
+ *
+ * The small radius is not revealed. I chose the empirically determined value of
+ * 0.125. New tests can give a better value for "small". I believe this value
+ * showed up because the optimization sharply penalize larger deviations.
+ */
+template<class T>
+Point<T> optimization_guess(Point<T> p)
+{
+    // See the value explanation in the function documentation.
+    T radius = 0.125;
+
+    T d[] = {
+        (T(std::rand()) / RAND_MAX) * radius * 2  - radius,
+        (T(std::rand()) / RAND_MAX) * radius * 2  - radius
+    };
+
+    return p + Point<T>(d[0], d[1]);
+}
+
+template<class T>
+std::vector< Point<T> > optimize(const std::vector< Point<T> > &path)
+{
+    typedef std::vector< Point<T> > Path;
+
+    Path ret = path;
+
+    /* The number of vertices not constrained by optimization */
+    unsigned n = 0;
+
+    /* Values chosen by test
+     * TODO: make values configurable via function parameters. */
+    const unsigned iterations = 4;
+    const unsigned nguess_per_iteration = 4;
+
+    for ( unsigned i = 0 ; i != iterations ; ++i ) {
+        n = 0;
+
+        /* This iteration bounds is not something to worry about, because the
+         * smallest path has size 4. */
+        for ( typename Path::size_type j = 0 ; j != ret.size() ; ++j ) {
+            Point<T> prev = ( j == 0 ) ? ret.back() : ret[j-1];
+            Point<T> next = ( j + 1 == ret.size() ) ? ret.front() : ret[j+1] ;
+
+            if ( !ret[j].visible || !ret[j].smooth )
+                continue;
+
+            {
+                typename Path::iterator it = ret.begin() + j;
+                typename Path::difference_type skip
+                    = border_detection<T>(it, ret.end());
+                j += skip;
+                if ( j == ret.size() )
+                    break;
+            }
+
+            ++n;
+
+            for ( unsigned k = 0 ; k != nguess_per_iteration ; ++k ) {
+                Point<T> guess = optimization_guess(ret[j]);
+
+                T s = smoothness_energy(prev, guess, next);
+                T p = positional_energy(guess, path[j]);
+
+                T e = s + p;
+
+                T prev_e = smoothness_energy(prev, ret[j], next)
+                    + positional_energy(ret[j], path[j]);
+
+                if ( prev_e > e ) {
+                    // We don't want to screw other metadata, then we manually
+                    // assign the new coords
+                    ret[j].x = guess.x;
+                    ret[j].y = guess.y;
+                }
+            }
+        }
+    }
+
+    return ret;
+}
+
+} // namespace Tracer
+
+#endif // LIBDEPIXELIZE_TRACER_OPTIMIZATION_KOPF2011_H
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/libdepixelize/priv/pixelgraph.h b/src/libdepixelize/priv/pixelgraph.h
index 32523d8ee..9e8c2124a 100644
--- a/src/libdepixelize/priv/pixelgraph.h
+++ b/src/libdepixelize/priv/pixelgraph.h
@@ -378,7 +378,7 @@ inline void PixelGraph::connectAllNeighbors()
 
                 it = nodeRight(it);
             }
-            // After the previous loop, 'it' is poiting to the last node from
+            // After the previous loop, 'it' is pointing to the last node from
             // the row.
             // Go south, then first node in the row (increment 'it' by 1)
             // Go to the second node in the line (increment 'it' by 1)
@@ -445,21 +445,34 @@ inline void PixelGraph::connectAllNeighbors()
     //  ...and the 4 corner nodes
     {
         Node *const top_left = &(*this)[0][0];
-        top_left->adj.right = 1;
-        top_left->adj.bottomright = 1;
-        top_left->adj.bottom = 1;
+
+        if ( _width > 1 )
+            top_left->adj.right = 1;
+
+        if ( _width > 1 && _height > 1 )
+            top_left->adj.bottomright = 1;
+
+        if ( _height > 1 )
+            top_left->adj.bottom = 1;
     }
     if ( _width > 1 ) {
         Node *const top_right = &(*this)[_width - 1][0];
-        top_right->adj.bottom = 1;
-        top_right->adj.bottomleft = 1;
+
+        if ( _height > 1 ) {
+            top_right->adj.bottom = 1;
+            top_right->adj.bottomleft = 1;
+        }
+
         top_right->adj.left = 1;
     }
     if ( _height > 1 ) {
         Node *const down_left = &(*this)[0][_height - 1];
         down_left->adj.top = 1;
-        down_left->adj.topright = 1;
-        down_left->adj.right = 1;
+
+        if ( _width > 1 ) {
+            down_left->adj.topright = 1;
+            down_left->adj.right = 1;
+        }
     }
     if ( _width > 1 && _height > 1 ) {
         Node *const down_right = &(*this)[_width - 1][_height - 1];
diff --git a/src/libdepixelize/priv/point.h b/src/libdepixelize/priv/point.h
index 6bea752ed..dc28a2b84 100644
--- a/src/libdepixelize/priv/point.h
+++ b/src/libdepixelize/priv/point.h
@@ -30,9 +30,9 @@ namespace Tracer {
 template<class T>
 struct Point
 {
-    Point() {}
-    Point(T x, T y) : x(x), y(y) {}
-    Point(T x, T y, bool smooth) : smooth(smooth), x(x), y(y) {}
+    Point() : visible(true) {}
+    Point(T x, T y) : visible(true), x(x), y(y) {}
+    Point(T x, T y, bool smooth) : smooth(smooth), visible(true), x(x), y(y) {}
 
     Point operator+(const Point &rhs) const
     {
@@ -44,21 +44,58 @@ struct Point
         return Point(x / foo, y / foo);
     }
 
+    Point invisible() const
+    {
+        Point p = *this;
+        p.visible = false;
+        return p;
+    }
+
     bool smooth;
 
+    /**
+     * By default, all points are visible, but the poor amount of information
+     * that B-Splines (libdepixelize-specific) allows us to represent forces us
+     * to create additional points. But... these additional points don't need to
+     * be visible.
+     */
+    bool visible;
+
     T x, y;
 };
 
 template<class T>
-inline bool operator==(const Point<T> &lhs, const Point<T> &rhs)
+Point<T> midpoint(const Point<T> &a, const Point<T> &b)
+{
+    return Point<T>((a.x + b.x) / 2, (a.y + b.y) / 2);
+}
+
+template<class T>
+bool operator==(const Point<T> &lhs, const Point<T> &rhs)
 {
     return
+        /*
+         * Will make a better job identifying which points can be eliminated by
+         * cells union.
+         */
 #ifndef LIBDEPIXELIZE_IS_VERY_WELL_TESTED
         lhs.smooth == rhs.smooth &&
 #endif // LIBDEPIXELIZE_IS_VERY_WELL_TESTED
         lhs.x == rhs.x && lhs.y == rhs.y;
 }
 
+template<class T>
+bool weakly_equal(const Point<T> &a, const Point<T> &b)
+{
+    return a.x == b.x && a.y == b.y;
+}
+
+template<class T>
+Geom::Point to_geom_point(Point<T> p)
+{
+    return Geom::Point(p.x, p.y);
+}
+
 } // namespace Tracer
 
 #endif // LIBDEPIXELIZE_TRACER_POINT_H
diff --git a/src/libdepixelize/priv/simplifiedvoronoi.h b/src/libdepixelize/priv/simplifiedvoronoi.h
index a33695ff7..d5ebc36e5 100644
--- a/src/libdepixelize/priv/simplifiedvoronoi.h
+++ b/src/libdepixelize/priv/simplifiedvoronoi.h
@@ -124,11 +124,6 @@ private:
     typedef void (*PointTransform)(Point<T> &p, T dx, T dy);
     typedef bool (*NodeTransform)(PixelGraph::const_iterator);
 
-    static Point<T> _midpoint(const Point<T> &p1, const Point<T> p2)
-    {
-        return Point<T>((p1.x + p2.x) / 2, (p1.y + p2.y) / 2);
-    }
-
     /**
      * Output is translated by -.5 in each axis. This function fixes this error.
      */
@@ -220,7 +215,7 @@ private:
      * code base 4 times smaller and bugs will be MUCH MUCH difficult to hide.
      *
      * Some maintainers may go mad because the level extra level of
-     * abstracation.
+     * abstraction.
      *
      * "All problems in computer science can be solved by another level of
      *  indirection, except for the problem of too many layers of indirection."
@@ -255,6 +250,11 @@ SimplifiedVoronoi<T>::SimplifiedVoronoi(const PixelGraph &graph) :
     if (!graph.size())
         return;
 
+    /*
+     * The insertion order of cells is not a big deal. Here I just follow
+     * the order of PixelGraph arrangement.
+     */
+
     // ...the "center" cells first...
     if ( _width > 2 && _height > 2 ) {
         PixelGraph::const_iterator graph_it = graph.begin() + _width + 1;
@@ -276,7 +276,7 @@ SimplifiedVoronoi<T>::SimplifiedVoronoi(const PixelGraph &graph) :
                 // Bottom-left
                 _complexBottomLeft(graph, graph_it, cells_it, j, i);
             }
-            // After the previous loop, 'it' is poiting to the last cell from
+            // After the previous loop, 'it' is pointing to the last cell from
             // the row.
             // Go south, then first node in the row (increment 'it' by 1)
             // Go to the second node in the line (increment 'it' by 1)
@@ -794,13 +794,13 @@ SimplifiedVoronoi<T>
                              PixelGraph::const_iterator d_it,
                              Cell *const cells_it, int x, int y,
                              PointTransform transform,
-                             NodeTransform top,
+                             NodeTransform,
                              NodeTransform topright,
                              NodeTransform right,
                              NodeTransform bottomright,
                              NodeTransform bottom,
                              NodeTransform bottomleft,
-                             NodeTransform left,
+                             NodeTransform,
                              NodeTransform topleft)
 {
     using colorspace::contour_edge;
@@ -808,31 +808,36 @@ SimplifiedVoronoi<T>
 
     const Point<T> initial(x, y);
 
+    /*
+     * The insertion order of points within the cell is very important. You must
+     * follow current practice: Clockwise order.
+     */
+
     if ( bottomright(a_it) ) {
         // this and bottom-right are connected
 
         bool smooth[2] = {
-            same_color(a_it->rgba, d_it->rgba) + right(a_it),
-            same_color(a_it->rgba, d_it->rgba) + bottom(a_it)
+            same_color(a_it->rgba, d_it->rgba) || right(a_it),
+            same_color(a_it->rgba, d_it->rgba) || bottom(a_it)
         };
 
         Point<T> borderMid = initial;
         {
             transform(borderMid, 1, 1);
-            borderMid = _midpoint(initial, borderMid);
+            borderMid = midpoint(initial, borderMid);
         }
 
         Point<T> vertices[2] = {initial, initial};
         {
             transform(vertices[0], 1, 0);
-            vertices[0] = _adjust(_midpoint(borderMid, vertices[0]), smooth[0]);
+            vertices[0] = _adjust(midpoint(borderMid, vertices[0]), smooth[0]);
 
             transform(vertices[1], 0, 1);
-            vertices[1] = _adjust(_midpoint(borderMid, vertices[1]), smooth[1]);
+            vertices[1] = _adjust(midpoint(borderMid, vertices[1]), smooth[1]);
         }
 
         if ( !smooth[0] ) {
-            cells_it->vertices.push_back(vertices[0]);
+            cells_it->vertices.push_back(vertices[0].invisible());
             {
                 Point<T> another = vertices[0];
                 transform(another,
@@ -841,7 +846,7 @@ SimplifiedVoronoi<T>
                           // y
                           - ( 0.5625
                               - ( topright(a_it) + topleft(b_it) ) * 0.1875 ));
-                cells_it->vertices.push_back(another);
+                cells_it->vertices.push_back(another.invisible());
             }
             {
                 Point<T> another = vertices[0];
@@ -862,7 +867,7 @@ SimplifiedVoronoi<T>
                           // y
                           0.0625
                           + ( bottomright(b_it) - topright(d_it) ) * 0.0625);
-                cells_it->vertices.push_back(another);
+                cells_it->vertices.push_back(another.invisible());
             }
             {
                 transform(vertices[0],
@@ -874,6 +879,7 @@ SimplifiedVoronoi<T>
                               + ( topright(d_it) - topright(a_it)
                                   - topleft(b_it) - bottomright(b_it) )
                               * 0.03125 ));
+                vertices[0].visible = false;
             }
         }
 
@@ -891,7 +897,7 @@ SimplifiedVoronoi<T>
                           0.0625
                           + ( bottomleft(a_it) - bottomleft(d_it)
                               - topleft(c_it) - bottomright(c_it) ) * 0.03125);
-                cells_it->vertices.push_back(another);
+                cells_it->vertices.push_back(another.invisible());
             }
             {
                 Point<T> another = vertices[1];
@@ -901,7 +907,7 @@ SimplifiedVoronoi<T>
                           // y
                           0.1875
                           - ( bottomright(c_it) + bottomleft(d_it) ) * 0.0625);
-                cells_it->vertices.push_back(another);
+                cells_it->vertices.push_back(another.invisible());
             }
             {
                 Point<T> another = vertices[1];
@@ -926,8 +932,9 @@ SimplifiedVoronoi<T>
                           - ( 0.1875
                               - ( bottomleft(a_it) - topleft(c_it) )
                               * 0.1875 ));
-                cells_it->vertices.push_back(another);
+                cells_it->vertices.push_back(another.invisible());
             }
+            vertices[1].visible = false;
         }
 
         cells_it->vertices.push_back(vertices[1]);
@@ -936,7 +943,7 @@ SimplifiedVoronoi<T>
 
         Point<T> vertex = initial;
         transform(vertex, 1, 1);
-        vertex = _adjust(_midpoint(_midpoint(initial, vertex), initial), true);
+        vertex = _adjust(midpoint(midpoint(initial, vertex), initial), true);
         cells_it->vertices.push_back(vertex);
     } else {
         // Connections don't affect the shape of this squared-like
@@ -944,7 +951,7 @@ SimplifiedVoronoi<T>
 
         Point<T> vertex = initial;
         transform(vertex, 1, 1);
-        vertex = _adjust(_midpoint(initial, vertex));
+        vertex = _adjust(midpoint(initial, vertex));
 
         // compute smoothness
         if ( right(a_it) ) {
@@ -969,13 +976,13 @@ SimplifiedVoronoi<T>
                                     - ( ( bottomright(c_it) + topleft(c_it) )
                                         * 0.03125 );
                                 transform(another, - amount, amount);
-                                cells_it->vertices.push_back(another);
+                                cells_it->vertices.push_back(another.invisible());
                             }
                             {
                                 Point<T> another = vertex;
                                 T amount = 0.0625 * bottomright(c_it);
                                 transform(another, amount, 0.25 - amount);
-                                cells_it->vertices.push_back(another);
+                                cells_it->vertices.push_back(another.invisible());
                             }
                             {
                                 Point<T> another = vertex;
@@ -990,8 +997,9 @@ SimplifiedVoronoi<T>
                                 T amount = 0.1875 * topleft(c_it);
                                 transform(another, - ( 0.75 - amount ),
                                           - amount);
-                                cells_it->vertices.push_back(another);
+                                cells_it->vertices.push_back(another.invisible());
                             }
+                            vertex.visible = false;
                         } else if ( twin_is_contour ) {
                             T amount = 0.125
                                 - ( ( bottomleft(d_it) + topright(d_it) )
@@ -1036,13 +1044,13 @@ SimplifiedVoronoi<T>
 
                     if ( !vertex.smooth ) {
                         if ( another_is_contour ) {
-                            cells_it->vertices.push_back(vertex);
+                            cells_it->vertices.push_back(vertex.invisible());
                             {
                                 Point<T> another = vertex;
                                 T amount = 0.1875 * topleft(b_it);
                                 transform(another, - amount,
                                           - ( 0.75 - amount ));
-                                cells_it->vertices.push_back(another);
+                                cells_it->vertices.push_back(another.invisible());
                             }
                             {
                                 Point<T> another = vertex;
@@ -1056,13 +1064,14 @@ SimplifiedVoronoi<T>
                                 Point<T> another = vertex;
                                 T amount = 0.0625 * bottomright(b_it);
                                 transform(another, 0.25 - amount, amount);
-                                cells_it->vertices.push_back(another);
+                                cells_it->vertices.push_back(another.invisible());
                             }
                             {
                                 T amount = 0.125
                                     - (bottomright(b_it) + topleft(b_it))
                                     * 0.03125;
                                 transform(vertex, amount, - amount);
+                                vertex.visible = false;
                             }
                         } else if ( twin_is_contour ) {
                             T amount = 0.125
@@ -1115,13 +1124,13 @@ SimplifiedVoronoi<T>
                                 - ( topleft(c_it) + bottomright(c_it) )
                                 * 0.03125;
                             transform(another, - amount, amount);
-                            cells_it->vertices.push_back(another);
+                            cells_it->vertices.push_back(another.invisible());
                         }
                         {
                             Point<T> another = vertex;
                             T amount = 0.0625 * bottomright(c_it);
                             transform(another, amount, 0.25 - amount);
-                            cells_it->vertices.push_back(another);
+                            cells_it->vertices.push_back(another.invisible());
                         }
                         {
                             Point<T> another = vertex;
@@ -1134,8 +1143,9 @@ SimplifiedVoronoi<T>
                             Point<T> another = vertex;
                             T amount = 0.1875 * topleft(c_it);
                             transform(another, - ( 0.75 - amount ), - amount);
-                            cells_it->vertices.push_back(another);
+                            cells_it->vertices.push_back(another.invisible());
                         }
+                        vertex.visible = false;
                     } else {
                         special = true;
                     }
@@ -1149,7 +1159,7 @@ SimplifiedVoronoi<T>
             }
 
             if ( special ) {
-                cells_it->vertices.push_back(vertex);
+                cells_it->vertices.push_back(vertex.invisible());
                 {
                     Point<T> another = vertex;
                     T amount = 0.1875;
@@ -1159,7 +1169,7 @@ SimplifiedVoronoi<T>
                               - ( 0.75
                                   - ( topleft(b_it) + topright(a_it) )
                                   * amount ));
-                    cells_it->vertices.push_back(another);
+                    cells_it->vertices.push_back(another.invisible());
                 }
                 {
                     Point<T> another = vertex;
@@ -1181,7 +1191,7 @@ SimplifiedVoronoi<T>
                               // y
                               0.25 - amount
                               * ( bottomleft(d_it) + bottomright(c_it) ));
-                    cells_it->vertices.push_back(another);
+                    cells_it->vertices.push_back(another.invisible());
                 }
                 {
                     transform(vertex,
@@ -1192,6 +1202,7 @@ SimplifiedVoronoi<T>
                               ( topleft(b_it) - bottomleft(d_it)
                                 + topright(a_it) - bottomright(c_it) )
                               * 0.03125);
+                    vertex.visible = false;
                 }
             }
         } else if ( right(c_it) ) {
@@ -1212,12 +1223,12 @@ SimplifiedVoronoi<T>
 
                 if ( !vertex.smooth ) {
                     if ( similar_neighbor_is_contour ) {
-                        cells_it->vertices.push_back(vertex);
+                        cells_it->vertices.push_back(vertex.invisible());
                         {
                             Point<T> another = vertex;
                             T amount = 0.1875 * topleft(b_it);
                             transform(another, - amount, - ( 0.75 - amount ));
-                            cells_it->vertices.push_back(another);
+                            cells_it->vertices.push_back(another.invisible());
                         }
                         {
                             Point<T> another = vertex;
@@ -1230,12 +1241,13 @@ SimplifiedVoronoi<T>
                             Point<T> another = vertex;
                             T amount = 0.0625 * bottomright(b_it);
                             transform(another, 0.25 - amount, amount);
-                            cells_it->vertices.push_back(another);
+                            cells_it->vertices.push_back(another.invisible());
                         }
                         {
                             T amount = 0.125
                                 - 0.03125 * (topleft(b_it) + bottomright(b_it));
                             transform(vertex, amount, - amount);
+                            vertex.visible = false;
                         }
                     } else {
                         special = true;
@@ -1261,7 +1273,7 @@ SimplifiedVoronoi<T>
                               - amount
                               * ( topleft(c_it) + topright(d_it)
                                   - bottomleft(a_it) - bottomright(b_it) ));
-                    cells_it->vertices.push_back(another);
+                    cells_it->vertices.push_back(another.invisible());
                 }
                 {
                     Point<T> another = vertex;
@@ -1272,7 +1284,7 @@ SimplifiedVoronoi<T>
                               // y
                               - amount
                               * ( topright(d_it) - bottomright(b_it) ));
-                    cells_it->vertices.push_back(another);
+                    cells_it->vertices.push_back(another.invisible());
                 }
                 {
                     Point<T> another = vertex;
@@ -1295,8 +1307,9 @@ SimplifiedVoronoi<T>
                               // y
                               -  amount
                               * ( topleft(c_it) - bottomleft(a_it) ));
-                    cells_it->vertices.push_back(another);
+                    cells_it->vertices.push_back(another.invisible());
                 }
+                vertex.visible = false;
             }
         } else {
             // there is a 4-color pattern, where the current node
diff --git a/src/libdepixelize/priv/splines.h b/src/libdepixelize/priv/splines-kopf2011.h
index c7ef2b492..eb84c3bfb 100644
--- a/src/libdepixelize/priv/splines.h
+++ b/src/libdepixelize/priv/splines-kopf2011.h
@@ -22,55 +22,65 @@
     Lesser General Public License for more details.
 */
 
-#ifndef LIBDEPIXELIZE_TRACER_SPLINES_PRIV_H
-#define LIBDEPIXELIZE_TRACER_SPLINES_PRIV_H
+#ifndef LIBDEPIXELIZE_TRACER_SPLINES_KOPF2011_H
+#define LIBDEPIXELIZE_TRACER_SPLINES_KOPF2011_H
 
 #include "../splines.h"
 #include "homogeneoussplines.h"
+#include "optimization-kopf2011.h"
 
 namespace Tracer {
 
+/**
+ * Maybe the pass-by-value and then move idiom should be more efficient. But all
+ * this is inlinable and we're not even in C++11 yet.
+ */
 template<class T>
-Geom::Point to_geom_point(Point<T> p)
-{
-    return Geom::Point(p.x, p.y);
-}
-
-template<class T>
-Geom::Path worker_helper(const std::vector< Point<T> > &source, bool optimize)
+Geom::Path worker_helper(const std::vector< Point<T> > &source1, bool optimize)
 {
     typedef Geom::LineSegment Line;
     typedef Geom::QuadraticBezier Quad;
     typedef typename std::vector< Point<T> >::const_iterator iterator;
 
-    iterator it = source.begin();
-    Geom::Path ret(to_geom_point((source.back() + *it) / 2));
+    std::vector< Point<T> > source;
 
-    for ( iterator end = --source.end() ; it != end ; ++it ) {
-        Point<T> next = *(it + 1);
-        Point<T> middle = (*it + next) / 2;
+    if ( optimize )
+        source = Tracer::optimize(source1);
+    else
+        source = source1;
 
-        if ( !it->smooth ) {
-            // TODO: remove redundant colinear points
-            ret.appendNew<Line>(Geom::Point(it->x, it->y));
-            ret.appendNew<Line>(Geom::Point(middle.x, middle.y));
-        } else {
-            ret.appendNew<Quad>(Geom::Point(it->x, it->y),
-                                Geom::Point(middle.x, middle.y));
+    iterator it = source.begin();
+    Point<T> prev = source.back();
+    Geom::Path ret(to_geom_point(midpoint(prev, *it)));
+
+    for ( iterator end = source.end() ; it != end ; ++it ) {
+#if LIBDEPIXELIZE_ENABLE_EXPERIMENTAL_FEATURES
+        // remove redundant points
+        if ( !it->visible ) {
+            prev = *it;
+            continue;
+        }
+#endif // LIBDEPIXELIZE_ENABLE_EXPERIMENTAL_FEATURES
+
+        if ( !prev.visible ) {
+            Geom::Point middle = to_geom_point(midpoint(prev, *it));
+            if ( ret.finalPoint() != middle ) {
+                // All generated invisible points are straight lines
+                ret.appendNew<Line>(middle);
+            }
         }
-    }
 
-    {
-        Point<T> next = source.front();
-        Point<T> middle = (*it + next) / 2;
+        Point<T> next = (it + 1 == end) ? source.front() : *(it + 1);
+        Point<T> middle = midpoint(*it, next);
 
         if ( !it->smooth ) {
-            ret.appendNew<Line>(Geom::Point(it->x, it->y));
-            ret.appendNew<Line>(Geom::Point(middle.x, middle.y));
+            ret.appendNew<Line>(to_geom_point(*it));
+            ret.appendNew<Line>(to_geom_point(middle));
         } else {
-            ret.appendNew<Quad>(Geom::Point(it->x, it->y),
-                                Geom::Point(middle.x, middle.y));
+            ret.appendNew<Quad>(to_geom_point(*it), to_geom_point(middle));
         }
+
+        prev = *it;
     }
 
     return ret;
@@ -98,7 +108,9 @@ void worker(const typename HomogeneousSplines<T>::Polygon &source,
 }
 
 template<typename T>
-Splines::Splines(const SimplifiedVoronoi<T> &diagram)
+Splines::Splines(const SimplifiedVoronoi<T> &diagram) :
+    _width(diagram.width()),
+    _height(diagram.height())
 {
     _paths.reserve(diagram.size());
 
@@ -141,7 +153,7 @@ Splines::Splines(const HomogeneousSplines<T> &homogeneousSplines,
 
 } // namespace Tracer
 
-#endif // LIBDEPIXELIZE_TRACER_SPLINES_PRIV_H
+#endif // LIBDEPIXELIZE_TRACER_SPLINES_KOPF2011_H
 
 /*
   Local Variables: