% Define paths to the folders containing masked TIFF stacks
parentFolder = 'PCA-6groups-final/Clustered_type6_results';

% Initialize a 4D array to store the pixel values
pixelValues4D = zeros(64, 682, 682, 218);

% Iterate through each folder
for i = 2:64
    % Folder name
    folderName = sprintf('Clustered_type6_S%03d', i);
    folderPath = fullfile(parentFolder, folderName);
    
    % Initialize a 3D array to store the masked TIFF stack
    maskedStack = zeros(682, 682, 218);
    
    % Read the TIFF stack
    for j = 1:218
        filename = sprintf('clustered_T%04d.tif', j);
        filePath = fullfile(folderPath, filename);
        maskedStack(:,:,j) = imread(filePath);
    end
    
    % Store the masked stack in the 4D array
    pixelValues4D(i,:,:,:) = maskedStack;
end

% Initialize arrays to store result
pixelsToKeep = zeros(64, 682, 682);
bestFitModel = cell(64, 682, 682);
rateOfChange = zeros(64, 682, 682);
pixelTypeArray = zeros(64, 682, 682);
startEDwithRate = zeros(64, 682, 682, 3);

% Counters for model fits
classicCount = 0;
classicLinearFitCount = 0;
classicQuadraticFitCount = 0;
classicSigmoidFitCount = 0;

mysteriousCount = 0;
mysteriousLinearFitCount = 0;
mysteriousQuadraticFitCount = 0;
mysteriousSigmoidFitCount = 0;

dehydrationHighTempCount = 0;
dehydrationHighTempLinearFitCount = 0;
dehydrationHighTempQuadraticFitCount = 0;
dehydrationHighTempSigmoidFitCount = 0;

dehydrationLowTempCount = 0;
dehydrationLowTempLinearFitCount = 0;
dehydrationLowTempQuadraticFitCount = 0;
dehydrationLowTempSigmoidFitCount = 0;

noTypeCount = 0;

% Initialize arrays to store rates of change for each type
classicRates = [];
mysteriousRates = [];
dehydrationHighTempRates = [];
dehydrationLowTempRates = [];
%% Analysing the rate

% Perform pixel analysis on the 4D array
for i = 2:64
    % Iterate through each pixel
    for row = 1:682
        for col = 1:682
            % Get the pixel values over time
             pixelValues = squeeze(pixelValues4D(i,row,col,:));
            
            % Check if pixel has 0 values across all slices
            if all(pixelValues == 0)
                continue;  % Skip this pixel
            end
        
        % Check for the three specific plateaus
        firstPlateau = mean(pixelValues(1:10));
        secondPlateau = mean(pixelValues(81:90));
        thirdPlateau = mean(pixelValues(172:181));
        fourthPlateau = mean(pixelValues(208:217));
        
        % Classify pixel based on plateau relationships
        if firstPlateau > 1000 && firstPlateau <= 4600 && fourthPlateau >= 7600 && fourthPlateau > thirdPlateau 
            pixelType = 'Classic';
            pixelTypeArray(i,row,col) = 1;
            classicCount = classicCount + 1;
        
        % Find the start and end of the period of change
        diffs = diff(pixelValues(181:204));  % Differences between adjacent values in the specified range
        [~, maxChangeIndex] = max(abs(diffs));  % Find the index with the largest change
        
        startOfChange = maxChangeIndex;
        endOfChange = maxChangeIndex;
        
        % Define criteria for the derivative
        derivativeThreshold = 5;  % Adjust this value as needed
        
        % Search for the start of the change
        for ind = maxChangeIndex:-1:2
            derivative = abs(pixelValues(ind+180) - pixelValues(ind+179));
            if derivative < derivativeThreshold || sign(diffs(ind)) ~= sign(diffs(ind-1))
                startOfChange = ind;
                break;
            end
        end

        % Search for the end of the change
        for ind = maxChangeIndex:numel(diffs)-1
            derivative = abs(pixelValues(ind+180) - pixelValues(ind+179));
            if derivative < derivativeThreshold || sign(diffs(ind)) ~= sign(diffs(ind+1))
                endOfChange = ind;
                break;
            end
        end
        
        % Adjust indices to match the original range
        startOfChange = startOfChange + 180;
        endOfChange = endOfChange + 180;
        periodOfChange = startOfChange:endOfChange;
        
        % Fit linear model (y = mx + c)
        linearFit = polyfit(periodOfChange, squeeze(pixelValues(periodOfChange)), 1);
        linearModel = @(x) linearFit(1) * x + linearFit(2);
        linearFitValues = linearModel(periodOfChange);
        linearError = norm(squeeze(pixelValues(periodOfChange)) - linearFitValues);
        
        % Fit quadratic model (y = ax^2 + bx + c)
         if numel(periodOfChange) >= 3
            quadraticFit = polyfit(periodOfChange, squeeze(pixelValues(periodOfChange)), 2);
            quadraticModel = @(x) quadraticFit(1) * x.^2 + quadraticFit(2) * x + quadraticFit(3);
            quadraticFitValues = quadraticModel(periodOfChange);
            quadraticError = norm(squeeze(pixelValues(periodOfChange)) - quadraticFitValues);  
        else
            % Set default values if there are not enough data points
            quadraticFit = NaN;
