A calculation of the one-loop contribution to the rare three-body flavor changing neutral current top quark decay $t\to c\gamma\gamma$ is presented in the framework of models with one or more scalar leptoquark $SU(2)$ doublets with hypercharge $7/6$. Analytical expressions for the invariant amplitude of the generic decay $f_i\to f_j\gamma\gamma$, with $f_{i,j}$ a lepton or quark, are presented in terms of Passarino-Veltman integral coefficients, from which the amplitudes for the processes $t\to c\gamma\gamma$ and $\ell_i\to \ell_j\gamma\gamma$ follow easily. An analysis of the current constraints on the parameter space is presented in the scenario with only one scalar LQ doublet and bounds on the LQ couplings are obtained from the muon $g-2$ anomaly, the lepton flavor violating (LFV) decay $\tau\to \mu\gamma$ and extra constraints meant to avoid tension between theory predictions and experimental data. For a LQ with a mass in the range of $1$--$1.5$ TeVs, the estimate ${\rm Br}(t\to c\gamma\gamma)\sim 10^{-11}$--$10^{-12}$ is obtained for the largest allowed values of the LQ coupling constants, which means that this decay would be below the reach of future experimental measurements. We also consider an scenario with three scalar doublets, which was recently proposed to explain the lepton flavor universality violation anomalies in $B$ decays as well as the muon $g-2$ anomaly. Although this scenario allows large LQ couplings to the tau lepton and the $c$ and $t$ quarks, the branching ratio of the $t\to c\gamma\gamma$ decay is also of the order of $10^{-11}$--$10^{-12}$ for LQ masses of 1.7 TeV.