In this paper, we propose a new algorithm for the fundamental problem of reconstructing surfaces from a large set of unorganized 3D data points. The local shapes of the surface are recovered by variational implicit surface represented as a weighted combination of radial basis functions. The variational implicit patches are then combined together to form the overall surface via a set of blending functions, which is also referred to as the partition-of-unity method. The reconstruction algorithm first partitions the input point set by octree subdivision and surface normal estimation is performed so as to orientate the local variational implicit patches. A new graph optimization scheme based on the belief propagation framework is proposed to determine the global consistent orientation for the entire set of data points. To achieve multi-scale reconstruction, we propose a novel progressive reconstruction algorithm which utilizes the Schur complement formula to reduce the computational cost of iteratively updating the radial basis function coefficients. Finally, we demonstrate the performance of the proposed algorithm by showing experimental results on some real-world 3D data sets.