(a) (b) Figure 1: (a) Gravitational field strength with harmonic coefficients of degree n and order m. (b) Vesta gravitational perturbations due to harmonics up to degree n = 3. This paper presents a pure vision based approach to solving for the gravitational field of extraterrestrial bodies with image data obtained by an orbiting spacecraft or satellite. Recovering a spacecraft's trajectory with modern day Structure from Motion approaches allows for further investigation for perturbations to accelerations due to variation in the strength of gravity. Understanding the variations of these forces, as well as developing a map, help to derive various models on the interior structure of the target planetary body or asteroid [1,4].Classical approaches for recovering the strength of a gravitational field study the motion of a satellite by tracking its position with Earth based telescopes. The basic principle behind this approach was developed in the field of satellite geodesy with the specific goal to define a highly accurate map of Earth's gravitational field. The same principle has not changed significantly, where the use of X-band Doppler and range measurements from a collection of Earth based radiometric tracking stations, known as the Deep Space Network, has been used to great effect.In this paper, we introduce method to recover an estimate of the gravitational field without any need of radiometric tracking. We formulate constraints on a set of spherical harmonic coefficients, which defines a map of gravitational variations on a sphere, as shown in Figure 1, that integrate with graphical models used in modern Structure from Motion techniques [2,3,6]. Our approach is a complete image-based pipeline based around a two-step optimization that recovers 3D structure, spacecraft kinematics, and a gravitational model.The basic process for gravity estimation is a two step iterative optimization. First, spacecraft pose and 3D landmark variables are estimated using batch bundle adjustment. The second step involves optimizing for the parameters of the gravitational field, in addition to camera pose velocities, using the local solutions found in step one. Here, tracking residuals are minimized with respect to global models.Development of two key error terms for recovering the gravitational field are presented. First, a dynamics based gravitational potential function is used to compare the error of a point mass' orbiting trajectory with the recovered camera positions given a set of spherical harmonic coefficients. The spherical harmonics, commonly referred to as Stokes coefficients, define a basis for the gravitational model, similar to a Fourier series but instead map to the surface of a unit sphere. This is the key error term behind recovering the gravitational perturbations. A second error term based upon the Kaula power law constrains the magnitude of the harmonic coefficients as a function of their degree.We evaluated our approach using camera data from the DAWN spacecraft's orbits around 4 Vesta, the second largest astero...