%           sigmoidFitValues = [];
            quadraticError = NaN;
         end
        
        % Check if there are enough data points
        if numel(periodOfChange) >= 3
            sigmoidModel = @(a, b, c, x) a ./ (1 + exp(-b * (x - c)));
            
            % Reshape pixelValues(periodOfChange) into a 2D array
            pixelValues2D = reshape(pixelValues(periodOfChange), [], numel(periodOfChange));

            % Fit sigmoid model
            sigmoidFit = fit(periodOfChange',pixelValues2D',sigmoidModel);                   % fit(x, y, sigmoidModel);
            sigmoidFitValues = sigmoidFit(periodOfChange);   %x
            sigmoidError = norm(pixelValues2D - sigmoidFitValues);
        else
            % Set default values if there are not enough data points
            sigmoidFit = NaN;
%           sigmoidFitValues = [];
            sigmoidError = NaN;
        end
        
        % Determine the best fit model
        errors = [linearError, quadraticError, sigmoidError];
        [~, bestModelIndex] = min(errors);
        
        startEDwithRate(i,row,col,1) = pixelValues(175);  %the 1st element in the 4th dimension will be the starting ED
        startEDwithRate(i,row,col,2) = bestModelIndex;  %the 2nd element in the 4th dimension will be the best model index
        
        switch bestModelIndex
            case 1
                bestModel = 'Linear';
                rateOfChange(i,row, col) = linearFit(1);
            case 2
                bestModel = 'Quadratic';
                rateOfChange(i,row, col) = quadraticFit(1) * 2;  % Rate of change for quadratic model
            case 3
                bestModel = 'Sigmoid';
                rateOfChange(i,row, col) = sigmoidFit.b;
        end
        
        startEDwithRate(i,row,col,3) = rateOfChange(i,row,col); 
              
        % Store the best fit model
        bestFitModel{i,row, col} = bestModel;
        
        % Update model fit counters
        switch pixelType
            case 'Classic'
                switch bestModel
                    case 'Linear'
                        classicLinearFitCount = classicLinearFitCount + 1;
                    case 'Quadratic'
                        classicQuadraticFitCount = classicQuadraticFitCount + 1;
                    case 'Sigmoid'
                        classicSigmoidFitCount = classicSigmoidFitCount + 1;
                end
            case 'Mysterious Phase'
                switch bestModel
                    case 'Linear'
                        mysteriousLinearFitCount = mysteriousLinearFitCount + 1;
                    case 'Quadratic'
                        mysteriousQuadraticFitCount = mysteriousQuadraticFitCount + 1;
                    case 'Sigmoid'
                        mysteriousSigmoidFitCount = mysteriousSigmoidFitCount + 1;
                end
            case 'Dehydration and Movement at High Temperatures'
                switch bestModel
                    case 'Linear'
                        dehydrationHighTempLinearFitCount = dehydrationHighTempLinearFitCount + 1;
                    case 'Quadratic'
                        dehydrationHighTempQuadraticFitCount = dehydrationHighTempQuadraticFitCount + 1;
                    case 'Sigmoid'
                        dehydrationHighTempSigmoidFitCount = dehydrationHighTempSigmoidFitCount + 1;
                end
            case 'Dehydration & Movement at Low Temperatures and Crystallization'
                switch bestModel
                    case 'Linear'
                        dehydrationLowTempLinearFitCount = dehydrationLowTempLinearFitCount + 1;
                    case 'Quadratic'
                        dehydrationLowTempQuadraticFitCount = dehydrationLowTempQuadraticFitCount + 1;
                    case 'Sigmoid'
                        dehydrationLowTempSigmoidFitCount = dehydrationLowTempSigmoidFitCount + 1;
                end
        end
        
        else
            noTypeCount = noTypeCount + 1;
            continue;  % Skip this pixel if it doesn't fit any category
        
        end
        end
    end
end

% Calculate average rate of change for each type
classicRates = rateOfChange(pixelTypeArray == 1 & classicCount > 0);
mysteriousRates = rateOfChange(pixelTypeArray == 2 & mysteriousCount > 0);
dehydrationHighTempRates = rateOfChange(pixelTypeArray == 3 & dehydrationHighTempCount > 0);
dehydrationLowTempRates = rateOfChange(pixelTypeArray == 4 & dehydrationLowTempCount > 0);

