In this paper, we give a complete description of efficient formulae for the numerical integration of fast oscillating functions of two variables. The focus is on the case encountered frequently in many engineering applications where an accurate value of the Lipschitz constant is not available. Using spline approximations, we demonstrate the main idea of our approach on the example of piecewise bilinear interpolation, and propose optimal-by-order (with a constant not exceeding two) cubature formulae that are applicable for a wide range of oscillatory patterns. This property makes the formulae indispensable in many engineering applications dealing with signal processing and image recognition. Illustrative results of numerical experiments are presented.