  Strong gravitational lensing is a powerful technique for probing galaxy mass
distributions and for measuring cosmological parameters. We present a pixelated
approach to modeling simultaneously the lens potential and source intensity of
strong gravitational lens systems with extended source-intensity distributions.
For systems with sources of sufficient extent such that the separate lensed
images are connected by intensity measurements, the accuracy in the
reconstructed potential is solely limited by the quality of the data. We apply
this potential reconstruction technique to deep HST observations of B1608+656,
a four-image gravitational lens system formed by a pair of interacting lens
galaxies. We present a comprehensive Bayesian analysis of the system that takes
into account the extended source-intensity distribution, dust extinction, and
the interacting lens galaxies. Our approach allows us to compare various models
of the components of the lens system, which include the point-spread function
(PSF), dust, lens galaxy light, source-intensity distribution, and lens
potential. Using optimal combinations of the PSF, dust, and lens galaxy light
models, we successfully reconstruct both the lens potential and the extended
source-intensity distribution of B1608+656. The resulting reconstruction can be
used as the basis of a measurement of the Hubble constant. We use our
reconstruction of the gravitational potential to study the relative
distribution of mass and light in the lensing galaxies. We find that the
mass-to-light ratio for the primary lens galaxy is (2.0+/-0.2)h M_{\sun}
L_{B,\sun}^{-1} within the Einstein radius 3.9 h^{-1} kpc, in agreement with
what is found for noninteracting lens galaxies at the same scales. (Abridged)
