We present a practical, timely, and effective radiative transfer algorithm, suitable for qualitative and quantitative analyses of high-resolution hyperspectral images of planetary surfaces in the visible and near-infrared domains. The bidirectional reflectance of a plane parallel, absorbing, scattering, and slightly stratified medium is generated. The local mean properties of scattering and absorption of such media are obtained apart, using semiempirical approaches. The functions which express the diffuse reflection and transmission behaviors of each homogeneous layer are then derived. For the multiple scattering term, we numerically resolve the equations appearing in the H, X and Y function method of radiative transfer, reducing the real phase function to a simplified one which can nevertheless be anisotropic. A better approach to the physical realism is obtained for the single and double scattering contributions, using their real analytical expressions. This contrasts with the Hapke model dedicated to homogeneous and semi-infinite media, where only an isotropic reduced phase function is adopted and the single scattering correction is applied. The bidirectional reflectance and the derived quantities (albedos) of an optically semi-infinite homogeneous medium are then easily derived from these quantities. For a stratified medium, a simple adding algorithm based on principles of invariance is presented. Compared to earlier and more complete theoretical developments, this model in most cases reproduces the dependence of the bidirectional reflectance according to the different geometrical and radiative parameters with a maximum of 10% relative error. It leads to important gains of computation time and significantly extends the validity of Hapke's or similar practical approaches. ¿ 1998 American Geophysical Union