% Define paths to the folders containing masked TIFF stacks
parentFolder = 'export/Upsampled_results/masked_64_slices';

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

% Iterate through each folder
for i = 2:64
    % Folder name
    folderName = sprintf('Upsampled_S%d', i);
    folderPath = fullfile(parentFolder, folderName);
    
    % Initialize a 3D array to store the masked TIFF stack
    maskedStack = zeros(682, 682, 98);
    
    % Read the TIFF stack
    for j = 1:98
        filename = sprintf('masked_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);

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

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

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

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

noTypeCount = 0;

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

% 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(51:60));
        thirdPlateau = mean(pixelValues(85:95));
        
        % Classify pixel based on plateau relationships
        if thirdPlateau > firstPlateau && thirdPlateau > secondPlateau && firstPlateau > secondPlateau
            pixelType = 'Classic';
            pixelTypeArray(i,row,col) = 1;
            classicCount = classicCount + 1;
        elseif thirdPlateau > secondPlateau && thirdPlateau > firstPlateau && secondPlateau > firstPlateau
            pixelType = 'Mysterious Phase';
            pixelTypeArray(i,row,col) = 2;
            mysteriousCount = mysteriousCount + 1;
        elseif firstPlateau > secondPlateau && secondPlateau > thirdPlateau && firstPlateau > thirdPlateau
            pixelType = 'Dehydration and Movement at High Temperatures';
            pixelTypeArray(i,row,col) = 3;
            dehydrationHighTempCount = dehydrationHighTempCount + 1;
        elseif thirdPlateau > firstPlateau && thirdPlateau > secondPlateau && secondPlateau > firstPlateau
            pixelType = 'Dehydration & Movement at Low Temperatures and Crystallization';
            pixelTypeArray(i,row,col) = 4;
            dehydrationLowTempCount = dehydrationLowTempCount + 1;
        else
            noTypeCount = noTypeCount + 1;
            continue;  % Skip this pixel if it doesn't fit any category
        end
        
        % Find the start and end of the period of change
        diffs = diff(pixelValues(11:50));  % 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+10) - pixelValues(ind+9));
            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+10) - pixelValues(ind+9));
            if derivative < derivativeThreshold || sign(diffs(ind)) ~= sign(diffs(ind+1))
                endOfChange = ind;
                break;
            end
        end
        
        % Adjust indices to match the original range
        startOfChange = startOfChange + 10;
        endOfChange = endOfChange + 10;
        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
       
%         % Fit logarithmic model only if there are enough data points
%         if numel(periodOfChange) >= 2
%             
% %             % Flatten the 3D array into a vector
% %             pixelValuesVector = pixelValues(periodOfChange);
% %             pixelValuesVector = pixelValuesVector(:);
% 
%           % Define an initial guess for the logarithmic model parameters
%             initialGuess = [max(pixelValues(periodOfChange)), 1];
%           % Fit logarithmic model using lsqcurvefit
%             logParams = lsqcurvefit(@logarithmicModel, initialGuess, periodOfChange, pixelValues(periodOfChange));   
%           % Calculate fitted values
%             logFitValues = logarithmicModel(logParams, periodOfChange);
%           % Calculate error
%             logError = norm(pixelValues(periodOfChange) - logFitValues);
%         else
%             % Set default values if there are not enough data points
%             logParams = [];
%             logFitValues = [];
%             logError = 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);
        
        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 = 'Logarithmic';
%                 rateOfChange(row, col) = logParams(1) * logParams(2);
%             case 4
                bestModel = 'Sigmoid';
                rateOfChange(i,row, col) = sigmoidFit.b;
        end
        
        % 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 'Logarithmic'
%                         classicLogFitCount = classicLogFitCount + 1;
                    case 'Sigmoid'
                        classicSigmoidFitCount = classicSigmoidFitCount + 1;
                end
            case 'Mysterious Phase'
                switch bestModel
                    case 'Linear'
                        mysteriousLinearFitCount = mysteriousLinearFitCount + 1;
                    case 'Quadratic'
                        mysteriousQuadraticFitCount = mysteriousQuadraticFitCount + 1;
%                     case 'Logarithmic'
%                         mysteriousLogFitCount = mysteriousLogFitCount + 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 'Logarithmic'
%                         dehydrationHighTempLogFitCount = dehydrationHighTempLogFitCount + 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 'Logarithmic'
%                         dehydrationLowTempLogFitCount = dehydrationLowTempLogFitCount + 1;
                    case 'Sigmoid'
                        dehydrationLowTempSigmoidFitCount = dehydrationLowTempSigmoidFitCount + 1;
                end

            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(:));

% 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 logarithmic fit in classic type: ', num2str(classicLogFitCount)]);
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 logarithmic fit in mysterious phase type: ', num2str(mysteriousLogFitCount)]);
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 logarithmic fit in dehydration and movement at high temperatures type: ', num2str(dehydrationHighTempLogFitCount)]);
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 logarithmic fit in dehydration & movement at low temperatures and crystallization type: ', num2str(dehydrationLowTempLogFitCount)]);
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)]);

   

    %% % Save all workspace variables to a .mat file
save('1st_4d_pixel_analysis_results.mat', 'pixelValues4D', 'bestFitModel', 'rateOfChange', 'pixelTypeArray', 'classicRates', 'mysteriousRates', 'dehydrationHighTempRates', 'dehydrationLowTempRates');;

%%  % 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

