Multi-resolution encoding with basis functions

In multibezier.py, an algorithm (adaptive_fit) is proposed to encode a (possibly multi-dimensional) function as a superposition of local Bernstein polynomials (which form the basis of Bézier curves) of increasingly higher resolution with recursive subdivision of the sample space.

This repository also includes demo_SDF01.py to generate signed distance functions (used here as test signals), and vis_sdf.py to visualize 3D SDFs, along with examples SDFs in .npy files.

Usage

Psi, weights = adaptive_fit(x, in_shape, nbOut, nbFct, iterations=1, threshold=0.0, C1=True)

where

• x is an array that contains the signal to fit
• in_shape is a tuple that contains the input dimensions
• nbOutis the number of output dimensions
• nbFctspecifies the degree + 1 of the polynomials used for the encoding (4 for cubic, 3 for quadratic)
• iterations is the maximum number of recursive subdivisions to perform
• threshold is the maximum root mean square error to allow in a subdivision
• C1 : if true, the reconstruction will be $C^1$ continuous. If false, it will be at least $C^0$ continuous

Tests using signed distance functions (SDFs) can be found at the bottom of multibezier.py as usage examples. These can be run with

python multibezier.py

assuming NumPy and matplotlib are installed.

Algorithm

A dictionary of low resolution basis polynomials $\Psi$ is first used to compte a coarse approximation of the signal by means of least squares estimation :

$$\hat{x} = \Psi w$$

with $$w = \Psi^\dagger x$$ The signal is split into $2^n$ regions, where $n$ is the number of dimensions.

Then, at each iteration and for each subdivision, the residual error is evaluated:

$$x = e - \Psi w$$

If the mean square error of the flattened region is bigger than some threshold, then the basis is extended to add more details in this region, i.e. the $\Psi$ matrix is augmented with additional columns that represent the detail basis functions: $$\Psi \gets \begin{bmatrix} \Psi & \Psi_{details} \end{bmatrix}$$ and new weights are re-computed based on the whole signal; effectively adding a correction to the approximation: $$ \hat{x}{i+1} = \begin{bmatrix} \Psi & \Psi{details} \end{bmatrix} \begin{bmatrix} \mathbf{w} \ \mathbf{w_{details}} \end{bmatrix} = \hat{x}{i} + \hat{x}{details} $$

The data structure resulting from the recursive subdivision is a binary tree in the case of a one-dimensional signal, a quadtree for a function with two input variables and an octree in the three variables case.

To preserve $C^1$ continuity when superposing Bézier curves, it is necessary to set the two first and two last control points of the additional curves to zero That way, the detail curves are null and flat at extremities, preserving $C^1$ continuity between subdivisions. This also requires the use of sub-segments (at least 2 for cubic curves and at least 3 for quadratic curves).

Limitations

• Not guaranteed to work for splines of higher order (only tested for quadratic and cubic)
• Computational complexity grows quickly in the multi-dimensional case