BSplineCurve interpolation currently assumes uniform B-splines?

[mandyxmq]

Hi, I was wondering whether the current BSplineCurve evaluation

mitsuba3/src/shapes/bsplinecurve.cpp

Line 1122 in 59d980a

cubic_interpolation(const Float v, const UInt32 prim_idx, Mask active) const {

assumes the control points are uniformly spaced (same as in the formula in the attached image). I tested using simple 1D examples that this evaluation is no longer accurate when we have repeated control points at the end (for clamped B-splines) or non-uniform spacing. I was wondering if a more general evaluation is desired in this case as control points are not guaranteed to be uniformly spaced. Thanks!

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BSpline

def cubic_interpolation(v, points):
    v2 = v * v
    v3 = v2 * v
    multiplier = 1.0 / 6.0
    c = (1-v) * (1-v) * (1-v) * points[0] + \
            (3.0 * v3 - 6.0 * v2 + 4.0) * points[1] + \
            (-3.0 * v3 + 3.0 * v2 + 3.0 * v + 1.0) * points[2] + \
            v3 * points[3]

    return c * multiplier

# uniform case
k = 3
t = [0, 1, 2, 3, 4, 5, 6, 7, 8]
basis_functions = [BSpline(t, (np.arange(len(t)-1) == i).astype(int), k) for i in range(len(t)-k-1)]

x = np.linspace(4, 5, 50)
basis_values = np.array([b(x) for b in basis_functions])

for i, b in enumerate(basis_values):
    plt.plot(x, b, label=f'B-spline {i+1}')
plt.legend()
plt.grid()

# clamped B-spline
k = 3
t = [0, 0, 0, 0, 1, 2, 2, 2, 2]
basis_functions = [BSpline(t, (np.arange(len(t)-1) == i).astype(int), k) for i in range(len(t)-k-1)]

x = np.linspace(0, 1, 50)
basis_values = np.array([b(x) for b in basis_functions])

for i, b in enumerate(basis_values):
    plt.plot(x, b, label=f'B-spline {i+1}')
plt.legend()
plt.grid()

# another non-uniform case
k = 3
t = [0, 1, 2, 3, 4, 5, 8, 9, 10]
basis_functions = [BSpline(t, (np.arange(len(t)-1) == i).astype(int), k) for i in range(len(t)-k-1)]

x = np.linspace(4, 5, 50)
basis_values = np.array([b(x) for b in basis_functions])

for i, b in enumerate(basis_values):
    plt.plot(x, b, label=f'B-spline {i+1}')
plt.legend()
plt.grid()

[ziyi-zhang]

Hi Mandy,

I vaguely recall that we are using the clamped formulation. But there is a conversion between v_local and v_global. For instance, if we want an endpoint p0 as the curve endpoint, (I think, haven't checked) we need to put control points like (2p0-p1, p0, p1, p2, ...).

I am not quite following this. What is the expected behavior and what is our current result?

this evaluation is no longer accurate when we have repeated control points at the end

[mandyxmq]

Hi Ziyi,

I think for enforcing the curve to go through a particular endpoint, we can use a clamped B-spline where we repeat the endpoint 4 times (including itself), so when interpolated, that point is guaranteed to be on the curve.

To evaluate a non-uniform B-spline (including the clamped case because towards the end the knots are not uniformly spaced), I believe we need to use Cox-de Boor recurrence to compute things https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.html.

I have compared Mitsuba's cubic_interpolation and Scipy's B-spline evaluations and they give slightly different results when we have clamped B-splines.

For simple illustrations, the code I attached in the first post shows how the basis functions look like on one segment. You can see that only in the first case they look the same as the image containing the formula Mitsuba uses, because there we have uniformly spaced knots. However in the other two cases the basis functions look different, thus using the same formula would result in inaccurate B-spline evaluation I think.

[wjakob]

Hi Mandy,

High level comment: I think that Mitsuba merely exposes what is (consistently) implemented by both Embree and OptiX. So there is not much freedom to change the definition, but perhaps we can better document the specifics. Could you take a look at the Bspline curve primitive in OptiX and Embree to see if this identifies the reason for why the resulting curves do not match your expectation?

[wjakob]

Hi Mandy,

High level comment: I think that Mitsuba merely exposes what is (consistently) implemented by both Embree and OptiX. So there is not much freedom to change the definition, but perhaps we can better document the specifics. Could you take a look at the Bspline curve primitive in OptiX and Embree to see if this identifies the reason for why the resulting curves do not match your expectation?

[ziyi-zhang]

Hi Mandy,
We are not following the scipy interpolation conventions since the preliminary ray-curve intersection is computed by Embree/Optix as Wenzel mentioned. Here is the optix documentation that explains the end point behavior (Sec 8.5).

We do not have a Mitsuba implementation of ray-curve intersection (for meshes, we still have the KD-tree if we don't use Optix/Embree). The related code you quoted about cubic interpolation is only to compute derived quantities (uv, normal, differentials) from the preliminary result Optix returns (v_global).

If we want the behavior to match the scipy one, an easy workaround would be to write a conversion script so Optix sees processed control points that follows its convention.

[mandyxmq]

Hi Ziyi and Wenzel,

Thanks for the pointer to the Optix documentation regarding their curve representation. They mentioned that their curve representation is defined by cubic uniform B-spline curve (or a quadratic uniform B-spline curve, a Catmull-Rom spline curve, a Bézier curve, or a series of linear segments). Therefore the control points need to be uniformly spaced in the cubic B-spline case, unless they process the control points input to Optix (I couldn't find that information in their guide). Similarly, in the cubic_interpolation function in Mitsuba, the B-spline curve is assumed to have uniformly spaced control points. This has nothing to do with Scipy's convention -- this means that we model a special kind of B-spline curves, instead of the more general ones. This is fine, as long as the control points given to Mitsuba are uniformly spaced, otherwise the inconsistency happens. I think maybe we could add that specification in the documentation, or more robustly handle this by preprocessing the control points before rendering.