Even though isogeometric analysis has a clear advantage regarding the number of degrees of freedom needed to achieve a certain level of accuracy, the time needed for matrix assembly (by means of numerical integration of stiffness or mass matrix entries) constitutes a severe bottleneck in the process. In order to overcome this difficulty, we propose the new approach of Integration by Interpolation and Lookup (IIL). Firstly, applying spline interpolation to the common factors in the occurring integrals approximately transforms them into integrals of piecewise polynomial functions, whose integrands are expressed in tensor-product B-spline form. The common factors represent the influence of the geometry mapping (i.e., the NURBS domain parameterization) and the contributions of possibly non-constant material coefficients. Secondly, these integrals are evaluated exactly using pre-computed look-up tables for integrals of tri-products of univariate B-splines and their derivatives. For the model case of elliptic problems, we perform a theoretical analysis to demonstrate that the IIL method maintains the overall approximation order of the Galerkin discretization, provided that the spline interpolation is sufficiently accurate. Moreover, we provide a comparison of the computational complexity of our method with that of a standard Gauss quadrature method. Finally, we present experimental results to illustrate the performance of the ILL method and to support our theoretical results.