Search code examples
matlabimage-processing3dgeometrysurface

Reconstructing faces of an upsampled icosahedron


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. 

Solution

  • 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.