In this study, a 4-D, computerized ionospheric tomography algorithm, IONOLAB-Fusion, is developed to reconstruct electron density using both actual and virtual vertical and horizontal paths for all ionospheric states. The user-friendly algorithm only requires the coordinates of the region of interest and range with the desired spatio-temporal resolutions. The model ionosphere is formed using spherical voxels in a lexicographical order so that a 4-D ionosphere can be mapped to a 2-D matrix. The model matrix is formed automatically using a background ionospheric model with an optimized retrospective or near-real time manner. The singular value decomposition is applied to extract a subset of significant singular values and corresponding signal subspace basis vectors. The measurement vector is filled automatically with the optimized number of ground-based and space-based paths. The reconstruction is obtained in closed form in the least squares sense. When the performance of IONOLAB-Fusion across Europe was compared with ionosonde profiles, a 26.51% and 32.33% improvement was observed over the background ionospheric model for quiet and disturbed days, respectively. When compared with GIM-TEC, the agreement of IONOLAB-Fusion was 37.89% and 31.58% better than those achieved with the background model for quiet and disturbed days, respectively.