averageClassicRate = mean(classicRates(:));
averageMysteriousRate = mean(mysteriousRates(:));
averageDehydrationHighTempRate = mean(dehydrationHighTempRates(:));
averageDehydrationLowTempRate = mean(dehydrationLowTempRates(:));
%% Save all workspace variables to a .mat file
save('mysterious_growth_rate_63slices.mat','-v7.3');

%% Plot figures
% Extract the 1st and 3rd elements along the 4th dimension
firstElement = startEDwithRate(:,:,:,1);
secondElement = startEDwithRate(:,:,:,2);
thirdElement = startEDwithRate(:,:,:,3);

% Reshape the extracted elements into 1D arrays
firstElement1D = reshape(firstElement, [], 1);
secondElement1D = reshape(secondElement, [], 1);
thirdElement1D = reshape(thirdElement, [], 1);

% Initialize a color array
colors = zeros(length(secondElement1D), 3);

% Assign specific colors based on the value of the 2nd element
colors(secondElement1D == 1, :) = repmat([0, 1, 0], sum(secondElement1D == 1), 1); % Green for value 1
colors(secondElement1D == 2, :) = repmat([1, 1, 0], sum(secondElement1D == 2), 1); % Yellow for value 2
colors(secondElement1D == 3, :) = repmat([1, 0, 0], sum(secondElement1D == 3), 1); % Red for value 3

% % Create a 2D histogram
% figure;
% histogram2(firstElement1D, thirdElement1D);
% title('2D Histogram of 1st vs 3rd Element in 4th Dimension');
% xlabel('1st Element Value');
% ylabel('3rd Element Value');
% zlabel('Frequency');

% Alternatively, create a scatter plot
figure;
scatter(firstElement1D, thirdElement1D, 10, colors, '.');
title('Scatter Plot of starting ED vs ACC growth rate');
xlabel('Starting electron density');
ylabel('ACC growth rate');


%% Display results
disp(['Number of pixels with classic fit: ', num2str(classicCount)]);
disp(['Number of pixels with linear fit in classic type: ', num2str(classicLinearFitCount)]);
disp(['Number of pixels with quadratic fit in classic type: ', num2str(classicQuadraticFitCount)]);
disp(['Number of pixels with sigmoid fit in classic type: ', num2str(classicSigmoidFitCount)]);
disp(['Average rate of change for classic type: ', num2str(averageClassicRate)]);

disp(['Number of pixels with mysterious phase fit: ', num2str(mysteriousCount)]);
disp(['Number of pixels with linear fit in mysterious phase type: ', num2str(mysteriousLinearFitCount)]);
disp(['Number of pixels with quadratic fit in mysterious phase type: ', num2str(mysteriousQuadraticFitCount)]);
disp(['Number of pixels with sigmoid fit in mysterious phase type: ', num2str(mysteriousSigmoidFitCount)]);
disp(['Average rate of change for mysterious phase type: ', num2str(averageMysteriousRate)]);

disp(['Number of pixels with dehydration and movement at high temperatures fit: ', num2str(dehydrationHighTempCount)]);
disp(['Number of pixels with linear fit in dehydration and movement at high temperatures type: ', num2str(dehydrationHighTempLinearFitCount)]);
disp(['Number of pixels with quadratic fit in dehydration and movement at high temperatures type: ', num2str(dehydrationHighTempQuadraticFitCount)]);
disp(['Number of pixels with sigmoid fit in dehydration and movement at high temperatures type: ', num2str(dehydrationHighTempSigmoidFitCount)]);
disp(['Average rate of change for dehydration and movement at high temperatures type: ', num2str(averageDehydrationHighTempRate)]);

disp(['Number of pixels with dehydration & movement at low temperatures and crystallization fit: ', num2str(dehydrationLowTempCount)]);
disp(['Number of pixels with linear fit in dehydration & movement at low temperatures and crystallization type: ', num2str(dehydrationLowTempLinearFitCount)]);
disp(['Number of pixels with quadratic fit in dehydration & movement at low temperatures and crystallization type: ', num2str(dehydrationLowTempQuadraticFitCount)]);
disp(['Number of pixels with sigmoid fit in dehydration & movement at low temperatures and crystallization type: ', num2str(dehydrationLowTempSigmoidFitCount)]);
disp(['Average rate of change for dehydration & movement at low temperatures and crystallization type: ', num2str(averageDehydrationLowTempRate)]);

disp(['Number of pixels with no type fit: ', num2str(noTypeCount)]);

%% % Filter the thirdElement1D values where secondElement1D == 2
thirdElement_filtered = thirdElement1D(secondElement1D == 3);

