We analyse the approximation and smoothness properties of fundamental and refinable functions that arise from interpolatory subdivision schemes in multidimensional spaces. In particular, we provide a general way for the construction of bivariate interpolatory refinement masks such that the corresponding fundamental and refinable functions attain the optimal approximation order and smoothness order. In addition, these interpolatory refinement masks are minimally supported and enjoy full symmetry. Several examples are explicitly computed.