We demonstrate highly accurate recovery of weak gravitational lensing shear using an implementation of the Bayesian Fourier Domain ( BFD ) method proposed by Bernstein & Armstrong ( 2014 , BA14 ) , extended to correct for selection biases . The BFD formalism is rigorously correct for Nyquist-sampled , background-limited , uncrowded image of background galaxies . BFD does not assign shapes to galaxies , instead compressing the pixel data D into a vector of moments M , such that we have an analytic expression for the probability P ( \mbox { \boldmath$M$ } | \mbox { \boldmath$g$ } ) of obtaining the observations with gravitational lensing distortion g along the line of sight . We implement an algorithm for conducting BFD ’ s integrations over the population of unlensed source galaxies which measures \approx 10 galaxies/second/core with good scaling properties . Initial tests of this code on \approx 10 ^ { 9 } simulated lensed galaxy images recover the simulated shear to a fractional accuracy of m = ( 2.1 \pm 0.4 ) \times 10 ^ { -3 } , substantially more accurate than has been demonstrated previously for any generally applicable method . Deep sky exposures generate a sufficiently accurate approximation to the noiseless , unlensed galaxy population distribution assumed as input to BFD . Potential extensions of the method include simultaneous measurement of magnification and shear ; multiple-exposure , multi-band observations ; and joint inference of photometric redshifts and lensing tomography .