Fluorescence molecular tomography in the near-infrared region is becoming a powerful modality for mapping the three-dimensional quantitative distributions of fluorochromes in live small animals. However, wider application of fluorescence molecular tomography still requires more accurate and stable reconstruction tools. We propose a shape-based reconstruction method that uses spherical harmonics parameterization, where fluorophores are assumed to be distributed as piecewise constants inside disjointed subdomains and the remaining background. The inverse problem is then formulated as a constrained nonlinear least-squares problem with respect to shape parameters, which decreases ill-posedness because of the significantly reduced number of unknowns. Since different shape parameters contribute differently to the boundary measurements, a two-step and modified block coordinate descent optimization algorithm is introduced to stabilize the reconstruction. We first evaluated our method using numerical simulations under various conditions for the noise level and fluorescent background; it showed significant superiority over conventional voxel-based methods in terms of the spatial resolution, reconstruction accuracy with regard to the morphology and intensity, and robustness against the initial estimated distribution. In our phantom experiment, our method again showed better spatial resolution and more accurate intensity reconstruction. Finally, the results of an in vivo experiment demonstrated its applicability to the imaging of mice.