Reconstructing gamma-ray spectra from detector measurements is an ill-posed inverse problem. Standard methods, such as Folding Iteration with Compton Subtraction (FICS), provide point estimates but lack calibrated uncertainties and may bias the spectrum. We introduce an unfolding framework based on regularized maximum-likelihood estimation (RMLE) that enforces nonnegativity and detector-response constraints while explicitly modeling background and contaminant contributions. Simulations and analytical results show that RMLE yields smoother reconstructions with well-calibrated confidence intervals and outperforms existing techniques for low-complexity spectra. Although high-complexity data remain challenging, the intervals produced by RMLE maintain correct coverage.