Surface that passes through points defined on grid

[snow54]

After almost a month of research and trials and errors, I have decided to take the following path.

1. Use either of the following methods to create upper and lower surfaces.
   • Create a bivariate spline surface that passes through all the input points in Python (Good in a sense that the surface passes through all the input points but high risk of overshoot)
   • Create a spline surface using control points (Not ideal but can avoid overshoot completely)
2. Trim the surface slightly inside of the boundary.
3. Connect the boundary curve defined by spline with the edge curves of surfaces by ruled surfaces.
4. Use Part Builder to do "Face from edges", "Shell from faces", and "Solid from shell" to create a solid.

I realized that use of surface filling or something that lets FreeCAD create a spline surface by itself produces much more control points than necessary when you have your inputs defined on a rectangular grid, which results in a significantly long computation time and sometimes a strange surface when you have more than 100x100 points. That's why I chose to find control points and knot vectors by other means and pass them to FreeCAD. I did lots of searches on this forum and on the internet but I could not find a sample code to enable it, so I have attached my sample code that creates surfaces by both methods. I hope someone trying to do something similar will find it helpful.

I still have some issues when creating the ruled surfaces to connect the boundary curve and the upper and lower surfaces. A very skewed surface is sometimes created when the boundary curve and the edge curves of upper or lower surfaces are complex. My current workaround is to cut those curves to let FreeCAD know which part of the curves should be joined by a straight line. It would be great if there is more control on the ruled surface.

Anyway, thank you very much for your help. I think I finally got what I wanted.

Code: Select all

import FreeCAD as App
from FreeCAD import Base
import Draft
import Part
import numpy as np
from scipy import interpolate

# Create a sample surface
x=np.linspace(-1,1,10)*np.pi
y=np.linspace(-1,1,10)*np.pi
X, Y=np.meshgrid(x,y)
Z=0.5*(np.sin(X)+2*np.cos(2*Y)+Y)

degree_u=3
degree_v=3

doc = App.newDocument()

#--Forcing the surface to pass through all input points--

# Create spline surface
f=interpolate.RectBivariateSpline(x,y,Z.T)
coeffs=f.get_coeffs()
coeffs=np.reshape(coeffs,Z.T.shape).T
knots=f.get_knots()
# Convert knot vector format
knot_u, mult_u=np.unique(knots[1],return_counts=True)
knot_v, mult_v=np.unique(knots[0],return_counts=True)
# Normalize knot vectors
knot_u=(knot_u-knot_u.min())/(knot_u.max()-knot_u.min())
knot_v=(knot_v-knot_v.min())/(knot_v.max()-knot_v.min())

# Calculate spline coefficients for x and y
fx=interpolate.make_interp_spline(x,x)
coeffx=fx.tck[1]
fy=interpolate.make_interp_spline(y,y)
coeffy=fy.tck[1]

# Create control point vectors
v=Base.Vector
ctrl=[]
for iZ,coeffsi in enumerate(coeffs):
    ctrl.append([])
    for jZ,_ in enumerate(coeffsi):
        ctrl[-1].append(v(coeffx[jZ],coeffy[iZ],coeffs[iZ,jZ]))

# Create spline surface in FreeCAD
periodic = False
bs=Part.BSplineSurface()
bs.buildFromPolesMultsKnots(ctrl, mult_u, mult_v, knot_u, knot_v, periodic, periodic, degree_u, degree_v)
s=bs.toShape()
Part.show(s)

# Show input points and control poins
for iZ, Zi in enumerate(Z):
    for jZ, _ in enumerate(Zi):
        # Input points
        Draft.make_point(v(X[iZ,jZ],Y[iZ,jZ],Z[iZ,jZ]))
        # Control points
        Draft.make_point(v(coeffx[jZ],coeffy[iZ],coeffs[iZ,jZ]))


#--Using input points as control points--

# Define control points
coeffs2=Z
coeffx2=X[0,:]
coeffy2=Y[:,0]

