Subdivision Curves¶
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
Subdvision methods progressively refine a discrete curve and converge to a smooth curve. This allows to perform an interpolation or approximation of a given coarse dataset.
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('toolbox_wavelet_meshes')
addpath('solutions/meshwav_1_subdivision_curves')
Curve Subdivision¶
Starting from an initial set of control points (which can be seen as a coarse curve), subdivision produces a smooth 2-D curve.
Shortcut to plot a periodic curve.
ms = 20; lw = 1.5;
myplot = @(f,c)plot(f([1:end 1]), c, 'LineWidth', lw, 'MarkerSize', ms);
myaxis = @(rho)axis([-rho 1+rho -rho 1+rho], 'off');
We represent a dicretized curve of $N$ points as a vector of complex numbers $f \in \CC^N$. Since we consider periodic boundary conditions, we assume the vectors have periodic boundary conditions.
Define the initial coarse set of control points $f_0 \in \CC^{N_0}$.
f0 = [0.11 0.18 0.26 0.36 0.59 0.64 0.80 0.89 0.58 0.22 0.18 0.30 0.58 0.43 0.42]' + ...
1i * [0.91 0.55 0.91 0.58 0.78 0.51 0.81 0.56 0.10 0.16 0.35 0.42 0.40 0.24 0.31]';
Rescale it to fit in a box.
f0 = rescale(real(f0),.01,.99) + 1i * rescale(imag(f0),.01,.99);
Display it.
clf; myplot(f0, 'k.-');
myaxis(0);
One subdivision step reads $$ f_{j+1} = (f_j \uparrow 2) \star h. $$ This produces discrete curves $f_j \in \CC^{N_j}$ where $N_j = N_0 2^j$.
Here $\uparrow 2$ is the up-sampling operator $$ (f \uparrow 2)_{2i}=f_i \qandq (f \uparrow 2)_{2i+1} = 0. $$
Recall that the periodic discrete convolution is defined as $$ (f \star h)_i = \sum_j f_j h_{i-j}, $$ where the filter $h$ is zero-padded to reach the same size as $f$.
The low pass filter (subdivision kernel) $h \in \CC^K$ should satisfies $$ \sum_i h_i = 2 . $$ This ensure that the center of gravity of the curve stays constant $$ \frac{1}{N_j} \sum_{i=1}^{N_j} f_{j,i} = \frac{1}{N_0} \sum_{i=1}^{N_0} f_{0,i}.$$
Define the subdivision operator that maps $f_j$ to $f_{j+1}$.
subdivide = @(f,h)cconvol( upsampling(f), h);
We use here the kernel $$ h = \frac{1}{8}[1, 4, 6, 4, 1]. $$ It produced a cubic B-spline interpolation.
h = [1 4 6 4 1];
h = 2*h/sum(h(:));
Initilize the subdivision with $f_0$ at scale $j=0$.
f = f0;
Perform one step.
f = subdivide(f,h);
Display the original and filtered discrete curves.
clf; hold on;
myplot(f, 'k.-');
myplot(f0, 'r.--');
myaxis(0);
Exercise 1
Perform several step of subdivision. Display the different curves.
exo1()
%% Insert your code here.
Under some restriction on the kernel $h$, one can show that these discrete curves converges (e.g. in Hausdorff distance) toward a smooth limit curve $f^\star : [0,1] \rightarrow \CC$.
We do not details here sufficient condition to ensure convergence and smoothness of the limitting curve. The interested reader can have a look at DynLevin02 for a review of theoritical guarantees for subdivision.
The limit curve $f^\star$ is a weighted average of the initial points $f_0 = (f_{0,i})_{i=0}^{N_0-1} \in \CC^{N_0}$ using a continuous scaling function $\phi : [0,1] \rightarrow \RR$ $$ f^\star(t) = \sum_{i=0}^{N_0-1} f_{0,i} \phi(t-i/N_0). $$ The continuous kernel $\phi$ is a low-pass function which as a compact support of width $K/N_0$. The control point $f_{0,i}$ thus only influences the final curve $f^\star$ around $t=i/N_0$.
The scaling function $\phi$ can be computed as the limit of the sub-division process $f_j$ when starting from $f_0 = \delta = [1,0,\ldots,0]$, which is the Dirac vector.
Exercise 2
Compute the scaling function $\phi$ associated to the subdivision.
exo2()
%% Insert your code here.
Exercise 3
Test with different configurations of control points.
exo3()
%% Insert your code here.
hcc = @(w)[1 w w 1]/(1+w);
For $w=3$, the scaling function $\phi$ is a quadratic B-spline.
Exercise 4
Test the corner-cutting for $w=3$.
exo4()
%% Insert your code here.
Exercise 5
Test the corner-cutting for vaious values of $w$.
exo5()