![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: fe9253d)
raw | patch | inline | side by side (parent: fe9253d)
| author | johanengelen <johanengelen@users.sourceforge.net> | |
| Mon, 1 Sep 2008 19:29:30 +0000 (19:29 +0000) | ||
| committer | johanengelen <johanengelen@users.sourceforge.net> | |
| Mon, 1 Sep 2008 19:29:30 +0000 (19:29 +0000) |
21 files changed:
index 48b990445ed3c0296378237811ed247750aade2b..334c23feaa50fcaea1010e8d1f3b4ffd699aee16 100644 (file)
+
+
+
#include <2geom/basic-intersection.h>
+#include <2geom/sbasis-to-bezier.h>
#include <2geom/exception.h>
+
+
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
-
+
unsigned intersect_steps = 0;
}
void split(double t, OldBezier &a, OldBezier &b) const;
Point operator()(double t) const;
-
+
~OldBezier() {}
- void bounds(double &minax, double &maxax,
+ void bounds(double &minax, double &maxax,
double &minay, double &maxay) {
// Compute bounding box for a
minax = p[0][X]; // These are the most likely to be extremal
}
}
-
+
};
static std::vector<std::pair<double, double> >
find_intersections( OldBezier a, OldBezier b);
-static std::vector<std::pair<double, double> >
+static std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A);
std::vector<std::pair<double, double> >
-find_intersections( vector<Geom::Point> const & A,
+find_intersections( vector<Geom::Point> const & A,
vector<Geom::Point> const & B) {
OldBezier a, b;
a.p = A;
return find_intersections(a,b);
}
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &/*Sb*/) {
THROW_NOTIMPLEMENTED();
}
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(D2<SBasis> const & A) {
OldBezier Sb;
- Sb.p = sbasis_to_bezier(A);
+ sbasis_to_bezier(Sb.p, A);
return find_self_intersections(Sb, A);
}
-static std::vector<std::pair<double, double> >
+static std::vector<std::pair<double, double> >
find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
-
+
vector<double> dr = roots(derivative(A[X]));
{
vector<double> dyr = roots(derivative(A[Y]));
// We want to be sure that we have no empty segments
sort(dr.begin(), dr.end());
unique(dr.begin(), dr.end());
-
+
std::vector<std::pair<double, double> > all_si;
-
+
vector<OldBezier> pieces;
{
OldBezier in = Sb, l, r;
}
for(unsigned i = 0; i < dr.size()-1; i++) {
for(unsigned j = i+1; j < dr.size()-1; j++) {
- std::vector<std::pair<double, double> > section =
+ std::vector<std::pair<double, double> > section =
find_intersections( pieces[i], pieces[j]);
for(unsigned k = 0; k < section.size(); k++) {
double l = section[k].first;
}
}
}
-
+
// Because i is in order, all_si should be roughly already in order?
//sort(all_si.begin(), all_si.end());
//unique(all_si.begin(), all_si.end());
-
+
return all_si;
}
/* The value of 1.0 / (1L<<14) is enough for most applications */
const double INV_EPS = (1L<<14);
-
+
/*
* split the curve at the midpoint, returning an array with the two parts
* Temporary storage is minimized by using part of the storage for the result
std::copy(p.begin(), p.end(), Vtemp[0]);
/* Triangle computation */
- for (unsigned i = 1; i < sz; i++) {
+ for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
}
}
-
+
left.p.resize(sz);
right.p.resize(sz);
for (unsigned j = 0; j < sz; j++)
right.p[j] = Vtemp[sz-1-j][j];
}
+#if 0
/*
* split the curve at the midpoint, returning an array with the two parts
* Temporary storage is minimized by using part of the storage for the result
std::copy(p.begin(), p.end(), Vtemp[0]);
/* Triangle computation */
- for (unsigned i = 1; i < sz; i++) {
+ for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
}
}
return Vtemp[sz-1][0];
}
+#endif
+
+// suggested by Sederberg.
+Point OldBezier::operator()(double t) const {
+ int n = p.size()-1;
+ double u, bc, tn, tmp;
+ int i;
+ Point r;
+ for(int dim = 0; dim < 2; dim++) {
+ u = 1.0 - t;
+ bc = 1;
+ tn = 1;
+ tmp = p[0][dim]*u;
+ for(i=1; i<n; i++){
+ tn = tn*t;
+ bc = bc*(n-i+1)/i;
+ tmp = (tmp + tn*bc*p[i][dim])*u;
+ }
+ r[dim] = (tmp + tn*t*p[n][dim]);
+ }
+ return r;
+}
/*
* Several observations:
* First, it is cheaper to compute the bounding box of the second curve
* and test its bounding box for interference than to use a more direct
- * approach of comparing all control points of the second curve with
+ * approach of comparing all control points of the second curve with
* the various edges of the bounding box of the first curve to test
* for interference.
* Second, after a few subdivisions it is highly probable that two corners
- * of the bounding box of a given Bezier curve are the first and last
+ * of the bounding box of a given Bezier curve are the first and last
* control point. Once this happens once, it happens for all subsequent
* subcurves. It might be worth putting in a test and then short-circuit
* code for further subdivision levels.
return not( ( minax > maxbx ) || ( minay > maxby )
|| ( minbx > maxax ) || ( minby > maxay ) );
}
-
-/*
- * Recursively intersect two curves keeping track of their real parameters
+
+/*
+ * Recursively intersect two curves keeping track of their real parameters
* and depths of intersection.
* The results are returned in a 2-D array of doubles indicating the parameters
* for which intersections are found. The parameters are in the order the
* intersections were found, which is probably not in sorted order.
- * When an intersection is found, the parameter value for each of the two
+ * When an intersection is found, the parameter value for each of the two
* is stored in the index elements array, and the index is incremented.
- *
+ *
* If either of the curves has subdivisions left before it is straight
* (depth > 0)
* that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
* the depth(s) is (are) decremented, and the parameter value(s) corresponding
* to the midpoints(s) is (are) computed.
- * Then each of the subcurves of one curve is intersected with each of the
+ * Then each of the subcurves of one curve is intersected with each of the
* subcurves of the other curve, first by testing the bounding boxes for
* interference. If there is any bounding box interference, the corresponding
* subcurves are recursively intersected.
- *
+ *
* If neither curve has subdivisions left, the line segments from the first
- * to last control point of each segment are intersected. (Actually the
+ * to last control point of each segment are intersected. (Actually the
* only the parameter value corresponding to the intersection point is found).
*
* The apriori flatness test is probably more efficient than testing at each
* level of recursion, although a test after three or four levels would
- * probably be worthwhile, since many curves become flat faster than their
+ * probably be worthwhile, since many curves become flat faster than their
* asymptotic rate for the first few levels of recursion.
*
* The bounding box test fails much more frequently than it succeeds, providing
* substantial pruning of the search space.
*
- * Each (sub)curve is subdivided only once, hence it is not possible that for
+ * Each (sub)curve is subdivided only once, hence it is not possible that for
* one final line intersection test the subdivision was at one level, while
* for another final line intersection test the subdivision (of the same curve)
* was at another. Since the line segments share endpoints, the intersection
}
inline double log4( double x ) { return log(x)/log(4.); }
-
+
/*
* Wang's theorem is used to estimate the level of subdivision required,
* but only if the bounding boxes interfere at the top level.
double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
double l0 = std::max(la1, la2);
unsigned ra;
- if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
+ if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
ra = 0;
else
ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
OldBezier &A;
OldBezier &B;
};
-
+
static int
intersect_polish_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
const double x0 = gsl_vector_get (x, 0);
const double x1 = gsl_vector_get (x, 1);
-
- Geom::Point dx = ((struct rparams *) params)->A(x0) -
+
+ Geom::Point dx = ((struct rparams *) params)->A(x0) -
((struct rparams *) params)->B(x1);
-
+
gsl_vector_set (f, 0, dx[0]);
gsl_vector_set (f, 1, dx[1]);
-
+
return GSL_SUCCESS;
}
+typedef union dbl_64{
+ long long i64;
+ double d64;
+};
+
+static double EpsilonBy(double value, int eps)
+{
+ dbl_64 s;
+ s.d64 = value;
+ s.i64 += eps;
+ return s.d64;
+}
+
+
static void intersect_polish_root (OldBezier &A, double &s,
OldBezier &B, double &t) {
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *sol;
-
+
int status;
size_t iter = 0;
-
+
const size_t n = 2;
struct rparams p = {A, B};
gsl_multiroot_function f = {&intersect_polish_f, n, &p};
-
+
double x_init[2] = {s, t};
gsl_vector *x = gsl_vector_alloc (n);
-
+
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
-
+
T = gsl_multiroot_fsolver_hybrids;
sol = gsl_multiroot_fsolver_alloc (T, 2);
gsl_multiroot_fsolver_set (sol, &f, x);
-
+
do
{
iter++;
status = gsl_multiroot_fsolver_iterate (sol);
-
+
if (status) /* check if solver is stuck */
break;
-
+
status =
gsl_multiroot_test_residual (sol->f, 1e-12);
}
while (status == GSL_CONTINUE && iter < 1000);
-
+
s = gsl_vector_get (sol->x, 0);
t = gsl_vector_get (sol->x, 1);
-
+
gsl_multiroot_fsolver_free (sol);
gsl_vector_free (x);
+
+ // This code does a neighbourhood search for minor improvements.
+ double best_v = L1(A(s) - B(t));
+ //std::cout << "------\n" << best_v << std::endl;
+ Point best(s,t);
+ while (true) {
+ Point trial = best;
+ double trial_v = best_v;
+ for(int nsi = -1; nsi < 2; nsi++) {
+ for(int nti = -1; nti < 2; nti++) {
+ Point n(EpsilonBy(best[0], nsi),
+ EpsilonBy(best[1], nti));
+ double c = L1(A(n[0]) - B(n[1]));
+ //std::cout << c << "; ";
+ if (c < trial_v) {
+ trial = n;
+ trial_v = c;
+ }
+ }
+ }
+ if(trial == best) {
+ //std::cout << "\n" << s << " -> " << s - best[0] << std::endl;
+ //std::cout << t << " -> " << t - best[1] << std::endl;
+ //std::cout << best_v << std::endl;
+ s = best[0];
+ t = best[1];
+ return;
+ } else {
+ best = trial;
+ best_v = trial_v;
+ }
+ }
}
@@ -432,8 +507,8 @@ std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezi
std::vector<std::pair<double, double> > parameters;
if( intersect_BB( a, b ) )
{
- recursively_intersect( a, 0., 1., wangs_theorem(a),
- b, 0., 1., wangs_theorem(b),
+ recursively_intersect( a, 0., 1., wangs_theorem(a),
+ b, 0., 1., wangs_theorem(b),
parameters);
}
for(unsigned i = 0; i < parameters.size(); i++)
index 0090b0305e72e0675401000bfb1f2f8cbfa08891..d464891f92be972d080b899b62d95919ad6e1ab5 100644 (file)
+
+
+#include <2geom/point.h>
#include <2geom/sbasis.h>
-#include <2geom/bezier-to-sbasis.h>
-#include <2geom/sbasis-to-bezier.h>
#include <2geom/d2.h>
+#include <vector>
+#include <utility>
+
+
namespace Geom {
-std::vector<std::pair<double, double> >
-find_intersections( D2<SBasis> const & A,
+std::vector<std::pair<double, double> >
+find_intersections( D2<SBasis> const & A,
D2<SBasis> const & B);
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(D2<SBasis> const & A);
// Bezier form
-std::vector<std::pair<double, double> >
-find_intersections( std::vector<Point> const & A,
+std::vector<std::pair<double, double> >
+find_intersections( std::vector<Point> const & A,
std::vector<Point> const & B);
-std::vector<std::pair<double, double> >
+std::vector<std::pair<double, double> >
find_self_intersections(std::vector<Point> const & A);
+
+/*
+ * find_intersection
+ *
+ * input: A, B - set of control points of two Bezier curve
+ * input: precision - required precision of computation
+ * output: xs - set of pairs of parameter values
+ * at which crossing happens
+ *
+ * This routine is based on the Bezier Clipping Algorithm,
+ * see: Sederberg - Computer Aided Geometric Design
+ */
+void find_intersections (std::vector< std::pair<double, double> > & xs,
+ std::vector<Point> const& A,
+ std::vector<Point> const& B,
+ double precision = 1e-5);
+
};
/*
index 8b421f2e7e872299ce7c2b246f11a243f51f6892..d5279a5702f0980a0e5cdf44c8fc0b45aa5e4a25 100644 (file)
#define _BEZIER_TO_SBASIS
#include <2geom/coord.h>
-
-#include <2geom/d2.h>
#include <2geom/point.h>
+#include <2geom/d2.h>
+#include <2geom/sbasis-to-bezier.h>
-namespace Geom{
+namespace Geom
+{
+#if 0
inline SBasis bezier_to_sbasis(Coord const *handles, unsigned order) {
if(order == 0)
return Linear(handles[0]);
template <typename T>
-inline D2<SBasis> handles_to_sbasis(T const &handles, unsigned order) {
+inline D2<SBasis> handles_to_sbasis(T const &handles, unsigned order)
+{
double v[2][order+1];
for(unsigned i = 0; i <= order; i++)
for(unsigned j = 0; j < 2; j++)
return D2<SBasis>(bezier_to_sbasis(v[0], order),
bezier_to_sbasis(v[1], order));
}
+#endif
+
+
+template <typename T>
+inline
+D2<SBasis> handles_to_sbasis(T const& handles, unsigned order)
+{
+ D2<SBasis> sbc;
+ size_t sz = order + 1;
+ std::vector<Point> v;
+ v.reserve(sz);
+ for (size_t i = 0; i < sz; ++i)
+ v.push_back(handles[i]);
+ bezier_to_sbasis(sbc, v);
+ return sbc;
+}
+
+
+} // end namespace Geom
+
+
+
-};
#endif
/*
Local Variables:
diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 79b15cd0364259270c4fcc2c1da3c78e5e983992..6e11ad58b3f94dd49b5d9a717eaa7b269d72ef6b 100644 (file)
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
#include <2geom/coord.h>
#include <valarray>
#include <2geom/isnan.h>
-#include <2geom/bezier-to-sbasis.h>
#include <2geom/d2.h>
#include <2geom/solver.h>
#include <boost/optional/optional.hpp>
std::copy(v, v+order+1, vtemp[0]);
/* Triangle computation */
- for (unsigned i = 1; i <= order; i++) {
+ for (unsigned i = 1; i <= order; i++) {
for (unsigned j = 0; j <= order - i; j++) {
vtemp[i][j] = lerp(t, vtemp[i-1][j], vtemp[i-1][j+1]);
}
public:
unsigned int order() const { return c_.size()-1;}
unsigned int size() const { return c_.size();}
-
- Bezier() :c_(0., 32) {}
+
+ Bezier() {}
Bezier(const Bezier& b) :c_(b.c_) {}
Bezier &operator=(Bezier const &other) {
if ( c_.size() != other.c_.size() ) {
c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
}
+ void resize (unsigned int n, Coord v = 0)
+ {
+ c_.resize (n, v);
+ }
+
+ void clear()
+ {
+ c_.resize(0);
+ }
+
inline unsigned degree() const { return order(); }
//IMPL: FragmentConcept
typedef Coord output_type;
- inline bool isZero() const {
+ inline bool isZero() const {
for(unsigned i = 0; i <= order(); i++) {
if(c_[i] != 0) return false;
}
inline Coord at0() const { return c_[0]; }
inline Coord at1() const { return c_[order()]; }
- inline Coord valueAt(double t) const {
- return subdivideArr(t, &c_[0], NULL, NULL, order());
+ inline Coord valueAt(double t) const {
+ return subdivideArr(t, &c_[0], NULL, NULL, order());
}
inline Coord operator()(double t) const { return valueAt(t); }
- inline SBasis toSBasis() const {
- return bezier_to_sbasis(&c_[0], order());
- }
+ SBasis toSBasis() const;
+// inline SBasis toSBasis() const {
+// SBasis sb;
+// bezier_to_sbasis(sb, (*this));
+// return sb;
+// //return bezier_to_sbasis(&c_[0], order());
+// }
//Only mutator
inline Coord &operator[](unsigned ix) { return c_[ix]; }
inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
-
+
/* This is inelegant, as it uses several extra stores. I think there might be a way to
* evaluate roughly in situ. */
-
+
std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
std::vector<Coord> val_n_der;
Coord d_[order()+1];
val_n_der.resize(n_derivs);
return val_n_der;
}
-
+
std::pair<Bezier, Bezier > subdivide(Coord t) const {
Bezier a(Bezier::Order(*this)), b(Bezier::Order(*this));
subdivideArr(t, &c_[0], &a.c_[0], &b.c_[0], order());
}
};
+
+void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
+
+
+inline
+SBasis Bezier::toSBasis() const {
+ SBasis sb;
+ bezier_to_sbasis(sb, (*this));
+ return sb;
+ //return bezier_to_sbasis(&c_[0], order());
+}
+
//TODO: implement others
inline Bezier operator+(const Bezier & a, double v) {
Bezier result = Bezier(Bezier::Order(a));
//if(a.order() == 1) return Bezier(0.0);
if(a.order() == 1) return Bezier(a.c_[1]-a.c_[0]);
Bezier der(Bezier::Order(a.order()-1));
-
+
for(unsigned i = 0; i < a.order(); i++) {
der.c_[i] = a.order()*(a.c_[i+1] - a.c_[i]);
}
inline Bezier integral(const Bezier & a) {
Bezier inte(Bezier::Order(a.order()+1));
-
+
inte[0] = 0;
for(unsigned i = 0; i < inte.order(); i++) {
inte[i+1] = inte[i] + a[i]/(inte.order());
diff --git a/src/2geom/crossing.cpp b/src/2geom/crossing.cpp
index a4d4f6ab1ac2e95afef6405a4ca022671ed74c90..d717a4ed560c6dbe4b3b17e37caf3ea3288efa84 100644 (file)
--- a/src/2geom/crossing.cpp
+++ b/src/2geom/crossing.cpp
}
}
}
-
+/*
CrossingGraph create_crossing_graph(std::vector<Path> const &p, Crossings const &crs) {
std::vector<Point> locs;
CrossingGraph ret;
}
return ret;
-
+ */
/* Various incoherent code bits
// list of sets of edges, each set corresponding to those emanating from the path
CrossingGraph ret;
}
}
*/
-}
+//}
void merge_crossings(Crossings &a, Crossings &b, unsigned i) {
Crossings n;
diff --git a/src/2geom/crossing.h b/src/2geom/crossing.h
index 546e33ebd41df747051a91379d3c11f83eee8d1e..b16c7e46d1a10a6e2b6d4634520c7433f0bddffb 100644 (file)
--- a/src/2geom/crossing.h
+++ b/src/2geom/crossing.h
};
+/*
struct Edge {
unsigned node, path;
double time;
}
};
-typedef std::vector<Crossing> Crossings;
typedef std::vector<CrossingNode> CrossingGraph;
class Path;
CrossingGraph create_crossing_graph(std::vector<Path> const &p, Crossings const &crs);
+*/
/*inline bool are_near(Crossing a, Crossing b) {
return are_near(a.ta, b.ta) && are_near(a.tb, b.tb);
}
};
+typedef std::vector<Crossing> Crossings;
typedef std::vector<Crossings> CrossingSet;
diff --git a/src/2geom/forward.h b/src/2geom/forward.h
index 9bd5dedbe7cdc18468a6b16792fde5e832399794..d4ac1af1ed514326add0c13e4234dcfad0658067 100644 (file)
--- a/src/2geom/forward.h
+++ b/src/2geom/forward.h
*
* Authors:
* Johan Engelen <goejendaagh@zonnet.nl>
- *
+ *
* Copyright 2008 authors
*
* This library is free software; you can redistribute it and/or
namespace Geom {
+class Bezier;
template <unsigned> class BezierCurve;
template<> class BezierCurve<0>;
typedef BezierCurve<2> QuadraticBezier;
diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 74e03aec9924c63234cbb01cc22b290a6eccad4d..2d9f4abd85ce01582edce7d679f7a34bf862ab31 100644 (file)
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
return contains(val._b[0]) || contains(val._b[1]) || val.contains(*this);
}
- inline bool operator==(Interval other) { return _b[0] == other._b[0] && _b[1] == other._b[1]; }
- inline bool operator!=(Interval other) { return _b[0] != other._b[0] || _b[1] != other._b[1]; }
+ inline bool operator==(Interval other) const { return _b[0] == other._b[0] && _b[1] == other._b[1]; }
+ inline bool operator!=(Interval other) const { return _b[0] != other._b[0] || _b[1] != other._b[1]; }
//IMPL: OffsetableConcept
//TODO: rename output_type to something else in the concept
index de8f24760aa91e8bad823aefd440fe5204038f90..564663cf78a1e188ac96fef6c6db4b9332d848bb 100644 (file)
namespace Geom { namespace NL {
+/*
+ * A model is an abstraction for an expression dependent from a parameter where
+ * the coefficients of this expression are the unknowns of the fitting problem.
+ * For a ceratain number of parameter values we know the related values
+ * the expression evaluates to: from each parameter value we get a row of
+ * the matrix of the fitting problem, from each expression value we get
+ * the related constant term.
+ * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get
+ * the equation a + b + c = 0, in this example the constant term is always
+ * the same for each parameter value.
+ *
+ * A model is required to implement 3 methods:
+ *
+ * - size : returns the number of unknown coefficients that appear in
+ * the expression of the fitting problem;
+ * - feed : its input is a parameter value and the related expression value,
+ * it generates a matrix row and a new entry of the constant vector
+ * of the fitting problem;
+ * - instance : it has an input parameter represented by the raw vector
+ * solution of the fitting problem and an output parameter
+ * of type InstanceType that return a specific object that is
+ * generated using the fitting problem solution, in the example
+ * above the object could be a Poly type.
+ */
/*
* completely unknown models must inherit from this template class;
index d589d86e39bcdd8a882ee8874c6a11e8f0dbcfe1..47c588cb22fcec488fdb91a336dbecc7a973eb41 100644 (file)
#include <vector>
+/*
+ * The least_square_fitter class represents a tool for solving a fitting
+ * problem with respect to a given model that represents an expression
+ * dependent from a parameter where the coefficients of this expression
+ * are the unknowns of the fitting problem.
+ * The minimizing solution is found by computing the pseudo-inverse
+ * of the problem matrix
+ */
+
+
namespace Geom { namespace NL {
namespace detail {
-
template< typename ModelT>
class lsf_base
{
diff --git a/src/2geom/path.h b/src/2geom/path.h
index 9dfc929754b5dfa3789f723a3701b289d92c0700..04b77862c25c715327b8e214cff6ec93964f1cfa 100644 (file)
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
size_type max_size() const { return get_curves().max_size()-1; }
- bool empty() const { return get_curves().size() == 1; }
+ bool empty() const { return (get_curves().size() == 1); }
bool closed() const { return closed_; }
void close(bool closed=true) { closed_ = closed; }
diff --git a/src/2geom/rect.h b/src/2geom/rect.h
index fa79201e80ad5d67c895e62a60202b4b3472a4df..21cdcf849a9cafc8d1dae73e6031c391df1caed5 100644 (file)
--- a/src/2geom/rect.h
+++ b/src/2geom/rect.h
inline double maxExtent() const { return std::max(f[X].extent(), f[Y].extent()); }
inline bool isEmpty() const {
- return f[X].isEmpty() && f[Y].isEmpty();
+ return f[X].isEmpty() || f[Y].isEmpty();
}
inline bool intersects(Rect const &r) const {
return f[X].intersects(r[X]) && f[Y].intersects(r[Y]);
index 861baf76d7637815ccff468811f5b17e656eff37..09a84050bbd8cb622f0a0e4b0ce28c74a5b13c39 100644 (file)
* multi-roots using bernstein method, one approach would be:
sort c
take median and find roots of that
- whenever a segment lies entirely on one side of the median,
+ whenever a segment lies entirely on one side of the median,
find the median of the half and recurse.
in essence we are implementing quicksort on a continuous function
going from f(a) up to C with slope M takes at least time (C-f(a))/M
From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
Do the same for b: compute some b' such that there are no roots in (b',b].
- -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
+ -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
making things tricky and unpleasant...
*/
double fa,
double b,
double fb){
-
+
if (f.size()==0){
int idx;
idx=upper_level(levels,0,vtol);
}
return;
}
-////usefull?
+////usefull?
// if (f.size()==1){
// int idxa=upper_level(levels,fa);
// int idxb=upper_level(levels,fb);
}
return;
}
-
+
int idxa=upper_level(levels,fa,vtol);
int idxb=upper_level(levels,fb,vtol);
if (bs.max()>0 && idxb>0)
tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
}
-
+
double t0,t1;
- t0=std::min(ta_hi,ta_lo);
+ t0=std::min(ta_hi,ta_lo);
t1=std::max(tb_hi,tb_lo);
//hum, rounding errors frighten me! so I add this +tol...
if (t0>t1+htol) return;//no root here.
std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
SBasis df=derivative(f);
- multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
+ multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
return(roots);
}
b = xk*b + p[j];
err = fabs(b) + abx*err;
}
-
+
err *= 1e-7; // magic epsilon for convergence, should be computed from tol
-
+
double px = b;
if(fabs(b) < err)
return xk;
// return xk;
double G = d / px;
double H = G*G - f / px;
-
+
//std::cout << "G = " << G << "H = " << H;
double radicand = (n - 1)*(n*H-G*G);
//assert(radicand.real() > 0);
#endif
void subdiv_sbasis(SBasis const & s,
- std::vector<double> & roots,
+ std::vector<double> & roots,
double left, double right) {
Interval bs = bounds_fast(s);
if(bs.min() > 0 || bs.max() < 0)
std::vector<double> roots(SBasis const & s) {
switch(s.size()) {
- case 0:
- return std::vector<double>();
- case 1:
- return roots1(s);
- default:
- return sbasis_to_bezier(s).roots();
+ case 0:
+ return std::vector<double>();
+ case 1:
+ return roots1(s);
+ default:
+ {
+ Bezier bz;
+ sbasis_to_bezier(bz, s);
+ return bz.roots();
+ }
}
}
index cbddccda86cabbd8771797aa1fa282d12cf24afa..27e3047fd1beae5cc226713fb8f60be36da6a40d 100644 (file)
+/*
+ * Symmetric Power Basis - Bernstein Basis conversion routines
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ * Nathan Hurst <njh@mail.csse.monash.edu.au>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/choose.h>
+#include <2geom/svg-path.h>
+#include <2geom/exception.h>
+
+#include <iostream>
+
+
+
+
+namespace Geom
+{
+
+/*
+ * Symmetric Power Basis - Bernstein Basis conversion routines
+ *
+ * some remark about precision:
+ * interval [0,1], subdivisions: 10^3
+ * - bezier_to_sbasis : up to degree ~72 precision is at least 10^-5
+ * up to degree ~87 precision is at least 10^-3
+ * - sbasis_to_bezier : up to order ~63 precision is at least 10^-15
+ * precision is at least 10^-14 even beyond order 200
+ *
+ * interval [-1,1], subdivisions: 10^3
+ * - bezier_to_sbasis : up to degree ~21 precision is at least 10^-5
+ * up to degree ~24 precision is at least 10^-3
+ * - sbasis_to_bezier : up to order ~11 precision is at least 10^-5
+ * up to order ~13 precision is at least 10^-3
+ *
+ * interval [-10,10], subdivisions: 10^3
+ * - bezier_to_sbasis : up to degree ~7 precision is at least 10^-5
+ * up to degree ~8 precision is at least 10^-3
+ * - sbasis_to_bezier : up to order ~3 precision is at least 10^-5
+ * up to order ~4 precision is at least 10^-3
+ *
+ * references:
+ * this implementation is based on the following article:
+ * J.Sanchez-Reyes - The Symmetric Analogue of the Polynomial Power Basis
+ */
+
+inline
+double binomial(unsigned int n, unsigned int k)
+{
+ return choose<double>(n, k);
+}
+
+inline
+int sgn(unsigned int j, unsigned int k)
+{
+ assert (j >= k);
+ // we are sure that j >= k
+ return ((j-k) & 1u) ? -1 : 1;
+}
+
+
+void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz)
+{
+ // if the degree is even q is the order in the symmetrical power basis,
+ // if the degree is odd q is the order + 1
+ // n is always the polynomial degree, i. e. the Bezier order
+ size_t q, n;
+ bool even;
+ if (sz == 0)
+ {
+ q = sb.size();
+ if (sb[q-1][0] == sb[q-1][1])
+ {
+ even = true;
+ --q;
+ n = 2*q;
+ }
+ else
+ {
+ even = false;
+ n = 2*q-1;
+ }
+ }
+ else
+ {
+ q = (sz > sb.size()) ? sb.size() : sz;
+ n = 2*sz-1;
+ even = false;
+ }
+ bz.clear();
+ bz.resize(n+1);
+ double Tjk;
+ for (size_t k = 0; k < q; ++k)
+ {
+ for (size_t j = k; j < n-k; ++j) // j <= n-k-1
+ {
+ Tjk = binomial(n-2*k-1, j-k);
+ bz[j] += (Tjk * sb[k][0]);
+ bz[n-j] += (Tjk * sb[k][1]); // n-k <-> [k][1]
+ }
+ }
+ if (even)
+ {
+ bz[q] += sb[q][0];
+ }
+ // the resulting coefficients are with respect to the scaled Bernstein
+ // basis so we need to divide them by (n, j) binomial coefficient
+ for (size_t j = 1; j < n; ++j)
+ {
+ bz[j] /= binomial(n, j);
+ }
+ bz[0] = sb[0][0];
+ bz[n] = sb[0][1];
+}
+
+void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz)
+{
+ Bezier bzx, bzy;
+ sbasis_to_bezier(bzx, sb[X], sz);
+ sbasis_to_bezier(bzy, sb[Y], sz);
+ size_t n = (bzx.size() >= bzy.size()) ? bzx.size() : bzy.size();
+
+ bz.resize(n, Point(0,0));
+ for (size_t i = 0; i < bzx.size(); ++i)
+ {
+ bz[i][X] = bzx[i];
+ }
+ for (size_t i = 0; i < bzy.size(); ++i)
+ {
+ bz[i][Y] = bzy[i];
+ }
+}
+
+
+void bezier_to_sbasis (SBasis & sb, Bezier const& bz)
+{
+ // if the degree is even q is the order in the symmetrical power basis,
+ // if the degree is odd q is the order + 1
+ // n is always the polynomial degree, i. e. the Bezier order
+ size_t n = bz.order();
+ size_t q = (n+1) / 2;
+ size_t even = (n & 1u) ? 0 : 1;
+ sb.clear();
+ sb.resize(q + even, Linear(0, 0));
+ double Tjk;
+ for (size_t k = 0; k < q; ++k)
+ {
+ for (size_t j = k; j < q; ++j)
+ {
+ Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
+ sb[j][0] += (Tjk * bz[k]);
+ sb[j][1] += (Tjk * bz[n-k]); // n-j <-> [j][1]
+ }
+ for (size_t j = k+1; j < q; ++j)
+ {
+ Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
+ sb[j][0] += (Tjk * bz[n-k]);
+ sb[j][1] += (Tjk * bz[k]); // n-j <-> [j][1]
+ }
+ }
+ if (even)
+ {
+ for (size_t k = 0; k < q; ++k)
+ {
+ Tjk = sgn(q,k) * binomial(n, k);
+ sb[q][0] += (Tjk * (bz[k] + bz[n-k]));
+ }
+ sb[q][0] += (binomial(n, q) * bz[q]);
+ sb[q][1] = sb[q][0];
+ }
+ sb[0][0] = bz[0];
+ sb[0][1] = bz[n];
+}
+
+
+void bezier_to_sbasis (D2<SBasis> & sb, std::vector<Point> const& bz)
+{
+ size_t n = bz.size() - 1;
+ size_t q = (n+1) / 2;
+ size_t even = (n & 1u) ? 0 : 1;
+ sb[X].clear();
+ sb[Y].clear();
+ sb[X].resize(q + even, Linear(0, 0));
+ sb[Y].resize(q + even, Linear(0, 0));
+ double Tjk;
+ for (size_t k = 0; k < q; ++k)
+ {
+ for (size_t j = k; j < q; ++j)
+ {
+ Tjk = sgn(j, k) * binomial(n-j-k, j-k) * binomial(n, k);
+ sb[X][j][0] += (Tjk * bz[k][X]);
+ sb[X][j][1] += (Tjk * bz[n-k][X]);
+ sb[Y][j][0] += (Tjk * bz[k][Y]);
+ sb[Y][j][1] += (Tjk * bz[n-k][Y]);
+ }
+ for (size_t j = k+1; j < q; ++j)
+ {
+ Tjk = sgn(j, k) * binomial(n-j-k-1, j-k-1) * binomial(n, k);
+ sb[X][j][0] += (Tjk * bz[n-k][X]);
+ sb[X][j][1] += (Tjk * bz[k][X]);
+ sb[Y][j][0] += (Tjk * bz[n-k][Y]);
+ sb[Y][j][1] += (Tjk * bz[k][Y]);
+ }
+ }
+ if (even)
+ {
+ for (size_t k = 0; k < q; ++k)
+ {
+ Tjk = sgn(q,k) * binomial(n, k);
+ sb[X][q][0] += (Tjk * (bz[k][X] + bz[n-k][X]));
+ sb[Y][q][0] += (Tjk * (bz[k][Y] + bz[n-k][Y]));
+ }
+ sb[X][q][0] += (binomial(n, q) * bz[q][X]);
+ sb[X][q][1] = sb[X][q][0];
+ sb[Y][q][0] += (binomial(n, q) * bz[q][Y]);
+ sb[Y][q][1] = sb[Y][q][0];
+ }
+ sb[X][0][0] = bz[0][X];
+ sb[X][0][1] = bz[n][X];
+ sb[Y][0][0] = bz[0][Y];
+ sb[Y][0][1] = bz[n][Y];
+}
+
+
+} // end namespace Geom
+
+namespace Geom{
+#if 0
+
/* From Sanchez-Reyes 1997
W_{j,k} = W_{n0j, n-k} = choose(n-2k-1, j-k)choose(2k+1,k)/choose(n,j)
for k=0,...,q-1; j = k, ...,n-k-1
W_{q,q} = 1 (n even)
*/
-#include <2geom/sbasis-to-bezier.h>
-#include <2geom/choose.h>
-#include <2geom/svg-path.h>
-#include <iostream>
-#include <2geom/exception.h>
-
-namespace Geom{
-
double W(unsigned n, unsigned j, unsigned k) {
unsigned q = (n+1)/2;
if((n & 1) == 0 && j == q && k == q)
choose<double>(n,j);
}
+
// this produces a degree 2q bezier from a degree k sbasis
Bezier
sbasis_to_bezier(SBasis const &B, unsigned q) {
return (i&1)?-1:1;
}
+// WARNING: this is wrong!
// this produces a degree k sbasis from a degree 2q bezier
SBasis
bezier_to_sbasis(Bezier const &B) {
return D2<Bezier<order> >(sbasis_to_bezier<order>(B[0]), sbasis_to_bezier<order>(B[1]));
}
*/
+#endif
#if 0 // using old path
//std::vector<Geom::Point>
@@ -135,7 +387,7 @@ subpath_from_sbasis(Geom::OldPathSetBuilder &pb, D2<SBasis> const &B, double tol
/*
* This version works by inverting a reasonable upper bound on the error term after subdividing the
* curve at $a$. We keep biting off pieces until there is no more curve left.
-*
+*
* Derivation: The tail of the power series is $a_ks^k + a_{k+1}s^{k+1} + \ldots = e$. A
* subdivision at $a$ results in a tail error of $e*A^k, A = (1-a)a$. Let this be the desired
* tolerance tol $= e*A^k$ and invert getting $A = e^{1/k}$ and $a = 1/2 - \sqrt{1/4 - A}$
@@ -146,7 +398,7 @@ subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, doubl
double te = B.tail_error(k);
assert(B[0].IS_FINITE());
assert(B[1].IS_FINITE());
-
+
//std::cout << "tol = " << tol << std::endl;
while(1) {
double A = std::sqrt(tol/te); // pow(te, 1./k)
@@ -169,10 +421,10 @@ subpath_from_sbasis_incremental(Geom::OldPathSetBuilder &pb, D2<SBasis> B, doubl
initial = false;
}
pb.push_cubic(bez[1], bez[2], bez[3]);
-
+
// move to next piece of curve
if(a >= 1) break;
- B = compose(B, Linear(a, 1));
+ B = compose(B, Linear(a, 1));
te = B.tail_error(k);
}
}
@@ -190,7 +442,8 @@ void build_from_sbasis(Geom::PathBuilder &pb, D2<SBasis> const &B, double tol, b
if( !only_cubicbeziers && (sbasis_size(B) <= 1) ) {
pb.lineTo(B.at1());
} else {
- std::vector<Geom::Point> bez = sbasis_to_bezier(B, 2);
+ std::vector<Geom::Point> bez;
+ sbasis_to_bezier(bez, B, 2);
pb.curveTo(bez[1], bez[2], bez[3]);
}
} else {
@@ -246,7 +499,7 @@ path_from_piecewise(Geom::Piecewise<Geom::D2<Geom::SBasis> > const &B, double to
return pb.peek();
}
-};
+}
/*
Local Variables:
index 2fedb83c74ccd3c418ae570beb62399b08525ed1..c9d5cbbbca1973b8288a88dbfb059ccfb1b0263f 100644 (file)
#include <2geom/d2.h>
#include <2geom/path.h>
-namespace Geom{
+#include <vector>
+
+namespace Geom {
+
+void sbasis_to_bezier (Bezier & bz, SBasis const& sb, size_t sz = 0);
+void sbasis_to_bezier (std::vector<Point> & bz, D2<SBasis> const& sb, size_t sz = 0);
+void bezier_to_sbasis (SBasis & sb, Bezier const& bz);
+void bezier_to_sbasis (D2<SBasis> & sb, std::vector<Point> const& bz);
+
+
+#if 0
// this produces a degree k bezier from a degree k sbasis
Bezier
sbasis_to_bezier(SBasis const &B, unsigned q = 0);
std::vector<Geom::Point>
sbasis_to_bezier(D2<SBasis> const &B, unsigned q = 0);
+#endif
+
std::vector<Path> path_from_piecewise(Piecewise<D2<SBasis> > const &B, double tol, bool only_cubicbeziers = false);
Path path_from_sbasis(D2<SBasis> const &B, double tol, bool only_cubicbeziers = false);
inline Path cubicbezierpath_from_sbasis(D2<SBasis> const &B, double tol)
- { return path_from_sbasis(B, tol, true); } ;
+ { return path_from_sbasis(B, tol, true); }
+
+} // end namespace Geom
+
+
-};
#endif
/*
index b922e970effd478165526de22624591774d68dff..9a7b5633fff5d37fba9ecc8e9115ec28e8b03e5a 100644 (file)
const unsigned MAXDEPTH = 23; // Maximum depth for recursion. Using floats means 23 bits precision max
-const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */
-
+const double BEPSILON = ldexp(1.0,(-MAXDEPTH-1)); /*Flatness control value */
+const double SECANT_EPSILON = 1e-13; // secant method converges much faster, get a bit more precision
/**
* This function is called _a lot_. We have included various manual memory management stuff to reduce the amount of mallocing that goes on. In the future it is possible that this will hurt performance.
**/
double *Vtemp;
unsigned N,degree;
std::vector<double> &solutions;
- Bernsteins(int degr, std::vector<double> &so) : N(degr+1), degree(degr),solutions(so) {
+ bool use_secant;
+ Bernsteins(int degr, std::vector<double> &so) : N(degr+1), degree(degr),solutions(so), use_secant(false) {
Vtemp = new double[N*2];
}
~Bernsteins() {
double *Left,
double *Right);
+ double horner(const double *b, double t);
+
unsigned
control_poly_flat_enough(double const *V);
unsigned degree, /* The degree of the polynomial */
std::vector<double> &solutions, /* RETURN candidate t-values */
unsigned depth, /* The depth of the recursion */
- double left_t, double right_t)
+ double left_t, double right_t, bool use_secant)
{
Bernsteins B(degree, solutions);
+ B.use_secant = use_secant;
B.find_bernstein_roots(w, depth, left_t, right_t);
}
}
// I thought secant method would be faster here, but it'aint. -- njh
-
- // The original code went to some effort to try and terminate early when the curve is sufficiently flat. However, it seems that in practice it almost always bottoms out, so leaving this code out actually speeds things up
- if(0) if (control_poly_flat_enough(w)) {
- //printf("flatten out %d\n", depth);
- const double Ax = right_t - left_t;
- const double Ay = w[degree] - w[0];
-
- solutions.push_back(left_t - Ax*w[0] / Ay);
- return;
+ // Actually, it was, I just was using the wrong method for bezier evaluation. Horner's rule results in a very efficient algorithm - 10* faster (20080816)
+ // Future work: try using brent's method
+ if(use_secant) { // false position
+ double s = 0;double t = 1;
+ double e = 1e-10;
+ int n,side=0;
+ double r,fr,fs = w[0],ft = w[degree];
+
+ for (n = 1; n <= 100; n++)
+ {
+ r = (fs*t - ft*s) / (fs - ft);
+ if (fabs(t-s) < e*fabs(t+s)) break;
+ fr = horner(w, r);
+
+ if (fr * ft > 0)
+ {
+ t = r; ft = fr;
+ if (side==-1) fs /= 2;
+ side = -1;
+ }
+ else if (fs * fr > 0)
+ {
+ s = r; fs = fr;
+ if (side==+1) ft /= 2;
+ side = +1;
+ }
+ else break;
}
-
+ solutions.push_back(r*right_t + (1-r)*left_t);
+ return;
+ }
}
/* Otherwise, solve recursively after subdividing control polygon */
return error < BEPSILON * a;
}
+// suggested by Sederberg.
+double Bernsteins::horner(const double *b, double t) {
+ int n = degree;
+ double u, bc, tn, tmp;
+ int i;
+ u = 1.0 - t;
+ bc = 1;
+ tn = 1;
+ tmp = b[0]*u;
+ for(i=1; i<n; i++){
+ tn = tn*t;
+ bc = bc*(n-i+1)/i;
+ tmp = (tmp + tn*bc*b[i])*u;
+ }
+ return (tmp + tn*t*b[n]);
+}
};
diff --git a/src/2geom/solver.h b/src/2geom/solver.h
index 87365ab336553926d8a51e016e4ed37e0c917980..13153c22cbae43568823ecfd037f8169f788fc97 100644 (file)
--- a/src/2geom/solver.h
+++ b/src/2geom/solver.h
unsigned degree, /* The degree of the polynomial */
std::vector<double> & solutions, /* RETURN candidate t-values */
unsigned depth, /* The depth of the recursion */
- double left_t=0, double right_t=1);
+ double left_t=0, double right_t=1, bool use_secant=true);
};
#endif
index 5ae92097e9c2ac22ab8449e3f98415dc249dec51..23645fe11b0d2f82dafca2947c9abd39a401b9cd 100644 (file)
}
-// D(E(t,C),t) = E(t+PI/2,O)
+// D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center
+// the derivative doesn't rotate the ellipse but there is a translation
+// of the parameter t by an angle of PI/2 so the ellipse points are shifted
+// of such an angle in the cw direction
Curve* SVGEllipticalArc::derivative() const
{
if (isDegenerate() && is_svg_compliant())
Curve* SVGEllipticalArc::portion(double f, double t) const
{
+ // fix input arguments
if (f < 0) f = 0;
if (f > 1) f = 1;
if (t < 0) t = 0;
if (t > 1) t = 1;
+
if ( are_near(f, t) )
{
SVGEllipticalArc* arc = new SVGEllipticalArc();
arc->m_sweep = m_sweep;
arc->m_large_arc = m_large_arc;
}
+
SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
arc->m_initial_point = pointAt(f);
arc->m_final_point = pointAt(t);
Coord cos_rot_angle = std::cos(rotation_angle());
Coord sin_rot_angle = std::sin(rotation_angle());
// order = 4 seems to be enough to get a perfect looking elliptical arc
- // should it be choosen in function of the arc length anyway ?
- // or maybe a user settable parameter: toSBasis(unsigned int order) ?
SBasis arc_x = ray(X) * cos(param,4);
SBasis arc_y = ray(Y) * sin(param,4);
arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
Curve* SVGEllipticalArc::transformed(Matrix const& m) const
{
- // return SBasisCurve(toSBasis()).transformed(m);
-
Ellipse e(center(X), center(Y), ray(X), ray(Y), rotation_angle());
Ellipse et = e.transformed(m);
Point inner_point = pointAt(0.5);
return ea.duplicate();
}
-
+/*
+ * helper routine to convert the parameter t value
+ * btw [0,1] and [0,2PI] domain and back
+ *
+ */
Coord SVGEllipticalArc::map_to_02PI(Coord t) const
{
Coord angle = start_angle();
namespace detail
{
-
+/*
+ * ellipse_equation
+ *
+ * this is an helper struct, it provides two routines:
+ * the first one evaluates the implicit form of an ellipse on a given point
+ * the second one computes the normal versor at a given point of an ellipse
+ * in implicit form
+ */
struct ellipse_equation
{
ellipse_equation(double a, double b, double c, double d, double e, double f)
{
}
+/*
+ * check that the coefficients computed by the fit method satisfy
+ * the tolerance parameters at the k-th sample point
+ */
bool
make_elliptical_arc::
bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
{
dist_err = std::fabs( ee(p[k]) );
dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );
+ // check that the angle btw the tangent versor to the input curve
+ // and the normal versor of the elliptical arc, both evaluate
+ // at the k-th sample point, are really othogonal
angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );
//angle_err *= angle_err;
return ( dist_err > dist_bound || angle_err > angle_tol );
}
+/*
+ * check that the coefficients computed by the fit method satisfy
+ * the tolerance parameters at each sample point
+ */
bool
make_elliptical_arc::
check_bound(double A, double B, double C, double D, double E, double F)
{
- // check error magnitude
detail::ellipse_equation ee(A, B, C, D, E, F);
+ // check error magnitude at the end-points
double e1x = (2*A + B) * tol_at_extr;
double e1y = (B + 2*C) * tol_at_extr;
double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr;
- if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+ if (bound_exceeded(0, ee, e1x, e1y, e2))
{
print_bound_error(0);
return false;
}
- if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+ if (bound_exceeded(0, ee, e1x, e1y, e2))
{
print_bound_error(last);
return false;
}
+ // e1x = derivative((ee(x,y), x) | x->tolerance, y->tolerance
e1x = (2*A + B) * tolerance;
+ // e1y = derivative((ee(x,y), y) | x->tolerance, y->tolerance
e1y = (B + 2*C) * tolerance;
+ // e2 = ee(tolerance, tolerance) - F;
e2 = ((D + E) + (A + B + C) * tolerance) * tolerance;
// std::cerr << "e1x = " << e1x << std::endl;
// std::cerr << "e1y = " << e1y << std::endl;
// std::cerr << "e2 = " << e2 << std::endl;
+ // check error magnitude at sample points
for ( unsigned int k = 1; k < last; ++k )
{
if ( bound_exceeded(k, ee, e1x, e1y, e2) )
return true;
}
+/*
+ * fit
+ *
+ * supply the samples to the fitter and compute
+ * the ellipse implicit equation coefficients
+ */
void make_elliptical_arc::fit()
{
for (unsigned int k = 0; k < N; ++k)
index 7d5087c48961b2f8485e84ba7b3f595c27b0da9f..5c42e6e1d1e57fecb3e88590b77a7850e20378b6 100644 (file)
: m_initial_point(Point(0,0)), m_final_point(Point(0,0)),
m_rx(0), m_ry(0), m_rot_angle(0),
m_large_arc(true), m_sweep(true),
- m_svg_compliant(_svg_compliant)
- {
- m_start_angle = m_end_angle = 0;
- m_center = Point(0,0);
- }
-
+ m_svg_compliant(_svg_compliant),
+ m_start_angle(0), m_end_angle(0),
+ m_center(Point(0,0))
+ {
+ }
+
+ /*
+ * constructor
+ *
+ * input parameters:
+ * _initial_point: initial arc end point;
+ * _rx: ellipse x-axis ray length
+ * _ry: ellipse y-axis ray length
+ * _rot_angle: ellipse x-axis rotation angle;
+ * _large_arc: if true the largest arc is chosen,
+ * if false the smallest arc is chosen;
+ * _sweep : if true the clockwise arc is chosen,
+ * if false the counter-clockwise arc is chosen;
+ * _final_point: final arc end point;
+ * _svg_compliant: if true the class behaviour follows the Standard
+ * SVG 1.1 implementation guidelines (see Appendix F.6)
+ * if false the class behavoiur is more strict
+ * on input parameter
+ *
+ * in case the initial and the final arc end-points overlaps
+ * a degenerate arc of zero length is generated
+ *
+ */
SVGEllipticalArc( Point _initial_point, double _rx, double _ry,
double _rot_angle, bool _large_arc, bool _sweep,
Point _final_point,
std::vector<double> roots(double v, Dim2 d) const;
+ /*
+ * find all the points on the curve portion between "from" and "to"
+ * at the same smallest distance from the point "p" the points are returned
+ * as their parameter t value;
+ */
std::vector<double>
allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;
+ /*
+ * find a point on the curve portion between "from" and "to"
+ * at the same smallest distance from the point "p";
+ * the point is returned as its parameter t value;
+ */
double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
{
if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )
D2<SBasis> toSBasis() const;
+ /*
+ * return true if the angle argument (in radiants) is contained
+ * in the range [start_angle(), end_angle() ]
+ */
bool containsAngle(Coord angle) const;
+ /*
+ * return the value of the d-dimensional coordinate related to "t"
+ * here t belongs to the [0,2PI] domain
+ */
double valueAtAngle(Coord t, Dim2 d) const;
+ /*
+ * return the point related to the parameter value "t"
+ * here t belongs to the [0,2PI] domain
+ */
Point pointAtAngle(Coord t) const
{
double sin_rot_angle = std::sin(rotation_angle());
return p * m;
}
+ /*
+ * return the value of the d-dimensional coordinate related to "t"
+ * here t belongs to the [0,1] domain
+ */
double valueAt(Coord t, Dim2 d) const
{
if (isDegenerate() && is_svg_compliant())
return valueAtAngle(tt, d);
}
+ /*
+ * return the point related to the parameter value "t"
+ * here t belongs to the [0,1] domain
+ */
Point pointAt(Coord t) const
{
if (isDegenerate() && is_svg_compliant())
return rarc;
}
+
double sweep_angle() const
{
Coord d = end_angle() - start_angle();
Point m_initial_point, m_final_point;
double m_rx, m_ry, m_rot_angle;
bool m_large_arc, m_sweep;
+ bool m_svg_compliant;
double m_start_angle, m_end_angle;
Point m_center;
- bool m_svg_compliant;
}; // end class SVGEllipticalArc
+
+/*
+ * useful for testing and debugging
+ */
template< class charT >
inline
std::basic_ostream<charT> &
+// forward declation
namespace detail
{
struct ellipse_equation;
}
-
+/*
+ * make_elliptical_arc
+ *
+ * convert a parametric polynomial curve given in symmetric power basis form
+ * into an SVGEllipticalArc type; in order to be successfull the input curve
+ * has to look like an actual elliptical arc even if a certain tolerance
+ * is allowed through an ad-hoc parameter.
+ * The conversion is performed through an interpolation on a certain amount of
+ * sample points computed on the input curve;
+ * the interpolation computes the coefficients of the general implicit equation
+ * of an ellipse (A*X^2 + B*XY + C*Y^2 + D*X + E*Y + F = 0), then from the
+ * implicit equation we compute the parametric form.
+ *
+ */
class make_elliptical_arc
{
public:
typedef D2<SBasis> curve_type;
+ /*
+ * constructor
+ *
+ * it doesn't execute the conversion but set the input and output parameters
+ *
+ * _ea: the output SVGEllipticalArc that will be generated;
+ * _curve: the input curve to be converted;
+ * _total_samples: the amount of sample points to be taken
+ * on the input curve for performing the conversion
+ * _tolerance: how much likelihood is required between the input curve
+ * and the generated elliptical arc; the smaller it is the
+ * the tolerance the higher it is the likelihood.
+ */
make_elliptical_arc( SVGEllipticalArc& _ea,
curve_type const& _curve,
unsigned int _total_samples,
}
public:
+ /*
+ * perform the actual conversion
+ * return true if the conversion is successfull, false on the contrary
+ */
bool operator()()
{
+ // initialize the reference
const NL::Vector & coeff = fitter.result();
fit();
if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )
return true;
}
+ /*
+ * you can set a boolean parameter to tell the conversion routine
+ * if the output elliptical arc has to be svg compliant or not;
+ * the default value is true
+ */
bool svg_compliant_flag() const
{
return svg_compliant;
}
private:
- SVGEllipticalArc& ea;
- const curve_type & curve;
- Piecewise<D2<SBasis> > dcurve;
- NL::LFMEllipse model;
+ SVGEllipticalArc& ea; // output elliptical arc
+ const curve_type & curve; // input curve
+ Piecewise<D2<SBasis> > dcurve; // derivative of the input curve
+ NL::LFMEllipse model; // model used for fitting
+ // perform the actual fitting task
NL::least_squeares_fitter<NL::LFMEllipse> fitter;
+ // tolerance: the user-defined tolerance parameter;
+ // tol_at_extr: the tolerance at end-points automatically computed
+ // on the value of "tolerance", and usually more strict;
+ // tol_at_center: tolerance at the center of the ellipse
+ // angle_tol: tolerance for the angle btw the input curve tangent
+ // versor and the ellipse normal versor at the sample points
double tolerance, tol_at_extr, tol_at_center, angle_tol;
- Point initial_point, final_point;
- unsigned int N;
+ Point initial_point, final_point; // initial and final end-points
+ unsigned int N; // total samples
unsigned int last; // N-1
double partitions; // N-1
- std::vector<Point> p; // sample points
+ std::vector<Point> p; // sample points
double dist_err, dist_bound, angle_err;
bool svg_compliant;
};
diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp
index bed04d3ffcfd1144dfcb6ac5a17290e6ad957844..898c72bf5c85b74625b72089b93801218c987668 100644 (file)
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
namespace Geom {
void output(Curve const &curve, SVGPathSink &sink) {
- std::vector<Point> pts = sbasis_to_bezier(curve.toSBasis(), 2); //TODO: use something better!
+ std::vector<Point> pts;
+ sbasis_to_bezier(pts, curve.toSBasis(), 2); //TODO: use something better!
sink.curveTo(pts[0], pts[1], pts[2]);
}
diff --git a/src/2geom/sweep.cpp b/src/2geom/sweep.cpp
index 1b881399036103801c60a7cfb02e44d757db34f4..227c822bdd1b417b08df92b8a1a5de92f5f10032 100644 (file)
--- a/src/2geom/sweep.cpp
+++ b/src/2geom/sweep.cpp
std::vector<std::vector<unsigned> > pairs(rs.size());
for(unsigned i = 0; i < rs.size(); i++) {
- events.push_back(Event(rs[i].left(), i, false));
- events.push_back(Event(rs[i].right(), i, true));
+ if(!rs[i].isEmpty()) {
+ events.push_back(Event(rs[i].left(), i, false));
+ events.push_back(Event(rs[i].right(), i, true));
+ }
}
std::sort(events.begin(), events.end());
unsigned sz = n ? b.size() : a.size();
events[n].reserve(sz*2);
for(unsigned i = 0; i < sz; i++) {
- events[n].push_back(Event(n ? b[i].left() : a[i].left(), i, false));
- events[n].push_back(Event(n ? b[i].right() : a[i].right(), i, true));
+ Rect r = n ? b[i] : a[i];
+ if(!r.isEmpty()) {
+ events[n].push_back(Event(r.left(), i, false));
+ events[n].push_back(Event(r.right(), i, true));
+ }
}
std::sort(events[n].begin(), events[n].end());
}