Current thick detectors used in medical imaging allow recording many attributes, such as the 3D location of interaction within the scintillation crystal and the amount of energy deposited. An efficient way of dealing with these data is by storing them in list-mode (LM). To reconstruct the data, maximum-likelihood expectation-maximization (MLEM) is efficiently applied to the list-mode data, resulting in the list-mode maximum-likelihood expectation-maximization (LMMLEM) reconstruction algorithm. In this work, we consider a PET system consisting of two thick detectors facing each other. PMT outputs are collected for each coincidence event and are used to perform 3D maximum-likelihood (ML) position estimation of location of interaction. The mathematical properties of the ML estimation allow accurate modeling of the detector blur and provide a theoretical framework for the subsequent estimation step, namely the LMM-LEM reconstruction. Indeed, a rigorous statistical model for the detector output can be obtained from calibration data and used in the calculation of the conditional probability density functions for the interaction location estimates. Our implementation of the 3D ML position estimation takes advantage of graphics processing unit (GPU) hardware and permits accurate real-time estimates of position of interaction. The LMMLEM algorithm is then applied to the list of position estimates, and the 3D radiotracer distribution is reconstructed on a voxel grid.