Background: I have recorded navigation trajectory positional x,y and its time point (50hz) data from a virtual reality navigation in an empty room.

Problem: I want to calculate how much the navigation deviates from it's closest point to the ideal path and must calculate to the correct segment they are still "on"/furthest along the path currently.

Inputs:

Current output using KNN() only (upper left subplot) of a good (upper figure) and bad (lower) navigation. Note how on bad navigations it goes crazy and compares to any random point closest to it when it needs to compare to the current path segment it is on:

The current method using function KNN() to find the nearest neighbor works well sometimes but it will not work in cases where navigations are completely wrong and it begins comparing deviations to a whole different path it is not even on yet. It also entirely skips corners as the hypotenuse euclidean distance always going to be the furthest... Code blocks to process both images and the calculation/plotting below:

% Single line shortest path analyses                    
shortstimImg = imread(shorteststimFile);
if size(stimImg, 3) == 3
    shortstimImg = rgb2gray(shortstimImg);
end
% Convert the image to double precision
shortstimImg = im2double(shortstimImg);                 
% Resize the image to 600x600
shortstimImgResized = imresize(shortstimImg, [600, 600], 'nearest');

% Convert to binary and skeletonize
if size(shortstimImgResized, 3) == 3
    shorteststimImgGray = rgb2gray(shortstimImgResized);
else
    shorteststimImgGray = shortstimImgResized;
end

% Binarize and skeletonize

% Convert to binary image
binary = imbinarize(shorteststimImgGray);  % Automatic thresholding
% binary = imbinarize(shorteststimImgBW, 'adaptive'); % Alternative adaptive threshold

% Ensure white lines on black background (check top-left corner)
if mean(binary(1:10,1:10), 'all') > 0.5
    binary = ~binary;  % Invert image
end

% Clean up with morphological closing
kernel = strel('square', 3);  % 3x3 structuring element
closed = imclose(binary, kernel);

% Skeletonization with branch cleaning
skeleton = bwskel(closed, 'MinBranchLength', 20);  % Remove small branches
skeleton = bwmorph(skeleton, 'bridge');  % Connect nearby segments

% Convert to uint8 image (255 for lines, 0 for background)
skeleton_uint8 = uint8(skeleton) * 255;

% Extract coordinates from 1-pixel skeleton BEFORE any dilation
[stimY, stimX] = find(skeleton);  % MATLAB returns [rows, cols]

% Convert coordinates to 600x600 cm space (preserve aspect ratio)
[height, width] = size(skeleton);
stimX_cm = (stimX / width) * 600;
stimY_cm = (stimY / height) * 600;
idealTrajectory = [stimX_cm, stimY_cm];

% Read the PNG image from the matfile
matImg = imread(rawTrajectory.png);
% Convert the image to grayscale if it's not already in grayscale format
if size(matImg, 3) == 3
    matImg = rgb2gray(matImg);
end
% Convert the image to double precision
matImg = im2double(matImg);
% Convert the matImg to uint8
matImgUint8 = im2uint8(matImg);  
% Extract trajectory coordinates
skeletonizedImg = bwmorph(matImgUint8 == 0, 'skel', Inf);
[trajY, trajX] = find(skeletonizedImg);
trajectory = [trajX, trajY];

function plotPathDeviation(metrics, idealTrajectory, trajectory_raw)
    % Add input validation
    if size(trajectory_raw, 2) ~= 2
        error('Trajectory must be Nx2 matrix');
    end
    
    % Clean trajectory data
    validIdx = all(~isnan(trajectory_raw), 2);
    trajectory_valid = trajectory_raw(validIdx,:);
    
    % Ensure idealTrajectory is Nx2
    idealTrajectory = validateTrajectoryDimensions(idealTrajectory);
    % Set white background
    set(gca, 'Color', 'w');
    
    % Plot ideal path
    scatter(idealTrajectory(:,1), idealTrajectory(:,2), 10, 'g', 'DisplayName', 'Ideal Path');
    hold on;
    
    % Plot trajectory
    validIdx = all(~isnan(trajectory_raw), 2);
    trajectory_valid = trajectory_raw(validIdx, :);
    scatter(trajectory_valid(:,1), trajectory_valid(:,2), 10, 'b', 'filled', 'DisplayName', 'Trajectory');
    
    % Sample points for deviation lines
    sample_points = 1:size(trajectory_valid, 1);
    
    % Compute KNN distances for visualization
    [nearestIndices, ~] = knnsearch(idealTrajectory, trajectory_valid);
    
    maxDev = 0;
    meanDev = 0;
    count = 0;
    
    for i = sample_points
        trajPoint = trajectory_valid(i,:);
        nearestPoint = idealTrajectory(nearestIndices(i),:);
        
        % Plot deviation line
        plot([trajPoint(1) nearestPoint(1)], [trajPoint(2) nearestPoint(2)], ...
            'Color', [1 0 0 0.2], ... % Red with 0.2 alpha (80% transparent)
            'LineWidth', 0.75, ...     % Slightly thicker to compensate for transparency
            'HandleVisibility', 'off');
        
        % Update deviation metrics
        dist = norm(trajPoint - nearestPoint);
        maxDev = max(maxDev, dist);
        meanDev = meanDev + dist;
        count = count + 1;
    end
    
    meanDev = meanDev / (count + eps);
    title(sprintf('Path Deviation\nMean Dev: %.1f cm, Max Dev: %.1f cm', meanDev, maxDev));
    grid on;
    configureAxes();
    hold off;
end

Suggestion perhaps: Expanding/growing radii circle on each pixel point on one of the images until it matches/touches the closest comparison image point? Not sure if it's best to expand it from the template ideal path or the recorded trajectory. Leaning towards recorded trajectory as that always varies person-to-person

Will need to implement a sequential pathing processing so comparisons do not jump ahead to compare to a path further ahead.

I have each segment composing the entire ideal path classified of each segment's endpoints, orientation, angle, length.

If there are any other methods you guys think it best for my scenario, please let me know! Many thanks all.