Three-dimensional optoacoustic imaging uses detection of laser-generated thermoelastic waves with an ultrasound sensor array. Integrated acoustic signals (velocity potentials) are back projected into the source volume to give a map of absorbed laser energy. Since the number of array elements and the receiving solid angle are limited, radial back projection produces artifacts such as back projection arcs. To solve this problem, we developed in this study an iterative method for image reconstruction. A first image estimate was generated by simple radial back projection. A model for signal generation from a volume containing arbitrary optoacoustic sources was then used to calculate acoustic wave propagation from this estimate. Calculated signals at the array elements were compared with the measured signals and the difference was used to improve the image. In simulations and experiments we used the algorithm to obtain three-dimensional images of multiple optoacoustic sources. With signals from an array of 3 x 3 detector elements a significant improvement was observed after about 10 iterations compared to the simple radial back projection. Although computationally more intensive, iterative reconstruction can minimize the time and instrumentation for signal acquisition because a small number of array elements already gives a good quality optoacoustic image.