Measurement of gamma-ray polarization is of interest in many areas of modern physics, both in fundamental and applied research. One interesting case is detection of polarized gamma radiation from positron annihilation. When the positronium decays to two gamma photons, they are created with orthogonal polarizations. It is possible to use coincidence measurements of the two photons to deduce azimuthal correlations via Compton scattering and estimate their initial polarization correlation. This information is of great interest in gamma imaging systems, such as Positron Emission Tomography, where it may be used as an additional handle to distinguish true coincidence events. In this study we used Monte Carlo simulations based on the Geant4 toolkit to model several multi-pixel detector configurations based on LYSO:Ce and GaGG:Ce scintillators. Each module consists of 64 crystals, set up in a single 8x8 matrix, where both the recoil electron and the Compton scattered photon are absorbed. We simulated positron annihilation by generating two back-to-back gamma photons of 511 keV with orthogonal polarizations, which were then detected by two detector modules. The results indicate that a pair of such detectors can observe the modulation of the relative azimuthal scattering angles. The characteristic X-rays following electron knock-out contribute to crosstalk between the crystals, hence we studied their influence on energy resolution for selected detector geometries and materials. The resulting polarimetric performance of LYSO:Ce and GaGG:Ce based configurations are compared and their application in experiments is discussed.