cgs algorithm not working with function defined system matrix

%% in following code i am trying to get forward diffraction using a point spread function given by psf
% whole pipeline is given as Au=uin, in which A=psf * original image.*u, u is to be estimated and uin is input .
% A is defined via a function which doe convolution between psf and elemsnt wise multiplication of original image and u.
%i tried to implement it but getting error as "Arrays have incompatile size for this opeartion , can someone hlep me rectify it.


%% this code contains only conjugate gradient based method for image reconstruction
clc
clear all;
close all;
% Load the original image
original_image = im2double(imread('cameraman.tif'));
N=size(original_image,1);
psf = green_f1(N);
% Normalize the PSF
psf = psf / sum(psf(:));
incident_field=ones(N,N);
I=eye(size(original_image));
I=I(:);
% Define the system matrix
A = @(x) convolution2(x.*original_image, psf);
At = @(x) convolution2(x, conj(flipud(fliplr(psf))));
% Define parameters for conjugate gradient method
max_iter = 100; % Maximum number of iterations
tol = 1e-8; % Tolerance for convergence
% Conjugate gradient method
b=incident_field(:);
%b = degraded_image(:);
x0=rand(N,N);
x0=x0(:);
%x0 = estimated_image(:);
[x, ~, relres, iter] = cgs(A, b, tol, max_iter, [], [], x0);
% Reshape the result to image size
estimated_image = reshape(x, size(degraded_image));
% Display the deconvolved image
figure;
imshow(estimated_image); title('Deconvolved Image');
% Helper function for 2D convolution
function y = convolution2(x, h)
    y = ifft2(fft2(x) .* fft2(h, size(x,1), size(x,2)));
end
function G=green_f1(N) ;
% Input parameters
% N=512; % Number of pixels
L=10;  %e-6; % Length of the field of view (in meters)
rad=30e-6; % Radius of aperture (in meters)
lambda=0.630;%e-9; % Wavelength of light (in meters)
z=100e-6; %Propagation distance (in meters)
% Coordinate vectors/grids
dx=L/N; % Length of one pixel
x=-L/2 : dx : L/2-dx; % Coordinate vector (in meters)
fx=-1/(2*dx):1/L:1/(2*dx)-1/L; % Spatial frequency vector (in 1/meters)
[X,Y]=meshgrid(x,x); % Coordinate grids (real space)
[FX,FY]=meshgrid(fx,fx); % Coordinate grids (Fourier space)
% Define aperture
R=sqrt(X.^2+Y.^2); % Radial coordinate
R(R==0) = 1; % Avoid division by zero
%Aperture=double(R<rad); % Aperture
nb=1.351;                % refractive index background
kg=2*pi/lambda;    % wavenumber
%p=parpool(4);
nn=size(R,2);
for i=1:size(R,1)
    %for i=1:size(r,1)
    for j=1:nn
        G(i,j)= (exp(1i*(kg*R(i,j)*nb)))/(4*pi*R(i,j));
    end
end
end

Answer 1

[Animesh]

Hey @asim asrar

The error message "Arrays have incompatible sizes for this operation" typically indicates that there is a size mismatch during operations that involve arrays. So, we need to verify that the sizes of x and original_image are compatible during their element-wise multiplication. Since x is vectorized before being passed to function A, it is necessary to reshape x back to the original image size within function A before performing the multiplication.

Additionally, the provided code snippet appears to have a few more potential issues that need attention:

1.The following line uses degraded_image, which is not defined within the provided context:

estimated_image = reshape(x, size(degraded_image)); 

This line should likely be corrected to ensure that x is reshaped according to the size of the original_image.

estimated_image = reshape(x, size(original_image));  

2. For effective use of the conjugate gradient method it is recommended to use adjoint operator. The function At (used to find the adjoint) is defined in the code but is not being used.

3. The function provided to the cgs method (in this case, a function handle for A) is expected to return a column vector. However, the current implementation of A is returning a 2D array. So, we need to ensure that the output of A is turned back into a column vector before it is returned.

