Biarc Curve Fitting 

Hakan Tiftikci
Turkey
hakan.tiftikci@yahoo.com.tr 

Introduction 

Biarc curve fitting, as the name implies, fits a curve consisting of two circular arcs through two points with specified tangents. For a given set points, biarc curve fitting may be applied for successive pairs of points with common given /computed tangent for each point, to generalize the curve fitting process. This worksheet  demonstrates biarc curve fitting, following the work of Jaroslaw and Aristides[1]. Another source for biarc curves is Choi [2] which contains only 2D specific formulation with less freedom and thus limited compared to formulation of [1]. 

 

Initialization 

In this part, auxiliary functions
 

 

• dot

 

 

 

• cross

 

 

 

• norm2

 

 

 

• threepoint2circlearc

 

 

 

• rationalArc
  
  are defined.

 

 

>

restart;

 

 

dot 

>	dot := (p,q) -> LinearAlgebra:-DotProduct(p,q,conjugate=false):

 

cross 

>	cross := (p,q) -> LinearAlgebra:-CrossProduct(p,q):

 

norm2  

>	norm2 := v -> sqrt(dot(v,v)):

 

 

Determination of circular arc parameters from arc end  points and apex point 

Biarc curve fitting, as described here, determines junction point and here-so-called apex points. Circular arc parameters (radius, angle, center, ...) from two end points and meeting point of tangents at are determined via following formulas 

 

 

 

 

 

where
 

 

threepoint2circlearc := proc(A,B,C)
description "determines the center, radius,angle extent and basis vectors for a arc end points 'A,B' and apex point 'C'";
local AB,BC,lAB,lBC,ab,bc,E1,E2,E3,cos2t,sin2t,cost,sint,tant,theta,r,bisector,center;
AB := B-A;
BC := C-B;
lAB := norm2(AB);
lBC := norm2(BC);
ab := AB/lAB;
bc := BC/lBC;

E3 := cross(ab,bc);
sin2t := norm2(E3);
E3 := E3/sin2t;

cos2t := dot(ab,bc);
cost := sqrt((1+cos2t)/2);
sint := sqrt((1-cos2t)/2);
tant := sint/cost;
theta := arctan(sin2t,cos2t)/2;
r := lAB/tant;  # tant (instead of tan(theta)) is intentially used since it is positive

bisector := bc-ab;
center := B + (r/cost)*bisector/norm2(bisector);

E1 := A-center;
E1 := E1/norm2(E1);

E2 := cross(E3,E1);

center,r,theta,E1,E2,E3;
end proc:  

 

 

Rational polynomial representation of circular arc 

Although not necessary, rational representation of circular arc is employed here. A circular arc centered at origin and having angle sweep is defiend as 

 

 

 

 

 

 

 

Following procedure generalizes this formula to general orientation by using basis vectors as input arguments 

 

 

>

rationalArc := proc(A,B,E3,theta,u) description "defines rational arc with origin 'A' base 'AB' normal 'E3', angle extent 'theta' and parameter 'u'"; local P0,P1,P2,b0,b1,b2,w0,w1,w2,E1,E2,N,D; E1 := B-A; E2 := cross(E3,E1); P0 := B; P1 := P0+tan(theta)*E2; P2 := A+cos(2*theta)*E1+sin(2*theta)*E2; w0 := 1; w2 := 1; w1 := cos(theta); b0 := (1-u)^2; b1 := 2*(1-u)*u; b2 := u^2; N := b0*w0*P0+b1*w1*P1+b2*w2*P2; D := b0*w0+b1*w1+b2*w2; N/D; end proc:

 

 

Test of routine "rationalArc" 

>

arc0 := rationalArc(<0,0,0>,<R,0,0>,<0,0,1>,theta,u): simplify(norm2(arc0)) assuming R::real; arc1 := subs(theta=+Pi/4/2,R=1,arc0): arc2 := subs(theta=-Pi/4/2,R=1,arc0): plots[display]({ plot([arc1[1],arc1[2],u=0..1],scaling=constrained,axes=boxed,color=red), plot([arc2[1],arc2[2],u=0..1],scaling=constrained,axes=boxed,color=blue)},view=[-1..1,-1..1],scaling=constrained);

 

 

 

 

 

