%% STEP 1: Generate synthetic image dataset
clear; close all; clc;

% Parameters
numImages = 50;   % number of images in the dataset
imgSize   = [50 50]; % image dimensions

% Preallocate array
images = zeros([imgSize, numImages]);

% Generate base features: a horizontal stripe + diagonal gradient
[x, y] = meshgrid(1:imgSize(2), 1:imgSize(1));

for k = 1:numImages
    % Feature 1: horizontal stripe
    stripe = double(y > 20 & y < 30);
    
    % Feature 2: diagonal gradient
    gradient = (x + y) / max(x(:) + y(:));
    
    % Random variation in stripe intensity
    stripeIntensity = 0.8 + 0.4*randn(1);
    
    % Random rotation of gradient strength
    gradIntensity = 0.5 + 0.2*randn(1);
    
    % Combine features + Gaussian noise
    images(:,:,k) = stripeIntensity*stripe + gradIntensity*gradient + 0.1*randn(imgSize);
end

%% STEP 2: Visualize some sample images
figure;
for i = 1:6
    subplot(2,3,i);
    imagesc(images(:,:,i));
    axis image off; colormap gray;
    title(sprintf('Image %d', i));
end
sgtitle('Sample Synthetic Images');

%% STEP 3: Reshape images into a 2D data matrix
% Each image is a row, each pixel is a column
X = reshape(images, [], numImages)';  % size: [numImages × numPixels]

%% STEP 4: Perform PCA
[coeff, score, latent, tsquared, explained] = pca(X);

% coeff  : principal component directions (eigenvectors) in pixel space
% score  : projection of each image onto the principal components
% latent : eigenvalues (variance captured by each PC)
% explained : percentage variance explained by each PC

%% STEP 5: Visualize variance explained
figure;
pareto(explained);
xlabel('Principal Component');
ylabel('Variance Explained (%)');
title('PCA Variance Explained');

%% STEP 6: View first few principal components as "eigenimages"
numPCsToShow = 4;
figure;
for i = 1:numPCsToShow
    pcImage = reshape(coeff(:,i), imgSize); % reshape back to image form
    subplot(1,numPCsToShow,i);
    imagesc(pcImage);
    axis image off; colormap gray;
    title(sprintf('PC %d', i));
end
sgtitle('First Principal Components (Eigenimages)');

%% STEP 7: View projection (scores) in PC space
figure;
scatter(score(:,1), score(:,2), 50, 'filled');
xlabel('PC 1 Score');
ylabel('PC 2 Score');
title('Images in Principal Component Space');
grid on;