Subdivision Curves

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.

Contents

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox, general toolbox, graph toolbox and wavelet_meshes toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal, toolbox_general, toolbox_graph and toolbox_wavelet_meshes in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');
getd('toolbox_graph/');
getd('toolbox_wavelet_meshes/');

Curve Subdivision

Starting from an initial set of control points (which can be seen as a coarse curve), subdivision produces a smooth 2D curve.

In the following, a 2D curve is discretized as a vector of complex numbers.

Define the low pass filter (subdivision kernel). It should sum to 2.

h = [1 4 6 4 1];
h = 2*h/sum(h(:));

Define the initial coarse set of control points.

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);

Initilize the subdivision.

f = f0;

One step of subdivision correspond to a step of inverse wavelet transform (applied to each X/Y channel of f0) but with vanishing wavelet coefficients.

First upsample each channel, which corresponds to inserting zeros. The filter using a circular convolution.

f = cconv( upsampling(f), h);

Display the original and filtered discrete curves.

ms = 20; lw = 1.5;
hold on;
hh = plot(f([1:end 1]), 'k.-');
set(hh, 'MarkerSize', ms); set(hh, 'LineWidth', lw);
hh = plot(f0([1:end 1]), 'r.--');
set(hh, 'LineWidth', lw);
axis([0 1 0 1]); axis off;

Exercice 1: (the solution is exo1.m) Perform several step of subdivision. Display the different curves. 

exo1;

Exercice 2: (the solution is exo2.m) Compute the scaling function (interpolation kernel) associated to the subdivision. Here this is a cubic B-spline interpolation. 

exo2;

Exercice 3: (the solution is exo3.m) Test with different configuration of control points. 

exo3;

Exercice 4: (the solution is exo4.m) Try with different filters h. Here we use h=[1 3 3 1]/4, which corresponds to the famous Chaikin "corner cutting" scheme. This corresponds to a quadratic B-spline interpolation. 

exo4;

Interpolating Subdivision

Interpolating schemes keeps unchange the set of point at the previous level, and only smooth the position of the added points.

It means that the filter h should satisfies h(0)=1 and h(2k)=0 for all non zero k.

This is a filter proposed by Dyn, Levin and Gregory.

h = [-1, 0, 9, 1, 9, 0, -1]/16;
h((end+1)/2)=1;

Exercice 5: (the solution is exo5.m) Perform the interpolating subdivision. 

exo5;

Exercice 6: (the solution is exo6.m) Compare the result of the quadratic B-spline, cubic B-spline, and Dyn Levin Gregory. 

exo6;

Curve Approximation

Given an input, complicated curve, it is possible to approximate is by sampling the curve, and then

Load an initial random curve.

options.bound = 'per';
n = 1024*2;
sigma = n/8;
F = perform_blurring(randn(n,1),sigma,options) + 1i*perform_blurring(randn(n,1),sigma,options);
F = rescale(real(F),.01,.99) + 1i * rescale(imag(F),.01,.99);

Display it.

clf;
hh = plot(F([1:end end]), 'k');
set(hh, 'LineWidth', lw);
axis([0 1 0 1]);

Load an interpolating subvision mask.

h = [-1, 0, 9, 1, 9, 0, -1]/16;
h((end+1)/2)=1;

Exercice 7: (the solution is exo7.m) Perform an approximation f of the curve using a uniform sampling with 20 points. 

exo7;

To quantify the quality of the approximation, we use the L2 Hausdorff distance.

Compute the pairwise distances between points.

D = abs( repmat(F, [1 length(f)]) - repmat(transpose(f), [length(F) 1]) );

Compute the Hausdorff distance.

d = sqrt( ( mean(min(D).^2) + mean(min(D').^2) )/2 );
disp(['Hausdorff distance: ' num2str(d, 3)]);

Hausdorff distance: 0.00943

Exercice 8: (the solution is exo8.m) Display the decay of the Hausdorff approximation error as the number of sampling points increases. 

exo8;

Exercice 9: (the solution is exo9.m) Perform an optimization of the position of the points to minimize the Hausdorff approximation error. 

exo9;

Exercice 10: (the solution is exo10.m) Test on other curves. 

exo10;

3D Curve Subdivision

It is possible to construct 3D curves by subdivision.

Exercice 11: (the solution is exo11.m) Perform curve subdivision in 3D space. 

exo11;


Copyright ® 2008 Gabriel Peyre