Generate K lines in a degraded image

Generate $K$ lines with parameters $(\theta_k,\alpha_k,\eta_k)$ where:

$\theta_k$ is the angle of the line
$\alpha_k$ is the amplitude of the line
$\eta_k$ is the horizontal offset of the line

Lines are blurred with a horizontal filter $g$ and a vertical one $h$ . The image of size $W\times H$ degraded by some noiselevel is called $y$ .

global h hadj hfourextend Hs M Gfour S;

S=ceil(spread*4)-1;            % half-length of the discrete filter h
Hs=H+2*S;                      % height after adding S pixels beyond border
M=(W-1)/2;                     % half-width of the image
h=exp(-((-S:S)'/spread).^2/2); % the Gaussian filter h
h=h/sum(h);                    % filter normalization
vstar=zeros((H+2*S),W);        % image of the lines horizontally blurred
g=h;                           % same blur g horizontally and h vertically

Compute $v^{\star}$ thanks to these two relations:

$v^\sharp(t_1,t_2)=\sum_{k=1}^K \frac{\alpha_k}{\cos(\theta_k)} \varphi_1\Big(t_1+\tan(\theta_k)t_2+\frac{\gamma_k}{\cos(\theta_k)}\Big)$

for k=1:K,
    Angle=t_k(k);           % the angle of line with respect to -y axis
    a_k(k)=a_k(k)/max(h);   % after blurring a_k*h is at most equal to 255
    amplitude=a_k(k);       % horizontal offset p_k from the image's center
    offset=p_k(k);          % may varying between -M and M
    for n2=1:(H+2*S),
        t=(0:W-1)-M+tan(Angle)*((H+1)/2+S-n2)-offset;
        vec=zeros(1,W);
        for n=-S:S,         % Dirichlet : sin(pi*(t-n))./sin(pi*(t-n)/W);
            vec=vec+g(n+S+1)*arrayfun(@(x) Dirichlet_Kernel(x,W),t-n);
        end
        vstar(n2,:)=vstar(n2,:)+amplitude/cos(Angle)*vec/W;
    end
end

xstar=conv2(vstar,h,'valid');  % vertical blur convolution
noise=randn(H,W)*noiselevel;   % additionnal Gaussian noise
y=xstar+noise;                 % degraded image

%imwrite(uint8(y),'degraded.tif');

We display the images $v^{\sharp}$ , $x^{\sharp}$ and $y$

if plotData
    f=figure;
    subplot(1,3,1);
    imagesc(vstar); colormap gray; axis image;
    title('$$v^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,3,2);
    imagesc(xstar); colormap gray; axis image;
    title('$$x^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,3,3);
    imagesc(y); colormap gray; axis image;
    title('$$y$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    warning off;
    truesize(f,[200 200]);
    tightfig;
end

We compute the FFT of these filters for a later use

gextend = [zeros(M-S,1) ; g ; zeros(M-S,1)]; % size 2M+1=W cause we need
gfour=fftshift(fft(ifftshift(gextend)));     % W Fourier coefficients ^g[m]

Hm=(Hs-1)/2;
hextend=[zeros(Hm-S,1) ; h ; zeros(Hm-S,1)];
hfourextend=fftshift(fft(ifftshift(hextend)));

Gfour=diag(gfour(M+1:end)); % matrix diag containing the Fourier coefs
hadj=conj(flipud(h));       % filter of the adjoint of convolution by h

At this stage, we have all we need to perform the minimization on the data $y$ which is the degraded image, in order to recover the solution $w^{\sharp}=\mathcal{F}_1 s$ and then the parameters of lines.

The exact solution $w^{\sharp}$ we want to recover is

$\hat w^{\sharp}[m,n_2]=\sum_{k=1}^K \frac{\alpha_k}{\cos \theta_k}e^{j2\pi \left(\tan (\theta_k)n_2+\gamma_k/\cos(\theta_k)\right)m/W}$

% Theoritical formula for wstar
wstar=zeros((H+2*S),W);
for n2=1:(H+2*S)
    for m=-M:M
        wstar(n2,m+M+1)=sum(a_k./cos(t_k).*exp(1i*2*pi*...
            (tan(t_k).*((H+1)/2+S-n2)-p_k).*m/W));
    end
end

We compute other images to compare with the future provided solution

Compare theoritical and empirical $\mathcal{F}_1 v^{\sharp}$

% Theoritical formula for F_1(v)
vfourtheo=zeros(Hs,W);
for m=1:W
    vfourtheo(:,m)=wstar(:,m).*gfour(m);
end

% Compute empirical F_1(v)
vfouremp=zeros((H+2*S),W);
for n2=1:(H+2*S)
    vfouremp(n2,:) = fftshift(fft(ifftshift(vstar(n2,:)))); % horizontal FFT
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(vfourtheo)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$\mathcal{F}_1 v^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(vfouremp)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$\mathcal{F}_1 v^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end

Compare theoritical and empirical $v^{\sharp}$

% Compute theoritical v
vtheo=zeros(Hs,W);
for n2=1:H+2*S
    vtheo(n2,:)=fftshift(ifft(ifftshift(vfourtheo(n2,:))));
end

% Explicit formula for v (for one vertical line at m=offset)
vtheo1=zeros((H+2*S),W);
for n2=1:(H+2*S)
    for m=-M:M
        vtheo1(n2,m+M+1) = 255*exp(-((m-offset)/spread)^2/2);
    end
end

% Compute empirical v
vemp=zeros(Hs,W);
for n2=1:H+2*S
    vemp(n2,:)=fftshift(ifft(ifftshift(vfouremp(n2,:))));
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(vtheo)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$v^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(vemp)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$v^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end

Compare theoritical and empirical $\hat w^{\sharp}$

% Compute empirical wstar
wemp=zeros((H+2*S),W);
for m=1:W
    wemp(:,m)=vfouremp(:,m)./gfour(m);
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(wstar)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$\hat w^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(wemp)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$\hat w^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end

Compare theoritical and empirical $s^{\sharp}$

% Compute theoritical sstar
stheo=zeros((H+2*S),W);
for n2=1:(H+2*S)
    stheo(n2,:) = fftshift(ifft(ifftshift(wstar(n2,:)))); % horizontal FFT
end

% Compute empirical sstar
semp=zeros((H+2*S),W);
for n2=1:(H+2*S)
    semp(n2,:) = fftshift(ifft(ifftshift(wemp(n2,:)))); % horizontal FFT
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(stheo)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$s^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(semp)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$s^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end

Compare theoritical and empirical $\mathcal{F}_1 x^{\sharp}$

% Compute theoritical F_1(x)
xfourtheo=conv2(vfourtheo,h,'valid');

% Compute empirical F_1(x)
xfour=zeros(H,W);
for n2=1:H
    xfour(n2,:)=fftshift(fft(ifftshift(xstar(n2,:)))); % horizontal FFT
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(xfourtheo)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$\mathcal{F}_1 x^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(xfour)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$\mathcal{F}_1 x^{\sharp}$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end

Compare theoritical and empirical $\mathcal{F}_1 y$

noisefour=zeros(H,W);
for n2=1:H
    noisefour(n2,:)=fftshift(fft(ifftshift(noise(n2,:)))); % horizontal FFT
end
yfourtheo=xfourtheo+noisefour;

yfour=zeros(H,W);
for n2=1:H
    yfour(n2,:)=fftshift(fft(ifftshift(y(n2,:)))); % horizontal FFT
end

if plotComp
    f=figure;
    subplot(1,2,1);
    imagesc(abs(yfourtheo)); colormap gray; axis image;
    title('Theoritical');
    xlabel('$$\mathcal{F}_1 y$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    subplot(1,2,2);
    imagesc(abs(yfour)); colormap gray; axis image;
    title('Empirical');
    xlabel('$$\mathcal{F}_1 y$$','Interpreter','latex');
    set(gca,'xtick',[],'ytick',[]);
    truesize(f,[300 300]);
end