I will discuss a new, open-source astronomical image-fitting program, specialized for galaxies, which is fast, flexible, and highly extensible. A key characteristic is an object-oriented design which allows new types of image components (2D surface-brightness functions) to be easily written and added to the program. Image functions provided with the program include the usual suspects for galaxy decompositions (Sersic, exponential, Gaussian), along with Core-Sersic and broken-exponential profiles, elliptical rings, and components which perform line-of-sight integration through 3D luminosity-density models of disks and rings seen at arbitrary inclinations. Minimization can be done using the standard chi^2 statistic (using either data or model values to estimate per-pixel errors) or the Cash statistic, which is appropriate for Poisson data in low-count regimes; different minimization algorithms allow trade-offs between speed and decreased sensitivity to local minima in the fit landscape. I will also show that fitting low-S/N galaxy images by minimizing chi^2 can lead to significant biases in fitted parameter values, which are avoided if the Cash statistic is used; this is true even when Gaussian read noise is present.