Computing a source-target map that yields integrable surface normal field is quite challenging for freeform illumination design. Here, we propose a least-squares ray mapping method to calculate a superior ray mapping by iteratively correcting an integrable map to approach the energy conservation and boundary condition. The process is implemented via solving three minimization problems. The first two problems can be figured out pointwise and the third can be converted to two decoupled Poisson equations with Robin boundary conditions. We demonstrate the robustness and high efficiency of the proposed method with several design examples.