# Fit bezier curve to mesh

Hi,
is there a oprtion via python to fit a curve (bezier) to mesh?
Blender has already a cruve fit algorithm (i guess - curve decimate) but is this algorithm accessable via python ?.

What do you mean by āfit a curve to meshā? Please explain the before/after that you want.

before:
There are a set of points from a mesh

after:
A Bezier segment where the curve fits the point position.

There is a function to get the interpolated bezier points

`````` mathutils.geometry.interpolate_bezier(knot1, handle1, handle2, knot2, resolution)
``````

but i need the reverse behaviour wehere you have the start and endpoint and also the interpolated points and get the handle vectors.

I hope that make sense.

1 Like

Thanks, that helps me answer you. Well, there is a library that you can easily build and import that was written by a Blender developer for, IIRC, internal use with the GP. It can import into Python. I havenāt tried it personally but Iāve been meaning to. I donāt know how portable it is, either. Iāve built it and imported it but I didnāt try to use it in my addon, because there was an easier, less precise way to do what I needed.

Within Blender, I donāt know of any curve fit function thatās exposed to the user, but I think internally it uses the library I linked for everything.

The easier, less precise way: Just make a Bezier curve from the points with all the control points set to āAutoā. Thereās an operator for doing this already, and then I suppose you can use the curve-decimate addon (I donāt know how that one works).

1 Like

Iāve found that this isnāt even a very helpful function. In my addon I convert the curve to a mesh and extract the data from it. Itās fast, reliable, and completely matches what the user expects. And itās simple, too. Let Blender do the complex math!

Oh, also, this may be worth looking at: A Primer on BĆ©zier Curves (thanks to @sybren for showing this to me )

2 Likes

Thanks for all your help

ok i tried the decimate option but the result is not so presice like in the curve fit library.
can you explain how i can use this library , do i have to compile blender?

No you donāt have to compile Blender, just the library. I think itās as simple as cloning the repo and `make` `make install`. But I havenāt done anything beyond building it, myself.

@floka
Perhaps the experimental script at the end of this post can help you.
It employs `numpy.linalg.lstsq` to fit Bezier spline curves to stand-alone mesh edge loops.
The standalone mesh edge loops should themselves be products of conversions of curves to meshes.
I imagine that in a workflow, one would start with a curve and shape it to oneās needs, perhaps using a shrink modifier to profile the surface of another object. I suppose, while in edit mode, one could also copy out an edge loop of an existing mesh (āShift-Dā), turn it into an independent mesh object (āPā) then convert that mesh object into a curve. See Howto. For the terminally curious, Rainy Day Reading discusses the scriptās design and limitations.

Caveat: This is not production code. It illustrates a use of `numpy` in Blender. Take care.

### TL;DR

Dataset: A mesh generated from a curve (Object -> Convert To -> Mesh) Previously, the shrink-wrap modifier plastered the curve onto another mesh, blue in the illustration.

Fitted Bezier spline: Estimated using `numpy.linalg.lstsq` as a Least Squares solver.

Comparison: The estimated Bezier curve overlaying the mesh.

Misfits: Arise if there are irregularities in the mesh.
Howto:

1. Copy the script; paste it into a Blender text buffer.
2. Use any kind of Blender curve. Call it āProbeā for discussion.
• Shape Probe or capture a shape with it - shrinkwraping Probe onto a surface, for example.
1. Convert Probe to a mesh - a free-standing edge loop.
• Caution: Donāt start with a mesh, bypassing Probe. See FIXME_1: in the script.
• If you DO copy out an edge loop from an existing mesh (Edit mode: Shift-D, P), creating a new mesh object, go to object mode, convert the new mesh object temporarily to a curve (Object -> Convert To -> Curve) then back to a mesh (Object -> Convert To -> Mesh).
1. Select the mesh-converted Probe; it should be the active object.
2. Run the script in the text buffer.
3. A fitted Bezier spline, āFitted_Bezier,ā curve appears. It should closely track Probe.
• Adjust TOL in the script, making it bigger for more loosely fitting curves.
1. Use āFitted_Bezierā in your further modeling schemes.
2. Come some idle day, peruse Rainy Day Reading.

