This paper proposes an algorithm of simulating spatially correlated polarimetric synthetic aperture radar (PolSAR) images based on the inverse transform method (ITM). Three flexible non-Gaussian models are employed as the underlying distributions of PolSAR images, including the KummerU, and models. Additionally, the spatial correlation of the texture component is considered, which is described by a parametric model called the anisotropic Gaussian function. In the algorithm, PolSAR images are simulated by multiplying two independent components, the speckle and texture, that are generated separately. There are two main contributions referring to two important aspects of the ITM. First, the inverse cumulative distribution functions of all the considered texture distributions are mathematically derived, including the Fisher, Beta, and inverse Beta models. Second, considering the high computational complexities the implicitly expressed correlation transfer functions of these texture distributions have, we develop an alternative fast scheme for their computation by using piecewise linear functions. The effectiveness of the proposed simulation algorithm is demonstrated with respect to both the probability density function and spatial correlation.