A = @(x) reshape(convolution2(reshape(x,size(original_image)).*original_image, psf), [], 1);  

This modification ensures that x is first reshaped to match the size of the original_image before performing the element-wise multiplication. Subsequently, the result of the convolution operation is reshaped back into a column vector, aligning with the expectations of the cgs method.

You can refer the following MathWorks documentation for more information on cgs function : https://www.mathworks.com/help/releases/R2022b/matlab/ref/cgs.html

[asim asrar]

 % thank you animesh for your kind response



% in following code i am trying to get forward diffraction using a point spread function given by psf
% whole pipeline is given as Au=uin, in which A=psf * original image.*u, u is to be estimated and uin is input .
% A is defined via a function which doe convolution between psf and elemsnt wise multiplication of original image and u.
%i tried to implement it but getting error as "Arrays have incompatile size for this opeartion , can someone hlep me rectify it.


%% this code contains only conjugate gradient based method for image reconstruction
clc
clear all;
close all;
% Load the original image
N=200;
original_image = im2double(imread('cameraman.tif'));
original_image=imresize(original_image,[N,N]);
f1=original_image(:);
% N=size(original_image,1);
psf = green_f1(N);
% Normalize the PSF
%psf = psf / sum(psf(:));
G=psf;
H=G;
incident_field=ones(N,N);
uin=incident_field(:);
% I=eye(size(original_image));
% I=I(:);
% Define the system matrix

[u1,relres,iter]=fwd_prop(f1,G,uin,N);

% Reshape the result to image size
estimated_image = reshape(u1,[N,N]);
% Display the deconvolved image
figure;
subplot 121
imagesc(abs(estimated_image)); title('total field');
colorbar
colormap gray
% Helper function for 2D convolution

vecOp = @(x) x(:);
usc=vecOp(conv2(reshape(f1.*u1,[N,N]),G,'same'));
scattered_field = reshape(usc,[N,N]);


subplot 122
imagesc(abs(scattered_field)); title('scattered field');
colorbar
colormap gray


function G=green_f1(N) 
% Input parameters
%N=256; % Number of pixels
L=40;  %e-6; % Length of the field of view (in meters)
rad=30e-6; % Radius of aperture (in meters)
lambda=0.630;%e-9; % Wavelength of light (in meters)
z=100e-6; %Propagation distance (in meters)
% Coordinate vectors/grids
dx=L/N; % Length of one pixel
x=-L/2 : dx : L/2-dx; % Coordinate vector (in meters)
fx=-1/(2*dx):1/L:1/(2*dx)-1/L; % Spatial frequency vector (in 1/meters)
[X,Y]=meshgrid(x,x); % Coordinate grids (real space)
[FX,FY]=meshgrid(fx,fx); % Coordinate grids (Fourier space)
% Define aperture
R=sqrt(X.^2+Y.^2+(15*lambda)^2); % Radial coordinate
%Aperture=double(R<rad); % Aperture
nb=1.351;                % refractive index background
kg=2*pi/lambda;    % wavenumber
%p=parpool(4);
nn=size(R,2);
for i=1:size(R,1)
    %for i=1:size(r,1)
    for j=1:nn
        G(i,j)= (exp(1i*(kg*R(i,j)*nb)))/(4*pi*R(i,j));
    end
end
%delete(p);
end

%%  version of frward prop using cgs with AT


function [u,relres,iter] = fwd_prop(uin,G,f,N)

A= @(u) A_forward(u,f,G,N);
AT = @(z) A_adjoint(z,f,G,N);

%incident_field=ones(N,N);
max_iter = 500; % Maximum number of iterations
tol = 1e-5; % Tolerance for convergence
% Conjugate gradient method

%b = degraded_image(:);
x0=rand(N,N);
x0=x0(:);

