linearized bregman algo 1

clear
% Load the original image and blur it
original_img = imread('cameraman.tif');
original_img = im2double(original_img);
g_orig = original_img;

sigma = 3;

filter_size = 7; % Size of the Gaussian filter (e.g., 7x7)
gaussian_filter = fspecial('gaussian', filter_size, sigma);
blurred_img = conv2(original_img, gaussian_filter, 'same');

% Pad the Gaussian filter with zeros to match the size of the image
image_size = size(original_img); % Assuming the image size is 256x256
padded_filter = zeros(image_size);
padded_filter(1:filter_size, 1:filter_size) = gaussian_filter;
% Create the circulant matrix using FFT
C = fft2(padded_filter);

theta = 1; % A small positive number
% Create the first-order difference matrix with Neumann boundary condition
G = eye(image_size) - circshift(eye(image_size), [1, 0]);

P = real(inv(C' * C + theta * (G' * G)));

% Initialize variables
max_iter = 500;
g = zeros(image_size);
delta = 0.1; % You may need to adjust this value
mu = 1e-6; % You may need to adjust this value

[A, D] = dualtree2(original_img, "Level", 4);

% Initialize u with the same structure as A and D
uA = zeros(size(A));
uD = cell(size(D));

for k = 1:numel(D)
    uD{k} = zeros(size(D{k}));
end

stopping_variance = 0.02;

for k = 1:max_iter

    F_u = idualtree2(uA, uD);
    gradient = blurred_img - conv2(F_u, gaussian_filter, 'same');

    % Stopping criterion
    gradient_norm_to_power_of_2 = norm(gradient)^2;
    %fprintf("iter: %d, gradient norm: %e\n", k, gradient_norm_to_power_of_2)
    if gradient_norm_to_power_of_2 <= stopping_variance
        break
    end


    %gradient_mean = mean(gradient(:)); % debugging purposes
    % Update g^(k+1)
    g = g + gradient;
    
    % Update u^(k+1)
    convolved_P_g = real(ifft2(P .* fft2(g)));
    Gaussian_blur_filter_applied = conv2(convolved_P_g, gaussian_filter, 'same');
    [new_uA, new_uD] = dualtree2(Gaussian_blur_filter_applied, "Level", 4);
    %mean_new_uA = mean(new_uA(:)); % debugging purposes

    %fprintf("iter: %d, gradient mean: %e, new_uA mean: %e\n", k, gradient_mean, mean_new_uA) % debugging purposes
    uA = delta * soft_threshold(new_uA, mu);
    for i = 1:numel(uD)
        uD{i} = delta * soft_threshold(new_uD{i}, mu);
    end
end

% Reconstruct the image from the framelet coefficients u
reconstructed_img = idualtree2(uA, uD);
reconstructed_img = reconstructed_img - min(reconstructed_img(:));
reconstructed_img = reconstructed_img / max(reconstructed_img(:));
psnr = calculate_psnr(original_img, reconstructed_img);
% Display the original, blurred, and reconstructed images
figure('Position', [0, 0, 1000, 600]);
subplot(1, 3, 1); imshow(original_img); title('Original Image');
subplot(1, 3, 2); imshow(blurred_img); title('Blurred Image');
subplot(1, 3, 3); imshow(reconstructed_img); title(join(['Reconstructed Image PSNR '  num2str(psnr)  "dB"]));