In particle-tracing radiosity algorithms, energy-carrying particles are traced through an environment for simulating global illumination. Illumination on a surface is reconstructed from particle “hit points” on the surface, which is a density estimation problem [11]. Several methods can be used to solve this problem, such as the adaptive meshing method [14], the kernel method [15], and the orthogonal series estimator [3]. In this paper, a new orthogonal series estimator is proposed to tackle the problem. In the new method, the appropriate number of terms that should be used in the series is determined adaptively and automatically. Moreover, a surface subdivision scheme is combined with the estimator to increase the accuracy of estimation. The new method has several advantages over other existing methods: (1) it requires less memory than the adaptive meshing method; (2) it does not store all the particle-hit points as in the kernel method; (3) it determines automatically how many terms should be used in the orthogonal series; (4) it incorporates surface subdivision to further in-crease the accuracy of estimation.