Matlab to Python Conversion

Question

I am trying to convert some code from Matlab to Python, but I am unfamiliar with a considerable amount of the Matlab syntax and functionality. I have managed to do a some of the conversion using the PIL and Numpy python package, but I was hoping someone would be able to explain what is going on with some elements of this code.

clear all;close all;clc;

% Set gray scale to 0 for color images. Will need more memory
GRAY_SCALE = 1


% The physical mask placed close to the sensor has 4 harmonics, therefore
% we will have 9 angular samples in the light field
nAngles = 9;
cAngles = (nAngles+1)/2;


% The fundamental frequency of the cosine in the mask in pixels
F1Y = 238;  F1X = 191;  %Cosine Frequency in Pixels from Calibration Image
F12X = floor(F1X/2);
F12Y = floor(F1Y/2);

%PhaseShift due to Mask In-Plane Translation wrt Sensor
phi1 = 300;  phi2 = 150;


%read 2D image
disp('Reading Input Image...');
I = double(imread('InputCones.png'));

if(GRAY_SCALE)
   %take green channel only
   I = I(:,:,2);
end


%make image odd size
I = I(1:end,1:end-1,:);

%find size of image
[m,n,CH] = size(I);


%Compute Spectral Tile Centers, Peak Strengths and Phase
for i = 1:nAngles
    for j = 1:nAngles
        CentY(i,j) = (m+1)/2 + (i-cAngles)*F1Y;
        CentX(i,j) = (n+1)/2 + (j-cAngles)*F1X;
        %Mat(i,j) = exp(-sqrt(-1)*((phi1*pi/180)*(i-cAngles) + (phi2*pi/180)*(j-cAngles)));
    end
end

Mat = ones(nAngles,nAngles);
% 20 is because we cannot have negative values in the mask. So strenght of
% DC component is 20 times that of harmonics
Mat(cAngles,cAngles) = Mat(cAngles,cAngles) * 20;


% Beginning of 4D light field computation
% do for all color channel

for ch = 1:CH

    disp('=================================');
    disp(sprintf('Processing channel %d',ch));

    % Find FFT of image
    disp('Computing FFT of 2D image');
    f = fftshift(fft2(I(:,:,ch)));


    %If you want to visaulize the FFT of input 2D image (Figure 8 of
    %paper), uncomment the next 2 lines
    %     figure;imshow(log10(abs(f)),[]);colormap gray;
    %     title('2D FFT of captured image (Figure 8 of paper). Note the spectral     replicas');


    %Rearrange Tiles of 2D FFT into 4D Planes to obtain FFT of 4D Light-Field
    disp('Rearranging 2D FFT into 4D');
    for i = 1: nAngles
        for j = 1: nAngles
            FFT_LF(:,:,i,j) =  f( CentY(i,j)-F12Y:CentY(i,j)+F12Y, CentX(i,j)-F12X:CentX(i,j)+F12X)/Mat(i,j);
        end
    end
    clear f

    k = sqrt(-1);
    for i = 1:nAngles
        for j = 1:nAngles
            shift = (phi1*pi/180)*(i-cAngles) + (phi2*pi/180)*(j-cAngles);
            FFT_LF(:,:,i,j) = FFT_LF(:,:,i,j)*exp(k*shift);
        end
    end



    disp('Computing inverse 4D FFT');
    LF     =    ifftn(ifftshift(FFT_LF)); %Compute Light-Field by 4D Inverse FFT
    clear FFT_LF

    if(ch==1)
        LF_R = LF;
    elseif(ch==2)
        LF_G = LF;
    elseif(ch==3)
        LF_B = LF;
    end
    clear LF

end
clear I

%Now we have 4D light fiel

disp('Light Field Computed. Done...');
disp('==========================================');




% Digital Refocusing Code
% Take a 2D slice of 4D light field
% For refocusing, we only need the FFT of light field, not the light field

disp('Synthesizing Refocused Images by taking 2D slice of 4D Light Field');

if(GRAY_SCALE)
    FFT_LF_R = fftshift(fftn(LF_R));
    clear LF_R
else
    FFT_LF_R = fftshift(fftn(LF_R));
    clear LF_R

    FFT_LF_G = fftshift(fftn(LF_G));
    clear LF_G

    FFT_LF_B = fftshift(fftn(LF_B));
    clear LF_B
end


