In this paper we propose using planar and spherical Bernstein polynomials over triangular domain for radiative transfer computations. In the planar domain, we propose using piecewise Bernstein basis functions and symmetric Gaussian quadrature formulas over triangular elements for high quality radiosity solution. In the spherical domain, we propose using piecewise Bernstein basis functions over a geodesic triangulation to represent the radiance function. The representation is intrinsic to the unit sphere, and may be efficiently stored, evaluated, and subdivided by the de Casteljau algorithm. The computation of other fundamental radiometric quantities such as vector irradiance and reflected radiance may be reduced to the integration of the piecewise Bernstein basis functions on the unit sphere. The key result of our work is a simple geometric integration algorithm based on adaptive domain subdivision for the Bernstein-Bézier polynomials over a geodesic triangle on the unit sphere.