Implement ability to construct B-Splines from Poles, Knots, Multiplicities, and Degree parameters.

[nomad-vagabond]

Existing spline construction methods assume that the user passes an array of points to interpolate/approximate. This approach, however, does not provide control over interpolation methods, and in some cases fails to satisfy expected smoothness.

Example. For the airfoil geometry, I often obtain wavy shapes due to overfitting (I use Selig format for the profile points). A possible solution would be to construct the desired shape using external spline approx libraries or services (SciPy, NurbsPython, SplineCloud) and pass spline-supporting data (control points, knots, degree) directly into CadQuery to end up with the shape of desired smoothness.

OCC supports this type of the spline creation via OCC.Geom.Geom_BSplineCurve class.

Links to OCCT documentation:

• https://liuxinwin_admin.gitee.io/pythonocc-docs/OCC.Geom.html#OCC.Geom.Geom_BSplineCurve

• https://dev.opencascade.org/doc/refman/html/class_geom___b_spline_curve.html#a00d969dca532b7cd71594f805608fc7d

[adam-urbanczyk]

First take a look at cq.Shape.makeSplineApprox and cq.Workplane.splineApprox.

[nomad-vagabond]

First take a look at cq.Shape.makeSplineApprox and cq.Workplane.splineApprox.

I see no way how these methods are related to the proposed functionality.

SplineApprox is supposed to be used for fitting data points. The proposed method is about enabling the construction of B-Splines based on control points, knot vector, and degree (ideally weights too).

[adam-urbanczyk]

Example. For the airfoil geometry, I often obtain wavy shapes due to overfitting (I use Selig format for the profile points). A possible solution would be to construct the desired shape using external spline approx libraries or services (SciPy, NurbsPython, SplineCloud) and pass spline-supporting data (control points, knots, degree) directly into CadQuery to end up with the shape of desired smoothness.

You mentioned external approximation libraries and airfoils yourself. I'd use the functions above in those contexts.

[nomad-vagabond]

Due to the makeSplineApprox and splineApprox limitations, I can not use them for this particular case (no option to add data point weights). Please, consider the described case as only one of many possible cases when users might want to construct B-Splines with given parameters.

Moreover - the requested method will enable users to write custom methods for spline construction. Imagine a pipeline connector or other similar cases when you need to define spline by control points without any data points to approximate/interpolate.

[adam-urbanczyk]

Can you describe your actual use case precisely? BTW you can always use OCCT methods directly.

[nomad-vagabond]

I do use OCCT methods right now: https://github.com/nomad-vagabond/cq-uav/blob/main/cquav/wing/profile.py#L15

However, I think it should be implemented as a CadQuery method.

My use case: I'm trying to build smooth airfoil shapes using data points in this format: https://splinecloud.com/compose/dfl_lBoz4QrH9wgJ/datasets/dst_h2IXUCbZS3l3

I can not use splineApprox due to: 1) lack of flexibility; 2) exception when trying to fit profile points in the given format.

Here is how I find the best fit with SciPy: https://github.com/nomad-vagabond/cq-uav/blob/main/cquav/wing/profile.py#L104

And weights parameter is crucial for obtaining the desired shape.

Then I reconstruct the spline with the OCCT methods as mentioned before.

[adam-urbanczyk]

This is what I get with splineApprox (see below). I don't see anything obviously wrong. I also don't understand the story with weights, your profile is supposed to be exact after all. If you want to use different profile, simply use a different profile.

What you say you need is more or less implemented here https://github.com/CadQuery/cadquery/blob/79e64e557e87d63b84c1c8a60c0df8e941a1a4a1/cadquery/occ_impl/importers/dxf.py#L61 but with the dxf import use case in mind. You could refactor it into e.g. occ_impl.importers.utils and enable it to accept nurbs control points knots and weights as e.g. iterators or sequences. You wrote some code already, so PRs are welcome. There will be some review process though.

import numpy as np
import cadquery as cq

data = np.loadtxt(r'naca23015_Subset_1.csv',skiprows=1,delimiter=',')
w = cq.Workplane().splineApprox(data[:,[1,2]]).close()