% height and width of refocused image
H = size(FFT_LF_R,1);
W = size(FFT_LF_R,2);


count = 0;
for theta = -14:14

    count = count + 1;

    disp('===============================================');
    disp(sprintf('Calculating New ReFocused Image: theta = %d',theta));

    if(GRAY_SCALE)
        RefocusedImage = Refocus2D(FFT_LF_R,[theta,theta]);
    else
        RefocusedImage = zeros(H,W,3);
        RefocusedImage(:,:,1) = Refocus2D(FFT_LF_R,[theta,theta]);
        RefocusedImage(:,:,2) = Refocus2D(FFT_LF_G,[theta,theta]);
        RefocusedImage(:,:,3) = Refocus2D(FFT_LF_B,[theta,theta]);
    end

    str = sprintf('RefocusedImage%03d.png',count);

    %Scale RefocusedImage in [0,255]
    RefocusedImage = RefocusedImage - min(RefocusedImage(:));
    RefocusedImage = 255*RefocusedImage/max(RefocusedImage(:));

    %write as png image
    clear tt
    for ii = 1:CH
        tt(:,:,ii) = fliplr(RefocusedImage(:,:,ii)');
    end
    imwrite(uint8(tt),str);

    disp(sprintf('Refocused image written as %s',str));

end

Here is the Refocus2d function:

function IOut = Refocus2D(FFTLF,theta)


[m,n,p,q] = size(FFTLF);
Theta1 = theta(1);
Theta2 = theta(2);

cTem = floor(size(FFTLF)/2)  +   1;


% find the coordinates of 2D slice
[XX,YY] = meshgrid(1:n,1:m);
cc = (XX - cTem(2))/size(FFTLF,2);
cc = Theta2*cc + cTem(4);
dd = (YY - cTem(1))/size(FFTLF,1);
dd = Theta1*dd + cTem(3);

% Resample 4D light field along the 2D slice
v = interpn(FFTLF,YY,XX,dd,cc,'cubic');

%set nan values to zero
idx = find(isnan(v)==1);
disp(sprintf('Number of Nans in sampling = %d',size(idx,1)))
v(isnan(v)) =   0;

% take inverse 2D FFT to get the image
IOut = real(ifft2(ifftshift(v)));

If anyone could help it would be greatly appreciated.

Thanks in advance

---

Apologies: Here is a brief description of what the code does:

The code reads in an image of a light field, and with prior knowledge of the plenoptic mask it we store the relevant nAngles and the fundamental frequencies of the mask and the phase shift, these are used to find multiple spectral replicas of the image.

Once the image is read in and the green channel is extracted we perform a Fast Fourier Transform on the image, and start taking slices from the image matrix that represent one of the spectral replicas.

We then take the Inverse Fourier Transform of all the spectral replicas to produce a the light field.

The Refocus2d function, then takes a 2 dimensional slice of the 4d data to recreate a refocused image.

The things I am struggling with specifically are:

FFT_LF(:,:,i,j) =  f( CentY(i,j)-F12Y:CentY(i,j)+F12Y, CentX(i,j)-F12X:CentX(i,j)+F12X)/Mat(i,j);

We are taking a slice from the Matrix f, but where is that data in FFT_LF? What does (:,:,i,j) mean? Is it a multidimensional array?

and what does the size function return:

[m,n,p,q] = size(FFTLF);

Just a brief explanation of how this translates to python would be a great help.

Thanks everyone so far :)

Answer 1

How about getting start with this page http://www.scipy.org/NumPy_for_Matlab_Users? Also if you have brief description of what this is supposed to do, that would be good

Answer 2

You're correct: FFT_LF(:,:,i,j) refers to a multidimensional array. In this case, FFT_LF is a 4-D array, but the calculations result in a 2-D array. The (:,:,i,j) tells MATLAB exactly how to place the 2-D results into the 4-D variable.

In effect, it is storing one MxN array for each pair of indices (i,j). The colons (:) effectively mean "get every element in that dimension."

What [m,n,p,q] = size(FFTLF) will do is return the length of each dimension in your array. So, if FFTLF ends up being a 5x5x3x2 array, you get:

 m=5, n=5, p=3, q=2.

If you have MATLAB available, typing "help size" should give a good explanation of what it does. The same can be said for most MATLAB functions: I've always been quite impressed with their documentation.

Hope that Helps