How to use nufft or nufftn to compute NonUniformFFT of an image with missing (not sampled) pixels?

Given an image x for which I have non sampled values i.e. some pixels are randomly (or given a pattern but lets say randomly) sampled, i.e. I have a nonuniformly sampled 2D signal. I would like to compute its NUFFT.

In the Matlab doc I see this example:

t = [1:10 11:2:29]';
x = t;
y = t';
z = reshape(t,[1 1 20]);
X = cos(2*pi*0.01*x) + sin(2*pi*0.02*y) + cos(2*pi*0.03*z);
Y = nufftn(X,{t,t,t});

Unfortunately this does not help me. The X is a 3D matrix of size 20x20x20 with just the nonuniformly sampled values however I would like to keep the dimenstions of my image e.g. 256x256. The idea I had was to put zeros in the non-sampled pixels and the greyscale value otherwise and give in t the coordinates of the sampled pixels but this does not work. I could just take the nonuniformly sampled pixels of my image and form a new matrix but its shape would be arbitrary, I could even just concatenate all the pixels in a long vector... I have no idea how to do it and how nufftn understands this.

So my question is: how to obtain the NUFFT (2D) of an undersampled 2D image? I would like to understand what happends with the reshaping/what shapes are understood/used by NUFFT... all this seems very obscure to me.

Answer 1

[Chris Turnes]

I'm not sure what format your data is in; however, for the purposes of showing how you might do this, suppose you had a list of pixel locations (stored as linear indices) called sampledLocations for which you have corresponding values. Then to compute the 2-D non-uniform FFT of your non-uniform image samples onto a uniform 2-D grid you would do:

% Just to get some data.
X = randn(256, 256);
% Suppose you pick 100 random pixels to sample.
sampledLocations = randperm(256*256, 100)';
% Convert these to x-y locations.
[i, j] = ind2sub([256 256], sampledLocations);
% Specify the uniformly spaced output frequencies you want.
funiform = (0:255) / 256;
% Non-uniform input points specified as M x 2 vector.
% Uniform output grid specified via cells.
Y = nufftn(X(sampledLocations), [i, j], {funiform, funiform});
% Reshape to your grid size.
Y = reshape(Y, [256 256]);

[Michael Thompson]

IMHO, the more "powerful"/complicated a function is the more important it is to verify its function before use. My attempts to verify the behavior of the nufftn function haven't gone swimmingly so I wanted to ask:

Has the nufftn been tested with fully sampled data sets, and shown to be equivalent to the fft functions?

Is there a location where Matlab provides access to verification code or results?

There is some sample code here and on the nufftn help page, but it's limited to pretty artificial scenarios; my opinion.

Answer 2

[Simon Großmann]

So I hope that I can help you with this issue which I had myself when trying to understand the concept of the 2D nufftn. After experimenting I found a way to get a basic 2D example - also with the inverse nufftn - out of the official 3D example that you stated in combination with the the answer from here (let's call it "frequency example"). It's a 1D example but could be expanded to 2D.

This reply/answer doesn't directly target your specific problem but in both of these examples there are also "missing" points.

In the following will be 3 code blocks: Block 1 is the "frquency example" with using nufftn instead of nufft, Block 2 is the expansion of the "frequency example" to 2D nufftn and Block 3 is the example you stated / which is on the website for 2D.

In each case also the inverse nufftn is applied (using nufftn with a minus in the points).

Remarks:

- The comment regarding the nyquist frequency is from the original "frequency example", it kinda makes sense to choose the frequencys like that, since it's symmetrical around zero but idk the details about that.

- The first choice of n=701 since the last sample was taken at 700.5s also stems from the original "frequency example", just a higher number worked great for me.

I hope I could help you in some way.

Code Block 1: "Frequency example" with nufftn:

% https://connections.mathworks.com/matlabcentral/answers/2032659-adjoint-inverse-nufft#toggle-comments

t = [0:300 500.5:700.5];

S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t);
X = abs(S + rand(size(t)));

n=701;      % The last sample taken at 700.5 sec

f=(-ceil((n-1)/2):floor((n-1)/2))/n; % Nyquist frequency = 0.5 Hz
% NUFFT
Y_NUFFT = nufftn(X(:),t(:),f(:));
S_NUFFT=Y_NUFFT/length(X);  % Power Spectrum
% EDFT

