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...