A numerical method is presented for calculating neutron transport problems in three-dimensional (x,y,z) geometry on the basis of a method of direct integration of the integral transport equation. Several new techniques are introduced to the method to make it well adapted to practical neutron transport calculations in three-dimensional geometry. A technique for evaluating the scattering source based on an estimated spectral shape in each material region allows use of coarse energy mesh intervals without reducing calculational accuracy as compared with the calculation with fine meshes. A quadratic function approximation for the source spatial distribution in each spatial mesh interval is found to improve the mathematical error in direct integration of the source term over the spatial variable as compared with the linear- or exponential-function approximation used in the original method. In addition, Lagrange's interpolation formula is applied instead of the linear interpolation used in the original method for more accurate estimation of both flux and source. Comparisons are made of the calculations with experiments for three neutron transport problems, the pool critical assembly experiment, the Winfrith iron benchmark experiment, and the annular duct neutron streaming experiment, and also with the three-dimensional Sn calculation to verify the validity of the present method for neutron transport calculations in (x,y,z) geometry.