### Rainy Day Reading

At itās heart this script:

(1) makes explicit what is implicit in a very large number of āeducatedā guesses we make, embodied in a matrix A - āeducatedā in so far as we have a theory - a statistical model - with a few choice parameters x. This theory predicts how a large quantity of plot points b occur in the locations weāve been given.

(2) tells us the liklihood of whether our theory is viable or a busted flush. If likely, then we have a solution with the x of Ax = b. That is, these x are very likely to be the correct control points for fitting a Bezier spline to the observed plot points constituting b.

What we donāt know: The explicit x - aka, our control points.

What we do know: An awful lot of locations b that - by gum! - look like a Bezier spline!

Our theory: A spline of cubic Bezier curves, each with a few choice parameters - their control points - can plot all the observed positions we have been given. Using just (a few) choice control points, we can predict other positions that could also be among the points we have been given.

Constructing our theory: For each spline segment, a cubic Bezier polynomial pairs a control point from x with a bezier basis function written in terms of an independent parameter u. Given u, the basis function returns a weight in the range [0ā¦1] that weights the control point, establishing its contribution to a computed location. These basis functions are linked together in the sense that all the weights they compute add up to 100% - unity. Such things are called affine combinations. All the parts - contributions - make up a whole - the computed location in space. No more. No less. Otherwise the weights do not make an affine combination and the shape of the curve changes relative to its position with respect to the origin. Cubic polynomials have four terms, four pairings of control points with their allied basis functions. This weighted sum of control point locations constitutes one plotted point.

Relating the known to the unknown: A is a set of educated guesses produced by our theory. Each row of A connects a plot point from b to our theory through a set of weights, corresponding to each purported control point. Some people call A an observation matrix because of the way it connects observational reality (plot points) to theoretical fantasy (educated guesses). These weights set the contribution of purported control points drawn from x to the location of a specific plot point drawn from b. We donāt know where those purported control points are, but we provide hints in the form of implied behavior: wherever those control points are, when we weight their contributions - using these particular weights! - we compute the location of a corresponding plot point in b. Behind the scenes, weāve used the Bezier basis functions by themselves to compute the weights. All of our myriad hints, one hint for each plot point, taken together imply a specific set of locations for the control points, or so we claim. For these hints are not grounded in any kind of reality. They are grounded in our theory, to wit: that a Bezier spline curve underlies this large list of locations. We are doing statistics here, not mathematics, and Bezier curves constitutes our statistical model.

Technique: Having composed our myriad educated guesses in A we ask a Least Squares Solver in numpy - `numpy.linalg.lstsq` to invert A and put on the right hand side of the equation, solving for the control points x. Iām not going into the LS technique here. I recommend Sec 15, Numerical Recipes in C, āLeast Squares as a Maximum Likelihood Estimatorā (Cambridge University Press, 2nd Ed. pp 657ff, Press, Teukolsky, Vetterling & Flannery). About a bazillion PDFās of this reference are floating around the Web. While your nose is in that reference, also take a peek at Singular Value Decomposition (SVD), section 2.6, a technique to separate matrices like A into three, more computationally tractible matrices.

Bookkeeping: Most of it involves cooking up the guesswork matrix: A.

`tryit()` contains the development run-up to A culminating in the request to invert A and provide (1) likey coefficients (control points): x and (2) an array of residues in each coordinate of the plot point. These tell us if we have a busted flush in any one or more dimensions. The larger the residue, the more likely we have a busted theory.

1. `plotpoints` is a (very long) vector of plot points.

2. `diff` is a (one less long) vector of difference from one plot point to the next.

3. `slen` are those differences expressed as lengths, found by way of the Pythagorean Theorem.

