Subdivision Surfaces

$\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 mesh and converge to a smooth surface. This allows to perform an interpolation or approximation of a given coarse dataset.

In [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('toolbox_wavelet_meshes')
addpath('solutions/meshwav_2_subdivision_surfaces')

Subdivision of a Regular Polyedra

Starting from a control mesh which is a regular polyhedra, one can construct a sequence of mesh that converge to a sphere by subdividing each edge into two edges, and each triangle into four smaller triangles. The position of the mid points are projected onto the sphere.

Compute two examples of initial base mesh.

In [3]:
[vertex1,face1] = compute_base_mesh('oct');
[vertex0,face0] = compute_base_mesh('ico');

Display it.

In [4]:
clf;
subplot(1,2,1);
plot_mesh(vertex1,face1); 
shading('faceted'); lighting('flat'); view(3); axis('tight');
subplot(1,2,2);
plot_mesh(vertex0,face0); 
shading('faceted'); lighting('flat'); view(3); axis('tight');

Initialize the subdivision.

In [5]:
face = face0;
vertex = vertex0;

Compute the set of edges.

In [6]:
edge = compute_edges(face);

Number of vertex and edges.

In [7]:
n = size(vertex,2); 
ne = size(edge,2);

Compute the number of the three edges associated to each face.

In [8]:
A = sparse([edge(1,:);edge(2,:)],[edge(2,:);edge(1,:)],[n+(1:ne);n+(1:ne)],n,n);
v12 = full( A( face(1,:) + (face(2,:)-1)*n ) );
v23 = full( A( face(2,:) + (face(3,:)-1)*n ) );
v31 = full( A( face(3,:) + (face(1,:)-1)*n ) );

Compute the new faces, each old face generates 4 faces.

In [9]:
face = [  cat(1,face(1,:),v12,v31),...
    cat(1,face(2,:),v23,v12),...
    cat(1,face(3,:),v31,v23),...
    cat(1,v12,v23,v31)   ];

Add new vertices at the edges center.

In [10]:
vertex = [vertex, (vertex(:,edge(1,:))+vertex(:,edge(2,:)))/2 ];

Project the points on the sphere.

In [11]:
d = sqrt( sum(vertex.^2,1) );
vertex = vertex ./ repmat( d, [size(vertex,1) 1]);

Display before/after subdivision.

In [12]:
clf;
subplot(1,2,1);
plot_mesh(vertex0,face0); 
shading('faceted'); lighting('flat'); view(3); axis('tight');
subplot(1,2,2);
plot_mesh(vertex,face); 
shading('faceted'); lighting('flat'); view(3); axis('tight');

Exercise 1

Perform the full subdivision.

In [13]:
exo1()
In [14]:
%% Insert your code here.

Exercise 2

Try with other control meshes.

In [15]:
exo2()
In [16]:
%% Insert your code here.

Triangulated Mesh Subdivision

The same method can be applied to an arbitrary control mesh, but without the projection on the sphere. More clever interpolations should be used to avoid having a simple piecewise linear surface.

Load the base control mesh.

In [17]:
name = 'mannequin';
[vertex0,face0] = read_mesh(name);

Display it.

In [18]:
options.name = name;
clf;
plot_mesh(vertex0,face0,options);
shading('faceted'); lighting('flat'); axis('tight');

Initialize.

In [19]:
face = face0;
vertex = vertex0;

Perform the subdivision.

In [20]:
edge = compute_edges(face);
n = size(vertex,2);
ne = size(edge,2);

Compute the number of the three edges associated to each face.

In [21]:
A = sparse([edge(1,:);edge(2,:)],[edge(2,:);edge(1,:)],[n+(1:ne);n+(1:ne)],n,n);
v12 = full( A( face(1,:) + (face(2,:)-1)*n ) );
v23 = full( A( face(2,:) + (face(3,:)-1)*n ) );
v31 = full( A( face(3,:) + (face(1,:)-1)*n ) );

Compute the new faces, each old face generates 4 faces.

In [22]:
face_old = face;
face = [  cat(1,face(1,:),v12,v31),...
    cat(1,face(2,:),v23,v12),...
    cat(1,face(3,:),v31,v23),...
    cat(1,v12,v23,v31)   ];

Compute the vertex and face ring.

In [23]:
global vring e2f fring facej;
vring = compute_vertex_ring(face);
e2f = compute_edge_face_ring(face_old);
fring = compute_face_ring(face_old);
facej = face_old;

Compute the interpolated position using

In [24]:
for k=n+1:n+ne
    [e,v,g] = compute_butterfly_neighbors(k, n);
    vertex(:,k) = 1/2*sum(vertex(:,e),2) + 1/8*sum(vertex(:,v),2) - 1/16*sum(vertex(:,g),2);
end

Display before/after subdivision.

In [25]:
clf;
subplot(1,2,1);
plot_mesh(vertex0,face0,options); 
shading('faceted'); lighting('flat'); axis('tight');
subplot(1,2,2);
plot_mesh(vertex,face,options); 
shading('faceted'); lighting('flat'); axis('tight');