diff options
Diffstat (limited to 'Build/source/texk/lcdf-typetools/lcdf-typetools-2.99/liblcdf/bezier.cc')
-rw-r--r-- | Build/source/texk/lcdf-typetools/lcdf-typetools-2.99/liblcdf/bezier.cc | 376 |
1 files changed, 376 insertions, 0 deletions
diff --git a/Build/source/texk/lcdf-typetools/lcdf-typetools-2.99/liblcdf/bezier.cc b/Build/source/texk/lcdf-typetools/lcdf-typetools-2.99/liblcdf/bezier.cc new file mode 100644 index 00000000000..95c8d8b9fb0 --- /dev/null +++ b/Build/source/texk/lcdf-typetools/lcdf-typetools-2.99/liblcdf/bezier.cc @@ -0,0 +1,376 @@ +// -*- related-file-name: "../include/lcdf/bezier.hh" -*- + +/* bezier.{cc,hh} -- cubic Bezier curves + * + * Copyright (c) 1998-2012 Eddie Kohler + * + * 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. + */ + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif +#include <lcdf/bezier.hh> + +// +// bounding box +// + +void +Bezier::make_bb() const throw () +{ + _bb = 0; + for (int i = 1; i < 4; i++) { + if (_p[i].x > bb_right_x()) + _bb = (_bb & ~0x03) | (i << 0); + else if (_p[i].x < bb_left_x()) + _bb = (_bb & ~0x0C) | (i << 2); + if (_p[i].y > bb_top_x()) + _bb = (_bb & ~0x30) | (i << 4); + else if (_p[i].y < bb_bottom_x()) + _bb = (_bb & ~0xC0) | (i << 6); + } +} + + +// +// is_flat, eval +// + +bool +Bezier::is_flat(double t) const throw () +{ + return (_p[2].on_segment(_p[0], _p[3], t) + && _p[1].on_segment(_p[0], _p[3], t)); +} + +static Point +eval_bezier(Point *b_in, int degree, double u) +{ + assert(degree < 4); + Point b[4]; + for (int i = 0; i <= degree; i++) + b[i] = b_in[i]; + + double m = 1.0 - u; + for (int i = 1; i <= degree; i++) + for (int j = 0; j <= degree - i; j++) + b[j] = b[j]*m + b[j+1]*u; + return b[0]; +} + +Point +Bezier::eval(double u) const throw () +{ + Bezier b = *this; + double m = 1.0 - u; + for (int i = 1; i < 4; i++) + for (int j = 0; j < 4 - i; j++) + b._p[j] = m * b._p[j] + u * b._p[j+1]; + return b._p[0]; +} + + +// +// halve +// + +void +Bezier::halve(Bezier &l, Bezier &r) const throw () +{ + Point half = Point::midpoint(_p[1], _p[2]); + l._p[0] = _p[0]; + l._p[1] = Point::midpoint(_p[0], _p[1]); + l._p[2] = Point::midpoint(l._p[1], half); + r._p[3] = _p[3]; + r._p[2] = Point::midpoint(_p[2], _p[3]); + r._p[1] = Point::midpoint(r._p[2], half); + r._p[0] = l._p[3] = Point::midpoint(l._p[2], r._p[1]); +} + + +// +// hit testing +// + +bool +Bezier::in_bb(const Point &p, double tolerance) const throw () +{ + ensure_bb(); + if (bb_right() + tolerance < p.x + || bb_left() - tolerance > p.x + || bb_top() + tolerance < p.y + || bb_bottom() - tolerance > p.y) + return false; + else + return true; +} + +double +Bezier::hit_recurse(const Point &p, double tolerance, double leftd, + double rightd, double leftt, double rightt) const throw () +{ + Bezier left, right; + double middled, resultt; + + if (is_flat(tolerance)) { + if (p.on_segment(_p[0], _p[3], tolerance)) + return (leftt + rightt) / 2; + else + return -1; + } + + if (leftd < tolerance * tolerance) + return leftt; + if (rightd < tolerance * tolerance) + return rightt; + + if (!in_bb(p, tolerance)) + return -1; + + halve(left, right); + middled = (right._p[0] - p).squared_length(); + resultt = left.hit_recurse + (p, tolerance, leftd, middled, leftt, (leftt + rightt) / 2); + if (resultt >= 0) + return resultt; + + return right.hit_recurse + (p, tolerance, middled, rightd, (leftt + rightt) / 2, rightt); +} + +bool +Bezier::hit(const Point &p, double tolerance) const throw () +{ + double leftd = (_p[0] - p).squared_length(); + double rightd = (_p[3] - p).squared_length(); + double resultt = hit_recurse(p, tolerance, leftd, rightd, 0, 1); + return resultt >= 0; +} + + +// +// segmentize to list of points +// +// uses recursive subdivision +// + +void +Bezier::segmentize(Vector<Point> &v, bool first) const +{ + if (is_flat(0.5)) { + if (first) + v.push_back(_p[0]); + v.push_back(_p[3]); + } else { + Bezier left, right; + halve(left, right); + left.segmentize(v, first); + right.segmentize(v, false); + } +} + + +// +// curve fitting +// +// code after Philip J. Schneider's algorithm described, with code, in the +// first Graphics Gems +// + +static void +chord_length_parameterize(const Point *d, int nd, Vector<double> &result) +{ + assert(result.size() == 0); + result.reserve(nd); + result.push_back(0); + for (int i = 1; i < nd; i++) + result.push_back(result.back() + Point::distance(d[i-1], d[i])); + double last_dist = result.back(); + for (int i = 1; i < nd; i++) + result[i] /= last_dist; +} + +static inline double +B0(double u) +{ + double m = 1.0 - u; + return m*m*m; +} + +static inline double +B1(double u) +{ + double m = 1.0 - u; + return 3*m*m*u; +} + +static inline double +B2(double u) +{ + double m = 1.0 - u; + return 3*m*u*u; +} + +static inline double +B3(double u) +{ + return u*u*u; +} + +static Bezier +generate_bezier(const Point *d, int nd, const Vector<double> ¶meters, + const Point &left_tangent, const Point &right_tangent) +{ + Point *a0 = new Point[nd]; + Point *a1 = new Point[nd]; + + for (int i = 0; i < nd; i++) { + a0[i] = left_tangent * B1(parameters[i]); + a1[i] = right_tangent * B2(parameters[i]); + } + + double c[2][2], x[2]; + c[0][0] = c[0][1] = c[1][0] = c[1][1] = x[0] = x[1] = 0.0; + + int last = nd - 1; + for (int i = 0; i < nd; i++) { + c[0][0] += Point::dot(a0[i], a0[i]); + c[0][1] += Point::dot(a0[i], a1[i]); + c[1][1] += Point::dot(a1[i], a1[i]); + + Point tmp = d[i] - (d[0] * (B0(parameters[i]) + B1(parameters[i])) + + d[last] * (B2(parameters[i]) + B3(parameters[i]))); + x[0] += Point::dot(a0[i], tmp); + x[1] += Point::dot(a1[i], tmp); + } + c[1][0] = c[0][1]; + + // compute determinants + double det_c0_c1 = c[0][0]*c[1][1] - c[1][0]*c[0][1]; + double det_c0_x = c[0][0]*x[1] - c[0][1]*x[0]; + double det_x_c1 = x[0]*c[1][1] - x[1]*c[0][1]; + + // finally, derive alpha values + if (det_c0_c1 == 0.0) + det_c0_c1 = c[0][0]*c[1][1] * 10e-12; + double alpha_l = det_x_c1 / det_c0_c1; + double alpha_r = det_c0_x / det_c0_c1; + + // if alpha negative, use the Wu/Barsky heuristic + if (alpha_l < 0.0 || alpha_r < 0.0) { + double distance = Point::distance(d[0], d[last]) / 3; + return Bezier(d[0], d[0] + left_tangent*distance, + d[last] + right_tangent*distance, d[last]); + } else + return Bezier(d[0], d[0] + left_tangent*alpha_l, + d[last] + right_tangent*alpha_r, d[last]); +} + +static double +newton_raphson_root_find(const Bezier &b, const Point &p, double u) +{ + const Point *b_pts = b.points(); + + Point b_det[3]; + for (int i = 0; i < 3; i++) + b_det[i] = (b_pts[i+1] - b_pts[i]) * 3; + + Point b_det_det[2]; + for (int i = 0; i < 2; i++) + b_det_det[i] = (b_det[i+1] - b_det[i]) * 2; + + Point b_u = b.eval(u); + Point b_det_u = eval_bezier(b_det, 2, u); + Point b_det_det_u = eval_bezier(b_det_det, 1, u); + + double numerator = Point::dot(b_u - p, b_det_u); + double denominator = Point::dot(b_det_u, b_det_u) + + Point::dot(b_u - p, b_det_det_u); + + return u - numerator/denominator; +} + +static void +reparameterize(const Point *d, int nd, Vector<double> ¶meters, + const Bezier &b) +{ + for (int i = 0; i < nd; i++) + parameters[i] = newton_raphson_root_find(b, d[i], parameters[i]); +} + +static double +compute_max_error(const Point *d, int nd, const Bezier &b, + const Vector<double> ¶meters, int *split_point) +{ + *split_point = nd/2; + double max_dist = 0.0; + for (int i = 1; i < nd - 1; i++) { + double dist = (b.eval(parameters[i]) - d[i]).squared_length(); + if (dist >= max_dist) { + max_dist = dist; + *split_point = i; + } + } + return max_dist; +} + +static void +fit0(const Point *d, int nd, Point left_tangent, Point right_tangent, + double error, Vector<Bezier> &result) +{ + // Use a heuristic for small regions (only two points) + if (nd == 2) { + double dist = Point::distance(d[0], d[1]) / 3; + result.push_back(Bezier(d[0], + d[0] + dist*left_tangent, + d[1] + dist*right_tangent, + d[1])); + return; + } + + // Parameterize points and attempt to fit curve + Vector<double> parameters; + chord_length_parameterize(d, nd, parameters); + Bezier b = generate_bezier(d, nd, parameters, left_tangent, right_tangent); + + // find max error + int split_point; + double max_error = compute_max_error(d, nd, b, parameters, &split_point); + if (max_error < error) { + result.push_back(b); + return; + } + + // if error not too large, try iteration and reparameterization + if (max_error < error*error) + for (int i = 0; i < 4; i++) { + reparameterize(d, nd, parameters, b); + b = generate_bezier(d, nd, parameters, left_tangent, right_tangent); + max_error = compute_max_error(d, nd, b, parameters, &split_point); + if (max_error < error) { + result.push_back(b); + return; + } + } + + // fitting failed -- split at max error point and fit again + Point center_tangent = ((d[split_point-1] - d[split_point+1])/2).normal(); + fit0(d, split_point+1, left_tangent, center_tangent, error, result); + fit0(d+split_point, nd-split_point, -center_tangent, right_tangent, error, result); +} + +void +Bezier::fit(const Vector<Point> &points, double error, Vector<Bezier> &result) +{ + int npoints = points.size(); + Point left_tangent = (points[1] - points[0]).normal(); + Point right_tangent = (points[npoints-2] - points[npoints-1]).normal(); + fit0(&points[0], npoints, left_tangent, right_tangent, error, result); +} |