[Python] Trouble implementing a super resolution algorithm with a frequency method

I'm currently trying to implement the following algorithm for image super resolution https://ieeexplore.ieee.org/document/56062 I have trouble in the last step : after having calculated the values of F^c (equation (10b)) which represent the values of the Fourier transform of the continuous image, I do not understand how to put them correctly in a matrix in order to get the HR image by calculating its inverse fft. Here is my code :

import numpy as np
import cv2
import matplotlib.pyplot as plt
import math


#Définition des constantes
N = 128 #Résolution de 128 pixels pour les g_k
p = 4 #Nombre d'images floues à disposition
shifts = [(0,0), (0,30), (30,0), (30,30)] #Shifts entre les images
x = [0, 0, 30, 30]
y = [0, 30, 0, 30]
T = 1 #Période d'échnatillonnage (égale selon x et y)
L = 1 #L défini tel que F(x,y) = 0 lorsque |x| > LT et |y| > LT
centre_HR = L*N #centre de l'image HR_fourier

def load_images(): #charges les images dans le dossier, en appliquant les shifts voulus
images = []
for t in shifts:
x,y = t
img = cv2.imread('jojo.jpg', cv2.IMREAD_GRAYSCALE)
images.append(cv2.resize(img[xx+226), yy+226)], (128,128)))

for k in range(len(images)):
cv2.imwrite('g_%d.jpg' % k, images[k])
return images

def coef_i(r):
return ((r-1) % (2*L) ) - L

def coef_l(r):
return math.floor((r-1)/(2*L)) - L

def calcule_Phi(m,n): #Calcule la matrice Phi
tab = [
[(1/(T*T)) * np.exp(1j*2*np.pi*(x[k]*( (m + N*coef_i(r))/(N*T) ) + y[k]*( (n+N*coef_l(r))/(N*T) )))
for r in range (1, 4*L*L +1)] for k in range (p)]
return np.matrix(tab)

def calcule_colonne_G(m,n): #Calcule la colonne associée à G
tab = [[G[k][m][n]] for k in range(p)]
return np.matrix(tab)

def calcule_solution(m,n): #Résout le problème matriciel avec un algo least squares
return np.linalg.lstsq(calcule_Phi(m,n), calcule_colonne_G(m,n), rcond=None)

def initialise_sol():
sol = [[0 + 0j for i in range(2*L*N)] for j in range(2*L*N)]
return np.array(sol)

if __name__ == "__main__":
noir_et_blanc = cv2.imread('jojo.jpg', cv2.IMREAD_GRAYSCALE)
cv2.imwrite('noir_et_blanc.jpg', noir_et_blanc)

HR_fourier = initialise_sol() #The continuous Fourier image
g = load_images()
G = [np.fft.fft2(im) for im in g]

for m in range(N):
for n in range(N):
X, residus, rang, singular_values = calcule_solution(m,n)
for r in range(1, 4*L*L +1):
HR_fourier[centre_HR + m + N*(coef_i(r))][centre_HR + n + N*(coef_l(r))] = X[r-1]

#To do : the F_c values are calculated : do the inverse fft to get the HR image
#res = np.fft.ifft2(HR_fourier)
#res = np.abs(res)
#cv2.imwrite('HR.jpg', res)


I also do not understand how to correctly chose the values L and T (I chose L_x = L_y and T_x = T_y to simplify). L is defined in the article by the equation (6), under the assumption that F_c is bandlimited, but how to choose such a value in practice ? For T, it is defined as the period of sampling, but I also don't understand how to choose it in practice (is it linked with the upscale factor wanted ?).

I am a beginner in this topic, sorry for the potential mistakes and thanks in advance for your help

Continue reading...