Formulation (Jaroslaw and Aristides) 

The figure below illustrates a a biarc curve passing through points with  tangents together with the dimensions used to formulate the problem..  

 

 

In this document, the intersection point of tangents for one of the arcs is called "apex". Apex points in the figure are given by 

 

 

 

 

 

Junction point of two arcs is 

 

Those points are subject to following length constraints 

 

 

 

 

 

 

 

 

 

which yield, after some manipulation, the general equation as 

 

 

 

If the ratio of lengths is specified as 

 

 

 

then quadratic equation 

 

 

is obtained. These formulas are implemented in "biarcBasic" and "biarc" procedures. Note that "biarc" procedure returns only positive solution if "sign" parameter is 0, otherwise sign should be (or may be left as a symbol to denote ) and the result will be dependent on the value of sign (for a generic solution) 

 

>	biarcBasic := proc(P1,T1,P2,T2,a1,a2) local S; S := P2-P1; print(dot(T1,T2)); print(dot(S,T2)); print(dot(S,T1)); print(dot(S,S)); a1*a2*(dot(T1,T2)-1) - a2*dot(S,T2) - a1*dot(S,T1) + dot(S,S)/2; end proc:

 

>

biarc := proc(P1,T1,P2,T2,rho,sign:=0) description "determines two G1-continuous arcs connecting P1 to P2 and having tangents T1,T2 at P1,P2."; local S,x,c2,c1,c0,a1s,a1,a2,B1,B2,B12,center1,center2,r1,r2,E1,E2,E3,F1,F2,F3,theta1,theta2,arc1,arc2; S := P2-P1; c2 := rho*(dot(T1,T2)-1); c1 := -dot(S,T1+rho*T2); c0 := dot(S,S)/2; if evalb(not type(sign,symbol) and  (sign=0)) then a1s := solve(c2*x^2+c1*x+c0,x); a1 := select( x -> verify( x, 0, greater_than), [a1s]); a1 := a1[1]; else a1 := (-c1+sign*sqrt(c1^2-4*c2*c0))/(2*c2); end if; a2 := rho*a1; B1 := P1+a1*T1; B2 := P2-a2*T2; B12 := (a2*B1+a1*B2)/(a1+a2); center1,r1,theta1,E1,E2,E3 := threepoint2circlearc(P1,B1,B12); center2,r2,theta2,F1,F2,F3 := threepoint2circlearc(B12,B2,P2); arc1 := rationalArc(center1,P1,E3,theta1,u); arc2 := rationalArc(center2,B12,F3,theta2,u); (arc1,arc2,[center1,P1,B1,B12],[center2,B12,B2,P2]); end proc:

 

>

plotArc := proc(g,arcpoints,u,color1,color2,trans,thick1,thick2,in2D:=false) description "plots arc 'g' using color 'color1' and thickness 'thick1' and plots arc points as transparent polygon using color 'color2' and thickness 'thick2'"; local arcplot,polyplot; if in2D then arcplot := plot([g[1],g[2],u=0..1],color=color1,thickness=thick1,adaptive=8); polyplot := plottools[polygon](map(x->convert(x[1..2],'list'),arcpoints),color=color2,transparency=trans,thickness=thick2,linestyle=DASHDOT); else arcplot := plots[spacecurve]([g[1],g[2],g[3]],u=0..1,color=color1,thickness=thick1); polyplot := plottools[polygon](map(x->convert(x,'list'),arcpoints),color=color2,transparency=trans,thickness=thick2,linestyle=DASHDOT); end if; arcplot,polyplot; end proc:

 

Formulation (Choi) 

Although not examplified, formulation given in Choi[2] is summarized here for completeness. Formulation given in Choi[2] has less freedom (as opposed to 2-parameter family of solutions of [1]) due to some geometric conditions imposed  on the biarc curve. It is based on radii , instead of lengths . Also it is discriminated for unimodal and  inflection curves. For unimodal curves with , radii are given by 

 

 

 

 

 