figure(1) %Plot Power Spectrums
plot(f,abs(S_NUFFT)), hold on
X0 = real(nufftn(Y_NUFFT(:),f(:),-t(:)))/n; % reconstructed signal NUFFT
figure(2),plot(t,X), hold on, plot(t,X0), hold off

% Conclusion:
% The nufftn function results in the same result as nufft and also the
% inversion works the same. It's important that when entering a vector, one
% has to put (:) at the vectors so Matlab understands what is meant.

Code Block 2: "Frequency example" in 2D with nufftn:

%% Test for Dimension 2: expand the upper example to 2D
disp("Start test 2D frequency example")

Start test 2D frequency example

t = [0:300 500.5:700.5];
S = 2*sin(2*pi*0.02*t) + sin(2*pi*0.1*t');
X = abs(S);% + rand(size(t)));

n=2001;      % The last sample taken at 700.5 sec -> Larger n, 
            % more detail in frequency spectrum, better reconstruction

f=(-ceil((n-1)/2):floor((n-1)/2))/n; % Nyquist frequency = 0.5 Hz
% NUFFT
Y_NUFFT = nufftn(X(:,:),{t(:),t(:)},{f(:),f(:)});
S_NUFFT=Y_NUFFT/length(X);  % Power Spectrum

% Plotting
S_NUFFT_matrix = real(reshape(S_NUFFT,length(f),length(f)));
[X1,Y1] = meshgrid(f,f);
figure(1) %Plot Power Spectrums
contourf(X1,Y1,abs(S_NUFFT_matrix))
colorbar

X0 = real(nufftn(Y_NUFFT(:,:),{f(:),f(:)},{-t(:),-t(:)}))/n^2; % reconstructed signal NUFFT
X0_matrix = reshape(X0, length(t),length(t));
[X2,Y2] = meshgrid(t,t);
figure(2)
contourf(X2,Y2,abs(X))
colorbar

figure(3)
contourf(X2,Y2,abs(X0_matrix))
colorbar

disp("Errornorm of nufftn-inverseNufftn: " + norm(X0_matrix-X));

Errornorm of nufftn-inverseNufftn: 1.2389

% Conclusion:
% Here figure(2) and figure(3) look alike and therefore the inversion
% seemingly works.
% Remark: The approximation is better when using more sampling points in
% the frequency spectrum, meaning larger n.

Code Block 2: Documentation/your example in 2D with nufftn:

%% Now for 2D nufftn with a different function and sampling
% example from https://de.mathworks.com/help/matlab/ref/double.nufftn.html

disp("Start test 2D Documentation example")

Start test 2D Documentation example

t = [1:10 11:2:29]';
x = t;
z = t;
Z = cos(2*pi*0.01*x) + cos(2*pi*0.03*z');


n=1001;      % just a high sampling number is better

f=(-ceil((n-1)/2):floor((n-1)/2))/n; % Nyquist frequency = 0.5 Hz
                                    % that was the reason in 1D, works here too
% NUFFT
Y_NUFFT = nufftn(Z(:,:),{t(:),t(:)},{f(:),f(:)});
S_NUFFT=Y_NUFFT/length(Z);  % Power Spectrum


% Plotting
S_NUFFT_matrix = real(reshape(S_NUFFT,length(f),length(f)));
[X1,Y1] = meshgrid(f,f);
figure(1) %Plot Power Spectrums
contourf(X1,Y1,abs(S_NUFFT_matrix))
colorbar

Z0 = real(nufftn(Y_NUFFT(:,:),{f(:),f(:)},{-t(:),-t(:)}))/n^2; % reconstructed signal NUFFT
Z0_matrix = reshape(Z0, length(t),length(t));
[X2,Y2] = meshgrid(t,t);
figure(2)
contourf(X2,Y2,abs(Z))
colorbar

figure(3)
contourf(X2,Y2,abs(Z0_matrix))
colorbar

disp("Errornorm of nufftn-inverseNufftn: " + norm(Z0_matrix-Z));

Errornorm of nufftn-inverseNufftn: 1.2187e-12