decompose B-Spline -> circles

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
http://opencascade.wikidot.com/projects
?

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?

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

Best,

AlexP

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

/* initialize part with context and folder */
#define INITPART\
UnitsAPI::SetLocalSystem(UnitsAPI_MDTV);\
int viscount=0;\
TopoDS_Compound folder;\
BRep_Builder B;\
B.MakeCompound(folder);\

/* create shape and visulize */
#define vis(name)\
if (!name.IsNull()){\
B.Add(folder,name);\
viscount++;}\

void biarc::update2()
{
//INITPART

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);

//visgeo(point1,SphereSurface(p1,1))
//visgeo(point2,SphereSurface(p2,1))
//visgeo(point3,SphereSurface(p3,1))
//visgeo(point4,SphereSurface(p4,1))
//visgeo(point5,SphereSurface(p5,1))

splineset << p1 << p2 << p3 << p4 << p5;
TopoDS_Shape spline = hsf::AddNewSplineSShape(splineset);
//vis(spline)
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(domain.at(0),arcratio);

TopoDS_Shape biarc1 = lineptarc(p1,j,t1,infinite);
TopoDS_Shape biarc2 = lineptarc(p2,j,t2,infinite);
QList returndata;
returndata << biarc1 << biarc2 << domain.at(0) << domain.at(1) << domain.at(2);
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)
{

INITPART
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(biarc1.at(0)) vis(biarc1.at(1)) vis(biarc1.at(2)) vis(biarc1.at(3)) vis(biarc1.at(4))
//vis(biarc2.at(0)) vis(biarc2.at(1)) vis(biarc2.at(2)) vis(biarc2.at(3))
//vis(biarc3.at(0)) vis(biarc3.at(1)) vis(biarc3.at(2)) vis(biarc3.at(3))
//vis(biarc4.at(0)) vis(biarc4.at(1)) vis(biarc4.at(2)) vis(biarc4.at(3))
//vis(biarc5.at(0)) vis(biarc5.at(1)) vis(biarc5.at(2)) vis(biarc5.at(3))
//vis(biarc6.at(0)) vis(biarc6.at(1)) vis(biarc6.at(2)) vis(biarc6.at(3))
//vis(biarc7.at(0)) vis(biarc7.at(1)) vis(biarc7.at(2)) vis(biarc7.at(3))
//vis(biarc8.at(0)) vis(biarc8.at(1)) vis(biarc8.at(2)) vis(biarc8.at(3))
//vis(biarc9.at(0)) vis(biarc9.at(1)) vis(biarc9.at(2)) vis(biarc9.at(3))
//vis(biarc10.at(0)) vis(biarc10.at(1)) vis(biarc10.at(2)) vis(biarc10.at(3))

double dis1a = hsf::GetMaxDis(biarc1.at(0),crv);
double dis2a = hsf::GetMaxDis(biarc2.at(0),crv);
double dis3a = hsf::GetMaxDis(biarc3.at(0),crv);
double dis4a = hsf::GetMaxDis(biarc4.at(0),crv);
double dis5a = hsf::GetMaxDis(biarc5.at(0),crv);
double dis6a = hsf::GetMaxDis(biarc6.at(0),crv);
double dis7a = hsf::GetMaxDis(biarc7.at(0),crv);
double dis8a = hsf::GetMaxDis(biarc8.at(0),crv);
double dis9a = hsf::GetMaxDis(biarc9.at(0),crv);
double dis10a = hsf::GetMaxDis(biarc10.at(0),crv);

double dis1b = hsf::GetMaxDis(biarc1.at(1),crv);
double dis2b = hsf::GetMaxDis(biarc2.at(1),crv);
double dis3b = hsf::GetMaxDis(biarc3.at(1),crv);
double dis4b = hsf::GetMaxDis(biarc4.at(1),crv);
double dis5b = hsf::GetMaxDis(biarc5.at(1),crv);
double dis6b = hsf::GetMaxDis(biarc6.at(1),crv);
double dis7b = hsf::GetMaxDis(biarc7.at(1),crv);
double dis8b = hsf::GetMaxDis(biarc8.at(1),crv);
double dis9b = hsf::GetMaxDis(biarc9.at(1),crv);
double dis10b = hsf::GetMaxDis(biarc10.at(1),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;
map.insert(dis1,biarc1);
map.insert(dis2,biarc2);
map.insert(dis3,biarc3);
map.insert(dis4,biarc4);
map.insert(dis5,biarc5);
map.insert(dis6,biarc6);
map.insert(dis7,biarc7);
map.insert(dis8,biarc8);
map.insert(dis9,biarc9);
map.insert(dis10,biarc10);

QList distances = map.keys();
qSort(distances);

double minval = distances.at(0);
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;
aCurve->D1(aP,p1,V1);

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;}
try
{
Result= BRepBuilderAPI_MakeEdge(mycircle);
}
catch(...)
{

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;
}
else
{
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();

}

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.

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