New upstream version 5.1.0upstream/5.1.0

12 files changed, 2486 insertions, 0 deletions

diff --git a/app/cornu/CMakeLists.txt b/app/cornu/CMakeLists.txt
new file mode 100644
index 0000000..b54fc80
--- /dev/null
+++ b/app/cornu/CMakeLists.txt
@@ -0,0 +1,27 @@
+PROJECT(cornu)
+
+FILE(GLOB HEADERS *.h)
+
+SET(SOURCES
+	bezctx.c
+	bezctx_xtrkcad.c
+	spiroentrypoints.c
+	spiro.c
+	)
+	
+SET(HEADERS 
+	spiro.h
+	)
+
+
+INCLUDE_DIRECTORIES(${XTrkCAD_BINARY_DIR})
+INCLUDE_DIRECTORIES(${XTrkCAD_SOURCE_DIR})
+INCLUDE_DIRECTORIES(${wlib_SOURCE_DIR}/include)
+INCLUDE_DIRECTORIES(${help_BINARY_DIR})
+
+ADD_LIBRARY(xtrkcad-cornu ${HEADERS} ${SOURCES})
+
+
+
+
+
diff --git a/app/cornu/bezctx.c b/app/cornu/bezctx.c
new file mode 100644
index 0000000..722f5db
--- /dev/null
+++ b/app/cornu/bezctx.c
@@ -0,0 +1,48 @@
+/*
+ppedit - A pattern plate editor for Spiro splines.
+Copyright (C) 2007 Raph Levien
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+#include "bezctx.h"
+
+void bezctx_moveto(bezctx *bc, double x, double y, int is_open)
+{
+    bc->moveto(bc, x, y, is_open);
+}
+
+void bezctx_lineto(bezctx *bc, double x, double y)
+{
+    bc->lineto(bc, x, y);
+}
+
+void bezctx_quadto(bezctx *bc, double x1, double y1, double x2, double y2)
+{
+    bc->quadto(bc, x1, y1, x2, y2);
+}
+
+void bezctx_curveto(bezctx *bc, double x1, double y1, double x2, double y2,
+		    double x3, double y3)
+{
+    bc->curveto(bc, x1, y1, x2, y2, x3, y3);
+}
+
+void bezctx_mark_knot(bezctx *bc, int knot_idx)
+{
+    if (bc->mark_knot)
+	bc->mark_knot(bc, knot_idx);
+}
diff --git a/app/cornu/bezctx.h b/app/cornu/bezctx.h
new file mode 100644
index 0000000..1216be9
--- /dev/null
+++ b/app/cornu/bezctx.h
@@ -0,0 +1,13 @@
+#ifndef _BEZCTX_H
+#define _BEZCTX_H
+#include "bezctx_intf.h"
+
+struct _bezctx {
+    void (*moveto)(bezctx *bc, double x, double y, int is_open);
+    void (*lineto)(bezctx *bc, double x, double y);
+    void (*quadto)(bezctx *bc, double x1, double y1, double x2, double y2);
+    void (*curveto)(bezctx *bc, double x1, double y1, double x2, double y2,
+		    double x3, double y3);
+    void (*mark_knot)(bezctx *bc, int knot_idx);
+};
+#endif
diff --git a/app/cornu/bezctx_intf.h b/app/cornu/bezctx_intf.h
new file mode 100644
index 0000000..4e488c6
--- /dev/null
+++ b/app/cornu/bezctx_intf.h
@@ -0,0 +1,23 @@
+#ifndef _BEZCTX_INTF_H
+#define _BEZCTX_INTF_H
+typedef struct _bezctx bezctx;
+
+bezctx *
+new_bezctx(void);
+
+void
+bezctx_moveto(bezctx *bc, double x, double y, int is_open);
+
+void
+bezctx_lineto(bezctx *bc, double x, double y);
+
+void
+bezctx_quadto(bezctx *bc, double x1, double y1, double x2, double y2);
+
+void
+bezctx_curveto(bezctx *bc, double x1, double y1, double x2, double y2,
+	       double x3, double y3);
+
+void
+bezctx_mark_knot(bezctx *bc, int knot_idx);
+#endif
diff --git a/app/cornu/bezctx_xtrkcad.c b/app/cornu/bezctx_xtrkcad.c
new file mode 100644
index 0000000..1b902b2
--- /dev/null
+++ b/app/cornu/bezctx_xtrkcad.c
@@ -0,0 +1,217 @@
+/*
+xtrkcad_spiro - An adapter for Spiro splines within XtrkCAD.
+
+Copyright (C) XtrkCad.org, based on the Spiro Toolkit of Ralph Levien.
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+
+#include "zmisc.h"
+#include "common.h"
+#include "bezctx.h"
+#include "bezctx_xtrkcad.h"
+#include "track.h"
+#include "tbezier.h"
+#include "i18n.h"
+#include "math.h"
+#include "utility.h"
+
+#define trkSeg(N) DYNARR_N(trkSeg_t,* bc->segsArray, N );
+
+
+typedef struct {
+    bezctx base;
+    dynArr_t * segsArray;
+    BOOL_T track;
+    BOOL_T is_open;
+    BOOL_T has_NAN;
+    BOOL_T draw_spots;
+    coOrd last_pos;						// For moveTo
+    int ends[2];						//Start and End knot number
+
+} bezctx_xtrkcad;
+
+static void
+bezctx_xtrkcad_moveto(bezctx *z, double x, double y, int is_open) {
+    bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+    bc->last_pos.x = x;
+    bc->last_pos.y = y;
+    if (!(isfinite(x) && isfinite(y))) {
+        	bc->has_NAN = TRUE;
+        	return;
+    }
+}
+
+static void
+bezctx_xtrkcad_lineto(bezctx *z, double x, double y) {
+
+	bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+	if (!(isfinite(x) && isfinite(y))) {
+	    bc->has_NAN = TRUE;
+	}
+	if (!bc->is_open || bc->has_NAN) {
+		bc->last_pos.x = x;
+		bc->last_pos.y = y;
+		return;
+	}
+    DYNARR_APPEND(trkSeg_t,* bc->segsArray,10);
+    trkSeg_p seg = &trkSeg(bc->segsArray->cnt-1);
+    seg->u.l.pos[0].x = bc->last_pos.x;
+    seg->u.l.pos[0].y = bc->last_pos.y;
+    seg->u.l.pos[1].x = x;
+    seg->u.l.pos[1].y = y;
+    seg->u.l.option = 0;
+    seg->width = 0.0;
+    seg->color = wDrawColorBlack;
+    seg->type = SEG_STRTRK;
+    if (seg->bezSegs.ptr) MyFree(seg->bezSegs.ptr);
+    seg->bezSegs.max =0;
+    seg->bezSegs.cnt = 0;
+    seg->bezSegs.ptr = NULL;
+    seg->u.l.angle = FindAngle(seg->u.l.pos[0],seg->u.l.pos[1]);
+    bc->last_pos.x = x;
+    bc->last_pos.y = y;
+
+
+}
+
+static void
+bezctx_xtrkcad_quadto(bezctx *z, double x1, double y1, double x2, double y2)
+{
+    bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+    if ((!isfinite(x1) || !isfinite(y1)
+    			|| !isfinite(x2) || !isfinite(y2))) {
+    	bc->has_NAN = TRUE;
+    }
+    if (!bc->is_open || bc->has_NAN) {
+    	 bc->last_pos.x = x2;
+    	 bc->last_pos.y = y2;
+    	return;
+    }
+    DYNARR_APPEND(trkSeg_t,* bc->segsArray,10);
+    trkSeg_p seg = &trkSeg(bc->segsArray->cnt-1);
+    seg->u.b.pos[0] = bc->last_pos;
+    seg->u.b.pos[1].x = x1;
+    seg->u.b.pos[1].y = y1;
+    seg->u.b.pos[2].x = x1;
+    seg->u.b.pos[2].y = y1;
+    seg->u.b.pos[3].x = x2;
+    seg->u.b.pos[3].y = y2;
+    seg->width = 0.0;
+    seg->color = wDrawColorBlack;
+    seg->type = SEG_BEZTRK;
+    if (seg->bezSegs.ptr) MyFree(seg->bezSegs.ptr);
+    seg->bezSegs.max =0;
+    seg->bezSegs.cnt = 0;
+    seg->bezSegs.ptr = NULL;
+    bc->last_pos.x = x2;
+    bc->last_pos.y = y2;
+
+    FixUpBezierSeg(seg->u.b.pos,seg,bc->track);
+}
+
+static void
+ bezctx_xtrkcad_curveto(bezctx *z, double x1, double y1, double x2, double y2,
+	       double x3, double y3)
+{
+	bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+	if (!(isfinite(x1) && isfinite(y1)
+			&& isfinite(x2) && isfinite(y2)
+			&& isfinite(x3) && isfinite(y3))) {
+		bc->has_NAN = TRUE;
+	}
+	if (!bc->is_open || bc->has_NAN) {
+		bc->last_pos.x = x3;
+		bc->last_pos.y = y3;
+		return;
+	}
+	DYNARR_APPEND(trkSeg_t,* bc->segsArray,10);
+	trkSeg_p seg = &trkSeg(bc->segsArray->cnt-1);
+		seg->u.b.pos[0].x = bc->last_pos.x;
+		seg->u.b.pos[0].y = bc->last_pos.y;
+	    seg->u.b.pos[1].x = x1;
+	    seg->u.b.pos[1].y = y1;
+	    seg->u.b.pos[2].x = x2;
+	    seg->u.b.pos[2].y = y2;
+	    seg->u.b.pos[3].x = x3;
+	    seg->u.b.pos[3].y = y3;
+	    seg->width = 0.0;
+	    seg->color = wDrawColorBlack;
+	    seg->type = SEG_BEZTRK;
+	    if (seg->bezSegs.ptr) MyFree(seg->bezSegs.ptr);
+	    seg->bezSegs.max = 0;
+	    seg->bezSegs.cnt = 0;
+	    seg->bezSegs.ptr = NULL;
+	    bc->last_pos.x = x3;
+	    bc->last_pos.y = y3;
+
+	    FixUpBezierSeg(seg->u.b.pos,seg,bc->track);
+
+	 if (bc->draw_spots) {
+		 DYNARR_APPEND(trkSeg_t,* bc->segsArray,10);
+		 seg = &trkSeg(bc->segsArray->cnt-1);
+	 	 seg->type=SEG_FILCRCL;
+	 	 seg->u.c.center.x = bc->last_pos.x;
+	 	 seg->u.c.center.y = bc->last_pos.y;
+	 	 seg->u.c.radius = 0.25;
+	 	 seg->width = 0.0;
+	 	 seg->color = wDrawColorBlack;
+	 }
+
+}
+
+void
+bezctx_xtrkcad_mark_knot(bezctx *z, int knot_idx) {
+
+	bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+	if (knot_idx >= bc->ends[0]) bc->is_open = TRUE;    //Only worry about segs inside our gap
+	if (knot_idx >= bc->ends[1]) bc->is_open = FALSE;
+
+}
+
+
+
+bezctx *
+new_bezctx_xtrkcad(dynArr_t * segArray, int ends[2], BOOL_T spots) {
+
+    bezctx_xtrkcad *result = znew(bezctx_xtrkcad, 1);
+
+    result->segsArray = segArray;
+    result->ends[0] = ends[0];
+    result->ends[1] = ends[1];
+
+    result->base.moveto = bezctx_xtrkcad_moveto;
+    result->base.lineto = bezctx_xtrkcad_lineto;
+    result->base.quadto = bezctx_xtrkcad_quadto;
+    result->base.curveto = bezctx_xtrkcad_curveto;
+    result->base.mark_knot = bezctx_xtrkcad_mark_knot;
+    result->is_open = FALSE;
+    result->has_NAN = FALSE;
+    result->draw_spots = spots;
+    result->track = TRUE;
+
+    return &result->base;
+}
+
+BOOL_T bezctx_xtrkcad_close(bezctx *z) {
+	bezctx_xtrkcad *bc = (bezctx_xtrkcad *)z;
+	if (bc->has_NAN) return FALSE;
+	return TRUE;
+}
+
+
diff --git a/app/cornu/bezctx_xtrkcad.h b/app/cornu/bezctx_xtrkcad.h
new file mode 100644
index 0000000..4117870
--- /dev/null
+++ b/app/cornu/bezctx_xtrkcad.h
@@ -0,0 +1,4 @@
+bezctx * new_bezctx_xtrkcad(dynArr_t * segs, int ends[2], BOOL_T spots);
+
+void bezctx_to_xtrkcad(bezctx *bc);
+BOOL_T bezctx_xtrkcad_close(bezctx *bc);
diff --git a/app/cornu/ppedit_gtk1.c b/app/cornu/ppedit_gtk1.c
new file mode 100644
index 0000000..b81299e
--- /dev/null
+++ b/app/cornu/ppedit_gtk1.c
@@ -0,0 +1,930 @@
+/*
+ppedit - A pattern plate editor for Spiro splines.
+Copyright (C) 2007 Raph Levien
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+#include <gdk/gdkkeysyms.h>
+#include <gtk/gtk.h>
+#include <libart_lgpl/libart.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "zmisc.h"
+#include "bezctx.h"
+#include "bezctx_libart.h"
+#include "bezctx_ps.h"
+#include "cornu.h"
+#include "spiro.h"
+#include "plate.h"
+#include "image.h"
+
+int n_iter = 10;
+
+typedef struct {
+    const char *description;
+    plate *p;
+} undo_record;
+
+typedef struct {
+    GtkWidget *da;
+    const char *description;
+    plate *p;
+    int undo_n;
+    int undo_i;
+    undo_record undo_buf[16];
+    int undo_xn_state;
+
+    GtkWidget *undo_me;
+    GtkWidget *redo_me;
+
+    GtkWidget *show_knots_me;
+    GtkWidget *show_bg_me;
+    int show_knots;
+    int show_bg;
+
+    image *bg_image;
+} plate_edit;
+
+int
+quit_func(GtkWidget *widget, gpointer dummy)
+{
+    gtk_main_quit();
+    return TRUE;
+}
+
+#define C1 0.55228
+static void
+draw_dot(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	 double x, double y, double r, guint32 rgba)
+{
+    ArtBpath bp[6];
+    ArtVpath *vp;
+    ArtSVP *svp;
+
+    bp[0].code = ART_MOVETO;
+    bp[0].x3 = x + r;
+    bp[0].y3 = y;
+    bp[1].code = ART_CURVETO;
+    bp[1].x1 = x + r;
+    bp[1].y1 = y - C1 * r;
+    bp[1].x2 = x + C1 * r;
+    bp[1].y2 = y - r;
+    bp[1].x3 = x;
+    bp[1].y3 = y - r;
+    bp[2].code = ART_CURVETO;
+    bp[2].x1 = x - C1 * r;
+    bp[2].y1 = y - r;
+    bp[2].x2 = x - r;
+    bp[2].y2 = y - C1 * r;
+    bp[2].x3 = x - r;
+    bp[2].y3 = y;
+    bp[3].code = ART_CURVETO;
+    bp[3].x1 = x - r;
+    bp[3].y1 = y + C1 * r;
+    bp[3].x2 = x - C1 * r;
+    bp[3].y2 = y + r;
+    bp[3].x3 = x;
+    bp[3].y3 = y + r;
+    bp[4].code = ART_CURVETO;
+    bp[4].x1 = x + C1 * r;
+    bp[4].y1 = y + r;
+    bp[4].x2 = x + r;
+    bp[4].y2 = y + C1 * r;
+    bp[4].x3 = x + r;
+    bp[4].y3 = y;
+    bp[5].code = ART_END;
+
+    vp = art_bez_path_to_vec(bp, 0.25);
+    svp = art_svp_from_vpath(vp);
+    art_free(vp);
+
+    art_rgb_svp_alpha(svp, x0, y0, x1, y1, rgba, buf, rowstride, NULL);
+    art_svp_free(svp);
+}
+
+static void
+draw_raw_rect(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	      double rx0, double ry0, double rx1, double ry1, guint32 rgba)
+{
+    ArtVpath vp[6];
+    ArtSVP *svp;
+
+    vp[0].code = ART_MOVETO;
+    vp[0].x = rx0;
+    vp[0].y = ry1;
+    vp[1].code = ART_LINETO;
+    vp[1].x = rx1;
+    vp[1].y = ry1;
+    vp[2].code = ART_LINETO;
+    vp[2].x = rx1;
+    vp[2].y = ry0;
+    vp[3].code = ART_LINETO;
+    vp[3].x = rx0;
+    vp[3].y = ry0;
+    vp[4].code = ART_LINETO;
+    vp[4].x = rx0;
+    vp[4].y = ry1;
+    vp[5].code = ART_END;
+
+    svp = art_svp_from_vpath(vp);
+
+    art_rgb_svp_alpha(svp, x0, y0, x1, y1, rgba, buf, rowstride, NULL);
+    art_svp_free(svp);
+}
+
+static void
+draw_rect(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	 double x, double y, double r, guint32 rgba)
+{
+    draw_raw_rect(buf, x0, y0, x1, y1, rowstride,
+		  x - r, y - r, x + r, y + r, rgba);
+}
+
+static void
+draw_half(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	 double x, double y, double r, double th, guint32 rgba)
+{
+    ArtBpath bp[6];
+    ArtVpath *vp;
+    ArtSVP *svp;
+    double c = cos(th);
+    double s = sin(th);
+
+    bp[0].code = ART_MOVETO;
+    bp[0].x3 = x + c * r;
+    bp[0].y3 = y + s * r;
+    bp[1].code = ART_CURVETO;
+    bp[1].x1 = x + c * r + C1 * s * r;
+    bp[1].y1 = y + s * r - C1 * c * r;
+    bp[1].x2 = x + s * r + C1 * c * r;
+    bp[1].y2 = y - c * r + C1 * s * r;
+    bp[1].x3 = x + s * r;
+    bp[1].y3 = y - c * r;
+    bp[2].code = ART_CURVETO;
+    bp[2].x1 = x + s * r - C1 * c * r;
+    bp[2].y1 = y - c * r - C1 * s * r;
+    bp[2].x2 = x - c * r + C1 * s * r;
+    bp[2].y2 = y - s * r - C1 * c * r;
+    bp[2].x3 = x - c * r;
+    bp[2].y3 = y - s * r;
+    bp[3].code = ART_LINETO;
+    bp[3].x3 = x + c * r;
+    bp[3].y3 = y + s * r;
+    bp[4].code = ART_END;
+
+    vp = art_bez_path_to_vec(bp, 0.25);
+    svp = art_svp_from_vpath(vp);
+    art_free(vp);
+
+    art_rgb_svp_alpha(svp, x0, y0, x1, y1, rgba, buf, rowstride, NULL);
+    art_svp_free(svp);
+}
+
+static ArtVpath *
+bezctx_to_vpath(bezctx *bc)
+{
+    ArtBpath *bp = bezctx_to_bpath(bc);
+    ArtVpath *vp = art_bez_path_to_vec(bp, .25);
+
+    g_free(bp);
+    if (vp[0].code == ART_END || vp[1].code == ART_END) {
+	g_free(vp);
+	vp = NULL;
+    }
+    return vp;
+}
+
+static void
+draw_plate(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	   plate_edit *pe)
+{
+    plate *p = pe->p;
+    int i, j;
+
+    /* find an existing point to select, if any */
+    for (i = 0; i < p->n_sp; i++) {
+	bezctx *bc = new_bezctx_libart();
+	subpath *sp = &p->sp[i];
+	spiro_seg *s = draw_subpath(sp, bc);
+	ArtVpath *vp = bezctx_to_vpath(bc);
+
+	if (vp != NULL) {
+	    ArtSVP *svp = art_svp_vpath_stroke(vp, ART_PATH_STROKE_JOIN_MITER,
+					       ART_PATH_STROKE_CAP_BUTT,
+					       1.5, 4.0, 0.25);
+
+	    art_free(vp);
+	    art_rgb_svp_alpha(svp, x0, y0, x1, y1, 0x000000ff, buf, rowstride,
+			      NULL);
+	    art_svp_free(svp);
+	}
+
+	for (j = 0; j < sp->n_kt; j++) {
+	    if (pe->show_knots) {
+		knot *kt = &sp->kt[j];
+		kt_flags kf = kt->flags;
+		if ((kf & KT_SELECTED) && (kf & KT_OPEN)) {
+		    draw_dot(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			     3, 0x000000ff);
+		    draw_dot(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			     1.5, 0xffffffff);
+		} else if ((kf & KT_SELECTED) && (kf & KT_CORNER)) {
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      3, 0x000000ff);
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      1.5, 0xffffffff);
+		} else if (!(kf & KT_SELECTED) && (kf & KT_CORNER)) {
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      2.5, 0x000080ff);
+		} else if ((kf & KT_SELECTED) && (kf & KT_CORNU)) {
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      3, 0xc000c0ff);
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      1.5, 0xffffffff);
+		} else if (!(kf & KT_SELECTED) && (kf & KT_CORNU)) {
+		    draw_rect(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			      2.5, 0x800080ff);
+		} else if ((kf & KT_LEFT) || (kf & KT_RIGHT)) {
+		    double th = 1.5708 + (s ? get_knot_th(s, j) : 0);
+		    if (kf & KT_LEFT)
+			th += 3.1415926;
+		    if (kf & KT_SELECTED) {
+			draw_half(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+				 4, th, 0x000000ff);
+			draw_half(buf, x0, y0, x1, y1, rowstride,
+				  kt->x + sin(th), kt->y - cos(th),
+				  2, th, 0xffffffff);
+		    } else {
+			draw_half(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+				  3, th, 0x000080ff);
+		    }
+		} else {
+		    draw_dot(buf, x0, y0, x1, y1, rowstride, kt->x, kt->y,
+			     2, 0x000080ff);
+		}
+	    }
+	}
+	free_spiro(s);
+    }
+}
+
+static void
+draw_selection(art_u8 *buf, int x0, int y0, int x1, int y1, int rowstride,
+	       plate_edit *pe)
+{
+    plate *p = pe->p;
+
+    if (p->motmode == MOTION_MODE_SELECT) {
+	double rx0 = p->sel_x0;
+	double ry0 = p->sel_y0;
+	double rx1 = p->x0;
+	double ry1 = p->y0;
+	if (rx0 > rx1) {
+	    double tmp = rx1;
+	    rx1 = rx0;
+	    rx0 = tmp;
+	}
+	if (ry0 > ry1) {
+	    double tmp = ry1;
+	    ry1 = ry0;
+	    ry0 = tmp;
+	}
+	if (rx1 > rx0 && ry1 > ry0)
+	    draw_raw_rect(buf, x0, y0, x1, y1, rowstride,
+			  rx0, ry0, rx1, ry1, 0x0000ff20);
+    }
+}
+
+static void
+render_bg_layer(guchar *buf, int rowstride, int x0, int y0, int x1, int y1,
+		plate_edit *pe)
+{
+    const double affine[6] = { 1, 0, 0, 1, 0, 0 };
+
+    if (pe->show_bg && pe->bg_image)
+	render_image(pe->bg_image, affine,
+		     buf, rowstride, x0, y0, x1, y1);
+    else
+	memset(buf, 255, (y1 - y0) * rowstride);
+}
+
+static gint
+data_expose (GtkWidget *widget, GdkEventExpose *event, void *data)
+{
+    plate_edit *pe = (plate_edit *)data;
+    int x0 = event->area.x;
+    int y0 = event->area.y;
+    int width = event->area.width;
+    int height = event->area.height;
+    guchar *rgb;
+    int rowstride = (width * 3 + 3) & -4;
+
+    rgb = g_new (guchar, event->area.height * rowstride);
+
+    render_bg_layer(rgb, rowstride, x0, y0, x0 + width, y0 + height, pe);
+
+    draw_plate(rgb, x0, y0, x0 + width, y0 + height, rowstride, pe);
+
+    draw_selection(rgb, x0, y0, x0 + width, y0 + height, rowstride, pe);
+
+    gdk_draw_rgb_image(widget->window,
+		       widget->style->black_gc,
+		       x0, y0, width, height,
+		       GDK_RGB_DITHER_NONE, rgb,
+		       rowstride);
+    g_free(rgb);
+    return FALSE;
+}
+
+/* Make sure there's room for at least one more undo record. */
+static void
+makeroom_undo(plate_edit *pe)
+{
+    const int undo_max = sizeof(pe->undo_buf) / sizeof(undo_record);
+
+    if (pe->undo_n == undo_max) {
+	free_plate(pe->undo_buf[0].p);
+	memmove(pe->undo_buf, pe->undo_buf + 1, (undo_max - 1) * sizeof(undo_record));
+	pe->undo_i--;
+	pe->undo_n--;
+    }
+}
+
+static void
+set_undo_menuitem(GtkWidget *me, const char *name, const char *desc)
+{
+    char str[256];
+
+    if (desc) {
+	sprintf(str, "%s %s", name, desc);
+    } else {
+	strcpy(str, name);
+    }
+    gtk_container_foreach(GTK_CONTAINER(me),
+			  (GtkCallback)gtk_label_set_text,
+			  str);
+    gtk_widget_set_sensitive(me, desc != NULL);
+}
+
+static void
+set_undo_state(plate_edit *pe, const char *undo_desc, const char *redo_desc)
+{
+    set_undo_menuitem(pe->undo_me, "Undo", undo_desc);
+    set_undo_menuitem(pe->redo_me, "Redo", redo_desc);
+}
+
+static void
+begin_undo_xn(plate_edit *pe)
+{
+    int i;
+
+    if (pe->undo_xn_state != 1) {
+	for (i = pe->undo_i; i < pe->undo_n; i++)
+	    free_plate(pe->undo_buf[i].p);
+	pe->undo_n = pe->undo_i;
+	makeroom_undo(pe);
+	i = pe->undo_i;
+	pe->undo_buf[i].description = pe->description;
+	pe->undo_buf[i].p = copy_plate(pe->p);
+	pe->undo_n = i + 1;
+	pe->undo_xn_state = 1;
+    }
+}
+
+static void
+dirty_undo_xn(plate_edit *pe, const char *description)
+{
+    if (pe->undo_xn_state == 0) {
+	g_warning("dirty_undo_xn: not in begin_undo_xn state");
+	begin_undo_xn(pe);
+    }
+    if (description == NULL)
+	description = pe->p->description;
+    if (pe->undo_xn_state == 1) {
+	pe->undo_i++;
+	pe->undo_xn_state = 2;
+	set_undo_state(pe, description, NULL);
+    }
+    pe->description = description;
+}
+
+static void
+begindirty_undo_xn(plate_edit *pe, const char *description)
+{
+    begin_undo_xn(pe);
+    dirty_undo_xn(pe, description);
+}
+ 
+static void
+end_undo_xn(plate_edit *pe)
+{
+    if (pe->undo_xn_state == 0) {
+	g_warning("end_undo_xn: not in undo xn");
+    }
+    pe->undo_xn_state = 0;
+}
+
+static int
+undo(plate_edit *pe)
+{
+    if (pe->undo_i == 0)
+	return 0;
+
+    if (pe->undo_i == pe->undo_n) {
+	makeroom_undo(pe);
+	pe->undo_buf[pe->undo_i].description = pe->description;
+	pe->undo_buf[pe->undo_i].p = pe->p;
+	pe->undo_n++;
+    } else {
+	free_plate(pe->p);
+    }
+    pe->undo_i--;
+    pe->description = pe->undo_buf[pe->undo_i].description;
+    set_undo_state(pe,
+		   pe->undo_i > 0 ? pe->description : NULL,
+		   pe->undo_buf[pe->undo_i + 1].description);
+    g_print("undo: %d of %d\n", pe->undo_i, pe->undo_n);
+    pe->p = copy_plate(pe->undo_buf[pe->undo_i].p);
+    return 1;
+}
+
+static int
+redo(plate_edit *pe)
+{
+    if (pe->undo_i >= pe->undo_n - 1)
+	return 0;
+    free_plate(pe->p);
+    pe->undo_i++;
+    set_undo_state(pe,
+		   pe->undo_buf[pe->undo_i].description,
+		   pe->undo_i < pe->undo_n - 1 ?
+		   pe->undo_buf[pe->undo_i + 1].description : NULL);
+    pe->description = pe->undo_buf[pe->undo_i].description;
+    pe->p = copy_plate(pe->undo_buf[pe->undo_i].p);
+    g_print("redo: %d of %d\n", pe->undo_i, pe->undo_n);
+    return 1;
+}
+
+static gint
+data_button_press (GtkWidget *widget, GdkEventButton *event, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+    plate *p = pe->p;
+    double x, y;
+    press_mod mods = 0;
+
+#define noVERBOSE
+#ifdef VERBOSE
+    g_print ("button press %f %f %f %d\n",
+	     event->x, event->y, event->pressure, event->type);
+
+#endif
+    x = event->x;
+    y = event->y;
+    if (event->state & GDK_SHIFT_MASK) mods |= PRESS_MOD_SHIFT;
+    if (event->state & GDK_CONTROL_MASK) mods |= PRESS_MOD_CTRL;
+    if (event->type == GDK_2BUTTON_PRESS) mods |= PRESS_MOD_DOUBLE;
+    if (event->type == GDK_3BUTTON_PRESS) mods |= PRESS_MOD_TRIPLE;
+
+    begin_undo_xn(pe);
+    p->description = NULL;
+    plate_press(p, x, y, mods);
+    if (p->description) dirty_undo_xn(pe, NULL);
+    gtk_widget_queue_draw(widget);
+
+    return TRUE;
+}
+
+static gint
+data_motion_move (GtkWidget *widget, GdkEventMotion *event, plate_edit *pe)
+{
+    double x, y;
+    x = event->x;
+    y = event->y;
+
+    plate_motion_move(pe->p, x, y);
+    dirty_undo_xn(pe, NULL);
+    
+    gtk_widget_queue_draw(widget);
+
+    return TRUE;
+}
+
+static gint
+data_motion_select (GtkWidget *widget, GdkEventMotion *event, plate_edit *pe)
+{
+    double x, y;
+
+    x = event->x;
+    y = event->y;
+    
+    plate_motion_select(pe->p, x, y);
+    gtk_widget_queue_draw(widget);
+
+    return TRUE;
+}
+
+static gint
+data_motion (GtkWidget *widget, GdkEventMotion *event, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+#ifdef VERBOSE
+    g_print ("motion %f %f %f\n", event->x, event->y, event->pressure);
+
+#endif
+    if (pe->p->motmode == MOTION_MODE_MOVE)
+	return data_motion_move(widget, event, pe);
+    else if (pe->p->motmode == MOTION_MODE_SELECT)
+	return data_motion_select(widget, event, pe);
+    return TRUE;
+}
+
+static gint
+data_button_release (GtkWidget *widget, GdkEventMotion *event, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+    int need_redraw;
+    
+    need_redraw = (pe->p->motmode == MOTION_MODE_SELECT);
+
+    plate_unpress(pe->p);
+    
+    gtk_widget_queue_draw(widget);
+
+    return TRUE;
+}
+
+static gboolean
+key_press(GtkWidget *widget, GdkEventKey *event, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+    int delta = event->state & 4 ? 10 : 1;
+    int old_n_iter = n_iter;
+    double dx = 0, dy = 0;
+    gboolean did_something = FALSE;
+
+    g_print("key press %d %s %d\n", event->keyval, event->string, event->state);
+
+    if (event->keyval == '<') {
+	did_something = TRUE;
+	n_iter -= delta;
+    } else if (event->keyval == '>') {
+	n_iter += delta;
+    }
+    if (n_iter < 0) n_iter = 0;
+    if (n_iter != old_n_iter)
+	g_print("n_iter = %d\n", n_iter);
+
+    if (event->keyval == GDK_Left)
+	dx = -1;
+    else if (event->keyval == GDK_Right)
+	dx = 1;
+    else if (event->keyval == GDK_Up)
+	dy = -1;
+    else if (event->keyval == GDK_Down)
+	dy = 1;
+    if (event->state & GDK_SHIFT_MASK) {
+	dx *= 10;
+	dy *= 10;
+    } else if (event->state & GDK_CONTROL_MASK) {
+	dx *= .1;
+	dy *= .1;
+    }
+    if (dx != 0 || dy != 0) {
+	begindirty_undo_xn(pe, "Keyboard move");
+	plate_motion_move(pe->p, pe->p->x0 + dx, pe->p->y0 + dy);
+	end_undo_xn(pe);
+	did_something = TRUE;
+    }
+
+    if (did_something) {
+	gtk_signal_emit_stop_by_name(GTK_OBJECT(widget), "key-press-event");
+	gtk_widget_queue_draw(widget);
+    }
+
+    return did_something;
+}
+
+static gint
+toggle_corner_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    begindirty_undo_xn(pe, "Toggle Corner");
+    plate_toggle_corner(pe->p);
+    end_undo_xn(pe);
+    gtk_widget_queue_draw(pe->da);
+
+    return TRUE;
+}
+
+static gint
+delete_pt_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    begindirty_undo_xn(pe, "Delete Point");
+    plate_delete_pt(pe->p);
+    end_undo_xn(pe);
+    gtk_widget_queue_draw(pe->da);
+
+    return TRUE;
+}
+
+static gint
+set_select_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_SELECT;
+    return TRUE;
+}
+
+static gint
+set_curve_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_ADD_CURVE;
+    pe->p->last_curve_mmode = pe->p->mmode;
+    return TRUE;
+}
+
+static gint
+set_corner_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_ADD_CORNER;
+    return TRUE;
+}
+
+static gint
+set_left_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_ADD_LEFT;
+    return TRUE;
+}
+
+static gint
+set_right_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_ADD_RIGHT;
+    return TRUE;
+}
+
+static gint
+set_cornu_mode_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->p->mmode = MOUSE_MODE_ADD_CORNU;
+    pe->p->last_curve_mmode = pe->p->mmode;
+    return TRUE;
+}
+
+static gint
+undo_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    undo(pe);
+    gtk_widget_queue_draw(pe->da);
+
+    return TRUE;
+}
+
+static gint
+redo_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    redo(pe);
+    gtk_widget_queue_draw(pe->da);
+
+    return TRUE;
+}
+
+static gint
+save_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    file_write_plate("plate", pe->p);
+
+    return TRUE;
+}
+
+static gint
+toggle_show_knots_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->show_knots = !pe->show_knots;
+    gtk_container_foreach(GTK_CONTAINER(pe->show_knots_me),
+			  (GtkCallback)gtk_label_set_text,
+			  pe->show_knots ? "Hide Knots" : "Show Knots");
+    gtk_widget_queue_draw(pe->da);
+    return TRUE;
+}
+
+static gint
+toggle_show_bg_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+
+    pe->show_bg = !pe->show_bg;
+    gtk_container_foreach(GTK_CONTAINER(pe->show_bg_me),
+			  (GtkCallback)gtk_label_set_text,
+			  pe->show_bg ? "Hide BG" : "Show BG");
+    gtk_widget_queue_draw(pe->da);
+    return TRUE;
+}
+
+static gint
+print_func(GtkWidget *widget, gpointer data)
+{
+    plate_edit *pe = (plate_edit *)data;
+    plate *p = pe->p;
+    int i;
+    FILE *f = fopen("/tmp/foo.ps", "w");
+    bezctx *bc = new_bezctx_ps(f);
+
+    fputs(ps_prolog, f);
+    for (i = 0; i < p->n_sp; i++) {
+	subpath *sp = &p->sp[i];
+	free_spiro(draw_subpath(sp, bc));
+    }
+    bezctx_ps_close(bc);
+    fputs(ps_postlog, f);
+    fclose(f);
+    return TRUE;
+}
+
+static GtkWidget *
+add_menuitem(GtkWidget *menu, const char *name, GtkSignalFunc callback,
+	     gpointer callback_data, GtkAccelGroup *ag, const char *accel)
+{
+    GtkWidget *menuitem;
+
+    menuitem = gtk_menu_item_new_with_label(name);
+    gtk_menu_append(GTK_MENU(menu), menuitem);
+    gtk_widget_show(menuitem);
+    if (accel != NULL) {
+	guint accel_key, accel_mods;
+
+	gtk_accelerator_parse(accel, &accel_key, &accel_mods);
+	gtk_widget_add_accelerator(menuitem, "activate", ag,
+				   accel_key, accel_mods, GTK_ACCEL_VISIBLE);
+    }
+    gtk_signal_connect(GTK_OBJECT(menuitem), "activate",
+		       (GtkSignalFunc)callback, callback_data);
+    return (menuitem);
+}
+
+static void
+create_mainwin(plate_edit *p)
+{
+    GtkWidget *mainwin;
+    GtkWidget *eb;
+    GtkWidget *da;
+    GtkWidget *vbox;
+    GtkWidget *menubar;
+    GtkWidget *menu;
+    GtkWidget *menuitem;
+    GtkAccelGroup *ag;
+    void *data = p;
+
+    mainwin = gtk_widget_new(gtk_window_get_type(),
+			  "GtkWindow::type", GTK_WINDOW_TOPLEVEL,
+			  "GtkWindow::title", "pattern plate editor",
+			  NULL);
+    gtk_signal_connect(GTK_OBJECT(mainwin), "destroy",
+		       (GtkSignalFunc)quit_func, NULL);
+
+    vbox = gtk_vbox_new(FALSE, 0);
+    gtk_container_add(GTK_CONTAINER(mainwin), vbox);
+
+    menubar = gtk_menu_bar_new();
+    ag = gtk_accel_group_new();
+    gtk_window_add_accel_group(GTK_WINDOW(mainwin), ag);
+    gtk_box_pack_start(GTK_BOX(vbox), menubar, FALSE, FALSE, 0);
+
+    menu = gtk_menu_new();
+    menuitem = gtk_menu_item_new_with_label("File");
+    gtk_menu_item_set_submenu(GTK_MENU_ITEM(menuitem), menu);
+    gtk_menu_bar_append(GTK_MENU_BAR(menubar), menuitem);
+    gtk_widget_show(menuitem);
+    gtk_menu_set_accel_group(GTK_MENU(menu), ag);
+    add_menuitem(menu, "Save", (GtkSignalFunc)save_func, data, ag, "<ctrl>S");
+    add_menuitem(menu, "Quit", (GtkSignalFunc)quit_func, data, ag, "<ctrl>Q");
+    add_menuitem(menu, "Print", (GtkSignalFunc)print_func, data, ag, "<ctrl>P");
+
+    menu = gtk_menu_new();
+    menuitem = gtk_menu_item_new_with_label("Edit");
+    gtk_menu_item_set_submenu(GTK_MENU_ITEM(menuitem), menu);
+    gtk_menu_bar_append(GTK_MENU_BAR(menubar), menuitem);
+    gtk_widget_show(menuitem);
+    gtk_menu_set_accel_group(GTK_MENU(menu), ag);
+    p->undo_me = add_menuitem(menu, "Undo", (GtkSignalFunc)undo_func, data, ag, "<ctrl>Z");
+    p->redo_me = add_menuitem(menu, "Redo", (GtkSignalFunc)redo_func, data, ag, "<ctrl>Y");
+    set_undo_state(p, NULL, NULL);
+    add_menuitem(menu, "Toggle Corner", (GtkSignalFunc)toggle_corner_func, data, ag, "<ctrl>T");
+    add_menuitem(menu, "Delete Point", (GtkSignalFunc)delete_pt_func, data, ag, "<ctrl>D");
+    add_menuitem(menu, "Selection Mode", (GtkSignalFunc)set_select_mode_func, data, ag, "1");
+    add_menuitem(menu, "Add Curve Mode", (GtkSignalFunc)set_curve_mode_func, data, ag, "2");
+    add_menuitem(menu, "Add Corner Mode", (GtkSignalFunc)set_corner_mode_func, data, ag, "3");
+    add_menuitem(menu, "Add Left Mode", (GtkSignalFunc)set_left_mode_func, data, ag, "4");
+    add_menuitem(menu, "Add Right Mode", (GtkSignalFunc)set_right_mode_func, data, ag, "5");
+    add_menuitem(menu, "Add Cornu Mode", (GtkSignalFunc)set_cornu_mode_func, data, ag, "6");
+
+
+    menu = gtk_menu_new();
+    menuitem = gtk_menu_item_new_with_label("View");
+    gtk_menu_item_set_submenu(GTK_MENU_ITEM(menuitem), menu);
+    gtk_menu_bar_append(GTK_MENU_BAR(menubar), menuitem);
+    gtk_widget_show(menuitem);
+    gtk_menu_set_accel_group(GTK_MENU(menu), ag);
+    p->show_knots_me = add_menuitem(menu, "Hide Knots",
+				    (GtkSignalFunc)toggle_show_knots_func,
+				    data, ag, "<ctrl>K");
+    p->show_bg_me = add_menuitem(menu, "Hide BG",
+				 (GtkSignalFunc)toggle_show_bg_func,
+				 data, ag, "<ctrl>B");
+
+    eb = gtk_event_box_new ();
+    GTK_WIDGET_SET_FLAGS(eb, GTK_CAN_FOCUS);
+    gtk_box_pack_start(GTK_BOX(vbox), eb, TRUE, TRUE, 0);
+    gtk_widget_set_extension_events (eb, GDK_EXTENSION_EVENTS_ALL);
+    gtk_signal_connect(GTK_OBJECT (eb), "button-press-event",
+		       (GtkSignalFunc) data_button_press, data);
+    gtk_signal_connect(GTK_OBJECT (eb), "motion-notify-event",
+		       (GtkSignalFunc) data_motion, data);
+    gtk_signal_connect(GTK_OBJECT (eb), "button-release-event",
+		       (GtkSignalFunc) data_button_release, data);
+    gtk_signal_connect(GTK_OBJECT(eb), "key-press-event",
+		       (GtkSignalFunc)key_press, data);
+
+    da = gtk_drawing_area_new();
+    p->da = da;
+    gtk_window_set_default_size(GTK_WINDOW(mainwin), 512, 512);
+    gtk_container_add(GTK_CONTAINER(eb), da);
+    gtk_signal_connect(GTK_OBJECT (da), "expose-event",
+		       (GtkSignalFunc) data_expose, data);
+#if 0
+    gtk_widget_set_double_buffered(da, FALSE);
+#endif
+
+    gtk_widget_grab_focus(eb);
+    gtk_widget_show(da);
+    gtk_widget_show(eb);
+    gtk_widget_show(menubar);
+    gtk_widget_show(vbox);
+    gtk_widget_show(mainwin);
+}
+
+int main(int argc, char **argv)
+{
+    plate_edit pe;
+    plate *p = NULL;
+    gtk_init(&argc, &argv);
+    gtk_widget_set_default_colormap(gdk_rgb_get_cmap());
+    gtk_widget_set_default_visual(gdk_rgb_get_visual());
+    char *reason;
+
+    if (argc > 1)
+	p = file_read_plate(argv[1]);
+    if (p == NULL)
+	p = new_plate();
+    pe.p = p;
+    pe.undo_n = 0;
+    pe.undo_i = 0;
+    pe.undo_xn_state = 0;
+    pe.show_knots = 1;
+    pe.show_bg = 1;
+    pe.bg_image = load_image_file("/tmp/foo.ppm", &reason);
+    create_mainwin(&pe);
+    gtk_main();
+    return 0;
+}
diff --git a/app/cornu/spiro.c b/app/cornu/spiro.c
new file mode 100644
index 0000000..5aae665
--- /dev/null
+++ b/app/cornu/spiro.c
@@ -0,0 +1,1099 @@
+/*
+ppedit - A pattern plate editor for Spiro splines.
+Copyright (C) 2007 Raph Levien
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+/* C implementation of third-order polynomial spirals. */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "bezctx_intf.h"
+#include "spiro.h"
+
+struct spiro_seg_s {
+    double x;
+    double y;
+    char ty;
+    double bend_th;
+    double ks[4];
+    double seg_ch;
+    double seg_th;
+    double l;
+};
+
+typedef struct {
+    double a[11]; /* band-diagonal matrix */
+    double al[5]; /* lower part of band-diagonal decomposition */
+} bandmat;
+
+#ifndef M_PI
+#define M_PI            3.14159265358979323846  /* pi */
+#endif
+
+int n = 4;
+
+#ifndef ORDER
+#define ORDER 12
+#endif
+
+/* Integrate polynomial spiral curve over range -.5 .. .5. */
+void
+integrate_spiro(const double ks[4], double xy[2])
+{
+#if 0
+    int n = 1024;
+#endif
+    double th1 = ks[0];
+    double th2 = .5 * ks[1];
+    double th3 = (1./6) * ks[2];
+    double th4 = (1./24) * ks[3];
+    double x, y;
+    double ds = 1. / n;
+    double ds2 = ds * ds;
+    double ds3 = ds2 * ds;
+    double k0 = ks[0] * ds;
+    double k1 = ks[1] * ds;
+    double k2 = ks[2] * ds;
+    double k3 = ks[3] * ds;
+    int i;
+    double s = .5 * ds - .5;
+
+    x = 0;
+    y = 0;
+
+    for (i = 0; i < n; i++) {
+
+#if ORDER > 2
+	double u, v;
+	double km0, km1, km2, km3;
+
+	if (n == 1) {
+	    km0 = k0;
+	    km1 = k1 * ds;
+	    km2 = k2 * ds2;
+	} else {
+	    km0 = (((1./6) * k3 * s + .5 * k2) * s + k1) * s + k0;
+	    km1 = ((.5 * k3 * s + k2) * s + k1) * ds;
+	    km2 = (k3 * s + k2) * ds2;
+	}
+	km3 = k3 * ds3;
+#endif
+
+	{
+
+#if ORDER == 4
+	double km0_2 = km0 * km0;
+	u = 24 - km0_2;
+	v = km1;
+#endif
+
+#if ORDER == 6
+	double km0_2 = km0 * km0;
+	double km0_4 = km0_2 * km0_2;
+	u = 24 - km0_2 + (km0_4 - 4 * km0 * km2 - 3 * km1 * km1) * (1./80);
+	v = km1 + (km3 - 6 * km0_2 * km1) * (1./80);
+#endif
+
+#if ORDER == 8
+	double t1_1 = km0;
+	double t1_2 = .5 * km1;
+	double t1_3 = (1./6) * km2;
+	double t1_4 = (1./24) * km3;
+	double t2_2 = t1_1 * t1_1;
+	double t2_3 = 2 * (t1_1 * t1_2);
+	double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
+	double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
+	double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
+	double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
+	double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
+	double t4_4 = t2_2 * t2_2;
+	double t4_5 = 2 * (t2_2 * t2_3);
+	double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
+	double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
+	double t6_6 = t4_4 * t2_2;
+	u = 1;
+	v = 0;
+	v += (1./12) * t1_2 + (1./80) * t1_4;
+	u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6;
+	v -= (1./480) * t3_4 + (1./2688) * t3_6;
+	u += (1./1920) * t4_4 + (1./10752) * t4_6;
+	v += (1./53760) * t5_6;
+	u -= (1./322560) * t6_6;
+#endif
+
+#if ORDER == 10
+	double t1_1 = km0;
+	double t1_2 = .5 * km1;
+	double t1_3 = (1./6) * km2;
+	double t1_4 = (1./24) * km3;
+	double t2_2 = t1_1 * t1_1;
+	double t2_3 = 2 * (t1_1 * t1_2);
+	double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
+	double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
+	double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
+	double t2_7 = 2 * (t1_3 * t1_4);
+	double t2_8 = t1_4 * t1_4;
+	double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
+	double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
+	double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
+	double t4_4 = t2_2 * t2_2;
+	double t4_5 = 2 * (t2_2 * t2_3);
+	double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
+	double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
+	double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
+	double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
+	double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
+	double t6_6 = t4_4 * t2_2;
+	double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
+	double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
+	double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
+	double t8_8 = t6_6 * t2_2;
+	u = 1;
+	v = 0;
+	v += (1./12) * t1_2 + (1./80) * t1_4;
+	u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
+	v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8;
+	u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8;
+	v += (1./53760) * t5_6 + (1./276480) * t5_8;
+	u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8;
+	v -= (1./1.16122e+07) * t7_8;
+	u += (1./9.28973e+07) * t8_8;
+#endif
+
+#if ORDER == 12
+	double t1_1 = km0;
+	double t1_2 = .5 * km1;
+	double t1_3 = (1./6) * km2;
+	double t1_4 = (1./24) * km3;
+	double t2_2 = t1_1 * t1_1;
+	double t2_3 = 2 * (t1_1 * t1_2);
+	double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
+	double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
+	double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
+	double t2_7 = 2 * (t1_3 * t1_4);
+	double t2_8 = t1_4 * t1_4;
+	double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
+	double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
+	double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
+	double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
+	double t4_4 = t2_2 * t2_2;
+	double t4_5 = 2 * (t2_2 * t2_3);
+	double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
+	double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
+	double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
+	double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
+	double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
+	double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
+	double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
+	double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
+	double t6_6 = t4_4 * t2_2;
+	double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
+	double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
+	double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
+	double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
+	double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
+	double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
+	double t8_8 = t6_6 * t2_2;
+	double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
+	double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
+	double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
+	double t10_10 = t8_8 * t2_2;
+	u = 1;
+	v = 0;
+	v += (1./12) * t1_2 + (1./80) * t1_4;
+	u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
+	v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10;
+	u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10;
+	v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10;
+	u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10;
+	v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10;
+	u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10;
+	v += (1./4.08748e+09) * t9_10;
+	u -= (1./4.08748e+10) * t10_10;
+#endif
+
+#if ORDER == 14
+	double t1_1 = km0;
+	double t1_2 = .5 * km1;
+	double t1_3 = (1./6) * km2;
+	double t1_4 = (1./24) * km3;
+	double t2_2 = t1_1 * t1_1;
+	double t2_3 = 2 * (t1_1 * t1_2);
+	double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
+	double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
+	double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
+	double t2_7 = 2 * (t1_3 * t1_4);
+	double t2_8 = t1_4 * t1_4;
+	double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
+	double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
+	double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
+	double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
+	double t3_12 = t2_8 * t1_4;
+	double t4_4 = t2_2 * t2_2;
+	double t4_5 = 2 * (t2_2 * t2_3);
+	double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
+	double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
+	double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
+	double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
+	double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
+	double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
+	double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
+	double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
+	double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
+	double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
+	double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
+	double t6_6 = t4_4 * t2_2;
+	double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
+	double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
+	double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
+	double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
+	double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
+	double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
+	double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
+	double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
+	double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
+	double t8_8 = t6_6 * t2_2;
+	double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
+	double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
+	double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
+	double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
+	double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
+	double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
+	double t10_10 = t8_8 * t2_2;
+	double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
+	double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
+	double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
+	double t12_12 = t10_10 * t2_2;
+	u = 1;
+	v = 0;
+	v += (1./12) * t1_2 + (1./80) * t1_4;
+	u -= (1./24) * t2_2 + (1./160) * t2_4 + (1./896) * t2_6 + (1./4608) * t2_8;
+	v -= (1./480) * t3_4 + (1./2688) * t3_6 + (1./13824) * t3_8 + (1./67584) * t3_10 + (1./319488) * t3_12;
+	u += (1./1920) * t4_4 + (1./10752) * t4_6 + (1./55296) * t4_8 + (1./270336) * t4_10 + (1./1.27795e+06) * t4_12;
+	v += (1./53760) * t5_6 + (1./276480) * t5_8 + (1./1.35168e+06) * t5_10 + (1./6.38976e+06) * t5_12;
+	u -= (1./322560) * t6_6 + (1./1.65888e+06) * t6_8 + (1./8.11008e+06) * t6_10 + (1./3.83386e+07) * t6_12;
+	v -= (1./1.16122e+07) * t7_8 + (1./5.67706e+07) * t7_10 + (1./2.6837e+08) * t7_12;
+	u += (1./9.28973e+07) * t8_8 + (1./4.54164e+08) * t8_10 + (1./2.14696e+09) * t8_12;
+	v += (1./4.08748e+09) * t9_10 + (1./1.93226e+10) * t9_12;
+	u -= (1./4.08748e+10) * t10_10 + (1./1.93226e+11) * t10_12;
+	v -= (1./2.12549e+12) * t11_12;
+	u += (1./2.55059e+13) * t12_12;
+#endif
+
+#if ORDER == 16
+	double t1_1 = km0;
+	double t1_2 = .5 * km1;
+	double t1_3 = (1./6) * km2;
+	double t1_4 = (1./24) * km3;
+	double t2_2 = t1_1 * t1_1;
+	double t2_3 = 2 * (t1_1 * t1_2);
+	double t2_4 = 2 * (t1_1 * t1_3) + t1_2 * t1_2;
+	double t2_5 = 2 * (t1_1 * t1_4 + t1_2 * t1_3);
+	double t2_6 = 2 * (t1_2 * t1_4) + t1_3 * t1_3;
+	double t2_7 = 2 * (t1_3 * t1_4);
+	double t2_8 = t1_4 * t1_4;
+	double t3_4 = t2_2 * t1_2 + t2_3 * t1_1;
+	double t3_6 = t2_2 * t1_4 + t2_3 * t1_3 + t2_4 * t1_2 + t2_5 * t1_1;
+	double t3_8 = t2_4 * t1_4 + t2_5 * t1_3 + t2_6 * t1_2 + t2_7 * t1_1;
+	double t3_10 = t2_6 * t1_4 + t2_7 * t1_3 + t2_8 * t1_2;
+	double t3_12 = t2_8 * t1_4;
+	double t4_4 = t2_2 * t2_2;
+	double t4_5 = 2 * (t2_2 * t2_3);
+	double t4_6 = 2 * (t2_2 * t2_4) + t2_3 * t2_3;
+	double t4_7 = 2 * (t2_2 * t2_5 + t2_3 * t2_4);
+	double t4_8 = 2 * (t2_2 * t2_6 + t2_3 * t2_5) + t2_4 * t2_4;
+	double t4_9 = 2 * (t2_2 * t2_7 + t2_3 * t2_6 + t2_4 * t2_5);
+	double t4_10 = 2 * (t2_2 * t2_8 + t2_3 * t2_7 + t2_4 * t2_6) + t2_5 * t2_5;
+	double t4_11 = 2 * (t2_3 * t2_8 + t2_4 * t2_7 + t2_5 * t2_6);
+	double t4_12 = 2 * (t2_4 * t2_8 + t2_5 * t2_7) + t2_6 * t2_6;
+	double t4_13 = 2 * (t2_5 * t2_8 + t2_6 * t2_7);
+	double t4_14 = 2 * (t2_6 * t2_8) + t2_7 * t2_7;
+	double t5_6 = t4_4 * t1_2 + t4_5 * t1_1;
+	double t5_8 = t4_4 * t1_4 + t4_5 * t1_3 + t4_6 * t1_2 + t4_7 * t1_1;
+	double t5_10 = t4_6 * t1_4 + t4_7 * t1_3 + t4_8 * t1_2 + t4_9 * t1_1;
+	double t5_12 = t4_8 * t1_4 + t4_9 * t1_3 + t4_10 * t1_2 + t4_11 * t1_1;
+	double t5_14 = t4_10 * t1_4 + t4_11 * t1_3 + t4_12 * t1_2 + t4_13 * t1_1;
+	double t6_6 = t4_4 * t2_2;
+	double t6_7 = t4_4 * t2_3 + t4_5 * t2_2;
+	double t6_8 = t4_4 * t2_4 + t4_5 * t2_3 + t4_6 * t2_2;
+	double t6_9 = t4_4 * t2_5 + t4_5 * t2_4 + t4_6 * t2_3 + t4_7 * t2_2;
+	double t6_10 = t4_4 * t2_6 + t4_5 * t2_5 + t4_6 * t2_4 + t4_7 * t2_3 + t4_8 * t2_2;
+	double t6_11 = t4_4 * t2_7 + t4_5 * t2_6 + t4_6 * t2_5 + t4_7 * t2_4 + t4_8 * t2_3 + t4_9 * t2_2;
+	double t6_12 = t4_4 * t2_8 + t4_5 * t2_7 + t4_6 * t2_6 + t4_7 * t2_5 + t4_8 * t2_4 + t4_9 * t2_3 + t4_10 * t2_2;
+	double t6_13 = t4_5 * t2_8 + t4_6 * t2_7 + t4_7 * t2_6 + t4_8 * t2_5 + t4_9 * t2_4 + t4_10 * t2_3 + t4_11 * t2_2;
+	double t6_14 = t4_6 * t2_8 + t4_7 * t2_7 + t4_8 * t2_6 + t4_9 * t2_5 + t4_10 * t2_4 + t4_11 * t2_3 + t4_12 * t2_2;
+	double t7_8 = t6_6 * t1_2 + t6_7 * t1_1;
+	double t7_10 = t6_6 * t1_4 + t6_7 * t1_3 + t6_8 * t1_2 + t6_9 * t1_1;
+	double t7_12 = t6_8 * t1_4 + t6_9 * t1_3 + t6_10 * t1_2 + t6_11 * t1_1;
+	double t7_14 = t6_10 * t1_4 + t6_11 * t1_3 + t6_12 * t1_2 + t6_13 * t1_1;
+	double t8_8 = t6_6 * t2_2;
+	double t8_9 = t6_6 * t2_3 + t6_7 * t2_2;
+	double t8_10 = t6_6 * t2_4 + t6_7 * t2_3 + t6_8 * t2_2;
+	double t8_11 = t6_6 * t2_5 + t6_7 * t2_4 + t6_8 * t2_3 + t6_9 * t2_2;
+	double t8_12 = t6_6 * t2_6 + t6_7 * t2_5 + t6_8 * t2_4 + t6_9 * t2_3 + t6_10 * t2_2;
+	double t8_13 = t6_6 * t2_7 + t6_7 * t2_6 + t6_8 * t2_5 + t6_9 * t2_4 + t6_10 * t2_3 + t6_11 * t2_2;
+	double t8_14 = t6_6 * t2_8 + t6_7 * t2_7 + t6_8 * t2_6 + t6_9 * t2_5 + t6_10 * t2_4 + t6_11 * t2_3 + t6_12 * t2_2;
+	double t9_10 = t8_8 * t1_2 + t8_9 * t1_1;
+	double t9_12 = t8_8 * t1_4 + t8_9 * t1_3 + t8_10 * t1_2 + t8_11 * t1_1;
+	double t9_14 = t8_10 * t1_4 + t8_11 * t1_3 + t8_12 * t1_2 + t8_13 * t1_1;
+	double t10_10 = t8_8 * t2_2;
+	double t10_11 = t8_8 * t2_3 + t8_9 * t2_2;
+	double t10_12 = t8_8 * t2_4 + t8_9 * t2_3 + t8_10 * t2_2;
+	double t10_13 = t8_8 * t2_5 + t8_9 * t2_4 + t8_10 * t2_3 + t8_11 * t2_2;
+	double t10_14 = t8_8 * t2_6 + t8_9 * t2_5 + t8_10 * t2_4 + t8_11 * t2_3 + t8_12 * t2_2;
+	double t11_12 = t10_10 * t1_2 + t10_11 * t1_1;
+	double t11_14 = t10_10 * t1_4 + t10_11 * t1_3 + t10_12 * t1_2 + t10_13 * t1_1;
+	double t12_12 = t10_10 * t2_2;
+	double t12_13 = t10_10 * t2_3 + t10_11 * t2_2;
+	double t12_14 = t10_10 * t2_4 + t10_11 * t2_3 + t10_12 * t2_2;
+	double t13_14 = t12_12 * t1_2 + t12_13 * t1_1;
+	double t14_14 = t12_12 * t2_2;
+	u = 1;
+	u -= 1./24 * t2_2 + 1./160 * t2_4 + 1./896 * t2_6 + 1./4608 * t2_8;
+	u += 1./1920 * t4_4 + 1./10752 * t4_6 + 1./55296 * t4_8 + 1./270336 * t4_10 + 1./1277952 * t4_12 + 1./5898240 * t4_14;
+	u -= 1./322560 * t6_6 + 1./1658880 * t6_8 + 1./8110080 * t6_10 + 1./38338560 * t6_12 + 1./176947200 * t6_14;
+	u += 1./92897280 * t8_8 + 1./454164480 * t8_10 + 4.6577500191e-10 * t8_12 + 1.0091791708e-10 * t8_14;
+	u -= 2.4464949595e-11 * t10_10 + 5.1752777990e-12 * t10_12 + 1.1213101898e-12 * t10_14;
+	u += 3.9206649992e-14 * t12_12 + 8.4947741650e-15 * t12_14;
+	u -= 4.6674583324e-17 * t14_14;
+	v = 0;
+	v += 1./12 * t1_2 + 1./80 * t1_4;
+	v -= 1./480 * t3_4 + 1./2688 * t3_6 + 1./13824 * t3_8 + 1./67584 * t3_10 + 1./319488 * t3_12;
+	v += 1./53760 * t5_6 + 1./276480 * t5_8 + 1./1351680 * t5_10 + 1./6389760 * t5_12 + 1./29491200 * t5_14;
+	v -= 1./11612160 * t7_8 + 1./56770560 * t7_10 + 1./268369920 * t7_12 + 8.0734333664e-10 * t7_14;
+	v += 2.4464949595e-10 * t9_10 + 5.1752777990e-11 * t9_12 + 1.1213101898e-11 * t9_14;
+	v -= 4.7047979991e-13 * t11_12 + 1.0193728998e-13 * t11_14;
+	v += 6.5344416654e-16 * t13_14;
+#endif
+
+	}
+
+	if (n == 1) {
+#if ORDER == 2
+	    x = 1;
+	    y = 0;
+#else
+	    x = u;
+	    y = v;
+#endif
+	} else {
+	    double th = (((th4 * s + th3) * s + th2) * s + th1) * s;
+	    double cth = cos(th);
+	    double sth = sin(th);
+
+#if ORDER == 2
+	    x += cth;
+	    y += sth;
+#else
+	    x += cth * u - sth * v;
+	    y += cth * v + sth * u;
+#endif
+	    s += ds;
+	}
+    }
+
+#if ORDER == 4 || ORDER == 6
+    xy[0] = x * (1./24 * ds);
+    xy[1] = y * (1./24 * ds);
+#else
+    xy[0] = x * ds;
+    xy[1] = y * ds;
+#endif
+}
+
+static double
+compute_ends(const double ks[4], double ends[2][4], double seg_ch)
+{
+    double xy[2];
+    double ch, th;
+    double l, l2, l3;
+    double th_even, th_odd;
+    double k0_even, k0_odd;
+    double k1_even, k1_odd;
+    double k2_even, k2_odd;
+
+    integrate_spiro(ks, xy);
+    ch = hypot(xy[0], xy[1]);
+    th = atan2(xy[1], xy[0]);
+    l = ch / seg_ch;
+
+    th_even = .5 * ks[0] + (1./48) * ks[2];
+    th_odd = .125 * ks[1] + (1./384) * ks[3] - th;
+    ends[0][0] = th_even - th_odd;
+    ends[1][0] = th_even + th_odd;
+    k0_even = l * (ks[0] + .125 * ks[2]);
+    k0_odd = l * (.5 * ks[1] + (1./48) * ks[3]);
+    ends[0][1] = k0_even - k0_odd;
+    ends[1][1] = k0_even + k0_odd;
+    l2 = l * l;
+    k1_even = l2 * (ks[1] + .125 * ks[3]);
+    k1_odd = l2 * .5 * ks[2];
+    ends[0][2] = k1_even - k1_odd;
+    ends[1][2] = k1_even + k1_odd;
+    l3 = l2 * l;
+    k2_even = l3 * ks[2];
+    k2_odd = l3 * .5 * ks[3];
+    ends[0][3] = k2_even - k2_odd;
+    ends[1][3] = k2_even + k2_odd;
+
+    return l;
+}
+
+static void
+compute_pderivs(const spiro_seg *s, double ends[2][4], double derivs[4][2][4],
+		int jinc)
+{
+    double recip_d = 2e6;
+    double delta = 1./ recip_d;
+    double try_ks[4];
+    double try_ends[2][4];
+    int i, j, k;
+
+    compute_ends(s->ks, ends, s->seg_ch);
+    for (i = 0; i < jinc; i++) {
+	for (j = 0; j < 4; j++)
+	    try_ks[j] = s->ks[j];
+	try_ks[i] += delta;
+	compute_ends(try_ks, try_ends, s->seg_ch);
+	for (k = 0; k < 2; k++)
+	    for (j = 0; j < 4; j++)
+		derivs[j][k][i] = recip_d * (try_ends[k][j] - ends[k][j]);
+    }
+}
+
+static double
+mod_2pi(double th)
+{
+    double u = th / (2 * M_PI);
+    return 2 * M_PI * (u - floor(u + 0.5));
+}
+
+static spiro_seg *
+setup_path(const spiro_cp *src, int n)
+{
+    int n_seg = src[0].ty == '{' ? n - 1 : n;
+    spiro_seg *r = (spiro_seg *)malloc((n_seg + 1) * sizeof(spiro_seg));
+    int i;
+    int ilast;
+
+    for (i = 0; i < n_seg; i++) {
+	r[i].x = src[i].x;
+	r[i].y = src[i].y;
+	r[i].ty = src[i].ty;
+	r[i].ks[0] = 0.;
+	r[i].ks[1] = 0.;
+	r[i].ks[2] = 0.;
+	r[i].ks[3] = 0.;
+    }
+    r[n_seg].x = src[n_seg % n].x;
+    r[n_seg].y = src[n_seg % n].y;
+    r[n_seg].ty = src[n_seg % n].ty;
+
+    for (i = 0; i < n_seg; i++) {
+	double dx = r[i + 1].x - r[i].x;
+	double dy = r[i + 1].y - r[i].y;
+	r[i].seg_ch = hypot(dx, dy);
+	r[i].seg_th = atan2(dy, dx);
+    }
+
+    ilast = n_seg - 1;
+    for (i = 0; i < n_seg; i++) {
+	if (r[i].ty == '{' || r[i].ty == '}' || r[i].ty == 'v')
+	    r[i].bend_th = 0.;
+	else
+	    r[i].bend_th = mod_2pi(r[i].seg_th - r[ilast].seg_th);
+	ilast = i;
+    }
+    return r;
+}
+
+static void
+bandec11(bandmat *m, int *perm, int n)
+{
+    int i, j, k;
+    int l;
+
+    /* pack top triangle to the left. */
+    for (i = 0; i < 5; i++) {
+	for (j = 0; j < i + 6; j++)
+	    m[i].a[j] = m[i].a[j + 5 - i];
+	for (; j < 11; j++)
+	    m[i].a[j] = 0.;
+    }
+    l = 5;
+    for (k = 0; k < n; k++) {
+	int pivot = k;
+	double pivot_val = m[k].a[0];
+	double pivot_scale;
+
+	l = l < n ? l + 1 : n;
+
+	for (j = k + 1; j < l; j++)
+	    if (fabs(m[j].a[0]) > fabs(pivot_val)) {
+		pivot_val = m[j].a[0];
+		pivot = j;
+	    }
+
+	perm[k] = pivot;
+	if (pivot != k) {
+	    for (j = 0; j < 11; j++) {
+		double tmp = m[k].a[j];
+		m[k].a[j] = m[pivot].a[j];
+		m[pivot].a[j] = tmp;
+	    }
+	}
+
+	if (fabs(pivot_val) < 1e-12) pivot_val = 1e-12;
+	pivot_scale = 1. / pivot_val;
+	for (i = k + 1; i < l; i++) {
+	    double x = m[i].a[0] * pivot_scale;
+	    m[k].al[i - k - 1] = x;
+	    for (j = 1; j < 11; j++)
+		m[i].a[j - 1] = m[i].a[j] - x * m[k].a[j];
+	    m[i].a[10] = 0.;
+	}
+    }
+}
+
+static void
+banbks11(const bandmat *m, const int *perm, double *v, int n)
+{
+    int i, k, l;
+
+    /* forward substitution */
+    l = 5;
+    for (k = 0; k < n; k++) {
+	i = perm[k];
+	if (i != k) {
+	    double tmp = v[k];
+	    v[k] = v[i];
+	    v[i] = tmp;
+	}
+	if (l < n) l++;
+	for (i = k + 1; i < l; i++)
+	    v[i] -= m[k].al[i - k - 1] * v[k];
+    }
+
+    /* back substitution */
+    l = 1;
+    for (i = n - 1; i >= 0; i--) {
+	double x = v[i];
+	for (k = 1; k < l; k++)
+	    x -= m[i].a[k] * v[k + i];
+	v[i] = x / m[i].a[0];
+	if (l < 11) l++;
+    }
+}
+
+int compute_jinc(char ty0, char ty1)
+{
+    if (ty0 == 'o' || ty1 == 'o' ||
+	ty0 == ']' || ty1 == '[')
+	return 4;
+    else if (ty0 == 'c' && ty1 == 'c')
+	return 2;
+    else if (((ty0 == '{' || ty0 == 'v' || ty0 == '[') && ty1 == 'c') ||
+	     (ty0 == 'c' && (ty1 == '}' || ty1 == 'v' || ty1 == ']')))
+	return 1;
+    else
+	return 0;
+}
+
+int count_vec(const spiro_seg *s, int nseg)
+{
+    int i;
+    int n = 0;
+
+    for (i = 0; i < nseg; i++)
+	n += compute_jinc(s[i].ty, s[i + 1].ty);
+    return n;
+}
+
+static void
+add_mat_line(bandmat *m, double *v,
+	     double derivs[4], double x, double y, int j, int jj, int jinc,
+	     int nmat)
+{
+    int k;
+
+    if (jj >= 0) {
+	int joff =  (j + 5 - jj + nmat) % nmat;
+	if (nmat < 6) {
+	    joff = j + 5 - jj;
+	} else if (nmat == 6) {
+	    joff = 2 + (j + 3 - jj + nmat) % nmat;
+	}
+#ifdef VERBOSE
+	printf("add_mat_line j=%d jj=%d jinc=%d nmat=%d joff=%d\n", j, jj, jinc, nmat, joff);
+#endif
+	v[jj] += x;
+	for (k = 0; k < jinc; k++)
+	    m[jj].a[joff + k] += y * derivs[k];
+    }
+}
+
+static double
+spiro_iter(spiro_seg *s, bandmat *m, int *perm, double *v, int n)
+{
+    int cyclic = s[0].ty != '{' && s[0].ty != 'v';
+    int i, j, jj;
+    int nmat = count_vec(s, n);
+    double norm;
+    int n_invert;
+
+    for (i = 0; i < nmat; i++) {
+	v[i] = 0.;
+	for (j = 0; j < 11; j++)
+	    m[i].a[j] = 0.;
+	for (j = 0; j < 5; j++)
+	    m[i].al[j] = 0.;
+    }
+
+    j = 0;
+    if (s[0].ty == 'o')
+	jj = nmat - 2;
+    else if (s[0].ty == 'c')
+	jj = nmat - 1;
+    else
+	jj = 0;
+    for (i = 0; i < n; i++) {
+	char ty0 = s[i].ty;
+	char ty1 = s[i + 1].ty;
+	int jinc = compute_jinc(ty0, ty1);
+	double th = s[i].bend_th;
+	double ends[2][4];
+	double derivs[4][2][4];
+	int jthl = -1, jk0l = -1, jk1l = -1, jk2l = -1;
+	int jthr = -1, jk0r = -1, jk1r = -1, jk2r = -1;
+
+	compute_pderivs(&s[i], ends, derivs, jinc);
+
+	/* constraints crossing left */
+	if (ty0 == 'o' || ty0 == 'c' || ty0 == '[' || ty0 == ']') {
+	    jthl = jj++;
+	    jj %= nmat;
+	    jk0l = jj++;
+	}
+	if (ty0 == 'o') {
+	    jj %= nmat;
+	    jk1l = jj++;
+	    jk2l = jj++;
+	}
+
+	/* constraints on left */
+	if ((ty0 == '[' || ty0 == 'v' || ty0 == '{' || ty0 == 'c') &&
+	    jinc == 4) {
+	    if (ty0 != 'c')
+		jk1l = jj++;
+	    jk2l = jj++;
+	}
+
+	/* constraints on right */
+	if ((ty1 == ']' || ty1 == 'v' || ty1 == '}' || ty1 == 'c') &&
+	    jinc == 4) {
+	    if (ty1 != 'c')
+		jk1r = jj++;
+	    jk2r = jj++;
+	}
+
+	/* constraints crossing right */
+	if (ty1 == 'o' || ty1 == 'c' || ty1 == '[' || ty1 == ']') {
+	    jthr = jj;
+	    jk0r = (jj + 1) % nmat;
+	}
+	if (ty1 == 'o') {
+	    jk1r = (jj + 2) % nmat;
+	    jk2r = (jj + 3) % nmat;
+	}
+
+	add_mat_line(m, v, derivs[0][0], th - ends[0][0], 1, j, jthl, jinc, nmat);
+	add_mat_line(m, v, derivs[1][0], ends[0][1], -1, j, jk0l, jinc, nmat);
+	add_mat_line(m, v, derivs[2][0], ends[0][2], -1, j, jk1l, jinc, nmat);
+	add_mat_line(m, v, derivs[3][0], ends[0][3], -1, j, jk2l, jinc, nmat);
+	add_mat_line(m, v, derivs[0][1], -ends[1][0], 1, j, jthr, jinc, nmat);
+	add_mat_line(m, v, derivs[1][1], -ends[1][1], 1, j, jk0r, jinc, nmat);
+	add_mat_line(m, v, derivs[2][1], -ends[1][2], 1, j, jk1r, jinc, nmat);
+	add_mat_line(m, v, derivs[3][1], -ends[1][3], 1, j, jk2r, jinc, nmat);
+	j += jinc;
+    }
+    if (cyclic) {
+	memcpy(m + nmat, m, sizeof(bandmat) * nmat);
+	memcpy(m + 2 * nmat, m, sizeof(bandmat) * nmat);
+	memcpy(v + nmat, v, sizeof(double) * nmat);
+	memcpy(v + 2 * nmat, v, sizeof(double) * nmat);
+	n_invert = 3 * nmat;
+	j = nmat;
+    } else {
+	n_invert = nmat;
+	j = 0;
+    }
+/*
+#define VERBOSE
+*/
+#ifdef VERBOSE
+    for (i = 0; i < n; i++) {
+	int k;
+	for (k = 0; k < 11; k++)
+	    printf(" %2.4f", m[i].a[k]);
+	printf(": %2.4f\n", v[i]);
+    }
+    printf("---\n");
+#endif
+    bandec11(m, perm, n_invert);
+    banbks11(m, perm, v, n_invert);
+    norm = 0.;
+    for (i = 0; i < n; i++) {
+	char ty0 = s[i].ty;
+	char ty1 = s[i + 1].ty;
+	int jinc = compute_jinc(ty0, ty1);
+	int k;
+
+	for (k = 0; k < jinc; k++) {
+	    double dk = v[j++];
+
+#ifdef VERBOSE
+	    printf("s[%d].ks[%d] += %f\n", i, k, dk);
+#endif
+	    s[i].ks[k] += dk;
+	    norm += dk * dk;
+	}
+    }
+    return norm;
+}
+
+int
+solve_spiro(spiro_seg *s, int nseg)
+{
+    bandmat *m;
+    double *v;
+    int *perm;
+    int nmat = count_vec(s, nseg);
+    int n_alloc = nmat;
+    double norm;
+    int i;
+
+    if (nmat == 0)
+	return 0;
+    if (s[0].ty != '{' && s[0].ty != 'v')
+	n_alloc *= 3;
+    if (n_alloc < 5)
+	n_alloc = 5;
+    m = (bandmat *)malloc(sizeof(bandmat) * n_alloc);
+    v = (double *)malloc(sizeof(double) * n_alloc);
+    perm = (int *)malloc(sizeof(int) * n_alloc);
+
+    for (i = 0; i < 10; i++) {
+	norm = spiro_iter(s, m, perm, v, nseg);
+#ifdef VERBOSE
+	printf("%% norm = %g\n", norm);
+#endif
+	if (norm < 1e-12) break;
+    }
+
+    free(m);
+    free(v);
+    free(perm);
+    return 0;
+}
+
+static void
+spiro_seg_to_bpath(const double ks[4],
+		   double x0, double y0, double x1, double y1,
+		   bezctx *bc, int depth)
+{
+    double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
+	fabs((1./48) * ks[3]);
+
+    if (!(bend > 1e-8)) {
+	bezctx_lineto(bc, x1, y1);
+    } else {
+	double seg_ch = hypot(x1 - x0, y1 - y0);
+	double seg_th = atan2(y1 - y0, x1 - x0);
+	double xy[2];
+	double ch, th;
+	double scale, rot;
+	double th_even, th_odd;
+	double ul, vl;
+	double ur, vr;
+
+	integrate_spiro(ks, xy);
+	ch = hypot(xy[0], xy[1]);
+	th = atan2(xy[1], xy[0]);
+	scale = seg_ch / ch;
+	rot = seg_th - th;
+	if (depth > 10 || bend < 0.5) {
+	    th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
+	    th_odd = (1./48) * ks[2] + .5 * ks[0];
+	    ul = (scale * (1./3)) * cos(th_even - th_odd);
+	    vl = (scale * (1./3)) * sin(th_even - th_odd);
+	    ur = (scale * (1./3)) * cos(th_even + th_odd);
+	    vr = (scale * (1./3)) * sin(th_even + th_odd);
+	    bezctx_curveto(bc, x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
+	} else {
+	    /* subdivide */
+	    double ksub[4];
+	    double thsub;
+	    double xysub[2];
+	    double xmid, ymid;
+	    double cth, sth;
+
+	    ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
+	    ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
+	    ksub[2] = .125 * ks[2] - (1./32) * ks[3];
+	    ksub[3] = (1./16) * ks[3];
+	    thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
+	    cth = .5 * scale * cos(thsub);
+	    sth = .5 * scale * sin(thsub);
+	    integrate_spiro(ksub, xysub);
+	    xmid = x0 + cth * xysub[0] - sth * xysub[1];
+	    ymid = y0 + cth * xysub[1] + sth * xysub[0];
+	    spiro_seg_to_bpath(ksub, x0, y0, xmid, ymid, bc, depth + 1);
+	    ksub[0] += .25 * ks[1] + (1./384) * ks[3];
+	    ksub[1] += .125 * ks[2];
+	    ksub[2] += (1./16) * ks[3];
+	    spiro_seg_to_bpath(ksub, xmid, ymid, x1, y1, bc, depth + 1);
+	}
+    }
+}
+
+spiro_seg *
+run_spiro(const spiro_cp *src, int n)
+{
+    int nseg = src[0].ty == '{' ? n - 1 : n;
+    spiro_seg *s = setup_path(src, n);
+    if (nseg > 1)
+	solve_spiro(s, nseg);
+    return s;
+}
+
+void
+free_spiro(spiro_seg *s)
+{
+    free(s);
+}
+
+void
+spiro_to_bpath(const spiro_seg *s, int n, bezctx *bc)
+{
+    int i;
+    int nsegs = s[n - 1].ty == '}' ? n - 1 : n;
+
+    for (i = 0; i < nsegs; i++) {
+	double x0 = s[i].x;
+	double y0 = s[i].y;
+	double x1 = s[i + 1].x;
+	double y1 = s[i + 1].y;
+
+	if (i == 0)
+	    bezctx_moveto(bc, x0, y0, s[0].ty == '{');
+	bezctx_mark_knot(bc, i);
+	spiro_seg_to_bpath(s[i].ks, x0, y0, x1, y1, bc, 0);
+    }
+}
+
+double
+get_knot_th(const spiro_seg *s, int i)
+{
+    double ends[2][4];
+
+    if (i == 0) {
+	compute_ends(s[i].ks, ends, s[i].seg_ch);
+	return s[i].seg_th - ends[0][0];
+    } else {
+	compute_ends(s[i - 1].ks, ends, s[i - 1].seg_ch);
+	return s[i - 1].seg_th + ends[1][0];
+    }
+}
+
+#ifdef UNIT_TEST
+#include <stdio.h>
+#include <sys/time.h> /* for gettimeofday */
+
+static double
+get_time (void)
+{
+  struct timeval tv;
+  struct timezone tz;
+
+  gettimeofday (&tv, &tz);
+
+  return tv.tv_sec + 1e-6 * tv.tv_usec;
+}
+
+int
+test_integ(void) {
+    double ks[] = {1, 2, 3, 4};
+    double xy[2];
+    double xynom[2];
+    double ch, th;
+    int i, j;
+    int nsubdiv;
+
+    n = ORDER < 6 ? 4096 : 1024;
+    integrate_spiro(ks, xynom);
+    nsubdiv = ORDER < 12 ? 8 : 7;
+    for (i = 0; i < nsubdiv; i++) {
+	double st, en;
+	double err;
+	int n_iter = (1 << (20 - i));
+
+	n = 1 << i;
+	st = get_time();
+	for (j = 0; j < n_iter; j++)
+	    integrate_spiro(ks, xy);
+	en = get_time();
+	err = hypot(xy[0] - xynom[0], xy[1] - xynom[1]);
+#ifdef VERBOSE
+	printf("%d %d %g %g\n", ORDER, n, (en - st) / n_iter, err);
+#endif
+	ch = hypot(xy[0], xy[1]);
+	th = atan2(xy[1], xy[0]);
+#if 0
+	printf("n = %d: integ(%g %g %g %g) = %g %g, ch = %g, th = %g\n", n,
+	       ks[0], ks[1], ks[2], ks[3], xy[0], xy[1], ch, th);
+	printf("%d: %g %g\n", n, xy[0] - xynom[0], xy[1] - xynom[1]);
+#endif
+    }
+    return 0;
+}
+
+void
+print_seg(const double ks[4], double x0, double y0, double x1, double y1)
+{
+    double bend = fabs(ks[0]) + fabs(.5 * ks[1]) + fabs(.125 * ks[2]) +
+	fabs((1./48) * ks[3]);
+
+    if (bend < 1e-8) {
+#ifdef VERBOSE
+	printf("%g %g lineto\n", x1, y1);
+#endif
+    } else {
+	double seg_ch = hypot(x1 - x0, y1 - y0);
+	double seg_th = atan2(y1 - y0, x1 - x0);
+	double xy[2];
+	double ch, th;
+	double scale, rot;
+	double th_even, th_odd;
+	double ul, vl;
+	double ur, vr;
+
+	integrate_spiro(ks, xy);
+	ch = hypot(xy[0], xy[1]);
+	th = atan2(xy[1], xy[0]);
+	scale = seg_ch / ch;
+	rot = seg_th - th;
+	if (bend < 1.) {
+	    th_even = (1./384) * ks[3] + (1./8) * ks[1] + rot;
+	    th_odd = (1./48) * ks[2] + .5 * ks[0];
+	    ul = (scale * (1./3)) * cos(th_even - th_odd);
+	    vl = (scale * (1./3)) * sin(th_even - th_odd);
+	    ur = (scale * (1./3)) * cos(th_even + th_odd);
+	    vr = (scale * (1./3)) * sin(th_even + th_odd);
+#ifdef VERBOSE
+	    printf("%g %g %g %g %g %g curveto\n",
+		   x0 + ul, y0 + vl, x1 - ur, y1 - vr, x1, y1);
+#endif
+	    
+	} else {
+	    /* subdivide */
+	    double ksub[4];
+	    double thsub;
+	    double xysub[2];
+	    double xmid, ymid;
+	    double cth, sth;
+
+	    ksub[0] = .5 * ks[0] - .125 * ks[1] + (1./64) * ks[2] - (1./768) * ks[3];
+	    ksub[1] = .25 * ks[1] - (1./16) * ks[2] + (1./128) * ks[3];
+	    ksub[2] = .125 * ks[2] - (1./32) * ks[3];
+	    ksub[3] = (1./16) * ks[3];
+	    thsub = rot - .25 * ks[0] + (1./32) * ks[1] - (1./384) * ks[2] + (1./6144) * ks[3];
+	    cth = .5 * scale * cos(thsub);
+	    sth = .5 * scale * sin(thsub);
+	    integrate_spiro(ksub, xysub);
+	    xmid = x0 + cth * xysub[0] - sth * xysub[1];
+	    ymid = y0 + cth * xysub[1] + sth * xysub[0];
+	    print_seg(ksub, x0, y0, xmid, ymid);
+	    ksub[0] += .25 * ks[1] + (1./384) * ks[3];
+	    ksub[1] += .125 * ks[2];
+	    ksub[2] += (1./16) * ks[3];
+	    print_seg(ksub, xmid, ymid, x1, y1);
+	}
+    }
+}
+
+void
+print_segs(const spiro_seg *segs, int nsegs)
+{
+    int i;
+
+    for (i = 0; i < nsegs; i++) {
+	double x0 = segs[i].x;
+	double y0 = segs[i].y;
+	double x1 = segs[i + 1].x;
+	double y1 = segs[i + 1].y;
+
+	if (i == 0)
+	    printf("%g %g moveto\n", x0, y0);
+	printf("%% ks = [ %g %g %g %g ]\n",
+	       segs[i].ks[0], segs[i].ks[1], segs[i].ks[2], segs[i].ks[3]);
+	print_seg(segs[i].ks, x0, y0, x1, y1);
+    }
+    printf("stroke\n");
+}
+
+int
+test_curve(void)
+{
+    spiro_cp path[] = {
+	{334, 117, 'v'},
+	{305, 176, 'v'},
+	{212, 142, 'c'},
+	{159, 171, 'c'},
+	{224, 237, 'c'},
+	{347, 335, 'c'},
+	{202, 467, 'c'},
+	{81, 429, 'v'},
+	{114, 368, 'v'},
+	{201, 402, 'c'},
+	{276, 369, 'c'},
+	{218, 308, 'c'},
+	{91, 211, 'c'},
+	{124, 111, 'c'},
+	{229, 82, 'c'}
+    };
+    spiro_seg *segs;
+    int i;
+
+    n = 1;
+    for (i = 0; i < 1000; i++) {
+	segs = setup_path(path, 15);
+	solve_spiro(segs, 15);
+    }
+    printf("100 800 translate 1 -1 scale 1 setlinewidth\n");
+    print_segs(segs, 15);
+    printf("showpage\n");
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+    return test_curve();
+}
+#endif
diff --git a/app/cornu/spiro.h b/app/cornu/spiro.h
new file mode 100644
index 0000000..1bcca37
--- /dev/null
+++ b/app/cornu/spiro.h
@@ -0,0 +1,22 @@
+#ifndef _SPIRO_H
+#define _SPIRO_H
+#include "bezctx_intf.h"
+typedef struct {
+    double x;
+    double y;
+    char ty;
+} spiro_cp;
+
+typedef struct spiro_seg_s spiro_seg;
+
+spiro_seg *
+run_spiro(const spiro_cp *src, int n);
+
+void
+free_spiro(spiro_seg *s);
+
+void
+spiro_to_bpath(const spiro_seg *s, int n, bezctx *bc);
+
+double get_knot_th(const spiro_seg *s, int i);
+#endif
diff --git a/app/cornu/spiroentrypoints.c b/app/cornu/spiroentrypoints.c
new file mode 100644
index 0000000..0128115
--- /dev/null
+++ b/app/cornu/spiroentrypoints.c
@@ -0,0 +1,60 @@
+/*
+libspiro - conversion between spiro control points and bezier's
+Copyright (C) 2007 Raph Levien
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+/* Interface routines to Raph's spiro package. */
+
+#include "spiroentrypoints.h"
+
+void
+SpiroCPsToBezier(spiro_cp *spiros,int n,int isclosed,bezctx *bc)
+{
+    spiro_seg *s;
+
+    if ( n<1 )
+return;
+    if ( !isclosed ) {
+	char oldty_start = spiros[0].ty;
+	char oldty_end   = spiros[n-1].ty;
+	spiros[0].ty = '{';
+	spiros[n-1].ty = '}';
+	s = run_spiro(spiros,n);
+	spiros[n-1].ty = oldty_end;
+	spiros[0].ty = oldty_start;
+    } else
+	s = run_spiro(spiros,n);
+    spiro_to_bpath(s,n,bc);
+    free_spiro(s);
+}
+
+void
+TaggedSpiroCPsToBezier(spiro_cp *spiros,bezctx *bc)
+{
+    spiro_seg *s;
+    int n;
+
+    for ( n=0; spiros[n].ty!='z' && spiros[n].ty!='}'; ++n );
+    if ( spiros[n].ty == '}' ) ++n;
+
+    if ( n<1 )
+return;
+    s = run_spiro(spiros,n);
+    spiro_to_bpath(s,n,bc);
+    free_spiro(s);
+}
diff --git a/app/cornu/spiroentrypoints.h b/app/cornu/spiroentrypoints.h
new file mode 100644
index 0000000..006804c
--- /dev/null
+++ b/app/cornu/spiroentrypoints.h
@@ -0,0 +1,31 @@
+#ifndef _SPIROENTRYPOINTS_H
+# define _SPIROENTRYPOINTS_H
+# include "bezctx_intf.h"
+# include "spiro.h"
+
+    /* Possible values of the "ty" field. */
+#define SPIRO_CORNER	'v'
+#define SPIRO_G4		'o'
+#define SPIRO_G2		'c'
+#define SPIRO_LEFT		'['
+#define SPIRO_RIGHT		']'
+
+    /* For a closed contour add an extra cp with a ty set to */
+#define SPIRO_END		'z'
+    /* For an open contour the first cp must have a ty set to*/
+#define SPIRO_OPEN_CONTOUR	'{'
+    /* For an open contour the last cp must have a ty set to */
+#define SPIRO_END_OPEN_CONTOUR	'}'
+
+/* The spiros array should indicate it's own end... So              */
+/* Open contours must have the ty field of the first cp set to '{'  */
+/*               and have the ty field of the last cp set to '}'    */
+/* Closed contours must have an extra cp at the end whose ty is 'z' */
+/*               the x&y values of this extra cp are ignored        */
+extern void TaggedSpiroCPsToBezier(spiro_cp *spiros,bezctx *bc);
+
+/* The first argument is an array of spiro control points.          */
+/* Open contours do not need to start with '{', nor to end with '}' */
+/* Close contours do not need to end with 'z'                       */
+extern void SpiroCPsToBezier(spiro_cp *spiros,int n,int isclosed,bezctx *bc);
+#endif
diff --git a/app/cornu/zmisc.h b/app/cornu/zmisc.h
new file mode 100644
index 0000000..ce3d4d3
--- /dev/null
+++ b/app/cornu/zmisc.h
@@ -0,0 +1,12 @@
+/**
+ * Misc portability and convenience macros.
+ **/
+
+#include <stdlib.h>
+
+#define zalloc malloc
+#define zrealloc realloc
+#define zfree free
+
+#define znew(type, n) (type *)zalloc(sizeof(type) * (n))
+#define zrenew(type, p, n) (type *)zrealloc((p), sizeof(type) * (n))