silgon - 1 month ago 8x

Python Question

I'm trying to do something like the following (image extracted from wikipedia)

`#!/usr/bin/env python`

from scipy import interpolate

import numpy as np

import matplotlib.pyplot as plt

# sampling

x = np.linspace(0, 10, 10)

y = np.sin(x)

# spline trough all the sampled points

tck = interpolate.splrep(x, y)

x2 = np.linspace(0, 10, 200)

y2 = interpolate.splev(x2, tck)

# spline with all the middle points as knots (not working yet)

# knots = x[1:-1] # it should be something like this

knots = np.array([x[1]]) # not working with above line and just seeing what this line does

weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))

tck = interpolate.splrep(x, y, t=knots, w=weights)

x3 = np.linspace(0, 10, 200)

y3 = interpolate.splev(x2, tck)

# plot

plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')

plt.show()

The first part of the code is the code extracted from the main reference but it's not explained how to use the points as control knots.

The result of this code is the following image.

The points are the samples, the blue line is the spline taking into account all the points. And the red line is the one that is not working for me. I'm trying to take into account

`knots=x[1:-1]`

Note: this last image is exactly what I need, and it's the difference between what I have (spline passing all the points) and what I need (spline with control knots). Any ideas?

Answer

I just found something really interesting with the answer that I need with a bézier in this link. Then I used the code to try on my own. It's working fine apparently. This is my implementation:

```
#! /usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import binom
def Bernstein(n, k):
"""Bernstein polynomial.
"""
coeff = binom(n, k)
def _bpoly(x):
return coeff * x ** k * (1 - x) ** (n - k)
return _bpoly
def Bezier(points, num=200):
"""Build Bézier curve from points.
"""
N = len(points)
t = np.linspace(0, 1, num=num)
curve = np.zeros((num, 2))
for ii in range(N):
curve += np.outer(Bernstein(N - 1, ii)(t), points[ii])
return curve
xp = np.array([2,3,4,5])
yp = np.array([2,1,4,0])
x, y = Bezier(list(zip(xp, yp))).T
plt.plot(x,y)
plt.plot(xp,yp,"ro")
plt.plot(xp,yp,"b--")
plt.show()
```

And an image for the example.

The red points represent the control points. That's it =)

Source (Stackoverflow)

Comments