I am attempting to up-sample an icosahedron in MATLAB by adding new vertices at the mid-points of each edge and recomputing the faces. I've managed to place the new vertices, but I can't figure out how to recompute the new faces.
I want each vertex to connect to its closest neighbours to form a new, up-sampled icosahedron.
Sample data:
S.vertices = ...
[0.5000 0.1625 0.8507;
0.5000 0.6882 0.5257;
0.8507 0.2764 0.4472;
0.8090 0.5878 0.0000;
0.3090 0.9511 0.0000;
0.0000 0.5257 0.8507;
-0.5000 0.6882 0.5257;
0.0000 0.8944 0.4472;
-0.3090 0.9511 0.0000;
-0.8090 0.5878 0.0000;
-0.5000 0.1625 0.8507;
-0.8090 -0.2629 0.5257;
-0.8507 0.2764 0.4472;
-1.0000 0.0000 0.0000;
-0.8090 -0.5878 0.0000;
-0.3090 -0.4253 0.8507;
-0.0000 -0.8507 0.5257;
-0.5257 -0.7236 0.4472;
-0.3090 -0.9511 0.0000;
0.3090 -0.9511 0.0000;
0.3090 -0.4253 0.8507;
0.8090 -0.2629 0.5257;
0.5257 -0.7236 0.4472;
0.8090 -0.5878 0.0000;
1.0000 -0.0000 0.0000;
0.0000 0.0000 1.0000;
0.5000 0.1625 -0.8507;
0.5000 0.6882 -0.5257;
0.8507 0.2764 -0.4472;
0.0000 0.5257 -0.8507;
-0.5000 0.6882 -0.5257;
0.0000 0.8944 -0.4472;
-0.5000 0.1625 -0.8507;
-0.8090 -0.2629 -0.5257;
-0.8507 0.2764 -0.4472;
-0.3090 -0.4253 -0.8507;
-0.0000 -0.8507 -0.5257;
-0.5257 -0.7236 -0.4472;
0.3090 -0.4253 -0.8507;
0.8090 -0.2629 -0.5257;
0.5257 -0.7236 -0.4472;
0.0000 0.0000 -1.0000];
S.faces = ...
[13 14 12;
12 15 18;
1 3 2;
25 22 24;
12 16 11;
22 23 24;
17 21 16;
18 17 16;
23 17 20;
18 16 12;
26 6 11;
11 13 12;
14 15 12;
19 18 15;
17 23 21;
8 9 7;
10 14 13;
7 9 10;
13 7 10;
11 7 13;
6 8 7;
21 26 16;
6 7 11;
2 8 6;
2 5 8;
8 5 9;
26 1 6;
4 5 2;
6 1 2;
21 23 22;
2 3 4;
1 21 22;
3 25 4;
16 26 11;
21 1 26;
19 17 18;
19 20 17;
23 20 24;
3 22 25;
3 1 22;
27 29 28;
34 36 33;
37 39 36;
38 37 36;
38 36 34;
42 30 33;
33 35 34;
37 41 39;
33 31 35;
30 32 31;
39 42 36;
30 31 33;
28 32 30;
42 27 30;
30 27 28;
39 41 40;
27 39 40;
36 42 33;
39 27 42;
29 27 40;
35 14 34;
34 15 38;
25 40 24;
40 41 24;
41 37 20;
14 15 34;
19 38 15;
32 9 31;
10 14 35;
31 9 10;
35 31 10;
28 5 32;
32 5 9;
4 5 28;
28 29 4;
29 25 4;
19 37 38;
19 20 37;
41 20 24;
29 40 25];
Code:
% Center icosahedron on the origin.
offset = mean(S.vertices);
S.vertices = S.vertices - offset;
% Calculate radius.
radius = mean(pdist2([0 0 0], S.vertices));
% Get all edges.
edges = [S.faces(:,1),S.faces(:,2); ...
S.faces(:,1),S.faces(:,3); ...
S.faces(:,2),S.faces(:,3)];
% Find midpoints of each edge.
V = (S.vertices(edges(:,1),:) + S.vertices(edges(:,2),:)) / 2;
% Convert to spherical, fix radius, and convert back to cartesian.
[theta,phi] = cart2sph(V(:,1),V(:,2),V(:,3));
[V(:,1),V(:,2),V(:,3)] = sph2cart(theta,phi,radius);
% Add new points, fix faces.
S.vertices = [S.vertices;V];
% Some clever trick to build the new faces matrix.
For a triangle {1, 2, 3}
, subdividing will create 3 new vertices {a, b, c}
:
1
o
/ \
/ \
a o-----o c
/ \ / \
/ \ / \
o-----o-----o
2 b 3
It will also create 4 new faces:
[c 1 a], [a 2 b], [b 3 c], and [a b c]
Your code creates the coordinates of the new vertices in V
such that the a
vertices for all of the triangles come first, followed by the c
vertices, and then the b
vertices.
First let's generate the indices of these vertices. The new indices start from size(S.vertices)+1
and there are a total of 3*size(S.faces)
of them. Let's also group the new vertices for a each existing face together.
% before updating S.vertices...
Vx = reshape(1:size(V,1), [], 3) + size(S.vertices, 1);
This gives us all of the new [a b c]
faces (actually, in [a c b]
order) in Vx(f,:)
with the a
vertices in Vx(:,1)
, the b
vertices in Vx(:,3)
<-- not 2!, and the c
vertices in Vx(:,2)
. Now we need to create the other 3 new faces for each existing face. You can do this using cat
, permute
and reshape
, but it gets really ugly and incomprehensible. Since we've only got 9 vertices that go into these faces, it's easier to just plug them in manually to create a m x 9
matrix and then reshape into m*3 x 3
:
% face format: [c 1 a, a 2 b, b 3 c]
Fx = [Vx(:,2) S.faces(:,1) Vx(:,1)
Vx(:,1) S.faces(:,2) Vx(:,3)
Vx(:,3) S.faces(:,3) Vx(:,2)];
% now reshape to [c 1 a; a 2 b; b 3 c]
Fx = reshape(Fx.', 3, []).';
Now all that's left is updating S
:
S.vertices = [S.vertices;V];
S.faces = [S.faces; Fx; Vx]; % oops... you probably don't want the old faces in there. they're being replaced.