Algorithm that enlarges the image to a resolution of $2N \times 2N$ using DFT operations

Question

I'm trying to solve the following question:

Given an image at a resolution of $N \times N$. Describe an algorithm that enlarges the image to a resolution of $2N \times 2N$ using DFT operations.

As I understand, we need to use the DFT, split it into four sections and add zero padding and then perform IDFT (it was mentioned by my lecturer). But I'm not sure if it's correct or why it's correct. Is it possible to explain the full algorithm (step-by-step) that is required here?

Answer 1

I am copying my answer from Applying 2D Sinc Interpolation in the Fourier Domain (DFT / FFT).

Given a Matrix $ A \in \mathbb{R}^{m \times n} $ in order to interpolate it into a grid of size $ k \times l $ where $ k \geq m $ and $ l \geq n $ one could use the "Separability" property of the DFT and apply it once over the rows and once over the columns.

Things might be trickier for the case of upsampling on one dimension and downsampling on the other. In higher dimensions it even get trickier.

One of the approaches to implement this would be applying DFT along a single dimension and then apply the same interpolation as in the link per slice (Which has N-1 dimensions).

In the case of 2D it is quite simple as each slice is basically a 1D signal which can be done easily as in the link.

I implemented [ mY ] = DFTUpSample2D( mX, vSizeO ) with the following test:

clear();
close('all');

numRowsI = 5000;
numColsI = 5200;

numRowsO = 10000;
numColsO = 10400;

sincRadius = 5;

mX      = GenTest([numRowsI, numColsI], sincRadius);
mYRef   = GenTest([numRowsO, numColsO], sincRadius);

mY = DFTUpSample2D(mX, [numRowsO, numColsO]);

figure();
imshow(mX, []);
figure();
imshow(mY, []);
figure();
imshow(mYRef, []);

max(abs(mYRef - mY), [], 'all')

function [ mX ] = GenTest( vSize, sincRadius )

vX = linspace(-sincRadius, sincRadius, vSize(2) + 1);
vX(end) = [];
vY = linspace(-sincRadius, sincRadius, vSize(1) + 1);
vY = vY(:);
vY(end) = [];

% mX = abs(vX) + abs(vY) + sinc(sqrt(vX .^2 + vY .^2));
mX = sinc(sqrt(vX .^2 + vY .^2));

end

The result is 7.2676e-05 which means the interpolation is valid.

The code is available at my StackExchange Codes Signal Processing Q81493 GitHub Repository (Look at the SignalProcessing\Q81493 folder).