decompose B-Spline -> circles


I'm trying to decompose a b-spline to a set of circle definitions. If you're working; yes I'm trying to generate G-Code ;)
Since a number of CNC / CAM apps are based on OCC I'm wondering whether you might be able to point me in the right direction.



Fabian Hachenberg's picture

Extract control points from BSplines? Raster a spline? Calculate derivatives on the spline? In order to point you in the right direction, we need to know exactly what you require to implement your decomposition algorithm using OCC...

jelle's picture

Sure, sorry I was a little unclear. What I'm trying to do is to generate a number of position & vectors along a TopoDS_Wire such that I can write ABB RAPID [ industrial robot ] code from it. The precision of the robot is ~ 0.1 mm. Much like a CNC router, a robot can use linear and circular movements. Since my controller is pretty old, I'd like to be economic in the description of the toolpaths. Hence, moving circular is the preferred option. Wrt the spline decomposition algorithm, honestly, I haven't given it too much thought... I figured that CPnts_UniformDeflection could be useful to get the circle points. For the moment simplicity of implementation outweighs sophistication of the algorithm [ uh, my robots are coming friday ;) ]



Fabian Hachenberg's picture

I would like to point out that your project sounds pretty interesting. For the case that you might end up with something useful to use OCC data as input source for CNC machines, maybe you could add your project to

And concerning the decomposition:
I'm not familiar with the ABB RADID code - could you describe your problem at bit more mathematically? Do you want to interpolate a general curve from OCC (spline, parabola, hyperbola, ...) using only arcs and lines (preferring arcs over lines)? Can you arbitrarily change the radius and central point of the arcs from one segment of the interpolation to the next? (can you approximate a "S"-line using 2 arcs with same radius and central points at (-x,a) and (x,b) ?

Another question: Aren't there any decomposition algorithms publicly avaible yet? So wouldn't it be easier to write a converter from OCC curve description into the respective format of the existing decomposition?

Forum supervisor's picture

Dear all,

OPEN CASCADE has a good experience in solving similar tasks for its customers working in CNC domain. We have good level of expertise and some re-usable solutions we can offer on a paid basis. Please, do not hesitate to contact us via for more detailed information.

Best regards,
Forum supervisor

jelle's picture

Hi Fabian,

>Do you want to interpolate a general curve from OCC (spline, parabola, hyperbola, ...) using only arcs and lines (preferring arcs over lines)?