pts = [cq.Vertex.makeVertex(*d,0) for d in data[:,[1,2]]]

show_object(w)
show_object(pts)

[nomad-vagabond]

Wow, thanks for digging deep into my problem!

Must confess now, that the CQ spline fitting methods do work with airfoil data (do not know why I got the exception last time). However, since we have this discussion, let me explain the drawbacks of these methods in more detail.

I've been into spline fitting for quite some time and even created a dedicated platform for spline fitting (SplineCloud). As with many other fitting methods, two issues hurt the most: overfitting and underfitting. And this is also the case for airfoil shape approximation. There are a bunch of research papers dedicated to this problem, here are just a few of them:

Universal Airfoil Parametrization Using B-Splines

Parametric Airfoil Catalog

A Comparison of Airfoil Shape Parameterization Techniques

Smoothing of the airfoil section having curvature graph oscillations

As I mentioned before, the CQ methods (spline of the Sketch class, makeSplineApprox'of the Edge class and splineApprox of the Workplane class) are quite generic and lack advanced fitting parameters. Moreover, there is no limit to ingenuity when it comes to development of advanced fitting techniques, and one framework just can not handle all possible options. In this regard, restricting users to only given functionality binds their hands when constructing complex shapes programmatically.

I deed some analysis to show the obvious issues of overfitting and underfitting that are induced by blindly using default interpolation and approximation with CadQuery.

The airfoil from your example show no visible signs of these issues, so I picked another one from the database that is more susceptible.

splineApprox defaults to interpolation, which causes overfitting:

Setting smoothing to some values makes no effect at all:

bug or it's meaning differs form makeSplineApprox?

makeSplineApprox actually generates in the smoothed spline, but leads to underfitting:

As you can see, by diving deep and trying to optimize airfoil shape one can find out the need to use some external libraries. And this is what I did with SciPy after trying out different other options. It's method splrep allows to set weights to datapoints. I used this trick to set the weight to the nose tip point at (0,0) to 20, forcing spline to pass extremely close to this point. This trick helped me to obtain nice smooth spline approximation.

You can find code with my analyses here: https://github.com/nomad-vagabond/airfoil-spline-approx

[adam-urbanczyk]

OK, your use case is clear and very specific. I can see how having this import functionality would be beneficial for some users. As I mentioned the code is there:

        degree = el.dxf.degree
        periodic = el.closed
        rational = False

        knots_unique = OrderedDict()
        for k in el.knots:
            if k in knots_unique:
                knots_unique[k] += 1
            else:
                knots_unique[k] = 1

        # assmble knots
        knots = TColStd_Array1OfReal(1, len(knots_unique))
        multiplicities = TColStd_Array1OfInteger(1, len(knots_unique))
        for i, (k, m) in enumerate(knots_unique.items()):
            knots.SetValue(i + 1, k)
            multiplicities.SetValue(i + 1, m)

        # assemble weights if present:
        if el.weights:
            rational = True

            weights = TColStd_Array1OfReal(1, len(el.weights))
            for i, w in enumerate(el.weights):
                weights.SetValue(i + 1, w)

        # assemble control points
        pts = TColgp_Array1OfPnt(1, len(el.control_points))
        for i, p in enumerate(el.control_points):
            pts.SetValue(i + 1, gp_Pnt(*p))

        if rational:
            spline = Geom_BSplineCurve(
                pts, weights, knots, multiplicities, degree, periodic
            )
        else:
            spline = Geom_BSplineCurve(pts, knots, multiplicities, degree, periodic)

        return (Edge(BRepBuilderAPI_MakeEdge(spline).Edge()),)

It just needs to be refactored into a separate function along the lines of def toSpline(knots, points, weights, order) -> cq.Edge.

[nomad-vagabond]

I can take this on in a week or two.

[nomad-vagabond]

In the meantime, can you please confirm if there is a bug in the smoothing parameter of the splineApprox. If so - it deserves a separate issue to be opened.

[adam-urbanczyk]

I don't think it is a bug. AFAICT the parameters are relative weights of certain penalty terms, try e.g (1,0,0) vs (0,1,0). I'm afraid OCCT documentation of this function is rather thin, but if you are really interested, you can always read the code.