Traditional methods of point spread function (PSF) modeling of pinhole SPECT systems, such as fitting the PSF with a 2D Gaussian, are generally sufficient in characterizing imaging systems in which the detector is the main source of PSF blur. However, when modeling the PSF of a high-resolution imaging system, these methods are too simplistic and fail to capture important PSF features, resulting in errors in later simulation studies. In this work, we present a method for parameterizing point spread functions that is then used to rapidly generate system matrices of gamma-ray imaging systems using high-resolution detectors with sub-millimeter spatial resolution. Our algorithm, which utilizes a GPU-based ray tracer to simulate a system's true PSF, accounts for the PSF blur due to not only the detector's depth of interaction, but also the penetration through the pinhole and the finite size of of the source. By considering all three blurring sources, we are able to better model the PSF, capturing features like the flat-top in the center and the tilt and drop-off towards the edges. Comparisons to the 2D Gaussian fitting model are examined and we report that as the ratio between the source size and pinhole radius increases, the two models converge to one another.