BRepAlgoAPI_Section and Shape-Shape Intersection Optimization

Hello, I was wondering if i could get some help/feedback/suggestions in a better understanding of the Mechanism by which BRepAlgoAPI_Section works, the philosophy, readings, who invented this technique?(Ivan Sutherland?, University of Utah?) is it documented on any siggraph papers? or any other proceedings.

the code runs fast, but i want to make it even faster, because I need to execute thousands of intersections and when summed up the accumulative time it takes will be in hours, would love to drop that total time to minutes.

so far i have been executing intersections between an imported IGS shell that may have more than one face, and another shell that may have more than one face (sphere).

I like Catia\'s vb CAA api, so if you have used the HybridShapeFactory, I built one for my self using Opencascade, Hence the HSF:: acronim.

here is a static function that returns the intersection:

TopoDS_Shape HSF::AddNewIntersectSrfW(TopoDS_Shape srf1,TopoDS_Shape srf2)
{
BRepAlgoAPI_Section asect(srf1,srf2,Standard_False);
asect.Approximation(Standard_True);
asect.Build();
TopoDS_Shape R = asect.Shape();
return R;
}

so as i said the code above runs fast, i just need to make it go faster.

So i looked into the way BRepAlgoAPI_Section works:

Stage 1(construction):
BRepAlgoAPI_Section subclasses BRepAlgoAPI_BooleanOperation
in the construction stage we feed the inputs to the base class BRepAlgoAPI_BooleanOperation(Sh1, Sh2, BOP_SECTION) letting the boolean operation know it is a BOP_SECTION.

If PerformNow is set to true it will build the curves at construction.

Stage2:
the base class BRepAlgoAPI_BooleanOperation calls preparefiller()
this creates a BOPTools_DSFiller if there is none
and it calls BOPTools_DSFiller::SetShapes(sh1,sh2)
this besides setting shape1 and shape 2 as references, checks their types and returns an error if the types are not compatible.

stage 3(build)

a call to BRepAlgoAPI_Section::Build()

after confirming there where no errors in the creation of the filler
a call to BOPTools_DSFiller::Perform passing the section attributes (computepcurveon1...aprox..)

in the perform function a BooleanOperations_ShapesDataStructure is created
then a reference is passed to BOPTools_InterferencePool.

then a BOPTools_PaveFiller object is created and the BOPTools_InterferencePool is passed with section attributes.

then the BOPTools_PaveFiller::perform is called

6 checks are performed:
vertex-vertex
vertex-edge
vertex-face
edge-edge
edge-face
face-face

and for each check, the respective geometry is created inside interference tables of the pavefiller.

back in the BRepAlgoAPI_Section::Build()

a BOP_Section object is created
BOP_SectionHistoryCollector object is created and passed the 2 shapes
then this is set as the history collector of the BOP_Section

a call to BOP_Section::dowithfiller passing a reference to the BOPTools_DSFiller made earlier

finally the shape can be extracted from the BOP_Section
by calling BOP_Section::Result()

**************************************************************************************

I have lots of questions about the whole process:

it seems the pseudo code is:

- check the type compatibility of the shapes to intersect.
- build a table of all the subshapes of shape 1 and shape 2
- check the shape ancestors and successors (dont know what this means)
- check the interference of the subshapes in shape1 with the subshapes in shape 2, based on ancestry.
(what is ment by interference? i noticed alot of mapping between a what and a with...)
- using this data gathered, calculate the apropiate intersection algorithm, vv,ve,vf,ee,ef,ff
- return results as a single shape.

so i thought if you already know before hand what shape 1 and shape 2 is you could jump to the corresponding ff or vv .. code

so i wrote this function, i had to make it work in the case shape1 has more than one face or shape2 has more than one face.

TopoDS_Shape HSF::AddNewIntersectFast(TopoDS_Shape srf1,TopoDS_Shape srf2)
{

TopoDS_Shape result;
BRepBuilderAPI_MakeWire thewire = BRepBuilderAPI_MakeWire();

TopExp_Explorer shapeX;
for (shapeX.Init(srf1,TopAbs_FACE); shapeX.More(); shapeX.Next())
{
TopoDS_Face aF1 = TopoDS::Face(shapeX.Current());

TopExp_Explorer shapeY;
for (shapeY.Init(srf2,TopAbs_FACE); shapeY.More(); shapeY.Next())
{
TopoDS_Face aF2 = TopoDS::Face(shapeY.Current());

double anApproxTol=1.e-7;

IntTools_FaceFace aFF;
aFF.SetParameters (true,
true,
true,
anApproxTol);

aFF.Perform(aF1, aF2);

if (aFF.IsDone()) {}
aFF.PrepareLines3D();
//
Standard_Integer aNbCurves, aNbPoints;
const IntTools_SequenceOfCurves& aCvs=aFF.Lines();
aNbCurves=aCvs.Length();

for (int i=1; i const IntTools_Curve& aIC=aCvs(i);
Handle(Geom_Curve) mycrv = aIC.Curve();
if (aIC.HasBounds())
{
double fp,lp;
gp_Pnt fpoint, lpoint;
aIC.Bounds(fp,lp,fpoint,lpoint);
TopoDS_Edge myedge = BRepBuilderAPI_MakeEdge(mycrv);
thewire.Add(myedge);
}
}
}
}
}

