Arbitrary law-based curve and surface modeling

by - 17:13

A recent case with enhancing CAD Exchanger ACIS importer to broaden support of ACIS primitives inspired me to write this post.

Like other modeling kernels, Open CASCADE comes with a finite set of supported types of curves and surfaces. For instance, for curves this includes lines and conic curves (circles, ellipses, parabolas and hyperbolas), B-Splines, Bezier curves, offset curves plus explicitly trimmed curves. OCC supports parametrics definition where 3D coordinate (x,y,z) is evaluated via parameter t, for surface it is calculated from a pair (u,v).

For instance, for line it is:
C(t) = O + Dt, where O is an origin and D is a unit vector.

Thus, a line in OCC is parametrized by its length.

ACIS and Parasolid additionally support so called procedural geometries (e.g. intersection curve or rolling ball surface) when there is no explicit formula but each point is still unambiguously calculated from parameters t or (u,v).

However, a limited set of explicitly supported types does not allow to express other possible types of curves and surfaces which could be easily represented in parametric definition. For instance, there were a few questions on the OCC forum how to model a helix curve with the help of OCC.

Let us consider how you could that indeed.

Say, a helix is parametrized as follows:

X(t) = O + Rcos(u) X
Y(t) = O + Rsin(u) Y
Z(t) = O + (pitch/(2*PI) * v Z

Where {O ,X, Y, Z} is a local axis system (origin and 3 unit vectors). Pitch – is a helix step along the Z axis. See the below screenshot (taken from ACIS documentation):


Now given this explicit definition (or law) how could you map this to Open CASCADE ?

One option would be to subclass Geom_Curve and implement respective methods (D0(), D1(), ..., Continuity(), etc). That would be fine if you only need a limited time span of such an object and won’t feed it into various OCC modeling algorithms. The OCC classes ShapeExtend_ComplexCurve and _CompositeSurface follow this approach.

However, although you could create an edge on such a curve, you would not be able to save it in a .brep file right away (you would need to enable a vehicle for saving/retrieving user-defined types). Some algorithms may throw an exception on an unrecognized type and so on.

Another option is to do one time approximation with B-Spline and keep the B-Spline representation after that. You would certainly lose original definition and evaluation of a point and derivative would be times more expensive. However you could keep this representation in persistent representation and confidently use it anywhere. By the way ACIS does support helix type and combines both the original definition and its B-Spline approximation.

To approximate with B-Spline, the GeomConvert_ApproxCurve class will do the job.

GeomConvert_ApproxCurve accepts a curve adaptor which essentially implements an Adaptor design pattern and produces a B-Spline curve. You would essentially need to create a subclass of Adaptor3d_Curve and implement respective methods (D0(), D1(), Intervals(), etc). You could possibly mix this approach with the former and create GeomAdaptor_Curve which would accept your Geom_Curve subclass.

Below is an excerpt from CAD Exchanger that approximates a helix with OCC B-Spline:
/*! Uses adaptor classes and invokes GeomConvert_ApproxCurve to approximate with a B-Spline.
    Created B-Spline is polynomial and is of C2-continuity.
    Returns true if the B-Spline has been successfully created and false otherwise.
*/
bool ACISAlgo_Helix::MakeHelix (const ACISGeom_HelixData& theSource,
    Handle_Geom_BSplineCurve& theTarget)
{
    Handle_ACISAlgo_HHelixCurveAdaptor anAdaptor = new ACISAlgo_HHelixCurveAdaptor (theSource);
    Standard_Real aTol = Precision::Confusion();
    GeomAbs_Shape aContinuity = GeomAbs_C2 /*highest supported continuity*/;
    Standard_Integer aMaxSeg = 10000, /*max number of spans*/
                     aMaxDeg = 9; /*max degree, consistent with settings in Algo*/
    GeomConvert_ApproxCurve anApprox (anAdaptor, aTol,
        aContinuity,
        aMaxSeg,
        aMaxDeg);
    if (anApprox.HasResult()) {
        theTarget = anApprox.Curve();
        Base_Debugger& aDebugger = Base_Debugger::GlobalInstance();
        aDebugger.Save (theTarget, "curve");
    }
    return !theTarget.IsNull();
 }
//! Defines data elements that can be reused for both common ACIS 'helix' and ASM 'helix_int_cur'.
struct ACISGeom_HelixData
{
    ACISGeom_HelixData() :
        myXRadius(0.),
        myYRadius(0.),
        myPitch(0.),
        myTaper(0.),
        myRangeMin (0.),
        myRangeMax (2 * M_PI),
        myScaleFactor (1.)
    {}

    __CADEX_DEFINE_PROPERTY(gp_Ax3,Position) //can be right- or left-handed
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,XRadius)
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,YRadius)
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,Pitch)  //must be >= 0, if = 0 then is planar
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,Taper) //if > 0, helix widens along the Z-axis, if 0 - then lies on a cylinder
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,RangeMin)
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,RangeMax)
    __CADEX_DEFINE_PRIMITIVE_PROPERTY(double,ScaleFactor)
};

