In this paper we propose a class of high order numerical methods to delta function integrals appearing in level set methods in the three dimensional case by extending the idea for designing high order methods to two dimensional delta function integrals in [X. Wen, J. Comput. Phys., 228 (2009), pp. 4273–4290]. The methods comprise approximating the mesh cell restrictions of the delta function integral. In each mesh cell the three dimensional delta function integral can be rewritten as a two dimensional ordinary integral with the smooth integrand being a one dimensional delta function integral. A key idea in designing high order methods in this paper is that the high order accuracy of the methods is not ensured in each mesh cell, and the methods are designed in a consistent way to foster error cancelations in mesh cells where high order accuracy is not ensured in the individual cell. This issue is essentially related to the construction of integral area for the two dimensional ordinary...