The 3D surface to be discretized is represented by a triangular control grid of arbitrary topology. The boundary edges of the grid form the boundary curves of the surface. The end nodes of each curve represent vertices. Nodes not representing vertices are called curve nodes if they form a curve, otherwise they are called surface nodes.
In order to obtain smooth triangulation (independent of the control grid e.g. due to different element size specifications), a smooth representation of the discrete surface (typically continuity is sufficient), called the limit surface, is derived. The limit surface is reconstructed using a suitable subdivision technique which is based on the hierarchical recursive refinement of triangular simplices forming the control grid (Fig. 1). Each step of the refinement consists of two stages -- splitting and averaging (Fig. 2). In the splitting stage, new nodes are introduced exactly in the middle of individual edges. During the averaging, the nodes are repositioned to a new location evaluated as a weighted average of nodes in the neighbourhood (according to the so called averaging mask). As the level of the refinement grows the resulting grid approaches the (smooth) limit surface. The averaging scheme may be either interpolating or approximating. While an interpolating scheme produces the limit surface which is passing through the nodes of the control grid (original and refined), the limit surface created by an approximating scheme does not generally interpolate any of the nodes. The other classification  of the averaging scheme may be uniform (if the same averaging mask is used everywhere on the surface) or non-uniform (if not), and stationary (if for a given node the same averaging mask is used on each level of the refinement) or non-stationary (if not).
In the presented approach, the recursive subdivision based on the modified Butterfly scheme  has been employed. It is the interpolating non-uniform stationary scheme, in which the existing nodes (on the current level of refinement) remain unchanged and the position of a new node on the next level (see Fig. 3a) is calculated as
Similarly, the limit boundary curves are recovered using a one-dimensional interpolating subdivision  producing continuous curves. The adopted 4-point (for a new node between two curve nodes) and 3-point (for a new node between vertex node and curve node) averaging masks are depicted in Figures 3b and 3c.
The final interpolating procedure evaluates the position of a new node according to the classification and regularity of the end nodes of its parent edge (a surface node of valence 6 is called regular, otherwise it is called irregular):