The intense magnetic fields of neutron stars (NSs) generate strong anisotropy and polarization of radiation emanating from their surfaces. These emission characteristics are sensitive to the emission locales and thus depend on principal geometry parameters, namely the angles between the magnetic axis and the observers viewing direction relative to the spin axis, as well as the sizes of emission regions. In this work, we present results from a Monte Carlo simulation of NS atmospheres that uses a complex electric field vector formalism to treat polarized radiative transfer due to magnetic Thomson scattering. General relativistic and QED effects on light propagation from the stellar surface to a distant observer are considered. We outline a range of theoretical predictions of phase-resolved intensity and polarization pulse profiles at different X-ray energies for magnetars, the most highly magnetized NSs, and NSs with lower magnetizations. By comparing these models with the corresponding observations for different neutron stars, we obtain constraints on the geometry angles and the emission region sizes. Finally, we highlight the critical impacts of magnetospheric vacuum birefringence on polarization properties.