%u = cgs(@(u) AT(A(u)), AT(uin), 'xinit',x0)
[u, ~, relres, iter] = cgs(@(u) AT(A(u)), AT(uin), tol,max_iter, [], [], x0);    % b = uin
%[u, ~, relres, iter] = cgs(A, uin, tol, max_iter, [], [], x0);

end
function z = A_forward(u, f, g,N)
%%% A*u = (I - G*diag(f))*u operator required for forward scattering
%%%
%%% Input:
%%% - u: (Ny x Nx) field
%%% - f: (Ny x Nx) object
%%% - g: 2*(Ny x Nx) Green's function
%%%
%%% Output:
%%% - z: (Ny x Nx) field
%%%
%%% U. S. Kamilov, MERL, 2017.

% assert(all(size(g)==2*size(u)), 'all(size(g)==2*size(u))');
% assert(all(size(f)==size(u)), 'all(size(f)==size(u))');

%z = u - dx*dy*conv2D(f.*u, g);
vecOp = @(x) x(:);

%z = u - dx*dy*conv2D(f.*u, g);
z =u - vecOp(conv2(reshape(f.*u,[N,N]), g,'same'));

end


function u = A_adjoint(z, f, g, N)
%%% Adjoint of A = I - G*diag(f) operator required for forward scattering
%%%
%%% Input:
%%% - z: (Ny x Nx) field
%%% - f: (Ny x Nx) object
%%% - g: 2*(Ny x Nx) Green's function
%%%
%%% Output:
%%% - u: (Ny x Nx) field
%%%
%%% U. S. Kamilov, MERL, 2017.

% assert(all(size(g)==2*size(z)), 'all(size(g)==2*size(z))');
% assert(all(size(f)==size(z)), 'all(size(f)==size(z))');
% z1=size(z)
% g1=size(g)
% f1=size(f)

vecOp = @(x) x(:);

% u = z - dx*dy*conj(f).*conv2DAdj(z, g);
u = z -conj(f).* vecOp(conv2(reshape(z,[N,N]), conj(g),'same'));

end

%% another eversion of forward prop with cgs  without AT
% 
% 
% function [x,relres,iter]=fwd_prop(f,G,uin,N)
% 
% A = @(x) A_forward(x,f,G,N);
% %At = @(x) convolution2(x, conj(flipud(fliplr(psf))));
% max_iter = 500; % Maximum number of iterations
% tol = 1e-5; % Tolerance for convergence
% x0=rand(N,N);
% x0=x0(:);
% %x0 = estimated_image(:);
% [x, ~, relres, iter] = cgs(A, uin, tol, max_iter, [], [], x0);
% 
% end
% 
% function z = A_forward(u, f, g,N)
% % A_forward is calculating Au i.e. (I - G*diag(f)*u) i.e.  u - G*diag(f)*u
% vecOp = @(x) x(:);
% z =u- vecOp(conv2(reshape(f.*u,[N,N]), g,'same'));
% %z =u- f.*vecOp(conv2(reshape(u,[N,N]), g,'same'));
% end
% 


%%
% here i have made two version of A_forward, on usign AT i.e. the adjoint
% operator and other without using adjoint operator. 
% what adwantage is there with usage of adjoint operator in accordance to
% the comment  "For effective use of the conjugate gradient method it is
% recommended to use adjoint operator. The function
%  At (used to find the adjoint) is defined in the code but is not being used"

[asim asrar]

thank you animesh for the kind response. but i hve a dbout regarding usage of AT i.e adjoint operator for

cgs algorithm. can you help me out a bt on it.

[asim asrar]

hey @Animesh can you give your valuable insights on this

[Animesh]

hey @asim asrar, the At function in the code symbolizes the adjoint (transpose) of the system matrix A, which is essential for certain mathematical operations and algorithms, especially in iterative solvers. It plays a crucial role in solving systems involving least square problems or systems where the matrix is not square. Additionally, it helps in improving the convergence of iterative solvers (cgs in this case).