%% 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%% 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 %% 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 optimization using cgs, matlab MATLAB Answers — New Questions

​