Adding upstream version 4:7.4.7.upstream/4%7.4.7upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

5 files changed, 2480 insertions, 0 deletions

diff --git a/basegfx/source/workbench/Makefile b/basegfx/source/workbench/Makefile
new file mode 100644
index 000000000..6218141da
--- /dev/null
+++ b/basegfx/source/workbench/Makefile
@@ -0,0 +1,23 @@
+#
+# This file is part of the LibreOffice project.
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#
+# This file incorporates work covered by the following license notice:
+#
+#   Licensed to the Apache Software Foundation (ASF) under one or more
+#   contributor license agreements. See the NOTICE file distributed
+#   with this work for additional information regarding copyright
+#   ownership. The ASF licenses this file to you under the Apache
+#   License, Version 2.0 (the "License"); you may not use this file
+#   except in compliance with the License. You may obtain a copy of
+#   the License at http://www.apache.org/licenses/LICENSE-2.0 .
+#
+
+prog : bezierclip.cxx convexhull.cxx
+	g++ -I. -Wall -g bezierclip.cxx convexhull.cxx -o bezierclip
+
+.cxx.o:
+	g++ -c $(LOCALDEFINES) $(CCFLAGS) $<
diff --git a/basegfx/source/workbench/bezierclip.cxx b/basegfx/source/workbench/bezierclip.cxx
new file mode 100644
index 000000000..676f239ef
--- /dev/null
+++ b/basegfx/source/workbench/bezierclip.cxx
@@ -0,0 +1,2006 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * This file incorporates work covered by the following license notice:
+ *
+ *   Licensed to the Apache Software Foundation (ASF) under one or more
+ *   contributor license agreements. See the NOTICE file distributed
+ *   with this work for additional information regarding copyright
+ *   ownership. The ASF licenses this file to you under the Apache
+ *   License, Version 2.0 (the "License"); you may not use this file
+ *   except in compliance with the License. You may obtain a copy of
+ *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
+ */
+
+#include <iostream>
+#include <cassert>
+#include <algorithm>
+#include <iterator>
+#include <vector>
+#include <utility>
+
+#include <math.h>
+
+#include <bezierclip.hxx>
+#include <gauss.hxx>
+
+// what to test
+#define WITH_ASSERTIONS
+//#define WITH_CONVEXHULL_TEST
+//#define WITH_MULTISUBDIVIDE_TEST
+//#define WITH_FATLINE_TEST
+//#define WITH_CALCFOCUS_TEST
+//#define WITH_SAFEPARAMBASE_TEST
+//#define WITH_SAFEPARAMS_TEST
+//#define WITH_SAFEPARAM_DETAILED_TEST
+//#define WITH_SAFEFOCUSPARAM_CALCFOCUS
+//#define WITH_SAFEFOCUSPARAM_TEST
+//#define WITH_SAFEFOCUSPARAM_DETAILED_TEST
+#define WITH_BEZIERCLIP_TEST
+
+/* Implementation of the so-called 'Fat-Line Bezier Clipping Algorithm' by Sederberg et al.
+ *
+ * Actual reference is: T. W. Sederberg and T Nishita: Curve
+ * intersection using Bezier clipping. In Computer Aided Design, 22
+ * (9), 1990, pp. 538--549
+ */
+
+/* Misc helper
+ * ===========
+ */
+int fallFac( int n, int k )
+{
+#ifdef WITH_ASSERTIONS
+    assert(n>=k); // "For factorials, n must be greater or equal k"
+    assert(n>=0); // "For factorials, n must be positive"
+    assert(k>=0); // "For factorials, k must be positive"
+#endif
+
+    int res( 1 );
+
+    while( k-- && n ) res *= n--;
+
+    return res;
+}
+
+int fac( int n )
+{
+    return fallFac(n, n);
+}
+
+/* Bezier fat line clipping part
+ * =============================
+ */
+
+void Impl_calcFatLine( FatLine& line, const Bezier& c )
+{
+    // Prepare normalized implicit line
+    // ================================
+
+    // calculate vector orthogonal to p1-p4:
+    line.a = -(c.p0.y - c.p3.y);
+    line.b = (c.p0.x - c.p3.x);
+
+    // normalize
+    const double len(std::hypot(line.a, line.b));
+    if( !tolZero(len) )
+    {
+        line.a /= len;
+        line.b /= len;
+    }
+
+    line.c = -(line.a*c.p0.x + line.b*c.p0.y);
+
+    // Determine bounding fat line from it
+    // ===================================
+
+    // calc control point distances
+    const double dP2( calcLineDistance(line.a, line.b, line.c, c.p1.x, c.p1.y ) );
+    const double dP3( calcLineDistance(line.a, line.b, line.c, c.p2.x, c.p2.y ) );
+
+    // calc approximate bounding lines to curve (tight bounds are
+    // possible here, but more expensive to calculate and thus not
+    // worth the overhead)
+    if( dP2 * dP3 > 0.0 )
+    {
+        line.dMin = 3.0/4.0 * std::min(0.0, std::min(dP2, dP3));
+        line.dMax = 3.0/4.0 * std::max(0.0, std::max(dP2, dP3));
+    }
+    else
+    {
+        line.dMin = 4.0/9.0 * std::min(0.0, std::min(dP2, dP3));
+        line.dMax = 4.0/9.0 * std::max(0.0, std::max(dP2, dP3));
+    }
+}
+
+void Impl_calcBounds( Point2D&          leftTop,
+                      Point2D&          rightBottom,
+                      const Bezier&     c1          )
+{
+    leftTop.x = std::min( c1.p0.x, std::min( c1.p1.x, std::min( c1.p2.x, c1.p3.x ) ) );
+    leftTop.y = std::min( c1.p0.y, std::min( c1.p1.y, std::min( c1.p2.y, c1.p3.y ) ) );
+    rightBottom.x = std::max( c1.p0.x, std::max( c1.p1.x, std::max( c1.p2.x, c1.p3.x ) ) );
+    rightBottom.y = std::max( c1.p0.y, std::max( c1.p1.y, std::max( c1.p2.y, c1.p3.y ) ) );
+}
+
+bool Impl_doBBoxIntersect( const Bezier& c1,
+                           const Bezier& c2 )
+{
+    // calc rectangular boxes from c1 and c2
+    Point2D lt1;
+    Point2D rb1;
+    Point2D lt2;
+    Point2D rb2;
+
+    Impl_calcBounds( lt1, rb1, c1 );
+    Impl_calcBounds( lt2, rb2, c2 );
+
+    if( std::min(rb1.x, rb2.x) < std::max(lt1.x, lt2.x) ||
+        std::min(rb1.y, rb2.y) < std::max(lt1.y, lt2.y) )
+    {
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+/* calculates two t's for the given bernstein control polygon: the first is
+ * the intersection of the min value line with the convex hull from
+ * the left, the second is the intersection of the max value line with
+ * the convex hull from the right.
+ */
+bool Impl_calcSafeParams( double&           t1,
+                          double&           t2,
+                          const Polygon2D&  rPoly,
+                          double            lowerYBound,
+                          double            upperYBound )
+{
+    // need the convex hull of the control polygon, as this is
+    // guaranteed to completely bound the curve
+    Polygon2D convHull( convexHull(rPoly) );
+
+    // init min and max buffers
+    t1 = 0.0 ;
+    double currLowerT( 1.0 );
+
+    t2 = 1.0;
+    double currHigherT( 0.0 );
+
+    if( convHull.size() <= 1 )
+        return false; // only one point? Then we're done with clipping
+
+    /* now, clip against lower and higher bounds */
+    Point2D p0;
+    Point2D p1;
+
+    bool bIntersection( false );
+
+    for( Polygon2D::size_type i=0; i<convHull.size(); ++i )
+    {
+        // have to check against convHull.size() segments, as the
+        // convex hull is, by definition, closed. Thus, for the
+        // last point, we take the first point as partner.
+        if( i+1 == convHull.size() )
+        {
+            // close the polygon
+            p0 = convHull[i];
+            p1 = convHull[0];
+        }
+        else
+        {
+            p0 = convHull[i];
+            p1 = convHull[i+1];
+        }
+
+        // is the segment in question within or crossing the
+        // horizontal band spanned by lowerYBound and upperYBound? If
+        // not, we've got no intersection. If yes, we maybe don't have
+        // an intersection, but we've got to update the permissible
+        // range, nevertheless. This is because inside lying segments
+        // leads to full range forbidden.
+        if( (tolLessEqual(p0.y, upperYBound) || tolLessEqual(p1.y, upperYBound)) &&
+            (tolGreaterEqual(p0.y, lowerYBound) || tolGreaterEqual(p1.y, lowerYBound)) )
+        {
+            // calc intersection of convex hull segment with
+            // one of the horizontal bounds lines
+            // to optimize a bit, r_x is calculated only in else case
+            const double r_y( p1.y - p0.y );
+
+            if( tolZero(r_y) )
+            {
+                // r_y is virtually zero, thus we've got a horizontal
+                // line. Now check whether we maybe coincide with lower or
+                // upper horizontal bound line.
+                if( tolEqual(p0.y, lowerYBound) ||
+                    tolEqual(p0.y, upperYBound) )
+                {
+                    // yes, simulate intersection then
+                    currLowerT = std::min(currLowerT, std::min(p0.x, p1.x));
+                    currHigherT = std::max(currHigherT, std::max(p0.x, p1.x));
+                }
+            }
+            else
+            {
+                // check against lower and higher bounds
+                // =====================================
+                const double r_x( p1.x - p0.x );
+
+                // calc intersection with horizontal dMin line
+                const double currTLow( (lowerYBound - p0.y) * r_x / r_y + p0.x );
+
+                // calc intersection with horizontal dMax line
+                const double currTHigh( (upperYBound - p0.y) * r_x / r_y + p0.x );
+
+                currLowerT = std::min(currLowerT, std::min(currTLow, currTHigh));
+                currHigherT = std::max(currHigherT, std::max(currTLow, currTHigh));
+            }
+
+            // set flag that at least one segment is contained or
+            // intersects given horizontal band.
+            bIntersection = true;
+        }
+    }
+
+#ifndef WITH_SAFEPARAMBASE_TEST
+    // limit intersections found to permissible t parameter range
+    t1 = std::max(0.0, currLowerT);
+    t2 = std::min(1.0, currHigherT);
+#endif
+
+    return bIntersection;
+}
+
+/* calculates two t's for the given bernstein polynomial: the first is
+ * the intersection of the min value line with the convex hull from
+ * the left, the second is the intersection of the max value line with
+ * the convex hull from the right.
+ *
+ * The polynomial coefficients c0 to c3 given to this method
+ * must correspond to t values of 0, 1/3, 2/3 and 1, respectively.
+ */
+bool Impl_calcSafeParams_clip( double&          t1,
+                               double&          t2,
+                               const FatLine&   bounds,
+                               double           c0,
+                               double           c1,
+                               double           c2,
+                               double           c3 )
+{
+    /* first of all, determine convex hull of c0-c3 */
+    Polygon2D poly(4);
+    poly[0] = Point2D(0,        c0);
+    poly[1] = Point2D(1.0/3.0,  c1);
+    poly[2] = Point2D(2.0/3.0,  c2);
+    poly[3] = Point2D(1,        c3);
+
+#ifndef WITH_SAFEPARAM_DETAILED_TEST
+
+    return Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax );
+
+#else
+    bool bRet( Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax ) );
+
+    Polygon2D convHull( convexHull( poly ) );
+
+    std::cout << "# convex hull testing" << std::endl
+         << "plot [t=0:1] ";
+    std::cout << " bez("
+         << poly[0].x << ","
+         << poly[1].x << ","
+         << poly[2].x << ","
+         << poly[3].x << ",t),bez("
+         << poly[0].y << ","
+         << poly[1].y << ","
+         << poly[2].y << ","
+         << poly[3].y << ",t), "
+         << "t, " << bounds.dMin << ", "
+         << "t, " << bounds.dMax << ", "
+         << t1 << ", t, "
+         << t2 << ", t, "
+         << "'-' using ($1):($2) title \"control polygon\" with lp, "
+         << "'-' using ($1):($2) title \"convex hull\" with lp" << std::endl;
+
+    unsigned int k;
+    for( k=0; k<poly.size(); ++k )
+    {
+        std::cout << poly[k].x << " " << poly[k].y << std::endl;
+    }
+    std::cout << poly[0].x << " " << poly[0].y << std::endl;
+    std::cout << "e" << std::endl;
+
+    for( k=0; k<convHull.size(); ++k )
+    {
+        std::cout << convHull[k].x << " " << convHull[k].y << std::endl;
+    }
+    std::cout << convHull[0].x << " " << convHull[0].y << std::endl;
+    std::cout << "e" << std::endl;
+
+    return bRet;
+#endif
+}
+
+void Impl_deCasteljauAt( Bezier&        part1,
+                         Bezier&        part2,
+                         const Bezier&  input,
+                         double         t        )
+{
+    // deCasteljau bezier arc, scheme is:
+
+    // First row is    C_0^n,C_1^n,...,C_n^n
+    // Second row is         P_1^n,...,P_n^n
+    // etc.
+    // with P_k^r = (1 - x_s)P_{k-1}^{r-1} + x_s P_k{r-1}
+
+    // this results in:
+
+    // P1  P2  P3  P4
+    // L1  P2  P3  R4
+    //     L2  H   R3
+    //         L3  R2
+    //             L4/R1
+    if( tolZero(t) )
+    {
+        // t is zero -> part2 is input curve, part1 is empty (input.p0, that is)
+        part1.p0.x = part1.p1.x = part1.p2.x = part1.p3.x = input.p0.x;
+        part1.p0.y = part1.p1.y = part1.p2.y = part1.p3.y = input.p0.y;
+        part2 = input;
+    }
+    else if( tolEqual(t, 1.0) )
+    {
+        // t is one -> part1 is input curve, part2 is empty (input.p3, that is)
+        part1 = input;
+        part2.p0.x = part2.p1.x = part2.p2.x = part2.p3.x = input.p3.x;
+        part2.p0.y = part2.p1.y = part2.p2.y = part2.p3.y = input.p3.y;
+    }
+    else
+    {
+        part1.p0.x = input.p0.x;                                    part1.p0.y = input.p0.y;
+        part1.p1.x = (1.0 - t)*part1.p0.x + t*input.p1.x;           part1.p1.y = (1.0 - t)*part1.p0.y + t*input.p1.y;
+        const double Hx ( (1.0 - t)*input.p1.x + t*input.p2.x ),    Hy ( (1.0 - t)*input.p1.y + t*input.p2.y );
+        part1.p2.x = (1.0 - t)*part1.p1.x + t*Hx;                   part1.p2.y = (1.0 - t)*part1.p1.y + t*Hy;
+        part2.p3.x = input.p3.x;                                    part2.p3.y = input.p3.y;
+        part2.p2.x = (1.0 - t)*input.p2.x + t*input.p3.x;           part2.p2.y = (1.0 - t)*input.p2.y + t*input.p3.y;
+        part2.p1.x = (1.0 - t)*Hx + t*part2.p2.x;                   part2.p1.y = (1.0 - t)*Hy + t*part2.p2.y;
+        part2.p0.x = (1.0 - t)*part1.p2.x + t*part2.p1.x;           part2.p0.y = (1.0 - t)*part1.p2.y + t*part2.p1.y;
+        part1.p3.x = part2.p0.x;                                    part1.p3.y = part2.p0.y;
+    }
+}
+
+void printCurvesWithSafeRange( const Bezier& c1, const Bezier& c2, double t1_c1, double t2_c1,
+                               const Bezier& c2_part, const FatLine& bounds_c2 )
+{
+    static int offset = 0;
+
+    std::cout << "# safe param range testing" << std::endl
+         << "plot [t=0.0:1.0] ";
+
+    // clip safe ranges off c1
+    Bezier c1_part1;
+    Bezier c1_part2;
+    Bezier c1_part3;
+
+    // subdivide at t1_c1
+    Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1_c1 );
+    // subdivide at t2_c1
+    Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, t2_c1 );
+
+    // output remaining segment (c1_part1)
+
+    std::cout << "bez("
+         << c1.p0.x+offset << ","
+         << c1.p1.x+offset << ","
+         << c1.p2.x+offset << ","
+         << c1.p3.x+offset << ",t),bez("
+         << c1.p0.y << ","
+         << c1.p1.y << ","
+         << c1.p2.y << ","
+         << c1.p3.y << ",t), bez("
+         << c2.p0.x+offset << ","
+         << c2.p1.x+offset << ","
+         << c2.p2.x+offset << ","
+         << c2.p3.x+offset << ",t),bez("
+         << c2.p0.y << ","
+         << c2.p1.y << ","
+         << c2.p2.y << ","
+         << c2.p3.y << ",t), "
+#if 1
+         << "bez("
+         << c1_part1.p0.x+offset << ","
+         << c1_part1.p1.x+offset << ","
+         << c1_part1.p2.x+offset << ","
+         << c1_part1.p3.x+offset << ",t),bez("
+         << c1_part1.p0.y << ","
+         << c1_part1.p1.y << ","
+         << c1_part1.p2.y << ","
+         << c1_part1.p3.y << ",t), "
+#endif
+#if 1
+         << "bez("
+         << c2_part.p0.x+offset << ","
+         << c2_part.p1.x+offset << ","
+         << c2_part.p2.x+offset << ","
+         << c2_part.p3.x+offset << ",t),bez("
+         << c2_part.p0.y << ","
+         << c2_part.p1.y << ","
+         << c2_part.p2.y << ","
+         << c2_part.p3.y << ",t), "
+#endif
+         << "linex("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c << ",t)+" << offset << ", liney("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c << ",t) title \"fat line (center)\", linex("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c-bounds_c2.dMin << ",t)+" << offset << ", liney("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c-bounds_c2.dMin << ",t) title \"fat line (min) \", linex("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c-bounds_c2.dMax << ",t)+" << offset << ", liney("
+         << bounds_c2.a << ","
+         << bounds_c2.b << ","
+         << bounds_c2.c-bounds_c2.dMax << ",t) title \"fat line (max) \"" << std::endl;
+
+    offset += 1;
+}
+
+void printResultWithFinalCurves( const Bezier& c1, const Bezier& c1_part,
+                                 const Bezier& c2, const Bezier& c2_part,
+                                 double t1_c1, double t2_c1 )
+{
+    static int offset = 0;
+
+    std::cout << "# final result" << std::endl
+         << "plot [t=0.0:1.0] ";
+
+    std::cout << "bez("
+         << c1.p0.x+offset << ","
+         << c1.p1.x+offset << ","
+         << c1.p2.x+offset << ","
+         << c1.p3.x+offset << ",t),bez("
+         << c1.p0.y << ","
+         << c1.p1.y << ","
+         << c1.p2.y << ","
+         << c1.p3.y << ",t), bez("
+         << c1_part.p0.x+offset << ","
+         << c1_part.p1.x+offset << ","
+         << c1_part.p2.x+offset << ","
+         << c1_part.p3.x+offset << ",t),bez("
+         << c1_part.p0.y << ","
+         << c1_part.p1.y << ","
+         << c1_part.p2.y << ","
+         << c1_part.p3.y << ",t), "
+         << " pointmarkx(bez("
+         << c1.p0.x+offset << ","
+         << c1.p1.x+offset << ","
+         << c1.p2.x+offset << ","
+         << c1.p3.x+offset << ","
+         << t1_c1 << "),t), "
+         << " pointmarky(bez("
+         << c1.p0.y << ","
+         << c1.p1.y << ","
+         << c1.p2.y << ","
+         << c1.p3.y << ","
+         << t1_c1 << "),t), "
+         << " pointmarkx(bez("
+         << c1.p0.x+offset << ","
+         << c1.p1.x+offset << ","
+         << c1.p2.x+offset << ","
+         << c1.p3.x+offset << ","
+         << t2_c1 << "),t), "
+         << " pointmarky(bez("
+         << c1.p0.y << ","
+         << c1.p1.y << ","
+         << c1.p2.y << ","
+         << c1.p3.y << ","
+         << t2_c1 << "),t), "
+
+         << "bez("
+         << c2.p0.x+offset << ","
+         << c2.p1.x+offset << ","
+         << c2.p2.x+offset << ","
+         << c2.p3.x+offset << ",t),bez("
+         << c2.p0.y << ","
+         << c2.p1.y << ","
+         << c2.p2.y << ","
+         << c2.p3.y << ",t), "
+         << "bez("
+         << c2_part.p0.x+offset << ","
+         << c2_part.p1.x+offset << ","
+         << c2_part.p2.x+offset << ","
+         << c2_part.p3.x+offset << ",t),bez("
+         << c2_part.p0.y << ","
+         << c2_part.p1.y << ","
+         << c2_part.p2.y << ","
+         << c2_part.p3.y << ",t)" << std::endl;
+
+    offset += 1;
+}
+
+/** determine parameter ranges [0,t1) and (t2,1] on c1, where c1 is guaranteed to lie outside c2.
+      Returns false, if the two curves don't even intersect.
+
+    @param t1
+    Range [0,t1) on c1 is guaranteed to lie outside c2
+
+    @param t2
+    Range (t2,1] on c1 is guaranteed to lie outside c2
+
+    @param c1_orig
+    Original curve c1
+
+    @param c1_part
+    Subdivided current part of c1
+
+    @param c2_orig
+    Original curve c2
+
+    @param c2_part
+    Subdivided current part of c2
+ */
+bool Impl_calcClipRange( double&        t1,
+                         double&        t2,
+                         const Bezier&  c1_orig,
+                         const Bezier&  c1_part,
+                         const Bezier&  c2_orig,
+                         const Bezier&  c2_part )
+{
+    // TODO: Maybe also check fat line orthogonal to P0P3, having P0
+    //       and P3 as the extremal points
+
+    if( Impl_doBBoxIntersect(c1_part, c2_part) )
+    {
+        // Calculate fat lines around c1
+        FatLine bounds_c2;
+
+        // must use the subdivided version of c2, since the fat line
+        // algorithm works implicitly with the convex hull bounding
+        // box.
+        Impl_calcFatLine(bounds_c2, c2_part);
+
+        // determine clip positions on c2. Can use original c1 (which
+        // is necessary anyway, to get the t's on the original curve),
+        // since the distance calculations work directly in the
+        // Bernstein polynomial parameter domain.
+        if( Impl_calcSafeParams_clip( t1, t2, bounds_c2,
+                                      calcLineDistance( bounds_c2.a,
+                                                        bounds_c2.b,
+                                                        bounds_c2.c,
+                                                        c1_orig.p0.x,
+                                                        c1_orig.p0.y    ),
+                                      calcLineDistance( bounds_c2.a,
+                                                        bounds_c2.b,
+                                                        bounds_c2.c,
+                                                        c1_orig.p1.x,
+                                                        c1_orig.p1.y    ),
+                                      calcLineDistance( bounds_c2.a,
+                                                        bounds_c2.b,
+                                                        bounds_c2.c,
+                                                        c1_orig.p2.x,
+                                                        c1_orig.p2.y    ),
+                                      calcLineDistance( bounds_c2.a,
+                                                        bounds_c2.b,
+                                                        bounds_c2.c,
+                                                        c1_orig.p3.x,
+                                                        c1_orig.p3.y    ) ) )
+        {
+            //printCurvesWithSafeRange(c1_orig, c2_orig, t1, t2, c2_part, bounds_c2);
+
+            // they do intersect
+            return true;
+        }
+    }
+
+    // they don't intersect: nothing to do
+    return false;
+}
+
+/* Tangent intersection part
+ * =========================
+ */
+
+void Impl_calcFocus( Bezier& res, const Bezier& c )
+{
+    // arbitrary small value, for now
+    // TODO: find meaningful value
+    const double minPivotValue( 1.0e-20 );
+
+    Point2D::value_type fMatrix[6];
+    Point2D::value_type fRes[2];
+
+    // calc new curve from hodograph, c and linear blend
+
+    // Coefficients for derivative of c are (C_i=n(C_{i+1} - C_i)):
+
+    // 3(P1 - P0), 3(P2 - P1), 3(P3 - P2) (bezier curve of degree 2)
+
+    // The hodograph is then (bezier curve of 2nd degree is P0(1-t)^2 + 2P1(1-t)t + P2t^2):
+
+    // 3(P1 - P0)(1-t)^2 + 6(P2 - P1)(1-t)t + 3(P3 - P2)t^2
+
+    // rotate by 90 degrees: x=-y, y=x and you get the normal vector function N(t):
+
+    // x(t) = -(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
+    // y(t) =   3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2
+
+    // Now, the focus curve is defined to be F(t)=P(t) + c(t)N(t),
+    // where P(t) is the original curve, and c(t)=c0(1-t) + c1 t
+
+    // This results in the following expression for F(t):
+
+    // x(t) =  P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
+    //          (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
+
+    // y(t) =  P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1.t)t^2 + P3.y t^3 +
+    //          (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2)
+
+    // As a heuristic, we set F(0)=F(1) (thus, the curve is closed and _tends_ to be small):
+
+    // For F(0), the following results:
+
+    // x(0) = P0.x - c0 3(P1.y - P0.y)
+    // y(0) = P0.y + c0 3(P1.x - P0.x)
+
+    // For F(1), the following results:
+
+    // x(1) = P3.x - c1 3(P3.y - P2.y)
+    // y(1) = P3.y + c1 3(P3.x - P2.x)
+
+    // Reorder, collect and substitute into F(0)=F(1):
+
+    // P0.x - c0 3(P1.y - P0.y) = P3.x - c1 3(P3.y - P2.y)
+    // P0.y + c0 3(P1.x - P0.x) = P3.y + c1 3(P3.x - P2.x)
+
+    // which yields
+
+    // (P0.y - P1.y)c0 + (P3.y - P2.y)c1 = (P3.x - P0.x)/3
+    // (P1.x - P0.x)c0 + (P2.x - P3.x)c1 = (P3.y - P0.y)/3
+
+    // so, this is what we calculate here (determine c0 and c1):
+    fMatrix[0] = c.p1.x - c.p0.x;
+    fMatrix[1] = c.p2.x - c.p3.x;
+    fMatrix[2] = (c.p3.y - c.p0.y)/3.0;
+    fMatrix[3] = c.p0.y - c.p1.y;
+    fMatrix[4] = c.p3.y - c.p2.y;
+    fMatrix[5] = (c.p3.x - c.p0.x)/3.0;
+
+    // TODO: determine meaningful value for
+    if( !solve(fMatrix, 2, 3, fRes, minPivotValue) )
+    {
+        // TODO: generate meaningful values here
+        // singular or nearly singular system -- use arbitrary
+        // values for res
+        fRes[0] = 0.0;
+        fRes[1] = 1.0;
+
+        std::cerr << "Matrix singular!" << std::endl;
+    }
+
+    // now, the reordered and per-coefficient collected focus curve is
+    // the following third degree bezier curve F(t):
+
+    // x(t) =  P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
+    //          (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
+    //      =  P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
+    //         (3c0P1.y(1-t)^3 - 3c0P0.y(1-t)^3 + 6c0P2.y(1-t)^2t - 6c0P1.y(1-t)^2t +
+    //          3c0P3.y(1-t)t^2 - 3c0P2.y(1-t)t^2 +
+    //          3c1P1.y(1-t)^2t - 3c1P0.y(1-t)^2t + 6c1P2.y(1-t)t^2 - 6c1P1.y(1-t)t^2 +
+    //          3c1P3.yt^3 - 3c1P2.yt^3)
+    //      =  (P0.x - 3 c0 P1.y + 3 c0 P0.y)(1-t)^3 +
+    //         3(P1.x - c1 P1.y + c1 P0.y - 2 c0 P2.y + 2 c0 P1.y)(1-t)^2t +
+    //         3(P2.x - 2 c1 P2.y + 2 c1 P1.y - c0 P3.y + c0 P2.y)(1-t)t^2 +
+    //         (P3.x - 3 c1 P3.y + 3 c1 P2.y)t^3
+    //      =  (P0.x - 3 c0(P1.y - P0.y))(1-t)^3 +
+    //         3(P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))(1-t)^2t +
+    //         3(P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))(1-t)t^2 +
+    //         (P3.x - 3 c1(P3.y - P2.y))t^3
+
+    // y(t) =  P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
+    //          (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2)
+    //      =  P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
+    //         3c0(P1.x - P0.x)(1-t)^3 + 6c0(P2.x - P1.x)(1-t)^2t + 3c0(P3.x - P2.x)(1-t)t^2 +
+    //         3c1(P1.x - P0.x)(1-t)^2t + 6c1(P2.x - P1.x)(1-t)t^2 + 3c1(P3.x - P2.x)t^3
+    //      =  (P0.y + 3 c0 (P1.x - P0.x))(1-t)^3 +
+    //         3(P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))(1-t)^2t +
+    //         3(P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))(1-t)t^2 +
+    //         (P3.y + 3 c1 (P3.x - P2.x))t^3
+
+    // Therefore, the coefficients F0 to F3 of the focus curve are:
+
+    // F0.x = (P0.x - 3 c0(P1.y - P0.y))                    F0.y = (P0.y + 3 c0 (P1.x - P0.x))
+    // F1.x = (P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))   F1.y = (P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))
+    // F2.x = (P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))  F2.y = (P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))
+    // F3.x = (P3.x - 3 c1(P3.y - P2.y))                    F3.y = (P3.y + 3 c1 (P3.x - P2.x))
+
+    res.p0.x = c.p0.x - 3*fRes[0]*(c.p1.y - c.p0.y);
+    res.p1.x = c.p1.x - fRes[1]*(c.p1.y - c.p0.y) - 2*fRes[0]*(c.p2.y - c.p1.y);
+    res.p2.x = c.p2.x - 2*fRes[1]*(c.p2.y - c.p1.y) - fRes[0]*(c.p3.y - c.p2.y);
+    res.p3.x = c.p3.x - 3*fRes[1]*(c.p3.y - c.p2.y);
+
+    res.p0.y = c.p0.y + 3*fRes[0]*(c.p1.x - c.p0.x);
+    res.p1.y = c.p1.y + 2*fRes[0]*(c.p2.x - c.p1.x) + fRes[1]*(c.p1.x - c.p0.x);
+    res.p2.y = c.p2.y + fRes[0]*(c.p3.x - c.p2.x) + 2*fRes[1]*(c.p2.x - c.p1.x);
+    res.p3.y = c.p3.y + 3*fRes[1]*(c.p3.x - c.p2.x);
+}
+
+bool Impl_calcSafeParams_focus( double&         t1,
+                                double&         t2,
+                                const Bezier&   curve,
+                                const Bezier&   focus )
+{
+    // now, we want to determine which normals of the original curve
+    // P(t) intersect with the focus curve F(t). The condition for
+    // this statement is P'(t)(P(t) - F) = 0, i.e. hodograph P'(t) and
+    // line through P(t) and F are perpendicular.
+    // If you expand this equation, you end up with something like
+
+    // (\sum_{i=0}^n (P_i - F)B_i^n(t))^T (\sum_{j=0}^{n-1} n(P_{j+1} - P_j)B_j^{n-1}(t))
+
+    // Multiplying that out (as the scalar product is linear, we can
+    // extract some terms) yields:
+
+    // (P_i - F)^T n(P_{j+1} - P_j) B_i^n(t)B_j^{n-1}(t) + ...
+
+    // If we combine the B_i^n(t)B_j^{n-1}(t) product, we arrive at a
+    // Bernstein polynomial of degree 2n-1, as
+
+    // \binom{n}{i}(1-t)^{n-i}t^i) \binom{n-1}{j}(1-t)^{n-1-j}t^j) =
+    // \binom{n}{i}\binom{n-1}{j}(1-t)^{2n-1-i-j}t^{i+j}
+
+    // Thus, with the defining equation for a 2n-1 degree Bernstein
+    // polynomial
+
+    // \sum_{i=0}^{2n-1} d_i B_i^{2n-1}(t)
+
+    // the d_i are calculated as follows:
+
+    // d_i = \sum_{j+k=i, j\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{j}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_j - F)
+
+    // Okay, but F is now not a single point, but itself a curve
+    // F(u). Thus, for every value of u, we get a different 2n-1
+    // bezier curve from the above equation. Therefore, we have a
+    // tensor product bezier patch, with the following defining
+    // equation:
+
+    // d(t,u) = \sum_{i=0}^{2n-1} \sum_{j=0}^m B_i^{2n-1}(t) B_j^{m}(u) d_{ij}, where
+    // d_{ij} = \sum_{k+l=i, l\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{l}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_l - F_j)
+
+    // as above, only that now F is one of the focus' control points.
+
+    // Note the difference in the binomial coefficients to the
+    // reference paper, these formulas most probably contained a typo.
+
+    // To determine, where D(t,u) is _not_ zero (these are the parts
+    // of the curve that don't share normals with the focus and can
+    // thus be safely clipped away), we project D(u,t) onto the
+    // (d(t,u), t) plane, determine the convex hull there and proceed
+    // as for the curve intersection part (projection is orthogonal to
+    // u axis, thus simply throw away u coordinate).
+
+    // \fallfac are so-called falling factorials (see Concrete
+    // Mathematics, p. 47 for a definition).
+
+    // now, for tensor product bezier curves, the convex hull property
+    // holds, too. Thus, we simply project the control points (t_{ij},
+    // u_{ij}, d_{ij}) onto the (t,d) plane and calculate the
+    // intersections of the convex hull with the t axis, as for the
+    // bezier clipping case.
+
+    // calc polygon of control points (t_{ij}, d_{ij}):
+
+    const int n( 3 ); // cubic bezier curves, as a matter of fact
+    const int i_card( 2*n );
+    const int j_card( n + 1 );
+    const int k_max( n-1 );
+    Polygon2D controlPolygon( i_card*j_card ); // vector of (t_{ij}, d_{ij}) in row-major order
+
+    int i, j, k, l; // variable notation from formulas above and Sederberg article
+    Point2D::value_type d;
+    for( i=0; i<i_card; ++i )
+    {
+        for( j=0; j<j_card; ++j )
+        {
+            // calc single d_{ij} sum:
+            for( d=0.0, k=std::max(0,i-n); k<=k_max && k<=i; ++k )
+            {
+                l = i - k; // invariant: k + l = i
+                assert(k>=0 && k<=n-1); // k \in {0,...,n-1}
+                assert(l>=0 && l<=n);   // l \in {0,...,n}
+
+                // TODO: find, document and assert proper limits for n and int's max_val.
+                // This becomes important should anybody wants to use
+                // this code for higher-than-cubic beziers
+                d += static_cast<double>(fallFac(n,l)*fallFac(n-1,k)*fac(i)) /
+                    static_cast<double>(fac(l)*fac(k) * fallFac(2*n-1,i)) * n *
+                    ( (curve[k+1].x - curve[k].x)*(curve[l].x - focus[j].x) +   // dot product here
+                      (curve[k+1].y - curve[k].y)*(curve[l].y - focus[j].y) );
+            }
+
+            // Note that the t_{ij} values are evenly spaced on the
+            // [0,1] interval, thus t_{ij}=i/(2n-1)
+            controlPolygon[ i*j_card + j ] = Point2D( i/(2.0*n-1.0), d );
+        }
+    }
+
+#ifndef WITH_SAFEFOCUSPARAM_DETAILED_TEST
+
+    // calc safe parameter range, to determine [0,t1] and [t2,1] where
+    // no zero crossing is guaranteed.
+    return Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 );
+
+#else
+    bool bRet( Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ) );
+
+    Polygon2D convHull( convexHull( controlPolygon ) );
+
+    std::cout << "# convex hull testing (focus)" << std::endl
+         << "plot [t=0:1] ";
+    std::cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
+         << "'-' using ($1):($2) title \"convex hull\" with lp" << std::endl;
+
+    unsigned int count;
+    for( count=0; count<controlPolygon.size(); ++count )
+    {
+        std::cout << controlPolygon[count].x << " " << controlPolygon[count].y << std::endl;
+    }
+    std::cout << controlPolygon[0].x << " " << controlPolygon[0].y << std::endl;
+    std::cout << "e" << std::endl;
+
+    for( count=0; count<convHull.size(); ++count )
+    {
+        std::cout << convHull[count].x << " " << convHull[count].y << std::endl;
+    }
+    std::cout << convHull[0].x << " " << convHull[0].y << std::endl;
+    std::cout << "e" << std::endl;
+
+    return bRet;
+#endif
+}
+
+/** Calc all values t_i on c1, for which safeRanges functor does not
+    give a safe range on c1 and c2.
+
+    This method is the workhorse of the bezier clipping. Because c1
+    and c2 must be alternatingly tested against each other (first
+    determine safe parameter interval on c1 with regard to c2, then
+    the other way around), we call this method recursively with c1 and
+    c2 swapped.
+
+    @param result
+    Output iterator where the final t values are added to. If curves
+    don't intersect, nothing is added.
+
+    @param delta
+    Maximal allowed distance to true critical point (measured in the
+    original curve's coordinate system)
+
+    @param safeRangeFunctor
+    Functor object, that must provide the following operator():
+    bool safeRangeFunctor( double& t1,
+                           double& t2,
+                           const Bezier& c1_orig,
+                           const Bezier& c1_part,
+                           const Bezier& c2_orig,
+                           const Bezier& c2_part );
+    This functor must calculate the safe ranges [0,t1] and [t2,1] on
+    c1_orig, where c1_orig is 'safe' from c2_part. If the whole
+    c1_orig is safe, false must be returned, true otherwise.
+ */
+template <class Functor> void Impl_applySafeRanges_rec( std::back_insert_iterator< std::vector< std::pair<double, double> > >&    result,
+                                                        double                                                                          delta,
+                                                        const Functor&                                                                  safeRangeFunctor,
+                                                        int                                                                             recursionLevel,
+                                                        const Bezier&                                                                   c1_orig,
+                                                        const Bezier&                                                                   c1_part,
+                                                        double                                                                          last_t1_c1,
+                                                        double                                                                          last_t2_c1,
+                                                        const Bezier&                                                                   c2_orig,
+                                                        const Bezier&                                                                   c2_part,
+                                                        double                                                                          last_t1_c2,
+                                                        double                                                                          last_t2_c2  )
+{
+    // check end condition
+    // ===================
+
+    // TODO: tidy up recursion handling. maybe put everything in a
+    // struct and swap that here at method entry
+
+    // TODO: Implement limit on recursion depth. Should that limit be
+    // reached, chances are that we're on a higher-order tangency. For
+    // this case, AW proposed to take the middle of the current
+    // interval, and to correct both curve's tangents at that new
+    // endpoint to be equal. That virtually generates a first-order
+    // tangency, and justifies to return a single intersection
+    // point. Otherwise, inside/outside test might fail here.
+
+    for( int i=0; i<recursionLevel; ++i ) std::cerr << " ";
+    if( recursionLevel % 2 )
+    {
+        std::cerr << std::endl << "level: " << recursionLevel
+             << " t: "
+             << last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0
+             << ", c1: " << last_t1_c2 << " " << last_t2_c2
+             << ", c2: " << last_t1_c1 << " " << last_t2_c1
+             << std::endl;
+    }
+    else
+    {
+        std::cerr << std::endl << "level: " << recursionLevel
+             << " t: "
+             << last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0
+             << ", c1: " << last_t1_c1 << " " << last_t2_c1
+             << ", c2: " << last_t1_c2 << " " << last_t2_c2
+             << std::endl;
+    }
+
+    // refine solution
+    // ===============
+
+    double t1_c1, t2_c1;
+
+    // Note: we first perform the clipping and only test for precision
+    // sufficiency afterwards, since we want to exploit the fact that
+    // Impl_calcClipRange returns false if the curves don't
+    // intersect. We would have to check that separately for the end
+    // condition, otherwise.
+
+    // determine safe range on c1_orig
+    if( safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ) )
+    {
+        // now, t1 and t2 are calculated on the original curve
+        // (but against a fat line calculated from the subdivided
+        // c2, namely c2_part). If the [t1,t2] range is outside
+        // our current [last_t1,last_t2] range, we're done in this
+        // branch - the curves no longer intersect.
+        if( tolLessEqual(t1_c1, last_t2_c1) && tolGreaterEqual(t2_c1, last_t1_c1) )
+        {
+            // As noted above, t1 and t2 are calculated on the
+            // original curve, but against a fat line
+            // calculated from the subdivided c2, namely
+            // c2_part. Our domain to work on is
+            // [last_t1,last_t2], on the other hand, so values
+            // of [t1,t2] outside that range are irrelevant
+            // here. Clip range appropriately.
+            t1_c1 = std::max(t1_c1, last_t1_c1);
+            t2_c1 = std::min(t2_c1, last_t2_c1);
+
+            // TODO: respect delta
+            // for now, end condition is just a fixed threshold on the t's
+
+            // check end condition
+            // ===================
+
+#if 1
+            if( fabs(last_t2_c1 - last_t1_c1) < 0.0001 &&
+                fabs(last_t2_c2 - last_t1_c2) < 0.0001  )
+#else
+            if( fabs(last_t2_c1 - last_t1_c1) < 0.01 &&
+                fabs(last_t2_c2 - last_t1_c2) < 0.01    )
+#endif
+            {
+                // done. Add to result
+                if( recursionLevel % 2 )
+                {
+                    // uneven level: have to swap the t's, since curves are swapped, too
+                    *result++ = std::make_pair( last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0,
+                                                  last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0 );
+                }
+                else
+                {
+                    *result++ = std::make_pair( last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0,
+                                                  last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0 );
+                }
+
+#if 0
+                //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, last_t1_c1, last_t2_c1 );
+                printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, t1_c1, t2_c1 );
+#else
+                // calc focus curve of c2
+                Bezier focus;
+                Impl_calcFocus(focus, c2_part); // need to use subdivided c2
+
+                safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
+
+                //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, t1_c1, t2_c1 );
+                printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, last_t1_c1, last_t2_c1 );
+#endif
+            }
+            else
+            {
+                // heuristic: if parameter range is not reduced by at least
+                // 20%, subdivide longest curve, and clip shortest against
+                // both parts of longest
+//                if( (last_t2_c1 - last_t1_c1 - t2_c1 + t1_c1) / (last_t2_c1 - last_t1_c1) < 0.2 )
+                if( false )
+                {
+                    // subdivide and descend
+                    // =====================
+
+                    Bezier part1;
+                    Bezier part2;
+
+                    double intervalMiddle;
+
+                    if( last_t2_c1 - last_t1_c1 > last_t2_c2 - last_t1_c2 )
+                    {
+                        // subdivide c1
+                        // ============
+
+                        intervalMiddle = last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0;
+
+                        // subdivide at the middle of the interval (as
+                        // we're not subdividing on the original
+                        // curve, this simply amounts to subdivision
+                        // at 0.5)
+                        Impl_deCasteljauAt( part1, part2, c1_part, 0.5 );
+
+                        // and descend recursively with swapped curves
+                        Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
+                                                  c2_orig, c2_part, last_t1_c2, last_t2_c2,
+                                                  c1_orig, part1, last_t1_c1, intervalMiddle );
+
+                        Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
+                                                  c2_orig, c2_part, last_t1_c2, last_t2_c2,
+                                                  c1_orig, part2, intervalMiddle, last_t2_c1 );
+                    }
+                    else
+                    {
+                        // subdivide c2
+                        // ============
+
+                        intervalMiddle = last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0;
+
+                        // subdivide at the middle of the interval (as
+                        // we're not subdividing on the original
+                        // curve, this simply amounts to subdivision
+                        // at 0.5)
+                        Impl_deCasteljauAt( part1, part2, c2_part, 0.5 );
+
+                        // and descend recursively with swapped curves
+                        Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
+                                                  c2_orig, part1, last_t1_c2, intervalMiddle,
+                                                  c1_orig, c1_part, last_t1_c1, last_t2_c1 );
+
+                        Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
+                                                  c2_orig, part2, intervalMiddle, last_t2_c2,
+                                                  c1_orig, c1_part, last_t1_c1, last_t2_c1 );
+                    }
+                }
+                else
+                {
+                    // apply calculated clip
+                    // =====================
+
+                    // clip safe ranges off c1_orig
+                    Bezier c1_part1;
+                    Bezier c1_part2;
+                    Bezier c1_part3;
+
+                    // subdivide at t1_c1
+                    Impl_deCasteljauAt( c1_part1, c1_part2, c1_orig, t1_c1 );
+
+                    // subdivide at t2_c1. As we're working on
+                    // c1_part2 now, we have to adapt t2_c1 since
+                    // we're no longer in the original parameter
+                    // interval. This is based on the following
+                    // assumption: t2_new = (t2-t1)/(1-t1), which
+                    // relates the t2 value into the new parameter
+                    // range [0,1] of c1_part2.
+                    Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2_c1-t1_c1)/(1.0-t1_c1) );
+
+                    // descend with swapped curves and c1_part1 as the
+                    // remaining (middle) segment
+                    Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
+                                              c2_orig, c2_part, last_t1_c2, last_t2_c2,
+                                              c1_orig, c1_part1, t1_c1, t2_c1 );
+                }
+            }
+        }
+    }
+}
+
+struct ClipBezierFunctor
+{
+    bool operator()( double& t1_c1,
+                     double& t2_c1,
+                     const Bezier& c1_orig,
+                     const Bezier& c1_part,
+                     const Bezier& c2_orig,
+                     const Bezier& c2_part ) const
+    {
+        return Impl_calcClipRange( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
+    }
+};
+
+struct BezierTangencyFunctor
+{
+    bool operator()( double& t1_c1,
+                     double& t2_c1,
+                     const Bezier& c1_orig,
+                     const Bezier& c1_part,
+                     const Bezier& c2_orig,
+                     const Bezier& c2_part ) const
+    {
+        // calc focus curve of c2
+        Bezier focus;
+        Impl_calcFocus(focus, c2_part); // need to use subdivided c2
+                                        // here, as the whole curve is
+                                        // used for focus calculation
+
+        // determine safe range on c1_orig
+        bool bRet( Impl_calcSafeParams_focus( t1_c1, t2_c1,
+                                              c1_orig, // use orig curve here, need t's on original curve
+                                              focus ) );
+
+        std::cerr << "range: " << t2_c1 - t1_c1 << ", ret: " << bRet << std::endl;
+
+        return bRet;
+    }
+};
+
+/** Perform a bezier clip (curve against curve)
+
+    @param result
+    Output iterator where the final t values are added to. This
+    iterator will remain empty, if there are no intersections.
+
+    @param delta
+    Maximal allowed distance to true intersection (measured in the
+    original curve's coordinate system)
+ */
+void clipBezier( std::back_insert_iterator< std::vector< std::pair<double, double> > >&   result,
+                 double                                                                         delta,
+                 const Bezier&                                                                  c1,
+                 const Bezier&                                                                  c2        )
+{
+#if 0
+    // first of all, determine list of collinear normals. Collinear
+    // normals typically separate two intersections, thus, subdivide
+    // at all collinear normal's t values beforehand. This will cater
+    // for tangent intersections, where two or more intersections are
+    // infinitesimally close together.
+
+    // TODO: evaluate effects of higher-than-second-order
+    // tangencies. Sederberg et al. state that collinear normal
+    // algorithm then degrades quickly.
+
+    std::vector< std::pair<double,double> > results;
+    std::back_insert_iterator< std::vector< std::pair<double, double> > > ii(results);
+
+    Impl_calcCollinearNormals( ii, delta, 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
+
+    // As Sederberg's collinear normal theorem is only sufficient, not
+    // necessary for two intersections left and right, we've to test
+    // all segments generated by the collinear normal algorithm
+    // against each other. In other words, if the two curves are both
+    // divided in a left and a right part, the collinear normal
+    // theorem does _not_ state that the left part of curve 1 does not
+    // e.g. intersect with the right part of curve 2.
+
+    // divide c1 and c2 at collinear normal intersection points
+    std::vector< Bezier > c1_segments( results.size()+1 );
+    std::vector< Bezier > c2_segments( results.size()+1 );
+    Bezier c1_remainder( c1 );
+    Bezier c2_remainder( c2 );
+    unsigned int i;
+    for( i=0; i<results.size(); ++i )
+    {
+        Bezier c1_part2;
+        Impl_deCasteljauAt( c1_segments[i], c1_part2, c1_remainder, results[i].first );
+        c1_remainder = c1_part2;
+
+        Bezier c2_part2;
+        Impl_deCasteljauAt( c2_segments[i], c2_part2, c2_remainder, results[i].second );
+        c2_remainder = c2_part2;
+    }
+    c1_segments[i] = c1_remainder;
+    c2_segments[i] = c2_remainder;
+
+    // now, c1/c2_segments contain all segments, then
+    // clip every resulting segment against every other
+    unsigned int c1_curr, c2_curr;
+    for( c1_curr=0; c1_curr<c1_segments.size(); ++c1_curr )
+    {
+        for( c2_curr=0; c2_curr<c2_segments.size(); ++c2_curr )
+        {
+            if( c1_curr != c2_curr )
+            {
+                Impl_clipBezier_rec(result, delta, 0,
+                                    c1_segments[c1_curr], c1_segments[c1_curr],
+                                    0.0, 1.0,
+                                    c2_segments[c2_curr], c2_segments[c2_curr],
+                                    0.0, 1.0);
+            }
+        }
+    }
+#else
+    Impl_applySafeRanges_rec( result, delta, BezierTangencyFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
+    //Impl_applySafeRanges_rec( result, delta, ClipBezierFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
+#endif
+    // that's it, boys'n'girls!
+}
+
+int main(int argc, const char *argv[])
+{
+    double curr_Offset( 0 );
+    unsigned int i,j,k;
+
+    Bezier someCurves[] =
+        {
+//            {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.0)},
+//            {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
+//            {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)}
+//            {Point2D(0.25+1,0.5),Point2D(0.25+1,0.708333),Point2D(0.423611+1,0.916667),Point2D(0.770833+1,0.980324)},
+//            {Point2D(0.0+1,0.0),Point2D(0.0+1,1.0),Point2D(1.0+1,1.0),Point2D(1.0+1,0.5)}
+
+// tangency1
+//            {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)},
+//            {Point2D(0.0,1.0),Point2D(0.1,1.0),Point2D(0.4,1.0),Point2D(0.5,1.0)}
+
+//            {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
+//            {Point2D(0.60114,0.933091),Point2D(0.69461,0.969419),Point2D(0.80676,0.992976),Point2D(0.93756,0.998663)}
+//            {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
+//            {Point2D(0.62712,0.828427),Point2D(0.76305,0.828507),Point2D(0.88555,0.77312),Point2D(0.95069,0.67325)}
+
+// clipping1
+//            {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
+//            {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)}
+
+// tangency2
+//            {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)},
+//            {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
+
+// tangency3
+//            {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
+//            {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)}
+
+// tangency4
+//            {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
+//            {Point2D(0.26,0.4),Point2D(0.25,0.5),Point2D(0.25,0.5),Point2D(0.26,0.6)}
+
+// tangency5
+//            {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
+//            {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
+
+// tangency6
+//            {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
+//            {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
+
+// tangency7
+//            {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)},
+//            {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
+
+// tangency Sederberg example
+            {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)},
+            {Point2D(5.33+9.311,0.0),Point2D(5.33+9.311-13.279,4.205),Point2D(5.33+9.311-10.681,9.119),Point2D(5.33+9.311-2.603,10.254)}
+
+// clipping2
+//            {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
+//            {Point2D(0.2575,0.4),Point2D(0.2475,0.5),Point2D(0.2475,0.5),Point2D(0.2575,0.6)}
+
+//            {Point2D(0.0,0.1),Point2D(0.2,3.5),Point2D(1.0,-2.5),Point2D(1.1,1.2)},
+//            {Point2D(0.0,1.0),Point2D(3.5,0.9),Point2D(-2.5,0.1),Point2D(1.1,0.2)}
+//            {Point2D(0.0,0.1),Point2D(0.2,3.0),Point2D(1.0,-2.0),Point2D(1.1,1.2)},
+//            {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)}
+//            {Point2D(0.0,1.0),Point2D(3.0,0.9),Point2D(-2.0,0.1),Point2D(1.1,0.2)}
+//            {Point2D(0.0,4.0),Point2D(0.1,5.0),Point2D(0.9,5.0),Point2D(1.0,4.0)},
+//            {Point2D(0.0,0.0),Point2D(0.1,0.5),Point2D(0.9,0.5),Point2D(1.0,0.0)},
+//            {Point2D(0.0,0.1),Point2D(0.1,1.5),Point2D(0.9,1.5),Point2D(1.0,0.1)},
+//            {Point2D(0.0,-4.0),Point2D(0.1,-5.0),Point2D(0.9,-5.0),Point2D(1.0,-4.0)}
+        };
+
+    // output gnuplot setup
+    std::cout << "#!/usr/bin/gnuplot -persist" << std::endl
+         << "#" << std::endl
+         << "# automatically generated by bezierclip, don't change!" << std::endl
+         << "#" << std::endl
+         << "set parametric" << std::endl
+         << "bez(p,q,r,s,t) = p*(1-t)**3+q*3*(1-t)**2*t+r*3*(1-t)*t**2+s*t**3" << std::endl
+         << "bezd(p,q,r,s,t) = 3*(q-p)*(1-t)**2+6*(r-q)*(1-t)*t+3*(s-r)*t**2" << std::endl
+         << "pointmarkx(c,t) = c-0.03*t" << std::endl
+         << "pointmarky(c,t) = c+0.03*t" << std::endl
+         << "linex(a,b,c,t) = a*-c + t*-b" << std::endl
+         << "liney(a,b,c,t) = b*-c + t*a" << std::endl << std::endl
+         << "# end of setup" << std::endl << std::endl;
+
+#ifdef WITH_CONVEXHULL_TEST
+    // test convex hull algorithm
+    const double convHull_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# convex hull testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Polygon2D aTestPoly(4);
+        aTestPoly[0] = someCurves[i].p0;
+        aTestPoly[1] = someCurves[i].p1;
+        aTestPoly[2] = someCurves[i].p2;
+        aTestPoly[3] = someCurves[i].p3;
+
+        aTestPoly[0].x += convHull_xOffset;
+        aTestPoly[1].x += convHull_xOffset;
+        aTestPoly[2].x += convHull_xOffset;
+        aTestPoly[3].x += convHull_xOffset;
+
+        std::cout << " bez("
+             << aTestPoly[0].x << ","
+             << aTestPoly[1].x << ","
+             << aTestPoly[2].x << ","
+             << aTestPoly[3].x << ",t),bez("
+             << aTestPoly[0].y << ","
+             << aTestPoly[1].y << ","
+             << aTestPoly[2].y << ","
+             << aTestPoly[3].y << ",t), '-' using ($1):($2) title \"convex hull " << i << "\" with lp";
+
+        if( i+1<sizeof(someCurves)/sizeof(Bezier) )
+            std::cout << ",\\" << std::endl;
+        else
+            std::cout << std::endl;
+    }
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Polygon2D aTestPoly(4);
+        aTestPoly[0] = someCurves[i].p0;
+        aTestPoly[1] = someCurves[i].p1;
+        aTestPoly[2] = someCurves[i].p2;
+        aTestPoly[3] = someCurves[i].p3;
+
+        aTestPoly[0].x += convHull_xOffset;
+        aTestPoly[1].x += convHull_xOffset;
+        aTestPoly[2].x += convHull_xOffset;
+        aTestPoly[3].x += convHull_xOffset;
+
+        Polygon2D convHull( convexHull(aTestPoly) );
+
+        for( k=0; k<convHull.size(); ++k )
+        {
+            std::cout << convHull[k].x << " " << convHull[k].y << std::endl;
+        }
+        std::cout << convHull[0].x << " " << convHull[0].y << std::endl;
+        std::cout << "e" << std::endl;
+    }
+#endif
+
+#ifdef WITH_MULTISUBDIVIDE_TEST
+    // test convex hull algorithm
+    const double multiSubdivide_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# multi subdivide testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Bezier c( someCurves[i] );
+        Bezier c1_part1;
+        Bezier c1_part2;
+        Bezier c1_part3;
+
+        c.p0.x += multiSubdivide_xOffset;
+        c.p1.x += multiSubdivide_xOffset;
+        c.p2.x += multiSubdivide_xOffset;
+        c.p3.x += multiSubdivide_xOffset;
+
+        const double t1( 0.1+i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
+        const double t2( 0.9-i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
+
+        // subdivide at t1
+        Impl_deCasteljauAt( c1_part1, c1_part2, c, t1 );
+
+        // subdivide at t2_c1. As we're working on
+        // c1_part2 now, we have to adapt t2_c1 since
+        // we're no longer in the original parameter
+        // interval. This is based on the following
+        // assumption: t2_new = (t2-t1)/(1-t1), which
+        // relates the t2 value into the new parameter
+        // range [0,1] of c1_part2.
+        Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
+
+        // subdivide at t2
+        Impl_deCasteljauAt( c1_part3, c1_part2, c, t2 );
+
+        std::cout << " bez("
+             << c1_part1.p0.x << ","
+             << c1_part1.p1.x << ","
+             << c1_part1.p2.x << ","
+             << c1_part1.p3.x << ",t), bez("
+             << c1_part1.p0.y+0.01 << ","
+             << c1_part1.p1.y+0.01 << ","
+             << c1_part1.p2.y+0.01 << ","
+             << c1_part1.p3.y+0.01 << ",t) title \"middle " << i << "\", "
+             << " bez("
+             << c1_part2.p0.x << ","
+             << c1_part2.p1.x << ","
+             << c1_part2.p2.x << ","
+             << c1_part2.p3.x << ",t), bez("
+             << c1_part2.p0.y << ","
+             << c1_part2.p1.y << ","
+             << c1_part2.p2.y << ","
+             << c1_part2.p3.y << ",t) title \"right " << i << "\", "
+             << " bez("
+             << c1_part3.p0.x << ","
+             << c1_part3.p1.x << ","
+             << c1_part3.p2.x << ","
+             << c1_part3.p3.x << ",t), bez("
+             << c1_part3.p0.y << ","
+             << c1_part3.p1.y << ","
+             << c1_part3.p2.y << ","
+             << c1_part3.p3.y << ",t) title \"left " << i << "\"";
+
+        if( i+1<sizeof(someCurves)/sizeof(Bezier) )
+            std::cout << ",\\" << std::endl;
+        else
+            std::cout << std::endl;
+    }
+#endif
+
+#ifdef WITH_FATLINE_TEST
+    // test fatline algorithm
+    const double fatLine_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# fat line testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Bezier c( someCurves[i] );
+
+        c.p0.x += fatLine_xOffset;
+        c.p1.x += fatLine_xOffset;
+        c.p2.x += fatLine_xOffset;
+        c.p3.x += fatLine_xOffset;
+
+        FatLine line;
+
+        Impl_calcFatLine(line, c);
+
+        std::cout << " bez("
+             << c.p0.x << ","
+             << c.p1.x << ","
+             << c.p2.x << ","
+             << c.p3.x << ",t), bez("
+             << c.p0.y << ","
+             << c.p1.y << ","
+             << c.p2.y << ","
+             << c.p3.y << ",t) title \"bezier " << i << "\", linex("
+             << line.a << ","
+             << line.b << ","
+             << line.c << ",t), liney("
+             << line.a << ","
+             << line.b << ","
+             << line.c << ",t) title \"fat line (center) on " << i << "\", linex("
+             << line.a << ","
+             << line.b << ","
+             << line.c-line.dMin << ",t), liney("
+             << line.a << ","
+             << line.b << ","
+             << line.c-line.dMin << ",t) title \"fat line (min) on " << i << "\", linex("
+             << line.a << ","
+             << line.b << ","
+             << line.c-line.dMax << ",t), liney("
+             << line.a << ","
+             << line.b << ","
+             << line.c-line.dMax << ",t) title \"fat line (max) on " << i << "\"";
+
+        if( i+1<sizeof(someCurves)/sizeof(Bezier) )
+            std::cout << ",\\" << std::endl;
+        else
+            std::cout << std::endl;
+    }
+#endif
+
+#ifdef WITH_CALCFOCUS_TEST
+    // test focus curve algorithm
+    const double focus_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# focus line testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Bezier c( someCurves[i] );
+
+        c.p0.x += focus_xOffset;
+        c.p1.x += focus_xOffset;
+        c.p2.x += focus_xOffset;
+        c.p3.x += focus_xOffset;
+
+        // calc focus curve
+        Bezier focus;
+        Impl_calcFocus(focus, c);
+
+        std::cout << " bez("
+             << c.p0.x << ","
+             << c.p1.x << ","
+             << c.p2.x << ","
+             << c.p3.x << ",t), bez("
+             << c.p0.y << ","
+             << c.p1.y << ","
+             << c.p2.y << ","
+             << c.p3.y << ",t) title \"bezier " << i << "\", bez("
+             << focus.p0.x << ","
+             << focus.p1.x << ","
+             << focus.p2.x << ","
+             << focus.p3.x << ",t), bez("
+             << focus.p0.y << ","
+             << focus.p1.y << ","
+             << focus.p2.y << ","
+             << focus.p3.y << ",t) title \"focus " << i << "\"";
+
+        if( i+1<sizeof(someCurves)/sizeof(Bezier) )
+            std::cout << ",\\" << std::endl;
+        else
+            std::cout << std::endl;
+    }
+#endif
+
+#ifdef WITH_SAFEPARAMBASE_TEST
+    // test safe params base method
+    double safeParamsBase_xOffset( curr_Offset );
+    std::cout << "# safe param base method testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Bezier c( someCurves[i] );
+
+        c.p0.x += safeParamsBase_xOffset;
+        c.p1.x += safeParamsBase_xOffset;
+        c.p2.x += safeParamsBase_xOffset;
+        c.p3.x += safeParamsBase_xOffset;
+
+        Polygon2D poly(4);
+        poly[0] = c.p0;
+        poly[1] = c.p1;
+        poly[2] = c.p2;
+        poly[3] = c.p3;
+
+        double t1, t2;
+
+        bool bRet( Impl_calcSafeParams( t1, t2, poly, 0, 1 ) );
+
+        Polygon2D convHull( convexHull( poly ) );
+
+        std::cout << " bez("
+             << poly[0].x << ","
+             << poly[1].x << ","
+             << poly[2].x << ","
+             << poly[3].x << ",t),bez("
+             << poly[0].y << ","
+             << poly[1].y << ","
+             << poly[2].y << ","
+             << poly[3].y << ",t), "
+             << "t+" << safeParamsBase_xOffset << ", 0, "
+             << "t+" << safeParamsBase_xOffset << ", 1, ";
+        if( bRet )
+        {
+            std::cout << t1+safeParamsBase_xOffset << ", t, "
+                 << t2+safeParamsBase_xOffset << ", t, ";
+        }
+        std::cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
+             << "'-' using ($1):($2) title \"convex hull\" with lp";
+
+        if( i+1<sizeof(someCurves)/sizeof(Bezier) )
+            std::cout << ",\\" << std::endl;
+        else
+            std::cout << std::endl;
+
+        safeParamsBase_xOffset += 2;
+    }
+
+    safeParamsBase_xOffset = curr_Offset;
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        Bezier c( someCurves[i] );
+
+        c.p0.x += safeParamsBase_xOffset;
+        c.p1.x += safeParamsBase_xOffset;
+        c.p2.x += safeParamsBase_xOffset;
+        c.p3.x += safeParamsBase_xOffset;
+
+        Polygon2D poly(4);
+        poly[0] = c.p0;
+        poly[1] = c.p1;
+        poly[2] = c.p2;
+        poly[3] = c.p3;
+
+        double t1, t2;
+
+        Impl_calcSafeParams( t1, t2, poly, 0, 1 );
+
+        Polygon2D convHull( convexHull( poly ) );
+
+        unsigned int k;
+        for( k=0; k<poly.size(); ++k )
+        {
+            std::cout << poly[k].x << " " << poly[k].y << std::endl;
+        }
+        std::cout << poly[0].x << " " << poly[0].y << std::endl;
+        std::cout << "e" << std::endl;
+
+        for( k=0; k<convHull.size(); ++k )
+        {
+            std::cout << convHull[k].x << " " << convHull[k].y << std::endl;
+        }
+        std::cout << convHull[0].x << " " << convHull[0].y << std::endl;
+        std::cout << "e" << std::endl;
+
+        safeParamsBase_xOffset += 2;
+    }
+    curr_Offset += 20;
+#endif
+
+#ifdef WITH_SAFEPARAMS_TEST
+    // test safe parameter range algorithm
+    const double safeParams_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# safe param range testing" << std::endl
+         << "plot [t=0.0:1.0] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
+        {
+            Bezier c1( someCurves[i] );
+            Bezier c2( someCurves[j] );
+
+            c1.p0.x += safeParams_xOffset;
+            c1.p1.x += safeParams_xOffset;
+            c1.p2.x += safeParams_xOffset;
+            c1.p3.x += safeParams_xOffset;
+            c2.p0.x += safeParams_xOffset;
+            c2.p1.x += safeParams_xOffset;
+            c2.p2.x += safeParams_xOffset;
+            c2.p3.x += safeParams_xOffset;
+
+            double t1, t2;
+
+            if( Impl_calcClipRange(t1, t2, c1, c1, c2, c2) )
+            {
+                // clip safe ranges off c1
+                Bezier c1_part1;
+                Bezier c1_part2;
+                Bezier c1_part3;
+
+                // subdivide at t1_c1
+                Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
+                // subdivide at t2_c1
+                Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
+
+                // output remaining segment (c1_part1)
+
+                std::cout << " bez("
+                     << c1.p0.x << ","
+                     << c1.p1.x << ","
+                     << c1.p2.x << ","
+                     << c1.p3.x << ",t),bez("
+                     << c1.p0.y << ","
+                     << c1.p1.y << ","
+                     << c1.p2.y << ","
+                     << c1.p3.y << ",t), bez("
+                     << c2.p0.x << ","
+                     << c2.p1.x << ","
+                     << c2.p2.x << ","
+                     << c2.p3.x << ",t),bez("
+                     << c2.p0.y << ","
+                     << c2.p1.y << ","
+                     << c2.p2.y << ","
+                     << c2.p3.y << ",t), bez("
+                     << c1_part1.p0.x << ","
+                     << c1_part1.p1.x << ","
+                     << c1_part1.p2.x << ","
+                     << c1_part1.p3.x << ",t),bez("
+                     << c1_part1.p0.y << ","
+                     << c1_part1.p1.y << ","
+                     << c1_part1.p2.y << ","
+                     << c1_part1.p3.y << ",t)";
+
+                if( i+2<sizeof(someCurves)/sizeof(Bezier) )
+                    std::cout << ",\\" << std::endl;
+                else
+                    std::cout << std::endl;
+            }
+        }
+    }
+#endif
+
+#ifdef WITH_SAFEPARAM_DETAILED_TEST
+    // test safe parameter range algorithm
+    const double safeParams2_xOffset( curr_Offset );
+    curr_Offset += 20;
+    if( sizeof(someCurves)/sizeof(Bezier) > 1 )
+    {
+        Bezier c1( someCurves[0] );
+        Bezier c2( someCurves[1] );
+
+        c1.p0.x += safeParams2_xOffset;
+        c1.p1.x += safeParams2_xOffset;
+        c1.p2.x += safeParams2_xOffset;
+        c1.p3.x += safeParams2_xOffset;
+        c2.p0.x += safeParams2_xOffset;
+        c2.p1.x += safeParams2_xOffset;
+        c2.p2.x += safeParams2_xOffset;
+        c2.p3.x += safeParams2_xOffset;
+
+        double t1, t2;
+
+        // output happens here
+        Impl_calcClipRange(t1, t2, c1, c1, c2, c2);
+    }
+#endif
+
+#ifdef WITH_SAFEFOCUSPARAM_TEST
+    // test safe parameter range from focus algorithm
+    const double safeParamsFocus_xOffset( curr_Offset );
+    curr_Offset += 20;
+    std::cout << "# safe param range from focus testing" << std::endl
+         << "plot [t=0.0:1.0] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
+        {
+            Bezier c1( someCurves[i] );
+            Bezier c2( someCurves[j] );
+
+            c1.p0.x += safeParamsFocus_xOffset;
+            c1.p1.x += safeParamsFocus_xOffset;
+            c1.p2.x += safeParamsFocus_xOffset;
+            c1.p3.x += safeParamsFocus_xOffset;
+            c2.p0.x += safeParamsFocus_xOffset;
+            c2.p1.x += safeParamsFocus_xOffset;
+            c2.p2.x += safeParamsFocus_xOffset;
+            c2.p3.x += safeParamsFocus_xOffset;
+
+            double t1, t2;
+
+            Bezier focus;
+#ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
+#if 0
+            {
+                // clip safe ranges off c1_orig
+                Bezier c1_part1;
+                Bezier c1_part2;
+                Bezier c1_part3;
+
+                // subdivide at t1_c1
+                Impl_deCasteljauAt( c1_part1, c1_part2, c2, 0.30204 );
+
+                // subdivide at t2_c1. As we're working on
+                // c1_part2 now, we have to adapt t2_c1 since
+                // we're no longer in the original parameter
+                // interval. This is based on the following
+                // assumption: t2_new = (t2-t1)/(1-t1), which
+                // relates the t2 value into the new parameter
+                // range [0,1] of c1_part2.
+                Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (0.57151-0.30204)/(1.0-0.30204) );
+
+                c2 = c1_part1;
+                Impl_calcFocus( focus, c2 );
+            }
+#else
+            Impl_calcFocus( focus, c2 );
+#endif
+#else
+            focus = c2;
+#endif
+            // determine safe range on c1
+            bool bRet( Impl_calcSafeParams_focus( t1, t2,
+                                                  c1, focus ) );
+
+            std::cerr << "t1: " << t1 << ", t2: " << t2 << std::endl;
+
+            // clip safe ranges off c1
+            Bezier c1_part1;
+            Bezier c1_part2;
+            Bezier c1_part3;
+
+            // subdivide at t1_c1
+            Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
+            // subdivide at t2_c1
+            Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
+
+            // output remaining segment (c1_part1)
+
+            std::cout << " bez("
+                 << c1.p0.x << ","
+                 << c1.p1.x << ","
+                 << c1.p2.x << ","
+                 << c1.p3.x << ",t),bez("
+                 << c1.p0.y << ","
+                 << c1.p1.y << ","
+                 << c1.p2.y << ","
+                 << c1.p3.y << ",t) title \"c1\", "
+#ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
+                 << "bez("
+                 << c2.p0.x << ","
+                 << c2.p1.x << ","
+                 << c2.p2.x << ","
+                 << c2.p3.x << ",t),bez("
+                 << c2.p0.y << ","
+                 << c2.p1.y << ","
+                 << c2.p2.y << ","
+                 << c2.p3.y << ",t) title \"c2\", "
+                 << "bez("
+                 << focus.p0.x << ","
+                 << focus.p1.x << ","
+                 << focus.p2.x << ","
+                 << focus.p3.x << ",t),bez("
+                 << focus.p0.y << ","
+                 << focus.p1.y << ","
+                 << focus.p2.y << ","
+                 << focus.p3.y << ",t) title \"focus\"";
+#else
+                 << "bez("
+                 << c2.p0.x << ","
+                 << c2.p1.x << ","
+                 << c2.p2.x << ","
+                 << c2.p3.x << ",t),bez("
+                 << c2.p0.y << ","
+                 << c2.p1.y << ","
+                 << c2.p2.y << ","
+                 << c2.p3.y << ",t) title \"focus\"";
+#endif
+            if( bRet )
+            {
+                std::cout << ", bez("
+                     << c1_part1.p0.x << ","
+                     << c1_part1.p1.x << ","
+                     << c1_part1.p2.x << ","
+                     << c1_part1.p3.x << ",t),bez("
+                     << c1_part1.p0.y+0.01 << ","
+                     << c1_part1.p1.y+0.01 << ","
+                     << c1_part1.p2.y+0.01 << ","
+                     << c1_part1.p3.y+0.01 << ",t) title \"part\"";
+            }
+
+            if( i+2<sizeof(someCurves)/sizeof(Bezier) )
+                std::cout << ",\\" << std::endl;
+            else
+                std::cout << std::endl;
+        }
+    }
+#endif
+
+#ifdef WITH_SAFEFOCUSPARAM_DETAILED_TEST
+    // test safe parameter range algorithm
+    const double safeParams3_xOffset( curr_Offset );
+    curr_Offset += 20;
+    if( sizeof(someCurves)/sizeof(Bezier) > 1 )
+    {
+        Bezier c1( someCurves[0] );
+        Bezier c2( someCurves[1] );
+
+        c1.p0.x += safeParams3_xOffset;
+        c1.p1.x += safeParams3_xOffset;
+        c1.p2.x += safeParams3_xOffset;
+        c1.p3.x += safeParams3_xOffset;
+        c2.p0.x += safeParams3_xOffset;
+        c2.p1.x += safeParams3_xOffset;
+        c2.p2.x += safeParams3_xOffset;
+        c2.p3.x += safeParams3_xOffset;
+
+        double t1, t2;
+
+        Bezier focus;
+#ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
+        Impl_calcFocus( focus, c2 );
+#else
+        focus = c2;
+#endif
+
+        // determine safe range on c1, output happens here
+        Impl_calcSafeParams_focus( t1, t2,
+                                   c1, focus );
+    }
+#endif
+
+#ifdef WITH_BEZIERCLIP_TEST
+    std::vector< std::pair<double, double> >                                result;
+    std::back_insert_iterator< std::vector< std::pair<double, double> > > ii(result);
+
+    // test full bezier clipping
+    const double bezierClip_xOffset( curr_Offset );
+    std::cout << std::endl << std::endl << "# bezier clip testing" << std::endl
+         << "plot [t=0:1] ";
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
+        {
+            Bezier c1( someCurves[i] );
+            Bezier c2( someCurves[j] );
+
+            c1.p0.x += bezierClip_xOffset;
+            c1.p1.x += bezierClip_xOffset;
+            c1.p2.x += bezierClip_xOffset;
+            c1.p3.x += bezierClip_xOffset;
+            c2.p0.x += bezierClip_xOffset;
+            c2.p1.x += bezierClip_xOffset;
+            c2.p2.x += bezierClip_xOffset;
+            c2.p3.x += bezierClip_xOffset;
+
+            std::cout << " bez("
+                 << c1.p0.x << ","
+                 << c1.p1.x << ","
+                 << c1.p2.x << ","
+                 << c1.p3.x << ",t),bez("
+                 << c1.p0.y << ","
+                 << c1.p1.y << ","
+                 << c1.p2.y << ","
+                 << c1.p3.y << ",t), bez("
+                 << c2.p0.x << ","
+                 << c2.p1.x << ","
+                 << c2.p2.x << ","
+                 << c2.p3.x << ",t),bez("
+                 << c2.p0.y << ","
+                 << c2.p1.y << ","
+                 << c2.p2.y << ","
+                 << c2.p3.y << ",t), '-' using (bez("
+                 << c1.p0.x << ","
+                 << c1.p1.x << ","
+                 << c1.p2.x << ","
+                 << c1.p3.x
+                 << ",$1)):(bez("
+                 << c1.p0.y << ","
+                 << c1.p1.y << ","
+                 << c1.p2.y << ","
+                 << c1.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << i << ")\", "
+                 << " '-' using (bez("
+                 << c2.p0.x << ","
+                 << c2.p1.x << ","
+                 << c2.p2.x << ","
+                 << c2.p3.x
+                 << ",$1)):(bez("
+                 << c2.p0.y << ","
+                 << c2.p1.y << ","
+                 << c2.p2.y << ","
+                 << c2.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << j << ")\"";
+
+            if( i+2<sizeof(someCurves)/sizeof(Bezier) )
+                std::cout << ",\\" << std::endl;
+            else
+                std::cout << std::endl;
+        }
+    }
+    for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
+    {
+        for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
+        {
+            result.clear();
+            Bezier c1( someCurves[i] );
+            Bezier c2( someCurves[j] );
+
+            c1.p0.x += bezierClip_xOffset;
+            c1.p1.x += bezierClip_xOffset;
+            c1.p2.x += bezierClip_xOffset;
+            c1.p3.x += bezierClip_xOffset;
+            c2.p0.x += bezierClip_xOffset;
+            c2.p1.x += bezierClip_xOffset;
+            c2.p2.x += bezierClip_xOffset;
+            c2.p3.x += bezierClip_xOffset;
+
+            clipBezier( ii, 0.00001, c1, c2 );
+
+            for( k=0; k<result.size(); ++k )
+            {
+                std::cout << result[k].first << std::endl;
+            }
+            std::cout << "e" << std::endl;
+
+            for( k=0; k<result.size(); ++k )
+            {
+                std::cout << result[k].second << std::endl;
+            }
+            std::cout << "e" << std::endl;
+        }
+    }
+#endif
+
+    return 0;
+}
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/basegfx/source/workbench/bezierclip.hxx b/basegfx/source/workbench/bezierclip.hxx
new file mode 100644
index 000000000..54cceb414
--- /dev/null
+++ b/basegfx/source/workbench/bezierclip.hxx
@@ -0,0 +1,84 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * This file incorporates work covered by the following license notice:
+ *
+ *   Licensed to the Apache Software Foundation (ASF) under one or more
+ *   contributor license agreements. See the NOTICE file distributed
+ *   with this work for additional information regarding copyright
+ *   ownership. The ASF licenses this file to you under the Apache
+ *   License, Version 2.0 (the "License"); you may not use this file
+ *   except in compliance with the License. You may obtain a copy of
+ *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
+ */
+
+#pragma once
+
+#include <vector>
+
+struct Point2D
+{
+    typedef double value_type;
+    Point2D( double _x, double _y ) : x(_x), y(_y) {}
+    Point2D() : x(), y() {}
+    double x;
+    double y;
+};
+
+struct Bezier
+{
+    Point2D p0;
+    Point2D p1;
+    Point2D p2;
+    Point2D p3;
+
+    Point2D&        operator[]( int i ) { return reinterpret_cast<Point2D*>(this)[i]; }
+    const Point2D&  operator[]( int i ) const { return reinterpret_cast<const Point2D*>(this)[i]; }
+};
+
+struct FatLine
+{
+    // line L through p1 and p4 in normalized implicit form
+    double a;
+    double b;
+    double c;
+
+    // the upper and lower distance from this line
+    double dMin;
+    double dMax;
+};
+
+template <typename DataType> DataType calcLineDistance( const DataType& a,
+                                                        const DataType& b,
+                                                        const DataType& c,
+                                                        const DataType& x,
+                                                        const DataType& y )
+{
+    return a*x + b*y + c;
+}
+
+typedef std::vector< Point2D > Polygon2D;
+
+/* little abs template */
+template <typename NumType> NumType absval( NumType x )
+{
+    return x<0 ? -x : x;
+}
+
+Polygon2D convexHull( const Polygon2D& rPoly );
+
+// TODO: find proper epsilon here (try std::numeric_limits<NumType>::epsilon()?)!
+constexpr auto DBL_EPSILON = 1.0e-100;
+
+/* little approximate comparisons */
+template <typename NumType> bool tolZero( NumType n ) { return fabs(n) < DBL_EPSILON; }
+template <typename NumType> bool tolEqual( NumType n1, NumType n2 ) { return tolZero(n1-n2); }
+template <typename NumType> bool tolLessEqual( NumType n1, NumType n2 ) { return tolEqual(n1,n2) || n1<n2; }
+template <typename NumType> bool tolGreaterEqual( NumType n1, NumType n2 ) { return tolEqual(n1,n2) || n1>n2; }
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/basegfx/source/workbench/convexhull.cxx b/basegfx/source/workbench/convexhull.cxx
new file mode 100644
index 000000000..ddcb04046
--- /dev/null
+++ b/basegfx/source/workbench/convexhull.cxx
@@ -0,0 +1,200 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * This file incorporates work covered by the following license notice:
+ *
+ *   Licensed to the Apache Software Foundation (ASF) under one or more
+ *   contributor license agreements. See the NOTICE file distributed
+ *   with this work for additional information regarding copyright
+ *   ownership. The ASF licenses this file to you under the Apache
+ *   License, Version 2.0 (the "License"); you may not use this file
+ *   except in compliance with the License. You may obtain a copy of
+ *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
+ */
+
+#include <algorithm>
+#include <vector>
+
+#include <bezierclip.hxx>
+
+/* Implements the theta function from Sedgewick: Algorithms in XXX, chapter 24 */
+template <class PointType> double theta( const PointType& p1, const PointType& p2 )
+{
+    typename PointType::value_type dx, dy, ax, ay;
+    double t;
+
+    dx = p2.x - p1.x; ax = absval( dx );
+    dy = p2.y - p1.y; ay = absval( dy );
+    t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
+    if( dx < 0 )
+        t = 2-t;
+    else if( dy < 0 )
+        t = 4+t;
+
+    return t*90.0;
+}
+
+/* Model of LessThanComparable for theta sort.
+ * Uses the theta function from Sedgewick: Algorithms in XXX, chapter 24
+ */
+template <class PointType> class ThetaCompare
+{
+public:
+    explicit ThetaCompare( const PointType& rRefPoint ) : maRefPoint( rRefPoint ) {}
+
+    bool operator() ( const PointType& p1, const PointType& p2 )
+    {
+        // return true, if p1 precedes p2 in the angle relative to maRefPoint
+        return theta(maRefPoint, p1) < theta(maRefPoint, p2);
+    }
+
+    double operator() ( const PointType& p ) const
+    {
+        return theta(maRefPoint, p);
+    }
+
+private:
+    PointType maRefPoint;
+};
+
+/* Implementation of the predicate 'counter-clockwise' for three points, from Sedgewick: Algorithms in XXX, chapter 24 */
+template <class PointType, class CmpFunctor> typename PointType::value_type ccw( const PointType& p0, const PointType& p1, const PointType& p2, CmpFunctor& thetaCmp )
+{
+    typename PointType::value_type dx1, dx2, dy1, dy2;
+    typename PointType::value_type theta0( thetaCmp(p0) );
+    typename PointType::value_type theta1( thetaCmp(p1) );
+    typename PointType::value_type theta2( thetaCmp(p2) );
+
+#if 0
+    if( theta0 == theta1 ||
+        theta0 == theta2 ||
+        theta1 == theta2 )
+    {
+        // cannot reliably compare, as at least two points are
+        // theta-equal. See bug description below
+        return 0;
+    }
+#endif
+
+    dx1 = p1.x - p0.x; dy1 = p1.y - p0.y;
+    dx2 = p2.x - p0.x; dy2 = p2.y - p0.y;
+
+    if( dx1*dy2 > dy1*dx2 )
+        return +1;
+
+    if( dx1*dy2 < dy1*dx2 )
+        return -1;
+
+    if( (dx1*dx2 < 0) || (dy1*dy2 < 0) )
+        return -1;
+
+    if( (dx1*dx1 + dy1*dy1) < (dx2*dx2 + dy2*dy2) )
+        return +1;
+
+    return 0;
+}
+
+/*
+ Bug
+ ===
+
+ Sometimes, the resulting polygon is not the convex hull (see below
+ for an edge configuration to reproduce that problem)
+
+ Problem analysis:
+ =================
+
+ The root cause of this bug is the fact that the second part of
+ the algorithm (the 'wrapping' of the point set) relies on the
+ previous theta sorting. Namely, it is required that the
+ generated point ordering, when interpreted as a polygon, is not
+ self-intersecting. This property, although, cannot be
+ guaranteed due to limited floating point accuracy. For example,
+ for two points very close together, and at the same time very
+ far away from the theta reference point p1, can appear on the
+ same theta value (because floating point accuracy is limited),
+ although on different rays to p1 when inspected locally.
+
+ Example:
+
+                    /
+             P3    /
+               |\ /
+               | /
+               |/ \
+             P2    \
+                    \
+      ...____________\
+     P1
+
+ Here, P2 and P3 are theta-equal relative to P1, but the local
+ ccw measure always says 'left turn'. Thus, the convex hull is
+ wrong at this place.
+
+ Solution:
+ =========
+
+ If two points are theta-equal and checked via ccw, ccw must
+ also classify them as 'equal'. Thus, the second stage of the
+ convex hull algorithm sorts the first one out, effectively
+ reducing a cluster of theta-equal points to only one. This
+ single point can then be treated correctly.
+*/
+
+/* Implementation of Graham's convex hull algorithm, see Sedgewick: Algorithms in XXX, chapter 25 */
+Polygon2D convexHull( const Polygon2D& rPoly )
+{
+    const Polygon2D::size_type N( rPoly.size() );
+    Polygon2D result( N + 1 );
+    std::copy(rPoly.begin(), rPoly.end(), result.begin()+1 );
+    Polygon2D::size_type min, i;
+
+    // determine safe point on hull (smallest y value)
+    for( min=1, i=2; i<=N; ++i )
+    {
+        if( result[i].y < result[min].y )
+            min = i;
+    }
+
+    // determine safe point on hull (largest x value)
+    for( i=1; i<=N; ++i )
+    {
+        if( result[i].y == result[min].y &&
+            result[i].x > result[min].x )
+        {
+            min = i;
+        }
+    }
+
+    // TODO: add inner elimination optimization from Sedgewick: Algorithms in XXX, chapter 25
+    // TODO: use radix sort instead of std::sort(), calc theta only once and store
+
+    // setup first point and sort
+    std::swap( result[1], result[min] );
+    ThetaCompare<Point2D> cmpFunc(result[1]);
+    std::sort( result.begin()+1, result.end(), cmpFunc );
+
+    // setup sentinel
+    result[0] = result[N];
+
+    // generate convex hull
+    Polygon2D::size_type M;
+    for( M=3, i=4; i<=N; ++i )
+    {
+        while( ccw(result[M], result[M-1], result[i], cmpFunc) >= 0 )
+            --M;
+
+        ++M;
+        std::swap( result[M], result[i] );
+    }
+
+    // copy range [1,M] to output
+    return Polygon2D( result.begin()+1, result.begin()+M+1 );
+}
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/basegfx/source/workbench/gauss.hxx b/basegfx/source/workbench/gauss.hxx
new file mode 100644
index 000000000..3605c1cac
--- /dev/null
+++ b/basegfx/source/workbench/gauss.hxx
@@ -0,0 +1,167 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * This file incorporates work covered by the following license notice:
+ *
+ *   Licensed to the Apache Software Foundation (ASF) under one or more
+ *   contributor license agreements. See the NOTICE file distributed
+ *   with this work for additional information regarding copyright
+ *   ownership. The ASF licenses this file to you under the Apache
+ *   License, Version 2.0 (the "License"); you may not use this file
+ *   except in compliance with the License. You may obtain a copy of
+ *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
+ */
+
+/** This method eliminates elements below main diagonal in the given
+    matrix by gaussian elimination.
+
+    @param matrix
+    The matrix to operate on. Last column is the result vector (right
+    hand side of the linear equation). After successful termination,
+    the matrix is upper triangular. The matrix is expected to be in
+    row major order.
+
+    @param rows
+    Number of rows in matrix
+
+    @param cols
+    Number of columns in matrix
+
+    @param minPivot
+    If the pivot element gets lesser than minPivot, this method fails,
+    otherwise, elimination succeeds and true is returned.
+
+    @return true, if elimination succeeded.
+ */
+
+#pragma once
+
+template <class Matrix, typename BaseType>
+bool eliminate(     Matrix&         matrix,
+                    int             rows,
+                    int             cols,
+                    const BaseType& minPivot    )
+{
+    BaseType    temp;
+
+    /* i, j, k *must* be signed, when looping like: j>=0 ! */
+    /* eliminate below main diagonal */
+    for(int i=0; i<cols-1; ++i)
+    {
+        /* find best pivot */
+        int max = i;
+        for(int j=i+1; j<rows; ++j)
+            if( fabs(matrix[ j*cols + i ]) > fabs(matrix[ max*cols + i ]) )
+                max = j;
+
+        /* check pivot value */
+        if( fabs(matrix[ max*cols + i ]) < minPivot )
+            return false;   /* pivot too small! */
+
+        /* interchange rows 'max' and 'i' */
+        for(int k=0; k<cols; ++k)
+        {
+            temp = matrix[ i*cols + k ];
+            matrix[ i*cols + k ] = matrix[ max*cols + k ];
+            matrix[ max*cols + k ] = temp;
+        }
+
+        /* eliminate column */
+        for(int j=i+1; j<rows; ++j)
+            for(int k=cols-1; k>=i; --k)
+                matrix[ j*cols + k ] -= matrix[ i*cols + k ] *
+                    matrix[ j*cols + i ] / matrix[ i*cols + i ];
+    }
+
+    /* everything went well */
+    return true;
+}
+
+/** Retrieve solution vector of linear system by substituting backwards.
+
+    This operation _relies_ on the previous successful
+    application of eliminate()!
+
+    @param matrix
+    Matrix in upper diagonal form, as e.g. generated by eliminate()
+
+    @param rows
+    Number of rows in matrix
+
+    @param cols
+    Number of columns in matrix
+
+    @param result
+    Result vector. Given matrix must have space for one column (rows entries).
+
+    @return true, if back substitution was possible (i.e. no division
+    by zero occurred).
+ */
+template <class Matrix, class Vector, typename BaseType>
+bool substitute(    const Matrix&   matrix,
+                    int             rows,
+                    int             cols,
+                    Vector&         result  )
+{
+    BaseType    temp;
+
+    /* j, k *must* be signed, when looping like: j>=0 ! */
+    /* substitute backwards */
+    for(int j=rows-1; j>=0; --j)
+    {
+        temp = 0.0;
+        for(int k=j+1; k<cols-1; ++k)
+            temp += matrix[ j*cols + k ] * result[k];
+
+        if( matrix[ j*cols + j ] == 0.0 )
+            return false;   /* imminent division by zero! */
+
+        result[j] = (matrix[ j*cols + cols-1 ] - temp) / matrix[ j*cols + j ];
+    }
+
+    /* everything went well */
+    return true;
+}
+
+/** This method determines solution of given linear system, if any
+
+    This is a wrapper for eliminate and substitute, given matrix must
+    contain right side of equation as the last column.
+
+    @param matrix
+    The matrix to operate on. Last column is the result vector (right
+    hand side of the linear equation). After successful termination,
+    the matrix is upper triangular. The matrix is expected to be in
+    row major order.
+
+    @param rows
+    Number of rows in matrix
+
+    @param cols
+    Number of columns in matrix
+
+    @param minPivot
+    If the pivot element gets lesser than minPivot, this method fails,
+    otherwise, elimination succeeds and true is returned.
+
+    @return true, if elimination succeeded.
+ */
+template <class Matrix, class Vector, typename BaseType>
+bool solve( Matrix&     matrix,
+            int         rows,
+            int         cols,
+            Vector&     result,
+            BaseType    minPivot    )
+{
+    if( eliminate<Matrix,BaseType>(matrix, rows, cols, minPivot) )
+        return substitute<Matrix,Vector,BaseType>(matrix, rows, cols, result);
+
+    return false;
+}
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */