NURBS样条曲线学习

搞3D几何内核算法研究，必须学习NURBS样条曲线曲面。

看《非均匀有理B样条  第2版》这本书，学习起来，事半功倍。

在《插件化算法研究平台》上，做了一个样条曲线研究功能，可以分析Bezier曲线、BSpline、NURBS曲线的各种性质，有直观的体验，能更好地理解。

一、基函数

参考示例代码：

// 计算基函数值
    double basis(int i, int k, double u) {
        if (k == 0) {
            if (u >= knots[i] && u < knots[i + 1]) {
                return 1.0;
            }
            return 0.0;
        }
        double a = 0.0, b = 0.0;
        if (knots[i + k] - knots[i] != 0.0) {
            a = (u - knots[i]) / (knots[i + k] - knots[i]);
        }
        if (knots[i + k + 1] - knots[i + 1] != 0.0) {
            b = (knots[i + k + 1] - u) / (knots[i + k + 1] - knots[i + 1]);
        }
        return a * basis(i, k - 1, u) + b * basis(i + 1, k - 1, u);
    }

二、B样条曲线

参考示例代码：

// 计算样条曲线上的点
   virtual Point3  calculatePoint(double u) {
        auto size = controlPoints3d.size();
        Point3 result;
        for (size_t i = 0; i < size; i++) {
            double factor = basis(i, degree, u) ;
            result.x += factor * controlPoints3d[i].x;
            result.y += factor * controlPoints3d[i].y;
            result.z += factor * controlPoints3d[i].z;
        }
        return result;
    }

三、NURBS曲线

参考示例代码：

 // 计算样条曲线上的点
    Point3  calculatePoint(double u)override
    {
        auto size = controlPoints3d.size();
        Point3 result;
        double denominator = 0.0;
        for (size_t i = 0; i < size; i++) {
            double factor = basis(i, degree, u) * weights[i];
            denominator += factor;
            result.x += factor * controlPoints3d[i].x;
            result.y += factor * controlPoints3d[i].y;
            result.z += factor * controlPoints3d[i].z;
        }
        result.x /= denominator;
        result.y /= denominator;
        result.z /= denominator;
        return result;
    }

四、程序效果

Bezier曲线

BSpline曲线

其中，Degree对曲线的影响，当Degree=1时，就成了线性插值。

NURBS曲线

     

权重对曲线的影响。

四、基于OCCT几何内核生成的样条曲线示例

参考示例代码：

//BSpline Bezier  Curve
void makeCurve(OccView* view)
{
    myObject3d.Clear();
    // Define points.
    gp_Pnt aPnt1(0.0, 0.0, 0.0);
    gp_Pnt aPnt2(5.0, 5.0, 0.0);
    gp_Pnt aPnt3(10.0, 5.0, 0.0);
    gp_Pnt aPnt4(15.0, 0.0, 5.0);

    // Add points to the curve poles array.
    TColgp_Array1OfPnt aPoles(1, 4);
    aPoles.SetValue(1, aPnt1);
    aPoles.SetValue(2, aPnt2);
    aPoles.SetValue(3, aPnt3);
    aPoles.SetValue(4, aPnt4);

    // Define BSpline weights.
    TColStd_Array1OfReal aBSplineWeights(1, 4);
    aBSplineWeights.SetValue(1, 1.0);
    aBSplineWeights.SetValue(2, 0.5);
    aBSplineWeights.SetValue(3, 0.5);
    aBSplineWeights.SetValue(4, 1.0);

    // Define knots.
    TColStd_Array1OfReal aKnots(1, 2);
    aKnots.SetValue(1, 0.0);
    aKnots.SetValue(2, 1.0);

    // Define multiplicities.
    TColStd_Array1OfInteger aMults(1, 2);
    aMults.SetValue(1, 4);
    aMults.SetValue(2, 4);

    // Define BSpline degree and periodicity.
    Standard_Integer aDegree = 3;
    Standard_Boolean aPeriodic = Standard_False;

    // Create a BSpline curve.
    Handle(Geom_BSplineCurve) aBSplineCurve = new Geom_BSplineCurve(
        aPoles, aBSplineWeights, aKnots, aMults, aDegree, aPeriodic);

    Handle(AIS_ColoredShape) anAisBSplineCurve = new AIS_ColoredShape(
        BRepBuilderAPI_MakeEdge(aBSplineCurve).Shape());
    anAisBSplineCurve->SetColor(Quantity_Color(Quantity_NOC_IVORY));

    // Compute BSpline surface bounding box.
    Bnd_Box aBndBox;
    BndLib_Add3dCurve::AddOptimal(GeomAdaptor_Curve(aBSplineCurve), Precision::Confusion(), aBndBox);
    Handle(AIS_ColoredShape) anAisBndBox = new AIS_ColoredShape(
        BRepPrimAPI_MakeBox(aBndBox.CornerMin(), aBndBox.CornerMax()).Shell());
    view->aisInterContext->SetDisplayMode(anAisBndBox, 0, false);
    myObject3d.Append(anAisBndBox);

    NCollection_List<Standard_Real> aParams;
    aParams.Append(0.75 * aBSplineCurve->FirstParameter() + 0.25 * aBSplineCurve->LastParameter());
    aParams.Append(0.50 * aBSplineCurve->FirstParameter() + 0.50 * aBSplineCurve->LastParameter());
    aParams.Append(0.25 * aBSplineCurve->FirstParameter() + 0.75 * aBSplineCurve->LastParameter());
    for (NCollection_List<Standard_Real>::Iterator anIt(aParams); anIt.More(); anIt.Next())
    {
        Standard_Real aParam = anIt.Value();
        gp_Pnt aPnt;
        gp_Vec aVec;
        aBSplineCurve->D1(aParam, aPnt, aVec);
        myObject3d.Append(new AIS_Point(new Geom_CartesianPoint(aPnt)));
        Handle(AIS_TextLabel) aisLabel = new AIS_TextLabel();
        Standard_SStream aSS;
        aSS << "P [" << aPnt.X() << ", " << aPnt.Y() << ", " << aPnt.Z() << "]" << std::endl;
        aSS << "D [" << aVec.X() << ", " << aVec.Y() << ", " << aVec.Z() << "]" << std::endl;
        aisLabel->SetHeight(fontHeight);
        aisLabel->SetText(aSS.str().c_str());
        aisLabel->SetPosition(aPnt);
        myObject3d.Append(aisLabel);
        Handle(AIS_Axis) anAisD = new AIS_Axis(new Geom_Axis1Placement(gp_Ax1(aPnt, aVec)));
        myObject3d.Append(anAisD);
    }


    // Define Bezier weights.
    TColStd_Array1OfReal aBezierWeights(1, 4);
    aBezierWeights.SetValue(1, 0.5);
    aBezierWeights.SetValue(2, 1.5);
    aBezierWeights.SetValue(3, 1.5);
    aBezierWeights.SetValue(4, 0.5);

    // Create Bezier curve.
    Handle(Geom_BezierCurve) aBezierCurve = new Geom_BezierCurve(aPoles, aBezierWeights);
    Handle(AIS_ColoredShape) anAisBezierCurve = new AIS_ColoredShape(BRepBuilderAPI_MakeEdge(aBezierCurve).Shape());
    anAisBezierCurve->SetColor(Quantity_Color(Quantity_NOC_BLUE));


    myObject3d.Append(anAisBSplineCurve);
    myObject3d.Append(anAisBezierCurve);
    myObject3d.Append(new AIS_Point(new Geom_CartesianPoint(aPnt1)));
    myObject3d.Append(new AIS_Point(new Geom_CartesianPoint(aPnt2)));
    myObject3d.Append(new AIS_Point(new Geom_CartesianPoint(aPnt3)));
    myObject3d.Append(new AIS_Point(new Geom_CartesianPoint(aPnt4)));

    Handle(AIS_TextLabel) aisBSplineLabel = new AIS_TextLabel();
    aisBSplineLabel->SetHeight(fontHeight);
    aisBSplineLabel->SetText("BSpline edge");
    aisBSplineLabel->SetPosition(aPnt4);
    aisBSplineLabel->SetColor(Quantity_Color(Quantity_NOC_IVORY));
    myObject3d.Append(aisBSplineLabel);

    Handle(AIS_TextLabel) aisLabel = new AIS_TextLabel();
    aisLabel->SetHeight(fontHeight);
    aisLabel->SetText("Bezier edge");
    aisLabel->SetPosition(aPnt3);
    aisLabel->SetColor(Quantity_Color(Quantity_NOC_BLUE));
    myObject3d.Append(aisLabel);

    gp_Ax2 anAxis2(gp_Pnt(0.0, 0.0, 0.0), gp_Dir(0.0, 0.0, 1.0));
    // Describes a parabola 抛物线 in 3D space. A parabola is defined by its focal length
  // (that is, the distance between its focus and apex) and positioned in space with
  // a coordinate system (a gp_Ax2 object)
    gp_Parab aParab(anAxis2.Translated(gp_Vec(0.0, 0.0, 20.0)), 2.0);
    Handle(Geom_Parabola) aGeomParabola = new Geom_Parabola(aParab);
    Handle(Geom_TrimmedCurve) aTrimmedParabola = new Geom_TrimmedCurve(aGeomParabola, 20.0, -20.0);
    Handle(AIS_ColoredShape) anAisParabola = new AIS_ColoredShape(BRepBuilderAPI_MakeEdge(aTrimmedParabola));
    myObject3d.Append(anAisParabola);

    // Describes a branch of a hyperbola 双曲线 in 3D space. A hyperbola is defined by its major
    // and minor radii and positioned in space with a coordinate system (a gp_Ax2 object)
    gp_Hypr aHypr(anAxis2.Translated(gp_Vec(0.0, 0.0, 30.0)), 20.0, 10.0);
    Handle(Geom_Hyperbola) aGeomHyperbola = new Geom_Hyperbola(aHypr);
    Handle(Geom_TrimmedCurve) aTrimmedHyperbola = new Geom_TrimmedCurve(aGeomHyperbola, 2.0, -2.0);
    Handle(AIS_ColoredShape) anAisHyperbola = new AIS_ColoredShape(BRepBuilderAPI_MakeEdge(aTrimmedHyperbola));
    myObject3d.Append(anAisHyperbola);

    auto context = view->getContext();
    NCollection_Vector<Handle(AIS_InteractiveObject)>::iterator it = myObject3d.begin();
    for (; it != myObject3d.end(); it++)
    {
        context->Display(*it, Standard_False);
    }
    view->fitAll();
}