clear;
% Define paths to binary videos
clusterVideoPath = 'FilePath\Clusters_Detected.avi';

largeParticlesVideoPath = 'FilePath\LargeParticles_Detected.avi';

videoDir='FilePath';

% Load binary videos
clusterReader = VideoReader(clusterVideoPath);
largeParticleReader = VideoReader(largeParticlesVideoPath);



% Ensure startFrame and endFrame are within the valid range of frames
startFrame = 1; %max(1, startFrame);  % At least the first frame
endFrame = clusterReader.NumFrames; % At most the total number of frames

% Set parameters for tracking
frameRate = 10;  % Assuming 10 frames per second (adjust as needed)
timeBetweenFrames = 1 / frameRate;  % Time between frames in seconds

minObjectSize = 40; % Minimum size to filter out noise
maxDistance = 15;  % Maximum distance to consider object tracked in next frame
areaChangeThreshold = 0.20;  % Area change threshold for splitting/merging
areaBinSize = 100;  % Bin size for area
aspectRatioBinSize = 0.1;  % Bin size for aspect ratio

% Initialize tracks structure with unique IDs
tracks = struct('id', {}, 'centroids', {}, 'areas', {}, 'aspectRatios', {},'frameStart', {}, 'frameEnd', {}, 'active', {}, 'velocities', {}, 'velocitiesX', {}, 'velocitiesY', {}, 'thetaCluster', {});

largeParticleTracks = struct('id', {}, 'centroids', {}, 'velocities', {}, 'velocitiesX', {}, 'velocitiesY', {}, 'thetaLarge', {}, 'frameStart', {}, 'frameEnd', {}, 'active', {});


uniqueID = 1;  % Start unique ID counter
uniqueLargeID = 1;  % Start unique ID counter
objectCount = 0; % Initialize object count

% Arrays for storing data for plots
velocities = [];
areas = [];
aspectRatios = [];
trackDurations = [];

% Sequential tracking loop
figure(1);  % Create a figure for visualization

% Pre-allocate an array for storing valid frame indices and their velocities
validFrameNumbers = [];
validFrameVelocities = [];

% Initialize velocity arrays
velocitiesX = [];
velocitiesY = [];

allFrameVelocities = cell(endFrame, 1); % Store velocities for each frame
frameVelocities = []; % Initialize frame-specific velocity storage

allLargeFrameVelocities = cell(endFrame, 1);



% Structured array to store per-frame cluster properties
framewiseClusterData = struct('FrameNumber', {}, 'ClusterID', {}, 'Centroid', {}, ...
                              'BoundaryFraction', {}, 'BoundaryParticleRatio', {}, ...
                              'NetBoundaryVelocity', {}, 'NetBoundaryAngle', {}, ...
                              'ThetaDiff', {}, 'Skewness', {}, 'ClusterVelocity', {});