result = thewire.Shape();
return result ;
}

this code is running way slower than the BRepAlgoAPI_Section code, and it is not robust, if i intersect a sphere with a surface, in some cases it wont return a full loop distorted circle, it returns partial sections of the circle, the BRepAlgoAPI_Section code always returns the correct full circle.

Any Suggestions are welcome...

I uploaded a sample application that should run SFMQTGUI.EXE if you have occ installed, and if you dont have occ installed, call occinstall.bat to make the enviroment variables. the interface is Derived from the QTOCC by Peter Dolby.

https://rcpt.yousendit.com/1134128785/3b0a3255ab628ebfd3b2fa44b4525593?cid=tr-cv-web-stnd-dl-literec-null-c-19073&s=19073

the Hcount slider sets how many spheres in the horizontal row, and the vcount sets how many spheres in the vertical direction.

Thanks for Listening.

Best,

Alex

the code that generates this example :

void mus::update()
{

QTime time;
time.start();

INITPART

double hcount = get(hcount)
double vcount = get(vcount)

try
{

QoccInputOutput* io_man = new QoccInputOutput();
QString srfpath = QoccApplication::instance()->applicationDirPath() + QString(\"/MPP_SOU_DRV_ENV_SRF1.igs\");

Handle(TopTools_HSequenceOfShape) famshape1 = io_man->importIGES(srfpath);
TopoDS_Shape SRF = famshape1->Value(1);

vis(SRF)

double Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
Bnd_Box Bnd;
BRepBndLib::Add(SRF, Bnd);
Bnd.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);

double width = abs(Xmax -Xmin);
double height = abs(Ymax -Ymin);

gp_Pnt centerp((Xmin+Xmax)/2,(Ymin+Ymax)/2,(Zmin+Zmax)/2);
gp_Vec up = gp::DZ();
gp_Pnt cog = HSF::GetCOG(SRF);

TopoDS_Shape rec = HSF::AddNewRectangle(cog,up,width,height);
TopoDS_Shape recfill = HSF::AddNewFillSurface(rec);
//vis (rec)
//vis(recfill)

TopoDS_Shape startcrv = HSF::AddNewIntersectSrfW(SRF,recfill);
vis(startcrv)

TopoDS_Shape startedge = HSF::getedgefromshape(startcrv);

double crvlen = HSF::GetLength(startedge);
crvlen = crvlen;
double minrad = floor(crvlen / hcount);
shapelist(spherelist1)

FOR(i,1,hcount)
try
{
gp_Pnt p1 = HSF::AddNewPointonCurve(startedge,(double) (i-1) / (hcount));
TopoDS_Shape sphere = HSF::AddNewSphereSurface(p1,minrad);
//vis(sphere);

TopoDS_Shape intsphere = HSF::AddNewIntersectSrfW(SRF,sphere);
//TopoDS_Shape intsphere = HSF::AddNewIntersectFast(SRF,sphere);

vis(intsphere);
addtolist(spherelist1,intsphere)
QString line;
QTextStream linestream(&line);
linestream ui::getInstance()->Statusbar(line);
}
catch(...)
{

}
END

shapelist(prevspheres)
prevspheres = spherelist1;
FOR(j,0,vcount)
shapelist(nextspheres)
FOR(i,0,prevspheres.size()-1)
try {
TopoDS_Shape sphere1 = prevspheres.at(i);
TopoDS_Shape sphere2 = prevspheres.at(i+1);

TopoDS_Shape p2 = HSF::AddNewIntersectionPointMulti(sphere1,sphere2);
vis(p2);
TopoDS_Shape point1 = HSF::GetLowestVertex(p2);
if (point1.IsNull()) continue;
gp_Pnt intp2 = BRep_Tool::Pnt(TopoDS::Vertex(point1));

TopoDS_Shape sphere = HSF::AddNewSphereSurface(intp2,minrad);
//vis(sphere);

TopoDS_Shape intsphere2 = HSF::AddNewIntersectSrfW(SRF,sphere);
//TopoDS_Shape intsphere2 = HSF::AddNewIntersectFast(SRF,sphere);
nextspheres vis(intsphere2);

QString line;
QTextStream linestream(&line);
linestream ui::getInstance()->Statusbar(line);

} catch (...) {}
END
prevspheres = nextspheres;
END

}

catch(...)
{

}

ENDPART

qGeomApp->myview->viewFront();
qGeomApp->myview->fitAll();

double milliseconds = time.elapsed();
double seconds = (milliseconds/1000);
double minutes = seconds/60;
double hours = minutes/60;

QString timetxt;
QTextStream linestream(&timetxt);
linestream

ui.time->setText(timetxt);
}

AP's picture

I have also tried this intersection code, but the intersection of a sphere and a a surface, even when they only have one face each, it doesnt produce a full circle, instead a partial shape.