/*! A few methods in OCC 6.9.0 have been made const.*/
#if OCC_VERSION_HEX < 0x060900
#define __CADEX_ADAPTOR3D_CURVE_CONST
#else
#define __CADEX_ADAPTOR3D_CURVE_CONST const
#endif

/* \class ACISAlgo_HelixCurveAdaptor
   \brief Defines an adaptor to represent a helix curve.

   Helix data is defined by ACISGeom_HelixData.
   Evaluation is performed in the Evaluator subclass which can either represent a helix lying on a
   cylinder (if taper is 0) or on a cone (if taper is not 0).

   Helix can have distinct radii along X and Y axes, i.e. to have an elliptical section.
*/
class ACISAlgo_HelixCurveAdaptor : public Adaptor3d_Curve
{
public:

    class Evaluator;

    //! Constructor.
    ACISAlgo_HelixCurveAdaptor (const ACISGeom_HelixData& theData);

    //! Constructor
    ACISAlgo_HelixCurveAdaptor (const std::shared_ptr& theEvaluator,
        Standard_Real theMin,
        Standard_Real theMax);

    virtual Standard_Real FirstParameter()const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual Standard_Real LastParameter() const __CADEX_OVERRIDE_ATTRIBUTE;

    virtual GeomAbs_Shape Continuity() const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual Standard_Integer NbIntervals (const GeomAbs_Shape S) __CADEX_ADAPTOR3D_CURVE_CONST
        __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void Intervals (TColStd_Array1OfReal& T, const GeomAbs_Shape S) __CADEX_ADAPTOR3D_CURVE_CONST
        __CADEX_OVERRIDE_ATTRIBUTE;
    virtual Handle(Adaptor3d_HCurve) Trim (const Standard_Real First,
        const Standard_Real Last,
        const Standard_Real Tol) const __CADEX_OVERRIDE_ATTRIBUTE;
  
    virtual Standard_Boolean IsClosed() const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual Standard_Boolean IsPeriodic() const __CADEX_OVERRIDE_ATTRIBUTE;

    virtual gp_Pnt Value (const Standard_Real U) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D0 (const Standard_Real U, gp_Pnt& P) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D1 (const Standard_Real U, gp_Pnt& P, gp_Vec& V) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D2 (const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D3 (const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2, gp_Vec& V3) const
        __CADEX_OVERRIDE_ATTRIBUTE;
    virtual gp_Vec DN (const Standard_Real U, const Standard_Integer N) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual Standard_Real Resolution (const Standard_Real R3d) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual GeomAbs_CurveType GetType() const __CADEX_OVERRIDE_ATTRIBUTE;

protected:
    std::shared_ptr  myEvaluator;
    Standard_Real               myMin;
    Standard_Real               myMax;
};
DEFINE_STANDARD_HANDLE(ACISAlgo_HHelixCurveAdaptor,Adaptor3d_HCurve)
class ACISAlgo_HHelixCurveAdaptor : public Adaptor3d_HCurve
{
public:

    //! Constructor.
    ACISAlgo_HHelixCurveAdaptor (const ACISGeom_HelixData& theData) : myAdaptor (theData) {}

    //! Constructor.
    ACISAlgo_HHelixCurveAdaptor (const std::shared_ptr& theEvaluator,
        Standard_Real theMin,
        Standard_Real theMax) : myAdaptor (theEvaluator, theMin, theMax) {}

    //! Returns the adaptor as Adaptor3d_Curve.
    /*! Return the internal ACISAlgo_HelixCurveAdaptor object.*/
    virtual const Adaptor3d_Curve& Curve() const __CADEX_OVERRIDE_ATTRIBUTE { return myAdaptor; }

