Seismic samples are typically designed on a perfect Cartesian grid. However, field constructions can disrupt the sampling geometry, resulting in the samples missing or off-the-grid. Our research goals are to simultaneously regularize off-the-grid samples and interpolate missing data for 3D seismic data under the framework of compressive sensing, which combines a 3D curvelet transform, a fast iterative threshold algorithm, and a merging sampling operator. The new sampling operator combines a binary mask with a barycentric Lagrangian operator for simultaneous interpolation and regularization. The fast iterative threshold algorithm is helpful to improve the interpolation accuracy and efficiency while solving the ill-posed problem. Finally, we demonstrate the effectiveness of the proposed approach by simulated and field datasets.