function roiData = RelativePos(roiData, Model)
%RELATIVEPOS Computes the relative position differences between ROIs.
%
%   DISTANCETABLE = RELATIVEPOS(ROIBASICFEATURES) calculates the 
%   triplet-wise differences in centroid coordinates (xcm, ycm, zcm) 
%   for each region of interest (ROI) in the input structure.
%
%   Input:
%       ROIBASICFEATURES - A struct where each field corresponds to an ROI
%                          and contains 'xcm', 'ycm', and 'zcm' values.
%
%   Output:
%       DISTANCETABLE - A table containing the coordinate differences
%                       between each pair of ROIs.
%
%   Example:
%       ROIBasicFeatures.ROI1 = struct('xcm', 10, 'ycm', 15, 'zcm', 20);
%       ROIBasicFeatures.ROI2 = struct('xcm', 5, 'ycm', 10, 'zcm', 15);
%       distanceTable = relativePos(ROIBasicFeatures);
       
    if nargin < 2
        Model = NaN;
    end

    roiDataTable = struct2table(roiData);

    roiDataTable.originalRowIndex = (1:height(roiDataTable))'; % this is done to keep track of the original index in the datatable later on

    if isstruct(Model)
        roiNames = fieldnames(Model); % calculate the relative Pos for the ROIs that are in the Model
    else
       roiNames = unique(roiDataTable.ROIName);
    end
    
    nRoi = numel(roiNames);
    nPatient = max([roiDataTable.PtNumber]);

    relPosTables = cell(height(roiDataTable), 1);
    relPosToRoiName = compose("relPos_%s", string(roiNames));

    % Extract coordinates or assign NaNs if missing
    for iPt = 1:nPatient
        ptMask = roiDataTable.PtNumber == iPt;
        ptData = roiDataTable(ptMask, ["xcm", "ycm", "zcm", "ROIName", "originalRowIndex"]);
        
        % Extract coordinates into a matrix
        coords = NaN(nRoi, 3);
        roiIdxMap = zeros(nRoi, 1);

        for kRoi = 1:nRoi
            roiMatch = strcmp(ptData.ROIName, roiNames{kRoi});
            if any(roiMatch)
                coords(kRoi, :) = ptData{roiMatch, {'xcm', 'ycm', 'zcm'}};
                roiIdxMap(kRoi) = find(roiMatch); % row in ptData
            else
                coords(kRoi, :) = NaN; % Fill with NaNs if ROI is missing
            end
        end
        
        % Compute relative positions
        relPos = cell(nRoi, 1);
        for iRoi = 1:nRoi
            diffs = num2cell(coords(iRoi, :) - coords, 2)';
            relPosTable = cell2table(diffs, 'VariableNames', relPosToRoiName);

            if roiIdxMap(iRoi) > 0
                relPosTables{ptData(strcmp(ptData.ROIName, roiNames{iRoi}), :).originalRowIndex} = relPosTable; 
            end
        end
    end
     
    % Combine all relPos tables into one big table
    flatRelPosTables = vertcat(relPosTables{:});
    roiDataTable = [roiDataTable, flatRelPosTables];
    roiData = table2struct(roiDataTable);
end
