In this thesis, we investigate the computational methods for both diffuse and general reflections in realistic image synthesis and propose two new approaches: the overrelaxation solution and the Bernstein polynomial solution. One of the major concerns with the radiosity method is its expensive computing time and memory requirements. In this thesis, we analyze the convergence behavior of the progressive refinement radiosity method and propose two overrelaxation algorithms: the gathering and shooting solution and the positive overshooting solution. We modify the conventional shooting method to make the optimal use of the visibility information computed in each iteration. Based on a concise record of the history of the unshot light energy distribution, a solid convergence speed-up is achieved. Though a great effort has been made to extend the radiosity method to accommodate general non-diffuse reflection, the current algorithms are still quite limited to simple environment settings. In this thesis, we propose using the piecewise spherical Bernstein basis functions over a geodesic triangulation to represent the radiance function. The representation is intrinsic to the unit sphere, and can be efficiently stored, evaluated, and subdivided by the numerically stable de Casteljau algorithm. We demonstrate that the computation of other fundamental radiometric quantities such as vector irradiance and reflected radiance can be reduced to the integration of the piecewise spherical Bernstein basis functions. A novel geometric integration algorithm based on adaptive domain subdivision is presented for the Bernstein-B´ezier polynomials over a geodesic triangle on the unit sphere.