videoFile = 'FilePath\OrginalMovie.avi'; outputMovieEvents = 'FilePath\movie1.avi'; outputMovieSplitMerge = 'FilePath\movie2.avi'; % Define the directory where the video and data files are located videoDir = 'FilePath'; videoReader = VideoReader(videoFile); % Continue with the rest of your code (initialization, segmentation, etc.) videoReader = VideoReader(videoFile); frameWidth = videoReader.Width; frameHeight = videoReader.Height; frameRate = videoReader.FrameRate; frameNumber = 1; duration = videoReader.Duration; % Get the duration of the video in seconds % Calculate the total number of frames (maxFrames) maxFrames = floor(duration * frameRate); % Use floor to get an integer frame count % Set the starting frame number startFrame = 1; endFrame=3600; frameNumber = startFrame; % Calculate the time corresponding to the start frame startTime = (startFrame - 1) / frameRate; % Subtract 1 to account for the 0-based index % Set the CurrentTime property of videoReader to start from frame 2000 videoReader.CurrentTime = startTime; % Create a VideoWriter object to save the movies videoWriter = VideoWriter(outputMovieEvents, 'Motion JPEG AVI'); videoWriter.FrameRate = frameRate; % Match the frame rate of the input video open(videoWriter); % Open the VideoWriter object for writing videoWriterTracking = VideoWriter(outputMovieSplitMerge, 'Motion JPEG AVI'); videoWriterTracking.FrameRate = frameRate; % Match the frame rate of the input video open(videoWriterTracking); % Open the VideoWriter object for tracking part % Initialize storage for previous frame labels, areas, and boundaries % Initialize two event queues eventQueue = {}; % For tracking and decision making (this can be modified and cleared) eventQueue2 = {}; % For storing data to populate event_finalized (this remains intact) binaryFrames = {}; % Store binary frames % Initialize counters for highlighted objects and tracks highlightedObjectsFirstLoop = 0; % Count objects highlighted in the first loop highlightedObjectsSecondLoop = 0; % Count objects tracked in the second loop uniqueTracksCount = 0; % Count the number of unique tracks found % Set distance thresholds and boundary proximity distanceThreshold = 200; % pixels for centroid distance threshold (for boundary distance calculation) distanceThresholdb = 7; % pixels for boundary distance (to detect event) minObjectSize = 40; % Minimum size to filter out noise manualThreshold = 0.548; % Threshold for binarization % Prepare for real-time visualization and save frames to the movie % Initialize the frame_event structure frame_event = struct('frame', {}, 'eventID', {}, 'split', {}, 'merge', {}); figure; % Set the figure size to match the video frame size figure('Position', [100, 100, frameWidth, frameHeight]); while hasFrame(videoReader) && frameNumber < endFrame % Read the frame frame = readFrame(videoReader); disp(['Current frame number for detection: ', num2str(frameNumber)]); grayFrame = rgb2gray(frame); % Convert to grayscale if needed % Set the intensity of the selected regions to 255 (maximum) in grayscale frames grayFrame(combinedRegion) = 255; smoothFrame=imgaussfilt(grayFrame,1); flatFieldFrame = imflatfield(smoothFrame, 50); % Normalize background binaryImage = imbinarize(flatFieldFrame, manualThreshold); % Apply threshold binaryImage = ~binaryImage; % Invert binary image if necessary binaryImage = imfill(binaryImage, "holes"); binaryImage = bwareaopen(binaryImage, minObjectSize); % Remove small objects % Store the binary frame for later use binaryFrames{frameNumber} = binaryImage; % Segment objects [B, L] = bwboundaries(binaryImage, 'noholes'); stats = regionprops(L, 'Area', 'Centroid', 'BoundingBox'); % Store the regionprops data propsFrames{frameNumber} = struct('stats', stats, 'boundaries', B); % Clear figure and show the frame imshow(binaryImage); hold on; % Loop over all objects to calculate minimum boundary distance between pairs for i = 1:length(stats) currentObject = stats(i); boundary1 = B{i}; boundaryo1 = B{i}; for j = i+1:length(stats) otherObject = stats(j); boundary2 = B{j}; boundaryo2 = B{j}; % Calculate centroid distance between the two objects centroidDistance = sqrt(sum((currentObject.Centroid - otherObject.Centroid).^2)); % If the centroid distance is below the threshold, calculate boundary distance if centroidDistance < distanceThreshold % Step 1: Calculate boundary-to-boundary minimum distance minBoundaryDistance = min(pdist2(boundary1, boundary2, 'euclidean'), [], 'all'); % Step 2: If the boundary distance is below the threshold, mark as an event if minBoundaryDistance < distanceThresholdb && abs(minBoundaryDistance) > 0 % Highlight both objects plot(boundary1(:,2), boundary1(:,1), 'b', 'LineWidth', 2); % Blue for first object plot(boundary2(:,2), boundary2(:,1), 'r', 'LineWidth', 2); % Red for second object % Increment the counter for highlighted objects in the first loop highlightedObjectsFirstLoop = highlightedObjectsFirstLoop + 2; % Store the event information for potential later analysis newEvent = struct('frame', frameNumber, ... 'centroid1', currentObject.Centroid, ... 'centroid2', otherObject.Centroid, ... 'boundary1', boundary1, ... 'boundary2', boundary2, ... 'boundaryo1', boundaryo1, ... 'boundaryo2', boundaryo2, ... 'object1', i, ... 'object2', j, ... 'trackingFrames', 5, ... % Set the tracking limit to 5 frames 'distanceTrend', minBoundaryDistance, ... % Initial distance 'distanceHistory', minBoundaryDistance, ... % Track distance changes 'tracks1', [], ... % Track of object 1 'tracks2', [], ... % Track of object 2 'areaHistory1', currentObject.Area, ... % Area history for object 1 'areaHistory2', otherObject.Area, ... % Area history for object 2 'pastTracks1', [], ... % For past tracking of object 1 'pastTracks2', [], ... % For past tracking of object 2 'pastAreaHistory1', [], ... % Past area history of object 1 'pastAreaHistory2', [], ... % Past area history of object 2 'pastAreaDiff1', 1, ... % Default past area difference for object 1 'pastAreaDiff2', 1, ... % Default past area difference for object 2 'splitDetected', false, ... 'mergeDetected', false, ... 'disappeared', false, ... 'trackAdded', false, ... 'processed', false, ... 'decisionFrame', NaN); % Initialize decisionFrame as NaN % Add the event to both eventQueue and eventQueue2 eventQueue{end+1} = newEvent; % Used for tracking and decision making eventQueue2{end+1} = newEvent; % This remains intact for finalization % Add the event to frame_event with default split and merge values (0) frame_event(end+1) = struct('frame', frameNumber, 'eventID', length(eventQueue2), 'split', 0, 'merge', 0); end end end end % Capture the frame for the movie frameCapture = getframe(gcf); % Get the current figure frame writeVideo(videoWriter, frameCapture); % Write the frame to the movie hold off; % Increment frame number frameNumber = frameNumber + 1; % Pause to visualize in real-time pause(0.1); % Adjust as needed for your frame rate end % Display the count of objects highlighted in the first loop disp(['Objects highlighted in the first loop: ', num2str(highlightedObjectsFirstLoop)]); disp(frame_event); % After event detection loop % This section will remove duplicate events from eventQueue and then copy the result to eventQueue2 % Step 1: Sort events by frame number [~, sortedIndices] = sort(cellfun(@(x) x.frame, eventQueue)); eventQueue = eventQueue(sortedIndices); % Step 2: Define parameters for spatial and time-based filtering frameDifferenceThreshold = 10; % Maximum allowed difference in frames between events proximityThreshold = 3; % Maximum distance between centroids to consider as duplicate (in pixels) areaDifferenceThreshold = 0.05; % Maximum allowed area difference (5%) % Step 3: Create a map to store events by frame number for fast lookup eventMap = containers.Map('KeyType', 'int32', 'ValueType', 'any'); for idx = 1:length(eventQueue) event = eventQueue{idx}; frame = event.frame; if isKey(eventMap, frame) eventMap(frame) = [eventMap(frame), idx]; else eventMap(frame) = [idx]; end end % Step 4: Initialize an array to keep track of events to keep eventsToKeep = false(1, length(eventQueue)); % Step 5: Loop through each event and check for nearby duplicates for idx = 1:length(eventQueue) if eventsToKeep(idx) continue; % Skip already kept events end event1 = eventQueue{idx}; frame1 = event1.frame; % Loop through a range of frames within the frame difference threshold for frameOffset = 0:frameDifferenceThreshold nearbyFrames = [frame1 - frameOffset, frame1 + frameOffset]; nearbyFrames = nearbyFrames(nearbyFrames >= 1 & nearbyFrames <= max(cellfun(@(x) x.frame, eventQueue))); for frame = nearbyFrames if isKey(eventMap, frame) nearbyEventIndices = eventMap(frame); for nearbyIdx = nearbyEventIndices if nearbyIdx == idx || eventsToKeep(nearbyIdx) continue; % Skip self or already marked events end event2 = eventQueue{nearbyIdx}; % Calculate distances and area differences between possible object pairs dist1 = sqrt(sum((event1.centroid1 - event2.centroid1).^2)); dist2 = sqrt(sum((event1.centroid2 - event2.centroid2).^2)); areaDiff1 = abs(event1.areaHistory1(end) - event2.areaHistory1(end)) / event1.areaHistory1(end); areaDiff2 = abs(event1.areaHistory2(end) - event2.areaHistory2(end)) / event1.areaHistory2(end); dist3 = sqrt(sum((event1.centroid1 - event2.centroid2).^2)); dist4 = sqrt(sum((event1.centroid2 - event2.centroid1).^2)); areaDiff3 = abs(event1.areaHistory1(end) - event2.areaHistory2(end)) / event1.areaHistory1(end); areaDiff4 = abs(event1.areaHistory2(end) - event2.areaHistory1(end)) / event1.areaHistory2(end); % Check if either pair meets the criteria for duplicates isDuplicatePair1 = (dist1 < proximityThreshold && dist2 < proximityThreshold && ... areaDiff1 < areaDifferenceThreshold && areaDiff2 < areaDifferenceThreshold); isDuplicatePair2 = (dist3 < proximityThreshold && dist4 < proximityThreshold && ... areaDiff3 < areaDifferenceThreshold && areaDiff4 < areaDifferenceThreshold); % If a duplicate is found, keep only the later event if isDuplicatePair1 || isDuplicatePair2 if event2.frame > event1.frame eventsToKeep(nearbyIdx) = true; % Keep the later event eventsToKeep(idx) = false; % Mark current event for removal else eventsToKeep(idx) = true; % Keep the current event eventsToKeep(nearbyIdx) = false; % Mark the nearby event for removal end end end end end end % Ensure the current event is marked to be kept if no duplicate was found if ~any(eventsToKeep(nearbyEventIndices)) eventsToKeep(idx) = true; end end % Step 6: Filter out the events that are marked as duplicates % Check the length of eventQueue and eventsToKeep before filtering if length(eventsToKeep) == length(eventQueue) eventQueue = eventQueue(eventsToKeep); else error('Mismatch in lengths: eventsToKeep and eventQueue must have the same length.'); end % Step 7: Create a copy of the filtered eventQueue to eventQueue2 eventQueue2 = eventQueue; disp(['Reduced events from ', num2str(length(sortedIndices)), ' to ', num2str(sum(eventsToKeep)), ' after removing duplicates.']); % Close the first VideoWriter object to finalize the first movie close(videoWriter); % % Save binaryFrames separately % save(fullfile(videoDir, 'binaryFrames.mat'), 'binaryFrames', '-v7.3'); % % % Save propsFrames separately % save(fullfile(videoDir, 'propsFrames.mat'), 'propsFrames', '-v7.3'); save(fullfile(videoDir, 'binaryFrames.mat'), 'binaryFrames', '-v7.3'); % save(fullfile(videoDir, 'propsFrames.mat'), 'propsFrames', '-v7.3'); % % Start tracking into the past (5 frames) for idx = 1:length(eventQueue) event = eventQueue{idx}; eventFrame = event.frame; % Track 5 frames into the past for pastFrame = eventFrame-1:-1:max(eventFrame-5, 1) % Ensure we don't go past frame 1 binaryImagePast = binaryFrames{pastFrame}; % Get the binary image from the past frame [Bpast, Lpast] = bwboundaries(binaryImagePast, 'noholes'); % Segment past objects statspast = regionprops(Lpast, 'Area', 'Centroid', 'BoundingBox'); % Find closest past objects to the ones involved in the event [pastObj1, obj1Dist] = findClosestObject(statspast, event.centroid1, 50); [pastObj2, obj2Dist] = findClosestObject(statspast, event.centroid2, 50); % Handle disappearances in the past and compute area differences if ~pastObj1.disappeared event.pastTracks1 = [event.pastTracks1; pastObj1.Centroid]; % Store past centroid of object 1 event.pastAreaHistory1 = [event.pastAreaHistory1, statspast(pastObj1.index).Area]; % Store past area of object 1 event.pastAreaDiff1 = (event.pastAreaHistory1(1) - event.areaHistory1(end))/event.areaHistory1(end); % Compute area difference for object 1 end if ~pastObj2.disappeared event.pastTracks2 = [event.pastTracks2; pastObj2.Centroid]; % Store past centroid of object 2 event.pastAreaHistory2 = [event.pastAreaHistory2, statspast(pastObj2.index).Area]; % Store past area of object 2 event.pastAreaDiff2 = (event.pastAreaHistory2(1) - event.areaHistory2(end))/event.areaHistory2(end); % Compute area difference for object 2 end end % Update the eventQueue with the past tracking data eventQueue{idx} = event; eventQueue2{idx} = event; % This remains intact for finalization end % Prepare tracking loop frameNumber = startFrame; disp('Starting tracking loop...'); while frameNumber < endFrame % Retrieve the binary frame from the saved list binaryImage = binaryFrames{frameNumber}; disp(['Current frame number for decision: ', num2str(frameNumber)]); % Segment objects in the current frame [B, L] = bwboundaries(binaryImage, 'noholes'); stats = regionprops(L, 'Area', 'Centroid', 'BoundingBox'); % % Clear figure and show the current binary image % imshow(binaryImage); % hold on; % Initialize an array to store the object IDs that were already tracked in this frame trackedObjectIndices = []; % Mark objects that were already processed in earlier events processedObjects = []; % Check event queue to track objects involved in potential split or merge events for idx = 1:length(eventQueue) event = eventQueue{idx}; event2 = eventQueue2{idx}; % Also retrieve the corresponding event2 if eventQueue{idx}.processed continue; end %disp(['merge detected in frame ', num2str(event2.frame)]); % Ensure that the objects haven't already been processed in another event if ismember(event.object1, processedObjects) || ismember(event.object2, processedObjects) continue; % Skip this event if the objects were processed end % Track for a maximum of `trackingFrames` frames after the event if frameNumber <= event.frame + event.trackingFrames % Find the closest objects to the tracked ones from the previous frame [nextObj1, obj1Dist] = findClosestObject(stats, event.centroid1, 50); % Adjusted threshold [nextObj2, obj2Dist] = findClosestObject(stats, event.centroid2, 50); % Adjusted threshold % Count unique tracks only once per event if ~event.trackAdded uniqueTracksCount = uniqueTracksCount + 1; event.trackAdded = true; % Mark the event so it's counted only once end % Only count objects if they haven't been tracked in this frame already if ~ismember(nextObj1.index, trackedObjectIndices) highlightedObjectsSecondLoop = highlightedObjectsSecondLoop + 1; trackedObjectIndices = [trackedObjectIndices, nextObj1.index]; % Mark as tracked processedObjects = [processedObjects, nextObj1.index]; % Mark as processed end if ~ismember(nextObj2.index, trackedObjectIndices) highlightedObjectsSecondLoop = highlightedObjectsSecondLoop + 1; trackedObjectIndices = [trackedObjectIndices, nextObj2.index]; % Mark as tracked processedObjects = [processedObjects, nextObj2.index]; % Mark as processed end % If both objects are still present, update their centroids and boundaries if ~nextObj1.disappeared && ~nextObj2.disappeared % Update event with the new centroids, boundaries, and area histories event.centroid1 = nextObj1.Centroid; event.centroid2 = nextObj2.Centroid; event.boundary1 = B{nextObj1.index}; event.boundary2 = B{nextObj2.index}; event.areaHistory1 = [event.areaHistory1, stats(nextObj1.index).Area]; event.areaHistory2 = [event.areaHistory2, stats(nextObj2.index).Area]; % Also update event2 % Update event2 with the new centroids, boundaries, and area histories event2.centroid1 = nextObj1.Centroid; event2.centroid2 = nextObj2.Centroid; event2.boundary1 = B{nextObj1.index}; event2.boundary2 = B{nextObj2.index}; event2.areaHistory1 = [event2.areaHistory1, stats(nextObj1.index).Area]; event2.areaHistory2 = [event2.areaHistory2, stats(nextObj2.index).Area]; % Update tracks in event2 for future velocity calculation event2.tracks1 = [event2.tracks1; nextObj1.Centroid]; % Store the centroid for future velocity event2.tracks2 = [event2.tracks2; nextObj2.Centroid]; % Store the centroid for future velocity % Calculate the new boundary-to-boundary minimum distance newBoundaryDistance = min(pdist2(event.boundary1, event.boundary2, 'euclidean'), [], 'all'); % Append the new distance to the history event.distanceHistory = [event.distanceHistory, newBoundaryDistance]; event2.distanceHistory = [event2.distanceHistory, newBoundaryDistance]; % Update event2 % Check if objects disappeared or at the final frame for this event if frameNumber == event.frame + event.trackingFrames || event.disappeared distanceChanges = diff(event.distanceHistory); % Calculate changes in distance numIncreases = sum(distanceChanges > 0.5); % Count increases numDecreases = sum(distanceChanges < -0.5); % Count decreases numStay = sum(-0.5 < distanceChanges & distanceChanges < 0.5); % Count stays % Check if area of the larger object has decreased over past frames if event.areaHistory1(end) > event.areaHistory2(end) pastAreaDiff = event.pastAreaDiff1; % Use object 1's area difference else pastAreaDiff = event.pastAreaDiff2; % Use object 2's area difference end % Determine if it's a split or merge based on the trend and area change if numIncreases/(numDecreases+numIncreases+numStay) > 0.9 && pastAreaDiff > 0.1 event.splitDetected = true; % It's a split event2.splitDetected = true; % It's a split % Highlight the objects in magenta for split plot(event.boundary1(:,2), event.boundary1(:,1), 'm', 'LineWidth', 2); plot(event.boundary2(:,2), event.boundary2(:,1), 'm', 'LineWidth', 2); event2.decisionFrame = frameNumber; % Store the frame where the decision is made % Update frame_event split status for this eventID frame_event([frame_event.eventID] == idx).split = 1; elseif numDecreases/(numDecreases+numIncreases) > 0.5 event.mergeDetected = true; % It's a merge event2.mergeDetected = true; % It's a merge % Highlight the objects in yellow for merge plot(event.boundary1(:,2), event.boundary1(:,1), 'g', 'LineWidth', 2); plot(event.boundary2(:,2), event.boundary2(:,1), 'g', 'LineWidth', 2); event2.decisionFrame = frameNumber; % Store the frame where the decision is made % Update frame_event merge status for this eventID frame_event([frame_event.eventID] == idx).merge = 1; % If a merge is detected, save the boundary at the decision frame if event.mergeDetected && frameNumber == event.decisionFrame event2.boundary1 = B{nextObj1.index}; % Save boundary for merge object 1 event2.boundary2 = B{nextObj2.index}; % Save boundary for merge object 2 end end % After processing, remove the event immediately from eventQueue %eventQueue{idx} = []; % Mark event for removal (but keep in event2) eventQueue{idx}.processed = true; end else % Mark as disappeared if either object is gone event.disappeared = true; event2.disappeared = true; % Update event2 end end % Update eventQueue2 with the changes to event2 eventQueue2{idx} = event2; end % Remove processed events from the queue eventQueue = eventQueue(~cellfun('isempty', eventQueue)); % Clean up the event queue % % Capture the frame for the tracking movie % frameCapture = getframe(gcf); % Get the current figure frame % writeVideo(videoWriterTracking, frameCapture); % Write the frame to the tracking movie % % hold off; % Increment frame number frameNumber = frameNumber + 1; % Pause to visualize in real-time pause(0.1); % Adjust as needed for your frame rate end % Preprocess events and group them by frames for optimized plotting eventMap = containers.Map('KeyType', 'int32', 'ValueType', 'any'); % Map to store events by frame % Group split and merge events by the frame they should be plotted for idx = 1:length(eventQueue2) event2 = eventQueue2{idx}; % Determine the frame where the event should be plotted if event2.splitDetected plotFrame = event2.frame; % Split is plotted at the event frame elseif event2.mergeDetected plotFrame = event2.decisionFrame; % Merge is plotted at the decision frame else continue; % No split or merge detected for this event end % Add event to the eventMap for the corresponding frame if isKey(eventMap, plotFrame) eventMap(plotFrame) = [eventMap(plotFrame), event2]; % Append to existing events for this frame else eventMap(plotFrame) = {event2}; % Add the first event for this frame end end % Optimized loop for plotting split and merge events frame-by-frame for frameNumber = startFrame:endFrame % Adjust the number of frames processed as needed % Retrieve the binary image for the current frame if frameNumber > 0 && frameNumber < endFrame binaryImage = binaryFrames{frameNumber}; % Segment the objects in the current frame (necessary to get updated boundaries) [B, L] = bwboundaries(binaryImage, 'noholes'); stats = regionprops(L, 'Centroid', 'Area'); % Retrieve centroids and areas % Display the binary image for each frame imshow(binaryImage); hold on; % Check if there are events (split or merge) for this frame if isKey(eventMap, frameNumber) % Retrieve the events for this frame eventsToPlot = eventMap(frameNumber); % Ensure eventsToPlot is treated as a cell array if ~iscell(eventsToPlot) eventsToPlot = {eventsToPlot}; % Wrap single events in a cell array end % Plot all events for this frame for i = 1:numel(eventsToPlot) event2 = eventsToPlot{i}; % Extract each event % For split, use boundary from event detection frame if isfield(event2, 'splitDetected') && event2.splitDetected boundary1 = event2.boundaryo1; boundary2 = event2.boundaryo2; % Plot in magenta for split plot(boundary1(:, 2), boundary1(:, 1), 'm', 'LineWidth', 4); plot(boundary2(:, 2), boundary2(:, 1), 'm', 'LineWidth', 4); % For merge, use the saved boundary from the decision frame elseif isfield(event2, 'mergeDetected') && event2.mergeDetected boundary1 = event2.boundary1; % Saved boundary for merge object 1 boundary2 = event2.boundary2; % Saved boundary for merge object 2 % Plot in green for merge plot(boundary1(:, 2), boundary1(:, 1), 'g', 'LineWidth', 4); plot(boundary2(:, 2), boundary2(:, 1), 'g', 'LineWidth', 4); end end end % Capture the frame for the movie frameCapture = getframe(gcf); % Get the current figure frame writeVideo(videoWriterTracking, frameCapture); % Write the frame to the tracking movie hold off; end % Increment frame number and reduce pause time to prevent slowness pause(0.01); % Adjust as needed for your frame rate end % Close the video writer close(videoWriterTracking); % Display the count of objects tracked in the second loop disp(['Objects tracked in the second loop: ', num2str(highlightedObjectsSecondLoop)]); % Display the number of unique tracks found disp(['Unique tracks found in the second loop: ', num2str(uniqueTracksCount)]); % Initialize new structure to store finalized events event_finalized = []; % Loop through eventQueue2 (instead of eventQueue) to finalize events for idx = 1:length(eventQueue2) event2 = eventQueue2{idx}; % Use event2 (unaltered data) for finalization eventFrame = event2.frame; disp(['frame: ', num2str(eventFrame)]); % Initialize finalized event structure for this event finalizedEvent = struct('frame', eventFrame, ... 'object1', event2.object1, ... 'object2', event2.object2, ... 'area1_before', [], ... 'area1_after', [], ... 'area1_at_event', event2.areaHistory1(end), ... 'circularity1_before', [], ... 'circularity1_after', [], ... 'circularity1_at_event', [], ... 'velocity1_before', [], ... 'velocity1_after', [], ... 'split', 0, ... 'merge', 0, ... 'area2_before', [], ... 'area2_after', [], ... 'area2_at_event', event2.areaHistory2(end), ... 'circularity2_before', [], ... 'circularity2_after', [], ... 'circularity2_at_event', [], ... 'velocity2_before', [], ... 'velocity2_after', []); %% Debugging Information disp(['Processing event ' num2str(idx) ' at frame ' num2str(eventFrame)]); % Calculate circularity for objects at event frame (current) if ~isempty(event2.boundary1) finalizedEvent.circularity1_at_event = 4 * pi * event2.areaHistory1(end) / (perimeter(event2.boundary1)^2); else disp('Warning: boundary1 is empty, cannot calculate circularity for object 1.'); end if ~isempty(event2.boundary2) finalizedEvent.circularity2_at_event = 4 * pi * event2.areaHistory2(end) / (perimeter(event2.boundary2)^2); else disp('Warning: boundary2 is empty, cannot calculate circularity for object 2.'); end %% Object 1 (Past Data) if ~isempty(event2.pastTracks1) area1_before = event2.pastAreaHistory1(end); % Last frame before the event finalizedEvent.area1_before = area1_before; % Use the event frame centroid and the last past track for velocity velocity1_before = sqrt(sum((event2.centroid1 - event2.pastTracks1(end,:)).^2)); % Compare event frame to the closest past frame finalizedEvent.velocity1_before = velocity1_before; else disp('Warning: Not enough pastTracks1 data for velocity1_before calculation.'); end %% Object 1 (Future Data) if ~isempty(event2.tracks1) area1_after = event2.areaHistory1(1); % First frame after the event finalizedEvent.area1_after = area1_after; % Use the event frame centroid and the first future track for velocity velocity1_after = sqrt(sum((event2.centroid1 - event2.tracks1(1,:)).^2)); % Compare event frame to the closest future frame finalizedEvent.velocity1_after = velocity1_after; else disp('Warning: Not enough tracks1 data for velocity1_after calculation.'); end %% Object 2 (Past Data) if ~isempty(event2.pastTracks2) area2_before = event2.pastAreaHistory2(end); % Last frame before the event finalizedEvent.area2_before = area2_before; % Use the event frame centroid and the last past track for velocity velocity2_before = sqrt(sum((event2.centroid2 - event2.pastTracks2(end,:)).^2)); % Compare event frame to the closest past frame finalizedEvent.velocity2_before = velocity2_before; else disp('Warning: Not enough pastTracks2 data for velocity2_before calculation.'); end %% Object 2 (Future Data) if ~isempty(event2.tracks2) area2_after = event2.areaHistory2(1); % First frame after the event finalizedEvent.area2_after = area2_after; % Use the event frame centroid and the first future track for velocity velocity2_after = sqrt(sum((event2.centroid2 - event2.tracks2(1,:)).^2)); % Compare event frame to the closest future frame finalizedEvent.velocity2_after = velocity2_after; else disp('Warning: Not enough tracks2 data for velocity2_after calculation.'); end % Set the split/merge flags based on event2's decision if event2.splitDetected finalizedEvent.split = 1; elseif event2.mergeDetected finalizedEvent.merge = 1; else % If neither split nor merge detected, set all properties to NaN finalizedEvent.area1_before = NaN; finalizedEvent.area1_after = NaN; finalizedEvent.area2_before = NaN; finalizedEvent.area2_after = NaN; finalizedEvent.circularity1_before = NaN; finalizedEvent.circularity1_after = NaN; finalizedEvent.circularity2_before = NaN; finalizedEvent.circularity2_after = NaN; finalizedEvent.velocity1_before = NaN; finalizedEvent.velocity1_after = NaN; finalizedEvent.velocity2_before = NaN; finalizedEvent.velocity2_after = NaN; end % Add the finalized event to the event_finalized list event_finalized = [event_finalized; finalizedEvent]; end % Define the maximum frame number from the frame_event structure maxFrame = max([frame_event.frame]); % Initialize arrays for counting split and merge events per exact frame splitCounts = zeros(1, maxFrame); mergeCounts = zeros(1, maxFrame); % Loop through frame_event to count split and merge events per frame for idx = 1:length(frame_event) frame = frame_event(idx).frame; % Ensure frame is within the valid range if frame > 0 && frame <= maxFrame if frame_event(idx).split == 1 splitCounts(frame) = splitCounts(frame) + 1; end if frame_event(idx).merge == 1 mergeCounts(frame) = mergeCounts(frame) + 1; end else disp(['Event at frame ', num2str(frame), ' is out of range.']); end end % Plot the split counts as a function of frame number figure; bar(1:maxFrame, splitCounts, 'FaceColor', 'blue'); xlabel('Frame Number'); ylabel('Split Event Count'); title('Count of Split Events per Frame'); grid on; % Save the figure in the same folder as the video file saveas(gcf, fullfile(videoDir, 'SplitCountsByFrame.png')); % Plot the merge counts as a function of frame number figure; bar(1:maxFrame, mergeCounts, 'FaceColor', 'red'); xlabel('Frame Number'); ylabel('Merge Event Count'); title('Count of Merge Events per Frame'); grid on; % Save the figure in the same folder as the video file saveas(gcf, fullfile(videoDir, 'MergeCountsByFrame.png')); % Save eventQueue, eventQueue2, event_finalized, and frame_event data to a .mat file save('FilePath\event_data.mat', ... 'eventQueue', 'eventQueue2', 'event_finalized', 'frame_event'); % Bin frame counts by 10 frames and calculate split, merge, and total event counts binSize = 10; numBins = ceil(maxFrame / binSize); binnedSplitCounts = zeros(1, numBins); binnedMergeCounts = zeros(1, numBins); binnedEventCounts = zeros(1, numBins); % Loop through each frame to accumulate counts into bins for frame = 1:maxFrame binIndex = ceil(frame / binSize); binnedSplitCounts(binIndex) = binnedSplitCounts(binIndex) + splitCounts(frame); binnedMergeCounts(binIndex) = binnedMergeCounts(binIndex) + mergeCounts(frame); binnedEventCounts(binIndex) = binnedSplitCounts(binIndex) + binnedMergeCounts(binIndex); end % Plot binned split counts figure; bar(1:numBins, binnedSplitCounts, 'FaceColor', 'blue'); xlabel('Binned Frame Number (x10)'); ylabel('Split Event Count'); title('Count of Split Events per 10 Frames'); grid on; % Save the figure saveas(gcf, fullfile(videoDir, 'BinnedSplitCounts.png')); % Plot binned merge counts figure; bar(1:numBins, binnedMergeCounts, 'FaceColor', 'red'); xlabel('Binned Frame Number (x10)'); ylabel('Merge Event Count'); title('Count of Merge Events per 10 Frames'); grid on; % Save the figure saveas(gcf, fullfile(videoDir, 'BinnedMergeCounts.png')); % Plot binned total event counts (split + merge) figure; bar(1:numBins, binnedEventCounts, 'FaceColor', 'magenta'); xlabel('Binned Frame Number (x10)'); ylabel('Total Event Count'); title('Count of Total Events per 10 Frames'); grid on; % Save the figure saveas(gcf, fullfile(videoDir, 'BinnedEventCounts.png')); % Helper function to calculate perimeter from boundary points function peri = perimeter(boundary) peri = sum(sqrt(sum(diff(boundary).^2, 2))); end % Helper function to find the closest object within a given threshold function [closestObj, distance] = findClosestObject(stats, centroid, threshold) closestObj = []; distance = inf; for i = 1:length(stats) currentDist = sqrt(sum((stats(i).Centroid - centroid).^2)); % Adjust the threshold based on object area if stats(i).Area > 1000 currentThreshold = 300; % Increase threshold for larger objects else currentThreshold = threshold; % Use the original threshold for smaller objects end if currentDist < currentThreshold && currentDist < distance closestObj = struct('index', i, 'Centroid', stats(i).Centroid, 'disappeared', false); distance = currentDist; end end % If no close object is found, mark it as disappeared if isempty(closestObj) closestObj = struct('index', -1, 'Centroid', [], 'disappeared', true); end end