% Calculate the average of the filtered thirdElement values
average_thirdElement = mean(thirdElement_filtered);

% Display the result
fprintf('The average of thirdElement where secondElement == 3 is: %.4f\n', average_thirdElement);
   
%% generate txt file
% Filter pairs of firstElement and thirdElement based on secondElement
linear_fit = [firstElement1D(secondElement1D == 1), thirdElement1D(secondElement1D == 1)];
% quadratic_fit = [firstElement1D(secondElement1D == 2), thirdElement1D(secondElement1D == 2)];
% sigmoid_fit = [firstElement1D(secondElement1D == 3), thirdElement1D(secondElement1D == 3)];

% % Find the maximum number of rows among the three sets
% max_rows = max([size(linear_fit, 1), size(quadratic_fit, 1), size(sigmoid_fit, 1)]);

% Pad the arrays with NaNs to ensure they all have the same number of rows
% linear_fit_padded = [linear_fit; NaN(max_rows - size(linear_fit, 1), 2)];
% quadratic_fit_padded = [quadratic_fit; NaN(max_rows - size(quadratic_fit, 1), 2)];
% sigmoid_fit_padded = [sigmoid_fit; NaN(max_rows - size(sigmoid_fit, 1), 2)];

% Combine all the data into one array with 6 columns
combined_array = [linear_fit];

% Define the output file name
output_file = 'ACC_crystal_growth_rate.txt';

% Open the file for writing
fid = fopen(output_file, 'w');

% Write a header for the columns
fprintf(fid, 'start ED\trate_for_linear\n');

% Write the data to the file
for i = 1:size(combined_array, 1)
    fprintf(fid, '%f\t%f\n', combined_array(i, 1), combined_array(i, 2));
end

% Close the file
fclose(fid);

% Display confirmation
fprintf('Data successfully written to %s\n', output_file);
%% % Histogram plot with average growth rate versus starting electron density
% Open the file
% filename = 'ACC_mysterious_growth_rate.txt';  % Replace with your actual file name
% fileID = fopen(filename, 'r');
% 
% % Read the data from the file (assuming it's numeric)
% data = textscan(fileID, '%f %f', 'Delimiter', ' ', 'HeaderLines', 1,'EmptyValue',NaN);  % Adjust delimiter as needed
% fclose(fileID);

linear_fit = [firstElement1D(secondElement1D == 1), thirdElement1D(secondElement1D == 1)];

% Extract the first and second columns
first_column = linear_fit(:,1);
second_column = linear_fit(:,2);

% Remove any NaN values (in case there are empty or invalid entries in the file)
% valid_rows = ~isnan(first_column) & ~isnan(second_column);
% first_column = first_column(valid_rows);
% second_column = second_column(valid_rows);

% Define the number of bins
num_bins = 20;

% Create the bin edges based on the range of the first column
bin_edges = linspace(min(first_column), max(first_column), num_bins+1);

% Initialize arrays to store the bin centers and the average of the second column
bin_centers = zeros(num_bins, 1);
bin_averages = zeros(num_bins, 1);

% Loop through each bin
for i = 1:num_bins
    % Find the indices of the data points that fall within the current bin
    in_bin = (first_column >= bin_edges(i)) & (first_column < bin_edges(i+1));
    
    % Calculate the bin center
    bin_centers(i) = (bin_edges(i) + bin_edges(i+1)) / 2;
    
    % Calculate the average of the second column values within this bin
    if sum(in_bin) > 0
        bin_averages(i) = mean(second_column(in_bin));
    else
        bin_averages(i) = NaN;  % If no data points in the bin, assign NaN
    end
end

bin_centers = bin_centers/10000;
bin_averages = bin_averages/10000;

% Plot the histogram (as a bar plot)
figure;
bar(bin_centers, bin_averages, 'FaceColor', [0.2, 0.7, 0.5]);
title('Histogram with Averages of ACC Growth Rate versus Starting Electron Density');
xlabel('Starting Electron density');
ylabel('Average of ACC Growth Rate');
%% Export histogram to a txt file

% Define the output file name
output_file = 'ACC _growth_rate_linear_model_histogram.txt';

% Open the file for writing
fid = fopen(output_file, 'w');

% Write a header for the columns
fprintf(fid, 'start ED\trate_for_linear\n');

% Write the data to the file
for i = 1:size(bin_centers, 1)
    fprintf(fid, '%f\t%f\n', bin_centers(i), bin_averages(i) );
       
end

% Close the file
fclose(fid);

% Display confirmation
fprintf('Data successfully written to %s\n', output_file);

% %%  % Define the logarithmic model function with two input arguments
%     function y_pred = logarithmicModel(params, x)
%         a = params(1);
%         b = params(2);
%         c = params(3);
%         y_pred = a * log(b * x(:)) + c;
%     end



