| /************************************************************** |
| * |
| * Licensed to the Apache Software Foundation (ASF) under one |
| * or more contributor license agreements. See the NOTICE file |
| * distributed with this work for additional information |
| * regarding copyright ownership. The ASF licenses this file |
| * to you under the Apache License, Version 2.0 (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.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, |
| * software distributed under the License is distributed on an |
| * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
| * KIND, either express or implied. See the License for the |
| * specific language governing permissions and limitations |
| * under the License. |
| * |
| *************************************************************/ |
| |
| |
| |
| // MARKER(update_precomp.py): autogen include statement, do not remove |
| #include "precompiled_basegfx.hxx" |
| |
| #include <algorithm> |
| #include <iterator> |
| #include <vector> |
| #include <utility> |
| |
| #include <math.h> |
| |
| #include "bezierclip.hxx" |
| #include "gauss.hxx" |
| |
| |
| |
| // what to test |
| #define WITH_ASSERTIONS |
| //#define WITH_CONVEXHULL_TEST |
| //#define WITH_MULTISUBDIVIDE_TEST |
| //#define WITH_FATLINE_TEST |
| //#define WITH_CALCFOCUS_TEST |
| //#define WITH_SAFEPARAMBASE_TEST |
| //#define WITH_SAFEPARAMS_TEST |
| //#define WITH_SAFEPARAM_DETAILED_TEST |
| //#define WITH_SAFEFOCUSPARAM_CALCFOCUS |
| //#define WITH_SAFEFOCUSPARAM_TEST |
| //#define WITH_SAFEFOCUSPARAM_DETAILED_TEST |
| #define WITH_BEZIERCLIP_TEST |
| |
| |
| |
| // ----------------------------------------------------------------------------- |
| |
| /* Implementation of the so-called 'Fat-Line Bezier Clipping Algorithm' by Sederberg et al. |
| * |
| * Actual reference is: T. W. Sederberg and T Nishita: Curve |
| * intersection using Bezier clipping. In Computer Aided Design, 22 |
| * (9), 1990, pp. 538--549 |
| */ |
| |
| // ----------------------------------------------------------------------------- |
| |
| /* Misc helper |
| * =========== |
| */ |
| int fallFac( int n, int k ) |
| { |
| #ifdef WITH_ASSERTIONS |
| assert(n>=k); // "For factorials, n must be greater or equal k" |
| assert(n>=0); // "For factorials, n must be positive" |
| assert(k>=0); // "For factorials, k must be positive" |
| #endif |
| |
| int res( 1 ); |
| |
| while( k-- && n ) res *= n--; |
| |
| return res; |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| int fac( int n ) |
| { |
| return fallFac(n, n); |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| /* Bezier fat line clipping part |
| * ============================= |
| */ |
| |
| // ----------------------------------------------------------------------------- |
| |
| void Impl_calcFatLine( FatLine& line, const Bezier& c ) |
| { |
| // Prepare normalized implicit line |
| // ================================ |
| |
| // calculate vector orthogonal to p1-p4: |
| line.a = -(c.p0.y - c.p3.y); |
| line.b = (c.p0.x - c.p3.x); |
| |
| // normalize |
| const double len( sqrt( line.a*line.a + line.b*line.b ) ); |
| if( !tolZero(len) ) |
| { |
| line.a /= len; |
| line.b /= len; |
| } |
| |
| line.c = -(line.a*c.p0.x + line.b*c.p0.y); |
| |
| |
| // Determine bounding fat line from it |
| // =================================== |
| |
| // calc control point distances |
| const double dP2( calcLineDistance(line.a, line.b, line.c, c.p1.x, c.p1.y ) ); |
| const double dP3( calcLineDistance(line.a, line.b, line.c, c.p2.x, c.p2.y ) ); |
| |
| // calc approximate bounding lines to curve (tight bounds are |
| // possible here, but more expensive to calculate and thus not |
| // worth the overhead) |
| if( dP2 * dP3 > 0.0 ) |
| { |
| line.dMin = 3.0/4.0 * ::std::min(0.0, ::std::min(dP2, dP3)); |
| line.dMax = 3.0/4.0 * ::std::max(0.0, ::std::max(dP2, dP3)); |
| } |
| else |
| { |
| line.dMin = 4.0/9.0 * ::std::min(0.0, ::std::min(dP2, dP3)); |
| line.dMax = 4.0/9.0 * ::std::max(0.0, ::std::max(dP2, dP3)); |
| } |
| } |
| |
| void Impl_calcBounds( Point2D& leftTop, |
| Point2D& rightBottom, |
| const Bezier& c1 ) |
| { |
| leftTop.x = ::std::min( c1.p0.x, ::std::min( c1.p1.x, ::std::min( c1.p2.x, c1.p3.x ) ) ); |
| leftTop.y = ::std::min( c1.p0.y, ::std::min( c1.p1.y, ::std::min( c1.p2.y, c1.p3.y ) ) ); |
| rightBottom.x = ::std::max( c1.p0.x, ::std::max( c1.p1.x, ::std::max( c1.p2.x, c1.p3.x ) ) ); |
| rightBottom.y = ::std::max( c1.p0.y, ::std::max( c1.p1.y, ::std::max( c1.p2.y, c1.p3.y ) ) ); |
| } |
| |
| bool Impl_doBBoxIntersect( const Bezier& c1, |
| const Bezier& c2 ) |
| { |
| // calc rectangular boxes from c1 and c2 |
| Point2D lt1; |
| Point2D rb1; |
| Point2D lt2; |
| Point2D rb2; |
| |
| Impl_calcBounds( lt1, rb1, c1 ); |
| Impl_calcBounds( lt2, rb2, c2 ); |
| |
| if( ::std::min(rb1.x, rb2.x) < ::std::max(lt1.x, lt2.x) || |
| ::std::min(rb1.y, rb2.y) < ::std::max(lt1.y, lt2.y) ) |
| { |
| return false; |
| } |
| else |
| { |
| return true; |
| } |
| } |
| |
| /* calculates two t's for the given bernstein control polygon: the first is |
| * the intersection of the min value line with the convex hull from |
| * the left, the second is the intersection of the max value line with |
| * the convex hull from the right. |
| */ |
| bool Impl_calcSafeParams( double& t1, |
| double& t2, |
| const Polygon2D& rPoly, |
| double lowerYBound, |
| double upperYBound ) |
| { |
| // need the convex hull of the control polygon, as this is |
| // guaranteed to completely bound the curve |
| Polygon2D convHull( convexHull(rPoly) ); |
| |
| // init min and max buffers |
| t1 = 0.0 ; |
| double currLowerT( 1.0 ); |
| |
| t2 = 1.0; |
| double currHigherT( 0.0 ); |
| |
| if( convHull.size() <= 1 ) |
| return false; // only one point? Then we're done with clipping |
| |
| /* now, clip against lower and higher bounds */ |
| Point2D p0; |
| Point2D p1; |
| |
| bool bIntersection( false ); |
| |
| for( Polygon2D::size_type i=0; i<convHull.size(); ++i ) |
| { |
| // have to check against convHull.size() segments, as the |
| // convex hull is, by definition, closed. Thus, for the |
| // last point, we take the first point as partner. |
| if( i+1 == convHull.size() ) |
| { |
| // close the polygon |
| p0 = convHull[i]; |
| p1 = convHull[0]; |
| } |
| else |
| { |
| p0 = convHull[i]; |
| p1 = convHull[i+1]; |
| } |
| |
| // is the segment in question within or crossing the |
| // horizontal band spanned by lowerYBound and upperYBound? If |
| // not, we've got no intersection. If yes, we maybe don't have |
| // an intersection, but we've got to update the permissible |
| // range, nevertheless. This is because inside lying segments |
| // leads to full range forbidden. |
| if( (tolLessEqual(p0.y, upperYBound) || tolLessEqual(p1.y, upperYBound)) && |
| (tolGreaterEqual(p0.y, lowerYBound) || tolGreaterEqual(p1.y, lowerYBound)) ) |
| { |
| // calc intersection of convex hull segment with |
| // one of the horizontal bounds lines |
| const double r_x( p1.x - p0.x ); |
| const double r_y( p1.y - p0.y ); |
| |
| if( tolZero(r_y) ) |
| { |
| // r_y is virtually zero, thus we've got a horizontal |
| // line. Now check whether we maybe coincide with lower or |
| // upper horizonal bound line. |
| if( tolEqual(p0.y, lowerYBound) || |
| tolEqual(p0.y, upperYBound) ) |
| { |
| // yes, simulate intersection then |
| currLowerT = ::std::min(currLowerT, ::std::min(p0.x, p1.x)); |
| currHigherT = ::std::max(currHigherT, ::std::max(p0.x, p1.x)); |
| } |
| } |
| else |
| { |
| // check against lower and higher bounds |
| // ===================================== |
| |
| // calc intersection with horizontal dMin line |
| const double currTLow( (lowerYBound - p0.y) * r_x / r_y + p0.x ); |
| |
| // calc intersection with horizontal dMax line |
| const double currTHigh( (upperYBound - p0.y) * r_x / r_y + p0.x ); |
| |
| currLowerT = ::std::min(currLowerT, ::std::min(currTLow, currTHigh)); |
| currHigherT = ::std::max(currHigherT, ::std::max(currTLow, currTHigh)); |
| } |
| |
| // set flag that at least one segment is contained or |
| // intersects given horizontal band. |
| bIntersection = true; |
| } |
| } |
| |
| #ifndef WITH_SAFEPARAMBASE_TEST |
| // limit intersections found to permissible t parameter range |
| t1 = ::std::max(0.0, currLowerT); |
| t2 = ::std::min(1.0, currHigherT); |
| #endif |
| |
| return bIntersection; |
| } |
| |
| |
| /* calculates two t's for the given bernstein polynomial: the first is |
| * the intersection of the min value line with the convex hull from |
| * the left, the second is the intersection of the max value line with |
| * the convex hull from the right. |
| * |
| * The polynomial coefficients c0 to c3 given to this method |
| * must correspond to t values of 0, 1/3, 2/3 and 1, respectively. |
| */ |
| bool Impl_calcSafeParams_clip( double& t1, |
| double& t2, |
| const FatLine& bounds, |
| double c0, |
| double c1, |
| double c2, |
| double c3 ) |
| { |
| /* first of all, determine convex hull of c0-c3 */ |
| Polygon2D poly(4); |
| poly[0] = Point2D(0, c0); |
| poly[1] = Point2D(1.0/3.0, c1); |
| poly[2] = Point2D(2.0/3.0, c2); |
| poly[3] = Point2D(1, c3); |
| |
| #ifndef WITH_SAFEPARAM_DETAILED_TEST |
| |
| return Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax ); |
| |
| #else |
| bool bRet( Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax ) ); |
| |
| Polygon2D convHull( convexHull( poly ) ); |
| |
| cout << "# convex hull testing" << endl |
| << "plot [t=0:1] "; |
| cout << " bez(" |
| << poly[0].x << "," |
| << poly[1].x << "," |
| << poly[2].x << "," |
| << poly[3].x << ",t),bez(" |
| << poly[0].y << "," |
| << poly[1].y << "," |
| << poly[2].y << "," |
| << poly[3].y << ",t), " |
| << "t, " << bounds.dMin << ", " |
| << "t, " << bounds.dMax << ", " |
| << t1 << ", t, " |
| << t2 << ", t, " |
| << "'-' using ($1):($2) title \"control polygon\" with lp, " |
| << "'-' using ($1):($2) title \"convex hull\" with lp" << endl; |
| |
| unsigned int k; |
| for( k=0; k<poly.size(); ++k ) |
| { |
| cout << poly[k].x << " " << poly[k].y << endl; |
| } |
| cout << poly[0].x << " " << poly[0].y << endl; |
| cout << "e" << endl; |
| |
| for( k=0; k<convHull.size(); ++k ) |
| { |
| cout << convHull[k].x << " " << convHull[k].y << endl; |
| } |
| cout << convHull[0].x << " " << convHull[0].y << endl; |
| cout << "e" << endl; |
| |
| return bRet; |
| #endif |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| void Impl_deCasteljauAt( Bezier& part1, |
| Bezier& part2, |
| const Bezier& input, |
| double t ) |
| { |
| // deCasteljau bezier arc, scheme is: |
| // |
| // First row is C_0^n,C_1^n,...,C_n^n |
| // Second row is P_1^n,...,P_n^n |
| // etc. |
| // with P_k^r = (1 - x_s)P_{k-1}^{r-1} + x_s P_k{r-1} |
| // |
| // this results in: |
| // |
| // P1 P2 P3 P4 |
| // L1 P2 P3 R4 |
| // L2 H R3 |
| // L3 R2 |
| // L4/R1 |
| if( tolZero(t) ) |
| { |
| // t is zero -> part2 is input curve, part1 is empty (input.p0, that is) |
| part1.p0.x = part1.p1.x = part1.p2.x = part1.p3.x = input.p0.x; |
| part1.p0.y = part1.p1.y = part1.p2.y = part1.p3.y = input.p0.y; |
| part2 = input; |
| } |
| else if( tolEqual(t, 1.0) ) |
| { |
| // t is one -> part1 is input curve, part2 is empty (input.p3, that is) |
| part1 = input; |
| part2.p0.x = part2.p1.x = part2.p2.x = part2.p3.x = input.p3.x; |
| part2.p0.y = part2.p1.y = part2.p2.y = part2.p3.y = input.p3.y; |
| } |
| else |
| { |
| part1.p0.x = input.p0.x; part1.p0.y = input.p0.y; |
| part1.p1.x = (1.0 - t)*part1.p0.x + t*input.p1.x; part1.p1.y = (1.0 - t)*part1.p0.y + t*input.p1.y; |
| const double Hx ( (1.0 - t)*input.p1.x + t*input.p2.x ), Hy ( (1.0 - t)*input.p1.y + t*input.p2.y ); |
| part1.p2.x = (1.0 - t)*part1.p1.x + t*Hx; part1.p2.y = (1.0 - t)*part1.p1.y + t*Hy; |
| part2.p3.x = input.p3.x; part2.p3.y = input.p3.y; |
| part2.p2.x = (1.0 - t)*input.p2.x + t*input.p3.x; part2.p2.y = (1.0 - t)*input.p2.y + t*input.p3.y; |
| part2.p1.x = (1.0 - t)*Hx + t*part2.p2.x; part2.p1.y = (1.0 - t)*Hy + t*part2.p2.y; |
| part2.p0.x = (1.0 - t)*part1.p2.x + t*part2.p1.x; part2.p0.y = (1.0 - t)*part1.p2.y + t*part2.p1.y; |
| part1.p3.x = part2.p0.x; part1.p3.y = part2.p0.y; |
| } |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| void printCurvesWithSafeRange( const Bezier& c1, const Bezier& c2, double t1_c1, double t2_c1, |
| const Bezier& c2_part, const FatLine& bounds_c2 ) |
| { |
| static int offset = 0; |
| |
| cout << "# safe param range testing" << endl |
| << "plot [t=0.0:1.0] "; |
| |
| // clip safe ranges off c1 |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| // subdivide at t1_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1_c1 ); |
| // subdivide at t2_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, t2_c1 ); |
| |
| // output remaining segment (c1_part1) |
| |
| cout << "bez(" |
| << c1.p0.x+offset << "," |
| << c1.p1.x+offset << "," |
| << c1.p2.x+offset << "," |
| << c1.p3.x+offset << ",t),bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",t), bez(" |
| << c2.p0.x+offset << "," |
| << c2.p1.x+offset << "," |
| << c2.p2.x+offset << "," |
| << c2.p3.x+offset << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t), " |
| #if 1 |
| << "bez(" |
| << c1_part1.p0.x+offset << "," |
| << c1_part1.p1.x+offset << "," |
| << c1_part1.p2.x+offset << "," |
| << c1_part1.p3.x+offset << ",t),bez(" |
| << c1_part1.p0.y << "," |
| << c1_part1.p1.y << "," |
| << c1_part1.p2.y << "," |
| << c1_part1.p3.y << ",t), " |
| #endif |
| #if 1 |
| << "bez(" |
| << c2_part.p0.x+offset << "," |
| << c2_part.p1.x+offset << "," |
| << c2_part.p2.x+offset << "," |
| << c2_part.p3.x+offset << ",t),bez(" |
| << c2_part.p0.y << "," |
| << c2_part.p1.y << "," |
| << c2_part.p2.y << "," |
| << c2_part.p3.y << ",t), " |
| #endif |
| << "linex(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c << ",t)+" << offset << ", liney(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c << ",t) title \"fat line (center)\", linex(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c-bounds_c2.dMin << ",t)+" << offset << ", liney(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c-bounds_c2.dMin << ",t) title \"fat line (min) \", linex(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c-bounds_c2.dMax << ",t)+" << offset << ", liney(" |
| << bounds_c2.a << "," |
| << bounds_c2.b << "," |
| << bounds_c2.c-bounds_c2.dMax << ",t) title \"fat line (max) \"" << endl; |
| |
| offset += 1; |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| void printResultWithFinalCurves( const Bezier& c1, const Bezier& c1_part, |
| const Bezier& c2, const Bezier& c2_part, |
| double t1_c1, double t2_c1 ) |
| { |
| static int offset = 0; |
| |
| cout << "# final result" << endl |
| << "plot [t=0.0:1.0] "; |
| |
| cout << "bez(" |
| << c1.p0.x+offset << "," |
| << c1.p1.x+offset << "," |
| << c1.p2.x+offset << "," |
| << c1.p3.x+offset << ",t),bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",t), bez(" |
| << c1_part.p0.x+offset << "," |
| << c1_part.p1.x+offset << "," |
| << c1_part.p2.x+offset << "," |
| << c1_part.p3.x+offset << ",t),bez(" |
| << c1_part.p0.y << "," |
| << c1_part.p1.y << "," |
| << c1_part.p2.y << "," |
| << c1_part.p3.y << ",t), " |
| << " pointmarkx(bez(" |
| << c1.p0.x+offset << "," |
| << c1.p1.x+offset << "," |
| << c1.p2.x+offset << "," |
| << c1.p3.x+offset << "," |
| << t1_c1 << "),t), " |
| << " pointmarky(bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << "," |
| << t1_c1 << "),t), " |
| << " pointmarkx(bez(" |
| << c1.p0.x+offset << "," |
| << c1.p1.x+offset << "," |
| << c1.p2.x+offset << "," |
| << c1.p3.x+offset << "," |
| << t2_c1 << "),t), " |
| << " pointmarky(bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << "," |
| << t2_c1 << "),t), " |
| |
| << "bez(" |
| << c2.p0.x+offset << "," |
| << c2.p1.x+offset << "," |
| << c2.p2.x+offset << "," |
| << c2.p3.x+offset << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t), " |
| << "bez(" |
| << c2_part.p0.x+offset << "," |
| << c2_part.p1.x+offset << "," |
| << c2_part.p2.x+offset << "," |
| << c2_part.p3.x+offset << ",t),bez(" |
| << c2_part.p0.y << "," |
| << c2_part.p1.y << "," |
| << c2_part.p2.y << "," |
| << c2_part.p3.y << ",t)" << endl; |
| |
| offset += 1; |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| /** determine parameter ranges [0,t1) and (t2,1] on c1, where c1 is guaranteed to lie outside c2. |
| Returns false, if the two curves don't even intersect. |
| |
| @param t1 |
| Range [0,t1) on c1 is guaranteed to lie outside c2 |
| |
| @param t2 |
| Range (t2,1] on c1 is guaranteed to lie outside c2 |
| |
| @param c1_orig |
| Original curve c1 |
| |
| @param c1_part |
| Subdivided current part of c1 |
| |
| @param c2_orig |
| Original curve c2 |
| |
| @param c2_part |
| Subdivided current part of c2 |
| */ |
| bool Impl_calcClipRange( double& t1, |
| double& t2, |
| const Bezier& c1_orig, |
| const Bezier& c1_part, |
| const Bezier& c2_orig, |
| const Bezier& c2_part ) |
| { |
| // TODO: Maybe also check fat line orthogonal to P0P3, having P0 |
| // and P3 as the extremal points |
| |
| if( Impl_doBBoxIntersect(c1_part, c2_part) ) |
| { |
| // Calculate fat lines around c1 |
| FatLine bounds_c2; |
| |
| // must use the subdivided version of c2, since the fat line |
| // algorithm works implicitely with the convex hull bounding |
| // box. |
| Impl_calcFatLine(bounds_c2, c2_part); |
| |
| // determine clip positions on c2. Can use original c1 (which |
| // is necessary anyway, to get the t's on the original curve), |
| // since the distance calculations work directly in the |
| // Bernstein polynom parameter domain. |
| if( Impl_calcSafeParams_clip( t1, t2, bounds_c2, |
| calcLineDistance( bounds_c2.a, |
| bounds_c2.b, |
| bounds_c2.c, |
| c1_orig.p0.x, |
| c1_orig.p0.y ), |
| calcLineDistance( bounds_c2.a, |
| bounds_c2.b, |
| bounds_c2.c, |
| c1_orig.p1.x, |
| c1_orig.p1.y ), |
| calcLineDistance( bounds_c2.a, |
| bounds_c2.b, |
| bounds_c2.c, |
| c1_orig.p2.x, |
| c1_orig.p2.y ), |
| calcLineDistance( bounds_c2.a, |
| bounds_c2.b, |
| bounds_c2.c, |
| c1_orig.p3.x, |
| c1_orig.p3.y ) ) ) |
| { |
| //printCurvesWithSafeRange(c1_orig, c2_orig, t1, t2, c2_part, bounds_c2); |
| |
| // they do intersect |
| return true; |
| } |
| } |
| |
| // they don't intersect: nothing to do |
| return false; |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| /* Tangent intersection part |
| * ========================= |
| */ |
| |
| // ----------------------------------------------------------------------------- |
| |
| void Impl_calcFocus( Bezier& res, const Bezier& c ) |
| { |
| // arbitrary small value, for now |
| // TODO: find meaningful value |
| const double minPivotValue( 1.0e-20 ); |
| |
| Point2D::value_type fMatrix[6]; |
| Point2D::value_type fRes[2]; |
| |
| // calc new curve from hodograph, c and linear blend |
| |
| // Coefficients for derivative of c are (C_i=n(C_{i+1} - C_i)): |
| // |
| // 3(P1 - P0), 3(P2 - P1), 3(P3 - P2) (bezier curve of degree 2) |
| // |
| // The hodograph is then (bezier curve of 2nd degree is P0(1-t)^2 + 2P1(1-t)t + P2t^2): |
| // |
| // 3(P1 - P0)(1-t)^2 + 6(P2 - P1)(1-t)t + 3(P3 - P2)t^2 |
| // |
| // rotate by 90 degrees: x=-y, y=x and you get the normal vector function N(t): |
| // |
| // x(t) = -(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) |
| // y(t) = 3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2 |
| // |
| // Now, the focus curve is defined to be F(t)=P(t) + c(t)N(t), |
| // where P(t) is the original curve, and c(t)=c0(1-t) + c1 t |
| // |
| // This results in the following expression for F(t): |
| // |
| // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - |
| // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) |
| // |
| // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1.t)t^2 + P3.y t^3 + |
| // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2) |
| // |
| // As a heuristic, we set F(0)=F(1) (thus, the curve is closed and _tends_ to be small): |
| // |
| // For F(0), the following results: |
| // |
| // x(0) = P0.x - c0 3(P1.y - P0.y) |
| // y(0) = P0.y + c0 3(P1.x - P0.x) |
| // |
| // For F(1), the following results: |
| // |
| // x(1) = P3.x - c1 3(P3.y - P2.y) |
| // y(1) = P3.y + c1 3(P3.x - P2.x) |
| // |
| // Reorder, collect and substitute into F(0)=F(1): |
| // |
| // P0.x - c0 3(P1.y - P0.y) = P3.x - c1 3(P3.y - P2.y) |
| // P0.y + c0 3(P1.x - P0.x) = P3.y + c1 3(P3.x - P2.x) |
| // |
| // which yields |
| // |
| // (P0.y - P1.y)c0 + (P3.y - P2.y)c1 = (P3.x - P0.x)/3 |
| // (P1.x - P0.x)c0 + (P2.x - P3.x)c1 = (P3.y - P0.y)/3 |
| // |
| |
| // so, this is what we calculate here (determine c0 and c1): |
| fMatrix[0] = c.p1.x - c.p0.x; |
| fMatrix[1] = c.p2.x - c.p3.x; |
| fMatrix[2] = (c.p3.y - c.p0.y)/3.0; |
| fMatrix[3] = c.p0.y - c.p1.y; |
| fMatrix[4] = c.p3.y - c.p2.y; |
| fMatrix[5] = (c.p3.x - c.p0.x)/3.0; |
| |
| // TODO: determine meaningful value for |
| if( !solve(fMatrix, 2, 3, fRes, minPivotValue) ) |
| { |
| // TODO: generate meaningful values here |
| // singular or nearly singular system -- use arbitrary |
| // values for res |
| fRes[0] = 0.0; |
| fRes[1] = 1.0; |
| |
| cerr << "Matrix singular!" << endl; |
| } |
| |
| // now, the reordered and per-coefficient collected focus curve is |
| // the following third degree bezier curve F(t): |
| // |
| // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - |
| // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) |
| // = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - |
| // (3c0P1.y(1-t)^3 - 3c0P0.y(1-t)^3 + 6c0P2.y(1-t)^2t - 6c0P1.y(1-t)^2t + |
| // 3c0P3.y(1-t)t^2 - 3c0P2.y(1-t)t^2 + |
| // 3c1P1.y(1-t)^2t - 3c1P0.y(1-t)^2t + 6c1P2.y(1-t)t^2 - 6c1P1.y(1-t)t^2 + |
| // 3c1P3.yt^3 - 3c1P2.yt^3) |
| // = (P0.x - 3 c0 P1.y + 3 c0 P0.y)(1-t)^3 + |
| // 3(P1.x - c1 P1.y + c1 P0.y - 2 c0 P2.y + 2 c0 P1.y)(1-t)^2t + |
| // 3(P2.x - 2 c1 P2.y + 2 c1 P1.y - c0 P3.y + c0 P2.y)(1-t)t^2 + |
| // (P3.x - 3 c1 P3.y + 3 c1 P2.y)t^3 |
| // = (P0.x - 3 c0(P1.y - P0.y))(1-t)^3 + |
| // 3(P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))(1-t)^2t + |
| // 3(P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))(1-t)t^2 + |
| // (P3.x - 3 c1(P3.y - P2.y))t^3 |
| // |
| // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 + |
| // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2) |
| // = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 + |
| // 3c0(P1.x - P0.x)(1-t)^3 + 6c0(P2.x - P1.x)(1-t)^2t + 3c0(P3.x - P2.x)(1-t)t^2 + |
| // 3c1(P1.x - P0.x)(1-t)^2t + 6c1(P2.x - P1.x)(1-t)t^2 + 3c1(P3.x - P2.x)t^3 |
| // = (P0.y + 3 c0 (P1.x - P0.x))(1-t)^3 + |
| // 3(P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))(1-t)^2t + |
| // 3(P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))(1-t)t^2 + |
| // (P3.y + 3 c1 (P3.x - P2.x))t^3 |
| // |
| // Therefore, the coefficients F0 to F3 of the focus curve are: |
| // |
| // F0.x = (P0.x - 3 c0(P1.y - P0.y)) F0.y = (P0.y + 3 c0 (P1.x - P0.x)) |
| // F1.x = (P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y)) F1.y = (P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x)) |
| // F2.x = (P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y)) F2.y = (P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x)) |
| // F3.x = (P3.x - 3 c1(P3.y - P2.y)) F3.y = (P3.y + 3 c1 (P3.x - P2.x)) |
| // |
| res.p0.x = c.p0.x - 3*fRes[0]*(c.p1.y - c.p0.y); |
| res.p1.x = c.p1.x - fRes[1]*(c.p1.y - c.p0.y) - 2*fRes[0]*(c.p2.y - c.p1.y); |
| res.p2.x = c.p2.x - 2*fRes[1]*(c.p2.y - c.p1.y) - fRes[0]*(c.p3.y - c.p2.y); |
| res.p3.x = c.p3.x - 3*fRes[1]*(c.p3.y - c.p2.y); |
| |
| res.p0.y = c.p0.y + 3*fRes[0]*(c.p1.x - c.p0.x); |
| res.p1.y = c.p1.y + 2*fRes[0]*(c.p2.x - c.p1.x) + fRes[1]*(c.p1.x - c.p0.x); |
| res.p2.y = c.p2.y + fRes[0]*(c.p3.x - c.p2.x) + 2*fRes[1]*(c.p2.x - c.p1.x); |
| res.p3.y = c.p3.y + 3*fRes[1]*(c.p3.x - c.p2.x); |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| bool Impl_calcSafeParams_focus( double& t1, |
| double& t2, |
| const Bezier& curve, |
| const Bezier& focus ) |
| { |
| // now, we want to determine which normals of the original curve |
| // P(t) intersect with the focus curve F(t). The condition for |
| // this statement is P'(t)(P(t) - F) = 0, i.e. hodograph P'(t) and |
| // line through P(t) and F are perpendicular. |
| // If you expand this equation, you end up with something like |
| // |
| // (\sum_{i=0}^n (P_i - F)B_i^n(t))^T (\sum_{j=0}^{n-1} n(P_{j+1} - P_j)B_j^{n-1}(t)) |
| // |
| // Multiplying that out (as the scalar product is linear, we can |
| // extract some terms) yields: |
| // |
| // (P_i - F)^T n(P_{j+1} - P_j) B_i^n(t)B_j^{n-1}(t) + ... |
| // |
| // If we combine the B_i^n(t)B_j^{n-1}(t) product, we arrive at a |
| // Bernstein polynomial of degree 2n-1, as |
| // |
| // \binom{n}{i}(1-t)^{n-i}t^i) \binom{n-1}{j}(1-t)^{n-1-j}t^j) = |
| // \binom{n}{i}\binom{n-1}{j}(1-t)^{2n-1-i-j}t^{i+j} |
| // |
| // Thus, with the defining equation for a 2n-1 degree Bernstein |
| // polynomial |
| // |
| // \sum_{i=0}^{2n-1} d_i B_i^{2n-1}(t) |
| // |
| // the d_i are calculated as follows: |
| // |
| // d_i = \sum_{j+k=i, j\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{j}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_j - F) |
| // |
| // |
| // Okay, but F is now not a single point, but itself a curve |
| // F(u). Thus, for every value of u, we get a different 2n-1 |
| // bezier curve from the above equation. Therefore, we have a |
| // tensor product bezier patch, with the following defining |
| // equation: |
| // |
| // d(t,u) = \sum_{i=0}^{2n-1} \sum_{j=0}^m B_i^{2n-1}(t) B_j^{m}(u) d_{ij}, where |
| // d_{ij} = \sum_{k+l=i, l\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{l}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_l - F_j) |
| // |
| // as above, only that now F is one of the focus' control points. |
| // |
| // Note the difference in the binomial coefficients to the |
| // reference paper, these formulas most probably contained a typo. |
| // |
| // To determine, where D(t,u) is _not_ zero (these are the parts |
| // of the curve that don't share normals with the focus and can |
| // thus be safely clipped away), we project D(u,t) onto the |
| // (d(t,u), t) plane, determine the convex hull there and proceed |
| // as for the curve intersection part (projection is orthogonal to |
| // u axis, thus simply throw away u coordinate). |
| // |
| // \fallfac are so-called falling factorials (see Concrete |
| // Mathematics, p. 47 for a definition). |
| // |
| |
| // now, for tensor product bezier curves, the convex hull property |
| // holds, too. Thus, we simply project the control points (t_{ij}, |
| // u_{ij}, d_{ij}) onto the (t,d) plane and calculate the |
| // intersections of the convex hull with the t axis, as for the |
| // bezier clipping case. |
| |
| // |
| // calc polygon of control points (t_{ij}, d_{ij}): |
| // |
| const int n( 3 ); // cubic bezier curves, as a matter of fact |
| const int i_card( 2*n ); |
| const int j_card( n + 1 ); |
| const int k_max( n-1 ); |
| Polygon2D controlPolygon( i_card*j_card ); // vector of (t_{ij}, d_{ij}) in row-major order |
| |
| int i, j, k, l; // variable notation from formulas above and Sederberg article |
| Point2D::value_type d; |
| for( i=0; i<i_card; ++i ) |
| { |
| for( j=0; j<j_card; ++j ) |
| { |
| // calc single d_{ij} sum: |
| for( d=0.0, k=::std::max(0,i-n); k<=k_max && k<=i; ++k ) |
| { |
| l = i - k; // invariant: k + l = i |
| assert(k>=0 && k<=n-1); // k \in {0,...,n-1} |
| assert(l>=0 && l<=n); // l \in {0,...,n} |
| |
| // TODO: find, document and assert proper limits for n and int's max_val. |
| // This becomes important should anybody wants to use |
| // this code for higher-than-cubic beziers |
| d += static_cast<double>(fallFac(n,l)*fallFac(n-1,k)*fac(i)) / |
| static_cast<double>(fac(l)*fac(k) * fallFac(2*n-1,i)) * n * |
| ( (curve[k+1].x - curve[k].x)*(curve[l].x - focus[j].x) + // dot product here |
| (curve[k+1].y - curve[k].y)*(curve[l].y - focus[j].y) ); |
| } |
| |
| // Note that the t_{ij} values are evenly spaced on the |
| // [0,1] interval, thus t_{ij}=i/(2n-1) |
| controlPolygon[ i*j_card + j ] = Point2D( i/(2.0*n-1.0), d ); |
| } |
| } |
| |
| #ifndef WITH_SAFEFOCUSPARAM_DETAILED_TEST |
| |
| // calc safe parameter range, to determine [0,t1] and [t2,1] where |
| // no zero crossing is guaranteed. |
| return Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ); |
| |
| #else |
| bool bRet( Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ) ); |
| |
| Polygon2D convHull( convexHull( controlPolygon ) ); |
| |
| cout << "# convex hull testing (focus)" << endl |
| << "plot [t=0:1] "; |
| cout << "'-' using ($1):($2) title \"control polygon\" with lp, " |
| << "'-' using ($1):($2) title \"convex hull\" with lp" << endl; |
| |
| unsigned int count; |
| for( count=0; count<controlPolygon.size(); ++count ) |
| { |
| cout << controlPolygon[count].x << " " << controlPolygon[count].y << endl; |
| } |
| cout << controlPolygon[0].x << " " << controlPolygon[0].y << endl; |
| cout << "e" << endl; |
| |
| for( count=0; count<convHull.size(); ++count ) |
| { |
| cout << convHull[count].x << " " << convHull[count].y << endl; |
| } |
| cout << convHull[0].x << " " << convHull[0].y << endl; |
| cout << "e" << endl; |
| |
| return bRet; |
| #endif |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| /** Calc all values t_i on c1, for which safeRanges functor does not |
| give a safe range on c1 and c2. |
| |
| This method is the workhorse of the bezier clipping. Because c1 |
| and c2 must be alternatingly tested against each other (first |
| determine safe parameter interval on c1 with regard to c2, then |
| the other way around), we call this method recursively with c1 and |
| c2 swapped. |
| |
| @param result |
| Output iterator where the final t values are added to. If curves |
| don't intersect, nothing is added. |
| |
| @param delta |
| Maximal allowed distance to true critical point (measured in the |
| original curve's coordinate system) |
| |
| @param safeRangeFunctor |
| Functor object, that must provide the following operator(): |
| bool safeRangeFunctor( double& t1, |
| double& t2, |
| const Bezier& c1_orig, |
| const Bezier& c1_part, |
| const Bezier& c2_orig, |
| const Bezier& c2_part ); |
| This functor must calculate the safe ranges [0,t1] and [t2,1] on |
| c1_orig, where c1_orig is 'safe' from c2_part. If the whole |
| c1_orig is safe, false must be returned, true otherwise. |
| */ |
| template <class Functor> void Impl_applySafeRanges_rec( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result, |
| double delta, |
| const Functor& safeRangeFunctor, |
| int recursionLevel, |
| const Bezier& c1_orig, |
| const Bezier& c1_part, |
| double last_t1_c1, |
| double last_t2_c1, |
| const Bezier& c2_orig, |
| const Bezier& c2_part, |
| double last_t1_c2, |
| double last_t2_c2 ) |
| { |
| // check end condition |
| // =================== |
| |
| // TODO: tidy up recursion handling. maybe put everything in a |
| // struct and swap that here at method entry |
| |
| // TODO: Implement limit on recursion depth. Should that limit be |
| // reached, chances are that we're on a higher-order tangency. For |
| // this case, AW proposed to take the middle of the current |
| // interval, and to correct both curve's tangents at that new |
| // endpoint to be equal. That virtually generates a first-order |
| // tangency, and justifies to return a single intersection |
| // point. Otherwise, inside/outside test might fail here. |
| |
| for( int i=0; i<recursionLevel; ++i ) cerr << " "; |
| if( recursionLevel % 2 ) |
| { |
| cerr << "level: " << recursionLevel |
| << " t: " |
| << last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0 |
| << ", c1: " << last_t1_c2 << " " << last_t2_c2 |
| << ", c2: " << last_t1_c1 << " " << last_t2_c1 |
| << endl; |
| } |
| else |
| { |
| cerr << "level: " << recursionLevel |
| << " t: " |
| << last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0 |
| << ", c1: " << last_t1_c1 << " " << last_t2_c1 |
| << ", c2: " << last_t1_c2 << " " << last_t2_c2 |
| << endl; |
| } |
| |
| // refine solution |
| // =============== |
| |
| double t1_c1, t2_c1; |
| |
| // Note: we first perform the clipping and only test for precision |
| // sufficiency afterwards, since we want to exploit the fact that |
| // Impl_calcClipRange returns false if the curves don't |
| // intersect. We would have to check that separately for the end |
| // condition, otherwise. |
| |
| // determine safe range on c1_orig |
| if( safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ) ) |
| { |
| // now, t1 and t2 are calculated on the original curve |
| // (but against a fat line calculated from the subdivided |
| // c2, namely c2_part). If the [t1,t2] range is outside |
| // our current [last_t1,last_t2] range, we're done in this |
| // branch - the curves no longer intersect. |
| if( tolLessEqual(t1_c1, last_t2_c1) && tolGreaterEqual(t2_c1, last_t1_c1) ) |
| { |
| // As noted above, t1 and t2 are calculated on the |
| // original curve, but against a fat line |
| // calculated from the subdivided c2, namely |
| // c2_part. Our domain to work on is |
| // [last_t1,last_t2], on the other hand, so values |
| // of [t1,t2] outside that range are irrelevant |
| // here. Clip range appropriately. |
| t1_c1 = ::std::max(t1_c1, last_t1_c1); |
| t2_c1 = ::std::min(t2_c1, last_t2_c1); |
| |
| // TODO: respect delta |
| // for now, end condition is just a fixed threshold on the t's |
| |
| // check end condition |
| // =================== |
| |
| #if 1 |
| if( fabs(last_t2_c1 - last_t1_c1) < 0.0001 && |
| fabs(last_t2_c2 - last_t1_c2) < 0.0001 ) |
| #else |
| if( fabs(last_t2_c1 - last_t1_c1) < 0.01 && |
| fabs(last_t2_c2 - last_t1_c2) < 0.01 ) |
| #endif |
| { |
| // done. Add to result |
| if( recursionLevel % 2 ) |
| { |
| // uneven level: have to swap the t's, since curves are swapped, too |
| *result++ = ::std::make_pair( last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0, |
| last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0 ); |
| } |
| else |
| { |
| *result++ = ::std::make_pair( last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0, |
| last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0 ); |
| } |
| |
| #if 0 |
| //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, last_t1_c1, last_t2_c1 ); |
| printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, t1_c1, t2_c1 ); |
| #else |
| // calc focus curve of c2 |
| Bezier focus; |
| Impl_calcFocus(focus, c2_part); // need to use subdivided c2 |
| |
| safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ); |
| |
| //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, t1_c1, t2_c1 ); |
| printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, last_t1_c1, last_t2_c1 ); |
| #endif |
| } |
| else |
| { |
| // heuristic: if parameter range is not reduced by at least |
| // 20%, subdivide longest curve, and clip shortest against |
| // both parts of longest |
| // if( (last_t2_c1 - last_t1_c1 - t2_c1 + t1_c1) / (last_t2_c1 - last_t1_c1) < 0.2 ) |
| if( false ) |
| { |
| // subdivide and descend |
| // ===================== |
| |
| Bezier part1; |
| Bezier part2; |
| |
| double intervalMiddle; |
| |
| if( last_t2_c1 - last_t1_c1 > last_t2_c2 - last_t1_c2 ) |
| { |
| // subdivide c1 |
| // ============ |
| |
| intervalMiddle = last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0; |
| |
| // subdivide at the middle of the interval (as |
| // we're not subdividing on the original |
| // curve, this simply amounts to subdivision |
| // at 0.5) |
| Impl_deCasteljauAt( part1, part2, c1_part, 0.5 ); |
| |
| // and descend recursively with swapped curves |
| Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, |
| c2_orig, c2_part, last_t1_c2, last_t2_c2, |
| c1_orig, part1, last_t1_c1, intervalMiddle ); |
| |
| Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, |
| c2_orig, c2_part, last_t1_c2, last_t2_c2, |
| c1_orig, part2, intervalMiddle, last_t2_c1 ); |
| } |
| else |
| { |
| // subdivide c2 |
| // ============ |
| |
| intervalMiddle = last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0; |
| |
| // subdivide at the middle of the interval (as |
| // we're not subdividing on the original |
| // curve, this simply amounts to subdivision |
| // at 0.5) |
| Impl_deCasteljauAt( part1, part2, c2_part, 0.5 ); |
| |
| // and descend recursively with swapped curves |
| Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, |
| c2_orig, part1, last_t1_c2, intervalMiddle, |
| c1_orig, c1_part, last_t1_c1, last_t2_c1 ); |
| |
| Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, |
| c2_orig, part2, intervalMiddle, last_t2_c2, |
| c1_orig, c1_part, last_t1_c1, last_t2_c1 ); |
| } |
| } |
| else |
| { |
| // apply calculated clip |
| // ===================== |
| |
| // clip safe ranges off c1_orig |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| // subdivide at t1_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c1_orig, t1_c1 ); |
| |
| // subdivide at t2_c1. As we're working on |
| // c1_part2 now, we have to adapt t2_c1 since |
| // we're no longer in the original parameter |
| // interval. This is based on the following |
| // assumption: t2_new = (t2-t1)/(1-t1), which |
| // relates the t2 value into the new parameter |
| // range [0,1] of c1_part2. |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2_c1-t1_c1)/(1.0-t1_c1) ); |
| |
| // descend with swapped curves and c1_part1 as the |
| // remaining (middle) segment |
| Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, |
| c2_orig, c2_part, last_t1_c2, last_t2_c2, |
| c1_orig, c1_part1, t1_c1, t2_c1 ); |
| } |
| } |
| } |
| } |
| } |
| |
| // ----------------------------------------------------------------------------- |
| |
| struct ClipBezierFunctor |
| { |
| bool operator()( double& t1_c1, |
| double& t2_c1, |
| const Bezier& c1_orig, |
| const Bezier& c1_part, |
| const Bezier& c2_orig, |
| const Bezier& c2_part ) const |
| { |
| return Impl_calcClipRange( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ); |
| } |
| }; |
| |
| // ----------------------------------------------------------------------------- |
| |
| struct BezierTangencyFunctor |
| { |
| bool operator()( double& t1_c1, |
| double& t2_c1, |
| const Bezier& c1_orig, |
| const Bezier& c1_part, |
| const Bezier& c2_orig, |
| const Bezier& c2_part ) const |
| { |
| // calc focus curve of c2 |
| Bezier focus; |
| Impl_calcFocus(focus, c2_part); // need to use subdivided c2 |
| // here, as the whole curve is |
| // used for focus calculation |
| |
| // determine safe range on c1_orig |
| bool bRet( Impl_calcSafeParams_focus( t1_c1, t2_c1, |
| c1_orig, // use orig curve here, need t's on original curve |
| focus ) ); |
| |
| cerr << "range: " << t2_c1 - t1_c1 << ", ret: " << bRet << endl; |
| |
| return bRet; |
| } |
| }; |
| |
| // ----------------------------------------------------------------------------- |
| |
| /** Perform a bezier clip (curve against curve) |
| |
| @param result |
| Output iterator where the final t values are added to. This |
| iterator will remain empty, if there are no intersections. |
| |
| @param delta |
| Maximal allowed distance to true intersection (measured in the |
| original curve's coordinate system) |
| */ |
| void clipBezier( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result, |
| double delta, |
| const Bezier& c1, |
| const Bezier& c2 ) |
| { |
| #if 0 |
| // first of all, determine list of collinear normals. Collinear |
| // normals typically separate two intersections, thus, subdivide |
| // at all collinear normal's t values beforehand. This will cater |
| // for tangent intersections, where two or more intersections are |
| // infinitesimally close together. |
| |
| // TODO: evaluate effects of higher-than-second-order |
| // tangencies. Sederberg et al. state that collinear normal |
| // algorithm then degrades quickly. |
| |
| ::std::vector< ::std::pair<double,double> > results; |
| ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(results); |
| |
| Impl_calcCollinearNormals( ii, delta, 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 ); |
| |
| // As Sederberg's collinear normal theorem is only sufficient, not |
| // necessary for two intersections left and right, we've to test |
| // all segments generated by the collinear normal algorithm |
| // against each other. In other words, if the two curves are both |
| // divided in a left and a right part, the collinear normal |
| // theorem does _not_ state that the left part of curve 1 does not |
| // e.g. intersect with the right part of curve 2. |
| |
| // divide c1 and c2 at collinear normal intersection points |
| ::std::vector< Bezier > c1_segments( results.size()+1 ); |
| ::std::vector< Bezier > c2_segments( results.size()+1 ); |
| Bezier c1_remainder( c1 ); |
| Bezier c2_remainder( c2 ); |
| unsigned int i; |
| for( i=0; i<results.size(); ++i ) |
| { |
| Bezier c1_part2; |
| Impl_deCasteljauAt( c1_segments[i], c1_part2, c1_remainder, results[i].first ); |
| c1_remainder = c1_part2; |
| |
| Bezier c2_part2; |
| Impl_deCasteljauAt( c2_segments[i], c2_part2, c2_remainder, results[i].second ); |
| c2_remainder = c2_part2; |
| } |
| c1_segments[i] = c1_remainder; |
| c2_segments[i] = c2_remainder; |
| |
| // now, c1/c2_segments contain all segments, then |
| // clip every resulting segment against every other |
| unsigned int c1_curr, c2_curr; |
| for( c1_curr=0; c1_curr<c1_segments.size(); ++c1_curr ) |
| { |
| for( c2_curr=0; c2_curr<c2_segments.size(); ++c2_curr ) |
| { |
| if( c1_curr != c2_curr ) |
| { |
| Impl_clipBezier_rec(result, delta, 0, |
| c1_segments[c1_curr], c1_segments[c1_curr], |
| 0.0, 1.0, |
| c2_segments[c2_curr], c2_segments[c2_curr], |
| 0.0, 1.0); |
| } |
| } |
| } |
| #else |
| Impl_applySafeRanges_rec( result, delta, BezierTangencyFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 ); |
| //Impl_applySafeRanges_rec( result, delta, ClipBezierFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 ); |
| #endif |
| // that's it, boys'n'girls! |
| } |
| |
| int main(int argc, const char *argv[]) |
| { |
| double curr_Offset( 0 ); |
| unsigned int i,j,k; |
| |
| Bezier someCurves[] = |
| { |
| // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.0)}, |
| // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)}, |
| // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)} |
| // {Point2D(0.25+1,0.5),Point2D(0.25+1,0.708333),Point2D(0.423611+1,0.916667),Point2D(0.770833+1,0.980324)}, |
| // {Point2D(0.0+1,0.0),Point2D(0.0+1,1.0),Point2D(1.0+1,1.0),Point2D(1.0+1,0.5)} |
| |
| // tangency1 |
| // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)}, |
| // {Point2D(0.0,1.0),Point2D(0.1,1.0),Point2D(0.4,1.0),Point2D(0.5,1.0)} |
| |
| // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)}, |
| // {Point2D(0.60114,0.933091),Point2D(0.69461,0.969419),Point2D(0.80676,0.992976),Point2D(0.93756,0.998663)} |
| // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)}, |
| // {Point2D(0.62712,0.828427),Point2D(0.76305,0.828507),Point2D(0.88555,0.77312),Point2D(0.95069,0.67325)} |
| |
| // clipping1 |
| // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)}, |
| // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)} |
| |
| // tangency2 |
| // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)}, |
| // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)} |
| |
| // tangency3 |
| // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)}, |
| // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)} |
| |
| // tangency4 |
| // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)}, |
| // {Point2D(0.26,0.4),Point2D(0.25,0.5),Point2D(0.25,0.5),Point2D(0.26,0.6)} |
| |
| // tangency5 |
| // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)}, |
| // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)} |
| |
| // tangency6 |
| // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)}, |
| // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)} |
| |
| // tangency7 |
| // {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)}, |
| // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)} |
| |
| // tangency Sederberg example |
| {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)}, |
| {Point2D(5.33+9.311,0.0),Point2D(5.33+9.311-13.279,4.205),Point2D(5.33+9.311-10.681,9.119),Point2D(5.33+9.311-2.603,10.254)} |
| |
| // clipping2 |
| // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)}, |
| // {Point2D(0.2575,0.4),Point2D(0.2475,0.5),Point2D(0.2475,0.5),Point2D(0.2575,0.6)} |
| |
| // {Point2D(0.0,0.1),Point2D(0.2,3.5),Point2D(1.0,-2.5),Point2D(1.1,1.2)}, |
| // {Point2D(0.0,1.0),Point2D(3.5,0.9),Point2D(-2.5,0.1),Point2D(1.1,0.2)} |
| // {Point2D(0.0,0.1),Point2D(0.2,3.0),Point2D(1.0,-2.0),Point2D(1.1,1.2)}, |
| // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)} |
| // {Point2D(0.0,1.0),Point2D(3.0,0.9),Point2D(-2.0,0.1),Point2D(1.1,0.2)} |
| // {Point2D(0.0,4.0),Point2D(0.1,5.0),Point2D(0.9,5.0),Point2D(1.0,4.0)}, |
| // {Point2D(0.0,0.0),Point2D(0.1,0.5),Point2D(0.9,0.5),Point2D(1.0,0.0)}, |
| // {Point2D(0.0,0.1),Point2D(0.1,1.5),Point2D(0.9,1.5),Point2D(1.0,0.1)}, |
| // {Point2D(0.0,-4.0),Point2D(0.1,-5.0),Point2D(0.9,-5.0),Point2D(1.0,-4.0)} |
| }; |
| |
| // output gnuplot setup |
| cout << "#!/usr/bin/gnuplot -persist" << endl |
| << "#" << endl |
| << "# automatically generated by bezierclip, don't change!" << endl |
| << "#" << endl |
| << "set parametric" << endl |
| << "bez(p,q,r,s,t) = p*(1-t)**3+q*3*(1-t)**2*t+r*3*(1-t)*t**2+s*t**3" << endl |
| << "bezd(p,q,r,s,t) = 3*(q-p)*(1-t)**2+6*(r-q)*(1-t)*t+3*(s-r)*t**2" << endl |
| << "pointmarkx(c,t) = c-0.03*t" << endl |
| << "pointmarky(c,t) = c+0.03*t" << endl |
| << "linex(a,b,c,t) = a*-c + t*-b" << endl |
| << "liney(a,b,c,t) = b*-c + t*a" << endl << endl |
| << "# end of setup" << endl << endl; |
| |
| #ifdef WITH_CONVEXHULL_TEST |
| // test convex hull algorithm |
| const double convHull_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# convex hull testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Polygon2D aTestPoly(4); |
| aTestPoly[0] = someCurves[i].p0; |
| aTestPoly[1] = someCurves[i].p1; |
| aTestPoly[2] = someCurves[i].p2; |
| aTestPoly[3] = someCurves[i].p3; |
| |
| aTestPoly[0].x += convHull_xOffset; |
| aTestPoly[1].x += convHull_xOffset; |
| aTestPoly[2].x += convHull_xOffset; |
| aTestPoly[3].x += convHull_xOffset; |
| |
| cout << " bez(" |
| << aTestPoly[0].x << "," |
| << aTestPoly[1].x << "," |
| << aTestPoly[2].x << "," |
| << aTestPoly[3].x << ",t),bez(" |
| << aTestPoly[0].y << "," |
| << aTestPoly[1].y << "," |
| << aTestPoly[2].y << "," |
| << aTestPoly[3].y << ",t), '-' using ($1):($2) title \"convex hull " << i << "\" with lp"; |
| |
| if( i+1<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Polygon2D aTestPoly(4); |
| aTestPoly[0] = someCurves[i].p0; |
| aTestPoly[1] = someCurves[i].p1; |
| aTestPoly[2] = someCurves[i].p2; |
| aTestPoly[3] = someCurves[i].p3; |
| |
| aTestPoly[0].x += convHull_xOffset; |
| aTestPoly[1].x += convHull_xOffset; |
| aTestPoly[2].x += convHull_xOffset; |
| aTestPoly[3].x += convHull_xOffset; |
| |
| Polygon2D convHull( convexHull(aTestPoly) ); |
| |
| for( k=0; k<convHull.size(); ++k ) |
| { |
| cout << convHull[k].x << " " << convHull[k].y << endl; |
| } |
| cout << convHull[0].x << " " << convHull[0].y << endl; |
| cout << "e" << endl; |
| } |
| #endif |
| |
| #ifdef WITH_MULTISUBDIVIDE_TEST |
| // test convex hull algorithm |
| const double multiSubdivide_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# multi subdivide testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Bezier c( someCurves[i] ); |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| c.p0.x += multiSubdivide_xOffset; |
| c.p1.x += multiSubdivide_xOffset; |
| c.p2.x += multiSubdivide_xOffset; |
| c.p3.x += multiSubdivide_xOffset; |
| |
| const double t1( 0.1+i/(3.0*sizeof(someCurves)/sizeof(Bezier)) ); |
| const double t2( 0.9-i/(3.0*sizeof(someCurves)/sizeof(Bezier)) ); |
| |
| // subdivide at t1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c, t1 ); |
| |
| // subdivide at t2_c1. As we're working on |
| // c1_part2 now, we have to adapt t2_c1 since |
| // we're no longer in the original parameter |
| // interval. This is based on the following |
| // assumption: t2_new = (t2-t1)/(1-t1), which |
| // relates the t2 value into the new parameter |
| // range [0,1] of c1_part2. |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) ); |
| |
| // subdivide at t2 |
| Impl_deCasteljauAt( c1_part3, c1_part2, c, t2 ); |
| |
| cout << " bez(" |
| << c1_part1.p0.x << "," |
| << c1_part1.p1.x << "," |
| << c1_part1.p2.x << "," |
| << c1_part1.p3.x << ",t), bez(" |
| << c1_part1.p0.y+0.01 << "," |
| << c1_part1.p1.y+0.01 << "," |
| << c1_part1.p2.y+0.01 << "," |
| << c1_part1.p3.y+0.01 << ",t) title \"middle " << i << "\", " |
| << " bez(" |
| << c1_part2.p0.x << "," |
| << c1_part2.p1.x << "," |
| << c1_part2.p2.x << "," |
| << c1_part2.p3.x << ",t), bez(" |
| << c1_part2.p0.y << "," |
| << c1_part2.p1.y << "," |
| << c1_part2.p2.y << "," |
| << c1_part2.p3.y << ",t) title \"right " << i << "\", " |
| << " bez(" |
| << c1_part3.p0.x << "," |
| << c1_part3.p1.x << "," |
| << c1_part3.p2.x << "," |
| << c1_part3.p3.x << ",t), bez(" |
| << c1_part3.p0.y << "," |
| << c1_part3.p1.y << "," |
| << c1_part3.p2.y << "," |
| << c1_part3.p3.y << ",t) title \"left " << i << "\""; |
| |
| |
| if( i+1<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| #endif |
| |
| #ifdef WITH_FATLINE_TEST |
| // test fatline algorithm |
| const double fatLine_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# fat line testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Bezier c( someCurves[i] ); |
| |
| c.p0.x += fatLine_xOffset; |
| c.p1.x += fatLine_xOffset; |
| c.p2.x += fatLine_xOffset; |
| c.p3.x += fatLine_xOffset; |
| |
| FatLine line; |
| |
| Impl_calcFatLine(line, c); |
| |
| cout << " bez(" |
| << c.p0.x << "," |
| << c.p1.x << "," |
| << c.p2.x << "," |
| << c.p3.x << ",t), bez(" |
| << c.p0.y << "," |
| << c.p1.y << "," |
| << c.p2.y << "," |
| << c.p3.y << ",t) title \"bezier " << i << "\", linex(" |
| << line.a << "," |
| << line.b << "," |
| << line.c << ",t), liney(" |
| << line.a << "," |
| << line.b << "," |
| << line.c << ",t) title \"fat line (center) on " << i << "\", linex(" |
| << line.a << "," |
| << line.b << "," |
| << line.c-line.dMin << ",t), liney(" |
| << line.a << "," |
| << line.b << "," |
| << line.c-line.dMin << ",t) title \"fat line (min) on " << i << "\", linex(" |
| << line.a << "," |
| << line.b << "," |
| << line.c-line.dMax << ",t), liney(" |
| << line.a << "," |
| << line.b << "," |
| << line.c-line.dMax << ",t) title \"fat line (max) on " << i << "\""; |
| |
| if( i+1<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| #endif |
| |
| #ifdef WITH_CALCFOCUS_TEST |
| // test focus curve algorithm |
| const double focus_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# focus line testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Bezier c( someCurves[i] ); |
| |
| c.p0.x += focus_xOffset; |
| c.p1.x += focus_xOffset; |
| c.p2.x += focus_xOffset; |
| c.p3.x += focus_xOffset; |
| |
| // calc focus curve |
| Bezier focus; |
| Impl_calcFocus(focus, c); |
| |
| cout << " bez(" |
| << c.p0.x << "," |
| << c.p1.x << "," |
| << c.p2.x << "," |
| << c.p3.x << ",t), bez(" |
| << c.p0.y << "," |
| << c.p1.y << "," |
| << c.p2.y << "," |
| << c.p3.y << ",t) title \"bezier " << i << "\", bez(" |
| << focus.p0.x << "," |
| << focus.p1.x << "," |
| << focus.p2.x << "," |
| << focus.p3.x << ",t), bez(" |
| << focus.p0.y << "," |
| << focus.p1.y << "," |
| << focus.p2.y << "," |
| << focus.p3.y << ",t) title \"focus " << i << "\""; |
| |
| |
| if( i+1<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| #endif |
| |
| #ifdef WITH_SAFEPARAMBASE_TEST |
| // test safe params base method |
| double safeParamsBase_xOffset( curr_Offset ); |
| cout << "# safe param base method testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Bezier c( someCurves[i] ); |
| |
| c.p0.x += safeParamsBase_xOffset; |
| c.p1.x += safeParamsBase_xOffset; |
| c.p2.x += safeParamsBase_xOffset; |
| c.p3.x += safeParamsBase_xOffset; |
| |
| Polygon2D poly(4); |
| poly[0] = c.p0; |
| poly[1] = c.p1; |
| poly[2] = c.p2; |
| poly[3] = c.p3; |
| |
| double t1, t2; |
| |
| bool bRet( Impl_calcSafeParams( t1, t2, poly, 0, 1 ) ); |
| |
| Polygon2D convHull( convexHull( poly ) ); |
| |
| cout << " bez(" |
| << poly[0].x << "," |
| << poly[1].x << "," |
| << poly[2].x << "," |
| << poly[3].x << ",t),bez(" |
| << poly[0].y << "," |
| << poly[1].y << "," |
| << poly[2].y << "," |
| << poly[3].y << ",t), " |
| << "t+" << safeParamsBase_xOffset << ", 0, " |
| << "t+" << safeParamsBase_xOffset << ", 1, "; |
| if( bRet ) |
| { |
| cout << t1+safeParamsBase_xOffset << ", t, " |
| << t2+safeParamsBase_xOffset << ", t, "; |
| } |
| cout << "'-' using ($1):($2) title \"control polygon\" with lp, " |
| << "'-' using ($1):($2) title \"convex hull\" with lp"; |
| |
| if( i+1<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| |
| safeParamsBase_xOffset += 2; |
| } |
| |
| safeParamsBase_xOffset = curr_Offset; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| Bezier c( someCurves[i] ); |
| |
| c.p0.x += safeParamsBase_xOffset; |
| c.p1.x += safeParamsBase_xOffset; |
| c.p2.x += safeParamsBase_xOffset; |
| c.p3.x += safeParamsBase_xOffset; |
| |
| Polygon2D poly(4); |
| poly[0] = c.p0; |
| poly[1] = c.p1; |
| poly[2] = c.p2; |
| poly[3] = c.p3; |
| |
| double t1, t2; |
| |
| Impl_calcSafeParams( t1, t2, poly, 0, 1 ); |
| |
| Polygon2D convHull( convexHull( poly ) ); |
| |
| unsigned int k; |
| for( k=0; k<poly.size(); ++k ) |
| { |
| cout << poly[k].x << " " << poly[k].y << endl; |
| } |
| cout << poly[0].x << " " << poly[0].y << endl; |
| cout << "e" << endl; |
| |
| for( k=0; k<convHull.size(); ++k ) |
| { |
| cout << convHull[k].x << " " << convHull[k].y << endl; |
| } |
| cout << convHull[0].x << " " << convHull[0].y << endl; |
| cout << "e" << endl; |
| |
| safeParamsBase_xOffset += 2; |
| } |
| curr_Offset += 20; |
| #endif |
| |
| #ifdef WITH_SAFEPARAMS_TEST |
| // test safe parameter range algorithm |
| const double safeParams_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# safe param range testing" << endl |
| << "plot [t=0.0:1.0] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j ) |
| { |
| Bezier c1( someCurves[i] ); |
| Bezier c2( someCurves[j] ); |
| |
| c1.p0.x += safeParams_xOffset; |
| c1.p1.x += safeParams_xOffset; |
| c1.p2.x += safeParams_xOffset; |
| c1.p3.x += safeParams_xOffset; |
| c2.p0.x += safeParams_xOffset; |
| c2.p1.x += safeParams_xOffset; |
| c2.p2.x += safeParams_xOffset; |
| c2.p3.x += safeParams_xOffset; |
| |
| double t1, t2; |
| |
| if( Impl_calcClipRange(t1, t2, c1, c1, c2, c2) ) |
| { |
| // clip safe ranges off c1 |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| // subdivide at t1_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 ); |
| // subdivide at t2_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) ); |
| |
| // output remaining segment (c1_part1) |
| |
| cout << " bez(" |
| << c1.p0.x << "," |
| << c1.p1.x << "," |
| << c1.p2.x << "," |
| << c1.p3.x << ",t),bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",t), bez(" |
| << c2.p0.x << "," |
| << c2.p1.x << "," |
| << c2.p2.x << "," |
| << c2.p3.x << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t), bez(" |
| << c1_part1.p0.x << "," |
| << c1_part1.p1.x << "," |
| << c1_part1.p2.x << "," |
| << c1_part1.p3.x << ",t),bez(" |
| << c1_part1.p0.y << "," |
| << c1_part1.p1.y << "," |
| << c1_part1.p2.y << "," |
| << c1_part1.p3.y << ",t)"; |
| |
| if( i+2<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| } |
| } |
| #endif |
| |
| #ifdef WITH_SAFEPARAM_DETAILED_TEST |
| // test safe parameter range algorithm |
| const double safeParams2_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| if( sizeof(someCurves)/sizeof(Bezier) > 1 ) |
| { |
| Bezier c1( someCurves[0] ); |
| Bezier c2( someCurves[1] ); |
| |
| c1.p0.x += safeParams2_xOffset; |
| c1.p1.x += safeParams2_xOffset; |
| c1.p2.x += safeParams2_xOffset; |
| c1.p3.x += safeParams2_xOffset; |
| c2.p0.x += safeParams2_xOffset; |
| c2.p1.x += safeParams2_xOffset; |
| c2.p2.x += safeParams2_xOffset; |
| c2.p3.x += safeParams2_xOffset; |
| |
| double t1, t2; |
| |
| // output happens here |
| Impl_calcClipRange(t1, t2, c1, c1, c2, c2); |
| } |
| #endif |
| |
| #ifdef WITH_SAFEFOCUSPARAM_TEST |
| // test safe parameter range from focus algorithm |
| const double safeParamsFocus_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << "# safe param range from focus testing" << endl |
| << "plot [t=0.0:1.0] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j ) |
| { |
| Bezier c1( someCurves[i] ); |
| Bezier c2( someCurves[j] ); |
| |
| c1.p0.x += safeParamsFocus_xOffset; |
| c1.p1.x += safeParamsFocus_xOffset; |
| c1.p2.x += safeParamsFocus_xOffset; |
| c1.p3.x += safeParamsFocus_xOffset; |
| c2.p0.x += safeParamsFocus_xOffset; |
| c2.p1.x += safeParamsFocus_xOffset; |
| c2.p2.x += safeParamsFocus_xOffset; |
| c2.p3.x += safeParamsFocus_xOffset; |
| |
| double t1, t2; |
| |
| Bezier focus; |
| #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS |
| #if 0 |
| { |
| // clip safe ranges off c1_orig |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| // subdivide at t1_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c2, 0.30204 ); |
| |
| // subdivide at t2_c1. As we're working on |
| // c1_part2 now, we have to adapt t2_c1 since |
| // we're no longer in the original parameter |
| // interval. This is based on the following |
| // assumption: t2_new = (t2-t1)/(1-t1), which |
| // relates the t2 value into the new parameter |
| // range [0,1] of c1_part2. |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (0.57151-0.30204)/(1.0-0.30204) ); |
| |
| c2 = c1_part1; |
| Impl_calcFocus( focus, c2 ); |
| } |
| #else |
| Impl_calcFocus( focus, c2 ); |
| #endif |
| #else |
| focus = c2; |
| #endif |
| // determine safe range on c1 |
| bool bRet( Impl_calcSafeParams_focus( t1, t2, |
| c1, focus ) ); |
| |
| cerr << "t1: " << t1 << ", t2: " << t2 << endl; |
| |
| // clip safe ranges off c1 |
| Bezier c1_part1; |
| Bezier c1_part2; |
| Bezier c1_part3; |
| |
| // subdivide at t1_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 ); |
| // subdivide at t2_c1 |
| Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) ); |
| |
| // output remaining segment (c1_part1) |
| |
| cout << " bez(" |
| << c1.p0.x << "," |
| << c1.p1.x << "," |
| << c1.p2.x << "," |
| << c1.p3.x << ",t),bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",t) title \"c1\", " |
| #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS |
| << "bez(" |
| << c2.p0.x << "," |
| << c2.p1.x << "," |
| << c2.p2.x << "," |
| << c2.p3.x << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t) title \"c2\", " |
| << "bez(" |
| << focus.p0.x << "," |
| << focus.p1.x << "," |
| << focus.p2.x << "," |
| << focus.p3.x << ",t),bez(" |
| << focus.p0.y << "," |
| << focus.p1.y << "," |
| << focus.p2.y << "," |
| << focus.p3.y << ",t) title \"focus\""; |
| #else |
| << "bez(" |
| << c2.p0.x << "," |
| << c2.p1.x << "," |
| << c2.p2.x << "," |
| << c2.p3.x << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t) title \"focus\""; |
| #endif |
| if( bRet ) |
| { |
| cout << ", bez(" |
| << c1_part1.p0.x << "," |
| << c1_part1.p1.x << "," |
| << c1_part1.p2.x << "," |
| << c1_part1.p3.x << ",t),bez(" |
| << c1_part1.p0.y+0.01 << "," |
| << c1_part1.p1.y+0.01 << "," |
| << c1_part1.p2.y+0.01 << "," |
| << c1_part1.p3.y+0.01 << ",t) title \"part\""; |
| } |
| |
| if( i+2<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| } |
| #endif |
| |
| #ifdef WITH_SAFEFOCUSPARAM_DETAILED_TEST |
| // test safe parameter range algorithm |
| const double safeParams3_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| if( sizeof(someCurves)/sizeof(Bezier) > 1 ) |
| { |
| Bezier c1( someCurves[0] ); |
| Bezier c2( someCurves[1] ); |
| |
| c1.p0.x += safeParams3_xOffset; |
| c1.p1.x += safeParams3_xOffset; |
| c1.p2.x += safeParams3_xOffset; |
| c1.p3.x += safeParams3_xOffset; |
| c2.p0.x += safeParams3_xOffset; |
| c2.p1.x += safeParams3_xOffset; |
| c2.p2.x += safeParams3_xOffset; |
| c2.p3.x += safeParams3_xOffset; |
| |
| double t1, t2; |
| |
| Bezier focus; |
| #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS |
| Impl_calcFocus( focus, c2 ); |
| #else |
| focus = c2; |
| #endif |
| |
| // determine safe range on c1, output happens here |
| Impl_calcSafeParams_focus( t1, t2, |
| c1, focus ); |
| } |
| #endif |
| |
| #ifdef WITH_BEZIERCLIP_TEST |
| ::std::vector< ::std::pair<double, double> > result; |
| ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(result); |
| |
| // test full bezier clipping |
| const double bezierClip_xOffset( curr_Offset ); |
| curr_Offset += 20; |
| cout << endl << endl << "# bezier clip testing" << endl |
| << "plot [t=0:1] "; |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j ) |
| { |
| Bezier c1( someCurves[i] ); |
| Bezier c2( someCurves[j] ); |
| |
| c1.p0.x += bezierClip_xOffset; |
| c1.p1.x += bezierClip_xOffset; |
| c1.p2.x += bezierClip_xOffset; |
| c1.p3.x += bezierClip_xOffset; |
| c2.p0.x += bezierClip_xOffset; |
| c2.p1.x += bezierClip_xOffset; |
| c2.p2.x += bezierClip_xOffset; |
| c2.p3.x += bezierClip_xOffset; |
| |
| cout << " bez(" |
| << c1.p0.x << "," |
| << c1.p1.x << "," |
| << c1.p2.x << "," |
| << c1.p3.x << ",t),bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",t), bez(" |
| << c2.p0.x << "," |
| << c2.p1.x << "," |
| << c2.p2.x << "," |
| << c2.p3.x << ",t),bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",t), '-' using (bez(" |
| << c1.p0.x << "," |
| << c1.p1.x << "," |
| << c1.p2.x << "," |
| << c1.p3.x |
| << ",$1)):(bez(" |
| << c1.p0.y << "," |
| << c1.p1.y << "," |
| << c1.p2.y << "," |
| << c1.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << i << ")\", " |
| << " '-' using (bez(" |
| << c2.p0.x << "," |
| << c2.p1.x << "," |
| << c2.p2.x << "," |
| << c2.p3.x |
| << ",$1)):(bez(" |
| << c2.p0.y << "," |
| << c2.p1.y << "," |
| << c2.p2.y << "," |
| << c2.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << j << ")\""; |
| |
| if( i+2<sizeof(someCurves)/sizeof(Bezier) ) |
| cout << ",\\" << endl; |
| else |
| cout << endl; |
| } |
| } |
| for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i ) |
| { |
| for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j ) |
| { |
| result.clear(); |
| Bezier c1( someCurves[i] ); |
| Bezier c2( someCurves[j] ); |
| |
| c1.p0.x += bezierClip_xOffset; |
| c1.p1.x += bezierClip_xOffset; |
| c1.p2.x += bezierClip_xOffset; |
| c1.p3.x += bezierClip_xOffset; |
| c2.p0.x += bezierClip_xOffset; |
| c2.p1.x += bezierClip_xOffset; |
| c2.p2.x += bezierClip_xOffset; |
| c2.p3.x += bezierClip_xOffset; |
| |
| clipBezier( ii, 0.00001, c1, c2 ); |
| |
| for( k=0; k<result.size(); ++k ) |
| { |
| cout << result[k].first << endl; |
| } |
| cout << "e" << endl; |
| |
| for( k=0; k<result.size(); ++k ) |
| { |
| cout << result[k].second << endl; |
| } |
| cout << "e" << endl; |
| } |
| } |
| #endif |
| |
| return 0; |
| } |
