We present a new construction of fractal interpolation surfaces defined on arbitrary rectangular lattices. We use this construction to form finite sets of fractal interpolation functions (FIFs) that generate multiresolution analyses of L(2)(R(2)) of multiplicity r. These multiresolution analyses are based on the dilation properties of the construction. The associated multi-wavelets are orthogonal and discontinuous functions. We give concrete examples to illustrate the method and generalize it to form multiresolution analyses of L(2)(R(d)), d > 2. To this end, we prove some results concerning the Holder exponent of FIFs defined on [0, 1](d).