%% LOAD DATA
clear all;  %close all;
N = 64;
 
if N ==32
    Power = 10;
else
    Power = 12;
end
 
file_path = "/Users/jon/Documents/PROJECTS/11_compressedSensing/CS_data/CS_data05.xlsx";
M = readmatrix(file_path,'sheet',6); % Read data from Excel file
 
file_path = "/Users/jon/Documents/PROJECTS/11_compressedSensing/CS_data/CS_data07.xlsx";
M = readmatrix(file_path,'sheet',4); % Read data from Excel file
 
% Convert to column vector
M = M(:); % Flatten in case data is multi-column
 
% Store the first value as the background
background = M(1);
 
% Extract the remaining values
M = M(2:end);
 
plotPatterns = 0;
H = hadamard(N);
 
M = M-background;
M = M./max(M(:));

%% Create Hadamard patterns
H_ordered = zeros(N,N);
for i = 1:N %Column
    counter = 0;
    temp = H(1,i);
    for j = 2:N %Row
        if temp ~= H(j,i)
            counter = counter + 1;
            temp = H(j,i);
        end
    end
    for k = 1:N
        H_ordered(k,counter+1) = H(k,i);
    end
end
 
%Plot the hadamard basis functions
L = zeros(N,N,N^2);  % Preallocate matrix
p=1;
 
for i = 1:N
    for j = 1:N
        temp = H_ordered(i,:)' * H_ordered(j,:);
        temp(temp==-1) = 0;
        L(:,:,p) = temp;
        p = p+1;
    end
end
 
if plotPatterns == 1
    figure;
    for p=1:N^2
        subplot(N,N,p)
        imagesc(L(:,:,p)); axis image;
        title(['Pattern ',num2str(p)]);
        % pause(0.5);
        
    end
end
 
%%
indStart =1;
 
PowerArr = linspace(Power,Power-3,4);
 
for p=1:4
    patternsUsed(p)= 2^PowerArr(p); %length(M); %length(M); % 2^11; % 1024; % 2^11
    percentPatternsUsed(p) = patternsUsed(p)/N^2;
    
    L2=L(:,:,indStart:patternsUsed(p));
    M2 = M(indStart:patternsUsed(p));
    
    %%% Reconstruction
    R = zeros(N,N);
    for m =1:N
        for n=1:N
            temp2 = squeeze(L2(m,n,:));
            R(m,n) = sum(M2.*temp2);
        end
    end
    
    R = imrotate(R,-90);
    R = fliplr(R);
    
    Rs(:,:,p) = flattenImage(R);
    Rs2(:,:,p) = imgaussfilt(Rs(:,:,p),0.5);
    
end
 
f=figure;
f.Position = [100 100 1000 600];
for p=1:4
    
    subplot(2,2,p);
    imagesc(Rs2(:,:,p)); axis image; colorbar; colormap('gray');
    title(['Patterns used ',num2str(patternsUsed(p))]);
    
end
 
f=figure('Color', 'k', 'Position', [100, 100, 2000, 900]); % Set figure size and background color
set(gcf, 'InvertHardCopy', 'off');
set(gcf,'Color',[0 0 0]);
 
imagesc(Rs2(:,:,1),[0.05 1]); axis image; colormap('gray');
ax = gca;
set(ax, 'XTick', [], 'YTick', []);
set(ax, 'XColor', 'w', 'YColor', 'w');
box on;
 
%%%%%%%%%%%%%%%
%% FUNCTIONS %%
%%%%%%%%%%%%%%%

function img_stretch = flattenImage(img)

% Flatten the image into a 1D array and sort it
sortedVals = sort(img(:));
numPixels = numel(sortedVals);
 
% Compute approximate percentiles
p1_idx = round(0.01 * numPixels);   % 1st percentile index
p99_idx = round(0.99 * numPixels);  % 99th percentile index
 
low = sortedVals(max(p1_idx, 1));    % Get 1st percentile value
high = sortedVals(min(p99_idx, numPixels)); % Get 99th percentile value
 
% Normalize and clip
img_stretch = (img - low) / (high - low);
img_stretch = max(0, min(1, img_stretch)); % Clip to range [0,1]

 
end