From d1a8285f818eb7e5c3d6a05709ea21a808490b8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 19 Mar 2018 19:55:58 +0100 Subject: New upstream version 5.1.0 --- app/cornu/CMakeLists.txt | 27 ++ app/cornu/bezctx.c | 48 ++ app/cornu/bezctx.h | 13 + app/cornu/bezctx_intf.h | 23 + app/cornu/bezctx_xtrkcad.c | 217 +++++++++ app/cornu/bezctx_xtrkcad.h | 4 + app/cornu/ppedit_gtk1.c | 930 +++++++++++++++++++++++++++++++++++ app/cornu/spiro.c | 1099 ++++++++++++++++++++++++++++++++++++++++++ app/cornu/spiro.h | 22 + app/cornu/spiroentrypoints.c | 60 +++ app/cornu/spiroentrypoints.h | 31 ++ app/cornu/zmisc.h | 12 + 12 files changed, 2486 insertions(+) create mode 100644 app/cornu/CMakeLists.txt create mode 100644 app/cornu/bezctx.c create mode 100644 app/cornu/bezctx.h create mode 100644 app/cornu/bezctx_intf.h create mode 100644 app/cornu/bezctx_xtrkcad.c create mode 100644 app/cornu/bezctx_xtrkcad.h create mode 100644 app/cornu/ppedit_gtk1.c create mode 100644 app/cornu/spiro.c create mode 100644 app/cornu/spiro.h create mode 100644 app/cornu/spiroentrypoints.c create mode 100644 app/cornu/spiroentrypoints.h create mode 100644 app/cornu/zmisc.h (limited to 'app/cornu') 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 +#include +#include +#include +#include +#include + +#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, "S"); + add_menuitem(menu, "Quit", (GtkSignalFunc)quit_func, data, ag, "Q"); + add_menuitem(menu, "Print", (GtkSignalFunc)print_func, data, ag, "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, "Z"); + p->redo_me = add_menuitem(menu, "Redo", (GtkSignalFunc)redo_func, data, ag, "Y"); + set_undo_state(p, NULL, NULL); + add_menuitem(menu, "Toggle Corner", (GtkSignalFunc)toggle_corner_func, data, ag, "T"); + add_menuitem(menu, "Delete Point", (GtkSignalFunc)delete_pt_func, data, ag, "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, "K"); + p->show_bg_me = add_menuitem(menu, "Hide BG", + (GtkSignalFunc)toggle_show_bg_func, + data, ag, "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 +#include +#include + +#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 +#include /* 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 + +#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)) -- cgit v1.2.3