We proposed an efficient and automatic method to reconstruct hierarchical tetrahedral meshes and build density models for bony anatomy from CT data sets. The mesh is reconstructed from contours extracted in CT images corresponding to the outer bone surfaces and the boundaries between compact bone, spongy bone and medullary cavity. Key problems, such as tiling problem, branching problem, correspondence problem, and constraint problem, have been solved. We then approximated bone density variations by means of continuous density functions in each tetrahedron. The density functions are written as smooth Bernstein polynomial spline expressed in terms of barycentric coordinates associated with each tetrahedron. We further performed the tetrahedral mesh simplification by edge collapsing and built hierarchical structure of multiple resolution meshes. We applied our density model to efficiently generate Digitally Reconstructed Radiographs. This research is part of our effort to build a bone density atlas.