Three-dimensional (3D) neutron imaging is the cutting-edge diagnostic method for pulsed radiation sources driven by laser or Z-pinch. Neutron images with different energies can accurately diagnose the spatial structure of the burning Deuterium-Tritium (DT) plasma and surrounding materials under high compression ratio, which is of great significance for improving the technical measures of compression symmetry, verifying the physical model and evaluating the ignition performance. Utilized the computed tomography method, 3D information of the target internal structure can be acquired by multi-view neutron imaging method. Along with the limitations of current technology, there are great difficulties in pixel level accurate matching, high-precision time correlation and camera sensitivity difference correction between different imaging axes due to the high radiation intensities, short durations, and wide energy ranges of the pulsed radiation sources. Thus, only a limited number of axis can be implemented.Reconstructing 3D images from severely incomplete two-dimensional data is an ill-posed problem, which suffers from large solution space, artifacts and inability to maintain edge features. JIANG Shaoen et al developed a 3D reconstruction algorithm based on algebraic reconstruction techniques. VOLEGOV P L et al used an iterative Expectation-Maximization (EM) algorithm to reconstruct 3D neutron or X-ray sources. The algorithm reconstructed the main geometrical features of the source, but introduced artifacts affected by the angular distribution of projection directions. VOLEGOV P L et al later adopted Spherical Harmonic Decomposition (SHD) and Cylindrical Harmonic Decomposition (CHD) methods to the reconstruction of 3D radiation sources. While these analytical methods based on basis functions expansion were fast and could acquire the unique solution under certain conditions, their representation abilities were limited by the number of imaging axes and would produce artifacts related to the truncation of expansion order.In this paper, we introduce and realize SHD reconstruction algorithm and EM iterative algorithm. Considering that the inertial confinement fusion radiation sources have some spherical symmetry and the EM algorithm can find local optimum solution, we exploit the SHD reconstruction results as the initial value of EM iterative algorithm. This method actually utilizes a physical prior based on the target physical characters overcoming the question of large solution space for the limited-view 3D reconstruction algorithm in some sense. To evaluate the performance of the proposed method, an ellipsoidal radiation source with a Gaussian intensity distribution is designed as the target source to be reconstructed. Compared with the results of SHD, OSEM and MF-OSEM algorithms, the reconstruction results of this method are closest to the ideal ellipsoid source. Our method also achieves the best results in terms of three quantitative evaluating indices, which are 36.545 4 (PSNR), 0.014 9 (RMSE) and 0.031 5 (D_KL), respectively. The results indicate that the EM iterative algorithm with initial values constrained by SHD algorithm exhibits better performances and better adaptability to noise compared with the SHD or EM algorithm. Our method can be used as a benchmark comparison algorithm for the reconstruction of similar 3D radiation sources with few-view projections.