I will present our pipeline for the surface brightness fitting of galaxies using optical and near infrared imaging data from large surveys, which we applied to ~13,000 nearby galaxies with z<0.08 from the Galaxy And Mass Assembly (GAMA) survey. We fit three models to each galaxy in each of our nine wavelength bands with a fully automated Markov-chain Monte Carlo analysis using the Bayesian two-dimensional profile fitting code ProFit. For the first time, we employ ProFits multi-frame fitting functionality, working with data at the pawprint level and fitting all exposures of the same galaxy in the same band simultaneously, thus avoiding point spread function (PSF) uncertainties due to stacking. All preparatory work, including image segmentation, background subtraction, PSF estimation, and obtaining initial guesses, is carried out using the complementary image analysis package ProFound; and we develop additional routines for post-processing, including model selection, extensive quality control and a detailed investigation into systematic uncertainties. The resulting catalogue of robust structural parameters for the stellar components of galaxies (bulges and disks) can be used to study a variety of properties of galaxies and their components such as colours, luminosity functions, mass-size relations and dust attenuation. At the same time, our work contributes to the advancement of image analysis, surface brightness fitting and post-processing routines for quality assurance in the context of automated large-scale bulge-disk decomposition studies. Such advancements are vital to fully exploit the high-quality data of current and upcoming large imaging surveys.