# Define knot vectors
knot_u2=np.hstack([np.zeros(degree_u),np.linspace(0,1,Z.shape[0]+1-degree_u),np.ones(degree_u)])
knot_v2=np.hstack([np.zeros(degree_v),np.linspace(0,1,Z.shape[1]+1-degree_v),np.ones(degree_v)])

# Convert knot vector format
knot_u2, mult_u2=np.unique(knot_u2,return_counts=True)
knot_v2, mult_v2=np.unique(knot_v2,return_counts=True)

# Normalize knot vectors
knot_u2=(knot_u2-knot_u2.min())/(knot_u2.max()-knot_u2.min())
knot_v2=(knot_v2-knot_v2.min())/(knot_v2.max()-knot_v2.min())

# Create control point vectors
v=Base.Vector
ctrl2=[]
for iZ,coeffsi in enumerate(coeffs2):
    ctrl2.append([])
    for jZ,_ in enumerate(coeffsi):
        ctrl2[-1].append(v(coeffx2[jZ],coeffy2[iZ],coeffs2[iZ,jZ]))

# Create spline surface in FreeCAD
periodic = False
bs2=Part.BSplineSurface()
bs2.buildFromPolesMultsKnots(ctrl2, mult_u2, mult_v2, knot_u2, knot_v2, periodic, periodic, degree_u, degree_v)
s2=bs2.toShape()
Part.show(s2)

doc.recompute()

[heda]

cool, and useful.

took the liberty to:

a) do some cosmetic changes to your script, and make the points with Points wb instead, makes the whole thing a bit snappier...

Code: Select all

import FreeCAD as App
import FreeCADGui as Gui
import Part, Points
import numpy as np
from scipy import interpolate

Vector = App.Vector

# Create a sample surface
x = np.linspace(-1, 1, 10) * np.pi
y = np.linspace(-1, 1, 10) * np.pi
X, Y = np.meshgrid(x, y)
Z = (np.sin(X) + 2*np.cos(2*Y) + Y) / 2

degree_u = degree_v = 3

doc = App.newDocument()

#--Forcing the surface to pass through all input points--

# Create spline surface
f = interpolate.RectBivariateSpline(x, y, Z.T)
coeffs = f.get_coeffs()
coeffs = np.reshape(coeffs, Z.T.shape).T
knots = f.get_knots()
# Convert knot vector format
knot_u, mult_u = np.unique(knots[1], return_counts=True)
knot_v, mult_v = np.unique(knots[0], return_counts=True)
# Normalize knot vectors
knot_u = (knot_u - knot_u.min()) / (knot_u.max() - knot_u.min())
knot_v = (knot_v - knot_v.min()) / (knot_v.max() - knot_v.min())

# Calculate spline coefficients for x and y
fx = interpolate.make_interp_spline(x, x)
coeffx = fx.tck[1]
fy = interpolate.make_interp_spline(y, y)
coeffy = fy.tck[1]

def RGB(*args):
    return tuple((float(i/255) for i in args))

def getVector(iz, jz, which='ctrl2'):
    if which == 'ctrl':
        cfx, cfy, cfz = coeffx, coeffy, coeffs
    elif which == 'spoints':
        cfx, cfy, cfz = X[iz], Y[iz], Z
    elif which == 'ctrl2':
        cfx, cfy, cfz = coeffx2, coeffy2, coeffs2
#    print(cfx[jz], cfy[iz], cfz[iz, jz])
    return Vector(cfx[jz], cfy[iz], cfz[iz, jz])

# Create control point vectors
ctrl = []
for iZ, coeffsi in enumerate(coeffs):
    ctrl.append([])
    for jZ, _ in enumerate(coeffsi):
        ctrl[-1].append(getVector(iZ, jZ, 'ctrl'))

# Create spline surface in FreeCAD
periodic = False
bs = Part.BSplineSurface()
bs.buildFromPolesMultsKnots(ctrl, mult_u, mult_v, knot_u, knot_v,
                            periodic, periodic, degree_u, degree_v)

