We present a method for the numerical evaluation of 6D singular integrals appearing in Volume Integral Equations. It is an extension of the Sauter-Schwab/Taylor-Duffy strategy for singular triangle-triangle interaction integrals to singular tetrahedron-tetrahedron interaction integrals. This general approach allows to use different kinds of kernel and basis functions. It also works on curvilinear domains. Our approach is based on relative coordinates and splitting the integration domain into subdomains for which quadrature rules can be constructed. Further, we show how to build these tensor-product quadrature rules economically using quadrature rules defined over 2D, 3D and 4D simplices. Compared to the existing approach where the integral is computed as a sequence of 1D integrations significant speedup can be achieved. The accuracy and convergence properties of the method are demonstrated by numerical experiments.