An easy-to-implement variance reduction technique based on the control variate method is proposed and demonstrated for Monte Carlo (MC) calculations of surface-to-surface (S-to-S) radiation exchange. The method is general, and is applicable to any S-to-S radiation exchange problem-gray or non-gray. The control function used in the proposed method is the Planck function evaluated at the minimum temperature of the domain. The Radiative Transfer Equation (RTE) and its boundary conditions are recast into a deviational form, and it is shown that only two minor changes are necessary within the standard S-to-S MC algorithm to implement the proposed variance reduction technique. No computational overheads are incurred is using the proposed technique. Three test cases are considered, and exact (deterministic) solutions are used to assess the accuracy of the MC results. It is found that in all cases, the statistical error in the computed heat flux, as measured by the variance, is reduced. The variance reduction can be substantial (almost an order of magnitude) if the problem at hand is dominated by temperatures sufficiently close to the control temperature. In the third test case, radiation is also iteratively coupled to convection, and it is shown that the variance reduced MC solutions for steady state temperatures are far more accurate than those produced by the standard MC method for the same number of rays. (C) 2018 Elsevier Ltd. All rights reserved.