We present a method for estimating the point spread function of a terahertz imaging system designed to operate in reflection mode. The method is based on imaging phantoms with known geometry, which have patterns with sharp edges at all orientations. The point spread functions are obtained by a deconvolution technique in the Fourier domain. We validate our results by using the estimated point spread functions to deblur several images of natural scenes and by direct comparison with a point source response. The estimations turn out to be robust and produce consistent deblurring quality over the entire depth of the focal region of the imaging system.