Stroke to path implementation

11 files changed, 2598 insertions, 0 deletions

diff --git a/BIBLIOGRAPHY b/BIBLIOGRAPHY
index 90a6cef20..e3081da1c 100644
--- a/BIBLIOGRAPHY
+++ b/BIBLIOGRAPHY
@@ -44,6 +44,31 @@ Now, "convolve" the "tracing" of the pen with the tracing of the path:
         24th IEEE Annual Symposium on Foundations of Computer Science
         (FOCS), November 1983, 100-111.
 
+An alternative method involves computing an 'offset path' for the
+convolution of the pen along the path.
+
+	A good introduction is given by:
+
+	"Comparing Offset Curve Approximation Methods"
+	Gershon Elber, In-Kwon Lee and Myung-Soo Kim
+        http://visualcomputing.yonsei.ac.kr/papers/1997/compare.pdf
+
+	The algorithm suggested by this paper is:
+
+	"Offsets of Two-Dimensional Profiles"
+	Wayne Tiller and Eric G. Hanson
+	IEEE Computer Graphics & Application, Vol 4. pp. 36-46
+	
+	The development of the symbolic global error
+        and iterative refinement is given by:
+
+	"Free Form Surface Analysis using a Hybrid of Symbolic and Numeric
+        Computation."
+	Gershon Elber, 1992
+	http://www.cs.utah.edu/gdc/publications/papers/elber_thesis.pdf
+
+
+
 The result of the convolution is a polygon that must be filled. A fill
 operations begins here. We use a very conventional Bentley-Ottmann
 pass for computing the intersections, informed by some hints on robust
diff --git a/src/Makefile.sources b/src/Makefile.sources
index b04d80fca..93f7df383 100644
--- a/src/Makefile.sources
+++ b/src/Makefile.sources
@@ -159,7 +159,9 @@ cairo_sources = \
 	cairo-slope.c \
 	cairo-spans.c \
 	cairo-spline.c \
+	cairo-spline-offset.c \
 	cairo-stroke-style.c \
+	cairo-stroke-to-path.c \
 	cairo-surface.c \
 	cairo-surface-fallback.c \
 	cairo-surface-clipper.c \
diff --git a/src/cairo-path-fixed.c b/src/cairo-path-fixed.c
index eea8630bd..44e6e50bf 100644
--- a/src/cairo-path-fixed.c
+++ b/src/cairo-path-fixed.c
@@ -177,6 +177,7 @@ _cairo_path_fixed_init_copy (cairo_path_fixed_t *path,
     return CAIRO_STATUS_SUCCESS;
 }
 
+
 unsigned long
 _cairo_path_fixed_hash (const cairo_path_fixed_t *path)
 {
@@ -919,6 +920,23 @@ _cairo_path_fixed_append (cairo_path_fixed_t		    *path,
 					&closure);
 }
 