> Can you arbitrarily change the radius and central point of the arcs from one segment of the interpolation to the next? (can you approximate a "S"-line using 2 arcs with same radius and central points at (-x,a) and (x,b) ?


>Aren't there any decomposition algorithms publicly avaible yet?

yes, there are, I suppose its a question of implementing these. pretty involved from what I've deduced from this [ excellent ] article:


AP's picture

if you have catia, this powercopy/user defined feature shows the principle, and it is well documented.



jelle's picture

cool, thanks for that reference Alex, the GT wiki is a pretty cool place...

AP's picture

Ok Jelle,

I got little bit obsesionado with this problem, I needed it not for NC of robots, but for architecture, the edge of slab of a building needs to be decomposed to arcs to rationalize the construction and simplify the shopdrawings makes it easier on the builder, less prone to errors.

I litterally translated the GTwiki Catia version included in Digital Project from the link i sent you earlier, so this is Alan Hudins method, not mine, im sure someone can point out best ways of getting to same result, but this one worked for me with only 10 biarcs.

I found literature on making more optimal configurations, this version is not the most optimal, im using abusing c++ and geometrical operations (intersections /mirrors) the literature suggests you can find the center points of the arcs by knowing the angles of the start and ending tangent and other minor pithagorean operations, I am not that good at math, so this is the geometrical version.

if you have ProQuest or any other research database, search for this paper:
Yuan-Jye Tseng & Yii-Der Chen (2000): Three dimensional
biarc approximation of freeform surfaces for machining tool path generation,
International Journal of Production Research, 38:4, 739-763



// *********** spaghetti code begins here

/* initialize part with context and folder */
#define INITPART\
int viscount=0;\
TopoDS_Compound folder;\
BRep_Builder B;\

/* create shape and visulize */
#define vis(name)\
if (!name.IsNull()){\

void biarc::update2()

double scale = 200;
QList splineset;
gp_Pnt p1(getdval(x1) * scale,getdval(y1)* scale,getdval(z1)* scale);
gp_Pnt p2(getdval(x2) * scale,getdval(y2) * scale,getdval(z2) * scale);
gp_Pnt(p3,getdval(x3) * scale,getdval(y3)* scale,getdval(z3) * scale);
gp_Pnt(p4,getdval(x4) * scale,getdval(y4) * scale,getdval(z4) * scale);
gp_Pnt(p5,getdval(x5) * scale,getdval(y5) * scale,getdval(z5) * scale);


splineset << p1 << p2 << p3 << p4 << p5;
TopoDS_Shape spline = hsf::AddNewSplineSShape(splineset);
double k = getdval(k)/100;

gp_Pnt sp = hsf::AddNewPointonCurve(spline,0);
gp_Pnt mp = hsf::AddNewPointonCurve(spline,k/2);
gp_Pnt ep = hsf::AddNewPointonCurve(spline,k);

gp_Vec t1 = hsf::getVectorTangentToCurveAtPoint(spline,0);
gp_Vec t2 = hsf::getVectorTangentToCurveAtPoint(spline,k);

double inf=0;
double sup=0;
double count = 30;
for(int i =1;i biarc::biarcjunctiondomain(gp_Pnt p1,gp_Pnt p2,gp_Vec t1,gp_Vec t2, double arcratio)

QList outputs;
gp_Pnt sp = p1;
gp_Pnt ep = p2;

gp_Pnt t1p = hsf::MovePointByVector(sp,t1,1);
//gp_Vec t2(p1,p2);

gp_Pln Spln(sp,t1);
gp_Pln Epln(ep,t2);

gp_Pln pln1 = hsf::AddNewPlane(sp,t1p,ep);
TopoDS_Shape axis = hsf::AddNewIntersectSrf(Spln,Epln);
gp_Vec axisvec = hsf::getVectorTangentToCurveAtPoint(axis,0.5);
gp_Pnt o1 = hsf::AddNewIntersectCrvPln(axis,pln1);

gp_Vec SVec(Spln.Axis().Direction());
gp_Vec EVec(Epln.Axis().Direction());

double angleoriented = (3.14159265 -SVec.AngleWithRef(EVec,axisvec))/2; // 180 deg in radians - angle between planes divided by 2

gp_Pln rotpln1 = Spln.Rotated(gp_Ax1(o1,axisvec),angleoriented); // rotate start plane
//gp_Pln rotpln2 = rotpln1.Rotated(gp_Ax1(o1,axisvec),90); //rotate 90 degs in radians to get perpendicular plane

gp_Pln rotpln2 = rotpln1.Rotated(gp_Ax1(o1,axisvec),1.57079633); //rotate 90 degs in radians

gp_Pnt p2mirror1 = ep.Mirrored(gp_Ax2(rotpln1.Location(),rotpln1.Axis().Direction(),rotpln1.XAxis().Direction()));
gp_Pnt p2mirror2 = ep.Mirrored(gp_Ax2(rotpln2.Location(),rotpln2.Axis().Direction(),rotpln2.XAxis().Direction()));
gp_Pnt farpoint;

double dis1 = p1.Distance(p2mirror1);
double dis2 = p1.Distance(p2mirror2);

if(dis1 > dis2){farpoint = p2mirror1;} else {farpoint = p2mirror2;}

TopoDS_Shape circ1 = hsf::AddNewCircle(ep,sp,farpoint);
Handle(Geom_Curve) circ1crv = hsf::convertshapetocurve(circ1);
TopoDS_Shape arc1 = BRepBuilderAPI_MakeEdge(circ1crv,ep,sp);

TopoDS_Shape point2m1 = hsf::AddNewPoint(p2mirror1);
TopoDS_Shape point2m2 = hsf::AddNewPoint(p2mirror2);
outputs << arc1 << point2m1 << point2m2 < biarc::biarconcurve(TopoDS_Shape crv, double ratio1,double ratio2,double arcratio,double infinite)

gp_Pnt p1 = hsf::AddNewPointonCurve(crv,ratio1);
gp_Pnt p2 = hsf::AddNewPointonCurve(crv,ratio2);

gp_Vec t1 = hsf::getVectorTangentToCurveAtPoint(crv,ratio1);
gp_Vec t2 = hsf::getVectorTangentToCurveAtPoint(crv,ratio2);

QList biarc = biarcptdir(p1,p2,t1,t2,arcratio,infinite);
return biarc;


QList biarc::biarcptdir(gp_Pnt p1,gp_Pnt p2,gp_Vec t1,gp_Vec t2, double arcratio,double infinite)

QList domain = biarcjunctiondomain(p1,p2,t1,t2,arcratio);
gp_Pnt j = hsf::AddNewPointonCurve(,arcratio);

TopoDS_Shape biarc1 = lineptarc(p1,j,t1,infinite);
TopoDS_Shape biarc2 = lineptarc(p2,j,t2,infinite);
QList returndata;
returndata << biarc1 << biarc2 << << <<;
return returndata;


//build 10 biarcs and get closest one to the crv.
TopoDS_Shape biarc::biarcoptimized(TopoDS_Shape crv, double ratio1,double ratio2,double percentinf,double percentsup,double infinite)

double r1 = percentinf;
double r2 = percentinf + (percentsup-percentinf)/9;
double r3 = percentinf + (percentsup-percentinf)* 2/9;
double r4 = percentinf + (percentsup-percentinf)* 3/9;
double r5 = percentinf + (percentsup-percentinf)* 4/9;
double r6 = percentinf + (percentsup-percentinf)* 5/9;
double r7 = percentinf + (percentsup-percentinf)* 6/9;
double r8 = percentinf + (percentsup-percentinf)* 7/9;
double r9 = percentinf + (percentsup-percentinf)* 8/9;
double r10 = percentsup;

QList biarc1 = biarconcurve(crv, ratio1, ratio2,r1, infinite);
QList biarc2 = biarconcurve(crv, ratio1, ratio2,r2, infinite);
QList biarc3 = biarconcurve(crv, ratio1, ratio2,r3, infinite);
QList biarc4 = biarconcurve(crv, ratio1, ratio2,r4, infinite);
QList biarc5 = biarconcurve(crv, ratio1, ratio2,r5, infinite);
QList biarc6 = biarconcurve(crv, ratio1, ratio2,r6, infinite);
QList biarc7 = biarconcurve(crv, ratio1, ratio2,r7, infinite);
QList biarc8 = biarconcurve(crv, ratio1, ratio2,r8, infinite);
QList biarc9 = biarconcurve(crv, ratio1, ratio2,r9, infinite);
QList biarc10 = biarconcurve(crv, ratio1, ratio2,r10, infinite);

int count = biarc1.length();
//vis( vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(
//vis( vis( vis( vis(

double dis1a = hsf::GetMaxDis(,crv);
double dis2a = hsf::GetMaxDis(,crv);
double dis3a = hsf::GetMaxDis(,crv);
double dis4a = hsf::GetMaxDis(,crv);
double dis5a = hsf::GetMaxDis(,crv);
double dis6a = hsf::GetMaxDis(,crv);
double dis7a = hsf::GetMaxDis(,crv);
double dis8a = hsf::GetMaxDis(,crv);
double dis9a = hsf::GetMaxDis(,crv);
double dis10a = hsf::GetMaxDis(,crv);

double dis1b = hsf::GetMaxDis(,crv);
double dis2b = hsf::GetMaxDis(,crv);
double dis3b = hsf::GetMaxDis(,crv);
double dis4b = hsf::GetMaxDis(,crv);
double dis5b = hsf::GetMaxDis(,crv);
double dis6b = hsf::GetMaxDis(,crv);
double dis7b = hsf::GetMaxDis(,crv);
double dis8b = hsf::GetMaxDis(,crv);
double dis9b = hsf::GetMaxDis(,crv);
double dis10b = hsf::GetMaxDis(,crv);

double dis1 = MAX(dis1a,dis1b);
double dis2 = MAX(dis2a,dis2b);
double dis3 = MAX(dis3a,dis3b);
double dis4 = MAX(dis4a,dis4b);
double dis5 = MAX(dis5a,dis5b);
double dis6 = MAX(dis6a,dis6b);
double dis7 = MAX(dis7a,dis7b);
double dis8 = MAX(dis8a,dis8b);
double dis9 = MAX(dis9a,dis9b);
double dis10 = MAX(dis10a,dis10b);

QMap> map;

QList distances = map.keys();

double minval =;
TopoDS_Shape arc1 = map.value(minval).at(0);
TopoDS_Shape arc2 = map.value(minval).at(1);

vis(arc1) vis(arc2)

return folder;


// utility functions from my Hybrid Shape Factory Library
// I am used to scripting Catia in Visual Basic
// so I use an Idiom that follows the Catia HybridShape Factory concept

const gp_Pnt& HSF::AddNewPointonCurve(TopoDS_Shape SupportEdge, Standard_Real uRatio)

const TopoDS_Edge& aEdge = TopoDS::Edge (SupportEdge);
Standard_Real aFP, aLP, aP;
Handle(Geom_Curve) aCurve = BRep_Tool::Curve(aEdge, aFP, aLP);
aP = aFP + (aLP - aFP) * uRatio;
gp_Pnt p1 = aCurve->Value(aP);

return p1;

const gp_Vec& HSF::getVectorTangentToCurveAtPoint(TopoDS_Shape SupportEdge, Standard_Real uRatio)

const TopoDS_Edge& aEdge = TopoDS::Edge (SupportEdge);
Standard_Real aFP, aLP, aP;
Handle(Geom_Curve) aCurve = BRep_Tool::Curve(aEdge, aFP, aLP);
aP = aFP + (aLP - aFP) * uRatio;
gp_Vec& V1 = gp_Vec() ;
gp_Pnt p1;

return V1;


Handle_Geom_Curve HSF::convertshapetocurve(TopoDS_Shape theedge)
double fp,lp;
Handle_Geom_Curve Result = BRep_Tool::Curve(TopoDS::Edge(theedge),fp,lp);
//TopoDS_Shape Result = revolved.Solid();
return Result;


TopoDS_Edge HSF::AddNewCircle(gp_Pnt p1,gp_Pnt p2,gp_Pnt p3)

Handle_Geom_Circle mycircle = GC_MakeCircle(p1,p2,p3);

TopoDS_Edge Result;
if (mycircle.IsNull()) {return Result;}
Result= BRepBuilderAPI_MakeEdge(mycircle);

return Result;
return Result;


TopoDS_Shape HSF::AddNewIntersectSrf(gp_Pln pln1,gp_Pln pln2)

Handle_Geom_Plane S1 = GC_MakePlane (pln1);
Handle_Geom_Plane S2 = GC_MakePlane (pln2);

//This class is instantiated as follows:
GeomAPI_IntSS Intersector(S1, S2, Precision::Intersection());

//Once the GeomAPI_IntSS object has been created, it can be interpreted.
//Calling the number of intersection curves
Standard_Integer nb = Intersector. NbLines();

//Calling the intersection curves
if (nb > 0) {
Handle(Geom_Curve) C = Intersector.Line(1);

//where Index is an integer between 1 and Nb.
Handle(Geom_TrimmedCurve) spinaxis = new Geom_TrimmedCurve (C, -200, 200);
TopoDS_Shape Result = BRepBuilderAPI_MakeEdge(spinaxis);

return Result;
TopoDS_Shape nullshape;
return nullshape;

gp_Pnt HSF::AddNewIntersectCrvPln(TopoDS_Shape crv1,gp_Pln pln1)
gp_Pnt intp;
Standard_Real af, al;
Handle(Geom_Curve) crv = BRep_Tool::Curve(TopoDS::Edge(crv1),af,al);

Handle_Geom_Plane aSurf = new Geom_Plane(pln1);

//This class is instantiated as follows:
GeomAPI_IntCS Intersector(crv, aSurf);

//Once the GeomAPI_IntSS object has been created, it can be interpreted.
//Calling the number of intersection curves
Standard_Integer nb = Intersector.NbPoints();

//Calling the intersection curves
if (nb > 0) {
intp = Intersector.Point(1);

return intp;


gp_Pnt HSF::MovePointByVector(gp_Pnt point, gp_Vec vector, double Offset)

gp_Vec aV = vector.Normalized();
aV = aV.Multiplied(Offset);

gp_Pnt newpoint;
newpoint.SetXYZ(point.XYZ() + aV.XYZ());

return newpoint;

TopoDS_Shape HSF::AddNewPoint(gp_Pnt point)
BRepBuilderAPI_MakeVertex mkVertex (point);
return mkVertex.Shape();


AP's picture

// Missed the Spline function :

TopoDS_Shape HSF::AddNewSplineSShape(QList Points)

TColgp_Array1OfPnt arrayval (1,Points.count()); // sizing array

QListIterator iT(Points);
int count = 0;
while (iT.hasNext())
gp_Pnt curpoint =;
count +=1;

Handle(Geom_BSplineCurve) SPL1 = GeomAPI_PointsToBSpline(arrayval).Curve();

/* GeomConvert_ApproxCurve aApprox(SPL1,Precision::Confusion(),GeomAbs_C1,3,3);
SPL1 = aApprox.Curve();*/

double fp,lp;
fp = SPL1->FirstParameter();
lp = SPL1->LastParameter();
TopoDS_Shape aShape = BRepBuilderAPI_MakeEdge(SPL1,fp,lp);
return aShape;


AP's picture

be wary of those "&gt&gt" they are suppossed to be >>
remove all the "&nbsp";

They must have been tabs and spaces in Visual Studio, who knows...



Nenad's picture

Hi Alex,
does your code work in Catia v5?

AP's picture

Solution Space before trimming

AP's picture

solution space after trimming

jelle's picture

some pretty gorgous screenshots Alex!
interesting; I'm an architect as well!
what about we catch up off list to talk OCC & architecture?
would be interesting!

jelleferinga youknowthecharacter gmail dot com

jelle's picture

btw [your?] HSF lib looks pretty neat!
is there a repo somewhere? I'd love to know more about it...

AP's picture

I guess its the forum:

remove the ampersand-n-b-s-p-semicolon from the the first post, if you see them.
and two ampersand-g-t-semicolon are suppossed to be ">>" if you see them.

jelle's picture

thanks for sharing that code Alex.
will add a pythonocc version of it to hb-robotics-code soon!

AP's picture

No Worries Jelle ...



Arthur Magill's picture

Hey Jelle,

In case you're looking for some code to work from, these guys have solved the same problem to convert SVG (from Inkscape) into G-code. I haven't used the plugin, but I think it does what you want. And you'll be pleased to see it's written in Python ;-)

jelle's picture

really cool! many thanks for the heads up Arthur!

JuryS's picture

Hello, Jelle!
It's cool think to decode approximation arcs, circles, B-spline curves with gl2ps module when we save the svg, ps or another vector format. Because when we search objects in frame buffer we have no token like GL_ARC or GL_CIRCLE.
And we can calculate the lines into Arc or Circle.
Now I'm test the svg, saved with Approximation lines and with Arcs. And the first svg file bigger for 10 times. Also more fast rendering for any vector formats and 1000 more reasons to use transformation of approximations to lines

some part of gl.h:
/* FeedBackToken */
#define GL_PASS_THROUGH_TOKEN 0x0700
#define GL_POINT_TOKEN 0x0701
#define GL_LINE_TOKEN 0x0702
#define GL_POLYGON_TOKEN 0x0703
#define GL_BITMAP_TOKEN 0x0704
#define GL_DRAW_PIXEL_TOKEN 0x0705
#define GL_COPY_PIXEL_TOKEN 0x0706
#define GL_LINE_RESET_TOKEN 0x0707

jelle's picture

hi Yuriy,

right, I see, so without knowing whether you're dealing with GL_ARC or GL_CIRCLE its hard to optimize / encode things efficiently?

JuryS's picture

Sorry for my bad English. I'm talk about optimize output of vector graphic files, no any OpenGL modification. I'm use the OCC HLR algorithm and all lines that was included in this algo after have approximations with StdPrs_DeflectionCurve. I'm know which type of object I'm input in HLR algo, but I don't know what I have after. Also I can save my scene without gl2ps module and OpenGL framebuffer, in my application all objects have drivers, that include information about geom and type of object.
And I don't know what I must to do: modify the HLR to save information about lines that I have input and after store this lines under driver of object, or modify gl2ps module for calculating the lines and found there some circles and arcs.

First things that I modify the HLR, because after use of this algorithm I don't have information about objects, I at all don't know witch color of lines was inserted in HLR algo.

alton's picture

I really like your way of expressing the opinion and sharing the information.

Zhongxing Xu's picture

I found a related paper that I have implemented its algorithm:

Les A. Piegl, Wayne Tiller, Biarc approximation of NURBS curves, Computer-Aided Design 34, 2002, 807-814