A novel approach is presented for obtaining fast robust three-dimensional (3-D) reconstructions of bioluminescent reporters buried deep inside animal subjects from multispectral images of surface bioluminescent photon densities. The proposed method iteratively acts upon the equations relating the multispectral data to the luminescent distribution with high computational efficiency to provide robust 3-D reconstructions. Unlike existing algebraic reconstruction techniques, the proposed method is designed to use adaptive projections that iteratively guide the updates to the solution with improved speed and robustness. Contrary to least-squares reconstruction methods, the proposed technique does not require parameter selection or optimization for optimal performance. Additionally, optimized schemes for thresholding, sampling, and ordering of the bioluminescence tomographic data used by the proposed method are presented. The performance of the proposed approach in reconstructing the shape, volume, flux, and depth of luminescent inclusions is evaluated in a multitude of phantom-based and dual-modality in vivo studies in which calibrated sources are implanted in animal subjects and imaged in a dual-modality optical/computed tomography platform. Statistical analysis of the errors in the depth and flux of the reconstructed inclusions and the convergence time of the proposed method is used to demonstrate its unbiased performance, low error variance, and computational efficiency.