A three-dimensional simulation of gas puff imaging (GPI) diagnostics is carried out by using a self-consistent multi-component model of the neutral-plasma interaction. The simulation, based on the drift-reduced Braginskii model for the plasma and a kinetic model for the neutrals, is performed in a toroidally limited plasma with gas puff sources located at the low field side equatorial midplane. In addition to electrons, the simulation evolves the turbulent dynamics of D+ and D-2(+) ions as well as D and D-2 neutral species. The D-alpha emission arising from the excitation of D atoms and the contributions from dissociation of D-2 molecules and D-2+ ions are considered. The simulation points out the importance of considering D-alpha emission due to molecular dissociation in the far scrape-off layer (SOL), since it is the dominant source of D-alpha emission at distances from the gas puff considerably smaller than the mean free path of D-2 molecules. The correlation functions between the D-alpha emission rate and the plasma and neutral quantities, namely, the electron density, n(e), electron temperature, T-e, and density of neutral atoms, n(D), are evaluated considering each contribution to D alpha emission and analyzing the correlation functions between these quantities. The correlation functions strongly depend on the location considered within the edge and SOL with an important impact on the interpretation of GPI measurements. The statistical moments and the turbulence properties computed for different components of the D-alpha emission as well as for the relevant plasma and neutral quantities are also investigated. While neglecting neutral density fluctuations is a reasonable approximation that is widely used in the analysis of GPI measurements, this work reveals a 20%-30% influence of neutral fluctuations on most of the quantities measured through the GPI diagnostics with a possibly larger impact for some quantities in specific regions. These results, therefore, suggest the importance of considering neutral fluctuations for the accurate quantitative interpretations of GPI measurements. (C )2022 Author(s)