Recent advances in materials, devices and fabrication technologies have motivated a strong momentum in developing solid-state sensors that can detect individual photons in space and time. It has been envisioned that such sensors can eventually achieve very high spatial resolutions (e.g., 109 pixels/chip) as well as high frame rates (e.g., 106 frames/sec). In this paper, we present an efficient algorithm to reconstruct images from the massive binary bit-streams generated by these sensors. Based on the concept of alternating direction method of multipliers (ADMM), we transform the computationally intensive optimization problem into a sequence of subproblems, each of which has efficient implementations in the form of polyphase-domain filtering or pixel-wise nonlinear mappings. Moreover, we reformulate the original maximum likelihood estimation as maximum a posterior estimation by introducing a total variation prior. Numerical results demonstrate the strong performance of the proposed method, which achieves several dB's of improvement in PSNR and requires a shorter runtime as compared to standard gradient-based approaches.