4. `ssum` are the relative distances of each plot point from the origin of the purported curve. The first plot gets zero; the last plot gets one; the interior plots get a relative value, these derived from sums of lengths along the curve from the origin to the plot at hand (arc length), normalized by the entire length of the curve. These become trial u values, one proposed for each plot point in `plotpoints`.

5. `obmx` The observation matrix, aka A. There is a bit of numpy `ufunc` magic here. For example, `bf0(ssum).reshape(ascol)` applies the first basis function computation, elementwise, to ssum, producing a new array which is rotated into a column vector by the numpy array method `reshape()`. The four column vectors are then horizontally stacked into the observation matrix.

Main loop: Divide and conquer by binary division. We manage this through an array of slices kept in `stack.` Typical game: First try for the whole curve as one Bezier segment. If that theory exhibits high residue, then divide the plot point set by two and attempt to find fits for the two halves. And so on. When the stack empties, we have a sequence of solutions: in aggregate, the purported Bezier spline.

This is an illustrative script, not a production one. Script comments discuss a few weaknesses. I would be interested if people find areas where it works and where it does not, perhaps flushing out other weaknesses. Thank you in advance for that.

``````#! /usr/bin/python
"""Fit a Bezier spline curve to a standalone edge loop generated from
a curve. For use with Blender 3D Python API. This code is for
illustrative purposes only. It is not production-ready. Available
under Creative Commons Attribution-NonCommercial-ShareAlike 3.0

import numpy            as np
import bpy
from bpy import data    as D
from bpy import context as C
from mathutils import Matrix, Vector

# In a production environment, these data should
# come from a user interface.
# Tolerance. Larger numbers give looser fits.
# TOL = 0.1    # This gives loose solutions - good for noisy data
# TOL = 0.01
# TOL = 0.001
# TOL = 0.0001
TOL = 0.00001  # This gives tight solutions - good for precision
RES = 8        # Sets mesh density for Blender curves.
NME = 'Fitted_Bezier' # Name of the fitted curve -
MINPTS = 4 # Don't fit slices smaller than four points.

# Lambdas for generating an observation matrix. These are
# Bezier basis functions which bakes in the idea that the
# observation data points lie along Bezier spline curves.
# Change these for other bases.
# ENHANCEMENT: Use Carl DeBoor's Algorithm instead of hard-
# wired basis functions.

bf0 = lambda u: u**3
bf1 = lambda u: 3*((u**2)*(1-u))
bf2 = lambda u: 3*(u*((1-u)**2))
bf3 = lambda u: (1-u)**3

def splitslice(stackslice):
"""From splice s return a list of two splices
0: from beginning to mid-element,
1: from mid-element to end"""

interval = (stackslice.stop - stackslice.start) >> 1
remainder = (stackslice.stop - stackslice.start) %  2
return [ \
slice(stackslice.start, stackslice.start + interval + remainder, 1), \
slice(stackslice.start + interval, stackslice.stop, 1)              \
]

def tryfit(plotpoints):
"""plotpoints: numpy array of 3D space points. shape (N, 3)
return: python list of numpy arrays, shape (4 x 3)
Cubic Bezier spline control points."""

dln = len(plotpoints)
ascol = (dln, 1)
diff = plotpoints[1:] - plotpoints[:-1]
slen = np.array([np.sqrt((p**2.0).sum()) for p in diff])
ssum = np.array([0.0] + [slen[:k].sum() for k in range(1, len(slen)+1)])/slen.sum()
obmx = np.hstack(                               \
(                             \
bf0(ssum).reshape(ascol), \
bf1(ssum).reshape(ascol), \
bf2(ssum).reshape(ascol), \
bf3(ssum).reshape(ascol)  \
)                             \
)
return np.linalg.lstsq(obmx, plotpoints, -1)

# Try to fit the whole dataset to a one-segment spline.
# If that fails, split and try halves. Recurse. Maintain a
# stack of to-do slices until the topmost slice on the
# stack fits. Pop it, save the fitted segment, and
# commence on the next slice on the stack. Recurse
# until the stack is empty.

aod = C.active_object.data

# FIXME_1: Generally, vertices read from a mesh data block aren't
# in any kind of order. vb assignment exploits the fact that an
# edge-loop generated from a curve 'happens-to-be' in geometric
# order. Replace this assignment below with an edge-loop-walker
# that reliably fetches vertices in geometric order - perhaps
# from a user-selected edge loop.

vb = np.array([v.co for v in aod.vertices])
stack = [slice(0, len(vb), 1)] # stack with one slice
sset = list() # Save set of fitted bezier curves

while bool(stack):
# Residue furnishes per-dimension quality of fit
# Larger residues indicate poorer fits. Not used -
# rank: If less than 4, results not reliable
# singular: As these approach zero, the corresponding
#           control point contributes less to the solution.
#           Small singulars suggest removing the
#           corresponding control point from the
#           solution.
# tryfit should be in a try - except block to
# catch fit failures.
solution, residue, rank, singular = tryfit(vb[stack[-1]])
ltol = np.sqrt((residue**2.0).sum())
if (ltol < TOL) or (len(vb[stack[-1]]) < MINPTS):
# Good fit - or slice too thin. Save
# solution and pop the slice.
stack.pop()
sset.insert(0, solution)
else:
# Too much residue. Pop top slice,
# divide it at mid-plotpoint and
# make a list of left and right
# slices.
lslice = stack.pop()
stack += splitslice(lslice)

# Construct Bezier curve data block from sset
DLEN = len(sset)
BINDX = 0
crv = D.curves.new(NME + '.Curve', type='CURVE')
crv.dimensions = '3D'
bseg = crv.splines.new('BEZIER')
bseg.resolution_u = RES
while bool(sset):
dp = sset.pop()
if not BINDX:
# At start of Bezier curve spline
bseg.bezier_points[BINDX].handle_left_type = 'FREE'
bseg.bezier_points[BINDX].handle_left = Vector(dp[0])
bseg.bezier_points[BINDX].co = Vector(dp[0])
bseg.bezier_points[BINDX].handle_right = Vector(dp[1])
bseg.bezier_points[BINDX].handle_right_type = 'ALIGNED'
bseg.bezier_points[BINDX + 1].handle_left = Vector(dp[2])
bseg.bezier_points[BINDX + 1].handle_left_type = 'ALIGNED'
elif bool(sset):
# In the midst of the Bezier curve spline
bseg.bezier_points[BINDX].co = Vector(dp[0])
bseg.bezier_points[BINDX].handle_right_type = 'ALIGNED'
bseg.bezier_points[BINDX].handle_right = Vector(dp[1])
bseg.bezier_points[BINDX + 1].handle_left = Vector(dp[2])
bseg.bezier_points[BINDX + 1].handle_left_type = 'ALIGNED'
else:
# At the end of the Bezier curve spline
bseg.bezier_points[BINDX].co = Vector(dp[0])
bseg.bezier_points[BINDX].handle_right = Vector(dp[1])
bseg.bezier_points[BINDX].handle_right_type = 'ALIGNED'
bseg.bezier_points[BINDX + 1].handle_left = Vector(dp[2])
bseg.bezier_points[BINDX + 1].handle_left_type = 'ALIGNED'
bseg.bezier_points[BINDX + 1].co = Vector(dp[3])
bseg.bezier_points[BINDX + 1].handle_right = Vector(dp[3])
bseg.bezier_points[BINDX + 1].handle_right_type = 'FREE'

BINDX += 1

# Construct Object data block. Link to Curve. Add to view layer
bobj = bpy.data.objects.new(NME, crv)
bobj.matrix_world = Matrix(bobj.matrix_world @ C.active_object.matrix_world)