Seismic samples are generally designed to be placed on perfect Cartesian coordinates, that is, on-the-grid. However, sampling geometry is disturbed by obstacles in field applications. Large obstacles result in missing samples. For small obstacles, geophones or sources are placed at an available off-the-grid location nearest to the designed grid. To achieve simultaneous off-the-grid regularization and missing data reconstruction for 3D seismic data, we develop a new mathematical model based on a new combined sampling operator, a 3D curvelet transform, and a fast projection onto convex sets (FPOCS) algorithm. The sampling operator is combined with a binary mask for on-the-grid samples reconstruction and a barycentric Lagrangian (BL) operator for off-the-grid samples regularization. A 2D BL operator is obtained using the tensor product of two 1D BL operators. The inversion problem is efficiently solved based on FPOCS. This method is tested on synthetic and field data sets. The reconstruction results outperform the methods based on the binary mask in terms of signal-to-noise ratio and visual effect.