where 

 

 

 

 

 

 

 

 

 

For inflection biarc with equal radii, formula is 

 

 

 

 

 

Examples 

This section contains three examples of biarc curve fitting 

Example-1 : Single 3D biarc curve segment 

>	theta1 := Pi/3: theta2 := -Pi/4: P1 := <0,0,0>: T1 := <cos(theta1),sin(theta1),0>: P2 := <1,0,1>: T2 := <cos(theta2),sin(theta2),0>:

 

>	arc1,arc2,arcnet1,arcnet2 := biarc(P1,T1,P2,T2,1):

 

>	plotBiarcSegment := proc(arc1,arc2,theta,phi) plots[display]({ plotArc(arc1,arcnet1,u,red,orange,0.5,3,2), plotArc(arc2,arcnet2,u,blue,cyan,0.5,3,2)}, axes=boxed,orientation=[theta,phi]); end proc:

 

>	plotBiarcSegment(arc1,arc2,60,110); plotBiarcSegment(arc1,arc2,20,-60);

 

 

 

 

 

 

Example-2: Multiple biarc curve segments 

In this example biarc curve fitting is performed on more than 2 points to yield more than 2 biarc curve segments. Computation of tangents is similar to that of Catmull-Rom splines, i.e.
 

 

 

same ratio is assumed for all segments. Not to exhaust the memory too much, results are returned in floating-point evaluated form. 

>

BiarcCurveFit := proc(P,rho) description "determines biarc curve fit to a sequence of points 'P'. Uses same length ratio 'rho' for  all biarc segments"; local i,n,T,A1,AN1,A2,AN2,ALLARCS,ALLARCNETS; n := nops(P); T := [0$n]; ALLARCS := []; ALLARCNETS := []; for i from 1 to n do if (i=1) then T[i] := (P[2]-P[1]); elif (i<n) then T[i] := (P[i+1]-P[i-1]); else T[i] := (P[n]-P[n-1]); end if; T[i] := evalf(T[i]/norm2(T[i])); if (i>1) then A1,A2,AN1,AN2 := biarc(P[i-1],T[i-1],P[i],T[i],rho); ALLARCS := [op(ALLARCS), evalf(A1),evalf(A2)]; ALLARCNETS := [op(ALLARCNETS), evalf(AN1),evalf(AN2)]; end if; end do; ALLARCS, ALLARCNETS; end proc:

 

Create sample points 

>	n := 11: P := [0$n]: for i from 1 to n do Rho := (i-1)/(n-1); P[i] := <(2+cos(Rho*2*Pi))*cos(3*2*Pi*Rho),(2+cos(Rho*2*Pi))*sin(3*2*Pi*Rho),0>; end do:

 

>	plots[display]( {plots[listplot]( map( x->convert(x[1..2],'list'), P )), plots[pointplot]( map( x->convert(x[1..2],'list'), P ), symbol=solidcircle,symbolsize=18 )}, axes=boxed );

 

 

Determine biarc curve fit and  display results 

>	arcs, arcpoints := BiarcCurveFit(P, 1):

 

>	twocolors :=[red,blue]; plots[display]({ seq( plotArc(arcs[i],arcpoints[i],u,twocolors[1+ modp(i-1,2)] ,blue,0.8,2,1,true), i=1..nops(arcs) ), plots[pointplot]( map( x->convert(x[1..2],'list'), P ), symbol=solidcircle,symbolsize=18, axes=boxed)}, axes=boxed,scaling=constrained,view=[-3..3,-3..3]);

 

 

 

 

Example-3 : Converting a parametric curve a sequence of circular arcs 

In this example procedure "curve2arcs" samples parametric curve at given parameter values for position and tangent values and determines the sequence of biarc curve fits for the computed curve points and tangents. The procedure returns both rational arc representations and arc points as lists. 

 

>

