Spherical harmonics transform plays an important role in research in shape description. Current methods to compute the spherical harmonics decomposition of the characteristic function of the intersection of a polyhedral solid with a sphere involve expensive voxelization, and are prone to numerical errors associated with the size of the voxels. This paper describes a fast and accurate technique for computing spherical harmonics coefficients directly from the description of the mesh. The algorithm runs in linear time O(pn), where n is the number of triangles of the mesh and p is the number of terms calculated, which is roughly linear.