s = Part.show(bs.toShape(), 'BSplineSurf')
s.Label2 = 'spline surface by points'
s.ViewObject.ShapeColor = RGB(0, 170, 255)
s.ViewObject.Transparency = 50

# Show input points and control poins
spts, cpts = [], []
for iZ, Zi in enumerate(Z):
    for jZ, _ in enumerate(Zi):
        # Input points
        spts.append(getVector(iZ, jZ, 'spoints'))
        # Control points
        cpts.append(getVector(iZ, jZ, 'ctrl'))

spoints = doc.addObject('Points::Feature')
spoints.Label2 = 'on surface points'
spoints.Points = Points.Points(spts)
spoints.ViewObject.ShapeColor = RGB(0, 170, 255)
spoints.ViewObject.PointSize = 5
cpoints = doc.addObject('Points::Feature')
cpoints.Label2 = 'control points'
cpoints.Points = Points.Points(cpts)
cpoints.ViewObject.ShapeColor = RGB(255, 170, 0)
cpoints.ViewObject.PointSize = 5


#--Using input points as control points--

# Define control points
coeffs2 = Z
coeffx2 = X[0, :]
coeffy2 = Y[:, 0]

# Define knot vectors
knot_u2 = np.hstack([np.zeros(degree_u),
                     np.linspace(0, 1, Z.shape[0] + 1 - degree_u),
                     np.ones(degree_u)])
knot_v2 = np.hstack([np.zeros(degree_v),
                     np.linspace(0, 1, Z.shape[1] + 1 - degree_v),
                     np.ones(degree_v)])

# Convert knot vector format
knot_u2, mult_u2 = np.unique(knot_u2, return_counts=True)
knot_v2, mult_v2 = np.unique(knot_v2, return_counts=True)

# Normalize knot vectors
knot_u2 = (knot_u2 - knot_u2.min()) / (knot_u2.max() - knot_u2.min())
knot_v2 = (knot_v2 - knot_v2.min()) / (knot_v2.max() - knot_v2.min())

# Create control point vectors
ctrl2 = []
for iZ, coeffsi in enumerate(coeffs2):
    ctrl2.append([])
    for jZ, _ in enumerate(coeffsi):
        ctrl2[-1].append(getVector(iZ, jZ, 'ctrl2'))

# Create spline surface in FreeCAD
periodic = False
bs2 = Part.BSplineSurface()
bs2.buildFromPolesMultsKnots(ctrl2, mult_u2, mult_v2, knot_u2, knot_v2,
                             periodic, periodic, degree_u, degree_v)

s2 = Part.show(bs2.toShape(), 'BSplineSurf')
s2.Label2 = 'spline surface by control polygon'
s2.ViewObject.ShapeColor = RGB(255, 170, 0)
s2.ViewObject.Transparency = 50


doc.recompute()

av = Gui.ActiveDocument.ActiveView
av.viewIsometric()
av.fitAll()

b) make Macro_BSurf_from_grid

[snow54]

Thank you for cleaning up my sample code. However, your code creates an error below.

<class 'AttributeError'>: 'NoneType' object has no attribute 'Label2'

Should you use the line below

Code: Select all

s = doc.addObject("Part::Feature", 'BSplineSurf')

instead of the one below?

Code: Select all

s = Part.show(bs.toShape(), 'BSplineSurf')

It applies to s2 as well. After replacing them, it works well.

[onekk]

What version if FC are you using?

Please put the info requested in the topic IMPORTANT: at top of the forum page.

Probably @heda is writing his code using a recent v0.20, as the way Part.show() is working has been modified, and now is returning an instance of the created "DocumentObject".

Direct creation of "DocumentObject" using the code you posted is somewhat the "old way" and the error seems to confirm that.

Regards

Carlo D.

[snow54]

Thank you for pointing that out. I was still using v0.19. That explains.

[heda]

indeed, it runs as is on 0.20 but not 0.19 - was not aware of that difference but good to know for the future