For high-spatial-resolution reconstructions in SPECT imaging, an accurate and fine system matrix (H matrix) should be constructed. In this study, the modified Gaussian interpolation method (GIM) and the modified GIM combined with geometric parameter estimations (GIMGPE) are developed to rapidly construct the complete H matrices of circular-orbit 1-pinhole/4-pinhole SPECT systems. The multi-module multi-resolution (M3R) pinhole SPECT system (in the Center for Gamma-Ray Imaging, University of Arizona) is utilized to verify the usefulness of the proposed interpolation methods. The modified GIM and GIMGPE apply a geometric calibration, a rough grid-scan experiment, and point response parameterization based on Gaussian fitting to construct the full H matrices with a 0.5-mm voxel size and 60 projection angles. The geometric calibration uses a 3-point calibration phantom projected at 60 projection angles to construct the geometric projection model. In the rough grid-scan experiment, a point source is stepped on a regular grid pattern with a 1-mm grid spacing to measure the point response functions (PRFs) at the 0° projection angle. Sequentially, the measured PRFs are parameterized into two-dimensional Gaussians. The PRFs of full H matrices are interpolated by the relations between the Gaussian coefficients and the geometric parameters. The processing time of constructing the complete H matrices with the modified GIM and GIMGPE are respectively reduced by about 347× and 336×, compared to that of the full three-dimensional grid-scan experiment. According to the results of geometric calibration, the magnification of M3R with 1- and 4-pinhole patterns are 1.9× and 2.6×, respectively. The reconstructions of a hot-rod phantom show that the spatial resolutions of M3R with 1- and 4-pinhole patterns are 1.4 mm and 1.0 mm, respectively.