A numerical method to the inverse problem in shape-based diffusion optical tomography is proposed to simultaneously recover the smooth boundaries of the tissue regions and the optical properties. The regions of different tissues are assumed to have piecewise constant optical properties,thus the forward problem of the light propagation of diffuse optical tomography (DOT) can be modeled as a set of coupled Helmholtz equations,and solved by the boundary element method (BEM). Fourier series expansion is used for the representation of the smooth complicated boundaries. For inverse problem the Levenberg-Marquardt optimization process is implemented here. The performance of the proposed approach is evaluated by the measurements at different noise levels. The numerical results illuminate that the methodology has faster convergence speed and global convergence,and the boundaries and the optical coefficients can both be recovered with the good accuracy from the noisy measurements.