    //! Returns the adaptor as Adaptor3d_Curve.
    /*! Return the internal ACISAlgo_HelixCurveAdaptor object.*/
    virtual Adaptor3d_Curve& GetCurve() __CADEX_OVERRIDE_ATTRIBUTE  { return myAdaptor; }
  
public:
    DEFINE_STANDARD_RTTI(ACISAlgo_HelixCurveAdaptor)

protected:
    ACISAlgo_HelixCurveAdaptor  myAdaptor;
};

__CADEX_IMPLEMENT_HANDLE(ACISAlgo_HHelixCurveAdaptor,Adaptor3d_HCurve)
/*********************************************************************************************/


/*! \class ACISAlgo_HelixCurveAdaptor::Evaluator
    \brief Base abstract class to evaluate helix.
*/
class ACISAlgo_HelixCurveAdaptor::Evaluator
{
public:
    __CADEX_DEFINE_MEMORY_MANAGEMENT

    Evaluator (const ACISGeom_HelixData& theData) : myData (theData), myVCoef (1.) {}
    virtual ~Evaluator() {}

    const ACISGeom_HelixData& Data() const { return myData; }

    double VParameter (Standard_Real U) const { return U * myVCoef; }
    virtual void D0 (Standard_Real U, gp_Pnt& P) const = 0;
    virtual void D1 (Standard_Real U, gp_Pnt& P, gp_Vec& V) const = 0;
    virtual void D2 (Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const = 0;
    virtual void D3 (Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2, gp_Vec& V3) const = 0;
    virtual gp_Vec DN (Standard_Real U, Standard_Integer N) const = 0;

protected:
    const gp_XYZ& Loc()  const { return myData.Position().Location().XYZ(); }
    const gp_XYZ& XDir() const { return myData.Position().XDirection().XYZ(); }
    const gp_XYZ& YDir() const { return myData.Position().YDirection().XYZ(); }
    const gp_XYZ& ZDir() const { return myData.Position().Direction().XYZ(); }

    ACISGeom_HelixData  myData;
    double              myVCoef; //coefficient to multiply U to get a V parameter on a respective surface
};
/*! \class ACISAlgo_HelixCurveAdaptor_CylinderEvaluator
    \brief Evaluates a helix lying on a cylinder.
*/
class ACISAlgo_HelixCurveAdaptor_CylinderEvaluator : public ACISAlgo_HelixCurveAdaptor::Evaluator
{
public:
    ACISAlgo_HelixCurveAdaptor_CylinderEvaluator (const ACISGeom_HelixData& theData);
    virtual void D0 (Standard_Real U, gp_Pnt& P) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D1 (Standard_Real U, gp_Pnt& P, gp_Vec& V) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D2 (Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual void D3 (Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2, gp_Vec& V3) const __CADEX_OVERRIDE_ATTRIBUTE;
    virtual gp_Vec DN (Standard_Real U, Standard_Integer N) const __CADEX_OVERRIDE_ATTRIBUTE;
};

ACISAlgo_HelixCurveAdaptor_CylinderEvaluator::ACISAlgo_HelixCurveAdaptor_CylinderEvaluator (
    const ACISGeom_HelixData& theData) :
    ACISAlgo_HelixCurveAdaptor::Evaluator (theData)
{
    myVCoef = theData.Pitch() * theData.ScaleFactor() / (2 * M_PI);
}


void ACISAlgo_HelixCurveAdaptor_CylinderEvaluator::D0 (Standard_Real U,
    gp_Pnt& P) const
{
    Standard_Real v    = VParameter (U);
    Standard_Real Rx   = myData.XRadius();
    Standard_Real Ry   = myData.YRadius();
    Standard_Real sinU = sin (U);
    Standard_Real cosU = cos (U);
    P = Rx * cosU * XDir() + Ry * sinU * YDir() + v * ZDir() + Loc();
}

void ACISAlgo_HelixCurveAdaptor_CylinderEvaluator::D1 (
    Standard_Real U,
    gp_Pnt& P,
    gp_Vec& V) const
{
    Standard_Real v    = VParameter (U);
    Standard_Real Rx   = myData.XRadius();
    Standard_Real Ry   = myData.YRadius();
    Standard_Real sinU = sin (U);
    Standard_Real cosU = cos (U);
    P = Rx * cosU * XDir() + Ry * sinU * YDir() + v * ZDir() + Loc();

    Standard_Real k    = myVCoef;
    V = -Rx * sinU * XDir() + Ry * cosU * YDir() + k * ZDir();
}

void ACISAlgo_HelixCurveAdaptor_CylinderEvaluator::D2 (
    Standard_Real U,
    gp_Pnt& P,
    gp_Vec& V1,
    gp_Vec& V2) const
{
    Standard_Real v    = VParameter (U);
    Standard_Real Rx   = myData.XRadius();
    Standard_Real Ry   = myData.YRadius();
    Standard_Real sinU = sin (U);
    Standard_Real cosU = cos (U);
    P = Rx * cosU * XDir() + Ry * sinU * YDir() + v * ZDir() + Loc();

    Standard_Real k    = myVCoef;
    V1 = -Rx * sinU * XDir() + Ry * cosU * YDir() + k * ZDir();
    V2 = -Rx * cosU * XDir() - Ry * sinU * YDir();
}
/*********************************************************************************************/

ACISAlgo_HelixCurveAdaptor::ACISAlgo_HelixCurveAdaptor (const ACISGeom_HelixData& theData) :
    myMin (theData.RangeMin()),
    myMax (theData.RangeMax())
{
    if (std::fabs (theData.Taper()) < Precision::Confusion()) { //cylinder
        myEvaluator.reset (new ACISAlgo_HelixCurveAdaptor_CylinderEvaluator (theData));

    } else { //cone
        myEvaluator.reset (new ACISAlgo_HelixCurveAdaptor_ConeEvaluator (theData));
    }
}

/*! Used when trimming*/
ACISAlgo_HelixCurveAdaptor::ACISAlgo_HelixCurveAdaptor (const std::shared_ptr& theEvaluator,
    Standard_Real theMin,
    Standard_Real theMax) :
    myEvaluator (theEvaluator),
    myMin (theMin),
    myMax (theMax)
{
    __CADEX_ASSERT_INVALID_VALUE(myEvaluator->Data().RangeMin() <= theMin);
    __CADEX_ASSERT_INVALID_VALUE(myEvaluator->Data().RangeMax() >= theMax);
}

Standard_Real ACISAlgo_HelixCurveAdaptor::FirstParameter() const
{
    return myMin;
}

Standard_Real ACISAlgo_HelixCurveAdaptor::LastParameter() const
{
    return myMax;
}

GeomAbs_Shape ACISAlgo_HelixCurveAdaptor::Continuity() const
{
    return GeomAbs_CN;
}

Standard_Integer ACISAlgo_HelixCurveAdaptor::NbIntervals (const GeomAbs_Shape /*S*/) __CADEX_ADAPTOR3D_CURVE_CONST
{
    return 1;
}

void ACISAlgo_HelixCurveAdaptor::Intervals (TColStd_Array1OfReal& T, const GeomAbs_Shape /*S*/) __CADEX_ADAPTOR3D_CURVE_CONST
{
    T (T.Lower()) = FirstParameter();
    T (T.Upper()) = LastParameter();
}

Handle(Adaptor3d_HCurve) ACISAlgo_HelixCurveAdaptor::Trim (const Standard_Real First,
    const Standard_Real Last,
    const Standard_Real /*Tol*/) const
{
    return new ACISAlgo_HHelixCurveAdaptor (myEvaluator, First, Last);
}

Standard_Boolean ACISAlgo_HelixCurveAdaptor::IsClosed() const
{
    return Standard_False;
}

Standard_Boolean ACISAlgo_HelixCurveAdaptor::IsPeriodic() const
{
    return Standard_False;
}

gp_Pnt ACISAlgo_HelixCurveAdaptor::Value (const Standard_Real U) const
{
    gp_Pnt aP;
    D0 (U, aP);
    return aP;
}

void ACISAlgo_HelixCurveAdaptor::D0 (const Standard_Real U, gp_Pnt& P) const
{
    myEvaluator->D0 (U, P);
}

void ACISAlgo_HelixCurveAdaptor::D1 (const Standard_Real U, gp_Pnt& P, gp_Vec& V) const
{
    myEvaluator->D1 (U, P, V);
}

void ACISAlgo_HelixCurveAdaptor::D2 (const Standard_Real U, gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) const
{
    myEvaluator->D2 (U, P, V1, V2);
}


Standard_Real ACISAlgo_HelixCurveAdaptor::Resolution (const Standard_Real R3d) const
{
    //see GeomAdaptor_Curve::Resolution()
    const auto& aData = myEvaluator->Data();
    Standard_Real R = std::max (aData.XRadius(), aData.YRadius());
    if (R3d < 2 * R)
        return 2 * ASin (R3d / (2 * R));
    else
        return 2 * M_PI;
}

GeomAbs_CurveType ACISAlgo_HelixCurveAdaptor::GetType() const
{
    return GeomAbs_OtherCurve;
}


Below are two screenshot of approximated helices in DRAW – the first one has different Rx and Ry radii and is lying on a cylindrical surface, the second is of equal Rx and Ry radii lying on a conical surface.



Hope this post will give you some hints on how to approximate arbitrary curves and surfaces using B-Spline approximation techniques.

If you have any interesting examples to share that would certainly be good to know.

Thanks,
Roman

 

You May Also Like

0 comments