TopoDS_Shape HSF::AddNewIntersectSrf(TopoDS_Shape srf1,TopoDS_Shape srf2)
{
if (srf1.ShapeType() != TopAbs_ShapeEnum::TopAbs_FACE)
{
srf1 = HSF::getfacefromshape(srf1);
}
if (srf2.ShapeType() != TopAbs_ShapeEnum::TopAbs_FACE)
{
srf2 = HSF::getfacefromshape(srf2);
}

Handle(Geom_Surface) S1 = BRep_Tool::Surface(TopoDS::Face(srf1));
Handle(Geom_Surface) S2 = BRep_Tool::Surface(TopoDS::Face(srf2));

//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.
BRepBuilderAPI_MakeEdge edgebuilder(C,C->FirstParameter(),C->LastParameter());

return edgebuilder.Shape();
}

this is how i build the sphere with one face:

TopoDS_Shape HSF::AddNewSphereSurface(gp_Pnt center, double radius)
{

Handle(Geom_SphericalSurface) SP = new Geom_SphericalSurface(gp_Ax3(center,gp::DZ()),radius);

Standard_Real u1, u2, v1, v2;
SP->Bounds(u1,u2,v1,v2);
fixParam(u1);
fixParam(u2);
fixParam(v1);
fixParam(v2);

TopoDS_Shape Result = BRepBuilderAPI_MakeFace(SP, u1, u2, v1, v2);

return Result;

}
else
{
TopoDS_Shape nullshape;
return nullshape;
}
}

with two faces:

TopoDS_Shape HSF::AddNewSphereSurface(gp_Pnt center, double radius)
{
Handle(Geom_SphericalSurface) SP = new Geom_SphericalSurface(gp_Ax3(center,gp::DZ()),radius);

Standard_Real u1, u2, v1, v2;
SP->Bounds(u1,u2,v1,v2);
fixParam(u1);
fixParam(u2);
fixParam(v1);
fixParam(v2);

TopoDS_Shape Result1 = BRepBuilderAPI_MakeFace(SP, u1, u2/2, v1, v2);
TopoDS_Shape Result2 = BRepBuilderAPI_MakeFace(SP, u2/2, u2, v1, v2);

TopoDS_Compound folder;
BRep_Builder B;
B.MakeCompound(folder);
B.Add(folder,Result1);
B.Add(folder,Result2);

return folder;
}

Evgeny Lodyzhehsky's picture

Dear Alex.
1. Let me try to analyse your function:
TopoDS_Shape HSF::AddNewIntersectFast(TopoDS_Shape srf1,TopoDS_Shape srf2)
{...}

1.1
The two for statements are used to obtain the pair of faces aF1/aF2 to intersect. OK. But what the performance will be when we have say 1000 faces in srf1, srf2? The things will be better if you use at least some preparing procedure to select the pairs of faces aF1/aF2. The selection can be done using the intersection between bounding boxes of source faces (The algorithm of unbalanced binary tree of overlapped bounding boxes: NCollection_UBTree, NCollection_UBTreeFiller). You can get to know about the details of the technique for e.g. in method
void NMTDS_Iterator::Intersect() /see open source SALOME, GEOM module/.

1.2.
aFF.SetParameters (true,true,true,anApproxTol);
As I understood you needs is the section. So, it is not necessary to compute p-curves on faces, i.e.:
aFF.SetParameters (Standard_True,Standard_False,Standard_False,anApproxTol);

1.3.
aFF.Perform(aF1, aF2);
The algorithm IntTools_FaceFace does not intersect aF1/aF2. It provides the intersection between the rectangular pathces of surfaces of aF1,aF2
//...
BRepTools::UVBounds(aF1, umin, umax, vmin, vmax);
//...
Please, have a look IntTools_FaceFace.cxx for more details. Thus the result of aFF.Perform(aF1, aF2) is just intersection curves Ci (i=1,2...nbC) and intersection points Pj (j=1,2...nbP) but not edges, vertices.
The curves Ci have the extensions that correspond to UV-bounds of aF1,aF2 but not correspond to the real bounds of aF1,aF2. Thus the edges you've made:
//...
TopoDS_Edge myedge = BRepBuilderAPI_MakeEdge(mycrv);
//...
will not be real section edges for the faces aF1, aF2.
In order to obtain correct result of the section operation it is necessary to
a) trim Ci by the real bounds of aF1,aF2 -> CTi;
b) make edges from CTi -> ESi;
c) provide sharing between the vertices of ESi where it is should be;
The points a),c) can be done if there is information about all intersecions (interferences) of order lower than F/F i.e. (V/V, V/E, V/F, E/E, E/F). That is why the interferences above are computed before F/F.

2.
"...who invented this technique?(Ivan Sutherland?, University of Utah?) is it documented on any siggraph papers? or any other proceedings..."

Try to find the same published information,docs,etc for Catia for e.g. and send me the link in sase of success, please.

AP's picture

Dear Lodyzhensky Evgeny Nicolaich

Thank You so much for the Response!

I will look into the suggestions.
I agree it will be difficult to find relevant information regarding the technique, especially if its from Dassault.

I will try my best, and keep you posted nonetheless...

Thanks for the hints.

Best,

Alexander Pena de Leon