// $Id: pathutil.cpp 1282 2006-06-09 09:46:49Z alex $ /* @@tag:xara-cn@@ DO NOT MODIFY THIS LINE ================================XARAHEADERSTART=========================== Xara LX, a vector drawing and manipulation program. Copyright (C) 1993-2006 Xara Group Ltd. Copyright on certain contributions may be held in joint with their respective authors. See AUTHORS file for details. LICENSE TO USE AND MODIFY SOFTWARE ---------------------------------- This file is part of Xara LX. Xara LX is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 2 as published by the Free Software Foundation. Xara LX and its component source files are distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Xara LX (see the file GPL in the root directory of the distribution); if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ADDITIONAL RIGHTS ----------------- Conditional upon your continuing compliance with the GNU General Public License described above, Xara Group Ltd grants to you certain additional rights. The additional rights are to use, modify, and distribute the software together with the wxWidgets library, the wxXtra library, and the "CDraw" library and any other such library that any version of Xara LX relased by Xara Group Ltd requires in order to compile and execute, including the static linking of that library to XaraLX. In the case of the "CDraw" library, you may satisfy obligation under the GNU General Public License to provide source code by providing a binary copy of the library concerned and a copy of the license accompanying it. Nothing in this section restricts any of the rights you have under the GNU General Public License. SCOPE OF LICENSE ---------------- This license applies to this program (XaraLX) and its constituent source files only, and does not necessarily apply to other Xara products which may in part share the same code base, and are subject to their own licensing terms. This license does not apply to files in the wxXtra directory, which are built into a separate library, and are subject to the wxWindows license contained within that directory in the file "WXXTRA-LICENSE". This license does not apply to the binary libraries (if any) within the "libs" directory, which are subject to a separate license contained within that directory in the file "LIBS-LICENSE". ARRANGEMENTS FOR CONTRIBUTION OF MODIFICATIONS ---------------------------------------------- Subject to the terms of the GNU Public License (see above), you are free to do whatever you like with your modifications. However, you may (at your option) wish contribute them to Xara's source tree. You can find details of how to do this at: http://www.xaraxtreme.org/developers/ Prior to contributing your modifications, you will need to complete our contributor agreement. This can be found at: http://www.xaraxtreme.org/developers/contribute/ Please note that Xara will not accept modifications which modify any of the text between the start and end of this header (marked XARAHEADERSTART and XARAHEADEREND). MARKS ----- Xara, Xara LX, Xara X, Xara X/Xtreme, Xara Xtreme, the Xtreme and Xara designs are registered or unregistered trademarks, design-marks, and/or service marks of Xara Group Ltd. All rights in these marks are reserved. Xara Group Ltd, Gaddesden Place, Hemel Hempstead, HP2 6EX, UK. http://www.xara.com/ =================================XARAHEADEREND============================ */ // Implementation file for functions which operate on sub path elements // ie lines and curves /* */ #include "camtypes.h" //#include "paths.h" - in camtypes.h [AUTOMATICALLY REMOVED] #include "pathutil.h" #include //#include "fixmem.h" - in camtypes.h [AUTOMATICALLY REMOVED] #include "macros.h" #define EPSILON (ldexp(1.0,-MAXDEPTH-1)) // Flatness control value #define DEGREE 3 // Cubic Bezier curve #define W_DEGREE 5 // Degree of eqn to find roots of const INT32 MAXDEPTH = 64; // Maximum depth for recursion /******************************************************************************************** > DocCoord RampItem::PointOnSemiCircle(const DocCoord& centre, const DocCoord& radialp, double t) Author: Mike_Kenny (Xara Group Ltd) Created: 05/03/97 Inputs: centre = the centre point of the circle radialp = another point (anywhere, but assumed to be somewhere on the radius of the circle. t = a parameter (0..1) defining the parametric locus about the centre from 0 to pi radians. Outputs: - Returns: A new doc coord which is a point on the circles perimeter which corresponds to the parameter t. At t==1 the function will evaluate to radialp At t==0 the function will evaluate to centre-(radialp-centre) Purpose: Find a point on the circumference of a semicircle. The semicircle is specified by two points, it's centre and a point on the circumference. If c=0,0 and p=(1,0) then the semicircle exists in the y positive half of the plane from sweeping from (-1,0) at t==0 to (1.0) at t==1.0 Errors: None. ********************************************************************************************/ DocCoord PathUtil::PointOnSemiCircle(const DocCoord& centre, const DocCoord& radialp, double t) { double X = (double)(radialp.x - centre.x); double Y = (double)(radialp.y - centre.y); double p = t * XS_PI; // spin clockwise. double s = sin(p); double c = -cos(p); DocCoord r; // circular rounding maybe? r.x = (INT32)(X*c - Y*s) + centre.x; r.y = (INT32)(X*s + Y*c) + centre.y; return r; } /******************************************************************************************** > double PathUtil::SqrDistanceToSemiCircle(const DocCoord* plist, const DocCoord& p1, double* param ) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: plist = a pointer to two doc coordinates describing a half circle p1 = the point to find the squared distance to Outputs: param = the parameter at which the closest point exists Returns: a double, the distance a given point is away from a line element Purpose: Calculates the distance p1 is away from a semi circle segment. plist[0] describes the centre of the semi circle plist[1] describes a point on the circumference of the circle. For instance say c the centre is (0,0) and cp the circumference point is (1,0) then the semicircle exists in the positive y half of the place and sweeps from (-1,0) to (1,0). Parameter space is param==0 at (-1,0) param==1 at ( 1,0) param==0.5 at ( 0,1) ********************************************************************************************/ double PathUtil::SqrDistanceToSemiCircle(const DocCoord* plist, const DocCoord& p, double *param) { double ex,ey,px,py,Px,Py,l,t; DocCoord c = plist[0]; // the centre point DocCoord e = plist[1]; // the end point DocCoord q; // translate to origin ex = e.x-c.x; ey = e.y-c.y; px = p.x-c.x; py = p.y-c.y; // build a (ex,ey) based transform l = sqrt(ex*ex+ey*ey); if (l!=0) { ex=ex/l; ey=ex/l; } // transform the click point to local canonical Px = px*ex+py*ey; Py = py*ex-px*ey; // are we below the half circle if (Py>0) { // no then calculate the dot product angle l = sqrt(Px*Px+Py*Py); if (l!=0) Px=Px/l; t = (XS_PI - acos(Px))/XS_PI; q = PathUtil::PointOnSemiCircle(c,p,t); } else { // yes then decide start point or endpoint if (Px>0) { t=1.0; q=e; } else { t=0.0; q=e-(c-e); } } // set the output parameter (*param)=t; // calculate the squared distance return PathUtil::SqrDistance(p,q); } /******************************************************************************************** > double PathUtil::SqrDistance(const Coord& Coord) Author: Mike_Kenny (Xara Group Ltd) Created: 20/8/94 Inputs: two coordinates Outputs: None Returns: the squared distance between the two coordinates Purpose: Accurate squared distance function. Errors: None. ********************************************************************************************/ double PathUtil::SqrDistance(const DocCoord& p1, const DocCoord& p2) { double dx = (double) (p2.x - p1.x); double dy = (double) (p2.y - p1.y); return ((dx*dx)+(dy*dy)); } /******************************************************************************************** > BOOL PathUtil::SplitLine(const double t, const DocCoord* plist, UINT32* NumElements, PathVerb* OutVerbs, DocCoord* OutCoords) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: t = parameter at which to split the line plist = pointer to 2 doc coords describing the line Outputs: NumElements = number of elements generated after split OutVerbs = pointer to verbs list to output data to OutCoords = pointer to coords list to output data to Returns: True if the line can be split, False if not. Purpose: Splits a line element into two lines, returning the lists of new coord points and verbs. ********************************************************************************************/ /* Technical notes: Two lineto coordinates will be returned as follows Verbs Coords Lineto OA Lineto OB = IA and can be inserted as follows Input Output Verbs Coords Verbs Coords .... IA .... IA Lineto IB Lineto OA ... XX Lineto OB ... XX OB is returned in the output list in order to allow removal of the split element and complete replacement with the generated elements ********************************************************************************************/ BOOL PathUtil::SplitLine(const double t, const DocCoord* plist, UINT32* NumElements, PathVerb* Verbs, DocCoord* Coords) { // check t is in range for this split. if (t < SPLIT_EPSILON) return FALSE; if (t > (1.0 - SPLIT_EPSILON)) return FALSE; // read the lines start and end points INT32 x0,y0,x1,y1; x0 = plist[0].x; y0 = plist[0].y; x1 = plist[1].x; y1 = plist[1].y; // fill in the output block details Verbs[0] = PT_LINETO; Coords[0].x = x0 + (INT32)(t*(x1-x0)+0.5); Coords[0].y = y0 + (INT32)(t*(y1-y0)+0.5); Verbs[1] = PT_LINETO; Coords[1].x = x1; Coords[1].y = y1; *NumElements = 2; return TRUE; } /******************************************************************************************** > BOOL PathUtil::SplitCurve(const double t, const DocCoord* plist, UINT32* NumElements, PathVerb* OutVerbs, DocCoord* OutCoords) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: t = parameter at which to split the curve plist = pointer to 4 curve control points to split Outputs: NumElements = number of elements generated OutVerbs = output verbs list OutCoords = output coords list Returns: True if the curve could be split, False if not. Purpose: Splits a curve element into two curves, returning the lists of new control points and verbs. ********************************************************************************************/ /* Technical notes: Six bezier control points will be returned as follows: Input Output IA OA IB OB IC OC ID OD OE OF = ID The resultant curve control lists should be encorperated by the caller as follows: Input Output Verbs Coords Verbs Coords .... IA .... IA Bez IB Bez OA Bez IC Bez OB Bez ID Bez OC ... XX Bez OD ... XX Bez OE Bez OF ... XX ... XX ********************************************************************************************/ BOOL PathUtil::SplitCurve(const double t, const DocCoord* plist, UINT32* NumElements, PathVerb* OutVerbs, DocCoord* OutCoords) { // check t is in range for this split. if (t < SPLIT_EPSILON) return FALSE; if (t > (1.0 - SPLIT_EPSILON)) return FALSE; // fill in the path verb array, these are all curveto's for (INT32 i=0; i<6; i++) OutVerbs[i] = PT_BEZIERTO; // read the four curve control points INT32 x0,y0,x1,y1,x2,y2,x3,y3; x0 = plist[0].x; // note this should be move,line,curve y0 = plist[0].y; x1 = plist[1].x; // these should all be curveto's y1 = plist[1].y; x2 = plist[2].x; y2 = plist[2].y; x3 = plist[3].x; y3 = plist[3].y; // now calculate six new coordinates from the given 4 double tx,ty; tx = x1+t*(x2-x1); ty = y1+t*(y2-y1); double Lx1,Ly1,Lx2,Ly2,Rx0,Ry0,Rx1,Ry1,Rx2,Ry2; Lx1 = x0+t*(x1-x0); Ly1 = y0+t*(y1-y0); Rx2 = x2+t*(x3-x2); Ry2 = y2+t*(y3-y2); Lx2 = Lx1+t*(tx-Lx1); Ly2 = Ly1+t*(ty-Ly1); Rx1 = tx+t*(Rx2-tx); Ry1 = ty+t*(Ry2-ty); Rx0 = Lx2+t*(Rx1-Lx2); Ry0 = Ly2+t*(Ry1-Ly2); // set the return values OutCoords[0].x = (INT32)(Lx1+0.5); OutCoords[0].y = (INT32)(Ly1+0.5); OutCoords[1].x = (INT32)(Lx2+0.5); OutCoords[1].y = (INT32)(Ly2+0.5); OutCoords[2].x = (INT32)(Rx0+0.5); OutCoords[2].y = (INT32)(Ry0+0.5); OutCoords[3].x = (INT32)(Rx1+0.5); OutCoords[3].y = (INT32)(Ry1+0.5); OutCoords[4].x = (INT32)(Rx2+0.5); OutCoords[4].y = (INT32)(Ry2+0.5); OutCoords[5].x = x3; OutCoords[5].y = y3; // Set the output number of elements. *NumElements = 6; return TRUE; } /******************************************************************************************** > void PathUtil::GetCurveCoefs(const DocCoord* coords, PtCoefs* coefs ) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: coords = a pointer to 4 control point doccoords defining a curve Outputs: coefs = a set of curve coefficients. Returns: None Purpose: Converts a curve from Bezier form to Canonical form ********************************************************************************************/ void PathUtil::GetCurveCoefs(const DocCoord* coords, PtCoefs* coefs ) { // Read the curve coordinates. INT32 X0,Y0,X1,Y1,X2,Y2,X3,Y3; X0 = coords->x; Y0 = coords->y; coords++; X1 = coords->x; Y1 = coords->y; coords++; X2 = coords->x; Y2 = coords->y; coords++; X3 = coords->x; Y3 = coords->y; // Calculate the curve coefficients. coefs->ax = X3-X0+3*(X1-X2); coefs->ay = Y3-Y0+3*(Y1-Y2); coefs->bx = 3*(X2-2*X1+X0); coefs->by = 3*(Y2-2*Y1+Y0); coefs->cx = 3*(X1-X0); coefs->cy = 3*(Y1-Y0); } /******************************************************************************************** > double PathUtil::SqrDistanceToLine(const DocCoord* plist, const DocCoord& p1, double* t ) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: plist = a pointer to two doc coordinates describing a line p1 = the point to find the squared distance to Outputs: t = the parameter at which the closest point exists Returns: a double, the distance a given point is away from a line element Purpose: Calculates the distance p1 is away from a line segment. The perpendicular distance is returned only when p1 is within the volume created by sweeping the line in the orthoganal direction. Otherwise the distance to the closest end point is returned. ********************************************************************************************/ double PathUtil::SqrDistanceToLine(const DocCoord* plist, const DocCoord& p1, double* t ) { // get hold of the lines end point coordinates INT32 x0,y0,x1,y1; x0 = plist[0].x; y0 = plist[0].y; x1 = plist[1].x; y1 = plist[1].y; INT32 pdx = (x1-x0); INT32 pdy = (y1-y0); // calculate the parameter at which the closest point exists // on an infinite line passing through (x0,y0),(x1,y1) // t=numdot/dendot. double numdot,dendot; numdot = (double)(p1.x - plist[0].x)*pdx + (double)(p1.y-plist[0].y)*pdy; dendot = (double)pdx*pdx + (double)pdy*pdy; // if t is less then zero then the closest point lies behind // the start point of the vector, hence return the squared // distance to the startpoint if (numdot<=0 || dendot<=0) { *t = 0.0; return (SqrDistance(p1, plist[0])); } // if t is greater than or equal to 1 then the closest point // lies beyond the line segment, ie on the projected line // beyond (x1,y1). So return the squared distance to that end // point. if (dendot<=numdot) { *t = 1.0; return (SqrDistance(p1,plist[1])); } // ok the closest point lies somewhere inside the line segment // (x0,y0),(x1,y1) so return its parameter and the squared // distance to it. *t = numdot/dendot; double c = (double)pdy*p1.x - (double)pdx*p1.y; double d = (double)y1*x0 - (double)y0*x1; double e = (c-d)*(c-d); return fabs(e/dendot); } /******************************************************************************************** > DocCoord PathUtil::PointOnLine(const double t, const DocCoord* startpt) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: t = parameter at which to evaluate point ptlist = pointer to two doc coords describing start and end of the line Outputs: Returns: DocCoord, an evaluation of a point on the specified line Purpose: Given a parameter t go and evaluate the point on the described line segment If t is out of range, ie t>1 or t<0 the end points of the line will be returned. ie the routine will not evaluate unbounded points. ********************************************************************************************/ DocCoord PathUtil::PointOnLine(const double t, const DocCoord* ptlist) { INT32 x0,y0,x1,y1; x0 = ptlist[0].x; y0 = ptlist[0].y; x1 = ptlist[1].x; y1 = ptlist[1].y; DocCoord d; if (t<=0.0) { d.x = x0; d.y = y0; return d; } if (t>=1.0) { d.x = x1; d.y = y1; return d; } d.x = x0 + (INT32)(t*(x1-x0)+0.5); d.y = y0 + (INT32)(t*(y1-y0)+0.5); return d; } /******************************************************************************************** > double PathUtil::SqrDistanceToCurve(const UINT32 step, const DocCoord* plist, const DocCoord& p1, double* t) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: step = value between 2 and 256 dictates the corseness of stepping internally across the curve. plist = pointer to 4 doccoord control points describing the curve p1 = coordinate reference used to find closest coordinate to. Outputs: t holds the parameter value of the closest point on the curve, 0<=t<=1 Returns: a double being the distance of the closest point on the curve. Purpose: This function returns the distance to the closest point on the curve and the parameter of this point ********************************************************************************************/ /* Technical notes: If Qx(t) = X0*(1-t)^3 + 3*X1*t*(1-t)^2 + 3*X2*t^2*(1-t) + X3*t^3 and Qy(t) is similar Then Qx(t) = (X0 - 3*X0*t + 3*X0*t^2 - X0*t^3) + (3*X1*t - 6*X1*t^2 + 3*X1*t^3 ) + (3*X2*t^2 - 3*X2*t^3) + X3*t^3 Collecting terms in t Qx(t) = t^3*(X3-X0+3*(X1-X2)) + t^2*(3*X0-6*X1+3*X2) + t *(3*X1-3*X0) + X0 For forward differencing Let Ax = X3-X0+3*(X1-X2) Bx = 3*(X0-2*X1+X2) Cx = 3*(X1-X0) Dx = X0 Then Q(t) = A*t^3 + B*t^2 + C*t + D = t*(C+t*(B+t*A)) + D Ok, having got this calculate the terms involved in D(t) = Q(t+delta)-Q(t) The terms involved in D(t) will give forward differences which can be used to walk the curve Q(t) ********************************************************************************************/ /* double PathUtil::SqrDistanceToCurve(const UINT32 step, const DocCoord* plist, const DocCoord& p1, double* t ) { ENSURE(step>1 && step<257,"curve step is out of range"); // Calculate the curve coefficients for this curve element PtCoefs CurveC; PathUtil::GetCurveCoefs(plist, &CurveC); double t0; double t1; double dist = PathUtil::CurveClosestRange(step, &CurveC, p1, plist[0].x, plist[0].y, t, &t0, &t1); // TRACE( _T("T0 = %f\n"),t0); // TRACE( _T("T = %f\n"),*t); // TRACE( _T("T1 = %f\n"),t1); if (t0 >= t1) return dist; if ((*t==t0) || (*t==t1)) return dist; if ((*t==0.0) || (*t==1.0)) return dist; INT32 x0,y0,x1,y1,x2,y2,x3,y3; x0 = plist[0].x; // note this should be move,line,curve y0 = plist[0].y; x1 = plist[1].x; // these should all be curveto's y1 = plist[1].y; x2 = plist[2].x; y2 = plist[2].y; x3 = plist[3].x; y3 = plist[3].y; // now calculate the coordinates of a bezier curve // lying between t0 and t1 double tx,ty; tx = x1+t0*(x2-x1); ty = y1+t0*(y2-y1); double Qx0,Qy0,Qx1,Qy1; double Rx1,Ry1,Rx2,Ry2,Lx1,Ly1,Lx2,Ly2; Lx1 = x0+t0*(x1-x0); Ly1 = y0+t0*(y1-y0); Rx2 = x2+t0*(x3-x2); Ry2 = y2+t0*(y3-y2); Lx2 = Lx1+t0*(tx-Lx1); Ly2 = Ly1+t0*(ty-Ly1); Qx1 = tx+t0*(Rx2-tx); Qy1 = ty+t0*(Ry2-ty); Qx0 = Lx2+t0*(Qx1-Lx2); Qy0 = Ly2+t0*(Qy1-Ly2); tx = x1+t1*(x2-x1); ty = y1+t1*(y2-y1); double Qx2,Qy2,Qx3,Qy3; Lx1 = x0+t1*(x1-x0); Ly1 = y0+t1*(y1-y0); Rx1 = x2+t1*(x3-x2); Ry1 = y2+t1*(y3-y2); Rx2 = tx+t1*(Rx1-tx); Ry2 = ty+t1*(Ry1-ty); Qx2 = Lx1+t1*(tx-Lx1); Qy2 = Ly1+t1*(ty-Ly1); Qx3 = Qx2+t1*(Rx2-Qx2); Qy3 = Qy2+t1*(Ry2-Qy2); // now home in a little bit closer CurveC.ax = Qx3-Qx0+3*(Qx1-Qx2); CurveC.ay = Qy3-Qy0+3*(Qy1-Qy2); CurveC.bx = 3*(Qx2-2*Qx1+Qx0); CurveC.by = 3*(Qy2-2*Qy1+Qy0); CurveC.cx = 3*(Qx1-Qx0); CurveC.cy = 3*(Qy1-Qy0); double t2; double t3; dist = PathUtil::CurveClosestRange(step, &CurveC, p1, Qx0, Qy0, t, &t2, &t3); *t = (1.0-(*t))*t0 + (*t)*t1; // TRACE( _T("T' = %f\n\n"),*t); return dist; } double PathUtil::CurveClosestRange(const UINT32 step, PtCoefs* CurveC, const DocCoord& p1, const double dx, const double dy, double* t, double* t0, double* t1) { // Calculate the step rate across the curve given the users // step value double Delta=(double) 1/step; // Calculate the forward difference factors for curve stepping double X0,Y0,X1,Y1,X2,Y2,X3,Y3; X3 = 6*Delta*Delta*Delta*CurveC->ax; Y3 = 6*Delta*Delta*Delta*CurveC->ay; X2 = 2*Delta*Delta*CurveC->bx+X3; Y2 = 2*Delta*Delta*CurveC->by+Y3; // Evaluate the first differenced point on the curve X1 = Delta*(CurveC->cx+Delta*(CurveC->bx+Delta*CurveC->ax)); Y1 = Delta*(CurveC->cy+Delta*(CurveC->by+Delta*CurveC->ay)); X0 = dx; Y0 = dy; // find the distance from the first point on the curve double Cdist=(p1.x - dx)*(p1.x - dx) + (p1.y - dy)*(p1.y - dy); double Fdist; INT32 Cpt=0; for (UINT32 i=1; i<=step; i++) { X0+=X1; Y0+=Y1; X1+=X2; Y1+=Y2; X2+=X3; Y2+=Y3; Fdist = (X0-p1.x)*(X0-p1.x); // Is this point closer? if (Cdist>Fdist) { Fdist += (Y0-p1.y)*(Y0-p1.y); if (Cdist>Fdist) { Cdist=Fdist; Cpt=i; } } } INT32 tl = Cpt-2; UINT32 tr = Cpt+2; if (tl<0) tl=0; if (tr>step) tr=step; *t1 = (double) tr/step; *t0 = (double) tl/step; *t = (double) Cpt/step; // set the point parameter return Cdist; // return the point distance } */ /******************************************************************************************** > DocCoord PathUtil::PointOnCurve(double t, const DocCoord* plist) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: t = parameter to evaluate curve at plist = pointer to list of 4 curve control points. Outputs: Returns: A DocCoord describing the point on the curve Purpose: Evaluate a bezier curve, given a pointer to a set of control points and a parameter value. ********************************************************************************************/ DocCoord PathUtil::PointOnCurve(double t, const DocCoord* plist) { PtCoefs p; DocCoord c; PathUtil::GetCurveCoefs(plist, &p); // calculate coefs block if (t<0) t=0; // clamp range of parameter if (t>1) t=1; // Eval bezier at specified parameter double dx = t*(p.cx+t*(p.bx+t*p.ax)); double dy = t*(p.cy+t*(p.by+t*p.ay)); c.x = plist->x + (INT32)(dx+0.5); c.y = plist->y + (INT32)(dy+0.5); return c; } /******************************************************************************************** Quick local vector functions to calculate various properties of vectors ********************************************************************************************/ Point2 V2ScaleII( Point2 *v, double s) { Point2 result; result.x = v->x * s; result.y = v->y * s; return (result); } double V2SquaredLength(Point2* a) { return( a->x*a->x + a->y*a->y ); } Point2 *V2Sub(Point2* a, Point2* b, Point2* c) { c->x = a->x - b->x; c->y = a->y - b->y; return(c); } double V2Dot(Point2* a, Point2* b) { return (a->x*b->x + a->y*b->y); } /******************************************************************************************** Point2 Bezier( Point2 *V, INT32 degree, double t, Point2 *Left, Point2 *Right) Author: Unattributed (Xara Group Ltd) Created: 22/08/94 Inputs: V = Control points of cubic Bezier degree = The degree of the polynomial t = Parameter value Outputs: Left = RETURN left half ctl pts (if NULL return none) Right = RETURN right half ctl pts (if NULL return none) Returns: Q(t), a point on the curve Purpose: Evaluate a Bezier curve at a particular parameter value Fill in control points for resulting sub-curves if "Left" and "Right" are non-null. ********************************************************************************************/ Point2 Bezier( Point2 *V, INT32 degree, double t, Point2 *Left, Point2 *Right) { INT32 i, j; // Index variables Point2 Vtemp[W_DEGREE+1][W_DEGREE+1]; // Copy control points for (j =0; j <= degree; j++) { Vtemp[0][j] = V[j]; } // Triangle computation for (i = 1; i <= degree; i++) { for (j =0 ; j <= degree - i; j++) { Vtemp[i][j].x = (1.0 - t) * Vtemp[i-1][j].x + t * Vtemp[i-1][j+1].x; Vtemp[i][j].y = (1.0 - t) * Vtemp[i-1][j].y + t * Vtemp[i-1][j+1].y; } } if (Left != NULL) { for (j = 0; j <= degree; j++) { Left[j] = Vtemp[j][0]; } } if (Right != NULL) { for (j = 0; j <= degree; j++) { Right[j] = Vtemp[degree-j][j]; } } return (Vtemp[degree][0]); } /******************************************************************************************** double ComputeXIntercept(Point2 *V, INT32 degree) Author: Created: 22/08/94 Inputs: V = Control points of cubic Bezier degree = The degree of the polynomial Outputs: Returns: Purpose: Compute intersection of chord from first control point to last with 0-axis. ********************************************************************************************/ double ComputeXIntercept(Point2 *V, INT32 degree) { double XLK, YLK, XNM, YNM, XMK, YMK; double det, detInv; double S, T; double X, Y; XLK = 1.0 - 0.0; YLK = 0.0 - 0.0; XNM = V[degree].x - V[0].x; YNM = V[degree].y - V[0].y; XMK = V[0].x - 0.0; YMK = V[0].y - 0.0; det = XNM*YLK - YNM*XLK; detInv = 1.0/det; S = (XNM*YMK - YNM*XMK) * detInv; T = (XLK*YMK - YLK*XMK) * detInv; X = 0.0 + XLK * S; Y = 0.0 + YLK * S; return X; } /******************************************************************************************** INT32 CrossingCount( Point2 *V, INT32 degree) Author: Created: 22/08/94 Inputs: V = Control points of cubic Bezier degree = The degree of the polynomial Outputs: Returns: number of crossings Purpose: Count the number of times a Bezier control polygon crosses the 0-axis. This number is >= the number of roots. *********************************************************************************************/ INT32 CrossingCount( Point2 *V, INT32 degree) { INT32 i; INT32 n_crossings = 0; // Number of zero-crossings INT32 sign, old_sign; // Sign of coefficients sign = old_sign = SGN(V[0].y); for (i = 1; i <= degree; i++) { sign = SGN(V[i].y); if (sign != old_sign) n_crossings++; old_sign = sign; } return n_crossings; } /******************************************************************************************** INT32 ControlPolygonFlatEnough(Point2 *V, INT32 degree) Author: Created: 22/08/94 Inputs: V = Control points of cubic Bezier degree = The degree of the polynomial Outputs: Returns: Purpose: Check if the control polygon of a Bezier curve is flat enough for recursive subdivision to bottom out. ********************************************************************************************/ INT32 ControlPolygonFlatEnough(Point2 *V, INT32 degree) { INT32 i; // Index variable double *distance; // Distances from pts to line double max_distance_above; // maximum of these double max_distance_below; double error; // Precision of root double intercept_1, intercept_2, left_intercept, right_intercept; double a, b, c; // Coefficients of implicit // eqn for line from V[0]-V[deg] /* Find the perpendicular distance from each interior control point to line connecting V[0] and V[degree] */ distance = (double *)CCMalloc((unsigned)(degree + 1) * sizeof(double)); { double abSquared; // Derive the implicit equation for line connecting first // and last control points a = V[0].y - V[degree].y; b = V[degree].x - V[0].x; c = V[0].x * V[degree].y - V[degree].x * V[0].y; abSquared = (a * a) + (b * b); for (i = 1; i < degree; i++) { // Compute distance from each of the points to that line distance[i] = a * V[i].x + b * V[i].y + c; if (distance[i] > 0.0) { distance[i] = (distance[i] * distance[i]) / abSquared; } if (distance[i] < 0.0) { distance[i] = -((distance[i] * distance[i]) / abSquared); } } } // Find the largest distance max_distance_above = 0.0; max_distance_below = 0.0; for (i = 1; i < degree; i++) { if (distance[i] < 0.0) { max_distance_below = MIN(max_distance_below, distance[i]); }; if (distance[i] > 0.0) { max_distance_above = MAX(max_distance_above, distance[i]); } } CCFree((char *)distance); { double det, dInv; double a1, b1, c1, a2, b2, c2; // Implicit equation for zero line a1 = 0.0; b1 = 1.0; c1 = 0.0; // Implicit equation for "above" line a2 = a; b2 = b; c2 = c + max_distance_above; det = a1 * b2 - a2 * b1; dInv = 1.0/det; intercept_1 = (b1 * c2 - b2 * c1) * dInv; // Implicit equation for "below" line a2 = a; b2 = b; c2 = c + max_distance_below; det = a1 * b2 - a2 * b1; dInv = 1.0/det; intercept_2 = (b1 * c2 - b2 * c1) * dInv; } // Compute intercepts of bounding box left_intercept = MIN(intercept_1, intercept_2); right_intercept = MAX(intercept_1, intercept_2); error = 0.5 * (right_intercept-left_intercept); if (error < EPSILON) { return 1; } else { return 0; } } /******************************************************************************************** INT32 FindRoots( Point2* w, INT32 degree, double* t, INT32 depth) Author: Created: 22/08/94 Inputs: w = Control points of cubic Bezier degree = The degree of the polynomial depth = The depth of the recursion Outputs: t = a list of candidate t-values Returns: Purpose: Given a 5th-degree equation in Bernstein-Bezier form, find all of the roots in the interval [0, 1]. Return the number of roots found. ********************************************************************************************/ INT32 FindRoots( Point2* w, INT32 degree, double* t, INT32 depth) { INT32 i; Point2 Left[W_DEGREE+1], // New left and right Right[W_DEGREE+1]; // control polygons INT32 left_count, // Solution count from children right_count; double left_t[W_DEGREE+1], // Solutions from kids right_t[W_DEGREE+1]; switch (CrossingCount(w, degree)) { case 0 : { // No solutions here return 0; break; } case 1 : { // Unique solution // Stop recursion when the tree is deep enough // if deep enough, return 1 solution at midpoint if (depth >= MAXDEPTH) { t[0] = (w[0].x + w[W_DEGREE].x) / 2.0; return 1; } if (ControlPolygonFlatEnough(w, degree)) { t[0] = ComputeXIntercept(w, degree); return 1; } break; } } // Otherwise, solve recursively after // subdividing control polygon Bezier(w, degree, 0.5, Left, Right); left_count = FindRoots(Left, degree, left_t, depth+1); right_count = FindRoots(Right, degree, right_t, depth+1); /* Gather solutions together */ for (i = 0; i < left_count; i++) { t[i] = left_t[i]; } for (i = 0; i < right_count; i++) { t[i+left_count] = right_t[i]; } /* Send back total number of solutions */ return (left_count+right_count); } /******************************************************************************************** Point2 *ConvertToBezierForm( Point2 P, Point2* V ) Author: Created: 22/08/94 Inputs: P = The point to find t for V = Control points of cubic Bezier Outputs: Returns: Purpose: Given a point and a Bezier curve, generate a 5th-degree Bezier-format equation whose solution finds the point on the curve nearest the user-defined point. *******************************************************************************************/ Point2 *ConvertToBezierForm( Point2 P, Point2* V ) { INT32 i, j, k, m, n, ub, lb; INT32 row, column; // Table indices Vector2 c[DEGREE+1]; // V(i)'s - P Vector2 d[DEGREE]; // V(i+1) - V(i) Point2 *w; // Ctrl pts of 5th-degree curve double cdTable[3][4]; // Dot product of c, d static double z[3][4] = { // Precomputed "z" for cubics {1.0, 0.6, 0.3, 0.1}, {0.4, 0.6, 0.6, 0.4}, {0.1, 0.3, 0.6, 1.0}, }; // Determine the c's -- these are vectors created by subtracting // point P from each of the control points for (i = 0; i <= DEGREE; i++) { V2Sub(&V[i], &P, &c[i]); } // Determine the d's -- these are vectors created by subtracting // each control point from the next for (i = 0; i <= DEGREE - 1; i++) { d[i] = V2ScaleII(V2Sub(&V[i+1], &V[i], &d[i]), 3.0); } // Create the c,d table -- this is a table of dot products of the // c's and d's for (row = 0; row <= DEGREE - 1; row++) { for (column = 0; column <= DEGREE; column++) { cdTable[row][column] = V2Dot(&d[row], &c[column]); } } // Now, apply the z's to the dot products, on the skew diagonal // Also, set up the x-values, making these "points" w = (Point2 *)CCMalloc((unsigned)(W_DEGREE+1) * sizeof(Point2)); for (i = 0; i <= W_DEGREE; i++) { w[i].y = 0.0; w[i].x = (double)(i) / W_DEGREE; } n = DEGREE; m = DEGREE-1; for (k = 0; k <= n + m; k++) { lb = MAX(0, k - m); ub = MIN(k, n); for (i = lb; i <= ub; i++) { j = k - i; w[i+j].y += cdTable[j][i] * z[j][i]; } } return (w); } /******************************************************************************************* > double PathUtil::SqrDistanceToCurve(DocCoord* V, DocCoord& P, double* mu ) Author: Mike_Kenny (Xara Group Ltd) Created: 22/08/94 Inputs: split = max level at which to split the curve (usually 64) P = The user-supplied point V = Control points of cubic Bezier Outputs: mu holds the parameter value of the closest point on the curve, 0<=mu<=1 Returns: The squared distance to the curve Purpose: Compute the parameter value of the point on a Bezier curve segment closest to some arbtitrary, user-input point. Return the squared distance to the curve. *******************************************************************************************/ double PathUtil::SqrDistanceToCurve(const UINT32 split, const DocCoord* V, const DocCoord& P, double* mu ) { Point2 *w; // Ctrl pts for 5th-degree equ Point2 lcurve[DEGREE+1]; // Local double version of control polygon double t_candidate[W_DEGREE]; // Possible roots INT32 n_solutions; // Number of roots found double t; // Parameter value of closest pt double dist; // Closest squared distance to curve // Convert control polygon to double type for (INT32 i=0; i < DEGREE+1; i++) { lcurve[i].x = (double) V[i].x; lcurve[i].y = (double) V[i].y; } Point2 p; p.x = (double) P.x; p.y = (double) P.y; // Convert problem to 5th-degree Bezier form w = ConvertToBezierForm(p, lcurve); // Find all possible roots of 5th-degree equation n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0); CCFree((char *)w); // Compare distances of P to all candidates, and to t=0, and t=1 { double new_dist; Point2 q; Vector2 v; // Check distance to beginning of curve, where t = 0 dist = V2SquaredLength(V2Sub(&p, lcurve, &v)); t = 0.0; // Find distances for candidate points for (INT32 i = 0; i < n_solutions; i++) { q = Bezier(lcurve, DEGREE, t_candidate[i], NULL, NULL); new_dist = V2SquaredLength(V2Sub(&p, &q, &v)); if (new_dist < dist) { dist = new_dist; t = t_candidate[i]; } } // Finally, look at distance to end point, where t = 1 new_dist = V2SquaredLength(V2Sub(&p, &(lcurve[DEGREE]), &v)); if (new_dist < dist) { dist = new_dist; t = 1.0; } } *mu = t; return dist; }