  Fitting parameterized models to images of galaxies has become the standard
for measuring galaxy morphology. This forward modelling technique allows one to
account for the PSF to effectively study semi-resolved galaxies. However, using
a specific parameterization for a galaxy's surface brightness profile can bias
measurements if it is not an accurate representation. Furthermore, it can be
difficult to assess systematic errors in parameterized profiles. To overcome
these issues we employ the Multi-Gaussian expansion (MGE) method of
representing a galaxy's profile together with a Bayesian framework for fitting
images. MGE flexibly represents a galaxy's profile using a series of Gaussians.
We introduce a novel Bayesian inference approach which uses pre-rendered
Gaussian components, which greatly speeds up computation time and makes it
feasible to run the fitting code on large samples of galaxies. We demonstrate
our method with a series of validation tests. By injecting galaxies, with
properties similar to those observed at $z\sim1.5$, into deep HST observations
we show that it can accurately recover total fluxes and effective radii of
realistic galaxies. Additionally we use degraded images of local galaxies to
show that our method can recover realistic galaxy surface brightness and color
profiles. Our implementation is available in an open source python package
$\texttt{imcascade}$, which contains all methods needed for the preparation of
images, fitting and analysis of results.
