extractImageFeatures.m

function [features,removedIndices] = extractImageFeatures(datastore,optional)

%% Description
% -------------------------------------------------------------------------
% A function that extracts image features from an image datastore. The
% function first resizes all images to the same size. Then, it extracts
% SIFT features and low-level RGB statistic features from each image. The
% images are divided into a grid, and features are extracted from each grid
% cell. Finally, the function applies PCA to reduce the number of features.

% Inputs:
% -------------------------------------------------------------------------
% => datastore: The imageDatastore that contains the raw data.
%
% => optional: A structure that contains optional parameters for the
%              function. It includes the number of grid rows and columns
%              for dividing the image, and the number of features to retain
%              after PCA.

% Outputs:
% -------------------------------------------------------------------------
% => features: A structure array that contains the extracted features for
%              each image. Each element of the array corresponds to one
%              image and contains the SIFT features, RGB statistic
%              features, and the reduced versions of these features after
%              PCA.
%
% Details:
% -------------------------------------------------------------------------
% The function works in several steps. First, it resizes all images to the
% same size. Then, it divides each image into a grid and extracts SIFT
% features from each grid cell. If no SIFT points are found in a cell, the
% function continues with the next cell. The SIFT features from all cells
% are concatenated vertically.
%
% In the next step, the function calculates low-level RGB statistic
% features for each grid cell. These include the mean, standard deviation,
% skewness, and kurtosis of the RGB values. The statistics for all cells
% are reshaped into a 1-by-96 vector.
%
% Finally, the function applies PCA to the SIFT features and RGB statistic
% features separately, and retains a specified number of principal
% components. The reduced features are added to the output structure.


arguments (Input)

    datastore          {mustBeUnderlyingType(datastore, ...
                                             ['matlab.io.datastore.' ...
                                             'ImageDatastore'])}



    optional.gridRows  {mustBePositive,mustBeNonempty,...
                        mustBeReal,mustBeInteger} = 2

    optional.gridCols  {mustBePositive,mustBeNonempty,...
                        mustBeReal,mustBeInteger} = 4

    optional.Feature_reduction {mustBePositive,mustBeNonempty,...
                                mustBeReal,mustBeInteger} = 50
end

%% First step is to resize all the images
i=0;
reset(datastore)
rows = zeros(length(datastore.Files),1);
cols = zeros(length(datastore.Files),1);

while hasdata(datastore)
    i=i+1;
    im = read(datastore);
    [rows(i),cols(i)] = size(im);
end

% Find the minimum dimensions of the images
minrows = min(rows);
mincols = min(cols);

imgSizeThresh = [minrows,mincols];

% Resize all the images based on the minimum dimensions using bilinear
% interpolation
resized_imds = transform(datastore,@(x)imresize(x,imgSizeThresh,"bilinear"));

reset(resized_imds)

%% Extract SIFT features from all the images

