Single-photon avalanche diodes (SPADs) are highly sensitive photodetectors that enable LIDAR imaging at extremely low photon flux levels. While conventional image formation methods require hundreds or thousands of photon detections per pixel to suppress noise, a recent computational approach achieves comparable results when forming reflectivity and depth images from on the order of 1 photon detection per pixel. This method uses the statistics underlying photon detections, along with the assumption that depth and reflectivity are spatially correlated in natural scenes, to perform noise censoring and regularized maximum-likelihood estimation. We expand on this research by adapting the method for use with SPAD arrays, accounting for the spatial non-uniformity of imaging parameters and the effects of crosstalk. We develop statistical models that incorporate these non-idealities, and present a statistical method for censoring crosstalk detections. We show results that demonstrate the performance of our method on simulated data with a range of imaging parameters.