Generating points on a face

I am trying to generate points on a face:

To this end I have written the following function:

TColgp_Array2OfPnt* generateFacePoints(TopoDS_Face face,int nu_p_iso){

TColgp_Array2OfPnt* points = new TColgp_Array2OfPnt(0,nu_p_iso-1,0,nu_p_iso-1);

BRepAdaptor_Surface *surface = new BRepAdaptor_Surface(face,true);
double u_start = surface->FirstUParameter();
double u_end = surface->LastUParameter();

double v_start = surface->FirstVParameter();
double v_end = surface->LastVParameter();

double du = (u_end - u_start)/(nu_p_iso-1);
double dv = (v_end - v_start)/(nu_p_iso-1);

for (int i=0; i for (int j=0; j double u = u_start + i*du;
double v = v_start + j*dv;

gp_Pnt point = surface->Value(u,v);
points->SetValue(i,j,point);
}
}
return points;
}

What it does is that it takes the face where the points are to be generated and then increments the surface in the parametric space in order to generate the points. The problem
with this implementation is that the points that are generated are often out of the surface.
This problem is observed for surfaces that are the starting and ending surfaces of a solid created by revolution of an initial face. How can I overcome this problem?

Aris Karagiannidis's picture

Looking at the problem further, the observation is that the case where the points
that are generated do not lie on the face is the case where the face is a plane. And not a curved surface.

jelle's picture

This could happen when working with trimmed surfaces. Than the parametric domain != to that of the topology. Do perhaps testing if that parameter represents a trimmed area is an idea

Aris Karagiannidis's picture

Hello jelle and thank you for your reply. This prety much happens on every flat surface.

Can you suggest a way to check that the parameter represents a trimmed area on the plane?

Mark Blome's picture

Dear Aris Karagiannidis,

you may use BRepClass_FaceClassifier to decide whether a point
evaluated in the surface of face actually is on the face:

def IsPointOnFace(aFace, P, aTol=0.0001):
from OCC.BRepClass import BRepClass_FaceClassifier
B=BRepClass_FaceClassifier()
B.Perform(TopoDS().face(aFace), P, aTol)
return (B.State()==TopAbs_ON)

There may be more efficient ways of doing this ...

Regards,
Mark

Aris Karagiannidis's picture

Dear Mark. Thank you very much for your response. Indeed the code you are giving here is working. However, while this solution can help filter out the points that are generated and which do not belong in the face that is examined, it does not facilitate the control I would like to have over the number of points I want to produce.

For example, If I want to generate a specific number of points I cannot know beforehand how many of them will be outside my face and therefore the amount of points that will be ultimately generated will be ambiguous.

In any case your help is really appreciated. I also attach a screenshot to clearly demonstrate the problem I am having, just in case someone else would like to pitch in.

Also from your code I can see that given that opencascade can "see" the correct U and V values that bind my surface but does not return them when
I use the methods:

BRepAdaptor_Surface::FirstUParameter()
BRepAdaptor_Surface::LastUParameter()

BRepAdaptor_Surface::FirstVParameter()
BRepAdaptor_Surface::LastVParameter()

Attachments: 
AP's picture

Hi Aris,

this might be overkill but here it goes:

Look into the class: StdPrs_WFRestrictedFace.cxx

check the StdPrs_WFRestrictedFace::Match function, in it, it shows how to use the hatch_hatcher class.

this file shows how to extract the iso curves on a restricted face. Technically it takes the unrestricted iso curves, then trims them against the boundaries of the face, then gets the new iso curves trimmed. in theory, you should be able to evaluate points on the trimmed iso curves and that could be a way to get your points:

post your results as soon as you get something. Best of luck.

Alex

adapt the following code to your purpose, take the anIso curve and evaluate points on them:

TopoDS_Shape theFace = yourface...;
bool theDrawUIso = true;
bool theDrawVIso = true;
integer theNbUIso = 10; //num of u iso curves
integer theNbVIso = 10; //num of v iso curves

StdPrs_ToolRFace aToolRst (theFace);

// Compute bounds of the restriction
Standard_Real anUMin,anUMax,aVMin,aVMax;
Standard_Real anU,aV,aStep;
Standard_Integer anI,anNbP = 10;
anUMin = aVMin = RealLast();
anUMax = aVMax = RealFirst();
gp_Pnt2d aPoint1,aPoint2;