% Loop through frames between startFrame and endFrame
for frameNumber = startFrame:endFrame
    % Read pre-loaded binary frame
     % Read binary frame from video
    clusterReader.CurrentTime = (frameNumber - 1) / clusterReader.FrameRate;
    binaryImage = im2bw(readFrame(clusterReader));

    % Remove small objects
    binaryImage = bwareaopen(binaryImage, minObjectSize);

    % Display the current frame number in the command window
    disp(['Processing Frame: ' num2str(frameNumber)]);

    % Get centroids, areas, and bounding boxes for aspect ratio calculation
    stats = regionprops(binaryImage, 'Centroid', 'Area', 'BoundingBox');
    centroids = cat(1, stats.Centroid);
    areasCurrent = [stats.Area];
    aspectRatiosCurrent = arrayfun(@(x) max(x.BoundingBox(3:4)) / min(x.BoundingBox(3:4)), stats);

    % Handle objects tracked from previous frames
    if frameNumber > startFrame
        frameVelocities = [];  % Store velocities for this frame
        
        % Get active tracks (objects being tracked)
        activeTrackIndices = find([tracks.active]);
        
        % Check for matches with centroids in the current frame
        if ~isempty(activeTrackIndices) && ~isempty(centroids)
            previousCentroids = cell2mat(arrayfun(@(x) x.centroids(end, :), tracks(activeTrackIndices), 'UniformOutput', false)');
            distances = pdist2(previousCentroids, centroids);  % Pairwise distances between old and new centroids
            
            % For each tracked object, find the closest object in the current frame
            [minDistances, closestIdx] = min(distances, [], 2);

            % Update tracks with matched centroids
            for i = 1:length(activeTrackIndices)
                trackIdx = activeTrackIndices(i);
                if minDistances(i) <= maxDistance
                    % Check area change
                    previousArea = tracks(trackIdx).areas(end);
                    currentArea = areasCurrent(closestIdx(i));
                    areaChange = abs(currentArea - previousArea) / previousArea;

                    if areaChange <= areaChangeThreshold
                        % Update track with new centroid, area, and aspect ratio
                        tracks(trackIdx).centroids = [tracks(trackIdx).centroids; centroids(closestIdx(i), :)];
                        tracks(trackIdx).areas = [tracks(trackIdx).areas; currentArea];
                        tracks(trackIdx).aspectRatios = [tracks(trackIdx).aspectRatios; aspectRatiosCurrent(closestIdx(i))];
                        tracks(trackIdx).frameEnd = frameNumber;

                        % Calculate velocity
                        if size(tracks(trackIdx).centroids, 1) > 1
                            displacement = tracks(trackIdx).centroids(end, :) - tracks(trackIdx).centroids(end-1, :);
                            velocityMagnitude = sqrt(sum(displacement.^2)) / timeBetweenFrames;
                            velocityX = displacement(1) / timeBetweenFrames;
                            velocityY = displacement(2) / timeBetweenFrames;

                            % Store frame-specific velocity for this cluster
                            frameVelocities = [frameVelocities; velocityMagnitude];
                            tracks(trackIdx).velocities = [tracks(trackIdx).velocities; velocityMagnitude];

                            % Store velocity components
                            velocities = [velocities; velocityMagnitude];
                            velocitiesX = [velocitiesX; velocityX];
                            velocitiesY = [velocitiesY; velocityY];

                            % [OPTIONAL] Also store velocity components in the track structure
                            tracks(trackIdx).velocitiesX = [tracks(trackIdx).velocitiesX; velocityX];
                            tracks(trackIdx).velocitiesY = [tracks(trackIdx).velocitiesY; velocityY];

                            areas = [areas; currentArea];
                            aspectRatios = [aspectRatios; aspectRatiosCurrent(closestIdx(i))];
                        end

                        % Mark centroid as processed
                        centroids(closestIdx(i), :) = NaN;
                        areasCurrent(closestIdx(i)) = NaN;
                    else
                        % Keep track inactive but retain its history
                        tracks(trackIdx).active = false;
                    end
                else
                    % Track is inactive, but retain its history
                    tracks(trackIdx).active = false;
                end
            end
        end
    end

    % Initialize new tracks for unmatched objects
    newObjectIndices = find(~isnan(centroids(:, 1)));  % Objects that are new or unmatched
    for i = 1:length(newObjectIndices)
        idx = newObjectIndices(i);
        tracks(end+1) = struct('id', uniqueID, 'centroids', centroids(idx, :), ...
                       'areas', areasCurrent(idx), 'aspectRatios', aspectRatiosCurrent(idx), ...
                       'frameStart', frameNumber, 'frameEnd', frameNumber, 'active', true, ...
                       'velocities', [], 'velocitiesX', [], 'velocitiesY', [], 'thetaCluster', []);

        uniqueID = uniqueID + 1;  % Increment unique ID for the next new object
    end

    % Store frame-specific velocity data
    allFrameVelocities{frameNumber} = frameVelocities;

    % Store average velocity for valid frames
    if frameNumber > startFrame && ~isempty(frameVelocities)
        validFrameNumbers = [validFrameNumbers, frameNumber];  % Store valid frame number
        validFrameVelocities = [validFrameVelocities, mean(frameVelocities)];  % Store the average velocity
    end

    % Visualization of active tracks
    imshow(binaryImage);  % Update the image data for current frame
    hold on;
    for j = 1:length(tracks)
        if tracks(j).active && ~isempty(tracks(j).centroids)  % Only visualize active tracks
            validCentroids = tracks(j).centroids(~isnan(tracks(j).centroids(:, 1)), :);
            if size(validCentroids, 1) > 1
                plot(validCentroids(:,1), validCentroids(:,2), '-', 'MarkerSize', 3, 'LineWidth', 1.5);
            end
        end
    end
    hold off;
    title(['Frame ' num2str(frameNumber)]);
    drawnow;
    pause(0.1);  % Optional pause to control visualization speed
end


% % Save the tracking data to a .mat file
% save(fullfile(videoDir, 'trackedData.mat'), 'tracks', 'velocities', 'areas', 'aspectRatios');
% disp('Tracking data saved successfully.');


maxDistance1=5;


%% 🟢 Loop 2: Large Particle Tracking (Fixed)
% Initialize large particle tracking structure
largeParticleTracks = struct('id', {}, 'centroids', {}, 'velocities', {}, ...
                             'velocitiesX', {}, 'velocitiesY', {}, 'frameStart', {}, ...
                             'frameEnd', {}, 'active', {});

uniqueLargeID = 1;

% Arrays to store velocities
largeVelocities = [];
largeVelocitiesX = [];
largeVelocitiesY = [];
validLargeFrameNumbers = [];
validLargeFrameVelocities = [];

% Per-frame velocity storage
allLargeFrameVelocities = cell(endFrame, 1);  % Store velocities for each frame

for frameNumber = startFrame:endFrame
    largeParticleReader.CurrentTime = (frameNumber - 1) / largeParticleReader.FrameRate;
    largeBinaryImage = im2bw(readFrame(largeParticleReader));

    disp(['Processing Large Particle Tracking: Frame ' num2str(frameNumber)]);

    stats1 = regionprops(largeBinaryImage, 'Centroid');

    % Print number of detected particles
    disp(['Frame ' num2str(frameNumber) ': Found ' num2str(length(stats1)) ' large particles']);

    if isempty(stats1)
        % If no particles detected, mark all active tracks with NaN velocities
        for i = 1:length(largeParticleTracks)
            if largeParticleTracks(i).active
                largeParticleTracks(i).velocities = [largeParticleTracks(i).velocities; NaN];
                largeParticleTracks(i).velocitiesX = [largeParticleTracks(i).velocitiesX; NaN];
                largeParticleTracks(i).velocitiesY = [largeParticleTracks(i).velocitiesY; NaN];
            end
        end
        continue; % Skip to next frame
    end

    largeCentroids = cat(1, stats1.Centroid);
    usedIndices = false(size(largeCentroids, 1));  % Track which centroids are assigned

    frameLargeVelocities = [];  % Store velocities per frame

    if frameNumber > startFrame
        activeLargeTrackIndices = find([largeParticleTracks.active]);

        if ~isempty(activeLargeTrackIndices) && ~isempty(largeCentroids)
            previousLargeCentroids = cell2mat(arrayfun(@(x) x.centroids(end, :), largeParticleTracks(activeLargeTrackIndices), 'UniformOutput', false)');
            distances = pdist2(previousLargeCentroids, largeCentroids);

            % Iterate over all active tracks
            for i = 1:length(activeLargeTrackIndices)
                trackIdx = activeLargeTrackIndices(i);
                [minDistance, closestIdx] = min(distances(i, :));

                if minDistance <= 20 && ~usedIndices(closestIdx) % Ensure no duplicate assignment
                    usedIndices(closestIdx) = true;  % Mark centroid as used
                    previousCentroid = largeParticleTracks(trackIdx).centroids(end, :);
                    newCentroid = largeCentroids(closestIdx, :);

                    % Compute velocity
                    displacement = newCentroid - previousCentroid;
                    velocityMagnitude = sqrt(sum(displacement.^2)) / timeBetweenFrames;
                    velocityX = displacement(1) / timeBetweenFrames;
                    velocityY = displacement(2) / timeBetweenFrames;

                    % Store velocity per particle per frame
                    frameLargeVelocities = [frameLargeVelocities; velocityMagnitude];

                    % Update track with velocity
                    largeParticleTracks(trackIdx).centroids = [largeParticleTracks(trackIdx).centroids; newCentroid];
                    % Ensure large particle centroids are stored correctly
%largeParticleTracks(trackIdx).centroids = [largeParticleTracks(trackIdx).centroids; newCentroid];

                    largeParticleTracks(trackIdx).velocities = [largeParticleTracks(trackIdx).velocities; velocityMagnitude];
                    largeParticleTracks(trackIdx).velocitiesX = [largeParticleTracks(trackIdx).velocitiesX; velocityX];
                    largeParticleTracks(trackIdx).velocitiesY = [largeParticleTracks(trackIdx).velocitiesY; velocityY];
                    largeParticleTracks(trackIdx).frameEnd = frameNumber;
                    largeParticleTracks(trackIdx).active = true;
                else
                    % If no valid match, store NaN velocity
                    largeParticleTracks(trackIdx).velocities = [largeParticleTracks(trackIdx).velocities; NaN];
                    largeParticleTracks(trackIdx).velocitiesX = [largeParticleTracks(trackIdx).velocitiesX; NaN];
                    largeParticleTracks(trackIdx).velocitiesY = [largeParticleTracks(trackIdx).velocitiesY; NaN];
                     % Track is inactive, but retain its history
                    largeParticleTracks(trackIdx).active = false;
                end
            end
        end
    end

    % Initialize new tracks for unmatched large particles
    newObjectIndices = find(~isnan(largeCentroids(:, 1)) & ~usedIndices);

    for i = 1:length(newObjectIndices)
        idx = newObjectIndices(i);
        if idx > size(largeCentroids, 1) || idx < 1
            continue;
        end
        largeParticleTracks(end+1) = struct('id', uniqueLargeID, 'centroids', largeCentroids(idx, :), ...
                                            'velocities', [], 'velocitiesX', [], 'velocitiesY', [], ...
                                            'frameStart', frameNumber, 'frameEnd', frameNumber, 'active', true);
        uniqueLargeID = uniqueLargeID + 1;
    end

    % Store per-frame velocity data
    allLargeFrameVelocities{frameNumber} = frameLargeVelocities;

    % Store average velocity for valid frames
    if frameNumber > startFrame && ~isempty(frameLargeVelocities)
        validLargeFrameNumbers = [validLargeFrameNumbers, frameNumber];
        validLargeFrameVelocities = [validLargeFrameVelocities, mean(frameLargeVelocities)];
    end

    % Visualization
    imshow(largeBinaryImage);
    hold on;
    for j = 1:length(largeParticleTracks)
        if largeParticleTracks(j).active && ~isempty(largeParticleTracks(j).centroids)
            validCentroids = largeParticleTracks(j).centroids(all(~isnan(largeParticleTracks(j).centroids), 2), :);
            if size(validCentroids, 1) > 1
                plot(validCentroids(:,1), validCentroids(:,2), '-', 'MarkerSize', 3, 'LineWidth', 1.5);
            end
        end
    end
    hold off;
    title(['Large Particle Tracking - Frame ' num2str(frameNumber)]);
    drawnow;
    pause(0.1);
end



%%% Loop 3: Cluster-Large Particle Association & Cluster-Level Calculations

% Initialize global arrays to store cluster-level metrics
allBoundaryFractions      = [];
allBoundaryParticleRatios = [];
allNetBoundaryVelocities  = [];
allNetBoundaryAngles      = [];
allThetaDiffs             = [];
allSkewnessValues         = [];
allClusterVelocities      = [];  % Store per-cluster velocities

for frameNumber = startFrame:endFrame
    disp(['Processing Cluster-Large Particle Association: Frame ' num2str(frameNumber)]);
    
    %% (A) Read current frames for clusters and large particles
    clusterReader.CurrentTime = (frameNumber - 1) / clusterReader.FrameRate;
    binaryImage = im2bw(readFrame(clusterReader));
    
    largeParticleReader.CurrentTime = (frameNumber - 1) / largeParticleReader.FrameRate;
    largeBinaryImage = im2bw(readFrame(largeParticleReader));
    
    %% (B) Label clusters so that each has a unique identifier
    labeledClusters = bwlabel(binaryImage);
    
    %% (C) Extract large particle centroids from the large-particle frame
    largeStats = regionprops(largeBinaryImage, 'Centroid');
    if ~isempty(largeStats)
        largeCentroids = cat(1, largeStats.Centroid);
    else
        largeCentroids = [];
    end
    
    % Associate each large particle with a cluster label
    numLarge = size(largeCentroids, 1);
    clusterLabelsForLarge = zeros(numLarge, 1);
    for p = 1:numLarge
        xIdx = round(largeCentroids(p,1));
        yIdx = round(largeCentroids(p,2));
        if xIdx >= 1 && xIdx <= size(labeledClusters,2) && yIdx >= 1 && yIdx <= size(labeledClusters,1)
            clusterLabelsForLarge(p) = labeledClusters(yIdx, xIdx);
        else
            clusterLabelsForLarge(p) = 0;
        end
    end
    
    %% (D) For each cluster in the current frame, compute cluster-level metrics.
    uniqueLabels = unique(labeledClusters);
    uniqueLabels(uniqueLabels == 0) = [];  % Remove background label
    
    % Loop over each unique cluster label
    for L = uniqueLabels'
        % For each track in 'tracks', check if its current centroid belongs to cluster L.
        for clusterIdx = 1:length(tracks)
            % Only process tracks active in this frame
            if frameNumber < tracks(clusterIdx).frameStart || frameNumber > tracks(clusterIdx).frameEnd
                continue;
            end
            velocityIdx = frameNumber - tracks(clusterIdx).frameStart + 1;
            if velocityIdx > 0 && velocityIdx <= size(tracks(clusterIdx).centroids,1)
                clusterCentroid = tracks(clusterIdx).centroids(velocityIdx, :);
                % Get the cluster label for this track from the labeled image
                xIdx = round(clusterCentroid(1));
                yIdx = round(clusterCentroid(2));
                if xIdx >= 1 && xIdx <= size(labeledClusters,2) && yIdx >= 1 && yIdx <= size(labeledClusters,1)
                    trackClusterLabel = labeledClusters(yIdx, xIdx);
                else
                    trackClusterLabel = 0;
                end
            else
                clusterCentroid = [NaN, NaN];
                trackClusterLabel = 0;
            end
            
            % Process this track only if it belongs to cluster L
            if trackClusterLabel ~= L
                continue;
            end
            
            % Retrieve the cluster's velocity (and its x & y components, if available)
            if velocityIdx > 0 && velocityIdx <= length(tracks(clusterIdx).velocities)
                clusterVelocity = tracks(clusterIdx).velocities(velocityIdx);
            else
                clusterVelocity = NaN;
            end
            if velocityIdx > 0 && velocityIdx <= length(tracks(clusterIdx).velocitiesX)
                clusterVelocityX = tracks(clusterIdx).velocitiesX(velocityIdx);
            else
                clusterVelocityX = NaN;
            end
            if velocityIdx > 0 && velocityIdx <= length(tracks(clusterIdx).velocitiesY)
                clusterVelocityY = tracks(clusterIdx).velocitiesY(velocityIdx);
            else
                clusterVelocityY = NaN;
            end
            
            %% Compute the cluster-specific boundary
            currentClusterMask = (labeledClusters == L);
            clusterBoundaryImage = bwperim(currentClusterMask);
            [Yboundary, Xboundary] = ind2sub(size(clusterBoundaryImage), find(clusterBoundaryImage));
            totalBoundaryPixels = length(Xboundary);
            
            %% Find large particles that belong to this cluster
            indicesThisCluster = find(clusterLabelsForLarge == L);
            largeCentroidsForCluster = largeCentroids(indicesThisCluster, :);
            
            %% Determine which of these large particles are at the boundary
            isBoundaryParticleLocal = false(size(largeCentroidsForCluster,1),1);
            for q = 1:size(largeCentroidsForCluster,1)
                dx = Xboundary - largeCentroidsForCluster(q,1);
                dy = Yboundary - largeCentroidsForCluster(q,2);
                distBoundary = sqrt(dx.^2 + dy.^2);
                if ~isempty(distBoundary) && min(distBoundary) <= 8  % 5 pixel threshold
                    isBoundaryParticleLocal(q) = true;
                end
            end
            
            % Count the number of boundary large particles for this cluster
            nBoundaryLarge = sum(isBoundaryParticleLocal);
            % Calculate boundary occupancy fraction using the alternative method:
            % Assume each boundary large particle occupies 10 pixels.
            occupancyFraction = (nBoundaryLarge * 10) / totalBoundaryPixels;
            
            % Calculate boundary particle ratio: boundary particles relative to total large particles in the cluster
            if ~isempty(largeCentroidsForCluster)
                boundaryParticleRatio = nBoundaryLarge / length(largeCentroidsForCluster);
            else
                boundaryParticleRatio = NaN;
            end
            
            %% Compute net velocity of boundary large particles for this cluster
            netBoundaryVx = 0;
            netBoundaryVy = 0;
            validBoundaryParticles = 0;
            boundaryAngles = [];
            % For each large particle (in this cluster) that is at the boundary,
            % use robust matching with stored large particle tracks (from Loop 2)
            for q = 1:length(indicesThisCluster)
                if ~isBoundaryParticleLocal(q)
                    continue;
                end
                % current large particle centroid for this cluster
                currentLargeCentroid = largeCentroidsForCluster(q, :);
                
                % Perform robust matching: gather active large particle tracks for the current frame
                activeLargeTrackIndices = [];
                lastCentroids = [];
                for lp = 1:length(largeParticleTracks)
                    if frameNumber >= largeParticleTracks(lp).frameStart && frameNumber <= largeParticleTracks(lp).frameEnd && largeParticleTracks(lp).active
                        activeLargeTrackIndices = [activeLargeTrackIndices, lp];
                        lastCentroids = [lastCentroids; largeParticleTracks(lp).centroids(end, :)];
                    end
                end
                
                particleVelocityX = NaN;
                particleVelocityY = NaN;
                if ~isempty(activeLargeTrackIndices)
                    dists = sqrt(sum((lastCentroids - currentLargeCentroid).^2, 2));
                    [minD, bestIdx] = min(dists);
                    matchingThreshold = 20;  % Adjust as needed
                    if minD <= matchingThreshold
                        trackIdxLP = activeLargeTrackIndices(bestIdx);
                        velocityIdxLP = frameNumber - largeParticleTracks(trackIdxLP).frameStart + 1;
                        if velocityIdxLP > 0 && velocityIdxLP <= length(largeParticleTracks(trackIdxLP).velocitiesX)
                            particleVelocityX = largeParticleTracks(trackIdxLP).velocitiesX(velocityIdxLP);
                            particleVelocityY = largeParticleTracks(trackIdxLP).velocitiesY(velocityIdxLP);
                        end
                    end
                end
                
                if ~isnan(particleVelocityX) && ~isnan(particleVelocityY)
                    netBoundaryVx = netBoundaryVx + particleVelocityX;
                    netBoundaryVy = netBoundaryVy + particleVelocityY;
                    validBoundaryParticles = validBoundaryParticles + 1;
                end
                
                % For skewness: compute the angle from the cluster centroid to this large particle
                dx = currentLargeCentroid(1) - clusterCentroid(1);
                dy = currentLargeCentroid(2) - clusterCentroid(2);
                angleParticle = atan2(dy, dx);
                boundaryAngles = [boundaryAngles; angleParticle];
            end
            
            if validBoundaryParticles > 0
                netBoundaryVelocity = sqrt(netBoundaryVx^2 + netBoundaryVy^2) / validBoundaryParticles;
                netBoundaryAngle = atan2(netBoundaryVy, netBoundaryVx);
            else
                netBoundaryVelocity = NaN;
                netBoundaryAngle = NaN;
            end
            
            %% Compute alignment: cosine of the difference between cluster velocity direction and net boundary velocity direction
            if ~isnan(clusterVelocityX) && ~isnan(clusterVelocityY) && ~isnan(netBoundaryAngle)
                thetaCluster = atan2(clusterVelocityY, clusterVelocityX);
                thetaDiff = cos(mod(thetaCluster - netBoundaryAngle + pi, 2*pi) - pi);
            else
                thetaDiff = NaN;
            end
            
            %% Compute skewness of boundary particle angles for this cluster
            if numel(boundaryAngles) > 1
                clusterSkewness = abs(skewness(boundaryAngles));
            else
                clusterSkewness = NaN;
            end
            
            %% Store the computed metrics for this cluster
            allBoundaryFractions      = [allBoundaryFractions; occupancyFraction];
            allBoundaryParticleRatios = [allBoundaryParticleRatios; boundaryParticleRatio];
            allNetBoundaryVelocities  = [allNetBoundaryVelocities; netBoundaryVelocity];
            allNetBoundaryAngles      = [allNetBoundaryAngles; netBoundaryAngle];
            allThetaDiffs             = [allThetaDiffs; thetaDiff];
            allSkewnessValues         = [allSkewnessValues; clusterSkewness];
            allClusterVelocities      = [allClusterVelocities; clusterVelocity];

            % Store the computed metrics for this cluster in the new structured dataset
framewiseClusterData(end+1) = struct('FrameNumber', frameNumber, ...
    'ClusterID', L, ...
    'Centroid', clusterCentroid, ...
    'BoundaryFraction', occupancyFraction, ...
    'BoundaryParticleRatio', boundaryParticleRatio, ...
    'NetBoundaryVelocity', netBoundaryVelocity, ...
    'NetBoundaryAngle', netBoundaryAngle, ...
    'ThetaDiff', thetaDiff, ...
    'Skewness', clusterSkewness, ...
    'ClusterVelocity', clusterVelocity);


        end % end loop over tracks
    end % end loop over unique cluster labels
end % end frame loop




%% 🟢 Save Data for Plots
%% Save Updated Data
save(fullfile(videoDir, 'trackedData.mat'), 'tracks', 'largeParticleTracks', ...
    'allBoundaryFractions', 'allBoundaryParticleRatios', 'allNetBoundaryVelocities', ...
    'allNetBoundaryAngles', 'allThetaDiffs', 'allClusterVelocities', 'allSkewnessValues', ...
    'framewiseClusterData');  % 🟢 NEW: Store framewise cluster data

disp('Tracking and association data saved successfully.');

%% 🟢 Create Plots Directory
plotsDir = fullfile(videoDir, 'Plots');
if ~exist(plotsDir, 'dir')
    mkdir(plotsDir);
end

%% Note:
% Instead of replacing NaNs with -9999 for export, we now filter out invalid data for plotting.
% We assume that valid values are > -1000.

%% 1️⃣ Average Velocity by Frame (Unbinned & Binned)
figure;
% Filter valid frame velocities (if needed)
validIdxFrames = ~isnan(validFrameNumbers) & ~isnan(validFrameVelocities);
plot(validFrameNumbers(validIdxFrames), validFrameVelocities(validIdxFrames), '-o');
xlabel('Frame Number'); ylabel('Average Velocity');
title('Average Velocity vs Frame Number');
saveas(gcf, fullfile(plotsDir, 'avg_velocity_vs_frame.jpg'));
writematrix([validFrameNumbers(validIdxFrames)', validFrameVelocities(validIdxFrames)'], fullfile(plotsDir, 'avg_velocity_vs_frame.xlsx'));

% Binned (50 frames, NaN-safe mean)
binSize = 50;
binnedFrameNumbers = [];
binnedVelocities = [];
for i = 1:binSize:length(validFrameNumbers)
    endIdx = min(i + binSize - 1, length(validFrameNumbers));
    frameRange = validFrameNumbers(i:endIdx);
    frameVelocitiesSubset = validFrameVelocities(i:endIdx);
    validVelocities = frameVelocitiesSubset(~isnan(frameVelocitiesSubset));
    if ~isempty(validVelocities)
        binnedFrameNumbers = [binnedFrameNumbers; mean(frameRange)];
        binnedVelocities = [binnedVelocities; mean(validVelocities)];
    end
end

figure;
plot(binnedFrameNumbers, binnedVelocities, '-o');
xlabel('Binned Frame Number'); ylabel('Binned Average Velocity');
title('Binned Average Velocity vs Frame Number');
saveas(gcf, fullfile(plotsDir, 'binned_avg_velocity_vs_frame.jpg'));
writematrix([binnedFrameNumbers, binnedVelocities], fullfile(plotsDir, 'binned_avg_velocity_vs_frame.xlsx'));