subdivideMesh

SUBDIVIDEMESH Subdivides each face of the mesh.

   [V2 F2] = subdivideMesh(V, F, N)
   Subdivides the mesh specified by (V,F) such that each face F is divided
   into N^2 smaller faces.

   Example
     [v, f] = createOctahedron;
     figure; drawMesh(v, f); view(3);
     [v2, f2] = subdivideMesh(v, f, 4);
     figure; drawMesh(v2, f2); view(3)

   See also
     meshes3d, drawMesh

0001 function [vertices2, faces2] = subdivideMesh(vertices, faces, n)
0002 %SUBDIVIDEMESH Subdivides each face of the mesh.
0003 %
0004 %   [V2 F2] = subdivideMesh(V, F, N)
0005 %   Subdivides the mesh specified by (V,F) such that each face F is divided
0006 %   into N^2 smaller faces.
0007 %
0008 %   Example
0009 %     [v, f] = createOctahedron;
0010 %     figure; drawMesh(v, f); view(3);
0011 %     [v2, f2] = subdivideMesh(v, f, 4);
0012 %     figure; drawMesh(v2, f2); view(3)
0013 %
0014 %   See also
0015 %     meshes3d, drawMesh
0016 %
0017 
0018 % ------
0019 % Author: David Legland
0020 % e-mail: david.legland@inra.fr
0021 % Created: 2013-08-22,    using Matlab 7.9.0.529 (R2009b)
0022 % Copyright 2013 INRA - Cepia Software Platform.
0023 
0024 
0025 %% Initialisations
0026 
0027 edges = [];
0028 faceEdgeIndices = [];
0029 
0030 if isstruct(vertices)
0031     % get relevant inputs
0032     mesh = vertices;
0033     n = faces;
0034     
0035     % parse fields from a mesh structure
0036     vertices = mesh.vertices;
0037     faces = mesh.faces;
0038     if isfield(mesh, 'edges')
0039         edges = mesh.edges;
0040     end
0041     if isfield(mesh, 'faceEdges')
0042         faceEdgeIndices = mesh.faceEdges;
0043     end
0044     
0045 end
0046 
0047 if ~isnumeric(faces) || size(faces, 2) ~= 3
0048     error('Requires a triangular mesh');
0049 end
0050 
0051 % compute the edge array
0052 if isempty(edges)
0053     edges = meshEdges(faces);
0054 end
0055 nEdges = size(edges, 1);
0056 
0057 % index of edges around each face
0058 if isempty(faceEdgeIndices)
0059     faceEdgeIndices = meshFaceEdges(vertices, edges, faces);
0060 end
0061 
0062 
0063 %% Create new vertices on edges
0064 
0065 % several interpolated positions
0066 t = linspace(0, 1, n + 1)';
0067 coef2 = t(2:end-1);
0068 coef1 = 1 - t(2:end-1);
0069 
0070 % initialise the array of new vertices
0071 vertices2 = vertices;
0072 
0073 % keep an array containing index of new vertices for each original edge
0074 edgeNewVertexIndices = zeros(nEdges, n-1);
0075 
0076 % create new vertices on each edge
0077 for iEdge = 1:nEdges
0078     % extract each extremity as a point
0079     v1 = vertices(edges(iEdge, 1), :);
0080     v2 = vertices(edges(iEdge, 2), :);
0081 
0082     % compute new points
0083     newPoints = coef1 * v1 + coef2 * v2;
0084     
0085     % add new vertices, and keep their indices
0086     edgeNewVertexIndices(iEdge,:) = size(vertices2, 1) + (1:n-1);
0087     vertices2 = [vertices2 ; newPoints]; %#ok<AGROW>
0088 end
0089 
0090 
0091 %% Process each face
0092 
0093 faces2 = zeros(0, 3);
0094 
0095 nFaces = size(faces, 1);
0096 for iFace = 1:nFaces
0097     % compute index of each corner vertex
0098     face = faces(iFace, :);
0099     iv1 = face(1);
0100     iv2 = face(2);
0101     iv3 = face(3);
0102     
0103     % compute index of each edge
0104     faceEdges = faceEdgeIndices{iFace};
0105     ie1 = faceEdges(1);
0106     ie2 = faceEdges(2);
0107     ie3 = faceEdges(3);
0108     
0109     % indices of new vertices on edges
0110     edge1NewVertexIndices = edgeNewVertexIndices(ie1, :);
0111     edge2NewVertexIndices = edgeNewVertexIndices(ie2, :);
0112     edge3NewVertexIndices = edgeNewVertexIndices(ie3, :);
0113     
0114     % keep vertex 1 as reference for edges 1 and 3
0115     if edges(ie1, 1) ~= iv1
0116         edge1NewVertexIndices = edge1NewVertexIndices(end:-1:1);
0117     end
0118     if edges(ie3, 1) ~= iv1
0119         edge3NewVertexIndices = edge3NewVertexIndices(end:-1:1);
0120     end
0121        
0122     % create the first new face, on 'top' of the original face
0123     topVertexInds = [edge1NewVertexIndices(1) edge3NewVertexIndices(1)];
0124     newFace = [iv1 topVertexInds];
0125     faces2 = [faces2; newFace]; %#ok<AGROW>
0126         
0127     % iterate over middle strips
0128     for iStrip = 2:n-1
0129         % index of extreme vertices of current row
0130         ivr1 = edge1NewVertexIndices(iStrip);
0131         ivr2 = edge3NewVertexIndices(iStrip);
0132         
0133         % extreme vertices as points
0134         v1 = vertices2(ivr1, :);
0135         v2 = vertices2(ivr2, :);
0136         
0137         % create additional vertices within the bottom row of the strip
0138         t = linspace(0, 1, iStrip+1)';
0139         coef2 = t(2:end-1);
0140         coef1 = 1 - t(2:end-1);
0141         newPoints = coef1 * v1 + coef2 * v2;
0142 
0143         % compute indices of new vertices in result array
0144         newInds = size(vertices2, 1) + (1:iStrip-1);
0145         botVertexInds = [ivr1 newInds ivr2];
0146         
0147         % add new vertices
0148         vertices2 = [vertices2 ; newPoints]; %#ok<AGROW>
0149         
0150         % create top faces of current strip
0151         for k = 1:iStrip-1
0152             newFace = [topVertexInds(k) botVertexInds(k+1) topVertexInds(k+1)];
0153             faces2 = [faces2; newFace]; %#ok<AGROW>
0154         end
0155         
0156         % create bottom faces of current strip
0157         for k = 1:iStrip
0158             newFace = [topVertexInds(k) botVertexInds(k) botVertexInds(k+1)];
0159             faces2 = [faces2; newFace]; %#ok<AGROW>
0160         end
0161         
0162         % bottom vertices of current strip are top vertices of next strip
0163         topVertexInds = botVertexInds;
0164     end
0165         
0166     % for edge 2, keep vertex 2 of the current face as reference
0167     if edges(ie2, 1) ~= iv2
0168         edge2NewVertexIndices = edge2NewVertexIndices(end:-1:1);
0169     end
0170     
0171     % consider new vertices together with extremities
0172     botVertexInds = [iv2 edge2NewVertexIndices iv3];
0173     
0174     % create top faces for last strip
0175     for k = 1:n-1
0176         newFace = [topVertexInds(k) botVertexInds(k+1) topVertexInds(k+1)];
0177         faces2 = [faces2; newFace]; %#ok<AGROW>
0178     end
0179     
0180     % create bottom faces for last strip
0181     for k = 1:n
0182         newFace = [topVertexInds(k) botVertexInds(k) botVertexInds(k+1)];
0183         faces2 = [faces2; newFace]; %#ok<AGROW>
0184     end
0185     
0186 end
0187 
0188 
0189 
0190 
0191 
0192 
0193 
0194 
0195 
0196 
0197 
0198 
0199 
0200 
0201