The purpose of this study was to develop and implement an accurate and computationally efficient method for determination of the mesh-domain system matrix including attenuation compensation for Ordered Subsets Expectation Maximization (OSEM) Single Photon Emission Computed Tomography (SPECT). The mesh-domain system matrix elements were estimated by first partitioning the object domain into strips parallel to detector face and with width not exceeding the size of a detector unit. This was followed by approximating the integration over the strip/mesh-element union. This approximation is product of: (i) strip width, (ii) intersection length of a ray central to strip with a mesh element, and (iii) the response and expansion function evaluated at midpoint of the intersection length. Reconstruction was performed using OSEM without regularization and with exact knowledge of the attenuation map. The method was evaluated using synthetic SPECT data generated using SIMIND Monte Carlo simulation software. Comparative quantitative and qualitative analysis included: bias, variance, standard deviation and line-profiles within three different regions of interest. We found that no more than two divisions per detector bin were needed for good quality reconstructed images when using a high resolution mesh.