totalFiles = length(resized_imds.UnderlyingDatastores{:}.Files);
loadingIcons = ['-', '\', '|', '/']; % Loading icon sequence

h = waitbar(0, 'Extracting SIFT features...'); warning("off")

i=0;
while hasdata(resized_imds)

    i=i+1;
    I = read(resized_imds);
    
    % Divide the image into 2-by-4 grid
    [rows, cols, ~] = size(I);
    gridRows = optional.gridRows;
    gridCols = optional.gridCols;

    regionHeight = floor(rows/gridRows);
    regionWidth = floor(cols/gridCols);

    points = detectSIFTFeatures(im2gray(I));
    Global_SIFT_features = extractFeatures(im2gray(I),points,"Method", ...
                                                                   "SIFT");

    % Find the SIFT features for each region
    region_SIFT_features = [];

    for r = 1:gridRows

        for c = 1:gridCols

            % Extract the region

            region = I((r-1)*regionHeight+1:r*regionHeight, ...
                       (c-1)*regionWidth+1:c*regionWidth, :);

            % Find the SIFT features for the region
            points = detectSIFTFeatures(im2gray(region));
            SIFT_features = extractFeatures(im2gray(region), ...
                points,"Method","SIFT");

            region_SIFT_features = [region_SIFT_features;SIFT_features];

            % Calculate percentage of completion
            percentComplete = i / totalFiles;

            % Update waitbar
            waitbar(percentComplete, h, sprintf(['Extracting regional SIFT ' ...
                    'features %s \n %.2f%% complete'],loadingIcons(mod(i-1, ...
                    numel(loadingIcons)) + 1),percentComplete*100));

        end

    end

    features(i).SIFT_Features = [Global_SIFT_features;
                                 region_SIFT_features];

end

% Close waitbar
close(h)

%% Extract low-level RGB statistic features from all the images

reset(resized_imds)


h = waitbar(0, 'Extracting regional RGB features...'); warning("off")

i=0;
while hasdata(resized_imds)

    i=i+1;
    I = im2double(read(resized_imds));

    % Divide the image into 2-by-4 grid
    [rows, cols, ~] = size(I);
    gridRows = optional.gridRows;
    gridCols = optional.gridCols;

    regionHeight = floor(rows/gridRows);
    regionWidth = floor(cols/gridCols);
    regionStats_reshaped = zeros(8,12);

    % Calculate statistics for each region
    regionIndex = 1;

    for r = 1:gridRows

        for c = 1:gridCols

            % Extract the region
            region = I((r-1)*regionHeight+1:r*regionHeight, (c-1)*regionWidth+1:c*regionWidth);
    
            % Calculate statistics for the region
            meanVal = mean(region, 'all');
            stdVal = std(double(region), 0, 'all');
            skewVal = skewness(double(region(:)));
            kurtVal = kurtosis(double(region(:)));
    
            % Replicate the statistics three times to mimic RGB statistics
            regionStats = repmat([meanVal, stdVal, skewVal, kurtVal], 1, 3);
    
            regionStats_reshaped(regionIndex,:) = regionStats;
    
            % Increment the region index
            regionIndex = regionIndex + 1;
            
        end

    end

    features(i).RGBfeatures = reshape(regionStats_reshaped,[1,96]);

    % Calculate percentage of completion
    percentComplete = i / totalFiles;

    % Update waitbar
    waitbar(percentComplete, h, sprintf(['Extracting RGB features %s \n ' ...
        '%.2f%% complete'],loadingIcons(mod(i-1, numel(loadingIcons)) + 1), ...
        percentComplete*100));
end

% Close waitbar
close(h)
%% Apply PCA to reduce to the optional number of features

for i = 1:length(features)

    [~,features(i).SIFT_scores] = pca(features(i).SIFT_Features);
    [~,features(i).RGB_scores] = pca(features(i).RGBfeatures');
    features(i).SIFT_scores_size = size(features(i).SIFT_scores,2);
    
end

%% Check that all data have at least the optional reduced features. If not delete them

% Due to the nature of the SIFT data, some of them couldn't provide the 
% optional reduced features (here with dimensionality of 50), so they had 
% to be removed. Please pay attention to the difference between the Trainds
% and the Training_features dimensionality!

% Get the field names
fields = fieldnames(features);

% Initialize an empty array to hold the indices of elements to remove
removedIndices = [];

% Loop over each element of features
for i = 1:length(features)

    % Check if SIFT_scores_size is in range [optional.Feature_reduction,128]

    if features(i).SIFT_scores_size < optional.Feature_reduction ...
            || features(i).SIFT_scores_size > 128

        removedIndices = [removedIndices, i];

        continue; % Skip to the next iteration, no need to check other
        % fields for this element
    end

    % Loop over each field of the current element
    for j = 1:length(fields)

        % Check if there are any NaN values in the current field
        if any(isnan(features(i).(fields{j})))

            % If there are, add the index to the removeIndices array
            removedIndices = [removedIndices, i];

            break; % No need to check other fields for this element
        end
    end
end

% Remove the elements with indices in removeIndices from the features array
features(removedIndices) = [];


for i = 1:length(features)

    % Select the minimum principal components as
    features(i).reduced_SIFT_features = features(i).SIFT_scores(:,1: ...
                                               optional.Feature_reduction);

    % Select the first 50 principal components
    features(i).reduced_RGB_features = features(i).RGB_scores(1: ...
                                            optional.Feature_reduction,:)';
end

features = rmfield(features,{'SIFT_scores_size','SIFT_scores','RGB_scores'});
end