+cairo_status_t
+_cairo_path_fixed_init_flat_copy (cairo_path_fixed_t *path,
+	const cairo_path_fixed_t *other,
+	double tolerance)
+{
+
+    _cairo_path_fixed_init (path);
+    return _cairo_path_fixed_interpret_flat (other,
+	                                      CAIRO_DIRECTION_FORWARD,
+					      _append_move_to,
+					      _append_line_to,
+					      _append_close_path,
+					      path,
+					      tolerance);
+}
+
+
 static void
 _cairo_path_fixed_offset_and_scale (cairo_path_fixed_t *path,
 				    cairo_fixed_t offx,
diff --git a/src/cairo-path-private.h b/src/cairo-path-private.h
index 61b4060fa..964fdd989 100644
--- a/src/cairo-path-private.h
+++ b/src/cairo-path-private.h
@@ -54,4 +54,8 @@ cairo_private cairo_status_t
 _cairo_path_append_to_context (const cairo_path_t	*path,
 			       cairo_t			*cr);
 
+cairo_private cairo_path_t *
+_cairo_stroked_path_create (cairo_path_fixed_t *path_fixed,
+			     cairo_gstate_t     *gstate);
+
 #endif /* CAIRO_PATH_DATA_PRIVATE_H */
diff --git a/src/cairo-path.c b/src/cairo-path.c
index 28182c0e4..f6e59a1e0 100644
--- a/src/cairo-path.c
+++ b/src/cairo-path.c
@@ -369,6 +369,58 @@ _cairo_path_create_internal (cairo_path_fixed_t *path_fixed,
     return path;
 }
 
+cairo_path_t *
+_cairo_stroked_path_create (cairo_path_fixed_t *path_fixed,
+			     cairo_gstate_t     *gstate)
+{
+    cairo_path_t *path;
+
+    path = malloc (sizeof (cairo_path_t));
+    if (unlikely (path == NULL)) {
+	_cairo_error_throw (CAIRO_STATUS_NO_MEMORY);
+	return (cairo_path_t*) &_cairo_path_nil;
+    }
+
+    cairo_path_fixed_t stroke_path;
+
+    _cairo_path_fixed_init (&stroke_path);
+
+    _cairo_path_stroke_to_path (path_fixed,
+	    &gstate->stroke_style,
+	    &stroke_path,
+	    &gstate->ctm,
+	    &gstate->ctm_inverse,
+	    gstate->tolerance);
+
+    path->num_data = _cairo_path_count (path, &stroke_path,
+					gstate->tolerance,
+					FALSE);
+    if (path->num_data < 0) {
+	free (path);
+	_cairo_path_fixed_fini (&stroke_path);
+	return (cairo_path_t*) &_cairo_path_nil;
+    }
+
+    if (path->num_data) {
+	path->data = _cairo_malloc_ab (path->num_data,
+				       sizeof (cairo_path_data_t));
+	if (unlikely (path->data == NULL)) {
+	    free (path);
+	    _cairo_error_throw (CAIRO_STATUS_NO_MEMORY);
+	    return (cairo_path_t*) &_cairo_path_nil;
+	}
+
+	path->status = _cairo_path_populate (path, &stroke_path,
+					     gstate, FALSE);
+    } else {
+	path->data = NULL;
+	path->status = CAIRO_STATUS_SUCCESS;
+    }
+    _cairo_path_fixed_fini (&stroke_path);
+
+    return path;
+}
+
 /**
  * cairo_path_destroy:
  * @path: a path previously returned by either cairo_copy_path() or
diff --git a/src/cairo-spline-offset.c b/src/cairo-spline-offset.c
new file mode 100644
index 000000000..283ed37db
--- /dev/null
+++ b/src/cairo-spline-offset.c
@@ -0,0 +1,945 @@
+/* -*- Mode: c; tab-width: 8; c-basic-offset: 4; indent-tabs-mode: t; -*- */
+/* cairo - a vector graphics library with display and print output
+ *
+ * Copyright © 2009 Mozilla Foundation
+ * Copyright © 2010 Jeff Muizelaar
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (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.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ * The Original Code is the cairo graphics library.
+ *
+ * The Initial Developer of the Original Code is Jeff Muizelaar
+ *
+ */
+
+#define _GNU_SOURCE
+
+#include "cairoint.h"
+
+#include <math.h>
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+
+#if _XOPEN_SOURCE >= 600 || defined (_ISOC99_SOURCE)
+#define ISFINITE(x) isfinite (x)
+#else
+#define ISFINITE(x) ((x) * (x) >= 0.) /* check for NaNs */
+#endif
+
+#include "cairoint.h"
+#include "cairo-spline-offset.h"
+
+static void
+_lerp (const point_t *a, const point_t *b, point_t *result, double t)
+{
+	result->x = (1-t)*a->x + t*b->x;
+	result->y = (1-t)*a->y + t*b->y;
+}
+
+/* returns false if there could be a discontinuity
+ * between the split curves */
+static cairo_bool_t
+_de_casteljau_t (knots_t *k1, knots_t *k2, double t)
+{
+	point_t ab, bc, cd;
+	point_t abbc, bccd;
+	point_t final;
+
+	_lerp (&k1->a, &k1->b, &ab, t);
+	_lerp (&k1->b, &k1->c, &bc, t);
+	_lerp (&k1->c, &k1->d, &cd, t);
+	_lerp (&ab, &bc, &abbc, t);
+	_lerp (&bc, &cd, &bccd, t);
+	_lerp (&abbc, &bccd, &final, t);
+
+	k2->a = final;
+	k2->b = bccd;
+	k2->c = cd;
+	k2->d = k1->d;
+
+	k1->b = ab;
+	k1->c = abbc;
+	k1->d = final;
+
+    /* we can test whether the split curve will have a cusp
+     * inbetween the two parts by checking the length of the
+     * line abbc-bccd. If this line has a non-zero length
+     * then the two curves will share a common normal
+     * at the points where they meet. Therefore the
+     * length must be zero for the normals to be different.
+     * However this condition is necessary but not sufficient.
+     * XXX: does this matter and how should we deal with it if it does? */
+    return abbc.x != final.x || abbc.y != final.y;
+}
+
+typedef point_t vector_t;
+
+static vector_t
+begin_tangent (knots_t c)
+{
+    vector_t t;
+    if (c.a.x != c.b.x || c.a.y != c.b.y) {
+	t.x = c.b.x - c.a.x;
+	t.y = c.b.y - c.a.y;
+	return t;
+    }
+    if (c.a.x != c.c.x || c.a.y != c.c.y) {
+	t.x = c.c.x - c.a.x;
+	t.y = c.c.y - c.a.y;
+	return t;
+    }
+    t.x = c.d.x - c.a.x;
+    t.y = c.d.y - c.a.y;
+    return t;
+}
+
+static vector_t
+end_tangent (knots_t c)
+{
+    vector_t t;
+    if (c.d.x != c.c.x || c.d.y != c.c.y) {
+	t.x = c.d.x - c.c.x;
+	t.y = c.d.y - c.c.y;
+	return t;
+    }
+    if (c.d.x != c.b.x || c.d.y != c.b.y) {
+	t.x = c.d.x - c.b.x;
+	t.y = c.d.y - c.b.y;
+	return t;
+    }
+    t.x = c.d.x - c.a.x;
+    t.y = c.d.y - c.a.y;
+    return t;
+}
+
+static void
+ensure_finite (knots_t k)
+{
+
+    if (!(
+	    isfinite(k.a.x) &&
+	    isfinite(k.a.y) &&
+	    isfinite(k.b.x) &&
+	    isfinite(k.b.y) &&
+	    isfinite(k.c.x) &&
+	    isfinite(k.c.y) &&
+	    isfinite(k.d.x) &&
+	    isfinite(k.d.y))) {
+	print_knot(k);
+	assert(0);
+    }
+}
+
+#if 0
+static void ensure_smoothness (knots_t a, knots_t b)
+{
+    double dist_len = hypot(a.d.x - b.a.x, a.d.y - b.a.y);
+    print_knot(a);
+    print_knot(b);
+    printf("dist: %f\n", dist_len);
+    vector_t in_norm = end_tangent(a);
+    vector_t out_norm = begin_tangent(b);
+    double in_a = atan2 (in_norm.y, in_norm.x);
+    double out_a = atan2 (out_norm.y, out_norm.x);
+    if (in_a < out_a)
+	in_a += 2*M_PI;
+    double diff = in_a - out_a;
+    if (diff > M_PI)
+	diff -= M_PI*2;
+    if (!(fabs(diff) < 0.001) && !(fabs(diff - M_PI) < 0.001) && !(fabs(diff + M_PI) < 0.001)) {
+	/* The angle difference between normals can be either 0 or +-M_PI for a change in direction.
+	 * The change in direction can occur in a region of high curvature with a large offset direction.
+	 * We move in reverse if the offset distance exceeds the radius of curvature. Note: that you need to
+	 * use signed curvature for this too work out properly. Because we only encounter this phenomenom in
+	 * the convex areas of the curve.
+	 *
+	 * XXX: explain why... */
+	printf("angle diff: %f %f %f\n", diff, in_a, out_a);
+	print_knot(a);
+	print_knot(b);
+	assert(0);
+    }
+    if (fabs(dist_len) > 0.01) {
+	printf("distlen: %f\n", dist_len);
+	print_knot(a);
+	print_knot(b);
+	assert(0);
+    }
+}
+
+static void
+ensure_smoothness2 (knots_t a, knots_t b, cairo_bool_t ccw) {
+    vector_t in_norm = end_tangent(a);
+    vector_t out_norm = begin_tangent(b);
+    if (!ccw) {
+	in_norm = end_tangent(a);
+	out_norm = begin_tangent(b);
+    }
+    double in_a = atan2 (in_norm.y, in_norm.x);
+    double out_a = atan2 (out_norm.y, out_norm.x);
+    if (in_a < out_a)
+	in_a += 2*M_PI;
+    double diff = in_a - out_a;
+    if (diff > M_PI)
+	diff -= M_PI*2;
+    if (fabs(diff) > 0.01) {
+	printf("gangle diff: %f\n", diff);
+	print_knot(a);
+	print_knot(b);
+    }
+}
+#endif
+
+static inline vector_t
+flip (vector_t a)
+{
+    vector_t ar;
+    ar.x = -a.x;
+    ar.y = -a.y;
+    return ar;
+}
+
+static vector_t
+perp (vector_t v)
+{
+    vector_t p = {-v.y, v.x};
+    return p;
+}
+
+static inline double
+dot (vector_t a, vector_t b)
+{
+    return a.x * b.x + a.y * b.y;
+}
+
+static vector_t
+normalize (vector_t v)
+{
+    double len = hypot(v.x, v.y);
+    v.x /= len;
+    v.y /= len;
+    return v;
+}
+
+/* given a normal rotate the vector 90 degrees to the right clockwise
+ * This function has a period of 4. e.g. swap(swap(swap(swap(x) == x */
+static vector_t
+swap (vector_t a)
+{
+    vector_t ar;
+    /* one of these needs to be negative. We choose a.x so that we rotate to the right instead of negating */
+    ar.x = a.y;
+    ar.y = -a.x;
+    return ar;
+}
+
+static vector_t
+unperp (vector_t a)
+{
+    return swap (a);
+}
+
+typedef struct {
+    double t1;
+    double t2;
+} inflection_points_t;
+
+/* Find the inflection points for a knots_t.
+ * The lower inflection point is returned in t1, the higher one in t2.
+ *
+ * This method for computing inflections points is from:
+ * "Fast, precise flattening of cubic Bezier path and offset curves", Hain et. al */
+static inflection_points_t
+inflection_points (knots_t s)
+{
+    inflection_points_t ret;
+
+    point_t a = {x:  -s.a.x + 3*s.b.x - 3*s.c.x + s.d.x,  -s.a.y + 3*s.b.y - 3*s.c.y + s.d.y};
+    point_t b = {x: 3*s.a.x - 6*s.b.x + 3*s.c.x,         3*s.a.y - 6*s.b.y + 3*s.c.y};
+    point_t c = {x:-3*s.a.x + 3*s.b.x,                  -3*s.a.y + 3*s.b.y};
+    /* we don't need 'd'
+    point_t d = {x:   s.a.x,                               s.a.y};
+    */
+
+    /* we're solving for the roots of this equation:
+        dx/dt * d2y/dt2 - d2x/dt2 * dy/dt
+        == 6*(a.y*b.x - a.x*b.y)*t^2 +6*(a.y*c.x-a.x*c.y)*t + 2*(b.y*c.x - b.x*c.y)
+    */
+    double t_cusp = (-1./2)*((a.y*c.x - a.x*c.y)/(a.y*b.x - a.x*b.y));
+    double t_inner = (t_cusp*t_cusp - ((b.y*c.x - b.x*c.y)/(a.y*b.x - a.x*b.y))/3);
+
+    if (t_inner > 0) {
+	double t_adj = sqrt (t_inner);
+	ret.t1 = t_cusp - t_adj;
+	ret.t2 = t_cusp + t_adj;
+    } else {
+	ret.t1 = ret.t2 = t_cusp;
+    }
+    return ret;
+}
+
+/* this function could use some tests.
+   I believe it only works for angles from 0 to pi. It will always use the smaller arc between the two vectors? */
+static knots_t
+arc_segment_cart (point_t start, point_t end, double radius, vector_t a, vector_t b)
+{
+    knots_t arc;
+    double r_sin_A = radius * a.y;
+    double r_cos_A = radius * a.x;
+    double r_sin_B = radius * b.y;
+    double r_cos_B = radius * b.x;
+
+    point_t mid, mid2;
+    double l, h;
+
+    mid.x = a.x + b.x;
+    mid.y = a.y + b.y;
+    if (dot (a, b) < 0) {
+	/* otherwise, we can flip a, add it
+	 * and then use the perpendicular of the result */
+	a = flip (a);
+	mid.x = a.x + b.x;
+	mid.y = a.y + b.y;
+	if (dot (a, perp(b)) <= 0) {
+	    mid = unperp (mid);
+	} else {
+	    mid = perp(mid);
+	}
+    }
+    /* normalize */
+    l = sqrt(mid.x*mid.x + mid.y*mid.y);
+    mid.x /= l;
+    mid.y /= l;
+
+    mid2.x = a.x + mid.x;
+    mid2.y = a.y + mid.y;
+
+    h = (4./3.)*dot(perp(a), mid2)/dot(a, mid2);
+
+    arc.a.x = start.x;
+    arc.a.y = start.y;
+
+    arc.b.x = start.x - h * r_sin_A;
+    arc.b.y = start.y + h * r_cos_A;
+
+    arc.c.x =   end.x + h * r_sin_B;
+    arc.c.y =   end.y - h * r_cos_B;
+
+    arc.d.x = end.x;
+    arc.d.y = end.y;
+
+    return arc;
+}
+
+static knots_t
+quarter_circle (double xc, double yc, double radius, vector_t normal)
+{
+    double r_sin_A, r_cos_A;
+    double r_sin_B, r_cos_B;
+    double h;
+    knots_t arc;
+
+    r_sin_A = radius * normal.y;
+    r_cos_A = radius * normal.x;
+    r_sin_B = radius * -normal.x; /* sin (angle_B); */
+    r_cos_B = radius * normal.y;  /* cos (angle_B); */
+
+    h = -0.55228474983079345; /* 4.0/3.0 * tan ((-M_PI/2) / 4.0); */
+
+    arc.a.x = xc + r_cos_A;
+    arc.a.y = yc + r_sin_A;
+
+    arc.b.x = xc + r_cos_A - h * r_sin_A;
+    arc.b.y = yc + r_sin_A + h * r_cos_A;
+
+    arc.c.x = xc + r_cos_B + h * r_sin_B;
+    arc.c.y = yc + r_sin_B - h * r_cos_B;
+
+    arc.d.x = xc + r_cos_B;
+    arc.d.y = yc + r_sin_B;
+
+    return arc;
+}
+
+static void semi_circle (double xc, double yc, double radius, vector_t out_normal, void (*curve_fn)(void *, knots_t), void *closure)
+{
+    double len = hypot (out_normal.y, out_normal.x);
+    vector_t mid_normal;
+
+    out_normal.x /= len;
+    out_normal.y /= len;
+
+    mid_normal.x = out_normal.y; /* we need to rotate by -M_PI/2 */
+    mid_normal.y = -out_normal.x;
+
+    curve_fn(closure, quarter_circle (xc, yc, radius, out_normal));
+    curve_fn(closure, quarter_circle (xc, yc, radius, mid_normal));
+}
+
+static cairo_bool_t
+is_knot_flat (knots_t c)
+{
+    /*
+       A well-known flatness test is both cheaper and more reliable
+       than the ones you have tried. The essential observation is that
+       when the curve is a uniform speed straight line from end to end,
+       the control points are evenly spaced from beginning to end.
+       Therefore, our measure of how far we deviate from that ideal
+       uses distance of the middle controls, not from the line itself,
+       but from their ideal *arrangement*. Point 2 should be halfway
+       between points 1 and 3; point 3 should be halfway between points
+       2 and 4.
+       This, too, can be improved. Yes, we can eliminate the square roots
+       in the distance tests, retaining the Euclidean metric; but the
+       taxicab metric is faster, and also safe. The length of displacement
+       (x,y) in this metric is |x|+|y|.
+    */
+
+    /* XXX: this isn't necessarily the best test because it
+       tests for flatness instead of smallness.
+       however, flat things will be offseted easily so we sort of implicitly do the right thing
+    */
+    /* this value should be much less than the tolerance because we only want to deal with situations where
+     * the curvature is extreme. */
+    double stop_value = 0.05;
+    double error = 0;
+    error += fabs((c.c.x - c.b.x) - (c.b.x - c.a.x));
+    error += fabs((c.c.y - c.b.y) - (c.b.y - c.a.y));
+    error += fabs((c.c.x - c.b.x) - (c.d.x - c.c.x));
+    error += fabs((c.c.y - c.b.y) - (c.d.y - c.c.y));
+    return error <= stop_value;
+}
+
+/* Finds the coefficient with the maximum error. This should be
+ * the coefficient closest to area of maximum error.
+ * This returns an upper bound of the error. */
+static cairo_bool_t
+is_curvy_t (double *bezier, int degree, double stop_value, double *t)
+{
+	double highest_error = .5;
+	double error = 0;
+	int i=0;
+	assert(degree == 6);
+	for (i=0; i<=degree; i++) {
+	    if (fabs(bezier[i]) > error) {
+		error = fabs(bezier[i]);
+		highest_error = i/(double)degree;
+	    }
+	}
+	*t = highest_error;
+
+	/* we shouldn't ever split at the edges because the approximation
+	 * should be perfect there. XXX: we might be able to take
+	 * advantage of this when computing the error polynomial. */
+#if 0
+	if (*t == 0 || *t == 1) {
+	    for (i=0; i<=degree; i++) {
+		error = fabs(bezier[i]);
+		printf("error: %f\n", error);
+	    }
+	    assert(0);
+	}
+#endif
+	/* *t = 0.5; */
+	return error >= stop_value;
+}
+
+static cairo_bool_t
+error_within_bounds_elber (knots_t bz1, knots_t bz_offset_approx,
+		double desired_offset,
+		double max_error,
+		double *t)
+{
+    double bzprod[3+3 +1] = {};
+
+    /* Elber and Cohnen propose the following method of computing the
+       error between an curve and it's offset:
+       e(t) = || C(t) - Coffset_approx(t) ||^2 - offset^2
+
+       It is important to note that this function doesn't represent a unique offset curve.
+       For example, it's possible for the offset curve to cross the original spline
+       provided it stays the appropriate distance away:
+
+                 xxx
+                x
+               ****
+             **x
+            *  x
+           *   x
+            *  x
+             **
+               ****
+               x
+                xxx
+
+      "Planar Curve Offset Base on Circle Approximation" suggests further problems
+      with this error measure. They give the following example: C(t) + (r, 0) will
+      have an error mesaure of 0, but is not even close to the actual offset curve.
+
+      This paper suggests using the angle between the difference vector and the normal
+      at point t, but doesn't go into much detail.
+     */
+
+    /* Now we compute this error function.
+     * First, subtract the approx */
+    bz1.a.x -= bz_offset_approx.a.x;
+    bz1.a.y -= bz_offset_approx.a.y;
+    bz1.b.x -= bz_offset_approx.b.x;
+    bz1.b.y -= bz_offset_approx.b.y;
+    bz1.c.x -= bz_offset_approx.c.x;
+    bz1.c.y -= bz_offset_approx.c.y;
+    bz1.d.x -= bz_offset_approx.d.x;
+    bz1.d.y -= bz_offset_approx.d.y;
+/*
+     Next we compute the dot product of bz1 with itself This involves squaring
+     the x and y polynomials and adding the result together.
+
+     However, we want to multiply the coefficients so that we can treat the
+     result as a higher degree bezier curve.
+
+     The following method is given in the Symbolic Computation section of
+     "Comparing Offset Curve Approximation Methods" by Elber, Lee and Kim
+
+  k     l       ⎛ k⎞  i     k-i⎛ l⎞  j       l-j
+ B (t) B (t) =  ⎜  ⎟ t (1-t)   ⎜  ⎟ t   (1-t)
+  i     j       ⎝ i⎠           ⎝ j⎠
+	        ⎛ k⎞ ⎛ l⎞  i+j     k+l-i-j
+	     =  ⎜  ⎟ ⎜  ⎟ t   (1-t)
+	        ⎝ i⎠ ⎝ j⎠
+
+	        ⎛ k⎞ ⎛ l⎞
+	     =  ⎝ i⎠ ⎝ j⎠  k+l
+	         ──────── B   (t)
+	        ⎛ k + l⎞   i+j
+	        ⎝ i + j⎠
+
+               m               n
+C1(t) C2(t) =  ∑ P_i B_i,m (t) ∑ P_j B_j,n (t)
+              i=0             j=0
+
+	       m   n
+            =  ∑   ∑ P_i P_j B_i,m (t) B_j,n (t)
+              i=0 j=0
+
+	       m+n       m + n
+	    =   ∑ Q_k B_k     (t)
+	       k=0
+
+Where Q_k is:
+		       ⎛ m⎞ ⎛ n⎞
+		       ⎝ i⎠ ⎝ j⎠
+       Q_k = P_i * P_j ─────────
+			⎛ m+n⎞
+			⎝ i+j⎠
+
+And for the case of degree 3 beziers
+
+		       ⎛ 3⎞ ⎛ 3⎞
+		       ⎝ i⎠ ⎝ j⎠
+       Q_k = P_i * P_j ─────────
+			⎛  6 ⎞
+			⎝ i+j⎠
+
+      This expands to:
+
+	i j
+	0 0 0 = (1 * 1) / 1  = 1
+
+	0 1 1 = (1 * 3) / 6  = .5
+	1 0 1 = (3 * 1) / 6  = .5
+
+	0 2 2 = (1 * 3) / 15 = .2
+	1 1 2 = (3 * 3) / 15 = .6
+	2 0 2 = (3 * 1) / 15 = .2
+
+	0 3 3 = (1 * 1) / 20 = .05
+	1 2 3 = (3 * 3) / 20 = .45
+	2 1 3 = (3 * 3) / 20 = .45
+	3 0 3 = (1 * 1) / 20 = .05
+
+	1 3 4 = (3 * 1) / 15 = .2
+	2 2 4 = (3 * 3) / 15 = .6
+	3 1 4 = (1 * 3) / 15 = .2
+
+	2 3 5 = (3 * 1) / 6  = .5
+	3 2 5 = (1 * 3) / 6  = .5
+
+	3 3 6 = (1 * 1) / 1  = 1
+
+Which simplifies to the following:
+	*/
+
+
+	bzprod[0] += bz1.a.x * bz1.a.x;
+	bzprod[0] += bz1.a.y * bz1.a.y;
+
+	bzprod[1] += bz1.a.x * bz1.b.x;
+	bzprod[1] += bz1.a.y * bz1.b.y;
+
+	bzprod[2] += bz1.a.x * bz1.c.x * 0.4;
+	bzprod[2] += bz1.b.x * bz1.b.x * 0.6;
+	bzprod[2] += bz1.a.y * bz1.c.y * 0.4;
+	bzprod[2] += bz1.b.y * bz1.b.y * 0.6;
+
+	bzprod[3] += bz1.d.x * bz1.a.x * 0.1;
+	bzprod[3] += bz1.b.x * bz1.c.x * 0.9;
+	bzprod[3] += bz1.d.y * bz1.a.y * 0.1;
+	bzprod[3] += bz1.b.y * bz1.c.y * 0.9;
+
+	bzprod[4] += bz1.b.x * bz1.d.x * 0.4;
+	bzprod[4] += bz1.c.x * bz1.c.x * 0.6;
+	bzprod[4] += bz1.d.y * bz1.b.y * 0.4;
+	bzprod[4] += bz1.c.y * bz1.c.y * 0.6;
+
+	bzprod[5] += bz1.c.x * bz1.d.x;
+	bzprod[5] += bz1.d.y * bz1.c.y;
+
+	bzprod[6] += bz1.d.x * bz1.d.x;
+	bzprod[6] += bz1.d.y * bz1.d.y;
+
+
+	/* the result is a bezier polynomial that represents the distance squared between
+	 * the two input polynomials */
+
+	/* we don't need to subtract the desired offset because our metric doesn't depend on being centered around 0 */
+#if 1
+	bzprod[0] -= (desired_offset*desired_offset);
+	bzprod[1] -= (desired_offset*desired_offset);
+	bzprod[2] -= (desired_offset*desired_offset);
+	bzprod[3] -= (desired_offset*desired_offset);
+	bzprod[4] -= (desired_offset*desired_offset);
+	bzprod[5] -= (desired_offset*desired_offset);
+	bzprod[6] -= (desired_offset*desired_offset);
+#endif
+	*t = 0.5;
+	return !is_curvy_t(bzprod, sizeof(bzprod)/sizeof(bzprod[0])-1, max_error, t);
+}
+
+void
+print_knot (knots_t c) {
+	/*printf("(%.13a %.13a) (%.13a %.13a) (%.13a %.13a) (%.13a %.13a)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);*/
+	printf("(%.5f %.5f) (%.5f %.5f) (%.5f %.5f) (%.5f %.5f)\n", c.a.x, c.a.y, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);
+}
+
+static knots_t
+convolve_with_circle (knots_t spline, double dist)
+{
+    vector_t in_normal  = normalize(perp(begin_tangent(spline)));
+    vector_t out_normal = normalize(perp(end_tangent(spline)));
+    /* find an arc the goes from the input normal to the output_normal */
+    knots_t arc_spline = arc_segment_cart(in_normal, out_normal, 1, in_normal, out_normal);
+
+    /* approximate convolving 'spline' with the arc spline (XXX: is convolve the right word?) */
+    spline.a.x += dist * arc_spline.a.x;
+    spline.a.y += dist * arc_spline.a.y;
+    spline.b.x += dist * arc_spline.b.x;
+    spline.b.y += dist * arc_spline.b.y;
+    spline.c.x += dist * arc_spline.c.x;
+    spline.c.y += dist * arc_spline.c.y;
+    spline.d.x += dist * arc_spline.d.x;
+    spline.d.y += dist * arc_spline.d.y;
+
+    return spline;
+}
+
+/* This approach is inspired by "Approximation of circular arcs and offset curves by Bezier curves of high degree"
+ * by YJ Ahn, Y soo Kim, Y Shin.
+ *
+ * This is a similar idea to:
+ *  "Offset approximation based on reparameterizaing the path of a moving point along the base circle", ZHAO, WANG
+ *  and a bunch of papers by Lee et al.
+ * "Planar curve offset based on circle approximation"
+ * "New approximation methods of planar offset and convolution curves"
+ * "Polynomial/rational approximation of Minkowski sum boundary curves"
+ *
+ * However these other approaches produce curves of higher degree than the input curve making them unsuitable
+ * for use here.
+ * */
+static void
+approximate_offset_curve_with_shin (knots_t self, double dist, void (*curve_fn)(void *, knots_t), void *closure)
+{
+    /* since we're going to approximating the offset curve with arcs
+       we want to ensure that we don't have any inflection points in our spline */
+    inflection_points_t t_inflect = inflection_points(self);
+
+    /* split the curve at the inflection points and convolve each with an approximation of
+     * a circle */
+    knots_t remaining = self;
+    double adjusted_inflect = t_inflect.t2;
+
+    /* check that the first inflection point is in range */
+    if (t_inflect.t1 > 0 && t_inflect.t1 < 1) {
+	    /* split at the inflection point */
+	    _de_casteljau_t (&self, &remaining, t_inflect.t1);
+
+	    /* approximate */
+	    curve_fn(closure, convolve_with_circle(self, dist));
+
+	    /* reparamaterize the second inflection point over the 'remaining' spline:
+	     * (domain between inflection points)/(domain of remaining) */
+	    adjusted_inflect = (t_inflect.t2-t_inflect.t1)/(1 - t_inflect.t1);
+    }
+
+    /* check that the second inflection point is in range and not equal to t1.
+     * we don't use the reparameterized value so that test remains simple */
+    if (t_inflect.t2 > t_inflect.t1 && t_inflect.t2 > 0 && t_inflect.t2 < 1) {
+	    /* split into one segment */
+	    knots_t self = remaining;
+	    _de_casteljau_t (&self, &remaining, adjusted_inflect);
+	    curve_fn(closure, convolve_with_circle(self, dist));
+    }
+
+    /* deal with what's left of the spline */
+    curve_fn(closure, convolve_with_circle(remaining, dist));
+}
+
+/* Computes the offset polygon for the control polygon of spline. We can
+ * use this is an approximation of the offset spline. This approximation is
+ * give by Tiller W, Hanson E G. in "Offsets of two-dimensional profile" */
+static knots_t
+knot_offset (knots_t self, double width, cairo_bool_t *is_parallel)
+{
+    /* this code is inspired by the libart mitreing code */
+
+    /* difference vectors between each point in the knot */
+    double dx0, dy0, dx1, dy1, dx2, dy2;
+    /* difference vector lengths */
+    double ld0, ld1, ld2;
+    /* temporary vector for dealing with degeneracies */
+    double last_x, last_y;
+
+    double scale;
+
+    /* normalized normal vectors for each difference vector */
+    double dlx0, dly0,  dlx1, dly1,  dlx2, dly2;
+
+    /* bisected/midpoint vectors */
+    double dm1x, dm1y,  dm2x, dm2y;
+
+    double dmr2_1, dmr2_2;
+
+    knots_t result;
+
+    dx0 = self.b.x - self.a.x;
+    dy0 = self.b.y - self.a.y;
+    ld0 = sqrt(dx0*dx0 + dy0*dy0);
+
+    dx1 = self.c.x - self.b.x;
+    dy1 = self.c.y - self.b.y;
+    ld1 = sqrt(dx1*dx1 + dy1*dy1);
+
+    dx2 = self.d.x - self.c.x;
+    dy2 = self.d.y - self.c.y;
+    ld2 = sqrt(dx2*dx2 + dy2*dy2);
+
+    /* compute normalized normals */
+    scale = 1. / ld0;
+    dlx0 = -dy0 * scale;
+    dly0 = dx0 * scale;
+
+    scale = 1. / ld1;
+    dlx1 = -dy1 * scale;
+    dly1 = dx1 * scale;
+
+    scale = 1. / ld2;
+    dlx2 = -dy2 * scale;
+    dly2 = dx2 * scale;
+
+    /* deal with any degeneracies in the spline by treating
+     * degenerate segments as having the same normal
+     * as their neighbour. If all of the segments are degenerate
+     * than we will fail the is_parallel test below */
+    last_x = 0;
+    last_y = 0;
+
+    /* first pass */
+    if (ld0) {
+	last_x = dlx0;
+	last_y = dly0;
+    }
+    if (ld1) {
+	last_x = dlx1;
+	last_y = dly1;
+    }
+    if (ld2) {
+	last_x = dlx2;
+	last_y = dly2;
+    }
+
+    /* second pass */
+    if (!ld2) {
+	dlx2 = last_x;
+	dly2 = last_y;
+    } else {
+	last_x = dlx2;
+	last_y = dly2;
+    }
+    if (!ld1) {
+	dlx1 = last_x;
+	dly1 = last_y;
+    } else {
+	last_x = dlx1;
+	last_y = dly1;
+    }
+    if (!ld0) {
+	dlx0 = last_x;
+	dly0 = last_y;
+    } else {
+	last_x = dlx0;
+	last_y = dly0;
+    }
+
+    /* mid-point vector sum */
+    dm1x = (dlx0 + dlx1);
+    dm1y = (dly0 + dly1);
+
+    dm2x = (dlx1 + dlx2);
+    dm2y = (dly1 + dly2);
+
+    /* length squared of the mid-point vector sum */
+    dmr2_1 = (dm1x * dm1x + dm1y * dm1y);
+    dmr2_2 = (dm2x * dm2x + dm2y * dm2y);
+
+    /* the choice of this EPSILON is arbitrary and could use further study */
+#define PARALLEL_EPSILON 1e-6
+    if (fabs(dmr2_1) < PARALLEL_EPSILON || fabs(dmr2_2) < PARALLEL_EPSILON)
+	*is_parallel = TRUE;
+    else
+	*is_parallel = FALSE;
+
+    scale = width * 2 / dmr2_1;
+    dm1x *= scale;
+    dm1y *= scale;
+
+    scale = width * 2 / dmr2_2;
+    dm2x *= scale;
+    dm2y *= scale;
+
+    /* write out the result */
+    result.a.x = self.a.x + dlx0 * width;
+    result.a.y = self.a.y + dly0 * width;
+
+    result.b.x = self.b.x + dm1x;
+    result.b.y = self.b.y + dm1y;
+
+    result.c.x = self.c.x + dm2x;
+    result.c.y = self.c.y + dm2y;
+
+    result.d.x = self.d.x + dlx2 * width;
+    result.d.y = self.d.y + dly2 * width;
+
+    return result;
+}
+
+/* plan:
+ * if we approximate the original bezier curve with line segments
+ * and then produce an offset to those lines. We will end up
+ * with sharp angles at places of high curvature.
+ *
+ * It is at these places of high curvature that we want to approximate
+ * the offset with an arc instead of continuing to subdivide.
+ *
+ * How do we find the areas of high curvature?
+ * - Ideas:
+ *   As we subdivide the original bezier
+ *   we can get an idea of how flat it is.
+ *   If the original bezier is flat but the offset approximation is not within our bounds
+ *   we should stop recursing and approximate the segment with an arc.
+ *
+ *   The degree of required flatness of the original curve should be related to the
+ *   the offset length. The greater the offset distance the less flat the original
+ *   bezier needs to be.
+ *
+ *   pros:
+ *   - this gives us an easier condition to reason about the termination of the
+ *     algorithm because we terminate the recursion when the input bezier has
+ *     been subdivided to 'flat' instead of when the offset is good.
+ *     This is valuable because it's possible to create an input curve that does not have a
+ *     good offset approximation down to the resolution of a double.
+ */
+void
+curve_offset (knots_t self, double dist, double tolerance, void (*curve_fn)(void *, knots_t), void *closure)
+{
+    cairo_bool_t recurse;
+    double break_point;
+    double max_error = tolerance;
+    knots_t offset = knot_offset(self, dist, &recurse);
+
+    if (!recurse) {
+	/* we need to make sure we have a finite offset before checking the error */
+	ensure_finite(offset);
+	recurse = !error_within_bounds_elber(self, offset, dist, max_error, &break_point);
+    } else {
+	/* TODO: we could probably choose a better place to split than halfway
+	 * for times when we know we're going to recurse */
+	break_point = .5;
+    }
+
+    if (recurse) {
+	/* error is too large. subdivide on max t and try again. */
+	knots_t left = self;
+	knots_t right;
+
+	cairo_bool_t is_continuous;
+
+	/* We need to do something to handle regions of high curvature
+	 *
+	 * skia uses the first and second derivitives at the points of 'max curvature' to check for
+	 * pinchynes.
+	 * qt tests whether the points are close and reverse direction.
+	 *   float l = (orig->x1 - orig->x2)*(orig->x1 - orig->x2) + (orig->y3 - orig->y4)*(orig->y3 - orig->y4);
+	 * float dot = (orig->x1 - orig->x2)*(orig->x3 - orig->x4) + (orig->y1 - orig->y2)*(orig->y3 - orig->y4);
+	 *
+	 * we just check if the knot is flat and has large error.
+	 */
+	if (is_knot_flat(self)) {
+	    /* we don't want to recurse anymore because we're pretty flat,
+	     * instead approximate the offset curve with an arc.
+	     *
+	     * The previously generated offset curve can become really bad if the the intersections
+	     * are far away from the original curve (self) */
+
+	    approximate_offset_curve_with_shin (self, dist, curve_fn, closure);
+	    return;
+	}
+
+	is_continuous = _de_casteljau_t(&left, &right, break_point);
+
+	curve_offset(left, dist, tolerance, curve_fn, closure);
+
+	if (!is_continuous) {
+	    /* add circle:
+	       we can use the exit the normal from the left side
+	       as the begining normal of our cusp */
+	    vector_t normal = perp (end_tangent (left));
+	    semi_circle (left.d.x, left.d.y, dist, normal, curve_fn, closure);
+	}
+
+	curve_offset(right, dist, tolerance, curve_fn, closure);
+    } else {
+	curve_fn(closure, offset);
+    }
+}
+
diff --git a/src/cairo-spline-offset.h b/src/cairo-spline-offset.h
new file mode 100644
index 000000000..9faa2cb49
--- /dev/null
+++ b/src/cairo-spline-offset.h
@@ -0,0 +1,18 @@
+typedef cairo_point_double_t point_t;
+typedef struct _knots {
+	    point_t a,b,c,d;
+} knots_t;
+
+void print_knot(knots_t k);
+
+struct curve_list {
+	knots_t curve;
+	double dist;
+	struct curve_list *next;
+	struct curve_list *prev;
+};
+
+typedef struct curve_list curve_list_t;
+
+void
+curve_offset(knots_t self, double offset, double tolerance, void (*curve_to_fn)(void *path, knots_t k), void *closure);
diff --git a/src/cairo-stroke-to-path.c b/src/cairo-stroke-to-path.c
new file mode 100644
index 000000000..8c0da821a
--- /dev/null
+++ b/src/cairo-stroke-to-path.c
@@ -0,0 +1,1508 @@
+/* -*- Mode: c; tab-width: 8; c-basic-offset: 4; indent-tabs-mode: t; -*- */
+/* cairo - a vector graphics library with display and print output
+ *
+ * Copyright © 2009 Mozilla Foundation
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (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.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ * The Original Code is the cairo graphics library.
+ *
+ * The Initial Developer of the Original Code is Jeff Muizelaar
+ *
+ */
+
+/* Takes a path as input and produces a stroked path as output.
+ * The resulting path may self intersect so it needs to be filled with
+ * a non-zero winding rule. */
+
+#define _BSD_SOURCE /* for hypot() */
+#include "cairoint.h"
+#include "cairo-path-fixed-private.h"
+#include "cairo-spline-offset.h"
+
+#define printf(a, ...)
+
+typedef struct {
+    cairo_path_fixed_t *path;
+    cairo_matrix_t *ctm;
+
+    /* we need to keep track of whether we've drawn anything yet
+       so that we can do a move_to to the intial point of a curve
+       before drawing it */
+    cairo_bool_t has_current_point;
+    cairo_int_status_t status;
+} path_output_t;
+
+static void line_to (path_output_t *path, double x, double y)
+{
+    if (path->ctm)
+	cairo_matrix_transform_point (path->ctm, &x, &y);
+    path->has_current_point = TRUE;
+
+    if (!path->status)
+	path->status = _cairo_path_fixed_line_to (path->path,
+	                                          _cairo_fixed_from_double(x),
+	                                          _cairo_fixed_from_double(y));
+}
+
+static void curve_to (path_output_t *path, double x0, double y0, double x1, double y1, double x2, double y2)
+{
+    if (path->ctm) {
+	cairo_matrix_transform_point (path->ctm, &x0, &y0);
+	cairo_matrix_transform_point (path->ctm, &x1, &y1);
+	cairo_matrix_transform_point (path->ctm, &x2, &y2);
+    }
+
+    if (!path->status)
+	path->status = _cairo_path_fixed_curve_to (path->path,
+	    _cairo_fixed_from_double (x0), _cairo_fixed_from_double (y0),
+	    _cairo_fixed_from_double (x1), _cairo_fixed_from_double (y1),
+	    _cairo_fixed_from_double (x2), _cairo_fixed_from_double (y2));
+
+}
+
+static void close_path (path_output_t *path)
+{
+    if (!path->status)
+	path->status = _cairo_path_fixed_close_path (path->path);
+    if (!path->status)
+	_cairo_path_fixed_new_sub_path (path->path);
+}
+
+typedef point_t vector_t;
+
+static inline double
+dot (vector_t a, vector_t b)
+{
+    return a.x * b.x + a.y * b.y;
+}
+
+static cairo_bool_t
+point_eq (point_t a, point_t b)
+{
+    return a.x == b.x && a.y == b.y;
+}
+
+static inline vector_t
+perp (vector_t v)
+{
+    vector_t p = {-v.y, v.x};
+    return p;
+}
+
+static inline vector_t
+flip (vector_t a)
+{
+    vector_t ar;
+    ar.x = -a.x;
+    ar.y = -a.y;
+    return ar;
+}
+
+/* given a normal rotate the vector 90 degrees to the right clockwise
+ * This function has a period of 4. e.g. swap(swap(swap(swap(x) == x */
+static vector_t
+swap (vector_t a)
+{
+    vector_t ar;
+    /* one of these needs to be negative. We choose a.x so that we rotate to the right instead of negating */
+    ar.x = a.y;
+    ar.y = -a.x;
+    return ar;
+}
+
+static vector_t
+unperp (vector_t a)
+{
+    return swap (a);
+}
+
+/* Compute a spline approximation of the arc
+   centered at xc, yc from the angle a to the angle b
+
+   The angle between a and b should not be more than a
+   quarter circle (pi/2)
+
+   The approximation is similar to an approximation given in:
+   "Approximation of a cubic bezier curve by circular arcs and vice versa"
+   by Alekas Riškus. However that approximation becomes unstable when the
+   angle of the arc approaches 0.
+
+   This approximation is inspired by a discusion with Boris Zbarsky
+   and essentially just computes:
+
+     h = 4.0/3.0 * tan ((angle_B - angle_A) / 4.0);
+
+   without converting to polar coordinates.
+
+   A different way to do this is covered in "Approximation of a cubic bezier
+   curve by circular arcs and vice versa" by Alekas Riškus. However, the method
+   presented there doesn't handle arcs with angles close to 0 because it
+   divides by the perp dot product of the two angle vectors.
+   */
+
+static void
+arc_segment (path_output_t *path,
+                    double   xc,
+                    double   yc,
+                    double   radius,
+		    vector_t a,
+		    vector_t b)
+{
+    double r_sin_A, r_cos_A;
+    double r_sin_B, r_cos_B;
+    double h;
+    vector_t mid, mid2;
+    double l;
+
+    r_sin_A = radius * a.y;
+    r_cos_A = radius * a.x;
+    r_sin_B = radius * b.y;
+    r_cos_B = radius * b.x;
+
+    /* bisect the angle between 'a' and 'b' with 'mid' */
+    mid.x = a.x + b.x;
+    mid.y = a.y + b.y;
+    l = sqrt(mid.x*mid.x + mid.y*mid.y);
+    mid.x /= l;
+    mid.y /= l;
+
+    /* bisect the angle between 'a' and 'mid' with 'mid2' this is parallel to a
+     * line with angle (B - A)/4 */
+    mid2.x = a.x + mid.x;
+    mid2.y = a.y + mid.y;
+
+    h = (4./3.)*dot(perp(a), mid2)/dot(a, mid2);
+    //h = (4./3.)*(-a.y*mid.x + a.x*mid.y)/(a.x*mid2.x + a.y*mid2.y);
+
+    curve_to (path,
+                    xc + r_cos_A - h * r_sin_A,
+                    yc + r_sin_A + h * r_cos_A,
+                    xc + r_cos_B + h * r_sin_B,
+                    yc + r_sin_B - h * r_cos_B,
+                    xc + r_cos_B,
+                    yc + r_sin_B);
+}
+
+
+static vector_t
+bisect (vector_t a, vector_t b)
+{
+    vector_t mid;
+    double len;
+
+    if (dot (a, b) >= 0) {
+	/* if the angle between a and b is accute, then we can
+	 * just add the vectors and normalize */
+	mid.x = a.x + b.x;
+	mid.y = a.y + b.y;
+    } else {
+	/* otherwise, we can flip a, add it
+	 * and then use the perpendicular of the result */
+	a = flip (a);
+	mid.x = a.x + b.x;
+	mid.y = a.y + b.y;
+	mid = perp (mid);
+    }
+
+    /* normalize */
+    /* because we assume that 'a' and 'b' are normalized, we can use
+     * sqrt instead of hypot because the range of mid is limited */
+    len = sqrt (mid.x*mid.x + mid.y*mid.y);
+    mid.x /= len;
+    mid.y /= len;
+    return mid;
+}
+
+static void
+arc (path_output_t *path, double xc, double yc, double radius, vector_t a, vector_t b)
+{
+    /* find a vector that bisects the angle between a and b */
+    vector_t mid_v = bisect (a, b);
+
+    /* construct the arc using two curve segments */
+    arc_segment (path, xc, yc, radius, a, mid_v);
+    arc_segment (path, xc, yc, radius, mid_v, b);
+}
+
+static inline void
+join_round (path_output_t *path, point_t center, vector_t a, vector_t b, double radius)
+{
+    /*
+    int ccw = dot (perp (b), a) >= 0; // XXX: is this always true?
+    yes, otherwise we have an interior angle.
+    assert (ccw);
+    */
+    arc (path, center.x, center.y, radius, a, b);
+}
+
+static void
+draw_circle (path_output_t *path, point_t center, double radius)
+{
+    vector_t a = {1, 0};
+    vector_t b = {-1, 0};
+
+    /* draw a line to the begining of the arc */
+    line_to (path, center.x + radius * a.x, center.y + radius * a.y);
+
+    arc (path, center.x, center.y, radius, a, b);
+    arc (path, center.x, center.y, radius, b, a);
+}
+
+static point_t
+c2p (cairo_matrix_t *ctm_inverse, cairo_point_t c)
+{
+    point_t p;
+    p.x = _cairo_fixed_to_double(c.x);
+    p.y = _cairo_fixed_to_double(c.y);
+    if (ctm_inverse)
+	cairo_matrix_transform_point (ctm_inverse, &p.x, &p.y);
+
+    return p;
+}
+
+#if 1
+static vector_t
+compute_normal (point_t p0, point_t p1)
+{
+    vector_t n;
+    /* u is parallel vector */
+    double ux = p1.x - p0.x;
+    double uy = p1.y - p0.y;
+
+    /* we actually want the inverse square root here. then we won't have to divide */
+    double ulen = hypot(ux, uy);
+    /* cairo uses hypot which does sqrt(ux*ux + uy*uy) without the accuracy problems
+       i.e. hypot is slower, but more accurate */
+    assert (ulen != 0);
+    // v is perpendicular *unit* vector
+    n.x = -uy/ulen;
+    n.y = ux/ulen;
+
+    return n;
+}
+#else
+#include <emmintrin.h>
+/* a quick and dirty sse implementation. This is much lower
+ * precision than the version above and not as fast as it
+ * could be. */
+static vector_t
+compute_normal (point_t p0, point_t p1)
+{
+    vector_t n;
+    /* u is parallel vector */
+    float ux = p1.x - p0.x;
+    float uy = p1.y - p0.y;
+
+
+    __m128 mux = _mm_load_ss (&ux);
+    __m128 muy = _mm_load_ss (&uy);
+
+    __m128 mux2 = _mm_mul_ss(mux, mux);
+    __m128 muy2 = _mm_mul_ss(muy, muy);
+
+    __m128 rlen = _mm_rsqrt_ss(_mm_add_ss(muy2, mux2));
+#if 1
+    mux = _mm_mul_ss(mux, rlen);
+    muy = _mm_sub_ps(_mm_setzero_ps(), _mm_mul_ss(muy, rlen)); //negate
+    __m128 results = _mm_shuffle_ps (mux, muy, _MM_SHUFFLE(0, 0, 0, 0)); // get bottom words
+    results = _mm_shuffle_ps (results, results, _MM_SHUFFLE(0, 0, 0, 2)); // x and y
+#else
+    rlen = _mm_shuffle_ps (rlen, rlen, _MM_SHUFFLE(0, 0, 0, 0)); // get bottom words
+    muy = _mm_sub_ps(_mm_setzero_ps(), muy);
+    __m128 results = _mm_shuffle_ps (mux, muy, _MM_SHUFFLE(0, 0, 0, 0)); // get bottom words
+    results = _mm_shuffle_ps (results, results, _MM_SHUFFLE(0, 0, 0, 2)); // x and y
+    results = _mm_mul_ps(mux, rlen);
+#endif
+
+
+    __m128d result = _mm_cvtps_pd (results);
+    _mm_storeu_pd(&n.x, result);
+#if 0
+    __m128 muxy = _mm_shuffle_ps (mux, muy, _MM_SHUFFLE(0, 0, 0, 0)); // get bottom words
+    /* we actually want the inverse square root here. then we won't have to divide */
+    muxy = _mm_mul_ps(muxy, muxy);
+    muxy = _mm_shuffle_ps (mux, muy, _MM_SHUFFLE(0, 0, 0, 0)); // get bottom words
+#endif
+    return n;
+}
+#endif
+
+
+
+/* computes the outgoing normal for curve given by a,b,c,d
+ * returning the original normal if the curve is point */
+static vector_t
+curve_normal (vector_t original_normal, point_t a, point_t b, point_t c, point_t d)
+{
+    vector_t n = original_normal;
+    if (point_eq (c, d)) {
+	if (point_eq (b, c)) {
+	    if (!point_eq (a, b)) {
+		n = compute_normal (a, b);
+	    }
+	} else {
+	    n = compute_normal (b, c);
+	}
+    } else {
+	n = compute_normal (c, d);
+    }
+    return n;
+}
+static void ensure_finite (knots_t k)
+{
+
+    if (!(
+	    isfinite(k.a.x) &&
+	    isfinite(k.a.y) &&
+	    isfinite(k.b.x) &&
+	    isfinite(k.b.y) &&
+	    isfinite(k.c.x) &&
+	    isfinite(k.c.y) &&
+	    isfinite(k.d.x) &&
+	    isfinite(k.d.y))) {
+	print_knot(k);
+	assert(0);
+    }
+}
+
+
+static void offset_curve_to(void *closure, knots_t c)
+{
+    path_output_t *path = closure;
+    ensure_finite(c);
+    //XXX: it would be nice if we didn't have to check this for every point
+    if (!path->has_current_point)
+	line_to(path, c.a.x, c.a.y);
+    curve_to(path, c.b.x, c.b.y, c.c.x, c.c.y, c.d.x, c.d.y);
+}
+
+static void
+draw_offset_curve (path_output_t *path,
+	point_t a, point_t b, point_t c, point_t d, double offset)
+{
+    knots_t curve = {a, b, c, d};
+    /* we want greater tolerance when the lines we are stroking are narrower
+     * because any deviations will be more noticeable.
+     * Fortunately, the narrow the lines are the easier it is to approximate
+     * the offset curve. It could be argued that this is a bad choice because
+     * one could draw two large but slightly different stroke widths on top of
+     * each other and the result wouldn't necessarily match. However,
+     * I'm not sure if tha's actually a problem. */
+    double tolerance = offset * 2 / 500.; /* the choice of 500 is pretty arbitrary */
+
+    curve_offset(curve, offset, tolerance, offset_curve_to, path);
+}
+
+/* Finds the intersection of two lines each defined by a point and a normal.
+   From "Example 2: Find the interesection of two lines" of
+   "The Pleasures of "Perp Dot" Products"
+   F. S. Hill, Jr. */
+static point_t
+line_intersection (point_t A, vector_t a_perp, point_t B, vector_t b_perp)
+{
+    point_t intersection;
+    double denom;
+    double t;
+
+    vector_t a = unperp(a_perp);
+    vector_t c = {B.x - A.x, B.y - A.y};
+    denom = dot(b_perp, a);
+    if (denom == 0.0) {
+	assert(0 && "TROUBLE");
+    }
+
+    t = dot (b_perp, c)/denom;
+
+    intersection.x = A.x+t*(a.x);
+    intersection.y = A.y+t*(a.y);
+
+    return intersection;
+}
+
+static cairo_bool_t
+is_interior_angle (vector_t a, vector_t b) {
+    /* angles of 180 and 0 degress will evaluate to 0, however
+     * we to treat 180 as an interior angle and 180 as an exterior angle */
+    return dot (perp (a), b) > 0 ||
+	(a.x == b.x && a.y == b.y); /* 0 degrees is interior */
+}
+
+static void
+join_segment_line (path_output_t *path, cairo_stroke_style_t *style, point_t p, vector_t s1_normal, vector_t s2_normal)
+{
+    double miter_limit = style->miter_limit;
+    cairo_line_join_t join = style->line_join;
+    double offset = style->line_width/2;
+
+    point_t start = {p.x + s1_normal.x*offset, p.y + s1_normal.y*offset};
+    point_t end   = {p.x + s2_normal.x*offset, p.y + s2_normal.y*offset};
+
+    if (is_interior_angle (s1_normal, s2_normal)) {
+	line_to(path, start.x, start.y);
+	//XXX: while you wouldn't think this point is necessary, it is when we the segments we are joining
+	//are shorter than the line width
+	/* qt does a check to avoid adding this additional point when possible */
+	line_to(path, p.x, p.y);
+
+	line_to(path, end.x, end.y);
+    } else {
+	switch (join) {
+	    case CAIRO_LINE_JOIN_ROUND:
+		line_to(path, start.x, start.y);
+		join_round(path, p, s1_normal, s2_normal, offset);
+		break;
+	    case CAIRO_LINE_JOIN_MITER:
+		{
+		    double in_dot_out = -s1_normal.x*s2_normal.x + -s1_normal.y*s2_normal.y;
+		    // XXX: we are going to have an extra colinear segment here */
+		    if (2 <= miter_limit*miter_limit * (1 - in_dot_out)) {
+			point_t intersection = line_intersection(start, s1_normal, end, s2_normal);
+			line_to(path, intersection.x, intersection.y);
+			break;
+		    }
+		}
+		// fall through
+	    case CAIRO_LINE_JOIN_BEVEL:
+		line_to(path, start.x, start.y);
+		line_to(path, end.x, end.y);
+		break;
+	}
+    }
+}
+
+
+static void
+join_segment_curve (path_output_t *path, cairo_stroke_style_t *style, point_t p, vector_t s1_normal, vector_t s2_normal)
+{
+    cairo_line_join_t join = style->line_join;
+    double offset = style->line_width/2;
+    double miter_limit = style->miter_limit;
+
+    point_t start = {p.x + s1_normal.x*offset, p.y + s1_normal.y*offset};
+    point_t end   = {p.x + s2_normal.x*offset, p.y + s2_normal.y*offset};
+    if (is_interior_angle (s1_normal, s2_normal)) {
+	//XXX: while you wouldn't think this point is necessary, it is when we the segments we are joining
+	//are shorter than the line width
+	/* qt does a check to avoid adding this additional point when possible */
+	line_to(path, p.x, p.y);
+	line_to(path, end.x, end.y);
+    } else {
+	switch (join) {
+	    case CAIRO_LINE_JOIN_ROUND:
+		join_round(path, p, s1_normal, s2_normal, offset);
+		break;
+	    case CAIRO_LINE_JOIN_MITER:
+	    {
+		double in_dot_out = -s1_normal.x*s2_normal.x + -s1_normal.y*s2_normal.y;
+		// XXX: we are going to have an extra colinear segment here */
+		if (2 <= miter_limit*miter_limit * (1 - in_dot_out)) {
+		    point_t intersection = line_intersection(start, s1_normal, end, s2_normal);
+		    line_to(path, intersection.x, intersection.y);
+		    break;
+		}
+	    }
+		// fall through
+	    case CAIRO_LINE_JOIN_BEVEL:
+		line_to(path, end.x, end.y);
+		break;
+	}
+    }
+}
+
+/* we have some different options as to what the semantics of cap_line should be:
+   - here's how skia's works:
+     - it uses a virtual function call instead of a switch statement.
+     - it assumes that the caller has added the start point of the cap
+     - here's an interesting bit: the square capper has different behaviour
+       depending on whether the last segment is a curve or not.
+       If it is a line it will adjust the last point of the path
+       and only set the outside point of the cap. This avoids the
+       the extra points when capping a square line.
+  - qt and mesa's openvg state tracker do this cute thing
+    where they treat caps as joins. i.e a butt cap maps
+    to a miter join, a round cap to a round join and
+    a square cap is special cased.
+    3 line segments are always emitted for a square cap
+    it seems the code assumes that the first point is implied?
+
+    I don't think the similarities of caps and joins outweigh the differences. Further
+    I think keeping them separate makes it easier to understand.
+*/
+static void
+cap_line (path_output_t *path, cairo_stroke_style_t *style, point_t end, vector_t normal)
+{
+    cairo_line_cap_t cap = style->line_cap;
+    double offset = style->line_width/2;
+
+    /* This function will add additional unnecessary lines when it is called after a curve
+     * segment. I'm not sure it's worth trying to eliminate them */
+    switch (cap) {
+	case CAIRO_LINE_CAP_ROUND:
+	    {
+		line_to (path, end.x + offset*normal.x, end.y + offset*normal.y);
+		arc(path, end.x, end.y, offset, normal, flip(normal));
+		break;
+	    }
+	case CAIRO_LINE_CAP_SQUARE:
+	    {
+		// parallel vector
+		double vx = normal.y;
+		double vy = -normal.x;
+		// move the end point down the line by offset
+		end.x += vx*offset;
+		end.y += vy*offset;
+
+		// offset the end point of last_line to draw the cap at the end
+		line_to(path, end.x + offset*normal.x, end.y + offset*normal.y);
+		line_to(path, end.x - offset*normal.x, end.y - offset*normal.y);
+		break;
+	    }
+	case CAIRO_LINE_CAP_BUTT:
+	    {
+		line_to (path, end.x + offset*normal.x, end.y + offset*normal.y);
+		// offset the end point of last_line to draw the cap at the end
+		line_to(path, end.x - offset*normal.x, end.y - offset*normal.y);
+		break;
+	    }
+    }
+}
+
+#define cairo_path_head(path__) (&(path__)->buf.base)
+#define cairo_path_tail(path__) cairo_path_buf_prev (cairo_path_head (path__))
+
+#define cairo_path_buf_next(pos__) \
+    cairo_list_entry ((pos__)->link.next, cairo_path_buf_t, link)
+#define cairo_path_buf_prev(pos__) \
+    cairo_list_entry ((pos__)->link.prev, cairo_path_buf_t, link)
+
+// index is the first segment we need to deal with
+typedef struct {
+    const cairo_path_buf_t *buf;
+    cairo_point_t *points;
+    unsigned int op_index;
+    const cairo_path_buf_t *first;
+} path_ittr;
+
+static const uint8_t num_args[] = {
+	1, /* cairo_path_move_to */
+	1, /* cairo_path_op_line_to */
+	3, /* cairo_path_op_curve_to */
+	0, /* cairo_path_op_close_path */
+    };
+
+static inline path_ittr
+ittr_next (path_ittr index)
+{
+    index.points += num_args[(int)index.buf->op[index.op_index]];
+    index.op_index++;
+    if (index.op_index == index.buf->num_ops) {
+	index.buf = cairo_path_buf_next (index.buf);
+#ifdef ITTR_DEBUG
+	if (index.buf == index.first) {
+	    assert(0);
+	    index.buf = NULL;
+	    index.points = NULL;
+	    index.op_index = 0xbaadf00d;
+	}
+#endif
+	index.op_index = 0;
+	index.points = index.buf->points;
+    }
+    return index;
+}
+
+static inline path_ittr
+ittr_prev (path_ittr index)
+{
+    if (index.op_index == 0) {
+	index.buf = cairo_path_buf_prev (index.buf);
+	index.op_index = index.buf->num_ops-1;
+	index.points = index.buf->points + (index.buf->num_points);
+    } else {
+	index.op_index--;
+    }
+    index.points -= num_args[(int)index.buf->op[index.op_index]];
+    return index;
+}
+
+static inline cairo_path_op_t
+ittr_op (path_ittr i)
+{
+    return i.buf->op[i.op_index];
+}
+
+static cairo_bool_t
+ittr_done(path_ittr i)
+{
+    return i.buf == NULL;
+}
+
+static cairo_bool_t
+ittr_last (path_ittr i)
+{
+    return i.op_index+1 == i.buf->num_ops && cairo_path_buf_next(i.buf) == i.first;
+}
+
+static cairo_bool_t
+ittr_eq (path_ittr a, path_ittr b)
+{
+    return a.buf == b.buf && a.op_index == b.op_index;
+}
+
+
+static void
+draw_degenerate_path (path_output_t *path_out, cairo_stroke_style_t *style, point_t point)
+{
+    if (style->line_cap == CAIRO_LINE_CAP_ROUND) {
+	draw_circle(path_out, point, style->line_width/2);
+    }
+}
+
+/* it might be interesting to see if we can just skip degenerate segments... */
+
+static cairo_status_t
+_cairo_subpath_stroke_to_path (path_ittr *ittr,
+			    cairo_stroke_style_t *stroke_style,
+			    path_output_t *path_out,
+			    cairo_matrix_t *ctm_inverse)
+{
+    double offset = stroke_style->line_width/2;
+    cairo_status_t status = CAIRO_STATUS_SUCCESS;
+    path_ittr i = *ittr;
+    cairo_matrix_t *ctm_i = ctm_inverse;
+
+    point_t begin = c2p (ctm_i, i.points[0]);
+    point_t end_point = begin;
+
+    cairo_path_op_t op;
+    cairo_point_t *points;
+    vector_t end_normal;
+    vector_t initial_normal;
+    point_t end;
+
+    int op_count = 0;
+    cairo_bool_t closed_path = FALSE;
+    point_t last_point;
+    vector_t closed_normal;
+
+    assert(ittr_op(i) == CAIRO_PATH_OP_MOVE_TO);
+
+    /* handle beginging of the line: we want to get an intial normal to work with
+     * before we start drawing */
+    while (1) {
+	i = ittr_next (i);
+	op = ittr_op (i);
+	points = i.points;
+
+	if (op == CAIRO_PATH_OP_MOVE_TO) {
+	    *ittr = i;
+	    draw_degenerate_path(path_out, stroke_style, begin);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+
+	if (op == CAIRO_PATH_OP_CLOSE_PATH) {
+	    *ittr = ittr_next (i);
+	    draw_degenerate_path(path_out, stroke_style, begin);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+
+	if (op == CAIRO_PATH_OP_LINE_TO) {
+	    point_t p = c2p(ctm_i, points[0]);
+	    //XXX: it would be better if we compared points in fixed point...
+	    if (!point_eq(begin, p)) {
+		initial_normal = compute_normal(begin, p);
+		break;
+	    }
+	} else {
+	    point_t p0 = c2p(ctm_i, points[0]);
+	    point_t p1 = c2p(ctm_i, points[1]);
+	    point_t p2 = c2p(ctm_i, points[2]);
+	    assert(ittr_op (i) == CAIRO_PATH_OP_CURVE_TO);
+	    //XXX: maybe use the helper function?
+	    if (!point_eq (begin, p0)) {
+		initial_normal = compute_normal (begin, p0);
+		break;
+	    } else {
+		if (!point_eq (begin, p1)) {
+			initial_normal = compute_normal (begin, p1);
+			break;
+		} else {
+		    if (!point_eq (begin, p2)) {
+			initial_normal = compute_normal (begin, p2);
+			break;
+		    }
+		}
+	    }
+	}
+
+	if (ittr_last (i)) {
+	    /* if we had a normal for the last op we would have broken out already */
+	    printf("degen\n");
+	    *ittr = i;
+	    draw_degenerate_path (path_out, stroke_style, begin);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+    }
+
+    end_normal = initial_normal;
+
+    /* end_normal needs to be the normal between i.points[] and the preceeding point */
+    while (!ittr_last (i)) {
+	    op_count++;
+
+	    i = ittr_next (i);
+	    cairo_path_op_t next_op = ittr_op (i);
+	    cairo_point_t *next_points = i.points;
+
+	    switch (next_op) {
+	    case CAIRO_PATH_OP_MOVE_TO:
+		*ittr = i;
+		i = ittr_prev (i);
+		op_count--;
+		// we are done with the current subpath
+		// cap it and return to the begining
+		//assert(0);
+		goto done;
+	    case CAIRO_PATH_OP_LINE_TO:
+	    case CAIRO_PATH_OP_CURVE_TO:
+		// get first point of next_op
+		end = c2p (ctm_i, next_points[0]);
+		break;
+	    case CAIRO_PATH_OP_CLOSE_PATH:
+		// close the path
+		// XXX: relies on the fact that CLOSE_PATH is never the last path_op (it is always followed by an OP_MOVE)
+		closed_path = TRUE;
+		// the first point of the first op
+		end = c2p (ctm_i, ittr->points[0]);
+		break;
+	    default:
+		ASSERT_NOT_REACHED;
+	    }
+	    if (unlikely (status))
+		return status;
+
+	    switch (op) {
+	    case CAIRO_PATH_OP_CURVE_TO:
+		{
+		    point_t p0 = c2p (ctm_i, points[0]);
+		    point_t p1 = c2p (ctm_i, points[1]);
+		    point_t p2 = c2p (ctm_i, points[2]);
+
+		    vector_t n2 = end_normal;
+		    if (!point_eq (p0, p1))
+			n2 = compute_normal (p0, p1);
+		    if (!point_eq (p1, p2))
+			n2 = compute_normal (p1, p2);
+		    end_normal = n2;
+
+		    draw_offset_curve (path_out, begin,
+			    p0,
+			    p1,
+			    p2,
+			    offset);
+
+		    if (!point_eq (p2, end)) {
+			vector_t n3 = compute_normal (p2, end);
+			join_segment_curve (path_out, stroke_style, p2, n2, n3);
+			end_normal = n3;
+		    }
+
+		    begin = p2;
+		    end_point = begin;
+		}
+		break;
+	    case CAIRO_PATH_OP_LINE_TO:
+		{
+		    point_t mid = c2p (ctm_i, points[0]);
+		    if (!point_eq (mid, end)) {
+			vector_t n1 = end_normal;
+			vector_t n2 = compute_normal (mid, end);
+			join_segment_line (path_out, stroke_style, mid, n1, n2);
+			end_normal = n2;
+			begin = mid;
+			end_point = mid;
+		    }
+		}
+		break;
+	    default:
+		ASSERT_NOT_REACHED;
+	    }
+	    if (closed_path)
+		break;
+	    op = next_op;
+	    points = next_points;
+    };
+
+    if (!closed_path)
+	*ittr = i;
+done:
+    // op should be set the last op in the path
+    if (closed_path) {
+	point_t first_point = c2p (ctm_i, ittr->points[0]);
+
+	assert(ittr_op (*ittr) == CAIRO_PATH_OP_MOVE_TO);
+	closed_normal = end_normal;
+	last_point = end_point;
+	join_segment_line (path_out, stroke_style, first_point, end_normal, initial_normal);
+
+	// XXX: is this right?
+	begin = first_point;
+	*ittr = ittr_next (i);
+	path_out->has_current_point = FALSE;
+	close_path (path_out);
+    } else {
+	if (op == CAIRO_PATH_OP_CURVE_TO) {
+	    point_t p0 = c2p (ctm_i, points[0]);
+	    point_t p1 = c2p (ctm_i, points[1]);
+	    point_t p2 = c2p (ctm_i, points[2]);
+
+	    /* draw capped curve segment */
+	    draw_offset_curve (path_out, end_point, p0, p1, p2, offset);
+	    cap_line (path_out, stroke_style, p2, curve_normal (end_normal, end_point, p0, p1, p2));
+	    draw_offset_curve (path_out, p2, p1, p0, end_point, offset);
+
+	    begin = p0;
+	} else if (op == CAIRO_PATH_OP_LINE_TO) {
+	    begin = c2p (ctm_i, points[0]);
+	    cap_line (path_out, stroke_style, begin, end_normal);
+	} else {
+	    assert(0);
+	}
+    }
+
+    end_normal = flip (end_normal);
+
+    // pre-conditions
+    i = ittr_prev(i);
+    op = ittr_op(i);
+    points = i.points;
+
+    /* traverse backwards to the begining drawing
+     * the other side of the stroked path */
+    while (op_count) {
+	i = ittr_prev (i);
+	cairo_path_op_t next_op = ittr_op (i);
+	cairo_point_t *next_points = i.points;
+
+	switch (next_op) {
+	    case CAIRO_PATH_OP_MOVE_TO:
+		end = c2p (ctm_i, next_points[0]);
+		break;
+	    case CAIRO_PATH_OP_LINE_TO:
+		end = c2p (ctm_i, next_points[0]);
+		break;
+	    case CAIRO_PATH_OP_CURVE_TO:
+		// get first point of next_op
+		end = c2p (ctm_i, next_points[2]);
+		break;
+	    case CAIRO_PATH_OP_CLOSE_PATH:
+		// close the path
+		assert(0);
+		break;
+	    default:
+		ASSERT_NOT_REACHED;
+	}
+	if (unlikely (status))
+	    return status;
+
+	switch (op) {
+	    case CAIRO_PATH_OP_CURVE_TO:
+		{
+		    point_t p2 = c2p (ctm_i, points[2]);
+		    point_t p1 = c2p (ctm_i, points[1]);
+		    point_t p0 = c2p (ctm_i, points[0]);
+		    if (!point_eq (p2, p1)) {
+			vector_t n0 = compute_normal (p2, p1);
+			// XXX whether we use join_segment_line or join_segment_curve here depends
+			// on whether we are coming from a line or from a curve. Does it matter?
+			// this used to be join_segment_curve but it was wrong for the following:
+			//         cairo_move_to(cr, 550, 400);
+			//	   cairo_curve_to(cr, 531.687500, 469.878906, 415.203125, 506.437500, 247.179688, 524.625000);
+			//         cairo_line_to(cr, 257.945312, 624.046875);
+			join_segment_line (path_out, stroke_style, p2, end_normal, n0);
+			end_normal = n0;
+		    }
+		    end_normal = curve_normal (end_normal, p2, p1, p0, end);
+		    draw_offset_curve (path_out,
+			    p2,
+			    p1,
+			    p0,
+			    end,
+			    offset);
+		    //XXX: this feels suspect...
+		    begin = p0;
+		    end_point = end;
+		}
+		break;
+	    case CAIRO_PATH_OP_LINE_TO:
+		{
+		    point_t mid = c2p (ctm_i, points[0]);
+		    vector_t n1 = end_normal;
+		    if (!point_eq (mid, end)) {
+			vector_t n2 = compute_normal (mid, end);
+			join_segment_line (path_out, stroke_style, mid, n1, n2);
+			end_normal = n2;
+			begin = mid;
+			end_point = end;
+		    }
+		}
+		break;
+	    default:
+		ASSERT_NOT_REACHED;
+	}
+
+	points = next_points;
+	op = next_op;
+	op_count--;
+
+    }
+
+    if (closed_path) {
+	vector_t n2 = end_normal;
+	printf("finish closed\n");
+	assert(ittr_op (i) == CAIRO_PATH_OP_MOVE_TO);
+	if (!point_eq (end_point, last_point)) {
+	    n2 = compute_normal (end_point, last_point);
+	    join_segment_line (path_out, stroke_style, c2p (ctm_i, i.points[0]), end_normal, n2);
+	}
+
+	join_segment_line (path_out, stroke_style, last_point, n2, flip (closed_normal));
+    } else {
+	cap_line (path_out, stroke_style, end_point, end_normal);
+    }
+    close_path (path_out);
+
+    return CAIRO_STATUS_SUCCESS;
+}
+
+typedef struct {
+        point_t point;
+        vector_t normal;
+} dash_point_t;
+
+static double
+distance (point_t p0, point_t p1)
+{
+    double ux = p1.x - p0.x;
+    double uy = p1.y - p0.y;
+    return hypot(ux, uy);
+}
+
+/* i points to the last op we've just finished */
+static dash_point_t
+skip_subpath (dash_point_t start_point, path_output_t *path_out, path_ittr *i, path_ittr subpath_start, path_ittr *subpath_end, cairo_matrix_t *ctm_i, cairo_bool_t *closed, double length)
+{
+    //XXX: review the naming of points to make sure it makes sense and is consistent
+    path_ittr index = *i;
+    path_ittr next_index;
+    dash_point_t end_point = start_point;
+    point_t next_point;
+    cairo_bool_t done = FALSE;
+
+    assert(!ittr_last(index));
+
+    next_index = ittr_next(index);
+
+    printf("skip_subpath: %d\n", index.op_index);
+
+    if (ittr_op(next_index) == CAIRO_PATH_OP_CLOSE_PATH) {
+	*closed = TRUE;
+	next_index = subpath_start;
+    }
+
+    assert(ittr_op(next_index) == CAIRO_PATH_OP_LINE_TO || ittr_eq(next_index, subpath_start));
+    next_point = c2p(ctm_i, next_index.points[0]);
+
+    while (length > 0) { //XXX: is this condition even needed?
+	double d = distance(end_point.point, next_point);
+	assert(ittr_op(next_index) == CAIRO_PATH_OP_LINE_TO || ittr_op(next_index) == CAIRO_PATH_OP_MOVE_TO);
+	printf("d(%d): %f of %f\n", index.op_index, d, length);
+	if (d >= length) {
+	    /* move the end_point a distance 'length' along the line with the normal end_point.normal */
+	    end_point.point.x += -end_point.normal.y*-length;
+	    end_point.point.y += end_point.normal.x*-length;
+	    length = 0;
+	    printf("skippath break\n");
+	    break;
+	} else {
+	    end_point.point = next_point;
+
+	    if (ittr_eq(next_index, subpath_start)) {
+	        /* we've passed the begining; we are done */
+		done = TRUE;
+		break;
+	    }
+
+	    if (ittr_last(next_index)) {
+		/* we're at the end of the path; we are done */
+		*subpath_end = next_index;
+		done = TRUE;
+		break;
+	    }
+
+	    index = next_index;
+	    next_index = ittr_next(next_index);
+
+	    if (ittr_op(next_index) == CAIRO_PATH_OP_CLOSE_PATH) {
+		/* we have a closed path so wrap around toward the begining (subpath_start) */
+		*closed = TRUE;
+		*subpath_end = ittr_next(next_index);
+		next_index = subpath_start;
+	    } else if (ittr_op(next_index) == CAIRO_PATH_OP_MOVE_TO) {
+		*subpath_end = next_index;
+		done = TRUE;
+		break;
+	    }
+
+	    next_point  = c2p(ctm_i, next_index.points[0]);
+	    if (!point_eq(end_point.point, next_point)) {
+		/* XXX: we might be able to avoid computing the normal for every point we vist
+		 * and instead only compute it for the last segment */
+		end_point.normal = compute_normal(end_point.point, next_point);
+	    }
+
+	    printf("next %f %f\n", d, length);
+	    length -= d;
+	}
+    }
+
+    *i = index;
+
+    if (done)
+	i->buf = NULL;
+
+    return end_point;
+}
+
+static dash_point_t
+dash_subpath (dash_point_t start_point, path_output_t *path_out, path_ittr *index, path_ittr subpath_start, path_ittr *subpath_done, cairo_matrix_t *ctm_i, cairo_stroke_style_t *style, double length, int begin_on, cairo_bool_t *closed, double start_length)
+{
+    dash_point_t end_point;
+    point_t first = start_point.point;
+    path_ittr i = *index;
+    path_ittr subpath_end = {};
+    int closed_path = 0;
+    int done = 0;
+    vector_t n1 = start_point.normal;
+    path_ittr start_index = i;
+    path_ittr next_i;
+    point_t p1;
+    // investigate more normal reuse.
+    // ideally we would only compute the normal for each segment once
+    // deal with short paths (count < 3)
+    assert(!ittr_last(i));
+    printf("dash_subpath(%f): %d: %f %f\n", length, i.op_index, first.x, first.y);
+
+    next_i = ittr_next(i);
+
+    if (ittr_op(next_i) == CAIRO_PATH_OP_CLOSE_PATH) {
+	subpath_end = i;
+	*subpath_done = ittr_next(next_i);
+	next_i = subpath_start;
+	*closed = TRUE;
+    }
+
+    assert(ittr_op(next_i) == CAIRO_PATH_OP_LINE_TO || ittr_eq(next_i, subpath_start));
+    //printf("off: %f %f, %f %f\n", end_point.point.x, end_point.point.y, path[index].x, path[index].y);
+    p1 = c2p(ctm_i, next_i.points[0]);
+    do {
+	double d;
+	vector_t n2;
+	point_t p2;
+	assert(ittr_op(next_i) == CAIRO_PATH_OP_LINE_TO || ittr_eq(next_i, subpath_start));
+	d = distance (first, p1);
+	if (d >= length) {
+	    end_point.point.x = first.x + length*n1.y;
+	    end_point.point.y = first.y + length*-n1.x;
+	    end_point.normal = n1;
+	    break;
+	}
+
+	if (ittr_last (next_i)) {
+	    /* we're at the end so we must be done */
+	    done = 1;
+	    *subpath_done = next_i;
+	    end_point.point = c2p(ctm_i, next_i.points[0]);
+	    end_point.normal = n1;
+	    break;
+	}
+
+	/* We don't have enough length in the current segment so get the next one
+	 * This means advancing i and next_i */
+
+	if (ittr_eq(next_i, subpath_start)) {
+	    /* we're already at the beginning of the path so we don't have anything left to dash */
+	    printf("done\n");
+	    //XXX done the path
+	    done = 1;
+	    if (!begin_on) {
+		/* we didn't skip the initial segment; so we're all done */
+		length = 0;
+		end_point.point = c2p(ctm_i, next_i.points[0]);
+		end_point.normal = n1;
+		break;
+	    } else {
+		/* the initial segment started on, so join with it and draw it */
+		length = start_length;
+		i = next_i;
+		next_i = ittr_next(next_i);
+	    }
+	} else {
+	    i = next_i;
+	    next_i = ittr_next(next_i);
+
+	    if (ittr_op(next_i) == CAIRO_PATH_OP_CLOSE_PATH) {
+		subpath_end = i;
+		*subpath_done = ittr_next(next_i);
+		next_i = subpath_start;
+		*closed = TRUE;
+	    } else if (ittr_op(next_i) == CAIRO_PATH_OP_MOVE_TO) {
+		done = 1;
+		end_point.point = c2p(ctm_i, i.points[0]);
+		end_point.normal = n1;
+		*subpath_done = next_i;
+		//XXX: it would be good if we could avoid doing this ittr_prev()
+		i = ittr_prev(i);
+		printf("dash: segment end break\n");
+		break;
+	    }
+	    length -= d;
+	}
+
+	p2 = c2p(ctm_i, next_i.points[0]);
+	if (!point_eq(p1, p2)) {
+	    n2 = compute_normal(p1, p2);
+	    end_point.normal = n2;
+	    end_point.point = p2;
+	}
+
+	if (ittr_eq(next_i, start_index)) {
+	    // we've reached the place where we started so we need to loop the path instead of
+	    // capping it
+	    done = 1;
+	    closed_path = 1;
+	    /*
+	    gcc complains that end_point can be used uninitialzed in this function because we don't set
+	    it here when we break out of the loop. I don't think we can get to here without setting end_point
+	    and here's why:
+
+	    p2 == start_index point
+	    p1 must have the same value as p2
+
+	    the only time that next_i can equal start_index is when start_index == subpath_start
+	    this means that the whole path needs to be here. However, unless the path is degenerate
+	    end_point will be set. But no degenerate paths will pass into here because they will be
+	    found else where.
+
+	    Can we make this more explict?
+	    */
+	    printf("closed path\n");
+	    // XXX: do we need to set end_point here? normally we would get an end_point from the code above.
+	    // but if we have a degenerate path we won't
+	    break;
+	}
+
+	// XXX: having to check this twice sort of sucks
+	if (!point_eq(p1, p2)) {
+	    join_segment_line(path_out, style, p1, n1, n2);
+	    n1 = n2;
+	    first = p1;
+	    p1 = p2;
+	}
+    } while (1);
+
+    /* at this point end_point must be initialized and 'i' must point to the point before the end_point */
+
+    *index = i;
+
+    point_t closed_p1;
+    point_t closed_p2;
+    vector_t closed_n1;
+    vector_t closed_n2;
+    vector_t closed_n3;
+    if (closed_path) {
+	    // join the line back with itself
+	    closed_p1 = c2p (ctm_i, i.points[0]);
+	    closed_p2 = end_point.point;
+	    closed_n1 = n1;
+	    closed_n2 = compute_normal (closed_p1, closed_p2);
+	    closed_n3 = compute_normal (closed_p2, c2p(ctm_i, ittr_next (next_i).points[0]));
+	    join_segment_line (path_out, style, closed_p1, closed_n1, closed_n2);
+	    join_segment_line (path_out, style, closed_p2, closed_n2, closed_n3);
+	    close_path (path_out);
+    } else {
+	// draw cap
+	cap_line (path_out, style, end_point.point, end_point.normal);
+    }
+
+    start_point.normal = flip (start_point.normal);
+
+    /* draw the other side of the path, returning to the begining */
+    n1 = flip (end_point.normal);
+    p1 = c2p (ctm_i, i.points[0]);
+    path_ittr prev_i = i;
+    while (!ittr_eq (prev_i, start_index)) {
+	i = prev_i;
+	if (ittr_eq (prev_i, subpath_start)) {
+	    prev_i = subpath_end;
+	} else {
+	    prev_i = ittr_prev (prev_i);
+	}
+	point_t p2 = c2p (ctm_i, prev_i.points[0]);
+	if (!point_eq (p1, p2)) {
+	    vector_t n2 = compute_normal (p1, p2);
+	    join_segment_line (path_out, style, p1, n1, n2);
+
+	    n1 = n2;
+	    p1 = p2;
+	}
+    }
+
+    /* finish up */
+
+    if (closed_path) {
+	join_segment_line (path_out, style, closed_p2, flip (closed_n3), flip (closed_n2));
+	join_segment_line (path_out, style, closed_p1, flip (closed_n2), flip (closed_n1));
+	close_path (path_out);
+    } else {
+	/* draw cap */
+	cap_line(path_out, style, start_point.point, start_point.normal);
+	close_path (path_out);
+    }
+
+    if (done)
+	index->buf = NULL;
+
+    return end_point;
+
+}
+
+static cairo_status_t
+_cairo_subpath_dash_to_path (path_ittr *ittr,
+			    cairo_stroke_style_t *style,
+			    path_output_t *path_out,
+			    cairo_matrix_t *ctm_i)
+{
+
+    cairo_bool_t on = TRUE, begin_on, closed = FALSE;
+
+    double length, start_length;
+
+    double *dash = style->dash;
+    double dash_init = style->dash_offset;
+    int dash_len = style->num_dashes;
+    int dash_state = 0;
+    path_ittr index = *ittr;
+    path_ittr subpath_end = {};
+
+    dash_point_t start_point;
+    dash_point_t end_point;
+    path_ittr subpath_start = index;
+
+    assert(ittr_op(index) == CAIRO_PATH_OP_MOVE_TO);
+
+    ///XXX: this code is different from cairo, figure out which is better
+    /* intiailize dash state */
+    length = dash[dash_state++ % dash_len];
+
+    /* decrease dash_init until we are in the correct dash_state */
+    /* the intervals are closed */
+    /* we make sure that we don't skip over 0 length segments */
+    while (length > 0 && dash_init - length >= 0) {
+	dash_init -= length;
+	on = !on;
+	length = dash[dash_state++ % dash_len];
+    }
+
+    begin_on = on;
+    length -= dash_init;
+    printf("subpath starts %s with length of %f\n", on ? "on" : "off", length);
+
+    start_length = length;
+
+    start_point.point = c2p(ctm_i, index.points[0]);
+    //XXX: to share this with the non dashed code we need to deal with dash on/off
+    while (1) {
+	if (ittr_op(ittr_next(index)) == CAIRO_PATH_OP_LINE_TO) {
+	    point_t p = c2p(ctm_i, ittr_next(index).points[0]);
+	    if (!point_eq(start_point.point, p)) {
+		/* we have a normal so we can start the regular stroking process */
+		start_point.normal = compute_normal(start_point.point, c2p(ctm_i, ittr_next(index).points[0]));
+		break;
+	    }
+	}
+	index = ittr_next(index);
+
+	if (ittr_op(index) == CAIRO_PATH_OP_MOVE_TO) {
+	    *ittr = index;
+	    if (on)
+		draw_degenerate_path (path_out, style, start_point.point);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+
+	if (ittr_op(index) == CAIRO_PATH_OP_CLOSE_PATH) {
+	    *ittr = ittr_next(index);
+	    if (on)
+		draw_degenerate_path (path_out, style, start_point.point);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+
+	if (ittr_last(index)) {
+	    *ittr = index;
+	    if (on)
+		draw_degenerate_path (path_out, style, start_point.point);
+	    return CAIRO_STATUS_SUCCESS;
+	}
+    }
+
+    end_point = start_point;
+
+    /* the first segment is special because we might need to join with it. However, we won't know
+     * what to do until we reach the end of the subpath
+     * Different renderers behave differently here:
+     *   No Join: skia, qt, mesa openvg, coregraphics
+     *   Join: opera, cairo, ghostscript, adobe */
+    if (begin_on) {
+	on = !on;
+	/* if the first dash segment is larger than the entire line we'll skip the whole thing */
+	end_point = skip_subpath (end_point, path_out, &index, subpath_start, &subpath_end, ctm_i, &closed, length);
+	length = dash[dash_state++ % dash_len];
+	printf("done skip: %d\n", index.op_index);
+	if (ittr_done(index)) {
+	    printf("reset\n");
+	    // flip to back on so that we don't skip everything again
+	    on = !on;
+	    index = *ittr;
+	    length = start_length;
+	    end_point = start_point;
+	}
+    }
+
+    do {
+	printf("index: %d\n", index.op_index);
+	if (on) {
+	    end_point = dash_subpath (end_point, path_out, &index, subpath_start, &subpath_end, ctm_i, style, length, begin_on, &closed, start_length);
+	} else {
+	    end_point = skip_subpath (end_point, path_out, &index, subpath_start, &subpath_end, ctm_i, &closed, length);
+	}
+	on = !on;
+	length = dash[dash_state++ % dash_len];
+
+	// cairo does an implicit move_to after a close_path so we need to start a new sub path
+	_cairo_path_fixed_new_sub_path (path_out->path);
+    } while (!ittr_done(index));
+
+    *ittr = index;
+    // draw the initial segment if we haven't yet.
+    if ((!closed && begin_on) || (closed && on)) {
+	printf("draw initial %f\n", start_length);
+	index = subpath_start;
+	end_point = start_point;
+	end_point = dash_subpath (end_point, path_out, &index, subpath_start, &subpath_start, ctm_i, style, start_length, begin_on, &closed, start_length);
+    }
+    *ittr = subpath_end;
+    return CAIRO_STATUS_SUCCESS;
+}
+
+static cairo_status_t
+_cairo_path_dash_to_path (const cairo_path_fixed_t		*path,
+			    cairo_stroke_style_t *stroke_style,
+			     cairo_path_fixed_t         *path_out,
+			     cairo_matrix_t *ctm,
+			     cairo_matrix_t *ctm_inverse)
+{
+    path_ittr i;
+    path_output_t output_path;
+
+    output_path.status = CAIRO_STATUS_SUCCESS;
+    output_path.ctm = ctm;
+    output_path.path = path_out;
+
+    i.buf = i.first = cairo_path_head(path);
+    i.op_index = 0;
+    i.points = i.buf->points;
+
+    while (!ittr_last(i)) {
+	_cairo_subpath_dash_to_path(&i, stroke_style, &output_path, ctm_inverse);
+    }
+
+    return output_path.status;
+}
+
+cairo_status_t
+_cairo_path_stroke_to_path (const cairo_path_fixed_t		*path,
+			    cairo_stroke_style_t *stroke_style,
+			     cairo_path_fixed_t         *path_out,
+			     cairo_matrix_t *ctm,
+			     cairo_matrix_t *ctm_inverse,
+			     double tolerance)
+{
+    path_ittr i;
+    path_output_t output_path;
+
+    cairo_status_t status = CAIRO_STATUS_SUCCESS;
+
+    /* if we have an empty path, we're alredy done */
+    if (cairo_path_head(path)->num_ops == 0)
+	return CAIRO_STATUS_SUCCESS;
+
+    if (stroke_style->num_dashes > 0) {
+	if (path->has_curve_to) {
+	    /* Dashing curved paths directly is tricky because I know of no easy way
+	       to split a bezier curve at a particular length. So, for now, flatten the
+	       path before dashing it */
+	    cairo_path_fixed_t flat_path;
+	    status = _cairo_path_fixed_init_flat_copy (&flat_path,
+		    path,
+		    tolerance);
+
+	    if (unlikely(status)) {
+		_cairo_path_fixed_fini (&flat_path);
+		return status;
+	    }
+
+	    status = _cairo_path_dash_to_path (&flat_path, stroke_style, path_out, ctm, ctm_inverse);
+	    _cairo_path_fixed_fini (&flat_path);
+	    return status;
+	} else {
+	    return _cairo_path_dash_to_path (path, stroke_style, path_out, ctm, ctm_inverse);
+	}
+    }
+
+    output_path.status = CAIRO_STATUS_SUCCESS;
+    output_path.ctm = ctm;
+    output_path.path = path_out;
+
+    i.buf = i.first = cairo_path_head (path);
+    i.op_index = 0;
+    i.points = i.buf->points;
+
+    while (!ittr_last (i)) {
+	output_path.has_current_point = FALSE;
+	_cairo_subpath_stroke_to_path (&i, stroke_style, &output_path, ctm_inverse);
+    }
+
+    return output_path.status;
+}
diff --git a/src/cairo.c b/src/cairo.c
index 157f898f4..8e488dd49 100644
--- a/src/cairo.c
+++ b/src/cairo.c
@@ -4156,6 +4156,15 @@ cairo_append_path (cairo_t		*cr,
 	_cairo_set_error (cr, status);
 }
 
+cairo_path_t *
+cairo_stroke_to_path (cairo_t *cr)
+{
+    if (unlikely (cr->status))
+	return _cairo_path_create_in_error (cr->status);
+
+    return _cairo_stroked_path_create (cr->path, cr->gstate);
+}
+
 /**
  * cairo_status:
  * @cr: a cairo context
diff --git a/src/cairo.h b/src/cairo.h
index 136c5dbf2..33883527d 100644
--- a/src/cairo.h
+++ b/src/cairo.h
@@ -1926,6 +1926,9 @@ cairo_append_path (cairo_t		*cr,
 cairo_public void
 cairo_path_destroy (cairo_path_t *path);
 
+cairo_public cairo_path_t *
+cairo_stroke_to_path (cairo_t *cr);
+
 /* Error status queries */
 
 cairo_public cairo_status_t
diff --git a/src/cairoint.h b/src/cairoint.h
index 346fc51ff..774c06c3c 100644
--- a/src/cairoint.h
+++ b/src/cairoint.h
@@ -1130,6 +1130,11 @@ cairo_private cairo_status_t
 _cairo_path_fixed_init_copy (cairo_path_fixed_t *path,
 			     const cairo_path_fixed_t *other);
 
+cairo_private cairo_status_t
+_cairo_path_fixed_init_flat_copy (cairo_path_fixed_t *path,
+				  const cairo_path_fixed_t *other,
+				  double tolerance);
+
 cairo_private cairo_bool_t
 _cairo_path_fixed_is_equal (const cairo_path_fixed_t *path,
 			    const cairo_path_fixed_t *other);
@@ -1326,6 +1331,15 @@ _cairo_path_fixed_stroke_to_traps (const cairo_path_fixed_t	*path,
 				   double		 tolerance,
 				   cairo_traps_t	*traps);
 
+/* cairo-stroke-to-path.c */
+cairo_status_t
+_cairo_path_stroke_to_path (const cairo_path_fixed_t            *path,
+			    cairo_stroke_style_t *stroke_style,
+			    cairo_path_fixed_t         *path_out,
+			    cairo_matrix_t *ctm,
+			    cairo_matrix_t *ctm_invserse,
+			    double tolerance); // XXX: the matrices could probably be const
+
 cairo_private cairo_status_t
 _cairo_path_fixed_stroke_to_shaper (cairo_path_fixed_t	*path,
 				   const cairo_stroke_style_t	*stroke_style,