The most known and used phase shifting interferometry (PSI) demodulation methods are one-dimensional temporal linear systems. These methods use the information of the interferogram sequence at a single pixel to recover the modulating phase. Accordingly, scanning all pixels, we obtain the two-dimensional (2-D) modulated phase sought. As PSI demodulation methods do not take into account spatial information, these methods cannot remove unwanted harmonics or noise from the interferogram image space (spatial domain). To remove these unwanted artifacts from the image space, spatial information must be included in the demodulation model. We are going to show that the well-known least-squares system for PSI can be used as a full-field 2-D linear system that uses the temporal and spatial information in conjunction in order to recover the modulating phase while removing noise, unwanted harmonics, and interpolating small empty sections of the image space all in the same process with a low computational time.