Abstract: Accurately modeling BRDF for real world materials is important and challenging for realistic image synthesis. For a majority of materials most of the incident light enters the material, undergoes multiple scattering under the surface before exiting the material's surface as reflection. Physically correct modeling of such BRDF must take into account of this subsurface volumetric light transport. Most of the accurate numerical solution methods (ex: Monte Carlo, Discrete Ordinate Methods (DOM)) for volumetric light transport compute radiance field for the whole volume, and are expensive. As BRDF ultimately relates only the outgoing radiation field at the boundary to the incident radiation, radiation field computed for the bulk of the material does not provide any useful information and hence the effort involved in computing them can be considered as wasteful. So for efficient BRDF computation any method that allows us to compute the radiance field only at the boundary would be a preferable choice. The search for such a method led us to the Ambartsumian's method [Sobolev 1975; Mishchenko et al. 1999].