In this paper we present a quasi-Newton projection method for image deblurring. The mathematical problem is a constrained minimization problem, where the objective function is a regularization function and the constraint is the positivity of the solution. The regularization function is a sum of the Kullback-Leibler divergence, used to minimize the error in the presence of Poisson noise, and of a Tikhonov term. The Hessian of the regularization function is approximated in order to invert it using Fast Fourier Transforms. The numerical experiments on some astronomical images blurred by simulated Point Spread Functions show that the method gives very good results both in terms of relative error and computational efficiency.