X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=lily%2Fbezier.cc;h=1eaf041b886e8c782efdd802301b0f2ebf78a60a;hb=6b876e65197de4fcc073ebc19a5ef744db92c1c9;hp=657cfea0d720e342f346f1fb314ea8cd7abbc6ad;hpb=234697355e387e3349018823d259c35f4ee49713;p=lilypond.git diff --git a/lily/bezier.cc b/lily/bezier.cc index 657cfea0d7..1eaf041b88 100644 --- a/lily/bezier.cc +++ b/lily/bezier.cc @@ -3,123 +3,265 @@ source file of the GNU LilyPond music typesetter - (c) 1998 Jan Nieuwenhuizen + (c) 1998--2007 Jan Nieuwenhuizen */ -#include -#include "offset.hh" #include "bezier.hh" -#include "dimen.hh" -#include "paper-def.hh" +#include "warn.hh" +#include "libc-extension.hh" -Bezier::Bezier (int steps_i) +Real binomial_coefficient_3[] = { + 1, 3, 3, 1 +}; + +void +scale (vector *array, Real x, Real y) { - steps_i_ = steps_i; - curve_ = new Offset [steps_i_ + 1]; + for (vsize i = 0; i < array->size (); i++) + { + (*array)[i][X_AXIS] = x * (*array)[i][X_AXIS]; + (*array)[i][Y_AXIS] = y * (*array)[i][Y_AXIS]; + } } -Bezier::~Bezier () +void +rotate (vector *array, Real phi) { - delete[] curve_; + Offset rot (complex_exp (Offset (0, phi))); + for (vsize i = 0; i < array->size (); i++) + (*array)[i] = complex_multiply (rot, (*array)[i]); } -//from GNU gs3.33: ega.c void -Bezier::calc (Offset control[4]) -{ - Real dt = 1.0 / steps_i_; - Real cx = 3.0 * (control[1].x() - control[0].x()); - Real bx = 3.0 * (control[2].x() - control[1].x()) - cx; - Real ax = control[3].x() - (control[0].x() + cx + bx); - Real cy = 3.0 * (control[1].y () - control[0].y ()); - Real by = 3.0 * (control[2].y () - control[1].y ()) - cy; - Real ay = control[3].y () - (control[0].y () + cy + by); - Real t = 0.0; - int i = 0; - while ( t <= 1.0 ) - { - curve_[i].x() = ((ax * t + bx) * t + cx) * t + control[0].x(); - curve_[i++].y () = ((ay * t + by) * t + cy) * t + control[0].y (); - t += dt; +translate (vector *array, Offset o) +{ + for (vsize i = 0; i < array->size (); i++) + (*array)[i] += o; +} + +/* + Formula of the bezier 3-spline + + sum_{j = 0}^3 (3 over j) z_j (1-t)^ (3-j) t^j + + + A is the axis of X coordinate. +*/ + +Real +Bezier::get_other_coordinate (Axis a, Real x) const +{ + Axis other = Axis ((a +1) % NO_AXES); + vector ts = solve_point (a, x); + + if (ts.size () == 0) + { + programming_error ("no solution found for Bezier intersection"); + return 0.0; } + +#ifdef PARANOID + Offset c = curve_point (ts[0]); + if (fabs (c[a] - x) > 1e-8) + programming_error ("bezier intersection not correct?"); +#endif + + return curve_coordinate (ts[0], other); } Real -Bezier::y (Real x) +Bezier::curve_coordinate (Real t, Axis a) const { - if (x <= curve_[0].x()) - return curve_[0].y (); - for (int i = 1; i < steps_i_; i++ ) + Real tj = 1; + Real one_min_tj[4]; + one_min_tj[0] = 1; + for (int i = 1; i < 4; i++) + one_min_tj[i] = one_min_tj[i - 1] * (1 - t); + + Real r = 0.0; + for (int j = 0; j < 4; j++) + { + r += control_[j][a] * binomial_coefficient_3[j] + * tj * one_min_tj[3 - j]; + + tj *= t; + } + + return r; +} + +Offset +Bezier::curve_point (Real t) const +{ + Real tj = 1; + Real one_min_tj[4]; + one_min_tj[0] = 1; + for (int i = 1; i < 4; i++) + one_min_tj[i] = one_min_tj[i - 1] * (1 - t); + + Offset o; + for (int j = 0; j < 4; j++) { - if (x < curve_[i].x()) - { - Real lin = (x - curve_[i-1].x()) / (curve_[i].x() - curve_[i-1].x()); - return curve_[i-1].y () + lin * (curve_[i].y () - curve_[i-1].y ()); - } + o += control_[j] * binomial_coefficient_3[j] + * tj * one_min_tj[3 - j]; + + tj *= t; } - return curve_[steps_i_-1].y (); + +#ifdef PARANOID + assert (fabs (o[X_AXIS] - polynomial (X_AXIS).eval (t)) < 1e-8); + assert (fabs (o[Y_AXIS] - polynomial (Y_AXIS).eval (t)) < 1e-8); +#endif + + return o; } +/* + Cache binom (3, j) t^j (1-t)^{3-j} +*/ +struct Polynomial_cache { + Polynomial terms_[4]; + Polynomial_cache () + { + for (int j = 0; j <= 3; j++) + terms_[j] + = binomial_coefficient_3[j] + * Polynomial::power (j, Polynomial (0, 1)) + * Polynomial::power (3 - j, Polynomial (1, -1)); + } +}; + +static Polynomial_cache poly_cache; -Bezier_bow::Bezier_bow (Paper_def* paper_l) - : Bezier(10) +Polynomial +Bezier::polynomial (Axis a) const { - paper_l_ = paper_l; + Polynomial p (0.0); + Polynomial q; + for (int j = 0; j <= 3; j++) + { + q = poly_cache.terms_[j]; + q *= control_[j][a]; + p += q; + } + + return p; } +/** + Remove all numbers outside [0, 1] from SOL +*/ +vector +filter_solutions (vector sol) +{ + for (vsize i = sol.size (); i--;) + if (sol[i] < 0 || sol[i] > 1) + sol.erase (sol.begin () + i); + return sol; +} + +/** + find t such that derivative is proportional to DERIV +*/ +vector +Bezier::solve_derivative (Offset deriv) const +{ + Polynomial xp = polynomial (X_AXIS); + Polynomial yp = polynomial (Y_AXIS); + xp.differentiate (); + yp.differentiate (); + + Polynomial combine = xp * deriv[Y_AXIS] - yp * deriv [X_AXIS]; -/* - from feta-sleur.ly + return filter_solutions (combine.solve ()); +} - slurheightlimit#:=staffsize#/2; - sluralpha:=slurheightlimit#*pi/2; - % slurratio:=1/3; - slurratio:=0.3333; - slurbeta:=3/4*pi*slurratio/sluralpha; +/* + Find t such that curve_point (t)[AX] == COORDINATE +*/ +vector +Bezier::solve_point (Axis ax, Real coordinate) const +{ + Polynomial p (polynomial (ax)); + p.coefs_[0] -= coordinate; - b#:=length(dx#,dy#); - % ugh: huh? 2/5 - indent#:=2/5*sluralpha*atan(slurbeta*b#); - height:=(indent+h)*d; - z1=(0,0); - z2=(b,0); - z3=(indent,height); - z4=(b-indent,height); + vector sol (p.solve ()); + return filter_solutions (sol); +} - boogje:=boogje rotated angle(dxs,dys); +/** + Compute the bounding box dimensions in direction of A. */ +Interval +Bezier::extent (Axis a) const +{ + int o = (a + 1)%NO_AXES; + Offset d; + d[Axis (o)] = 1.0; + Interval iv; + vector sols (solve_derivative (d)); + sols.push_back (1.0); + sols.push_back (0.0); + for (vsize i = sols.size (); i--;) + { + Offset o (curve_point (sols[i])); + iv.unite (Interval (o[a], o[a])); + } + return iv; +} + +Interval +Bezier::control_point_extent (Axis a) const +{ + Interval ext; + for (int i = CONTROL_COUNT; i--;) + ext.add_point (control_[i][a]); + + return ext; +} + +/** + Flip around axis A +*/ void -Bezier_bow::calc (Real dx, Real dy, Real h, Real d) -{ - // ugh - Real pi = M_PI; - // ugh - Real staffsize_f = paper_l_->get_var ("barsize"); - Real height_limit = staffsize_f; - Real alpha = height_limit * pi / 2.0; - Real ratio = 1.0/3.0; - Real beta = 3.0/4.0 * pi * ratio/alpha; - - Real b = sqrt (dx * dx + dy * dy); - Real indent = 2.0/5.0 * alpha * atan (beta * b); - // ugh, ugly height hack, see lily-ps-defs.tex - Real height = (indent + h) * d; - - Offset control[4]; - control[0] = Offset(0, 0); - control[1] = Offset(indent, height); - control[2] = Offset(b - indent, height); - control[3] = Offset( b, 0 ); - - Real phi = dx ? atan (dy/dx) : sign (dy) * pi / 2.0; - Real sphi = sin (phi); - Real cphi = cos (phi); - for (int i = 1; i < 4; i++) +Bezier::scale (Real x, Real y) +{ + for (int i = CONTROL_COUNT; i--;) { - control[i].x() = cphi * control[i].x() - sphi * control[i].y (); - control[i].y () = sphi * control[i].x() + cphi * control[i].y (); + control_[i][X_AXIS] = x * control_[i][X_AXIS]; + control_[i][Y_AXIS] = y * control_[i][Y_AXIS]; } - Bezier::calc (control); } +void +Bezier::rotate (Real phi) +{ + Offset rot (complex_exp (Offset (0, phi))); + for (int i = 0; i < CONTROL_COUNT; i++) + control_[i] = complex_multiply (rot, control_[i]); +} + +void +Bezier::translate (Offset o) +{ + for (int i = 0; i < CONTROL_COUNT; i++) + control_[i] += o; +} + +void +Bezier::assert_sanity () const +{ + for (int i = 0; i < CONTROL_COUNT; i++) + assert (!isnan (control_[i].length ()) + && !isinf (control_[i].length ())); +} + +void +Bezier::reverse () +{ + Bezier b2; + for (int i = 0; i < CONTROL_COUNT; i++) + b2.control_[CONTROL_COUNT - i - 1] = control_[i]; + *this = b2; +}