recursive-bezier-intersection.cpp revision 07717cb695cd79993e74908804d04e658a95482f
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
unsigned intersect_steps = 0;
namespace Geom {
class OldBezier {
public:
OldBezier() {
}
Point operator()(double t) const;
~OldBezier() {}
// Compute bounding box for a
minax = p[0][X]; // These are the most likely to be extremal
if( p[i][X] < minax )
minax = p[i][X];
else if( p[i][X] > maxax )
maxax = p[i][X];
}
minay = p[0][Y]; // These are the most likely to be extremal
if( p[i][Y] < minay )
minay = p[i][Y];
else if( p[i][Y] > maxay )
maxay = p[i][Y];
}
}
};
static void
OldBezier a,
OldBezier b);
void
double precision) {
OldBezier a, b;
a.p = A;
b.p = B;
return find_intersections_bezier_recursive(xs, a,b);
}
/* The value of 1.0 / (1L<<14) is enough for most applications */
/*
* 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
* to hold an intermediate value until it is no longer needed.
*/
/* Copy control points */
/* Triangle computation */
for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
}
}
for (unsigned j = 0; j < sz; j++)
for (unsigned j = 0; j < sz; 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
* to hold an intermediate value until it is no longer needed.
*/
/* Copy control points */
/* Triangle computation */
for (unsigned i = 1; i < sz; i++) {
for (unsigned j = 0; j < sz - i; j++) {
}
}
}
#endif
// suggested by Sederberg.
int n = p.size()-1;
int i;
Point r;
u = 1.0 - t;
bc = 1;
tn = 1;
for(i=1; i<n; i++){
}
}
return r;
}
/*
* Test the bounding boxes of two OldBezier curves for interference.
* 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
* 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
* 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.
* Third, in the final comparison (the interference test) the comparisons
* should both permit equality. We want to find intersections even if they
* occur at the ends of segments.
* Finally, there are tighter bounding boxes that can be derived. It isn't
* clear whether the higher probability of rejection (and hence fewer
* subdivisions and tests) is worth the extra work.
*/
// Test bounding box of b against bounding box of a
// Not >= : need boundary case
}
/*
* 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
* 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
* 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
* 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
* 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
* 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
* is robust: a near-tangential intersection will yield zero or two
* intersections.
*/
{
intersect_steps ++;
//std::cout << deptha << std::endl;
if( deptha > 0 )
{
OldBezier A[2];
deptha--;
if( depthb > 0 )
{
OldBezier B[2];
depthb--;
if( intersect_BB( A[0], B[0] ) )
parameters );
if( intersect_BB( A[1], B[0] ) )
parameters );
if( intersect_BB( A[0], B[1] ) )
parameters );
parameters );
}
else
{
if( intersect_BB( A[0], b ) )
parameters );
if( intersect_BB( A[1], b ) )
parameters );
}
}
else
if( depthb > 0 )
{
OldBezier B[2];
depthb--;
if( intersect_BB( a, B[0] ) )
parameters );
if( intersect_BB( a, B[1] ) )
parameters );
}
else // Both segments are fully subdivided; now do line segments
{
double xmk = b.p[0][X] - a.p[0][X];
double ymk = b.p[0][Y] - a.p[0][Y];
return;
else
{
if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
return;
}
}
}
/*
* Wang's theorem is used to estimate the level of subdivision required,
* but only if the bounding boxes interfere at the top level.
* Assuming there is a possible intersection, recursively_intersect is
* used to find all the parameters corresponding to intersection points.
* these are then sorted and returned in an array.
*/
}
unsigned wangs_theorem(OldBezier a) {
return 6; // seems a good approximation!
unsigned ra;
ra = 0;
else
//std::cout << ra << std::endl;
return ra;
}
struct rparams
{
OldBezier &A;
OldBezier &B;
};
static int
gsl_vector * f)
{
const double x0 = gsl_vector_get (x, 0);
gsl_vector_set (f, 0, dx[0]);
return GSL_SUCCESS;
}
union dbl_64{
long long i64;
double d64;
};
{
dbl_64 s;
return s.d64;
}
static void intersect_polish_root (OldBezier &A, double &s,
OldBezier &B, double &t) {
const gsl_multiroot_fsolver_type *T;
int status;
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_multiroot_fsolver_set (sol, &f, x);
do
{
iter++;
if (status) /* check if solver is stuck */
break;
status =
}
s = gsl_vector_get (sol->x, 0);
gsl_vector_free (x);
// This code does a neighbourhood search for minor improvements.
//std::cout << "------\n" << best_v << std::endl;
while (true) {
double c = L1(A(n[0]) - B(n[1]));
//std::cout << c << "; ";
if (c < trial_v) {
trial = n;
trial_v = c;
}
}
}
//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 {
}
}
}
{
if( intersect_BB( a, b ) )
{
b, 0., 1., wangs_theorem(b),
xs);
}
/*for(unsigned i = 0; i < xs.size(); i++)
intersect_polish_root(a, xs[i].first,
b, xs[i].second);*/
}
};
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :