Realistic image simulation is useful understanding artifacts introduced by lens aberrations. Simple simulations which convolve the system point spread function (PSF) with a scene are typically fast because Fast Fourier transform (FFT) techniques are used to calculate the convolution in the Fourier domain. This technique, however, inherently assumes that the PSF is shift invariant. In general, optical systems have a shift variant PSF and the speed of FFT is lost. To simulate shift variant cases, the scene is often broken down into a set of sub-regions over which the PSF is approximately isoplanatic. The FFT methods can then be employed over each sub-regions and then the sub-regions are recombined to create the image simulation. There is an obvious tradeoff between the number of sub-regions and the fidelity of the final image. This fidelity is dependent upon how quickly the PSF changes between adjacent sub-regions. Here, a different strategy is employed where PSFs at different points in the field are sampled and decomposed into a novel set of basis functions. The PSF at locations between the sampled points are estimated by interpolating the expansion coefficients of the decomposition. In this manner, the image simulation is built up by combining the interpolated PSFs across the scene. The technique is verified by generating a grid of PSFs in raytracing software and determining the interpolated PSFs at various points in the field. These interpolated PSFs are compared to the PSFs calculated directly for the same field point.