for (aToolRst.Init(); aToolRst.More(); aToolRst.Next())
{
Adaptor2d_Curve2dPtr aRCurve = aToolRst.Value();
anU = aRCurve->FirstParameter();
aV = aRCurve->LastParameter();
if (aRCurve->GetType() != GeomAbs_Line)
{
aStep = ( aV - anU) / anNbP;
for (anI = 0; anI <= anNbP; ++anI)
{
gp_Pnt2d aRCurvePoint = aRCurve->Value(anU);
if (aRCurvePoint.X() < anUMin) anUMin = aRCurvePoint.X();
if (aRCurvePoint.X() > anUMax) anUMax = aRCurvePoint.X();
if (aRCurvePoint.Y() < aVMin) aVMin = aRCurvePoint.Y();
if (aRCurvePoint.Y() > aVMax) aVMax = aRCurvePoint.Y();
anU += aStep;
}
}
else
{
aPoint1 = aRCurve->Value(anU);
if (aPoint1.X() < anUMin) anUMin = aPoint1.X();
if (aPoint1.X() > anUMax) anUMax = aPoint1.X();
if (aPoint1.Y() < aVMin) aVMin = aPoint1.Y();
if (aPoint1.Y() > aVMax) aVMax = aPoint1.Y();

aPoint2 = aRCurve->Value(aV);
if (aPoint2.X() < anUMin) anUMin = aPoint2.X();
if (aPoint2.X() > anUMax) anUMax = aPoint2.X();
if (aPoint2.Y() < aVMin) aVMin = aPoint2.Y();
if (aPoint2.Y() > aVMax) aVMax = aPoint2.Y();
}
}

// Load the isos
Hatch_Hatcher anIsoBuild(1.e-5,aToolRst.IsOriented());
Standard_Boolean anUClosed = theFace->IsUClosed();
Standard_Boolean aVClosed = theFace->IsVClosed();

if ( ! anUClosed )
{
anUMin = anUMin + ( anUMax - anUMin) /1000.;
anUMax = anUMax - ( anUMax - anUMin) /1000.;
}

if ( ! aVClosed )
{
aVMin = aVMin + ( aVMax - aVMin) /1000.;
aVMax = aVMax - ( aVMax - aVMin) /1000.;
}

if (theDrawUIso)
{
if (theNbUIso > 0)
{
anUClosed = Standard_False;
Standard_Real du= anUClosed ? (anUMax-anUMin)/theNbUIso : (anUMax-anUMin)/(1+theNbUIso);
for (anI=1; anI<=theNbUIso;anI++){
anIsoBuild.AddXLine(anUMin+du*anI);
}
}
}
if (theDrawVIso){
if ( theNBVIso > 0) {
aVClosed = Standard_False;
Standard_Real dv= aVClosed ?(aVMax-aVMin)/theNBVIso : (aVMax-aVMin)/(1+theNBVIso);
for (anI=1; anI<=theNBVIso;anI++){
anIsoBuild.AddYLine(aVMin+dv*anI);
}
}
}

// Trim the isos
Standard_Real anU1, anU2, aDU;

for (aToolRst.Init(); aToolRst.More(); aToolRst.Next())
{
TopAbs_Orientation Orient = aToolRst.Orientation();
if (Orient == TopAbs_FORWARD || Orient == TopAbs_REVERSED)
{
Adaptor2d_Curve2dPtr aRCurve = aToolRst.Value();
anU1 = aRCurve->FirstParameter();
anU2 = aRCurve->LastParameter();
if (aRCurve->GetType() != GeomAbs_Line) {
aDU = (anU2-anU1)/(aNbPoints-1);
aPoint2 = aRCurve->Value(anU1);
for (anI = 2; anI <= aNbPoints; anI++) {
anU = anU1 + (anI-1)*aDU;
aPoint1 = aPoint2;
aPoint2 = aRCurve->Value(anU);
if(Orient == TopAbs_FORWARD )
anIsoBuild.Trim(aPoint1,aPoint2);
else
anIsoBuild.Trim(aPoint2,aPoint1);
}
}
else {
aPoint1 = aRCurve->Value(anU1);
aPoint2 = aRCurve->Value(anU2);
if(Orient == TopAbs_FORWARD )
anIsoBuild.Trim(aPoint1,aPoint2);
else
anIsoBuild.Trim(aPoint2,aPoint1);
}
}
}

// Draw the isos

Adaptor3d_IsoCurve anIso;
anIso.Load(theFace);
Standard_Integer aNbLines = anIsoBuild.NbLines();

for (anI = 1; anI <= aNbLines; anI++)
{
Standard_Integer aNbIntervals = anIsoBuild.NbIntervals(anI);
Standard_Real aCoord = anIsoBuild.Coordinate(anI);
for (Standard_Integer j = 1; j <= aNbIntervals; j++)
{
Standard_Real anIsoStart=anIsoBuild.Start(anI,j),anIsoEnd=anIsoBuild.End(anI,j);

anIsoStart = anIsoStart == RealFirst() ? - aLimit : anIsoStart;
anIsoEnd = anIsoEnd == RealLast() ? aLimit : anIsoEnd;

if (anIsoBuild.IsXLine(anI))
anIso.Load(GeomAbs_IsoU,aCoord,anIsoStart,anIsoEnd); // decide which one to keep
anIso.Load(GeomAbs_IsoV,aCoord,anIsoStart,anIsoEnd); // decide which one to keep

// HERE ADD CODE to load points on anIso curve.

}
}

AP's picture
Aris Karagiannidis's picture

Dear Alex. I cannot thank you enough. The code sample you posted was very helpful.