curve2arcs := proc(g,t,tvals,rho) description "determines biarc curve fit to a parametric curve 'g(t)' sampled at parameter values 'tvals'. Uses same biarc curve parameter 'rho=a2/a1' for all segments"; local i,n,P1,T1,P2,T2,A1,A2,AN1,AN2,dg,ALLARCS,ALLARCNETS; n := nops(tvals); dg := map(diff,g,t); ALLARCS := []; ALLARCNETS := []; for i from 1 to n-1 do P1 := eval(g,t=tvals[i]); T1 := eval(dg,t=tvals[i]); P2 := eval(g,t=tvals[i+1]); T2 := eval(dg,t=tvals[i+1]); T1 := T1/norm2(T1); T2 := T2/norm2(T2); A1,A2,AN1,AN2 := biarc(P1,T1,P2,T2,rho); ALLARCS := [op(ALLARCS), A1,A2]; ALLARCNETS := [op(ALLARCNETS), AN1,AN2]; end do; ALLARCS,ALLARCNETS; end proc:

 

This procedure is applied to two curves 

 

 

>	g1 := <cos(5*t),2*sin(3*t),sin(t)>: plots[spacecurve]([g1[1],g1[2],g1[3]],t=0..2.4,numpoints=132,axes=boxed,scaling=constrained,thickness=2,title=gamma[1]);

 

 

>	g2 := <t*cos(5*t),t*sin(5*t),sin(3*t)>: plots[spacecurve]([g2[1],g2[2],g2[3]],t=0..2.4,numpoints=132,axes=boxed,scaling=constrained,thickness=2,title=gamma[2]);

 

 

 

>	arcsg1,arcpointsg1 := curve2arcs(g1,t,[0,0.4,0.8,1.2,1.6,2.0,2.4],1): arcsg2,arcpointsg2 := curve2arcs(g2,t,[0,0.4,0.8,1.2,1.6,2.0,2.4],1):

 

Resulting biarc curves are rendered below.  Note that the result for the first curve shown below contains a segment which does not fit to original curve, which is possibly due to non-adaptive (uniform) sampling of curve and/or non-adaptive selection of length ratio . Biarc curve fit to second curve fits quite well 

>

plots[display]({ seq( plotArc(arcsg1[i],arcpointsg1[i],u,red,blue,0.8,3,1)[1], i=1..nops(arcsg1) ), plots[spacecurve]([g1[1],g1[2],g1[3]],t=0..2.4,thickness=3,color=green)}, axes=boxed,scaling=constrained); plots[display]({ seq( plotArc(arcsg1[i],arcpointsg1[i],u,red,blue,0.8,3,1), i=1..nops(arcsg1) )}, view=[-1..1,-2..2,0..1],axes=boxed,scaling=constrained);

 

 

 

 

>

plots[display]({ seq( plotArc(arcsg2[i],arcpointsg2[i],u,red,blue,0.8,3,1)[1], i=1..nops(arcsg1) ), plots[spacecurve]([g2[1],g2[2],g2[3]],t=0..2.4,thickness=3,color=green)}, axes=boxed,scaling=constrained); plots[display]({ seq( plotArc(arcsg2[i],arcpointsg2[i],u,red,blue,0.8,3,1), i=1..nops(arcsg1) )}, axes=boxed,scaling=constrained);

 

 

 

 

 

Conclusions 

 

• Biarc curve fitting technique is implemented and illustrated on three examples.

 

 

 

• In [1], biarc curve fitting is applied to model b-rep surfaces (to represent bounadries in parameter domain), but in practice they may be applied to any problem that requires curve fitting. A peculiar application is the generation of CNC circular G-codes which are better than line G-codes.

 

 

References 

[1] Jaroslaw R.Rossignac, Aristides A.G. Requicha, "Piecewise circular curves for geometric modeling", , IBM J. RES. DEVELOP.VOL.3, NO.3 May 1987 

[2] B.K. Choi, "Surface Modeling for CAD/CAM". Elsevier, Amsterdam, 1991 

 

All worksheets must contain the following disclaimer: 

Legal Notice: ? Maplesoft, a division of Waterloo Maple Inc. 2009. Maplesoft and Maple are trademarks of Waterloo Maple Inc. Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting from the use of this material.  This application is intended for non-commercial, non-profit use only. Contact the authors for permission if you wish to use this application in for-profit activities.