We present a calculation of the NLO QCD corrections to the loop-induced production of a photon pair through gluon fusion, including massive top quarks at two loops. The two-loop virtual amplitudes are computed via projecting onto a set of Lorentz structures related to linear polarisation states of external gauge bosons. All two-loop master integrals dependent on the top quark mass are evaluated numerically. Matching the fixed order NLO result to a calculation with threshold resummation, we obtain an accurate description of the diphoton invariant mass spectrum around the threshold of top quark pair production. We analyse how the distribution of the diphoton invariant mass is affected by the top quark loop-induced corrections, especially around the top quark pair production threshold.