The purpose of this work is to generalize pulse shape modelling to include the effect of light bending due to general relativity. A simple analytical formula is found to describe the light-bending effect accurately, and is incorporated into the emission model for calculation of the pulse shapes from rotating accreting neutron stars. Least-squares fittings are performed to a subset of pulsars studied by Leahy. The fittings with light bending show significant improvement over the fittings without light bending, and the resulting derived emission region geometry gives narrower rings, as expected on theoretical grounds. The model is also applied to the transient pulsar EXO 2030+375. The relation between the luminosity and the properties of the emission rings is studied. For this pulsar we conclude that some effect other than changing emission ring size is responsible for the change in pulse shape.