For simplicity I will show an approach Id' use on 1D signal (A row of real world image).
You will be able to extend it and I will add few remarks on how you can even gain from having 2D data.
The general idea is as sketched in Estimate the Discrete Fourier Series of a Signal with Missing Samples. The trick here is to exploit prior information.
In our case we'll use the smoothness prior of images by utilizing the Total Variation prior.
Problem Formulation
So let't define our problem. We have the data to estimate $ x \in \mathbb{R}^{n} $ and the DFT samples $ y \in \mathbb{C}^{m} $ where $ n > m $. Let $ F \in \mathbb{C}^{m \times n} $ be the partial DFT matrix such that $ F x = y $, namely it transforms the to be estimated data into the given DFT samples.
This allows us to write the optimization problem as following:
$$ \arg \min_{x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
Where $ D \in \mathbb{R}^{\left( n - 1 \right) \times n} $ is the Forward Finite Differences Operator (Numerical approximation of the Derivative).
By intuition, we're looking for a vector $ x $ which its partial DFT is similar to $ y $ and it obeys local smoothness as real world images do. The parameter $ \lambda $ balances between the two.
Solving this is pretty easy using Proximal Gradient Method or ADMM. But even straight forward Sub Gradient Descent Method will do the trick here with the Gradient given by:
$$ \frac{\partial}{\partial x} \frac{1}{2} {\left\| F x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} = {F}^{H} \left( F x - y \right) + \lambda {D}^{T} \operatorname{sign} \left( D x \right) $$
Results
I used the Lena image and chose the 130th row:

I removed 3 random from the DFT of the row and estimated the row using that as a starting point.
The estimation after 1000 iterations:

The RMSE is about the value of 1 pixel Namely it is usually go unnoticed.
MATLAB Code
The MATLAB code is available at my StackExchange Signal Processing Q60119 GitHub Repository.
Pay attention one need to handle the fact the data to be estimated is Real while operations create Complex data. I used abs()
in MATLAB to generate real data from the complex data.
How to Extend to 2D Data
Well, you could just do as above per row of the image. but since you have 2D information you should use that for better results.
The trick is to apply the Total Variation prior in both directions.
By the way, on noisy image, if you don't want to apply any smoothing solve this as Least Squares problem (Just do $ \hat{x} = {\left( {F}^{H} F \right)}^{-1} {F}^{H} y $).